Coverity fixes
[u/mrichter/AliRoot.git] / 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   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 }