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