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   float vtxf[3] = {vtx->GetX(),vtx->GetY(),vtx->GetZ()};
 394   FindTracklets(vtxf);
 395   //
 396   CreateMultiplicityObject();
 397 }
 398
 399 //____________________________________________________________________
 400 void AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) {
 401   //
 402   // RS NOTE - this is old reconstructor invocation, to be used from VertexFinder
 403   //
 404   // - calls LoadClusterArray that finds the position of the clusters
 405   //   (in global coord)
 406   // - convert the cluster coordinates to theta, phi (seen from the
 407   //   interaction vertex).
 408   // - makes an array of tracklets
 409   //
 410   // After this method has been called, the clusters of the two layers
 411   // and the tracklets can be retrieved by calling the Get'er methods.
 412   if (fMult) delete fMult; fMult = 0;
 413   fNClustersLay1 = 0;
 414   fNClustersLay2 = 0;
 415   fNTracklets = 0;
 416   fNSingleCluster = 0;
 417   //
 418   if (!clusterTree) { AliError(" Invalid ITS cluster tree !\n"); return; }
 419   //
 420   fESDEvent = 0;
 421   fTreeRP = clusterTree;
 422   //
 423   FindTracklets(vtx);
 424   //
 425 }
 426
 427 //____________________________________________________________________
 428 void AliITSMultReconstructor::FindTracklets(const Float_t *vtx)
 429 {
 430   // Find tracklets converging to vertex
 431   //
 432   LoadClusterArrays(fTreeRP);
 433   // flag clusters used by ESD tracks
 434   if (fESDEvent) ProcessESDTracks();
 435
 436   if (!vtx) return;
 437
 438   const Double_t pi = TMath::Pi();
 439
 440   // dPhi shift is field dependent
 441   // get average magnetic field
 442   Float_t bz = 0;
 443   AliMagF* field = 0;
 444   if (TGeoGlobalMagField::Instance()) field = dynamic_cast<AliMagF*>(TGeoGlobalMagField::Instance()->GetField());
 445   if (!field)
 446   {
 447     AliError("Could not retrieve magnetic field. Assuming no field. Delta Phi shift will be deactivated in AliITSMultReconstructor.")
 448   }
 449   else
 450     bz = TMath::Abs(field->SolenoidField());
 451
 452   const Double_t dPhiShift = fPhiShift / 5 * bz;
 453   AliDebug(1, Form("Using phi shift of %f", dPhiShift));
 454
 455   const Double_t dPhiWindow2 = fPhiWindow * fPhiWindow;
 456   const Double_t dThetaWindow2 = fThetaWindow * fThetaWindow;
 457
 458   Int_t* partners = new Int_t[fNClustersLay2];
 459   Float_t* minDists = new Float_t[fNClustersLay2];
 460   Int_t* associatedLay1 = new Int_t[fNClustersLay1];
 461   TArrayI** blacklist = new TArrayI*[fNClustersLay1];
 462
 463   for (Int_t i=0; i<fNClustersLay2; i++) {
 464     partners[i] = -1;
 465     minDists[i] = 2;
 466   }
 467   for (Int_t i=0; i<fNClustersLay1; i++)
 468     associatedLay1[i] = 0;
 469   for (Int_t i=0; i<fNClustersLay1; i++)
 470     blacklist[i] = 0;
 471
 472   // find the tracklets
 473   AliDebug(1,"Looking for tracklets... ");
 474
 475   //###########################################################
 476   // Loop on layer 1 : finding theta, phi and z
 477   for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
 478     float *clPar = GetClusterLayer1(iC1);
 479     Float_t x = clPar[kClTh] - vtx[0];
 480     Float_t y = clPar[kClPh] - vtx[1];
 481     Float_t z = clPar[kClZ]  - vtx[2];
 482
 483     Float_t r    = TMath::Sqrt(x*x + y*y + z*z);
 484
 485     clPar[kClTh] = TMath::ACos(z/r);                   // Store Theta
 486     clPar[kClPh] = TMath::Pi() + TMath::ATan2(-y,-x);  // Store Phi
 487
 488     if (fHistOn) {
 489       Float_t eta = clPar[kClTh];
 490       eta= TMath::Tan(eta/2.);
 491       eta=-TMath::Log(eta);
 492       fhetaClustersLay1->Fill(eta);
 493       fhphiClustersLay1->Fill(clPar[kClPh]);
 494     }
 495   }
 496
 497   // Loop on layer 2 : finding theta, phi and r
 498   for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
 499     float *clPar = GetClusterLayer2(iC2);
 500     Float_t x = clPar[kClTh] - vtx[0];
 501     Float_t y = clPar[kClPh] - vtx[1];
 502     Float_t z = clPar[kClZ]  - vtx[2];
 503
 504     Float_t r    = TMath::Sqrt(x*x + y*y + z*z);
 505
 506     clPar[kClTh] = TMath::ACos(z/r);                   // Store Theta
 507     clPar[kClPh] = TMath::Pi() + TMath::ATan2(-y,-x);  // Store Phi
 508   }
 509
 510   //###########################################################
 511   Int_t found = 1;
 512   while (found > 0) {
 513     found = 0;
 514
 515     // Step1: find all tracklets allowing double assocation
 516     // Loop on layer 1
 517     for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
 518
 519       // already used or in the overlap ?
 520       if (associatedLay1[iC1] != 0 || fOverlapFlagClustersLay1[iC1]) continue;
 521
 522       found++;
 523
 524       // reset of variables for multiple candidates
 525       Int_t  iC2WithBestDist = -1;   // reset
 526       Double_t minDist       =  2;   // reset
 527       float* clPar1 = GetClusterLayer1(iC1);
 528
 529       // Loop on layer 2
 530       for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
 531
 532         // in the overlap ?
 533         if (fOverlapFlagClustersLay2[iC2]) continue;
 534         float* clPar2 = GetClusterLayer2(iC2);
 535
 536         if (blacklist[iC1]) {
 537           Bool_t blacklisted = kFALSE;
 538           for (Int_t i=blacklist[iC1]->GetSize(); i--;) {
 539             if (blacklist[iC1]->At(i) == iC2) {
 540               blacklisted = kTRUE;
 541               break;
 542             }
 543           }
 544           if (blacklisted) continue;
 545         }
 546
 547         // find the difference in angles
 548         Double_t dTheta = TMath::Abs(clPar2[kClTh] - clPar1[kClTh]);
 549         Double_t dPhi   = TMath::Abs(clPar2[kClPh] - clPar1[kClPh]);
 550         // take into account boundary condition
 551         if (dPhi>pi) dPhi=2.*pi-dPhi;
 552
 553         if (fHistOn) {
 554           fhClustersDPhiAll->Fill(dPhi);
 555           fhClustersDThetaAll->Fill(dTheta);
 556           fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
 557         }
 558
 559         dPhi -= dPhiShift;
 560
 561         // make "elliptical" cut in Phi and Theta!
 562         Float_t d = dPhi*dPhi/dPhiWindow2 + dTheta*dTheta/dThetaWindow2;
 563
 564         // look for the minimum distance: the minimum is in iC2WithBestDist
 565         if (d<1 && d<minDist) {
 566           minDist=d;
 567           iC2WithBestDist = iC2;
 568         }
 569       } // end of loop over clusters in layer 2
 570
 571       if (minDist<1) { // This means that a cluster in layer 2 was found that matches with iC1
 572
 573         if (minDists[iC2WithBestDist] > minDist) {
 574           Int_t oldPartner = partners[iC2WithBestDist];
 575           partners[iC2WithBestDist] = iC1;
 576           minDists[iC2WithBestDist] = minDist;
 577
 578           // mark as assigned
 579           associatedLay1[iC1] = 1;
 580
 581           if (oldPartner != -1) {
 582             // redo partner search for cluster in L0 (oldPartner), putting this one (iC2WithBestDist) on its blacklist
 583             if (blacklist[oldPartner] == 0) {
 584               blacklist[oldPartner] = new TArrayI(1);
 585             } else blacklist[oldPartner]->Set(blacklist[oldPartner]->GetSize()+1);
 586
 587             blacklist[oldPartner]->AddAt(iC2WithBestDist, blacklist[oldPartner]->GetSize()-1);
 588
 589             // mark as free
 590             associatedLay1[oldPartner] = 0;
 591           }
 592         } else {
 593           // try again to find a cluster without considering iC2WithBestDist
 594           if (blacklist[iC1] == 0) {
 595             blacklist[iC1] = new TArrayI(1);
 596           }
 597           else
 598             blacklist[iC1]->Set(blacklist[iC1]->GetSize()+1);
 599
 600           blacklist[iC1]->AddAt(iC2WithBestDist, blacklist[iC1]->GetSize()-1);
 601         }
 602
 603       } else // cluster has no partner; remove
 604         associatedLay1[iC1] = 2;
 605     } // end of loop over clusters in layer 1
 606   }
 607
 608   // Step2: store tracklets; remove used clusters
 609   for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
 610
 611     if (partners[iC2] == -1) continue;
 612
 613     if (fRemoveClustersFromOverlaps) FlagClustersInOverlapRegions (partners[iC2],iC2);
 614
 615
 616     if (fOverlapFlagClustersLay1[partners[iC2]] || fOverlapFlagClustersLay2[iC2]) continue;
 617
 618     float* clPar2 = GetClusterLayer2(iC2);
 619     float* clPar1 = GetClusterLayer1(partners[iC2]);
 620
 621     Float_t* tracklet = fTracklets[fNTracklets] = new Float_t[kTrNPar]; // RS Add also the cluster id's
 622
 623     // use the theta from the clusters in the first layer
 624     tracklet[kTrTheta] = clPar1[kClTh];
 625     // use the phi from the clusters in the first layer
 626     tracklet[kTrPhi] = clPar1[kClPh];
 627     // store the difference between phi1 and phi2
 628     tracklet[kTrDPhi] = clPar1[kClPh] - clPar2[kClPh];
 629
 630     // define dphi in the range [0,pi] with proper sign (track charge correlated)
 631     if (tracklet[kTrDPhi] > TMath::Pi())   tracklet[kTrDPhi] = tracklet[kTrDPhi]-2.*TMath::Pi();
 632     if (tracklet[kTrDPhi] < -TMath::Pi())  tracklet[kTrDPhi] = tracklet[kTrDPhi]+2.*TMath::Pi();
 633
 634     // store the difference between theta1 and theta2
 635     tracklet[kTrDTheta] = clPar1[kClTh] - clPar2[kClTh];
 636
 637     if (fHistOn) {
 638       fhClustersDPhiAcc->Fill(tracklet[kTrDPhi]);
 639       fhClustersDThetaAcc->Fill(tracklet[kTrDTheta]);
 640       fhDPhiVsDThetaAcc->Fill(tracklet[kTrDTheta],tracklet[kTrDPhi]);
 641     }
 642
 643     // find label
 644     // if equal label in both clusters found this label is assigned
 645     // if no equal label can be found the first labels of the L1 AND L2 cluster are assigned
 646     Int_t label1 = 0;
 647     Int_t label2 = 0;
 648     while (label2 < 3) {
 649       if ((Int_t) clPar1[kClMC0+label1] != -2 && (Int_t) clPar1[kClMC0+label1] == (Int_t) clPar2[kClMC0+label2])
 650         break;
 651       label1++;
 652       if (label1 == 3) {
 653         label1 = 0;
 654         label2++;
 655       }
 656     }
 657     if (label2 < 3) {
 658       AliDebug(AliLog::kDebug, Form("Found label %d == %d for tracklet candidate %d\n", (Int_t) clPar1[kClMC0+label1], (Int_t) clPar1[kClMC0+label2], fNTracklets));
 659       tracklet[kTrLab1] = clPar1[kClMC0+label1];
 660       tracklet[kTrLab2] = clPar2[kClMC0+label2];
 661     } else {
 662       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));
 663       tracklet[kTrLab1] = clPar1[kClMC0];
 664       tracklet[kTrLab2] = clPar2[kClMC0];
 665     }
 666
 667     if (fHistOn) {
 668       Float_t eta = tracklet[kTrTheta];
 669       eta= TMath::Tan(eta/2.);
 670       eta=-TMath::Log(eta);
 671       fhetaTracklets->Fill(eta);
 672       fhphiTracklets->Fill(tracklet[kTrPhi]);
 673     }
 674     //
 675     tracklet[kClID1] = partners[iC2];
 676     tracklet[kClID2] = iC2;
 677     //
 678     AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
 679     AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", partners[iC2], iC2));
 680     fNTracklets++;
 681
 682     associatedLay1[partners[iC2]] = 1;
 683   }
 684
 685   // Delete the following else if you do not want to save Clusters!
 686   // store the cluster
 687   for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
 688
 689     float* clPar1 = GetClusterLayer1(iC1);
 690
 691     if (associatedLay1[iC1]==2||associatedLay1[iC1]==0) {
 692       fSClusters[fNSingleCluster] = new Float_t[kClNPar];
 693       fSClusters[fNSingleCluster][kSCTh] = clPar1[kClTh];
 694       fSClusters[fNSingleCluster][kSCPh] = clPar1[kClPh];
 695       fSClusters[fNSingleCluster][kSCID] = iC1;
 696       AliDebug(1,Form(" Adding a single cluster %d (cluster %d  of layer 1)",
 697                 fNSingleCluster, iC1));
 698       fNSingleCluster++;
 699     }
 700   }
 701
 702   delete[] partners;
 703   delete[] minDists;
 704
 705   for (Int_t i=0; i<fNClustersLay1; i++)
 706     if (blacklist[i])
 707       delete blacklist[i];
 708   delete[] blacklist;
 709
 710   AliDebug(1,Form("%d tracklets found", fNTracklets));
 711 }
 712
 713 //____________________________________________________________________
 714 void AliITSMultReconstructor::CreateMultiplicityObject()
 715 {
 716   // create AliMultiplicity object and store it in the ESD event
 717   //
 718   TBits fastOrFiredMap,firedChipMap;
 719   if (fDetTypeRec) {
 720    fastOrFiredMap  = fDetTypeRec->GetFastOrFiredMap();
 721    firedChipMap    = fDetTypeRec->GetFiredChipMap(fTreeRP);
 722   }
 723   //
 724   fMult = new AliMultiplicity(fNTracklets,fNSingleCluster,fNFiredChips[0],fNFiredChips[1],fastOrFiredMap);
 725   fMult->SetFiredChipMap(firedChipMap);
 726   AliITSRecPointContainer* rcont = AliITSRecPointContainer::Instance();
 727   fMult->SetITSClusters(0,rcont->GetNClustersInLayer(1,fTreeRP));
 728   for(Int_t kk=2;kk<=6;kk++) fMult->SetITSClusters(kk-1,rcont->GetNClustersInLayerFast(kk));
 729   //
 730   for (int i=fNTracklets;i--;)  {
 731     float* tlInfo = fTracklets[i];
 732     fMult->SetTrackletData(i,tlInfo, fUsedClusLay1[int(tlInfo[kClID1])]|fUsedClusLay2[int(tlInfo[kClID2])]);
 733   }
 734   //
 735   for (int i=fNSingleCluster;i--;) {
 736     float* clInfo = fSClusters[i];
 737     fMult->SetSingleClusterData(i,clInfo,fUsedClusLay1[int(clInfo[kSCID])]);
 738   }
 739   fMult->CompactBits();
 740   //
 741 }
 742
 743
 744 //____________________________________________________________________
 745 void AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree)
 746 {
 747   // This method
 748   // - gets the clusters from the cluster tree
 749   // - convert them into global coordinates
 750   // - store them in the internal arrays
 751   // - count the number of cluster-fired chips
 752   //
 753   // RS: This method was strongly modified wrt original by Jan Fiete. In order to have the same numbering
 754   // of clusters as in the ITS reco I had to introduce sorting in Z
 755   // Also note that now the clusters data are stored not in float[6] attached to float**, but in 1-D array
 756
 757   AliDebug(1,"Loading clusters and cluster-fired chips ...");
 758
 759   fNClustersLay1 = 0;
 760   fNClustersLay2 = 0;
 761   fNFiredChips[0] = 0;
 762   fNFiredChips[1] = 0;
 763
 764   AliITSsegmentationSPD seg;
 765
 766   AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
 767   TClonesArray* itsClusters=rpcont->FetchClusters(0,itsClusterTree);
 768   if(!rpcont->IsSPDActive()){
 769     AliWarning("No SPD rec points found, multiplicity not calculated");
 770     return;
 771   }
 772   //
 773   // count clusters
 774   // loop over the SPD subdetectors
 775   TObjArray clArr(100);
 776   for (int il=0;il<2;il++) {
 777     int nclLayer = 0;
 778     int detMin = AliITSgeomTGeo::GetModuleIndex(il+1,1,1);
 779     int detMax = AliITSgeomTGeo::GetModuleIndex(il+2,1,1);
 780     for (int idt=detMin;idt<detMax;idt++) {
 781       itsClusters=rpcont->UncheckedGetClusters(idt);
 782       int nClusters = itsClusters->GetEntriesFast();
 783       if (!nClusters) continue;
 784       Int_t nClustersInChip[5] = {0,0,0,0,0};
 785       while(nClusters--) {
 786         AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
 787         if (!cluster) continue;
 788         clArr.AddAtAndExpand(cluster,nclLayer++);
 789         nClustersInChip[ seg.GetChipFromLocal(0,cluster->GetDetLocalZ()) ]++;
 790       }
 791       for(Int_t ifChip=5;ifChip--;) if (nClustersInChip[ifChip]) fNFiredChips[il]++;
 792     }
 793     // sort the clusters in Z (to have the same numbering as in ITS reco
 794     Float_t *z    = new Float_t[nclLayer];
 795     Int_t * index = new Int_t[nclLayer];
 796     for (int ic=0;ic<nclLayer;ic++) z[ic] = ((AliITSRecPoint*)clArr[ic])->GetZ();
 797     TMath::Sort(nclLayer,z,index,kFALSE);
 798     Float_t*   clustersLay              = new Float_t[nclLayer*kClNPar];
 799     Int_t*     detectorIndexClustersLay = new Int_t[nclLayer];
 800     Bool_t*    overlapFlagClustersLay   = new Bool_t[nclLayer];
 801     Char_t*    usedClusLay              = new Char_t[nclLayer];
 802     //
 803     for (int ic=0;ic<nclLayer;ic++) {
 804       AliITSRecPoint* cluster = (AliITSRecPoint*)clArr[index[ic]];
 805       float* clPar = &clustersLay[ic*kClNPar];
 806       //
 807       cluster->GetGlobalXYZ( clPar );
 808       detectorIndexClustersLay[ic] = cluster->GetDetectorIndex();
 809       overlapFlagClustersLay[ic]   = kFALSE;
 810       usedClusLay[ic]              = 0;
 811       for (Int_t i=3;i--;) clPar[kClMC0+i] = cluster->GetLabel(i);
 812     }
 813     clArr.Clear();
 814     delete[] z;
 815     delete[] index;
 816     //
 817     if (il==0) {
 818       fClustersLay1              = clustersLay;
 819       fOverlapFlagClustersLay1   = overlapFlagClustersLay;
 820       fDetectorIndexClustersLay1 = detectorIndexClustersLay;
 821       fUsedClusLay1              = usedClusLay;
 822       fNClustersLay1             = nclLayer;
 823     }
 824     else {
 825       fClustersLay2              = clustersLay;
 826       fOverlapFlagClustersLay2   = overlapFlagClustersLay;
 827       fDetectorIndexClustersLay2 = detectorIndexClustersLay;
 828       fUsedClusLay2              = usedClusLay;
 829       fNClustersLay2             = nclLayer;
 830     }
 831   }
 832   //
 833   // no double association allowed
 834   int nmaxT                  = TMath::Min(fNClustersLay1, fNClustersLay2);
 835   fTracklets                 = new Float_t*[nmaxT];
 836   fSClusters                 = new Float_t*[fNClustersLay1];
 837   for (Int_t i=nmaxT;i--;) fTracklets[i] = 0;
 838   //
 839   AliDebug(1,Form("(clusters in layer 1 : %d,  layer 2: %d)",fNClustersLay1,fNClustersLay2));
 840   AliDebug(1,Form("(cluster-fired chips in layer 1 : %d,  layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
 841 }
 842 //____________________________________________________________________
 843 void
 844 AliITSMultReconstructor::LoadClusterFiredChips(TTree* itsClusterTree) {
 845   // This method
 846   // - gets the clusters from the cluster tree
 847   // - counts the number of (cluster)fired chips
 848
 849   AliDebug(1,"Loading cluster-fired chips ...");
 850
 851   fNFiredChips[0] = 0;
 852   fNFiredChips[1] = 0;
 853
 854   AliITSsegmentationSPD seg;
 855   AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
 856   TClonesArray* itsClusters=rpcont->FetchClusters(0,itsClusterTree);
 857   if(!rpcont->IsSPDActive()){
 858     AliWarning("No SPD rec points found, multiplicity not calculated");
 859     return;
 860   }
 861
 862   // loop over the its subdetectors
 863   Int_t nSPDmodules=AliITSgeomTGeo::GetModuleIndex(3,1,1);
 864   for (Int_t iIts=0; iIts < nSPDmodules; iIts++) {
 865     itsClusters=rpcont->UncheckedGetClusters(iIts);
 866     Int_t nClusters = itsClusters->GetEntriesFast();
 867
 868     // number of clusters in each chip of the current module
 869     Int_t nClustersInChip[5] = {0,0,0,0,0};
 870     Int_t layer = 0;
 871
 872     // loop over clusters
 873     while(nClusters--) {
 874       AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
 875
 876       layer = cluster->GetLayer();
 877       if (layer>1) continue;
 878
 879       // find the chip for the current cluster
 880       Float_t locz = cluster->GetDetLocalZ();
 881       Int_t iChip = seg.GetChipFromLocal(0,locz);
 882       nClustersInChip[iChip]++;
 883
 884     }// end of cluster loop
 885
 886     // get number of fired chips in the current module
 887     for(Int_t ifChip=0; ifChip<5; ifChip++) {
 888       if(nClustersInChip[ifChip] >= 1)  fNFiredChips[layer]++;
 889     }
 890
 891   } // end of its "subdetector" loop
 892
 893
 894   AliDebug(1,Form("(cluster-fired chips in layer 1 : %d,  layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
 895 }
 896 //____________________________________________________________________
 897 void
 898 AliITSMultReconstructor::SaveHists() {
 899   // This method save the histograms on the output file
 900   // (only if fHistOn is TRUE).
 901
 902   if (!fHistOn)
 903     return;
 904
 905   fhClustersDPhiAll->Write();
 906   fhClustersDThetaAll->Write();
 907   fhDPhiVsDThetaAll->Write();
 908
 909   fhClustersDPhiAcc->Write();
 910   fhClustersDThetaAcc->Write();
 911   fhDPhiVsDThetaAcc->Write();
 912
 913   fhetaTracklets->Write();
 914   fhphiTracklets->Write();
 915   fhetaClustersLay1->Write();
 916   fhphiClustersLay1->Write();
 917 }
 918
 919 //____________________________________________________________________
 920 void
 921 AliITSMultReconstructor::FlagClustersInOverlapRegions (Int_t iC1, Int_t iC2WithBestDist) {
 922
 923   Float_t distClSameMod=0.;
 924   Float_t distClSameModMin=0.;
 925   Int_t   iClOverlap =0;
 926   Float_t meanRadiusLay1 = 3.99335; // average radius inner layer
 927   Float_t meanRadiusLay2 = 7.37935; // average radius outer layer;
 928
 929   Float_t zproj1=0.;
 930   Float_t zproj2=0.;
 931   Float_t deZproj=0.;
 932   Float_t* clPar1  = GetClusterLayer1(iC1);
 933   Float_t* clPar2B = GetClusterLayer2(iC2WithBestDist);
 934   // Loop on inner layer clusters
 935   for (Int_t iiC1=0; iiC1<fNClustersLay1; iiC1++) {
 936     if (!fOverlapFlagClustersLay1[iiC1]) {
 937       // only for adjacent modules
 938       if ((TMath::Abs(fDetectorIndexClustersLay1[iC1]-fDetectorIndexClustersLay1[iiC1])==4)||
 939          (TMath::Abs(fDetectorIndexClustersLay1[iC1]-fDetectorIndexClustersLay1[iiC1])==76)) {
 940         Float_t *clPar11 = GetClusterLayer1(iiC1);
 941         Float_t dePhi=TMath::Abs(clPar11[kClPh]-clPar1[kClPh]);
 942         if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
 943
 944         zproj1=meanRadiusLay1/TMath::Tan(clPar1[kClTh]);
 945         zproj2=meanRadiusLay1/TMath::Tan(clPar11[kClTh]);
 946
 947         deZproj=TMath::Abs(zproj1-zproj2);
 948
 949         distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
 950         if (distClSameMod<=1.) fOverlapFlagClustersLay1[iiC1]=kTRUE;
 951
 952 //        if (distClSameMod<=1.) {
 953 //          if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
 954 //            distClSameModMin=distClSameMod;
 955 //            iClOverlap=iiC1;
 956 //          }
 957 //        }
 958
 959
 960       } // end adjacent modules
 961     }
 962   } // end Loop on inner layer clusters
 963
 964 //  if (distClSameModMin!=0.) fOverlapFlagClustersLay1[iClOverlap]=kTRUE;
 965
 966   distClSameMod=0.;
 967   distClSameModMin=0.;
 968   iClOverlap =0;
 969   // Loop on outer layer clusters
 970   for (Int_t iiC2=0; iiC2<fNClustersLay2; iiC2++) {
 971     if (!fOverlapFlagClustersLay2[iiC2]) {
 972       // only for adjacent modules
 973       Float_t *clPar2 = GetClusterLayer2(iiC2);
 974       if ((TMath::Abs(fDetectorIndexClustersLay2[iC2WithBestDist]-fDetectorIndexClustersLay2[iiC2])==4) ||
 975          (TMath::Abs(fDetectorIndexClustersLay2[iC2WithBestDist]-fDetectorIndexClustersLay2[iiC2])==156)) {
 976         Float_t dePhi=TMath::Abs(clPar2[kClPh]-clPar2B[kClPh]);
 977         if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
 978
 979         zproj1=meanRadiusLay2/TMath::Tan(clPar2B[kClTh]);
 980         zproj2=meanRadiusLay2/TMath::Tan(clPar2[kClTh]);
 981
 982         deZproj=TMath::Abs(zproj1-zproj2);
 983         distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
 984         if (distClSameMod<=1.) fOverlapFlagClustersLay2[iiC2]=kTRUE;
 985
 986 //        if (distClSameMod<=1.) {
 987 //          if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
 988 //            distClSameModMin=distClSameMod;
 989 //            iClOverlap=iiC2;
 990 //          }
 991 //        }
 992
 993       } // end adjacent modules
 994     }
 995   } // end Loop on outer layer clusters
 996
 997 //  if (distClSameModMin!=0.) fOverlapFlagClustersLay2[iClOverlap]=kTRUE;
 998
 999 }
1000
1001 //____________________________________________________________________
1002 void AliITSMultReconstructor::ProcessESDTracks()
1003 {
1004   // Flag the clusters used by ESD tracks
1005   // Flag primary tracks to be used for multiplicity counting
1006   //
1007   if (!fESDEvent) return;
1008   AliESDVertex* vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexTracks();
1009   if (!vtx || vtx->GetNContributors()<1) vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexSPD();
1010   if (!vtx || vtx->GetNContributors()<1) {
1011     AliDebug(1,"No primary vertex: cannot flag primary tracks");
1012     return;
1013   }
1014   Int_t ntracks = fESDEvent->GetNumberOfTracks();
1015   for(Int_t itr=0; itr<ntracks; itr++) {
1016     AliESDtrack* track = fESDEvent->GetTrack(itr);
1017     if (!track->IsOn(AliESDtrack::kITSin)) continue; // use only tracks propagated in ITS to vtx
1018     FlagTrackClusters(track);
1019     FlagIfSecondary(track,vtx);
1020   }
1021   FlagV0s(vtx);
1022   //
1023 }
1024
1025 //____________________________________________________________________
1026 void AliITSMultReconstructor::FlagTrackClusters(const AliESDtrack* track)
1027 {
1028   // RS: flag the SPD clusters of the track if it is useful for the multiplicity estimation
1029   //
1030   Int_t idx[12];
1031   if ( track->GetITSclusters(idx)<3 ) return; // at least 3 clusters must be used in the fit
1032   //
1033   char mark = track->IsOn(AliESDtrack::kITSpureSA) ? kITSSAPBit : kITSTPCBit;
1034   char *uClus[2] = {fUsedClusLay1,fUsedClusLay2};
1035   for (int i=AliESDfriendTrack::kMaxITScluster;i--;) {
1036     // note: i>=6 is for extra clusters
1037     if (idx[i]<0) continue;
1038     int layID= (idx[i] & 0xf0000000) >> 28;
1039     if (layID>1) continue; // SPD only
1040     int clID = (idx[i] & 0x0fffffff);
1041     uClus[layID][clID] |= mark;
1042   }
1043   //
1044 }
1045
1046 //____________________________________________________________________
1047 void AliITSMultReconstructor::FlagIfSecondary(AliESDtrack* track, const AliVertex* vtx)
1048 {
1049   // RS: check if the track is primary and set the flag
1050   double cut = (track->HasPointOnITSLayer(0)||track->HasPointOnITSLayer(1)) ? fCutPxDrSPDin:fCutPxDrSPDout;
1051   float xz[2];
1052   track->GetDZ(vtx->GetX(),vtx->GetY(),vtx->GetZ(), fESDEvent->GetMagneticField(), xz);
1053   if (TMath::Abs(xz[0]*track->P())>cut || TMath::Abs(xz[1]*track->P())>fCutPxDz ||
1054       TMath::Abs(xz[0])>fCutDCArz   || TMath::Abs(xz[1])>fCutDCArz)
1055     track->SetStatus(AliESDtrack::kMultSec);
1056   else track->ResetStatus(AliESDtrack::kMultSec);
1057 }
1058
1059 //____________________________________________________________________
1060 void AliITSMultReconstructor::FlagV0s(const AliESDVertex *vtx)
1061 {
1062   // flag tracks belonging to v0s
1063   //
1064   const double kK0Mass = 0.4976;
1065   //
1066   AliV0 pvertex;
1067   AliKFVertex vertexKF;
1068   AliKFParticle epKF0,epKF1,pipmKF0,piKF0,piKF1,gammaKF,k0KF;
1069   Double_t mass,massErr,chi2c;
1070   enum {kKFIni=BIT(14)};
1071   //
1072   double recVtx[3];
1073   float recVtxF[3];
1074   vtx->GetXYZ(recVtx);
1075   for (int i=3;i--;) recVtxF[i] = recVtx[i];
1076   //
1077   int ntracks = fESDEvent->GetNumberOfTracks();
1078   if (ntracks<2) return;
1079   //
1080   vertexKF.X() = recVtx[0];
1081   vertexKF.Y() = recVtx[1];
1082   vertexKF.Z() = recVtx[2];
1083   vertexKF.Covariance(0,0) = vtx->GetXRes()*vtx->GetXRes();
1084   vertexKF.Covariance(1,2) = vtx->GetYRes()*vtx->GetYRes();
1085   vertexKF.Covariance(2,2) = vtx->GetZRes()*vtx->GetZRes();
1086   //
1087   AliESDtrack *trc0,*trc1;
1088   for (int it0=0;it0<ntracks;it0++) {
1089     trc0 = fESDEvent->GetTrack(it0);
1090     if (trc0->IsOn(AliESDtrack::kMultInV0)) continue;
1091     if (!trc0->IsOn(AliESDtrack::kITSin)) continue;
1092     Bool_t isSAP = trc0->IsPureITSStandalone();
1093     Int_t  q0 = trc0->Charge();
1094     Bool_t testGamma = CanBeElectron(trc0);
1095     epKF0.ResetBit(kKFIni);
1096     piKF0.ResetBit(kKFIni);
1097     double bestChi2=1e16;
1098     int bestID = -1;
1099     //
1100     for (int it1=it0+1;it1<ntracks;it1++) {
1101       trc1 = fESDEvent->GetTrack(it1);
1102       if (trc1->IsOn(AliESDtrack::kMultInV0)) continue;
1103       if (!trc1->IsOn(AliESDtrack::kITSin)) continue;
1104       if (trc1->IsPureITSStandalone() != isSAP) continue; // pair separately ITS_SA_Pure tracks and TPC/ITS+ITS_SA
1105       if ( (q0+trc1->Charge())!=0 ) continue;             // don't pair like signs
1106       //
1107       pvertex.SetParamN(q0<0 ? *trc0:*trc1);
1108       pvertex.SetParamP(q0>0 ? *trc0:*trc1);
1109       pvertex.Update(recVtxF);
1110       if (pvertex.P()<fCutMinP) continue;
1111       if (pvertex.GetV0CosineOfPointingAngle()<fCutMinPointAngle) continue;
1112       if (pvertex.GetDcaV0Daughters()>fCutMaxDCADauther) continue;
1113       double d = pvertex.GetD(recVtx[0],recVtx[1],recVtx[2]);
1114       if (d>fCutMaxDCA) continue;
1115       double dx=recVtx[0]-pvertex.Xv(), dy=recVtx[1]-pvertex.Yv();
1116       double rv = TMath::Sqrt(dx*dx+dy*dy);
1117       //
1118       // check gamma conversion hypothesis ----------------------------------------------------------->>>
1119       Bool_t gammaOK = kFALSE;
1120       while (testGamma && CanBeElectron(trc1)) {
1121         if (rv<fCutMinRGamma) break;
1122         if (!epKF0.TestBit(kKFIni)) {
1123           new(&epKF0) AliKFParticle(*trc0,q0>0 ? kPositron:kElectron);
1124           epKF0.SetBit(kKFIni);
1125         }
1126         new(&epKF1) AliKFParticle(*trc1,q0<0 ? kPositron:kElectron);
1127         gammaKF.Initialize();
1128         gammaKF += epKF0;
1129         gammaKF += epKF1;
1130         gammaKF.SetProductionVertex(vertexKF);
1131         gammaKF.GetMass(mass,massErr);
1132         if (mass>fCutMassGamma || (massErr>0&&(mass>massErr*fCutMassGammaNSigma))) break;
1133         if (gammaKF.GetS()<fCutGammaSFromDecay) break;
1134         gammaKF.SetMassConstraint(0.,0.001);
1135         chi2c = (gammaKF.GetNDF()!=0) ? gammaKF.GetChi2()/gammaKF.GetNDF() : 1000;
1136         if (chi2c>fCutChi2cGamma) break;
1137         gammaOK = kTRUE;
1138         if (chi2c>bestChi2) break;
1139         bestChi2 = chi2c;
1140         bestID = it1;
1141         break;
1142       }
1143       if (gammaOK) continue;
1144       // check gamma conversion hypothesis -----------------------------------------------------------<<<
1145       // check K0 conversion hypothesis    ----------------------------------------------------------->>>
1146       while (1) {
1147         if (rv<fCutMinRK0) break;
1148         if (!piKF0.TestBit(kKFIni)) {
1149           new(&piKF0) AliKFParticle(*trc0,q0>0 ? kPiPlus:kPiMinus);
1150           piKF0.SetBit(kKFIni);
1151         }
1152         new(&piKF1) AliKFParticle(*trc1,q0<0 ? kPiPlus:kPiMinus);
1153         k0KF.Initialize();
1154         k0KF += piKF0;
1155         k0KF += piKF1;
1156         k0KF.SetProductionVertex(vertexKF);
1157         k0KF.GetMass(mass,massErr);
1158         mass -= kK0Mass;
1159         if (TMath::Abs(mass)>fCutMassK0 || (massErr>0&&(abs(mass)>massErr*fCutMassK0NSigma))) break;
1160         if (k0KF.GetS()<fCutK0SFromDecay) break;
1161         k0KF.SetMassConstraint(kK0Mass,0.001);
1162         chi2c = (k0KF.GetNDF()!=0) ? k0KF.GetChi2()/k0KF.GetNDF() : 1000;
1163         if (chi2c>fCutChi2cK0) break;
1164         if (chi2c>bestChi2) break;
1165         bestChi2 = chi2c;
1166         bestID = it1;
1167         break;
1168       }
1169       // check K0 conversion hypothesis    -----------------------------------------------------------<<<
1170     }
1171     //
1172     if (bestID>=0) {
1173       trc0->SetStatus(AliESDtrack::kMultInV0);
1174       fESDEvent->GetTrack(bestID)->SetStatus(AliESDtrack::kMultInV0);
1175     }
1176   }
1177   //
1178 }
1179
1180 //____________________________________________________________________
1181 Bool_t AliITSMultReconstructor::CanBeElectron(const AliESDtrack* trc) const
1182 {
1183   // check if the track can be electron
1184   Double_t pid[AliPID::kSPECIES];
1185   if (!trc->IsOn(AliESDtrack::kESDpid)) return kTRUE;
1186   trc->GetESDpid(pid);
1187   return (trc->IsOn(AliESDtrack::kTPCpid)) ?
1188     pid[AliPID::kElectron]>fCutMinElectronProbTPC :
1189     pid[AliPID::kElectron]>fCutMinElectronProbESD;
1190   //
1191 }