PWGPP/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] = {static_cast<float>(vtx->GetX()),static_cast<float>(vtx->GetY()),static_cast<float>(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   delete[] associatedLay1;
 692
 693   for (Int_t i=0; i<fNClustersLay1; i++)
 694     if (blacklist[i])
 695       delete blacklist[i];
 696   delete[] blacklist;
 697 //  Printf("exiting...");
 698   AliDebug(1,Form("%d tracklets found", fNTracklets));
 699 }
 700
 701 //____________________________________________________________________
 702 void AliTrackletAlg::CreateMultiplicityObject()
 703 {
 704   // create AliMultiplicity object and store it in the ESD event
 705   //
 706   TBits fastOrFiredMap,firedChipMap;
 707   if (fDetTypeRec) {
 708    fastOrFiredMap  = fDetTypeRec->GetFastOrFiredMap();
 709    firedChipMap    = fDetTypeRec->GetFiredChipMap(fTreeRP);
 710   }
 711   //
 712   fMult = new AliMultiplicity(fNTracklets,fNSingleCluster,fNFiredChips[0],fNFiredChips[1],fastOrFiredMap);
 713   fMult->SetFiredChipMap(firedChipMap);
 714   AliITSRecPointContainer* rcont = AliITSRecPointContainer::Instance();
 715   fMult->SetITSClusters(0,rcont->GetNClustersInLayer(1,fTreeRP));
 716   for(Int_t kk=2;kk<=6;kk++) fMult->SetITSClusters(kk-1,rcont->GetNClustersInLayerFast(kk));
 717   //
 718   for (int i=fNTracklets;i--;)  {
 719     float* tlInfo = fTracklets[i];
 720     fMult->SetTrackletData(i,tlInfo, fUsedClusLay1[int(tlInfo[kClID1])],fUsedClusLay2[int(tlInfo[kClID2])]);
 721   }
 722   //
 723   for (int i=fNSingleCluster;i--;) {
 724     float* clInfo = fSClusters[i];
 725     fMult->SetSingleClusterData(i,clInfo,fUsedClusLay1[int(clInfo[kSCID])]);
 726   }
 727   fMult->CompactBits();
 728   //
 729 }
 730
 731
 732 //____________________________________________________________________
 733 void AliTrackletAlg::LoadClusterArrays(TTree* itsClusterTree)
 734 {
 735   // This method
 736   // - gets the clusters from the cluster tree
 737   // - convert them into global coordinates
 738   // - store them in the internal arrays
 739   // - count the number of cluster-fired chips
 740   //
 741   // RS: This method was strongly modified wrt original. In order to have the same numbering
 742   // of clusters as in the ITS reco I had to introduce sorting in Z
 743   // Also note that now the clusters data are stored not in float[6] attached to float**, but in 1-D array
 744   AliDebug(1,"Loading clusters and cluster-fired chips ...");
 745
 746   fNClustersLay1 = 0;
 747   fNClustersLay2 = 0;
 748   fNFiredChips[0] = 0;
 749   fNFiredChips[1] = 0;
 750
 751   AliITSsegmentationSPD seg;
 752
 753 //  AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
 754 //  TClonesArray* itsClusters=rpcont->FetchClusters(0,itsClusterTree);
 755 //  if(!rpcont->IsSPDActive()){
 756 //    AliWarning("No SPD rec points found, multiplicity not calculated");
 757 //    return;
 758 //  }
 759
 760   TClonesArray statITSrec("AliITSRecPoint");
 761   TClonesArray* itsClusters= &statITSrec;
 762
 763   TBranch* branch=itsClusterTree->GetBranch("ITSRecPoints");
 764   if(!branch) {
 765     printf("NO itsClusterTree branch available. Skipping...\n");
 766     return;
 767   }
 768
 769   branch->SetAddress(&itsClusters);
 770   // count clusters
 771   // loop over the SPD subdetectors
 772   static TClonesArray clArr("AliITSRecPoint",100);
 773 //  Float_t cluGlo[3] = {0.,0.,0.};
 774   for (int il=0;il<2;il++) {
 775     int nclLayer = 0;
 776     int detMin = AliITSgeomTGeo::GetModuleIndex(il+1,1,1);
 777     int detMax = AliITSgeomTGeo::GetModuleIndex(il+2,1,1);
 778     for (int idt=detMin;idt<detMax;idt++) {
 779       branch->GetEvent(idt);
 780       int nClusters = itsClusters->GetEntriesFast();
 781 //      itsClusters=rpcont->UncheckedGetClusters(idt);
 782       if (!nClusters) continue;
 783       Int_t nClustersInChip[5] = {0,0,0,0,0};
 784       while(nClusters--) {
 785         AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
 786         if (!cluster) continue;
 787 /*        cluster->GetGlobalXYZ(cluGlo);
 788         if (idt==0)
 789                    Printf("First Cl1 LoadClArr %f %f %f ",cluGlo[0],cluGlo[1],cluGlo[2]);
 790         if (idt==80)
 791                    Printf("First Cl2 LoadClArr %f %f %f ",cluGlo[0],cluGlo[1],cluGlo[2]);*/
 792         new (clArr[nclLayer++]) AliITSRecPoint(*cluster);
 793         // find the chip for the current cluster
 794         Float_t locz = cluster->GetDetLocalZ();
 795         Int_t iChip = seg.GetChipFromLocal(0,locz);
 796         if (iChip>=0)
 797           nClustersInChip[iChip]++;
 798       }
 799       for(Int_t ifChip=5;ifChip--;) if (nClustersInChip[ifChip]) fNFiredChips[il]++;
 800     }
 801     // sort the clusters in Z (to have the same numbering as in ITS reco
 802     Float_t *z    = new Float_t[nclLayer];
 803     Int_t * index = new Int_t[nclLayer];
 804     for (int ic=0;ic<nclLayer;ic++) z[ic] = ((AliITSRecPoint*)clArr[ic])->GetZ();
 805     TMath::Sort(nclLayer,z,index,kFALSE);
 806     Float_t*   clustersLay              = new Float_t[nclLayer*kClNPar];
 807     Int_t*     detectorIndexClustersLay = new Int_t[nclLayer];
 808     Bool_t*    overlapFlagClustersLay   = new Bool_t[nclLayer];
 809     UInt_t*    usedClusLay              = new UInt_t[nclLayer];
 810     //
 811     for (int ic=0;ic<nclLayer;ic++) {
 812       AliITSRecPoint* cluster = (AliITSRecPoint*)clArr[index[ic]];
 813       float* clPar = &clustersLay[ic*kClNPar];
 814       //
 815       cluster->GetGlobalXYZ( clPar );
 816 //      if (ic==0 && detMin==0) Printf("First Cl1 LoadClArrSorted %f %f %f ",clPar[kClTh],clPar[kClPh],clPar[kClZ]);
 817 //      if (ic==0 && detMin==80) Printf("First Cl2 LoadClArrSorted %f %f %f ",clPar[kClTh],clPar[kClPh],clPar[kClZ]);
 818 //      if (ic==0 && detMin==0) Printf("First Cl1 LoadClArrSorted %f %f %f ",clPar[0],clPar[1],clPar[2]);
 819 //      if (ic==0 && detMin==80) Printf("First Cl2 LoadClArrSorted %f %f %f ",clPar[0],clPar[1],clPar[2]);
 820       detectorIndexClustersLay[ic] = cluster->GetDetectorIndex();
 821       overlapFlagClustersLay[ic]   = kFALSE;
 822       usedClusLay[ic]              = 0;
 823       for (Int_t i=3;i--;) clPar[kClMC0+i] = cluster->GetLabel(i);
 824     }
 825     clArr.Clear();
 826     delete[] z;
 827     delete[] index;
 828     //
 829     if (il==0) {
 830       fClustersLay1              = clustersLay;
 831       fOverlapFlagClustersLay1   = overlapFlagClustersLay;
 832       fDetectorIndexClustersLay1 = detectorIndexClustersLay;
 833       fUsedClusLay1              = usedClusLay;
 834       fNClustersLay1             = nclLayer;
 835     }
 836     else {
 837       fClustersLay2              = clustersLay;
 838       fOverlapFlagClustersLay2   = overlapFlagClustersLay;
 839       fDetectorIndexClustersLay2 = detectorIndexClustersLay;
 840       fUsedClusLay2              = usedClusLay;
 841       fNClustersLay2             = nclLayer;
 842     }
 843   }
 844 //  Printf("First Cl1 %f %f %f ",fClustersLay1[0],fClustersLay1[1],fClustersLay1[2]);
 845 //  Printf("First Cl2 %f %f %f ",fClustersLay2[0],fClustersLay2[1],fClustersLay2[2]);
 846   Printf("LoadClusterArr: N cl1 %d",fNClustersLay1);
 847   Printf("LoadClusterArrN: N cl2 %d",fNClustersLay2);
 848   //
 849   // no double association allowed // ??
 850   int nmaxT                  = TMath::Min(fNClustersLay1, fNClustersLay2);
 851   fTracklets                 = new Float_t*[nmaxT];
 852   fSClusters                 = new Float_t*[fNClustersLay1];
 853   for (Int_t i=nmaxT;i--;) fTracklets[i] = 0;
 854   //
 855   AliDebug(1,Form("(clusters in layer 1 : %d,  layer 2: %d)",fNClustersLay1,fNClustersLay2));
 856   AliDebug(1,Form("(cluster-fired chips in layer 1 : %d,  layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
 857 }
 858 //____________________________________________________________________
 859 void
 860 AliTrackletAlg::LoadClusterFiredChips(TTree* itsClusterTree) {
 861   // This method
 862   // - gets the clusters from the cluster tree
 863   // - counts the number of (cluster)fired chips
 864
 865   AliDebug(1,"Loading cluster-fired chips ...");
 866
 867   fNFiredChips[0] = 0;
 868   fNFiredChips[1] = 0;
 869
 870   AliITSsegmentationSPD seg;
 871   AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
 872   TClonesArray* itsClusters=NULL;
 873   rpcont->FetchClusters(0,itsClusterTree);
 874   if(!rpcont->IsSPDActive()){
 875     AliWarning("No SPD rec points found, multiplicity not calculated");
 876     return;
 877   }
 878
 879   // loop over the its subdetectors
 880   Int_t nSPDmodules=AliITSgeomTGeo::GetModuleIndex(3,1,1);
 881   for (Int_t iIts=0; iIts < nSPDmodules; iIts++) {
 882     itsClusters=rpcont->UncheckedGetClusters(iIts);
 883     Int_t nClusters = itsClusters->GetEntriesFast();
 884
 885     // number of clusters in each chip of the current module
 886     Int_t nClustersInChip[5] = {0,0,0,0,0};
 887     Int_t layer = 0;
 888
 889     // loop over clusters
 890     while(nClusters--) {
 891       AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
 892
 893       layer = cluster->GetLayer();
 894       if (layer>1) break; // no point in further check for this module
 895
 896       // find the chip for the current cluster
 897       Float_t locz = cluster->GetDetLocalZ();
 898       Int_t iChip = seg.GetChipFromLocal(0,locz);
 899       if (iChip>=0)
 900         nClustersInChip[iChip]++;
 901
 902     }// end of cluster loop
 903     if (layer>1) continue; // make sure we are in the SPD layer
 904     // get number of fired chips in the current module
 905     for(Int_t ifChip=0; ifChip<5; ifChip++) {
 906       if(nClustersInChip[ifChip] >= 1)  fNFiredChips[layer]++;
 907     }
 908
 909   } // end of its "subdetector" loop
 910
 911
 912   AliDebug(1,Form("(cluster-fired chips in layer 1 : %d,  layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
 913 }
 914 //____________________________________________________________________
 915 void
 916 AliTrackletAlg::SaveHists() {
 917   // This method save the histograms on the output file
 918   // (only if fHistOn is TRUE).
 919
 920   if (!fHistOn)
 921     return;
 922
 923   fhClustersDPhiAll->Write();
 924   fhClustersDThetaAll->Write();
 925   fhDPhiVsDThetaAll->Write();
 926
 927   fhClustersDPhiAcc->Write();
 928   fhClustersDThetaAcc->Write();
 929   fhDPhiVsDThetaAcc->Write();
 930
 931   fhetaTracklets->Write();
 932   fhphiTracklets->Write();
 933   fhetaClustersLay1->Write();
 934   fhphiClustersLay1->Write();
 935 }
 936
 937 //____________________________________________________________________
 938 void
 939 AliTrackletAlg::FlagClustersInOverlapRegions (Int_t iC1, Int_t iC2WithBestDist) {
 940
 941   Float_t distClSameMod=0.;
 942   Float_t distClSameModMin=0.;
 943   Int_t   iClOverlap =0;
 944   Float_t meanRadiusLay1 = 3.99335; // average radius inner layer
 945   Float_t meanRadiusLay2 = 7.37935; // average radius outer layer;
 946
 947   Float_t zproj1=0.;
 948   Float_t zproj2=0.;
 949   Float_t deZproj=0.;
 950   Float_t* clPar1  = GetClusterLayer1(iC1);
 951   Float_t* clPar2B = GetClusterLayer2(iC2WithBestDist);
 952   // Loop on inner layer clusters
 953   for (Int_t iiC1=0; iiC1<fNClustersLay1; iiC1++) {
 954     if (!fOverlapFlagClustersLay1[iiC1]) {
 955       // only for adjacent modules
 956       if ((TMath::Abs(fDetectorIndexClustersLay1[iC1]-fDetectorIndexClustersLay1[iiC1])==4)||
 957          (TMath::Abs(fDetectorIndexClustersLay1[iC1]-fDetectorIndexClustersLay1[iiC1])==76)) {
 958         Float_t *clPar11 = GetClusterLayer1(iiC1);
 959         Float_t dePhi=TMath::Abs(clPar11[kClPh]-clPar1[kClPh]);
 960         if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
 961
 962         zproj1=meanRadiusLay1/TMath::Tan(clPar1[kClTh]);
 963         zproj2=meanRadiusLay1/TMath::Tan(clPar11[kClTh]);
 964
 965         deZproj=TMath::Abs(zproj1-zproj2);
 966
 967         distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
 968         if (distClSameMod<=1.) fOverlapFlagClustersLay1[iiC1]=kTRUE;
 969
 970 //        if (distClSameMod<=1.) {
 971 //          if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
 972 //            distClSameModMin=distClSameMod;
 973 //            iClOverlap=iiC1;
 974 //          }
 975 //        }
 976
 977
 978       } // end adjacent modules
 979     }
 980   } // end Loop on inner layer clusters
 981
 982 //  if (distClSameModMin!=0.) fOverlapFlagClustersLay1[iClOverlap]=kTRUE;
 983
 984   distClSameMod=0.;
 985   distClSameModMin=0.;
 986   iClOverlap =0;
 987   // Loop on outer layer clusters
 988   for (Int_t iiC2=0; iiC2<fNClustersLay2; iiC2++) {
 989     if (!fOverlapFlagClustersLay2[iiC2]) {
 990       // only for adjacent modules
 991       Float_t *clPar2 = GetClusterLayer2(iiC2);
 992       if ((TMath::Abs(fDetectorIndexClustersLay2[iC2WithBestDist]-fDetectorIndexClustersLay2[iiC2])==4) ||
 993          (TMath::Abs(fDetectorIndexClustersLay2[iC2WithBestDist]-fDetectorIndexClustersLay2[iiC2])==156)) {
 994         Float_t dePhi=TMath::Abs(clPar2[kClPh]-clPar2B[kClPh]);
 995         if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
 996
 997         zproj1=meanRadiusLay2/TMath::Tan(clPar2B[kClTh]);
 998         zproj2=meanRadiusLay2/TMath::Tan(clPar2[kClTh]);
 999
1000         deZproj=TMath::Abs(zproj1-zproj2);
1001         distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
1002         if (distClSameMod<=1.) fOverlapFlagClustersLay2[iiC2]=kTRUE;
1003
1004 //        if (distClSameMod<=1.) {
1005 //          if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
1006 //            distClSameModMin=distClSameMod;
1007 //            iClOverlap=iiC2;
1008 //          }
1009 //        }
1010
1011       } // end adjacent modules
1012     }
1013   } // end Loop on outer layer clusters
1014
1015 //  if (distClSameModMin!=0.) fOverlapFlagClustersLay2[iClOverlap]=kTRUE;
1016
1017 }
1018
1019 //____________________________________________________________________
1020 void AliTrackletAlg::ProcessESDTracks()
1021 {
1022   // Flag the clusters used by ESD tracks
1023   // Flag primary tracks to be used for multiplicity counting
1024   //
1025   if (!fESDEvent) return;
1026   AliESDVertex* vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexTracks();
1027   if (!vtx || vtx->GetNContributors()<1) vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexSPD();
1028   if (!vtx || vtx->GetNContributors()<1) {
1029     AliDebug(1,"No primary vertex: cannot flag primary tracks");
1030     return;
1031   }
1032   Int_t ntracks = fESDEvent->GetNumberOfTracks();
1033   for(Int_t itr=0; itr<ntracks; itr++) {
1034     AliESDtrack* track = fESDEvent->GetTrack(itr);
1035     if (!track->IsOn(AliESDtrack::kITSin)) continue; // use only tracks propagated in ITS to vtx
1036     FlagTrackClusters(itr);
1037     FlagIfSecondary(track,vtx);
1038   }
1039   FlagV0s(vtx);
1040   //
1041 }
1042
1043 //____________________________________________________________________
1044 void AliTrackletAlg::FlagTrackClusters(Int_t id)
1045 {
1046   // RS: flag the SPD clusters of the track if it is useful for the multiplicity estimation
1047   //
1048   const UInt_t   kMaskL = 0x0000ffff;
1049   const UInt_t   kMaskH = 0xffff0000;
1050   const UInt_t   kMaxTrID = kMaskL - 1; // max possible track id
1051   if (UInt_t(id)>kMaxTrID) return;
1052   const AliESDtrack* track = fESDEvent->GetTrack(id);
1053   Int_t idx[12];
1054   if ( track->GetITSclusters(idx)<3 ) return; // at least 3 clusters must be used in the fit
1055   UInt_t *uClus[2] = {fUsedClusLay1,fUsedClusLay2};
1056   //
1057   UInt_t mark = id+1;
1058   if (track->IsOn(AliESDtrack::kITSpureSA)) mark <<= 16;
1059   //
1060   for (int i=AliESDfriendTrack::kMaxITScluster;i--;) {
1061     // note: i>=6 is for extra clusters
1062     if (idx[i]<0) continue;
1063     int layID= (idx[i] & 0xf0000000) >> 28;
1064     if (layID>1) continue; // SPD only
1065     int clID = (idx[i] & 0x0fffffff);
1066     //
1067     if ( track->IsOn(AliESDtrack::kITSpureSA) ) {
1068       if (uClus[layID][clID]&kMaskH) {
1069         AliWarning(Form("Tracks %5d and %5d share cluster %6d of lr%d",id,int(uClus[layID][clID]>>16)-1,clID,layID));
1070         uClus[layID][clID] &= kMaskL;
1071       }
1072     }
1073     else if (uClus[layID][clID]&kMaskL) {
1074       AliWarning(Form("Tracks %5d and %5d share cluster %6d of lr%d",id,int(uClus[layID][clID]&kMaskL)-1,clID,layID));
1075       uClus[layID][clID] &= kMaskH;
1076     }
1077     uClus[layID][clID] |= mark;
1078   }
1079   //
1080 }
1081
1082 //____________________________________________________________________
1083 void AliTrackletAlg::FlagIfSecondary(AliESDtrack* track, const AliVertex* vtx)
1084 {
1085   // RS: check if the track is primary and set the flag
1086   double cut = (track->HasPointOnITSLayer(0)||track->HasPointOnITSLayer(1)) ? fCutPxDrSPDin:fCutPxDrSPDout;
1087   float xz[2];
1088   track->GetDZ(vtx->GetX(),vtx->GetY(),vtx->GetZ(), fESDEvent->GetMagneticField(), xz);
1089   if (TMath::Abs(xz[0]*track->P())>cut || TMath::Abs(xz[1]*track->P())>fCutPxDz ||
1090       TMath::Abs(xz[0])>fCutDCArz   || TMath::Abs(xz[1])>fCutDCArz)
1091     track->SetStatus(AliESDtrack::kMultSec);
1092   else track->ResetStatus(AliESDtrack::kMultSec);
1093 }
1094
1095 //____________________________________________________________________
1096 void AliTrackletAlg::FlagV0s(const AliESDVertex *vtx)
1097 {
1098   // flag tracks belonging to v0s
1099   //
1100   const double kK0Mass = 0.4976;
1101   //
1102   AliV0 pvertex;
1103   AliKFVertex vertexKF;
1104   AliKFParticle epKF0,epKF1,pipmKF0,piKF0,piKF1,gammaKF,k0KF;
1105   Double_t mass,massErr,chi2c;
1106   enum {kKFIni=BIT(14)};
1107   //
1108   double recVtx[3];
1109   float recVtxF[3];
1110   vtx->GetXYZ(recVtx);
1111   for (int i=3;i--;) recVtxF[i] = recVtx[i];
1112   //
1113   int ntracks = fESDEvent->GetNumberOfTracks();
1114   if (ntracks<2) return;
1115   //
1116   vertexKF.X() = recVtx[0];
1117   vertexKF.Y() = recVtx[1];
1118   vertexKF.Z() = recVtx[2];
1119   vertexKF.Covariance(0,0) = vtx->GetXRes()*vtx->GetXRes();
1120   vertexKF.Covariance(1,2) = vtx->GetYRes()*vtx->GetYRes();
1121   vertexKF.Covariance(2,2) = vtx->GetZRes()*vtx->GetZRes();
1122   //
1123   AliESDtrack *trc0,*trc1;
1124   for (int it0=0;it0<ntracks;it0++) {
1125     trc0 = fESDEvent->GetTrack(it0);
1126     if (trc0->IsOn(AliESDtrack::kMultInV0)) continue;
1127     if (!trc0->IsOn(AliESDtrack::kITSin)) continue;
1128     Bool_t isSAP = trc0->IsPureITSStandalone();
1129     Int_t  q0 = trc0->Charge();
1130     Bool_t testGamma = CanBeElectron(trc0);
1131     epKF0.ResetBit(kKFIni);
1132     piKF0.ResetBit(kKFIni);
1133     double bestChi2=1e16;
1134     int bestID = -1;
1135     //
1136     for (int it1=it0+1;it1<ntracks;it1++) {
1137       trc1 = fESDEvent->GetTrack(it1);
1138       if (trc1->IsOn(AliESDtrack::kMultInV0)) continue;
1139       if (!trc1->IsOn(AliESDtrack::kITSin)) continue;
1140       if (trc1->IsPureITSStandalone() != isSAP) continue; // pair separately ITS_SA_Pure tracks and TPC/ITS+ITS_SA
1141       if ( (q0+trc1->Charge())!=0 ) continue;             // don't pair like signs
1142       //
1143       pvertex.SetParamN(q0<0 ? *trc0:*trc1);
1144       pvertex.SetParamP(q0>0 ? *trc0:*trc1);
1145       pvertex.Update(recVtxF);
1146       if (pvertex.P()<fCutMinP) continue;
1147       if (pvertex.GetV0CosineOfPointingAngle()<fCutMinPointAngle) continue;
1148       if (pvertex.GetDcaV0Daughters()>fCutMaxDCADauther) continue;
1149       double d = pvertex.GetD(recVtx[0],recVtx[1],recVtx[2]);
1150       if (d>fCutMaxDCA) continue;
1151       double dx=recVtx[0]-pvertex.Xv(), dy=recVtx[1]-pvertex.Yv();
1152       double rv = TMath::Sqrt(dx*dx+dy*dy);
1153       //
1154       // check gamma conversion hypothesis ----------------------------------------------------------->>>
1155       Bool_t gammaOK = kFALSE;
1156       while (testGamma && CanBeElectron(trc1)) {
1157         if (rv<fCutMinRGamma) break;
1158         if (!epKF0.TestBit(kKFIni)) {
1159           new(&epKF0) AliKFParticle(*trc0,q0>0 ? kPositron:kElectron);
1160           epKF0.SetBit(kKFIni);
1161         }
1162         new(&epKF1) AliKFParticle(*trc1,q0<0 ? kPositron:kElectron);
1163         gammaKF.Initialize();
1164         gammaKF += epKF0;
1165         gammaKF += epKF1;
1166         gammaKF.SetProductionVertex(vertexKF);
1167         gammaKF.GetMass(mass,massErr);
1168         if (mass>fCutMassGamma || (massErr>0&&(mass>massErr*fCutMassGammaNSigma))) break;
1169         if (gammaKF.GetS()<fCutGammaSFromDecay) break;
1170         gammaKF.SetMassConstraint(0.,0.001);
1171         chi2c = (gammaKF.GetNDF()!=0) ? gammaKF.GetChi2()/gammaKF.GetNDF() : 1000;
1172         if (chi2c>fCutChi2cGamma) break;
1173         gammaOK = kTRUE;
1174         if (chi2c>bestChi2) break;
1175         bestChi2 = chi2c;
1176         bestID = it1;
1177         break;
1178       }
1179       if (gammaOK) continue;
1180       // check gamma conversion hypothesis -----------------------------------------------------------<<<
1181       // check K0 conversion hypothesis    ----------------------------------------------------------->>>
1182       while (1) {
1183         if (rv<fCutMinRK0) break;
1184         if (!piKF0.TestBit(kKFIni)) {
1185           new(&piKF0) AliKFParticle(*trc0,q0>0 ? kPiPlus:kPiMinus);
1186           piKF0.SetBit(kKFIni);
1187         }
1188         new(&piKF1) AliKFParticle(*trc1,q0<0 ? kPiPlus:kPiMinus);
1189         k0KF.Initialize();
1190         k0KF += piKF0;
1191         k0KF += piKF1;
1192         k0KF.SetProductionVertex(vertexKF);
1193         k0KF.GetMass(mass,massErr);
1194         mass -= kK0Mass;
1195         if (TMath::Abs(mass)>fCutMassK0 || (massErr>0&&(abs(mass)>massErr*fCutMassK0NSigma))) break;
1196         if (k0KF.GetS()<fCutK0SFromDecay) break;
1197         k0KF.SetMassConstraint(kK0Mass,0.001);
1198         chi2c = (k0KF.GetNDF()!=0) ? k0KF.GetChi2()/k0KF.GetNDF() : 1000;
1199         if (chi2c>fCutChi2cK0) break;
1200         if (chi2c>bestChi2) break;
1201         bestChi2 = chi2c;
1202         bestID = it1;
1203         break;
1204       }
1205       // check K0 conversion hypothesis    -----------------------------------------------------------<<<
1206     }
1207     //
1208     if (bestID>=0) {
1209       trc0->SetStatus(AliESDtrack::kMultInV0);
1210       fESDEvent->GetTrack(bestID)->SetStatus(AliESDtrack::kMultInV0);
1211     }
1212   }
1213   //
1214 }
1215
1216 //____________________________________________________________________
1217 Bool_t AliTrackletAlg::CanBeElectron(const AliESDtrack* trc) const
1218 {
1219   // check if the track can be electron
1220   Double_t pid[AliPID::kSPECIES];
1221   if (!trc->IsOn(AliESDtrack::kESDpid)) return kTRUE;
1222   trc->GetESDpid(pid);
1223   return (trc->IsOn(AliESDtrack::kTPCpid)) ?
1224     pid[AliPID::kElectron]>fCutMinElectronProbTPC :
1225     pid[AliPID::kElectron]>fCutMinElectronProbESD;
1226   //
1227 }