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