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