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