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