ITS/AliITSMultReconstructor.cxx

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