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