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