ITS/AliITSMultReconstructor.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //_________________________________________________________________________
  17 //
  18 //        Implementation of the ITS-SPD trackleter class
  19 //
  20 // It retrieves clusters in the pixels (theta and phi) and finds tracklets.
  21 // These can be used to extract charged particle multiplicity from the ITS.
  22 //
  23 // A tracklet consists of two ITS clusters, one in the first pixel layer and
  24 // one in the second. The clusters are associated if the differences in
  25 // Phi (azimuth) and Theta (polar angle) are within fiducial windows.
  26 // In case of multiple candidates the candidate with minimum
  27 // distance is selected.
  28 //
  29 // Two methods return the number of tracklets and the number of unassociated
  30 // clusters (i.e. not used in any tracklet) in the first SPD layer
  31 // (GetNTracklets and GetNSingleClusters)
  32 //
  33 // The cuts on phi and theta depend on the interacting system (p-p or Pb-Pb)
  34 // and can be set via AliITSRecoParam class
  35 // (SetPhiWindow and SetThetaWindow)
  36 //
  37 // Origin: Tiziano Virgili
  38 //
  39 // Current support and development:
  40 //         Domenico Elia, Maria Nicassio (INFN Bari)
  41 //         Domenico.Elia@ba.infn.it, Maria.Nicassio@ba.infn.it
  42 //
  43 // Most recent updates:
  44 //     - multiple association forbidden (fOnlyOneTrackletPerC2 = kTRUE)
  45 //     - phi definition changed to ALICE convention (0,2*TMath::pi())
  46 //     - cluster coordinates taken with GetGlobalXYZ()
  47 //     - fGeometry removed
  48 //     - number of fired chips on the two layers
  49 //     - option to cut duplicates in the overlaps
  50 //     - options and fiducial cuts via AliITSRecoParam
  51 //     - move from DeltaZeta to DeltaTheta cut
  52 //     - update to the new algorithm by Mariella and Jan Fiete
  53 //     - store also DeltaTheta in the ESD
  54 //     - less new and delete calls when creating the needed arrays
  55 //
  56 //     - RS: to decrease the number of new/deletes the clusters data are stored
  57 //           not in float[6] attached to float**, but in 1-D array.
  58 //     - RS: Clusters are sorted in Z in roder to have the same numbering as in the ITS reco
  59 //     - RS: Clusters used by ESDtrack are flagged, this information is passed to AliMulitiplicity object
  60 //           when storing the tracklets and single cluster info
  61 //     - MN: first MC label of single clusters stored
  62 //_________________________________________________________________________
  63
  64 #include <TClonesArray.h>
  65 #include <TH1F.h>
  66 #include <TH2F.h>
  67 #include <TTree.h>
  68 #include <TBits.h>
  69 #include <TArrayI.h>
  70 #include <string.h>
  71
  72 #include "AliITSMultReconstructor.h"
  73 #include "AliITSReconstructor.h"
  74 #include "AliITSRecPoint.h"
  75 #include "AliITSRecPointContainer.h"
  76 #include "AliITSgeom.h"
  77 #include "AliITSgeomTGeo.h"
  78 #include "AliITSDetTypeRec.h"
  79 #include "AliESDEvent.h"
  80 #include "AliESDVertex.h"
  81 #include "AliESDtrack.h"
  82 #include "AliMultiplicity.h"
  83 #include "AliLog.h"
  84 #include "TGeoGlobalMagField.h"
  85 #include "AliMagF.h"
  86 #include "AliESDv0.h"
  87 #include "AliV0.h"
  88 #include "AliKFParticle.h"
  89 #include "AliKFVertex.h"
  90 #include "AliRefArray.h"
  91
  92 //____________________________________________________________________
  93 ClassImp(AliITSMultReconstructor)
  94
  95
  96 //____________________________________________________________________
  97 AliITSMultReconstructor::AliITSMultReconstructor():
  98 fDetTypeRec(0),fESDEvent(0),fTreeRP(0),fTreeRPMix(0),
  99 fTracklets(0),
 100 fSClusters(0),
 101 fNTracklets(0),
 102 fNSingleCluster(0),
 103 fDPhiWindow(0),
 104 fDThetaWindow(0),
 105 fPhiShift(0),
 106 fRemoveClustersFromOverlaps(0),
 107 fPhiOverlapCut(0),
 108 fZetaOverlapCut(0),
 109 fPhiRotationAngle(0),
 110 fScaleDTBySin2T(0),
 111 fNStdDev(1.0),
 112 fNStdDevSq(1.0),
 113 //
 114 fCutPxDrSPDin(0.1),
 115 fCutPxDrSPDout(0.15),
 116 fCutPxDz(0.2),
 117 fCutDCArz(0.5),
 118 fCutMinElectronProbTPC(0.5),
 119 fCutMinElectronProbESD(0.1),
 120 fCutMinP(0.05),
 121 fCutMinRGamma(2.),
 122 fCutMinRK0(1.),
 123 fCutMinPointAngle(0.98),
 124 fCutMaxDCADauther(0.5),
 125 fCutMassGamma(0.03),
 126 fCutMassGammaNSigma(5.),
 127 fCutMassK0(0.03),
 128 fCutMassK0NSigma(5.),
 129 fCutChi2cGamma(2.),
 130 fCutChi2cK0(2.),
 131 fCutGammaSFromDecay(-10.),
 132 fCutK0SFromDecay(-10.),
 133 fCutMaxDCA(1.),
 134 //
 135 fHistOn(0),
 136 fhClustersDPhiAcc(0),
 137 fhClustersDThetaAcc(0),
 138 fhClustersDPhiAll(0),
 139 fhClustersDThetaAll(0),
 140 fhDPhiVsDThetaAll(0),
 141 fhDPhiVsDThetaAcc(0),
 142 fhetaTracklets(0),
 143 fhphiTracklets(0),
 144 fhetaClustersLay1(0),
 145 fhphiClustersLay1(0),
 146 //
 147   fDPhiShift(0),
 148   fDPhiWindow2(0),
 149   fDThetaWindow2(0),
 150   fPartners(0),
 151   fAssociatedLay1(0),
 152   fMinDists(0),
 153   fBlackList(0),
 154 //
 155   fCreateClustersCopy(0),
 156   fClustersLoaded(0),
 157   fRecoDone(0),
 158   fBuildRefs(kTRUE),
 159   fSPDSeg()
 160 {
 161   // default c-tor
 162   for (int i=0;i<2;i++) {
 163     fNFiredChips[i] = 0;
 164     fClArr[i] = 0;
 165     for (int j=0;j<2;j++) fUsedClusLay[i][j] = 0;
 166     fDetectorIndexClustersLay[i] = 0;
 167     fClusterCopyIndex[i] = 0;
 168     fOverlapFlagClustersLay[i] = 0;
 169     fNClustersLay[i] = 0;
 170     fClustersLay[i] = 0;
 171   }
 172   // Method to reconstruct the charged particles multiplicity with the
 173   // SPD (tracklets).
 174
 175   SetHistOn();
 176
 177   if (AliITSReconstructor::GetRecoParam()) {
 178     SetPhiWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiWindow());
 179     SetThetaWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterThetaWindow());
 180     SetPhiShift(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiShift());
 181     SetRemoveClustersFromOverlaps(AliITSReconstructor::GetRecoParam()->GetTrackleterRemoveClustersFromOverlaps());
 182     SetPhiOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiOverlapCut());
 183     SetZetaOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterZetaOverlapCut());
 184     SetPhiRotationAngle(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiRotationAngle());
 185     SetNStdDev(AliITSReconstructor::GetRecoParam()->GetTrackleterNStdDevCut());
 186     SetScaleDThetaBySin2T(AliITSReconstructor::GetRecoParam()->GetTrackleterScaleDThetaBySin2T());
 187     SetBuildRefs(AliITSReconstructor::GetRecoParam()->GetTrackleterBuildCl2TrkRefs());
 188     //
 189     SetCutPxDrSPDin(AliITSReconstructor::GetRecoParam()->GetMultCutPxDrSPDin());
 190     SetCutPxDrSPDout(AliITSReconstructor::GetRecoParam()->GetMultCutPxDrSPDout());
 191     SetCutPxDz(AliITSReconstructor::GetRecoParam()->GetMultCutPxDz());
 192     SetCutDCArz(AliITSReconstructor::GetRecoParam()->GetMultCutDCArz());
 193     SetCutMinElectronProbTPC(AliITSReconstructor::GetRecoParam()->GetMultCutMinElectronProbTPC());
 194     SetCutMinElectronProbESD(AliITSReconstructor::GetRecoParam()->GetMultCutMinElectronProbESD());
 195     SetCutMinP(AliITSReconstructor::GetRecoParam()->GetMultCutMinP());
 196     SetCutMinRGamma(AliITSReconstructor::GetRecoParam()->GetMultCutMinRGamma());
 197     SetCutMinRK0(AliITSReconstructor::GetRecoParam()->GetMultCutMinRK0());
 198     SetCutMinPointAngle(AliITSReconstructor::GetRecoParam()->GetMultCutMinPointAngle());
 199     SetCutMaxDCADauther(AliITSReconstructor::GetRecoParam()->GetMultCutMaxDCADauther());
 200     SetCutMassGamma(AliITSReconstructor::GetRecoParam()->GetMultCutMassGamma());
 201     SetCutMassGammaNSigma(AliITSReconstructor::GetRecoParam()->GetMultCutMassGammaNSigma());
 202     SetCutMassK0(AliITSReconstructor::GetRecoParam()->GetMultCutMassK0());
 203     SetCutMassK0NSigma(AliITSReconstructor::GetRecoParam()->GetMultCutMassK0NSigma());
 204     SetCutChi2cGamma(AliITSReconstructor::GetRecoParam()->GetMultCutChi2cGamma());
 205     SetCutChi2cK0(AliITSReconstructor::GetRecoParam()->GetMultCutChi2cK0());
 206     SetCutGammaSFromDecay(AliITSReconstructor::GetRecoParam()->GetMultCutGammaSFromDecay());
 207     SetCutK0SFromDecay(AliITSReconstructor::GetRecoParam()->GetMultCutK0SFromDecay());
 208     SetCutMaxDCA(AliITSReconstructor::GetRecoParam()->GetMultCutMaxDCA());
 209     //
 210   } else {
 211     SetPhiWindow();
 212     SetThetaWindow();
 213     SetPhiShift();
 214     SetRemoveClustersFromOverlaps();
 215     SetPhiOverlapCut();
 216     SetZetaOverlapCut();
 217     SetPhiRotationAngle();
 218
 219     //
 220     SetCutPxDrSPDin();
 221     SetCutPxDrSPDout();
 222     SetCutPxDz();
 223     SetCutDCArz();
 224     SetCutMinElectronProbTPC();
 225     SetCutMinElectronProbESD();
 226     SetCutMinP();
 227     SetCutMinRGamma();
 228     SetCutMinRK0();
 229     SetCutMinPointAngle();
 230     SetCutMaxDCADauther();
 231     SetCutMassGamma();
 232     SetCutMassGammaNSigma();
 233     SetCutMassK0();
 234     SetCutMassK0NSigma();
 235     SetCutChi2cGamma();
 236     SetCutChi2cK0();
 237     SetCutGammaSFromDecay();
 238     SetCutK0SFromDecay();
 239     SetCutMaxDCA();
 240   }
 241   //
 242   fTracklets                 = 0;
 243   fSClusters                 = 0;
 244   //
 245   // definition of histograms
 246   Bool_t oldStatus = TH1::AddDirectoryStatus();
 247   TH1::AddDirectory(kFALSE);
 248
 249   fhClustersDPhiAcc   = new TH1F("dphiacc",  "dphi",  100,-0.1,0.1);
 250   fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
 251
 252   fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,-0.1,0.1);
 253
 254   fhClustersDPhiAll   = new TH1F("dphiall",  "dphi",  100,0.0,0.5);
 255   fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,0.0,0.5);
 256
 257   fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,0.,0.5,100,0.,0.5);
 258
 259   fhetaTracklets  = new TH1F("etaTracklets",  "eta",  100,-2.,2.);
 260   fhphiTracklets  = new TH1F("phiTracklets",  "phi",  100, 0., 2*TMath::Pi());
 261   fhetaClustersLay1  = new TH1F("etaClustersLay1",  "etaCl1",  100,-2.,2.);
 262   fhphiClustersLay1  = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
 263   for (int i=2;i--;) fStoreRefs[i][0] =  fStoreRefs[i][1] = kFALSE;
 264   TH1::AddDirectory(oldStatus);
 265 }
 266
 267 //______________________________________________________________________
 268 AliITSMultReconstructor::AliITSMultReconstructor(const AliITSMultReconstructor &mr) :
 269 AliTrackleter(mr),
 270 fDetTypeRec(0),fESDEvent(0),fTreeRP(0),fTreeRPMix(0),
 271 fTracklets(0),
 272 fSClusters(0),
 273 fNTracklets(0),
 274 fNSingleCluster(0),
 275 fDPhiWindow(0),
 276 fDThetaWindow(0),
 277 fPhiShift(0),
 278 fRemoveClustersFromOverlaps(0),
 279 fPhiOverlapCut(0),
 280 fZetaOverlapCut(0),
 281 fPhiRotationAngle(0),
 282 fScaleDTBySin2T(0),
 283 fNStdDev(1.0),
 284 fNStdDevSq(1.0),
 285 //
 286 fCutPxDrSPDin(0.1),
 287 fCutPxDrSPDout(0.15),
 288 fCutPxDz(0.2),
 289 fCutDCArz(0.5),
 290 fCutMinElectronProbTPC(0.5),
 291 fCutMinElectronProbESD(0.1),
 292 fCutMinP(0.05),
 293 fCutMinRGamma(2.),
 294 fCutMinRK0(1.),
 295 fCutMinPointAngle(0.98),
 296 fCutMaxDCADauther(0.5),
 297 fCutMassGamma(0.03),
 298 fCutMassGammaNSigma(5.),
 299 fCutMassK0(0.03),
 300 fCutMassK0NSigma(5.),
 301 fCutChi2cGamma(2.),
 302 fCutChi2cK0(2.),
 303 fCutGammaSFromDecay(-10.),
 304 fCutK0SFromDecay(-10.),
 305 fCutMaxDCA(1.),
 306 //
 307 fHistOn(0),
 308 fhClustersDPhiAcc(0),
 309 fhClustersDThetaAcc(0),
 310 fhClustersDPhiAll(0),
 311 fhClustersDThetaAll(0),
 312 fhDPhiVsDThetaAll(0),
 313 fhDPhiVsDThetaAcc(0),
 314 fhetaTracklets(0),
 315 fhphiTracklets(0),
 316 fhetaClustersLay1(0),
 317 fhphiClustersLay1(0),
 318 fDPhiShift(0),
 319 fDPhiWindow2(0),
 320 fDThetaWindow2(0),
 321 fPartners(0),
 322 fAssociatedLay1(0),
 323 fMinDists(0),
 324 fBlackList(0),
 325 //
 326 fCreateClustersCopy(0),
 327 fClustersLoaded(0),
 328 fRecoDone(0),
 329 fBuildRefs(kTRUE),
 330 fSPDSeg()
 331  {
 332   // Copy constructor :!!! RS ATTENTION: old c-tor reassigned the pointers instead of creating a new copy -> would crash on delete
 333    AliError("May not use");
 334 }
 335
 336 //______________________________________________________________________
 337 AliITSMultReconstructor& AliITSMultReconstructor::operator=(const AliITSMultReconstructor& mr){
 338   // Assignment operator
 339   if (this != &mr) {
 340     this->~AliITSMultReconstructor();
 341     new(this) AliITSMultReconstructor(mr);
 342   }
 343   return *this;
 344 }
 345
 346 //______________________________________________________________________
 347 AliITSMultReconstructor::~AliITSMultReconstructor(){
 348   // Destructor
 349
 350   // delete histograms
 351   delete fhClustersDPhiAcc;
 352   delete fhClustersDThetaAcc;
 353   delete fhClustersDPhiAll;
 354   delete fhClustersDThetaAll;
 355   delete fhDPhiVsDThetaAll;
 356   delete fhDPhiVsDThetaAcc;
 357   delete fhetaTracklets;
 358   delete fhphiTracklets;
 359   delete fhetaClustersLay1;
 360   delete fhphiClustersLay1;
 361   //
 362   // delete arrays
 363   for(Int_t i=0; i<fNTracklets; i++) delete [] fTracklets[i];
 364
 365   for(Int_t i=0; i<fNSingleCluster; i++) delete [] fSClusters[i];
 366
 367   //
 368   for (int i=0;i<2;i++) {
 369     delete[] fClustersLay[i];
 370     delete[] fDetectorIndexClustersLay[i];
 371     delete[] fClusterCopyIndex[i];
 372     delete[] fOverlapFlagClustersLay[i];
 373     delete   fClArr[i];
 374     for (int j=0;j<2;j++) delete fUsedClusLay[i][j];
 375   }
 376   delete [] fTracklets;
 377   delete [] fSClusters;
 378   //
 379   delete[] fPartners;      fPartners = 0;
 380   delete[] fMinDists;      fMinDists = 0;
 381   delete   fBlackList;     fBlackList = 0;
 382   //
 383 }
 384
 385 //____________________________________________________________________
 386 void AliITSMultReconstructor::Reconstruct(AliESDEvent* esd, TTree* treeRP)
 387 {
 388   if (!treeRP) { AliError(" Invalid ITS cluster tree !\n"); return; }
 389   if (!esd) {AliError("ESDEvent is not available, use old reconstructor"); return;}
 390   // reset counters
 391   if (fMult) delete fMult; fMult = 0;
 392   fNClustersLay[0] = 0;
 393   fNClustersLay[1] = 0;
 394   fNTracklets = 0;
 395   fNSingleCluster = 0;
 396   //
 397   fESDEvent = esd;
 398   fTreeRP = treeRP;
 399   //
 400   // >>>> RS: this part is equivalent to former AliITSVertexer::FindMultiplicity
 401   //
 402   // see if there is a SPD vertex
 403   Bool_t isVtxOK=kTRUE, isCosmics=kFALSE;
 404   AliESDVertex* vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexSPD();
 405   if (!vtx || vtx->GetNContributors()<1) isVtxOK = kFALSE;
 406   if (vtx && strstr(vtx->GetTitle(),"cosmics")) {
 407     isVtxOK = kFALSE;
 408     isCosmics = kTRUE;
 409   }
 410   //
 411   if (!isVtxOK) {
 412     if (!isCosmics) {
 413       AliDebug(1,"Tracklets multiplicity not determined because the primary vertex was not found");
 414       AliDebug(1,"Just counting the number of cluster-fired chips on the SPD layers");
 415     }
 416     vtx = 0;
 417   }
 418   if(vtx){
 419     float vtxf[3] = {vtx->GetX(),vtx->GetY(),vtx->GetZ()};
 420     FindTracklets(vtxf);
 421   }
 422   else {
 423     FindTracklets(0);
 424   }
 425   //
 426   CreateMultiplicityObject();
 427 }
 428
 429 //____________________________________________________________________
 430 void AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) {
 431   //
 432   // RS NOTE - this is old reconstructor invocation, to be used from VertexFinder and in analysis mode
 433
 434   if (fMult) delete fMult; fMult = 0;
 435   fNClustersLay[0] = 0;
 436   fNClustersLay[1] = 0;
 437   fNTracklets = 0;
 438   fNSingleCluster = 0;
 439   //
 440   if (!clusterTree) { AliError(" Invalid ITS cluster tree !\n"); return; }
 441   //
 442   fESDEvent = 0;
 443   SetTreeRP(clusterTree);
 444   //
 445   FindTracklets(vtx);
 446   //
 447 }
 448
 449
 450 //____________________________________________________________________
 451 void AliITSMultReconstructor::ReconstructMix(TTree* clusterTree, TTree* clusterTreeMix, const Float_t* vtx, Float_t*)
 452 {
 453   //
 454   // RS NOTE - this is old reconstructor invocation, to be used from VertexFinder and in analysis mode
 455
 456   if (fMult) delete fMult; fMult = 0;
 457   fNClustersLay[0] = 0;
 458   fNClustersLay[1] = 0;
 459   fNTracklets = 0;
 460   fNSingleCluster = 0;
 461   //
 462   if (!clusterTree) { AliError(" Invalid ITS cluster tree !\n"); return; }
 463   if (!clusterTreeMix) { AliError(" Invalid ITS cluster tree 2nd event !\n"); return; }
 464   //
 465   fESDEvent = 0;
 466   SetTreeRP(clusterTree);
 467   SetTreeRPMix(clusterTreeMix);
 468   //
 469   FindTracklets(vtx);
 470   //
 471 }
 472
 473
 474 //____________________________________________________________________
 475 void AliITSMultReconstructor::FindTracklets(const Float_t *vtx)
 476 {
 477   // - calls LoadClusterArrays that finds the position of the clusters
 478   //   (in global coord)
 479
 480   // - convert the cluster coordinates to theta, phi (seen from the
 481   //   interaction vertex). Clusters in the inner layer can be now
 482   //   rotated for combinatorial studies
 483   // - makes an array of tracklets
 484   //
 485   // After this method has been called, the clusters of the two layers
 486   // and the tracklets can be retrieved by calling the Get'er methods.
 487
 488
 489   // Find tracklets converging to vertex
 490   //
 491   LoadClusterArrays(fTreeRP,fTreeRPMix);
 492   // flag clusters used by ESD tracks
 493   if (fESDEvent) ProcessESDTracks();
 494   fRecoDone = kTRUE;
 495
 496   if (!vtx) return;
 497
 498   InitAux();
 499
 500   // find the tracklets
 501   AliDebug(1,"Looking for tracklets... ");
 502
 503   ClusterPos2Angles(vtx); // convert cluster position to angles wrt vtx
 504   //
 505   // Step1: find all tracklets allowing double assocation:
 506   int found = 1;
 507   while (found > 0) {
 508     found = 0;
 509     for (Int_t iC1=0; iC1<fNClustersLay[0]; iC1++) found += AssociateClusterOfL1(iC1);
 510   }
 511   //
 512   // Step2: store tracklets; remove used clusters
 513   for (Int_t iC2=0; iC2<fNClustersLay[1]; iC2++) StoreTrackletForL2Cluster(iC2);
 514   //
 515   // store unused single clusters of L1
 516   StoreL1Singles();
 517   //
 518   AliDebug(1,Form("%d tracklets found", fNTracklets));
 519 }
 520
 521 //____________________________________________________________________
 522 void AliITSMultReconstructor::CreateMultiplicityObject()
 523 {
 524   // create AliMultiplicity object and store it in the ESD event
 525   //
 526   TBits fastOrFiredMap,firedChipMap;
 527   if (fDetTypeRec) {
 528    fastOrFiredMap  = fDetTypeRec->GetFastOrFiredMap();
 529    firedChipMap    = fDetTypeRec->GetFiredChipMap(fTreeRP);
 530   }
 531   //
 532   fMult = new AliMultiplicity(fNTracklets,fNSingleCluster,fNFiredChips[0],fNFiredChips[1],fastOrFiredMap);
 533   fMult->SetMultTrackRefs( fBuildRefs );
 534   // store some details of reco:
 535   fMult->SetScaleDThetaBySin2T(fScaleDTBySin2T);
 536   fMult->SetDPhiWindow2(fDPhiWindow2);
 537   fMult->SetDThetaWindow2(fDThetaWindow2);
 538   fMult->SetDPhiShift(fDPhiShift);
 539   fMult->SetNStdDev(fNStdDev);
 540   //
 541   fMult->SetFiredChipMap(firedChipMap);
 542   AliITSRecPointContainer* rcont = AliITSRecPointContainer::Instance();
 543   fMult->SetITSClusters(0,rcont->GetNClustersInLayer(1,fTreeRP));
 544   for(Int_t kk=2;kk<=6;kk++) fMult->SetITSClusters(kk-1,rcont->GetNClustersInLayerFast(kk));
 545   //
 546   UInt_t shared[100];
 547   AliRefArray *refs[2][2] = {{0,0},{0,0}};
 548   if (fBuildRefs) {
 549     for (int il=2;il--;)
 550       for (int it=2;it--;)  // tracklet_clusters->track references to stor
 551         if (fStoreRefs[il][it]) refs[il][it] = new AliRefArray(fNTracklets,0);
 552   }
 553   //
 554   for (int i=fNTracklets;i--;)  {
 555     float* tlInfo = fTracklets[i];
 556     fMult->SetTrackletData(i,tlInfo);
 557     //
 558     if (!fBuildRefs) continue; // do we need references?
 559     for (int itp=0;itp<2;itp++) {
 560       for (int ilr=0;ilr<2;ilr++) {
 561         if (!fStoreRefs[ilr][itp]) continue; // nothing to store
 562         int clID = int(tlInfo[ilr ? kClID2:kClID1]);
 563         int nref = fUsedClusLay[ilr][itp]->GetReferences(clID,shared,100);
 564         if (!nref) continue;
 565         else if (nref==1) refs[ilr][itp]->AddReference(i,shared[0]);
 566         else refs[ilr][itp]->AddReferences(i,shared,nref);
 567       }
 568     }
 569   }
 570   if (fBuildRefs) fMult->AttachTracklet2TrackRefs(refs[0][0],refs[0][1],refs[1][0],refs[1][1]);
 571   //
 572   AliRefArray *refsc[2] = {0,0};
 573   if (fBuildRefs) for (int it=2;it--;) if (fStoreRefs[0][it]) refsc[it] = new AliRefArray(fNClustersLay[0]);
 574   for (int i=fNSingleCluster;i--;) {
 575     float* clInfo = fSClusters[i];
 576     fMult->SetSingleClusterData(i,clInfo);
 577     //
 578     if (!fBuildRefs) continue; // do we need references?
 579     int clID = int(clInfo[kSCID]);
 580     for (int itp=0;itp<2;itp++) {
 581       if (!fStoreRefs[0][itp]) continue;
 582       int nref = fUsedClusLay[0][itp]->GetReferences(clID,shared,100);
 583       if (!nref) continue;
 584       else if (nref==1) refsc[itp]->AddReference(i,shared[0]);
 585       else refsc[itp]->AddReferences(i,shared,nref);
 586     }
 587   }
 588   if (fBuildRefs) fMult->AttachCluster2TrackRefs(refsc[0],refsc[1]);
 589   fMult->CompactBits();
 590   //
 591 }
 592
 593
 594 //____________________________________________________________________
 595 void AliITSMultReconstructor::LoadClusterArrays(TTree* tree, TTree* treeMix)
 596 {
 597   // load cluster info and prepare tracklets arrays
 598   //
 599   if (AreClustersLoaded()) {AliInfo("Clusters are already loaded"); return;}
 600   LoadClusterArrays(tree,0);
 601   LoadClusterArrays(treeMix ? treeMix:tree,1);
 602   int nmaxT = TMath::Min(fNClustersLay[0], fNClustersLay[1]);
 603   if (fTracklets) delete[] fTracklets;
 604   fTracklets = new Float_t*[nmaxT];
 605   memset(fTracklets,0,nmaxT*sizeof(Float_t*));
 606   //
 607   if (fSClusters) delete[] fSClusters;
 608   fSClusters = new Float_t*[fNClustersLay[0]];
 609   memset(fSClusters,0,fNClustersLay[0]*sizeof(Float_t*));
 610   //
 611   AliDebug(1,Form("(clusters in layer 1 : %d,  layer 2: %d)",fNClustersLay[0],fNClustersLay[1]));
 612   AliDebug(1,Form("(cluster-fired chips in layer 1 : %d,  layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
 613   SetClustersLoaded();
 614 }
 615
 616 //____________________________________________________________________
 617 void AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree, int il)
 618 {
 619   // This method
 620   // - gets the clusters from the cluster tree for layer il
 621   // - convert them into global coordinates
 622   // - store them in the internal arrays
 623   // - count the number of cluster-fired chips
 624   //
 625   // RS: This method was strongly modified wrt original. In order to have the same numbering
 626   // of clusters as in the ITS reco I had to introduce sorting in Z
 627   // Also note that now the clusters data are stored not in float[6] attached to float**, but in 1-D array
 628   AliDebug(1,Form("Loading clusters and cluster-fired chips for layer %d",il));
 629   //
 630   fNClustersLay[il] = 0;
 631   fNFiredChips[il]  = 0;
 632   for (int i=2;i--;) fStoreRefs[il][i] = kFALSE;
 633   //
 634   AliITSRecPointContainer* rpcont = 0;
 635   static TClonesArray statITSrec("AliITSRecPoint");
 636   static TObjArray clArr(100);
 637   TBranch* branch = 0;
 638   TClonesArray* itsClusters = 0;
 639   //
 640   if (!fCreateClustersCopy) {
 641     rpcont=AliITSRecPointContainer::Instance();
 642     itsClusters = rpcont->FetchClusters(0,itsClusterTree);
 643     if(!rpcont->IsSPDActive()){
 644       AliWarning("No SPD rec points found, multiplicity not calculated");
 645       return;
 646     }
 647   }
 648   else {
 649     itsClusters = &statITSrec;
 650     branch = itsClusterTree->GetBranch("ITSRecPoints");
 651     branch->SetAddress(&itsClusters);
 652     if (!fClArr[il]) fClArr[il] = new TClonesArray("AliITSRecPoint",100);
 653     delete[] fClusterCopyIndex[il];
 654   }
 655   //
 656   // count clusters
 657   // loop over the SPD subdetectors
 658   int nclLayer = 0;
 659   int detMin = TMath::Max(0,AliITSgeomTGeo::GetModuleIndex(il+1,1,1));
 660   int detMax = AliITSgeomTGeo::GetModuleIndex(il+2,1,1);
 661   for (int idt=detMin;idt<detMax;idt++) {
 662     if (!fCreateClustersCopy) itsClusters = rpcont->UncheckedGetClusters(idt);
 663     else                      branch->GetEvent(idt);
 664     int nClusters = itsClusters->GetEntriesFast();
 665     if (!nClusters) continue;
 666     Int_t nClustersInChip[5] = {0,0,0,0,0};
 667     while(nClusters--) {
 668       AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
 669       if (!cluster) continue;
 670       if (fCreateClustersCopy)  cluster = new ((*fClArr[il])[nclLayer]) AliITSRecPoint(*cluster);
 671       clArr.AddAtAndExpand(cluster,nclLayer++);
 672       Int_t chipNo = fSPDSeg.GetChipFromLocal(0,cluster->GetDetLocalZ());
 673       if(chipNo>=0)nClustersInChip[ chipNo ]++;
 674     }
 675     for(Int_t ifChip=5;ifChip--;) if (nClustersInChip[ifChip]) fNFiredChips[il]++;
 676   }
 677   // sort the clusters in Z (to have the same numbering as in ITS reco
 678   Float_t *z     = new Float_t[nclLayer];
 679   Int_t   *index = new Int_t[nclLayer];
 680   for (int ic=0;ic<nclLayer;ic++) z[ic] = ((AliITSRecPoint*)clArr[ic])->GetZ();
 681   TMath::Sort(nclLayer,z,index,kFALSE);
 682   Float_t*   clustersLay              = new Float_t[nclLayer*kClNPar];
 683   Int_t*     detectorIndexClustersLay = new Int_t[nclLayer];
 684   Bool_t*    overlapFlagClustersLay   = new Bool_t[nclLayer];
 685   if (fCreateClustersCopy) fClusterCopyIndex[il] = new Int_t[nclLayer];
 686   //
 687   for (int ic=0;ic<nclLayer;ic++) {
 688     AliITSRecPoint* cluster = (AliITSRecPoint*)clArr[index[ic]];
 689     float* clPar = &clustersLay[ic*kClNPar];
 690     //
 691     cluster->GetGlobalXYZ( clPar );
 692     detectorIndexClustersLay[ic] = cluster->GetDetectorIndex();
 693     overlapFlagClustersLay[ic]   = kFALSE;
 694     for (Int_t i=3;i--;) clPar[kClMC0+i] = cluster->GetLabel(i);
 695     if (fCreateClustersCopy) fClusterCopyIndex[il][ic] = index[ic];
 696   }
 697   clArr.Clear();
 698   delete[] z;
 699   delete[] index;
 700   //
 701   if (fOverlapFlagClustersLay[il]) delete[] fOverlapFlagClustersLay[il];
 702   fOverlapFlagClustersLay[il]   = overlapFlagClustersLay;
 703   //
 704   if (fDetectorIndexClustersLay[il]) delete[] fDetectorIndexClustersLay[il];
 705   fDetectorIndexClustersLay[il] = detectorIndexClustersLay;
 706   //
 707   if (fBuildRefs) {
 708     for (int it=0;it<2;it++) {
 709       if (fUsedClusLay[il][it]) delete fUsedClusLay[il][it];
 710       fUsedClusLay[il][it] = new AliRefArray(nclLayer);
 711     }
 712   }
 713   //
 714   if (fClustersLay[il]) delete[] fClustersLay[il];
 715   fClustersLay[il] = clustersLay;
 716   fNClustersLay[il] = nclLayer;
 717   //
 718 }
 719
 720 //____________________________________________________________________
 721 void AliITSMultReconstructor::LoadClusterFiredChips(TTree* itsClusterTree) {
 722   // This method
 723   // - gets the clusters from the cluster tree
 724   // - counts the number of (cluster)fired chips
 725
 726   AliDebug(1,"Loading cluster-fired chips ...");
 727
 728   fNFiredChips[0] = 0;
 729   fNFiredChips[1] = 0;
 730
 731   AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
 732   TClonesArray* itsClusters=NULL;
 733   rpcont->FetchClusters(0,itsClusterTree);
 734   if(!rpcont->IsSPDActive()){
 735     AliWarning("No SPD rec points found, multiplicity not calculated");
 736     return;
 737   }
 738
 739   // loop over the its subdetectors
 740   Int_t nSPDmodules=AliITSgeomTGeo::GetModuleIndex(3,1,1);
 741   for (Int_t iIts=0; iIts < nSPDmodules; iIts++) {
 742     itsClusters=rpcont->UncheckedGetClusters(iIts);
 743     Int_t nClusters = itsClusters->GetEntriesFast();
 744
 745     // number of clusters in each chip of the current module
 746     Int_t nClustersInChip[5] = {0,0,0,0,0};
 747     Int_t layer = 0;
 748     Int_t ladder=0;
 749     Int_t det=0;
 750     AliITSgeomTGeo::GetModuleId(iIts,layer,ladder,det);
 751     --layer;  // layer is from 1 to 6 in AliITSgeomTGeo, but from 0 to 5 here
 752     if(layer<0 || layer >1)continue;
 753
 754     // loop over clusters
 755     while(nClusters--) {
 756       AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
 757
 758       // find the chip for the current cluster
 759       Float_t locz = cluster->GetDetLocalZ();
 760       Int_t iChip = fSPDSeg.GetChipFromLocal(0,locz);
 761       if (iChip>=0) nClustersInChip[iChip]++;
 762
 763     }// end of cluster loop
 764
 765     // get number of fired chips in the current module
 766     for(Int_t ifChip=0; ifChip<5; ifChip++) {
 767       if(nClustersInChip[ifChip] >= 1)  fNFiredChips[layer]++;
 768     }
 769
 770   } // end of its "subdetector" loop
 771
 772
 773   AliDebug(1,Form("(cluster-fired chips in layer 1 : %d,  layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
 774 }
 775 //____________________________________________________________________
 776 void
 777 AliITSMultReconstructor::SaveHists() {
 778   // This method save the histograms on the output file
 779   // (only if fHistOn is TRUE).
 780
 781   if (!fHistOn)
 782     return;
 783
 784   fhClustersDPhiAll->Write();
 785   fhClustersDThetaAll->Write();
 786   fhDPhiVsDThetaAll->Write();
 787
 788   fhClustersDPhiAcc->Write();
 789   fhClustersDThetaAcc->Write();
 790   fhDPhiVsDThetaAcc->Write();
 791
 792   fhetaTracklets->Write();
 793   fhphiTracklets->Write();
 794   fhetaClustersLay1->Write();
 795   fhphiClustersLay1->Write();
 796 }
 797
 798 //____________________________________________________________________
 799 void AliITSMultReconstructor::FlagClustersInOverlapRegions (Int_t iC1, Int_t iC2WithBestDist)
 800 {
 801   // Flags clusters in the overlapping regions
 802   Float_t distClSameMod=0.;
 803   Float_t distClSameModMin=0.;
 804   Int_t   iClOverlap =0;
 805   Float_t meanRadiusLay1 = 3.99335; // average radius inner layer
 806   Float_t meanRadiusLay2 = 7.37935; // average radius outer layer;
 807
 808   Float_t zproj1=0.;
 809   Float_t zproj2=0.;
 810   Float_t deZproj=0.;
 811   Float_t* clPar1  = GetClusterLayer1(iC1);
 812   Float_t* clPar2B = GetClusterLayer2(iC2WithBestDist);
 813   // Loop on inner layer clusters
 814   for (Int_t iiC1=0; iiC1<fNClustersLay[0]; iiC1++) {
 815     if (!fOverlapFlagClustersLay[0][iiC1]) {
 816       // only for adjacent modules
 817       if ((TMath::Abs(fDetectorIndexClustersLay[0][iC1]-fDetectorIndexClustersLay[0][iiC1])==4)||
 818          (TMath::Abs(fDetectorIndexClustersLay[0][iC1]-fDetectorIndexClustersLay[0][iiC1])==76)) {
 819         Float_t *clPar11 = GetClusterLayer1(iiC1);
 820         Float_t dePhi=TMath::Abs(clPar11[kClPh]-clPar1[kClPh]);
 821         if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
 822
 823         zproj1=meanRadiusLay1/TMath::Tan(clPar1[kClTh]);
 824         zproj2=meanRadiusLay1/TMath::Tan(clPar11[kClTh]);
 825
 826         deZproj=TMath::Abs(zproj1-zproj2);
 827
 828         distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
 829         if (distClSameMod<=1.) fOverlapFlagClustersLay[0][iiC1]=kTRUE;
 830
 831 //        if (distClSameMod<=1.) {
 832 //          if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
 833 //            distClSameModMin=distClSameMod;
 834 //            iClOverlap=iiC1;
 835 //          }
 836 //        }
 837
 838
 839       } // end adjacent modules
 840     }
 841   } // end Loop on inner layer clusters
 842
 843 //  if (distClSameModMin!=0.) fOverlapFlagClustersLay[0][iClOverlap]=kTRUE;
 844
 845   distClSameMod=0.;
 846   distClSameModMin=0.;
 847   iClOverlap =0;
 848   // Loop on outer layer clusters
 849   for (Int_t iiC2=0; iiC2<fNClustersLay[1]; iiC2++) {
 850     if (!fOverlapFlagClustersLay[1][iiC2]) {
 851       // only for adjacent modules
 852       Float_t *clPar2 = GetClusterLayer2(iiC2);
 853       if ((TMath::Abs(fDetectorIndexClustersLay[1][iC2WithBestDist]-fDetectorIndexClustersLay[1][iiC2])==4) ||
 854          (TMath::Abs(fDetectorIndexClustersLay[1][iC2WithBestDist]-fDetectorIndexClustersLay[1][iiC2])==156)) {
 855         Float_t dePhi=TMath::Abs(clPar2[kClPh]-clPar2B[kClPh]);
 856         if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
 857
 858         zproj1=meanRadiusLay2/TMath::Tan(clPar2B[kClTh]);
 859         zproj2=meanRadiusLay2/TMath::Tan(clPar2[kClTh]);
 860
 861         deZproj=TMath::Abs(zproj1-zproj2);
 862         distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
 863         if (distClSameMod<=1.) fOverlapFlagClustersLay[1][iiC2]=kTRUE;
 864
 865 //        if (distClSameMod<=1.) {
 866 //          if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
 867 //            distClSameModMin=distClSameMod;
 868 //            iClOverlap=iiC2;
 869 //          }
 870 //        }
 871
 872       } // end adjacent modules
 873     }
 874   } // end Loop on outer layer clusters
 875
 876 //  if (distClSameModMin!=0.) fOverlapFlagClustersLay[1][iClOverlap]=kTRUE;
 877
 878 }
 879
 880 //____________________________________________________________________
 881 void AliITSMultReconstructor::InitAux()
 882 {
 883   // init arrays/parameters for tracklet reconstruction
 884
 885   // dPhi shift is field dependent, get average magnetic field
 886   Float_t bz = 0;
 887   AliMagF* field = 0;
 888   if (TGeoGlobalMagField::Instance()) field = dynamic_cast<AliMagF*>(TGeoGlobalMagField::Instance()->GetField());
 889   if (!field) {
 890     AliError("Could not retrieve magnetic field. Assuming no field. Delta Phi shift will be deactivated in AliITSMultReconstructor.");
 891   }
 892   else bz = TMath::Abs(field->SolenoidField());
 893   fDPhiShift = fPhiShift / 5 * bz;
 894   AliDebug(1, Form("Using phi shift of %f", fDPhiShift));
 895   //
 896   if (fPartners) delete[] fPartners; fPartners = new Int_t[fNClustersLay[1]];
 897   if (fMinDists) delete[] fMinDists; fMinDists = new Float_t[fNClustersLay[1]];
 898   if (fAssociatedLay1) delete[] fAssociatedLay1; fAssociatedLay1 = new Int_t[fNClustersLay[0]];
 899   //
 900   if (fBlackList) delete fBlackList; fBlackList = new AliRefArray(fNClustersLay[0]);
 901   //
 902   //  Printf("Vertex in find tracklets...%f %f %f",vtx[0],vtx[1],vtx[2]);
 903   for (Int_t i=0; i<fNClustersLay[1]; i++) {
 904     fPartners[i] = -1;
 905     fMinDists[i] = 2*fNStdDev;
 906   }
 907   memset(fAssociatedLay1,0,fNClustersLay[0]*sizeof(Int_t));
 908   //
 909 }
 910
 911 //____________________________________________________________________
 912 void AliITSMultReconstructor::ClusterPos2Angles(const Float_t *vtx)
 913 {
 914   // convert cluster coordinates to angles wrt vertex
 915   for (int ilr=0;ilr<2;ilr++) {
 916     for (Int_t iC=0; iC<fNClustersLay[ilr]; iC++) {
 917       float* clPar = GetClusterOfLayer(ilr,iC);
 918       CalcThetaPhi(clPar[kClTh]-vtx[0],clPar[kClPh]-vtx[1],clPar[kClZ]-vtx[2],clPar[kClTh],clPar[kClPh]);
 919       if (ilr==0) {
 920         clPar[kClPh] = clPar[kClPh] + fPhiRotationAngle;   // rotation of inner layer for comb studies
 921         if (fHistOn) {
 922           Float_t eta = clPar[kClTh];
 923           eta= TMath::Tan(eta/2.);
 924           eta=-TMath::Log(eta);
 925           fhetaClustersLay1->Fill(eta);
 926           fhphiClustersLay1->Fill(clPar[kClPh]);
 927         }
 928       }
 929     }
 930   }
 931   //
 932 }
 933
 934 //____________________________________________________________________
 935 Int_t AliITSMultReconstructor::AssociateClusterOfL1(Int_t iC1)
 936 {
 937   // search association of cluster iC1 of L1 with all clusters of L2
 938   if (fAssociatedLay1[iC1] != 0) return 0;
 939   Int_t  iC2WithBestDist = -1;   // reset
 940   Double_t minDist       =  2*fNStdDev;   // reset
 941   float* clPar1 = GetClusterLayer1(iC1);
 942   for (Int_t iC2=0; iC2<fNClustersLay[1]; iC2++) {
 943     //
 944     if (fBlackList->IsReferred(iC1,iC2)) continue;
 945     float* clPar2 = GetClusterLayer2(iC2);
 946     //
 947     // find the difference in angles
 948     Double_t dTheta = TMath::Abs(clPar2[kClTh] - clPar1[kClTh]);
 949     Double_t dPhi   = TMath::Abs(clPar2[kClPh] - clPar1[kClPh]);
 950     //        Printf("detheta %f  dephi %f", dTheta,dPhi);
 951     //
 952     if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;     // take into account boundary condition
 953     //
 954     if (fHistOn) {
 955       fhClustersDPhiAll->Fill(dPhi);
 956       fhClustersDThetaAll->Fill(dTheta);
 957       fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
 958     }
 959     Float_t d = CalcDist(dPhi,dTheta,clPar1[kClTh]);     // make "elliptical" cut in Phi and Theta!
 960     // look for the minimum distance: the minimum is in iC2WithBestDist
 961     if (d<fNStdDev && d<minDist) { minDist=d; iC2WithBestDist = iC2; }
 962   }
 963   //
 964   if (minDist<fNStdDev) { // This means that a cluster in layer 2 was found that matches with iC1
 965     //
 966     if (fMinDists[iC2WithBestDist] > minDist) {
 967       Int_t oldPartner = fPartners[iC2WithBestDist];
 968       fPartners[iC2WithBestDist] = iC1;
 969       fMinDists[iC2WithBestDist] = minDist;
 970       //
 971       fAssociatedLay1[iC1] = 1;      // mark as assigned
 972       //
 973       if (oldPartner != -1) {
 974         // redo partner search for cluster in L0 (oldPartner), putting this one (iC2WithBestDist) on its fBlackList
 975         fBlackList->AddReference(oldPartner,iC2WithBestDist);
 976         fAssociatedLay1[oldPartner] = 0;       // mark as free
 977       }
 978     } else {
 979       // try again to find a cluster without considering iC2WithBestDist
 980       fBlackList->AddReference(iC1,iC2WithBestDist);
 981     }
 982     //
 983   }
 984   else fAssociatedLay1[iC1] = 2;// cluster has no partner; remove
 985   //
 986   return 1;
 987 }
 988
 989 //____________________________________________________________________
 990 Int_t AliITSMultReconstructor::StoreTrackletForL2Cluster(Int_t iC2)
 991 {
 992   // build tracklet for cluster iC2 of layer 2
 993   if (fPartners[iC2] == -1) return 0;
 994   if (fRemoveClustersFromOverlaps) FlagClustersInOverlapRegions (fPartners[iC2],iC2);
 995   // Printf("saving tracklets");
 996   if (fOverlapFlagClustersLay[0][fPartners[iC2]] || fOverlapFlagClustersLay[1][iC2]) return 0;
 997   float* clPar2 = GetClusterLayer2(iC2);
 998   float* clPar1 = GetClusterLayer1(fPartners[iC2]);
 999   //
1000   Float_t* tracklet = fTracklets[fNTracklets] = new Float_t[kTrNPar]; // RS Add also the cluster id's
1001   //
1002   tracklet[kTrTheta] = clPar1[kClTh];    // use the theta from the clusters in the first layer
1003   tracklet[kTrPhi]   = clPar1[kClPh];    // use the phi from the clusters in the first layer
1004   tracklet[kTrDPhi] = clPar1[kClPh] - clPar2[kClPh];  // store the difference between phi1 and phi2
1005   //
1006   // define dphi in the range [0,pi] with proper sign (track charge correlated)
1007   if (tracklet[kTrDPhi] > TMath::Pi())   tracklet[kTrDPhi] = tracklet[kTrDPhi]-2.*TMath::Pi();
1008   if (tracklet[kTrDPhi] < -TMath::Pi())  tracklet[kTrDPhi] = tracklet[kTrDPhi]+2.*TMath::Pi();
1009   //
1010   tracklet[kTrDTheta] = clPar1[kClTh] - clPar2[kClTh]; // store the theta1-theta2
1011   //
1012   if (fHistOn) {
1013     fhClustersDPhiAcc->Fill(tracklet[kTrDPhi]);
1014     fhClustersDThetaAcc->Fill(tracklet[kTrDTheta]);
1015     fhDPhiVsDThetaAcc->Fill(tracklet[kTrDTheta],tracklet[kTrDPhi]);
1016   }
1017   //
1018   // find label
1019   // if equal label in both clusters found this label is assigned
1020   // if no equal label can be found the first labels of the L1 AND L2 cluster are assigned
1021   Int_t label1=0,label2=0;
1022   while (label2 < 3) {
1023     if ( int(clPar1[kClMC0+label1])!=-2 && int(clPar1[kClMC0+label1])==int(clPar2[kClMC0+label2])) break;
1024     if (++label1 == 3) { label1 = 0; label2++; }
1025   }
1026   if (label2 < 3) {
1027     AliDebug(AliLog::kDebug, Form("Found label %d == %d for tracklet candidate %d\n",
1028                                   (Int_t) clPar1[kClMC0+label1], (Int_t) clPar1[kClMC0+label2], fNTracklets));
1029     tracklet[kTrLab1] = tracklet[kTrLab2] = clPar1[kClMC0+label1];
1030   } else {
1031     AliDebug(AliLog::kDebug, Form("Did not find label %d %d %d %d %d %d for tracklet candidate %d\n",
1032                                   (Int_t) clPar1[kClMC0], (Int_t) clPar1[kClMC1], (Int_t) clPar1[kClMC2],
1033                                   (Int_t) clPar2[kClMC0], (Int_t) clPar2[kClMC1], (Int_t) clPar2[kClMC2], fNTracklets));
1034     tracklet[kTrLab1] = clPar1[kClMC0];
1035     tracklet[kTrLab2] = clPar2[kClMC0];
1036   }
1037   //
1038   if (fHistOn) {
1039     Float_t eta = tracklet[kTrTheta];
1040     eta= TMath::Tan(eta/2.);
1041     eta=-TMath::Log(eta);
1042     fhetaTracklets->Fill(eta);
1043     fhphiTracklets->Fill(tracklet[kTrPhi]);
1044   }
1045   //
1046   tracklet[kClID1] = fPartners[iC2];
1047   tracklet[kClID2] = iC2;
1048   //
1049   // Printf("Adding tracklet candidate");
1050   AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
1051   AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", fPartners[iC2], iC2));
1052   fNTracklets++;
1053   fAssociatedLay1[fPartners[iC2]] = 1;
1054   //
1055   return 1;
1056 }
1057
1058 //____________________________________________________________________
1059 void AliITSMultReconstructor::StoreL1Singles()
1060 {
1061   // Printf("saving single clusters...");
1062   for (Int_t iC1=0; iC1<fNClustersLay[0]; iC1++) {
1063     float* clPar1 = GetClusterLayer1(iC1);
1064     if (fAssociatedLay1[iC1]==2||fAssociatedLay1[iC1]==0) {
1065       fSClusters[fNSingleCluster] = new Float_t[kClNPar];
1066       fSClusters[fNSingleCluster][kSCTh] = clPar1[kClTh];
1067       fSClusters[fNSingleCluster][kSCPh] = clPar1[kClPh];
1068       fSClusters[fNSingleCluster][kSCLab] = clPar1[kClMC0];
1069       fSClusters[fNSingleCluster][kSCID] = iC1;
1070       AliDebug(1,Form(" Adding a single cluster %d (cluster %d  of layer 1)",
1071                       fNSingleCluster, iC1));
1072       fNSingleCluster++;
1073     }
1074   }
1075   //
1076 }
1077
1078 //____________________________________________________________________
1079 void AliITSMultReconstructor::ProcessESDTracks()
1080 {
1081   // Flag the clusters used by ESD tracks
1082   // Flag primary tracks to be used for multiplicity counting
1083   //
1084   if (!fESDEvent || !fBuildRefs) return;
1085   AliESDVertex* vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexTracks();
1086   if (!vtx || vtx->GetNContributors()<1) vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexSPD();
1087   if (!vtx || vtx->GetNContributors()<1) {
1088     AliDebug(1,"No primary vertex: cannot flag primary tracks");
1089     return;
1090   }
1091   Int_t ntracks = fESDEvent->GetNumberOfTracks();
1092   for(Int_t itr=0; itr<ntracks; itr++) {
1093     AliESDtrack* track = fESDEvent->GetTrack(itr);
1094     if (!track->IsOn(AliESDtrack::kITSin)) continue; // use only tracks propagated in ITS to vtx
1095     FlagTrackClusters(itr);
1096     FlagIfSecondary(track,vtx);
1097   }
1098   FlagV0s(vtx);
1099   //
1100 }
1101
1102 //____________________________________________________________________
1103 void AliITSMultReconstructor::FlagTrackClusters(Int_t id)
1104 {
1105   // RS: flag the SPD clusters of the track if it is useful for the multiplicity estimation
1106   //
1107   const AliESDtrack* track = fESDEvent->GetTrack(id);
1108   Int_t idx[12];
1109   if ( track->GetITSclusters(idx)<3 ) return; // at least 3 clusters must be used in the fit
1110   Int_t itsType = track->IsOn(AliESDtrack::kITSpureSA) ? 1:0;
1111
1112   for (int i=6/*AliESDfriendTrack::kMaxITScluster*/;i--;) { // ignore extras: note: i>=6 is for extra clusters
1113     if (idx[i]<0) continue;
1114     int layID= (idx[i] & 0xf0000000) >> 28;
1115     if (layID>1) continue; // SPD only
1116     int clID = (idx[i] & 0x0fffffff);
1117     fUsedClusLay[layID][itsType]->AddReference(clID,id);
1118     fStoreRefs[layID][itsType] = kTRUE;
1119   }
1120   //
1121 }
1122
1123 //____________________________________________________________________
1124 void AliITSMultReconstructor::FlagIfSecondary(AliESDtrack* track, const AliVertex* vtx)
1125 {
1126   // RS: check if the track is primary and set the flag
1127   double cut = (track->HasPointOnITSLayer(0)||track->HasPointOnITSLayer(1)) ? fCutPxDrSPDin:fCutPxDrSPDout;
1128   float xz[2];
1129   track->GetDZ(vtx->GetX(),vtx->GetY(),vtx->GetZ(), fESDEvent->GetMagneticField(), xz);
1130   if (TMath::Abs(xz[0]*track->P())>cut || TMath::Abs(xz[1]*track->P())>fCutPxDz ||
1131       TMath::Abs(xz[0])>fCutDCArz   || TMath::Abs(xz[1])>fCutDCArz)
1132     track->SetStatus(AliESDtrack::kMultSec);
1133   else track->ResetStatus(AliESDtrack::kMultSec);
1134 }
1135
1136 //____________________________________________________________________
1137 void AliITSMultReconstructor::FlagV0s(const AliESDVertex *vtx)
1138 {
1139   // flag tracks belonging to v0s
1140   //
1141   const double kK0Mass = 0.4976;
1142   //
1143   AliV0 pvertex;
1144   AliKFVertex vertexKF;
1145   AliKFParticle epKF0,epKF1,pipmKF0,piKF0,piKF1,gammaKF,k0KF;
1146   Double_t mass,massErr,chi2c;
1147   enum {kKFIni=BIT(14)};
1148   //
1149   double recVtx[3];
1150   float recVtxF[3];
1151   vtx->GetXYZ(recVtx);
1152   for (int i=3;i--;) recVtxF[i] = recVtx[i];
1153   //
1154   int ntracks = fESDEvent->GetNumberOfTracks();
1155   if (ntracks<2) return;
1156   //
1157   vertexKF.X() = recVtx[0];
1158   vertexKF.Y() = recVtx[1];
1159   vertexKF.Z() = recVtx[2];
1160   vertexKF.Covariance(0,0) = vtx->GetXRes()*vtx->GetXRes();
1161   vertexKF.Covariance(1,2) = vtx->GetYRes()*vtx->GetYRes();
1162   vertexKF.Covariance(2,2) = vtx->GetZRes()*vtx->GetZRes();
1163   //
1164   AliESDtrack *trc0,*trc1;
1165   for (int it0=0;it0<ntracks;it0++) {
1166     trc0 = fESDEvent->GetTrack(it0);
1167     if (trc0->IsOn(AliESDtrack::kMultInV0)) continue;
1168     if (!trc0->IsOn(AliESDtrack::kITSin)) continue;
1169     Bool_t isSAP = trc0->IsPureITSStandalone();
1170     Int_t  q0 = trc0->Charge();
1171     Bool_t testGamma = CanBeElectron(trc0);
1172     epKF0.ResetBit(kKFIni);
1173     piKF0.ResetBit(kKFIni);
1174     double bestChi2=1e16;
1175     int bestID = -1;
1176     //
1177     for (int it1=it0+1;it1<ntracks;it1++) {
1178       trc1 = fESDEvent->GetTrack(it1);
1179       if (trc1->IsOn(AliESDtrack::kMultInV0)) continue;
1180       if (!trc1->IsOn(AliESDtrack::kITSin)) continue;
1181       if (trc1->IsPureITSStandalone() != isSAP) continue; // pair separately ITS_SA_Pure tracks and TPC/ITS+ITS_SA
1182       if ( (q0+trc1->Charge())!=0 ) continue;             // don't pair like signs
1183       //
1184       pvertex.SetParamN(q0<0 ? *trc0:*trc1);
1185       pvertex.SetParamP(q0>0 ? *trc0:*trc1);
1186       pvertex.Update(recVtxF);
1187       if (pvertex.P()<fCutMinP) continue;
1188       if (pvertex.GetV0CosineOfPointingAngle()<fCutMinPointAngle) continue;
1189       if (pvertex.GetDcaV0Daughters()>fCutMaxDCADauther) continue;
1190       double d = pvertex.GetD(recVtx[0],recVtx[1],recVtx[2]);
1191       if (d>fCutMaxDCA) continue;
1192       double dx=recVtx[0]-pvertex.Xv(), dy=recVtx[1]-pvertex.Yv();
1193       double rv = TMath::Sqrt(dx*dx+dy*dy);
1194       //
1195       // check gamma conversion hypothesis ----------------------------------------------------------->>>
1196       Bool_t gammaOK = kFALSE;
1197       while (testGamma && CanBeElectron(trc1)) {
1198         if (rv<fCutMinRGamma) break;
1199         if (!epKF0.TestBit(kKFIni)) {
1200           new(&epKF0) AliKFParticle(*trc0,q0>0 ? kPositron:kElectron);
1201           epKF0.SetBit(kKFIni);
1202         }
1203         new(&epKF1) AliKFParticle(*trc1,q0<0 ? kPositron:kElectron);
1204         gammaKF.Initialize();
1205         gammaKF += epKF0;
1206         gammaKF += epKF1;
1207         gammaKF.SetProductionVertex(vertexKF);
1208         gammaKF.GetMass(mass,massErr);
1209         if (mass>fCutMassGamma || (massErr>0&&(mass>massErr*fCutMassGammaNSigma))) break;
1210         if (gammaKF.GetS()<fCutGammaSFromDecay) break;
1211         gammaKF.SetMassConstraint(0.,0.001);
1212         chi2c = (gammaKF.GetNDF()!=0) ? gammaKF.GetChi2()/gammaKF.GetNDF() : 1000;
1213         if (chi2c>fCutChi2cGamma) break;
1214         gammaOK = kTRUE;
1215         if (chi2c>bestChi2) break;
1216         bestChi2 = chi2c;
1217         bestID = it1;
1218         break;
1219       }
1220       if (gammaOK) continue;
1221       // check gamma conversion hypothesis -----------------------------------------------------------<<<
1222       // check K0 conversion hypothesis    ----------------------------------------------------------->>>
1223       while (1) {
1224         if (rv<fCutMinRK0) break;
1225         if (!piKF0.TestBit(kKFIni)) {
1226           new(&piKF0) AliKFParticle(*trc0,q0>0 ? kPiPlus:kPiMinus);
1227           piKF0.SetBit(kKFIni);
1228         }
1229         new(&piKF1) AliKFParticle(*trc1,q0<0 ? kPiPlus:kPiMinus);
1230         k0KF.Initialize();
1231         k0KF += piKF0;
1232         k0KF += piKF1;
1233         k0KF.SetProductionVertex(vertexKF);
1234         k0KF.GetMass(mass,massErr);
1235         mass -= kK0Mass;
1236         if (TMath::Abs(mass)>fCutMassK0 || (massErr>0&&(abs(mass)>massErr*fCutMassK0NSigma))) break;
1237         if (k0KF.GetS()<fCutK0SFromDecay) break;
1238         k0KF.SetMassConstraint(kK0Mass,0.001);
1239         chi2c = (k0KF.GetNDF()!=0) ? k0KF.GetChi2()/k0KF.GetNDF() : 1000;
1240         if (chi2c>fCutChi2cK0) break;
1241         if (chi2c>bestChi2) break;
1242         bestChi2 = chi2c;
1243         bestID = it1;
1244         break;
1245       }
1246       // check K0 conversion hypothesis    -----------------------------------------------------------<<<
1247     }
1248     //
1249     if (bestID>=0) {
1250       trc0->SetStatus(AliESDtrack::kMultInV0);
1251       fESDEvent->GetTrack(bestID)->SetStatus(AliESDtrack::kMultInV0);
1252     }
1253   }
1254   //
1255 }
1256
1257 //____________________________________________________________________
1258 Bool_t AliITSMultReconstructor::CanBeElectron(const AliESDtrack* trc) const
1259 {
1260   // check if the track can be electron
1261   Double_t pid[AliPID::kSPECIES];
1262   if (!trc->IsOn(AliESDtrack::kESDpid)) return kTRUE;
1263   trc->GetESDpid(pid);
1264   return (trc->IsOn(AliESDtrack::kTPCpid)) ?
1265     pid[AliPID::kElectron]>fCutMinElectronProbTPC :
1266     pid[AliPID::kElectron]>fCutMinElectronProbESD;
1267   //
1268 }
1269
1270 //____________________________________________________________________
1271 AliITSRecPoint* AliITSMultReconstructor::GetRecPoint(Int_t lr, Int_t n) const
1272 {
1273   // return a cluster of lr corresponding to orderer cluster index n
1274   if (fClArr[lr] && fClusterCopyIndex[lr] && n<fNClustersLay[lr])
1275     return (AliITSRecPoint*) fClArr[lr]->At(fClusterCopyIndex[lr][n]);
1276   else {
1277     AliError("To access the clusters SetCreateClustersCopy should have been called");
1278     return 0;
1279   }
1280 }