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