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