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