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