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] = {vtx->GetX(),vtx->GetY(),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 distClSameModMin=0.;
 817   Int_t   iClOverlap =0;
 818   Float_t meanRadiusLay1 = 3.99335; // average radius inner layer
 819   Float_t meanRadiusLay2 = 7.37935; // average radius outer layer;
 820
 821   Float_t zproj1=0.;
 822   Float_t zproj2=0.;
 823   Float_t deZproj=0.;
 824   Float_t* clPar1  = GetClusterLayer1(iC1);
 825   Float_t* clPar2B = GetClusterLayer2(iC2WithBestDist);
 826   // Loop on inner layer clusters
 827   for (Int_t iiC1=0; iiC1<fNClustersLay[0]; iiC1++) {
 828     if (!fOverlapFlagClustersLay[0][iiC1]) {
 829       // only for adjacent modules
 830       if ((TMath::Abs(fDetectorIndexClustersLay[0][iC1]-fDetectorIndexClustersLay[0][iiC1])==4)||
 831          (TMath::Abs(fDetectorIndexClustersLay[0][iC1]-fDetectorIndexClustersLay[0][iiC1])==76)) {
 832         Float_t *clPar11 = GetClusterLayer1(iiC1);
 833         Float_t dePhi=TMath::Abs(clPar11[kClPh]-clPar1[kClPh]);
 834         if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
 835
 836         zproj1=meanRadiusLay1/TMath::Tan(clPar1[kClTh]);
 837         zproj2=meanRadiusLay1/TMath::Tan(clPar11[kClTh]);
 838
 839         deZproj=TMath::Abs(zproj1-zproj2);
 840
 841         distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
 842         if (distClSameMod<=1.) fOverlapFlagClustersLay[0][iiC1]=kTRUE;
 843
 844 //        if (distClSameMod<=1.) {
 845 //          if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
 846 //            distClSameModMin=distClSameMod;
 847 //            iClOverlap=iiC1;
 848 //          }
 849 //        }
 850
 851
 852       } // end adjacent modules
 853     }
 854   } // end Loop on inner layer clusters
 855
 856 //  if (distClSameModMin!=0.) fOverlapFlagClustersLay[0][iClOverlap]=kTRUE;
 857
 858   distClSameMod=0.;
 859   distClSameModMin=0.;
 860   iClOverlap =0;
 861   // Loop on outer layer clusters
 862   for (Int_t iiC2=0; iiC2<fNClustersLay[1]; iiC2++) {
 863     if (!fOverlapFlagClustersLay[1][iiC2]) {
 864       // only for adjacent modules
 865       Float_t *clPar2 = GetClusterLayer2(iiC2);
 866       if ((TMath::Abs(fDetectorIndexClustersLay[1][iC2WithBestDist]-fDetectorIndexClustersLay[1][iiC2])==4) ||
 867          (TMath::Abs(fDetectorIndexClustersLay[1][iC2WithBestDist]-fDetectorIndexClustersLay[1][iiC2])==156)) {
 868         Float_t dePhi=TMath::Abs(clPar2[kClPh]-clPar2B[kClPh]);
 869         if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
 870
 871         zproj1=meanRadiusLay2/TMath::Tan(clPar2B[kClTh]);
 872         zproj2=meanRadiusLay2/TMath::Tan(clPar2[kClTh]);
 873
 874         deZproj=TMath::Abs(zproj1-zproj2);
 875         distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
 876         if (distClSameMod<=1.) fOverlapFlagClustersLay[1][iiC2]=kTRUE;
 877
 878 //        if (distClSameMod<=1.) {
 879 //          if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
 880 //            distClSameModMin=distClSameMod;
 881 //            iClOverlap=iiC2;
 882 //          }
 883 //        }
 884
 885       } // end adjacent modules
 886     }
 887   } // end Loop on outer layer clusters
 888
 889 //  if (distClSameModMin!=0.) fOverlapFlagClustersLay[1][iClOverlap]=kTRUE;
 890
 891 }
 892
 893 //____________________________________________________________________
 894 void AliITSMultReconstructor::InitAux()
 895 {
 896   // init arrays/parameters for tracklet reconstruction
 897
 898   // dPhi shift is field dependent, get average magnetic field
 899   Float_t bz = 0;
 900   AliMagF* field = 0;
 901   if (TGeoGlobalMagField::Instance()) field = dynamic_cast<AliMagF*>(TGeoGlobalMagField::Instance()->GetField());
 902   if (!field) {
 903     AliError("Could not retrieve magnetic field. Assuming no field. Delta Phi shift will be deactivated in AliITSMultReconstructor.");
 904   }
 905   else bz = TMath::Abs(field->SolenoidField());
 906   fDPhiShift = fPhiShift / 5 * bz;
 907   AliDebug(1, Form("Using phi shift of %f", fDPhiShift));
 908   //
 909   if (fPartners) delete[] fPartners; fPartners = new Int_t[fNClustersLay[1]];
 910   if (fMinDists) delete[] fMinDists; fMinDists = new Float_t[fNClustersLay[1]];
 911   if (fAssociatedLay1) delete[] fAssociatedLay1; fAssociatedLay1 = new Int_t[fNClustersLay[0]];
 912   //
 913   if (fBlackList) delete fBlackList; fBlackList = new AliRefArray(fNClustersLay[0]);
 914   //
 915   //  Printf("Vertex in find tracklets...%f %f %f",vtx[0],vtx[1],vtx[2]);
 916   for (Int_t i=0; i<fNClustersLay[1]; i++) {
 917     fPartners[i] = -1;
 918     fMinDists[i] = 2*fNStdDev;
 919   }
 920   memset(fAssociatedLay1,0,fNClustersLay[0]*sizeof(Int_t));
 921   //
 922 }
 923
 924 //____________________________________________________________________
 925 void AliITSMultReconstructor::ClusterPos2Angles(const Float_t *vtx)
 926 {
 927   // convert cluster coordinates to angles wrt vertex
 928   for (int ilr=0;ilr<2;ilr++) {
 929     for (Int_t iC=0; iC<fNClustersLay[ilr]; iC++) {
 930       float* clPar = GetClusterOfLayer(ilr,iC);
 931       CalcThetaPhi(clPar[kClTh]-vtx[0],clPar[kClPh]-vtx[1],clPar[kClZ]-vtx[2],clPar[kClTh],clPar[kClPh]);
 932       if (ilr==0) {
 933         clPar[kClPh] = clPar[kClPh] + fPhiRotationAngle;   // rotation of inner layer for comb studies
 934         if (fHistOn) {
 935           Float_t eta = clPar[kClTh];
 936           eta= TMath::Tan(eta/2.);
 937           eta=-TMath::Log(eta);
 938           fhetaClustersLay1->Fill(eta);
 939           fhphiClustersLay1->Fill(clPar[kClPh]);
 940         }
 941       }
 942     }
 943   }
 944   //
 945 }
 946
 947 //____________________________________________________________________
 948 Int_t AliITSMultReconstructor::AssociateClusterOfL1(Int_t iC1)
 949 {
 950   // search association of cluster iC1 of L1 with all clusters of L2
 951   if (fAssociatedLay1[iC1] != 0) return 0;
 952   Int_t  iC2WithBestDist = -1;   // reset
 953   Double_t minDist       =  2*fNStdDev;   // reset
 954   float* clPar1 = GetClusterLayer1(iC1);
 955   for (Int_t iC2=0; iC2<fNClustersLay[1]; iC2++) {
 956     //
 957     if (fBlackList->IsReferred(iC1,iC2)) continue;
 958     float* clPar2 = GetClusterLayer2(iC2);
 959     //
 960     // find the difference in angles
 961     Double_t dTheta = TMath::Abs(clPar2[kClTh] - clPar1[kClTh]);
 962     Double_t dPhi   = TMath::Abs(clPar2[kClPh] - clPar1[kClPh]);
 963     //        Printf("detheta %f  dephi %f", dTheta,dPhi);
 964     //
 965     if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi;     // take into account boundary condition
 966     //
 967     if (fHistOn) {
 968       fhClustersDPhiAll->Fill(dPhi);
 969       fhClustersDThetaAll->Fill(dTheta);
 970       fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
 971     }
 972     Float_t d = CalcDist(dPhi,dTheta,clPar1[kClTh]);     // make "elliptical" cut in Phi and Theta!
 973     // look for the minimum distance: the minimum is in iC2WithBestDist
 974     if (d<fNStdDev && d<minDist) { minDist=d; iC2WithBestDist = iC2; }
 975   }
 976   //
 977   if (minDist<fNStdDev) { // This means that a cluster in layer 2 was found that matches with iC1
 978     //
 979     if (fMinDists[iC2WithBestDist] > minDist) {
 980       Int_t oldPartner = fPartners[iC2WithBestDist];
 981       fPartners[iC2WithBestDist] = iC1;
 982       fMinDists[iC2WithBestDist] = minDist;
 983       //
 984       fAssociatedLay1[iC1] = 1;      // mark as assigned
 985       //
 986       if (oldPartner != -1) {
 987         // redo partner search for cluster in L0 (oldPartner), putting this one (iC2WithBestDist) on its fBlackList
 988         fBlackList->AddReference(oldPartner,iC2WithBestDist);
 989         fAssociatedLay1[oldPartner] = 0;       // mark as free
 990       }
 991     } else {
 992       // try again to find a cluster without considering iC2WithBestDist
 993       fBlackList->AddReference(iC1,iC2WithBestDist);
 994     }
 995     //
 996   }
 997   else fAssociatedLay1[iC1] = 2;// cluster has no partner; remove
 998   //
 999   return 1;
1000 }
1001
1002 //____________________________________________________________________
1003 Int_t AliITSMultReconstructor::StoreTrackletForL2Cluster(Int_t iC2)
1004 {
1005   // build tracklet for cluster iC2 of layer 2
1006   if (fPartners[iC2] == -1) return 0;
1007   if (fRemoveClustersFromOverlaps) FlagClustersInOverlapRegions (fPartners[iC2],iC2);
1008   // Printf("saving tracklets");
1009   if (fOverlapFlagClustersLay[0][fPartners[iC2]] || fOverlapFlagClustersLay[1][iC2]) return 0;
1010   float* clPar2 = GetClusterLayer2(iC2);
1011   float* clPar1 = GetClusterLayer1(fPartners[iC2]);
1012   //
1013   Float_t* tracklet = fTracklets[fNTracklets] = new Float_t[kTrNPar]; // RS Add also the cluster id's
1014   //
1015   tracklet[kTrTheta] = clPar1[kClTh];    // use the theta from the clusters in the first layer
1016   tracklet[kTrPhi]   = clPar1[kClPh];    // use the phi from the clusters in the first layer
1017   tracklet[kTrDPhi] = clPar1[kClPh] - clPar2[kClPh];  // store the difference between phi1 and phi2
1018   //
1019   // define dphi in the range [0,pi] with proper sign (track charge correlated)
1020   if (tracklet[kTrDPhi] > TMath::Pi())   tracklet[kTrDPhi] = tracklet[kTrDPhi]-2.*TMath::Pi();
1021   if (tracklet[kTrDPhi] < -TMath::Pi())  tracklet[kTrDPhi] = tracklet[kTrDPhi]+2.*TMath::Pi();
1022   //
1023   tracklet[kTrDTheta] = clPar1[kClTh] - clPar2[kClTh]; // store the theta1-theta2
1024   //
1025   if (fHistOn) {
1026     fhClustersDPhiAcc->Fill(tracklet[kTrDPhi]);
1027     fhClustersDThetaAcc->Fill(tracklet[kTrDTheta]);
1028     fhDPhiVsDThetaAcc->Fill(tracklet[kTrDTheta],tracklet[kTrDPhi]);
1029   }
1030   //
1031   // find label
1032   // if equal label in both clusters found this label is assigned
1033   // if no equal label can be found the first labels of the L1 AND L2 cluster are assigned
1034   Int_t label1=0,label2=0;
1035   while (label2 < 3) {
1036     if ( int(clPar1[kClMC0+label1])!=-2 && int(clPar1[kClMC0+label1])==int(clPar2[kClMC0+label2])) break;
1037     if (++label1 == 3) { label1 = 0; label2++; }
1038   }
1039   if (label2 < 3) {
1040     AliDebug(AliLog::kDebug, Form("Found label %d == %d for tracklet candidate %d\n",
1041                                   (Int_t) clPar1[kClMC0+label1], (Int_t) clPar1[kClMC0+label2], fNTracklets));
1042     tracklet[kTrLab1] = tracklet[kTrLab2] = clPar1[kClMC0+label1];
1043   } else {
1044     AliDebug(AliLog::kDebug, Form("Did not find label %d %d %d %d %d %d for tracklet candidate %d\n",
1045                                   (Int_t) clPar1[kClMC0], (Int_t) clPar1[kClMC1], (Int_t) clPar1[kClMC2],
1046                                   (Int_t) clPar2[kClMC0], (Int_t) clPar2[kClMC1], (Int_t) clPar2[kClMC2], fNTracklets));
1047     tracklet[kTrLab1] = clPar1[kClMC0];
1048     tracklet[kTrLab2] = clPar2[kClMC0];
1049   }
1050   //
1051   if (fHistOn) {
1052     Float_t eta = tracklet[kTrTheta];
1053     eta= TMath::Tan(eta/2.);
1054     eta=-TMath::Log(eta);
1055     fhetaTracklets->Fill(eta);
1056     fhphiTracklets->Fill(tracklet[kTrPhi]);
1057   }
1058   //
1059   tracklet[kClID1] = fPartners[iC2];
1060   tracklet[kClID2] = iC2;
1061   //
1062   // Printf("Adding tracklet candidate");
1063   AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
1064   AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", fPartners[iC2], iC2));
1065   fNTracklets++;
1066   fAssociatedLay1[fPartners[iC2]] = 1;
1067   //
1068   return 1;
1069 }
1070
1071 //____________________________________________________________________
1072 void AliITSMultReconstructor::StoreL1Singles()
1073 {
1074   // Printf("saving single clusters...");
1075   for (Int_t iC1=0; iC1<fNClustersLay[0]; iC1++) {
1076     float* clPar1 = GetClusterLayer1(iC1);
1077     if (fAssociatedLay1[iC1]==2||fAssociatedLay1[iC1]==0) {
1078       fSClusters[fNSingleCluster] = new Float_t[kClNPar];
1079       fSClusters[fNSingleCluster][kSCTh] = clPar1[kClTh];
1080       fSClusters[fNSingleCluster][kSCPh] = clPar1[kClPh];
1081       fSClusters[fNSingleCluster][kSCLab] = clPar1[kClMC0];
1082       fSClusters[fNSingleCluster][kSCID] = iC1;
1083       AliDebug(1,Form(" Adding a single cluster %d (cluster %d  of layer 1)",
1084                       fNSingleCluster, iC1));
1085       fNSingleCluster++;
1086     }
1087   }
1088   //
1089   if (GetStoreSPD2SingleCl()) {
1090     for (Int_t iC2=0; iC2<fNClustersLay[1]; iC2++) {
1091       if (fPartners[iC2]<0 || (fOverlapFlagClustersLay[0][fPartners[iC2]] || fOverlapFlagClustersLay[1][iC2])) {
1092         float* clPar2 = GetClusterLayer2(iC2);
1093         fSClusters[fNSingleCluster] = new Float_t[kClNPar];
1094         fSClusters[fNSingleCluster][kSCTh] = clPar2[kClTh];
1095         fSClusters[fNSingleCluster][kSCPh] = clPar2[kClPh];
1096         fSClusters[fNSingleCluster][kSCLab] = clPar2[kClMC0];
1097         fSClusters[fNSingleCluster][kSCID] = iC2;
1098         AliDebug(1,Form(" Adding a single cluster %d (cluster %d  of layer 2)",
1099                         fNSingleCluster, iC2));
1100         fNSingleCluster++;
1101         fNSingleClusterSPD2++;
1102       }
1103     }
1104   }
1105   //
1106 }
1107
1108 //____________________________________________________________________
1109 void AliITSMultReconstructor::ProcessESDTracks()
1110 {
1111   // Flag the clusters used by ESD tracks
1112   // Flag primary tracks to be used for multiplicity counting
1113   //
1114   if (!fESDEvent || !fBuildRefs) return;
1115   AliESDVertex* vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexTracks();
1116   if (!vtx || vtx->GetNContributors()<1) vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexSPD();
1117   if (!vtx || vtx->GetNContributors()<1) {
1118     AliDebug(1,"No primary vertex: cannot flag primary tracks");
1119     return;
1120   }
1121   Int_t ntracks = fESDEvent->GetNumberOfTracks();
1122   for(Int_t itr=0; itr<ntracks; itr++) {
1123     AliESDtrack* track = fESDEvent->GetTrack(itr);
1124     if (!track->IsOn(AliESDtrack::kITSin)) continue; // use only tracks propagated in ITS to vtx
1125     FlagTrackClusters(itr);
1126     FlagIfSecondary(track,vtx);
1127   }
1128   FlagV0s(vtx);
1129   //
1130 }
1131
1132 //____________________________________________________________________
1133 void AliITSMultReconstructor::FlagTrackClusters(Int_t id)
1134 {
1135   // RS: flag the SPD clusters of the track if it is useful for the multiplicity estimation
1136   //
1137   const AliESDtrack* track = fESDEvent->GetTrack(id);
1138   Int_t idx[12];
1139   if ( track->GetITSclusters(idx)<3 ) return; // at least 3 clusters must be used in the fit
1140   Int_t itsType = track->IsOn(AliESDtrack::kITSpureSA) ? 1:0;
1141
1142   for (int i=6/*AliESDfriendTrack::kMaxITScluster*/;i--;) { // ignore extras: note: i>=6 is for extra clusters
1143     if (idx[i]<0) continue;
1144     int layID= (idx[i] & 0xf0000000) >> 28;
1145     if (layID>1) continue; // SPD only
1146     int clID = (idx[i] & 0x0fffffff);
1147     fUsedClusLay[layID][itsType]->AddReference(clID,id);
1148     fStoreRefs[layID][itsType] = kTRUE;
1149   }
1150   //
1151 }
1152
1153 //____________________________________________________________________
1154 void AliITSMultReconstructor::FlagIfSecondary(AliESDtrack* track, const AliVertex* vtx)
1155 {
1156   // RS: check if the track is primary and set the flag
1157   double cut = (track->HasPointOnITSLayer(0)||track->HasPointOnITSLayer(1)) ? fCutPxDrSPDin:fCutPxDrSPDout;
1158   float xz[2];
1159   track->GetDZ(vtx->GetX(),vtx->GetY(),vtx->GetZ(), fESDEvent->GetMagneticField(), xz);
1160   if (TMath::Abs(xz[0]*track->P())>cut || TMath::Abs(xz[1]*track->P())>fCutPxDz ||
1161       TMath::Abs(xz[0])>fCutDCArz   || TMath::Abs(xz[1])>fCutDCArz)
1162     track->SetStatus(AliESDtrack::kMultSec);
1163   else track->ResetStatus(AliESDtrack::kMultSec);
1164 }
1165
1166 //____________________________________________________________________
1167 void AliITSMultReconstructor::FlagV0s(const AliESDVertex *vtx)
1168 {
1169   // flag tracks belonging to v0s
1170   //
1171   const double kK0Mass = 0.4976;
1172   //
1173   AliV0 pvertex;
1174   AliKFVertex vertexKF;
1175   AliKFParticle epKF0,epKF1,pipmKF0,piKF0,piKF1,gammaKF,k0KF;
1176   Double_t mass,massErr,chi2c;
1177   enum {kKFIni=BIT(14)};
1178   //
1179   double recVtx[3];
1180   float recVtxF[3];
1181   vtx->GetXYZ(recVtx);
1182   for (int i=3;i--;) recVtxF[i] = recVtx[i];
1183   //
1184   int ntracks = fESDEvent->GetNumberOfTracks();
1185   if (ntracks<2) return;
1186   //
1187   vertexKF.X() = recVtx[0];
1188   vertexKF.Y() = recVtx[1];
1189   vertexKF.Z() = recVtx[2];
1190   vertexKF.Covariance(0,0) = vtx->GetXRes()*vtx->GetXRes();
1191   vertexKF.Covariance(1,2) = vtx->GetYRes()*vtx->GetYRes();
1192   vertexKF.Covariance(2,2) = vtx->GetZRes()*vtx->GetZRes();
1193   //
1194   AliESDtrack *trc0,*trc1;
1195   for (int it0=0;it0<ntracks;it0++) {
1196     trc0 = fESDEvent->GetTrack(it0);
1197     if (trc0->IsOn(AliESDtrack::kMultInV0)) continue;
1198     if (!trc0->IsOn(AliESDtrack::kITSin)) continue;
1199     Bool_t isSAP = trc0->IsPureITSStandalone();
1200     Int_t  q0 = trc0->Charge();
1201     Bool_t testGamma = CanBeElectron(trc0);
1202     epKF0.ResetBit(kKFIni);
1203     piKF0.ResetBit(kKFIni);
1204     double bestChi2=1e16;
1205     int bestID = -1;
1206     //
1207     for (int it1=it0+1;it1<ntracks;it1++) {
1208       trc1 = fESDEvent->GetTrack(it1);
1209       if (trc1->IsOn(AliESDtrack::kMultInV0)) continue;
1210       if (!trc1->IsOn(AliESDtrack::kITSin)) continue;
1211       if (trc1->IsPureITSStandalone() != isSAP) continue; // pair separately ITS_SA_Pure tracks and TPC/ITS+ITS_SA
1212       if ( (q0+trc1->Charge())!=0 ) continue;             // don't pair like signs
1213       //
1214       pvertex.SetParamN(q0<0 ? *trc0:*trc1);
1215       pvertex.SetParamP(q0>0 ? *trc0:*trc1);
1216       pvertex.Update(recVtxF);
1217       if (pvertex.P()<fCutMinP) continue;
1218       if (pvertex.GetV0CosineOfPointingAngle()<fCutMinPointAngle) continue;
1219       if (pvertex.GetDcaV0Daughters()>fCutMaxDCADauther) continue;
1220       double d = pvertex.GetD(recVtx[0],recVtx[1],recVtx[2]);
1221       if (d>fCutMaxDCA) continue;
1222       double dx=recVtx[0]-pvertex.Xv(), dy=recVtx[1]-pvertex.Yv();
1223       double rv = TMath::Sqrt(dx*dx+dy*dy);
1224       //
1225       // check gamma conversion hypothesis ----------------------------------------------------------->>>
1226       Bool_t gammaOK = kFALSE;
1227       while (testGamma && CanBeElectron(trc1)) {
1228         if (rv<fCutMinRGamma) break;
1229         if (!epKF0.TestBit(kKFIni)) {
1230           new(&epKF0) AliKFParticle(*trc0,q0>0 ? kPositron:kElectron);
1231           epKF0.SetBit(kKFIni);
1232         }
1233         new(&epKF1) AliKFParticle(*trc1,q0<0 ? kPositron:kElectron);
1234         gammaKF.Initialize();
1235         gammaKF += epKF0;
1236         gammaKF += epKF1;
1237         gammaKF.SetProductionVertex(vertexKF);
1238         gammaKF.GetMass(mass,massErr);
1239         if (mass>fCutMassGamma || (massErr>0&&(mass>massErr*fCutMassGammaNSigma))) break;
1240         if (gammaKF.GetS()<fCutGammaSFromDecay) break;
1241         gammaKF.SetMassConstraint(0.,0.001);
1242         chi2c = (gammaKF.GetNDF()!=0) ? gammaKF.GetChi2()/gammaKF.GetNDF() : 1000;
1243         if (chi2c>fCutChi2cGamma) break;
1244         gammaOK = kTRUE;
1245         if (chi2c>bestChi2) break;
1246         bestChi2 = chi2c;
1247         bestID = it1;
1248         break;
1249       }
1250       if (gammaOK) continue;
1251       // check gamma conversion hypothesis -----------------------------------------------------------<<<
1252       // check K0 conversion hypothesis    ----------------------------------------------------------->>>
1253       while (1) {
1254         if (rv<fCutMinRK0) break;
1255         if (!piKF0.TestBit(kKFIni)) {
1256           new(&piKF0) AliKFParticle(*trc0,q0>0 ? kPiPlus:kPiMinus);
1257           piKF0.SetBit(kKFIni);
1258         }
1259         new(&piKF1) AliKFParticle(*trc1,q0<0 ? kPiPlus:kPiMinus);
1260         k0KF.Initialize();
1261         k0KF += piKF0;
1262         k0KF += piKF1;
1263         k0KF.SetProductionVertex(vertexKF);
1264         k0KF.GetMass(mass,massErr);
1265         mass -= kK0Mass;
1266         if (TMath::Abs(mass)>fCutMassK0 || (massErr>0&&(abs(mass)>massErr*fCutMassK0NSigma))) break;
1267         if (k0KF.GetS()<fCutK0SFromDecay) break;
1268         k0KF.SetMassConstraint(kK0Mass,0.001);
1269         chi2c = (k0KF.GetNDF()!=0) ? k0KF.GetChi2()/k0KF.GetNDF() : 1000;
1270         if (chi2c>fCutChi2cK0) break;
1271         if (chi2c>bestChi2) break;
1272         bestChi2 = chi2c;
1273         bestID = it1;
1274         break;
1275       }
1276       // check K0 conversion hypothesis    -----------------------------------------------------------<<<
1277     }
1278     //
1279     if (bestID>=0) {
1280       trc0->SetStatus(AliESDtrack::kMultInV0);
1281       fESDEvent->GetTrack(bestID)->SetStatus(AliESDtrack::kMultInV0);
1282     }
1283   }
1284   //
1285 }
1286
1287 //____________________________________________________________________
1288 Bool_t AliITSMultReconstructor::CanBeElectron(const AliESDtrack* trc) const
1289 {
1290   // check if the track can be electron
1291   Double_t pid[AliPID::kSPECIES];
1292   if (!trc->IsOn(AliESDtrack::kESDpid)) return kTRUE;
1293   trc->GetESDpid(pid);
1294   return (trc->IsOn(AliESDtrack::kTPCpid)) ?
1295     pid[AliPID::kElectron]>fCutMinElectronProbTPC :
1296     pid[AliPID::kElectron]>fCutMinElectronProbESD;
1297   //
1298 }
1299
1300 //____________________________________________________________________
1301 AliITSRecPoint* AliITSMultReconstructor::GetRecPoint(Int_t lr, Int_t n) const
1302 {
1303   // return a cluster of lr corresponding to orderer cluster index n
1304   if (fClArr[lr] && fClusterCopyIndex[lr] && n<fNClustersLay[lr])
1305     return (AliITSRecPoint*) fClArr[lr]->At(fClusterCopyIndex[lr][n]);
1306   else {
1307     AliError("To access the clusters SetCreateClustersCopy should have been called");
1308     return 0;
1309   }
1310 }