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