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