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