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