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