Added protection against missing ESD pointer
[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 //_________________________________________________________________________
62
63 #include <TClonesArray.h>
64 #include <TH1F.h>
65 #include <TH2F.h>
66 #include <TTree.h>
67 #include <TBits.h>
68 #include <TArrayI.h>
69
70 #include "AliITSMultReconstructor.h"
71 #include "AliITSReconstructor.h"
72 #include "AliITSsegmentationSPD.h"
73 #include "AliITSRecPoint.h"
74 #include "AliITSRecPointContainer.h"
75 #include "AliITSgeom.h"
76 #include "AliITSgeomTGeo.h"
77 #include "AliITSDetTypeRec.h"
78 #include "AliESDEvent.h"
79 #include "AliESDVertex.h"
80 #include "AliESDtrack.h"
81 #include "AliMultiplicity.h"
82 #include "AliLog.h"
83 #include "TGeoGlobalMagField.h"
84 #include "AliMagF.h"
85
86 //____________________________________________________________________
87 ClassImp(AliITSMultReconstructor)
88
89
90 //____________________________________________________________________
91 AliITSMultReconstructor::AliITSMultReconstructor():
92 fDetTypeRec(0),fESDEvent(0),fTreeRP(0),fUsedClusLay1(0),fUsedClusLay2(0),
93 fClustersLay1(0),
94 fClustersLay2(0),
95 fDetectorIndexClustersLay1(0),
96 fDetectorIndexClustersLay2(0),
97 fOverlapFlagClustersLay1(0),
98 fOverlapFlagClustersLay2(0),
99 fTracklets(0),
100 fSClusters(0),
101 fNClustersLay1(0),
102 fNClustersLay2(0),
103 fNTracklets(0),
104 fNSingleCluster(0),
105 fPhiWindow(0),
106 fThetaWindow(0),
107 fPhiShift(0),
108 fRemoveClustersFromOverlaps(0),
109 fPhiOverlapCut(0),
110 fZetaOverlapCut(0),
111 fHistOn(0),
112 fhClustersDPhiAcc(0),
113 fhClustersDThetaAcc(0),
114 fhClustersDPhiAll(0),
115 fhClustersDThetaAll(0),
116 fhDPhiVsDThetaAll(0),
117 fhDPhiVsDThetaAcc(0),
118 fhetaTracklets(0),
119 fhphiTracklets(0),
120 fhetaClustersLay1(0),
121 fhphiClustersLay1(0){
122
123   fNFiredChips[0] = 0;
124   fNFiredChips[1] = 0;
125   // Method to reconstruct the charged particles multiplicity with the 
126   // SPD (tracklets).
127
128   SetHistOn();
129
130   if(AliITSReconstructor::GetRecoParam()) { 
131     SetPhiWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiWindow());
132     SetThetaWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterThetaWindow());
133     SetPhiShift(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiShift());
134     SetRemoveClustersFromOverlaps(AliITSReconstructor::GetRecoParam()->GetTrackleterRemoveClustersFromOverlaps());
135     SetPhiOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiOverlapCut());
136     SetZetaOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterZetaOverlapCut());
137   } else {
138     SetPhiWindow();
139     SetThetaWindow();
140     SetPhiShift();
141     SetRemoveClustersFromOverlaps();
142     SetPhiOverlapCut();
143     SetZetaOverlapCut();
144   } 
145   
146   fClustersLay1              = 0;
147   fClustersLay2              = 0;
148   fDetectorIndexClustersLay1 = 0;
149   fDetectorIndexClustersLay2 = 0;
150   fOverlapFlagClustersLay1   = 0;
151   fOverlapFlagClustersLay2   = 0;
152   fTracklets                 = 0;
153   fSClusters                 = 0;
154
155   // definition of histograms
156   Bool_t oldStatus = TH1::AddDirectoryStatus();
157   TH1::AddDirectory(kFALSE);
158   
159   fhClustersDPhiAcc   = new TH1F("dphiacc",  "dphi",  100,-0.1,0.1);
160   fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
161
162   fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,-0.1,0.1);
163
164   fhClustersDPhiAll   = new TH1F("dphiall",  "dphi",  100,0.0,0.5);
165   fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,0.0,0.5);
166
167   fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,0.,0.5,100,0.,0.5);
168
169   fhetaTracklets  = new TH1F("etaTracklets",  "eta",  100,-2.,2.);
170   fhphiTracklets  = new TH1F("phiTracklets",  "phi",  100, 0., 2*TMath::Pi());
171   fhetaClustersLay1  = new TH1F("etaClustersLay1",  "etaCl1",  100,-2.,2.);
172   fhphiClustersLay1  = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
173   
174   TH1::AddDirectory(oldStatus);
175 }
176
177 //______________________________________________________________________
178 AliITSMultReconstructor::AliITSMultReconstructor(const AliITSMultReconstructor &mr) : 
179 AliTrackleter(mr),
180 fDetTypeRec(0),fESDEvent(0),fTreeRP(0),fUsedClusLay1(0),fUsedClusLay2(0),
181 fClustersLay1(0),
182 fClustersLay2(0),
183 fDetectorIndexClustersLay1(0),
184 fDetectorIndexClustersLay2(0),
185 fOverlapFlagClustersLay1(0),
186 fOverlapFlagClustersLay2(0),
187 fTracklets(0),
188 fSClusters(0),
189 fNClustersLay1(0),
190 fNClustersLay2(0),
191 fNTracklets(0),
192 fNSingleCluster(0),
193 fPhiWindow(0),
194 fThetaWindow(0),
195 fPhiShift(0),
196 fRemoveClustersFromOverlaps(0),
197 fPhiOverlapCut(0),
198 fZetaOverlapCut(0),
199 fHistOn(0),
200 fhClustersDPhiAcc(0),
201 fhClustersDThetaAcc(0),
202 fhClustersDPhiAll(0),
203 fhClustersDThetaAll(0),
204 fhDPhiVsDThetaAll(0),
205 fhDPhiVsDThetaAcc(0),
206 fhetaTracklets(0),
207 fhphiTracklets(0),
208 fhetaClustersLay1(0),
209 fhphiClustersLay1(0)
210  {
211   // Copy constructor :!!! RS ATTENTION: old c-tor reassigned the pointers instead of creating a new copy -> would crash on delete
212    AliError("May not use");
213 }
214
215 //______________________________________________________________________
216 AliITSMultReconstructor& AliITSMultReconstructor::operator=(const AliITSMultReconstructor& mr){
217   // Assignment operator
218   if (this != &mr) {
219     this->~AliITSMultReconstructor();
220     new(this) AliITSMultReconstructor(mr);
221   }
222   return *this;
223 }
224
225 //______________________________________________________________________
226 AliITSMultReconstructor::~AliITSMultReconstructor(){
227   // Destructor
228
229   // delete histograms
230   delete fhClustersDPhiAcc;
231   delete fhClustersDThetaAcc;
232   delete fhClustersDPhiAll;
233   delete fhClustersDThetaAll;
234   delete fhDPhiVsDThetaAll;
235   delete fhDPhiVsDThetaAcc;
236   delete fhetaTracklets;
237   delete fhphiTracklets;
238   delete fhetaClustersLay1;
239   delete fhphiClustersLay1;
240   delete[] fUsedClusLay1;
241   delete[] fUsedClusLay2;
242   // delete arrays    
243   for(Int_t i=0; i<fNTracklets; i++)
244     delete [] fTracklets[i];
245     
246   for(Int_t i=0; i<fNSingleCluster; i++)
247     delete [] fSClusters[i];
248     
249   delete [] fClustersLay1;
250   delete [] fClustersLay2;
251   delete [] fDetectorIndexClustersLay1;
252   delete [] fDetectorIndexClustersLay2;
253   delete [] fOverlapFlagClustersLay1;
254   delete [] fOverlapFlagClustersLay2;
255   delete [] fTracklets;
256   delete [] fSClusters;
257 }
258
259 //____________________________________________________________________
260 void AliITSMultReconstructor::Reconstruct(AliESDEvent* esd, TTree* treeRP) 
261 {
262   //
263   // - calls LoadClusterArrays that finds the position of the clusters
264   //   (in global coord) 
265   // - convert the cluster coordinates to theta, phi (seen from the
266   //   interaction vertex). 
267   // - makes an array of tracklets 
268   //   
269   // After this method has been called, the clusters of the two layers
270   // and the tracklets can be retrieved by calling the Get'er methods.
271
272   if (!treeRP) { AliError(" Invalid ITS cluster tree !\n"); return; }
273   if (!esd) {AliError("ESDEvent is not available, use old reconstructor"); return;}
274   // reset counters
275   if (fMult) delete fMult; fMult = 0;
276   fNClustersLay1 = 0;
277   fNClustersLay2 = 0;
278   fNTracklets = 0; 
279   fNSingleCluster = 0;
280   //
281   fESDEvent = esd;
282   fTreeRP = treeRP;
283   //
284   // >>>> RS: this part is equivalent to former AliITSVertexer::FindMultiplicity
285   //
286   // see if there is a SPD vertex 
287   Bool_t isVtxOK=kTRUE, isCosmics=kFALSE;
288   AliESDVertex* vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexSPD();
289   if (!vtx && vtx->GetNContributors()<0) isVtxOK = kFALSE;
290   if (vtx && strstr(vtx->GetTitle(),"cosmics")) {
291     isVtxOK = kFALSE;
292     isCosmics = kTRUE;
293   }
294   //
295   if (!isVtxOK) {
296     if (!isCosmics) {
297       AliDebug(1,"Tracklets multiplicity not determined because the primary vertex was not found");
298       AliDebug(1,"Just counting the number of cluster-fired chips on the SPD layers");
299     }
300     vtx = 0;
301   }
302   float vtxf[3] = {vtx->GetX(),vtx->GetY(),vtx->GetZ()};
303   FindTracklets(vtxf);
304   //
305   CreateMultiplicityObject();
306 }
307
308 //____________________________________________________________________
309 void AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) {
310   //
311   // RS NOTE - this is old reconstructor invocation, to be used from VertexFinder
312   //
313   // - calls LoadClusterArray that finds the position of the clusters
314   //   (in global coord) 
315   // - convert the cluster coordinates to theta, phi (seen from the
316   //   interaction vertex). 
317   // - makes an array of tracklets 
318   //   
319   // After this method has been called, the clusters of the two layers
320   // and the tracklets can be retrieved by calling the Get'er methods.
321   if (fMult) delete fMult; fMult = 0;
322   fNClustersLay1 = 0;
323   fNClustersLay2 = 0;
324   fNTracklets = 0; 
325   fNSingleCluster = 0;
326   //
327   if (!clusterTree) { AliError(" Invalid ITS cluster tree !\n"); return; }
328   //
329   fESDEvent = 0;
330   fTreeRP = clusterTree;
331   //
332   FindTracklets(vtx);
333   //
334 }
335
336 //____________________________________________________________________
337 void AliITSMultReconstructor::FindTracklets(const Float_t *vtx) 
338 {
339   // Find tracklets converging to vertex
340   //
341   LoadClusterArrays(fTreeRP);
342   // flag clusters used by ESD tracks
343   if (fESDEvent) ProcessESDTracks();
344
345   if (!vtx) return;
346
347   const Double_t pi = TMath::Pi();
348   
349   // dPhi shift is field dependent
350   // get average magnetic field
351   Float_t bz = 0;
352   AliMagF* field = 0;
353   if (TGeoGlobalMagField::Instance()) field = dynamic_cast<AliMagF*>(TGeoGlobalMagField::Instance()->GetField());
354   if (!field)
355   {
356     AliError("Could not retrieve magnetic field. Assuming no field. Delta Phi shift will be deactivated in AliITSMultReconstructor.")
357   }
358   else
359     bz = TMath::Abs(field->SolenoidField());
360   
361   const Double_t dPhiShift = fPhiShift / 5 * bz; 
362   AliDebug(1, Form("Using phi shift of %f", dPhiShift));
363   
364   const Double_t dPhiWindow2 = fPhiWindow * fPhiWindow;
365   const Double_t dThetaWindow2 = fThetaWindow * fThetaWindow;
366   
367   Int_t* partners = new Int_t[fNClustersLay2];
368   Float_t* minDists = new Float_t[fNClustersLay2];
369   Int_t* associatedLay1 = new Int_t[fNClustersLay1];
370   TArrayI** blacklist = new TArrayI*[fNClustersLay1];
371
372   for (Int_t i=0; i<fNClustersLay2; i++) {
373     partners[i] = -1;
374     minDists[i] = 2;
375   }
376   for (Int_t i=0; i<fNClustersLay1; i++)
377     associatedLay1[i] = 0; 
378   for (Int_t i=0; i<fNClustersLay1; i++)
379     blacklist[i] = 0;
380
381   // find the tracklets
382   AliDebug(1,"Looking for tracklets... ");  
383   
384   //###########################################################
385   // Loop on layer 1 : finding theta, phi and z 
386   for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {    
387     float *clPar = GetClusterLayer1(iC1);
388     Float_t x = clPar[kClTh] - vtx[0];
389     Float_t y = clPar[kClPh] - vtx[1];
390     Float_t z = clPar[kClZ]  - vtx[2];
391
392     Float_t r    = TMath::Sqrt(x*x + y*y + z*z);
393     
394     clPar[kClTh] = TMath::ACos(z/r);                   // Store Theta
395     clPar[kClPh] = TMath::Pi() + TMath::ATan2(-y,-x);  // Store Phi
396     
397     if (fHistOn) {
398       Float_t eta = clPar[kClTh];
399       eta= TMath::Tan(eta/2.);
400       eta=-TMath::Log(eta);
401       fhetaClustersLay1->Fill(eta);    
402       fhphiClustersLay1->Fill(clPar[kClPh]);
403     }      
404   }
405   
406   // Loop on layer 2 : finding theta, phi and r   
407   for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {    
408     float *clPar = GetClusterLayer2(iC2);
409     Float_t x = clPar[kClTh] - vtx[0];
410     Float_t y = clPar[kClPh] - vtx[1];
411     Float_t z = clPar[kClZ]  - vtx[2];
412    
413     Float_t r    = TMath::Sqrt(x*x + y*y + z*z);
414
415     clPar[kClTh] = TMath::ACos(z/r);                   // Store Theta
416     clPar[kClPh] = TMath::Pi() + TMath::ATan2(-y,-x);  // Store Phi    
417   }  
418   
419   //###########################################################
420   Int_t found = 1;
421   while (found > 0) {
422     found = 0;
423
424     // Step1: find all tracklets allowing double assocation
425     // Loop on layer 1 
426     for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {    
427
428       // already used or in the overlap ?
429       if (associatedLay1[iC1] != 0 || fOverlapFlagClustersLay1[iC1]) continue;
430
431       found++;
432
433       // reset of variables for multiple candidates
434       Int_t  iC2WithBestDist = -1;   // reset
435       Double_t minDist       =  2;   // reset
436       float* clPar1 = GetClusterLayer1(iC1);
437
438       // Loop on layer 2 
439       for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {      
440
441         // in the overlap ?
442         if (fOverlapFlagClustersLay2[iC2]) continue;
443         float* clPar2 = GetClusterLayer2(iC2);
444
445         if (blacklist[iC1]) {
446           Bool_t blacklisted = kFALSE;
447           for (Int_t i=blacklist[iC1]->GetSize(); i--;) {
448             if (blacklist[iC1]->At(i) == iC2) {
449               blacklisted = kTRUE;
450               break;
451             }
452           }
453           if (blacklisted) continue;
454         }
455
456         // find the difference in angles
457         Double_t dTheta = TMath::Abs(clPar2[kClTh] - clPar1[kClTh]);
458         Double_t dPhi   = TMath::Abs(clPar2[kClPh] - clPar1[kClPh]);
459         // take into account boundary condition
460         if (dPhi>pi) dPhi=2.*pi-dPhi;
461         
462         if (fHistOn) {
463           fhClustersDPhiAll->Fill(dPhi);
464           fhClustersDThetaAll->Fill(dTheta);    
465           fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
466         }
467         
468         dPhi -= dPhiShift;
469                 
470         // make "elliptical" cut in Phi and Theta! 
471         Float_t d = dPhi*dPhi/dPhiWindow2 + dTheta*dTheta/dThetaWindow2;
472
473         // look for the minimum distance: the minimum is in iC2WithBestDist
474         if (d<1 && d<minDist) {
475           minDist=d;
476           iC2WithBestDist = iC2;
477         }
478       } // end of loop over clusters in layer 2 
479     
480       if (minDist<1) { // This means that a cluster in layer 2 was found that matches with iC1
481
482         if (minDists[iC2WithBestDist] > minDist) {
483           Int_t oldPartner = partners[iC2WithBestDist];
484           partners[iC2WithBestDist] = iC1;
485           minDists[iC2WithBestDist] = minDist;
486
487           // mark as assigned
488           associatedLay1[iC1] = 1;
489           
490           if (oldPartner != -1) {
491             // redo partner search for cluster in L0 (oldPartner), putting this one (iC2WithBestDist) on its blacklist
492             if (blacklist[oldPartner] == 0) {
493               blacklist[oldPartner] = new TArrayI(1);
494             } else blacklist[oldPartner]->Set(blacklist[oldPartner]->GetSize()+1);
495
496             blacklist[oldPartner]->AddAt(iC2WithBestDist, blacklist[oldPartner]->GetSize()-1);
497
498             // mark as free
499             associatedLay1[oldPartner] = 0;
500           }
501         } else {
502           // try again to find a cluster without considering iC2WithBestDist 
503           if (blacklist[iC1] == 0) {
504             blacklist[iC1] = new TArrayI(1);
505           }
506           else 
507             blacklist[iC1]->Set(blacklist[iC1]->GetSize()+1);
508    
509           blacklist[iC1]->AddAt(iC2WithBestDist, blacklist[iC1]->GetSize()-1);
510         } 
511
512       } else // cluster has no partner; remove
513         associatedLay1[iC1] = 2;   
514     } // end of loop over clusters in layer 1
515   }  
516  
517   // Step2: store tracklets; remove used clusters
518   for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {
519
520     if (partners[iC2] == -1) continue;
521
522     if (fRemoveClustersFromOverlaps) FlagClustersInOverlapRegions (partners[iC2],iC2);
523
524
525     if (fOverlapFlagClustersLay1[partners[iC2]] || fOverlapFlagClustersLay2[iC2]) continue;
526
527     float* clPar2 = GetClusterLayer2(iC2);
528     float* clPar1 = GetClusterLayer1(partners[iC2]);
529
530     Float_t* tracklet = fTracklets[fNTracklets] = new Float_t[kTrNPar]; // RS Add also the cluster id's
531   
532     // use the theta from the clusters in the first layer
533     tracklet[kTrTheta] = clPar1[kClTh];
534     // use the phi from the clusters in the first layer
535     tracklet[kTrPhi] = clPar1[kClPh];
536     // store the difference between phi1 and phi2
537     tracklet[kTrDPhi] = clPar1[kClPh] - clPar2[kClPh];
538
539     // define dphi in the range [0,pi] with proper sign (track charge correlated)
540     if (tracklet[kTrDPhi] > TMath::Pi())   tracklet[kTrDPhi] = tracklet[kTrDPhi]-2.*TMath::Pi();
541     if (tracklet[kTrDPhi] < -TMath::Pi())  tracklet[kTrDPhi] = tracklet[kTrDPhi]+2.*TMath::Pi();
542
543     // store the difference between theta1 and theta2
544     tracklet[kTrDTheta] = clPar1[kClTh] - clPar2[kClTh];
545
546     if (fHistOn) {
547       fhClustersDPhiAcc->Fill(tracklet[kTrDPhi]); 
548       fhClustersDThetaAcc->Fill(tracklet[kTrDTheta]);    
549       fhDPhiVsDThetaAcc->Fill(tracklet[kTrDTheta],tracklet[kTrDPhi]);
550     }
551
552     // find label
553     // if equal label in both clusters found this label is assigned
554     // if no equal label can be found the first labels of the L1 AND L2 cluster are assigned
555     Int_t label1 = 0;
556     Int_t label2 = 0;
557     while (label2 < 3) {
558       if ((Int_t) clPar1[kClMC0+label1] != -2 && (Int_t) clPar1[kClMC0+label1] == (Int_t) clPar2[kClMC0+label2])
559         break;
560       label1++;
561       if (label1 == 3) {
562         label1 = 0;
563         label2++;
564       }
565     }
566     if (label2 < 3) {
567       AliDebug(AliLog::kDebug, Form("Found label %d == %d for tracklet candidate %d\n", (Int_t) clPar1[kClMC0+label1], (Int_t) clPar1[kClMC0+label2], fNTracklets));
568       tracklet[kTrLab1] = clPar1[kClMC0+label1];
569       tracklet[kTrLab2] = clPar2[kClMC0+label2];
570     } else {
571       AliDebug(AliLog::kDebug, Form("Did not find label %d %d %d %d %d %d for tracklet candidate %d\n", (Int_t) clPar1[kClMC0], (Int_t) clPar1[kClMC1], (Int_t) clPar1[kClMC2], (Int_t) clPar2[kClMC0], (Int_t) clPar2[kClMC1], (Int_t) clPar2[kClMC2], fNTracklets));
572       tracklet[kTrLab1] = clPar1[kClMC0];
573       tracklet[kTrLab2] = clPar2[kClMC0];
574     }
575
576     if (fHistOn) {
577       Float_t eta = tracklet[kTrTheta];
578       eta= TMath::Tan(eta/2.);
579       eta=-TMath::Log(eta);
580       fhetaTracklets->Fill(eta);
581       fhphiTracklets->Fill(tracklet[kTrPhi]);
582     }
583     //
584     tracklet[kClID1] = partners[iC2];
585     tracklet[kClID2] = iC2;
586     //
587     AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
588     AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", partners[iC2], iC2));
589     fNTracklets++;
590
591     associatedLay1[partners[iC2]] = 1;
592   }
593   
594   // Delete the following else if you do not want to save Clusters! 
595   // store the cluster
596   for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {
597
598     float* clPar1 = GetClusterLayer1(iC1);
599
600     if (associatedLay1[iC1]==2||associatedLay1[iC1]==0) { 
601       fSClusters[fNSingleCluster] = new Float_t[kClNPar];
602       fSClusters[fNSingleCluster][kSCTh] = clPar1[kClTh];
603       fSClusters[fNSingleCluster][kSCPh] = clPar1[kClPh];
604       fSClusters[fNSingleCluster][kSCID] = iC1;
605       AliDebug(1,Form(" Adding a single cluster %d (cluster %d  of layer 1)",
606                 fNSingleCluster, iC1));
607       fNSingleCluster++;
608     }
609   }
610
611   delete[] partners;
612   delete[] minDists;
613
614   for (Int_t i=0; i<fNClustersLay1; i++)
615     if (blacklist[i])
616       delete blacklist[i];
617   delete[] blacklist;
618
619   AliDebug(1,Form("%d tracklets found", fNTracklets));
620 }
621
622 //____________________________________________________________________
623 void AliITSMultReconstructor::CreateMultiplicityObject()
624 {
625   // create AliMultiplicity object and store it in the ESD event
626   //
627   TBits fastOrFiredMap,firedChipMap;
628   if (fDetTypeRec) {
629    fastOrFiredMap  = fDetTypeRec->GetFastOrFiredMap();
630    firedChipMap    = fDetTypeRec->GetFiredChipMap(fTreeRP);
631   }
632   //
633   fMult = new AliMultiplicity(fNTracklets,fNSingleCluster,fNFiredChips[0],fNFiredChips[1],fastOrFiredMap);
634   fMult->SetFiredChipMap(firedChipMap);
635   AliITSRecPointContainer* rcont = AliITSRecPointContainer::Instance();
636   fMult->SetITSClusters(0,rcont->GetNClustersInLayer(1,fTreeRP));
637   for(Int_t kk=2;kk<=6;kk++) fMult->SetITSClusters(kk-1,rcont->GetNClustersInLayerFast(kk));
638   //
639   for (int i=fNTracklets;i--;)  {
640     float* tlInfo = fTracklets[i];
641     fMult->SetTrackletData(i,tlInfo, fUsedClusLay1[int(tlInfo[kClID1])]|fUsedClusLay2[int(tlInfo[kClID2])]);
642   }
643   //  
644   for (int i=fNSingleCluster;i--;) {
645     float* clInfo = fSClusters[i];
646     fMult->SetSingleClusterData(i,clInfo,fUsedClusLay1[int(clInfo[kSCID])]);
647   }
648   fMult->CompactBits();
649   //
650 }
651
652
653 //____________________________________________________________________
654 void AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) 
655 {
656   // This method
657   // - gets the clusters from the cluster tree 
658   // - convert them into global coordinates 
659   // - store them in the internal arrays
660   // - count the number of cluster-fired chips
661   //
662   // RS: This method was strongly modified wrt original by Jan Fiete. In order to have the same numbering
663   // of clusters as in the ITS reco I had to introduce sorting in Z
664   // Also note that now the clusters data are stored not in float[6] attached to float**, but in 1-D array
665   
666   AliDebug(1,"Loading clusters and cluster-fired chips ...");
667   
668   fNClustersLay1 = 0;
669   fNClustersLay2 = 0;
670   fNFiredChips[0] = 0;
671   fNFiredChips[1] = 0;
672   
673   AliITSsegmentationSPD seg;
674
675   AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
676   TClonesArray* itsClusters=rpcont->FetchClusters(0,itsClusterTree);
677   if(!rpcont->IsSPDActive()){
678     AliWarning("No SPD rec points found, multiplicity not calculated");
679     return;
680   } 
681   //
682   // count clusters
683   // loop over the SPD subdetectors
684   TObjArray clArr(100);
685   for (int il=0;il<2;il++) {
686     int nclLayer = 0;
687     int detMin = AliITSgeomTGeo::GetModuleIndex(il+1,1,1);
688     int detMax = AliITSgeomTGeo::GetModuleIndex(il+2,1,1);
689     for (int idt=detMin;idt<detMax;idt++) {
690       itsClusters=rpcont->UncheckedGetClusters(idt);
691       int nClusters = itsClusters->GetEntriesFast();
692       if (!nClusters) continue;
693       Int_t nClustersInChip[5] = {0,0,0,0,0};
694       while(nClusters--) {
695         AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
696         if (!cluster) continue;
697         clArr[nclLayer++] = cluster;
698         nClustersInChip[ seg.GetChipFromLocal(0,cluster->GetDetLocalZ()) ]++; 
699       }
700       for(Int_t ifChip=5;ifChip--;) if (nClustersInChip[ifChip]) fNFiredChips[il]++;
701     }
702     // sort the clusters in Z (to have the same numbering as in ITS reco
703     Float_t *z    = new Float_t[nclLayer];
704     Int_t * index = new Int_t[nclLayer];
705     for (int ic=0;ic<nclLayer;ic++) z[ic] = ((AliITSRecPoint*)clArr[ic])->GetZ();
706     TMath::Sort(nclLayer,z,index,kFALSE);
707     Float_t*   clustersLay              = new Float_t[nclLayer*kClNPar];
708     Int_t*     detectorIndexClustersLay = new Int_t[nclLayer];
709     Bool_t*    overlapFlagClustersLay   = new Bool_t[nclLayer];
710     Char_t*    usedClusLay              = new Char_t[nclLayer];
711     //
712     for (int ic=0;ic<nclLayer;ic++) {
713       AliITSRecPoint* cluster = (AliITSRecPoint*)clArr[index[ic]];
714       float* clPar = &clustersLay[ic*kClNPar];
715       //      
716       cluster->GetGlobalXYZ( clPar );
717       detectorIndexClustersLay[ic] = cluster->GetDetectorIndex(); 
718       overlapFlagClustersLay[ic]   = kFALSE;
719       usedClusLay[ic]              = 0;
720       for (Int_t i=3;i--;) clPar[kClMC0+i] = cluster->GetLabel(i);
721     }
722     clArr.Clear();
723     delete[] z;
724     delete[] index;
725     //
726     if (il==0) {
727       fClustersLay1              = clustersLay;
728       fOverlapFlagClustersLay1   = overlapFlagClustersLay;
729       fDetectorIndexClustersLay1 = detectorIndexClustersLay;
730       fUsedClusLay1              = usedClusLay;
731       fNClustersLay1             = nclLayer;
732     }
733     else {
734       fClustersLay2              = clustersLay;
735       fOverlapFlagClustersLay2   = overlapFlagClustersLay;
736       fDetectorIndexClustersLay2 = detectorIndexClustersLay;
737       fUsedClusLay2              = usedClusLay;
738       fNClustersLay2             = nclLayer;
739     }
740   }
741   //
742   // no double association allowed
743   int nmaxT                  = TMath::Min(fNClustersLay1, fNClustersLay2);
744   fTracklets                 = new Float_t*[nmaxT];
745   fSClusters                 = new Float_t*[fNClustersLay1]; 
746   for (Int_t i=nmaxT;i--;) fTracklets[i] = 0;
747   //
748   AliDebug(1,Form("(clusters in layer 1 : %d,  layer 2: %d)",fNClustersLay1,fNClustersLay2));
749   AliDebug(1,Form("(cluster-fired chips in layer 1 : %d,  layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
750 }
751 //____________________________________________________________________
752 void
753 AliITSMultReconstructor::LoadClusterFiredChips(TTree* itsClusterTree) {
754   // This method
755   // - gets the clusters from the cluster tree 
756   // - counts the number of (cluster)fired chips
757   
758   AliDebug(1,"Loading cluster-fired chips ...");
759   
760   fNFiredChips[0] = 0;
761   fNFiredChips[1] = 0;
762   
763   AliITSsegmentationSPD seg;
764   AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
765   TClonesArray* itsClusters=rpcont->FetchClusters(0,itsClusterTree);
766   if(!rpcont->IsSPDActive()){
767     AliWarning("No SPD rec points found, multiplicity not calculated");
768     return;
769   } 
770
771   // loop over the its subdetectors
772   Int_t nSPDmodules=AliITSgeomTGeo::GetModuleIndex(3,1,1);
773   for (Int_t iIts=0; iIts < nSPDmodules; iIts++) {
774     itsClusters=rpcont->UncheckedGetClusters(iIts);
775     Int_t nClusters = itsClusters->GetEntriesFast();
776
777     // number of clusters in each chip of the current module
778     Int_t nClustersInChip[5] = {0,0,0,0,0};
779     Int_t layer = 0;
780     
781     // loop over clusters
782     while(nClusters--) {
783       AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
784       
785       layer = cluster->GetLayer();
786       if (layer>1) continue;            
787
788       // find the chip for the current cluster
789       Float_t locz = cluster->GetDetLocalZ();
790       Int_t iChip = seg.GetChipFromLocal(0,locz);
791       nClustersInChip[iChip]++; 
792       
793     }// end of cluster loop
794
795     // get number of fired chips in the current module
796     for(Int_t ifChip=0; ifChip<5; ifChip++) {
797       if(nClustersInChip[ifChip] >= 1)  fNFiredChips[layer]++;
798     }
799
800   } // end of its "subdetector" loop  
801   
802
803   AliDebug(1,Form("(cluster-fired chips in layer 1 : %d,  layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
804 }
805 //____________________________________________________________________
806 void
807 AliITSMultReconstructor::SaveHists() {
808   // This method save the histograms on the output file
809   // (only if fHistOn is TRUE). 
810   
811   if (!fHistOn)
812     return;
813
814   fhClustersDPhiAll->Write();
815   fhClustersDThetaAll->Write();
816   fhDPhiVsDThetaAll->Write();
817
818   fhClustersDPhiAcc->Write();
819   fhClustersDThetaAcc->Write();
820   fhDPhiVsDThetaAcc->Write();
821
822   fhetaTracklets->Write();
823   fhphiTracklets->Write();
824   fhetaClustersLay1->Write();
825   fhphiClustersLay1->Write();
826 }
827
828 //____________________________________________________________________
829 void 
830 AliITSMultReconstructor::FlagClustersInOverlapRegions (Int_t iC1, Int_t iC2WithBestDist) {
831
832   Float_t distClSameMod=0.;
833   Float_t distClSameModMin=0.;
834   Int_t   iClOverlap =0;
835   Float_t meanRadiusLay1 = 3.99335; // average radius inner layer
836   Float_t meanRadiusLay2 = 7.37935; // average radius outer layer;
837
838   Float_t zproj1=0.;
839   Float_t zproj2=0.;
840   Float_t deZproj=0.;
841   Float_t* clPar1  = GetClusterLayer1(iC1);
842   Float_t* clPar2B = GetClusterLayer2(iC2WithBestDist);
843   // Loop on inner layer clusters
844   for (Int_t iiC1=0; iiC1<fNClustersLay1; iiC1++) {
845     if (!fOverlapFlagClustersLay1[iiC1]) {
846       // only for adjacent modules
847       if ((TMath::Abs(fDetectorIndexClustersLay1[iC1]-fDetectorIndexClustersLay1[iiC1])==4)||
848          (TMath::Abs(fDetectorIndexClustersLay1[iC1]-fDetectorIndexClustersLay1[iiC1])==76)) {
849         Float_t *clPar11 = GetClusterLayer1(iiC1);
850         Float_t dePhi=TMath::Abs(clPar11[kClPh]-clPar1[kClPh]);
851         if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
852
853         zproj1=meanRadiusLay1/TMath::Tan(clPar1[kClTh]);
854         zproj2=meanRadiusLay1/TMath::Tan(clPar11[kClTh]);
855
856         deZproj=TMath::Abs(zproj1-zproj2);
857
858         distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
859         if (distClSameMod<=1.) fOverlapFlagClustersLay1[iiC1]=kTRUE;
860
861 //        if (distClSameMod<=1.) {
862 //          if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
863 //            distClSameModMin=distClSameMod;
864 //            iClOverlap=iiC1;
865 //          } 
866 //        }
867
868
869       } // end adjacent modules
870     } 
871   } // end Loop on inner layer clusters
872
873 //  if (distClSameModMin!=0.) fOverlapFlagClustersLay1[iClOverlap]=kTRUE;
874
875   distClSameMod=0.;
876   distClSameModMin=0.;
877   iClOverlap =0;
878   // Loop on outer layer clusters
879   for (Int_t iiC2=0; iiC2<fNClustersLay2; iiC2++) {
880     if (!fOverlapFlagClustersLay2[iiC2]) {
881       // only for adjacent modules
882       Float_t *clPar2 = GetClusterLayer2(iiC2);
883       if ((TMath::Abs(fDetectorIndexClustersLay2[iC2WithBestDist]-fDetectorIndexClustersLay2[iiC2])==4) ||
884          (TMath::Abs(fDetectorIndexClustersLay2[iC2WithBestDist]-fDetectorIndexClustersLay2[iiC2])==156)) {
885         Float_t dePhi=TMath::Abs(clPar2[kClPh]-clPar2B[kClPh]);
886         if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
887
888         zproj1=meanRadiusLay2/TMath::Tan(clPar2B[kClTh]);
889         zproj2=meanRadiusLay2/TMath::Tan(clPar2[kClTh]);
890
891         deZproj=TMath::Abs(zproj1-zproj2);
892         distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
893         if (distClSameMod<=1.) fOverlapFlagClustersLay2[iiC2]=kTRUE;
894
895 //        if (distClSameMod<=1.) {
896 //          if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
897 //            distClSameModMin=distClSameMod;
898 //            iClOverlap=iiC2;
899 //          }
900 //        }
901
902       } // end adjacent modules
903     }
904   } // end Loop on outer layer clusters
905
906 //  if (distClSameModMin!=0.) fOverlapFlagClustersLay2[iClOverlap]=kTRUE;
907
908 }
909
910 //____________________________________________________________________
911 void AliITSMultReconstructor::ProcessESDTracks()
912 {
913   // Flag the clusters used by ESD tracks
914   // Flag primary tracks to be used for multiplicity counting 
915   //
916   if (!fESDEvent) return;
917   AliESDVertex* vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexTracks();
918   if (!vtx) vtx = (AliESDVertex*)fESDEvent->GetPrimaryVertexSPD();
919   if (!vtx) {
920     AliDebug(1,"No primary vertex: cannot flag primary tracks");
921     return;
922   }
923   Int_t ntracks = fESDEvent->GetNumberOfTracks();
924   for(Int_t itr=0; itr<ntracks; itr++) {
925     AliESDtrack* track = fESDEvent->GetTrack(itr);
926     if (!track->IsOn(AliESDtrack::kITSin)) continue; // use only tracks propagated in ITS to vtx
927     FlagTrackClusters(track);
928     FlagIfPrimary(track,vtx);
929   }
930   //
931 }
932
933 //____________________________________________________________________
934 void AliITSMultReconstructor::FlagTrackClusters(const AliESDtrack* track)
935 {
936   // RS: flag the SPD clusters of the track if it is useful for the multiplicity estimation
937   //
938   Int_t idx[12];
939   if ( track->GetITSclusters(idx)<3 ) return; // at least 3 clusters must be used in the fit
940   //
941   char mark = track->IsOn(AliESDtrack::kITSpureSA) ? kITSSAPBit : kITSTPCBit;
942   char *uClus[2] = {fUsedClusLay1,fUsedClusLay2};
943   for (int i=AliESDfriendTrack::kMaxITScluster;i--;) {
944     // note: i>=6 is for extra clusters
945     if (idx[i]<0) continue;
946     int layID= (idx[i] & 0xf0000000) >> 28; 
947     if (layID>1) continue; // SPD only
948     int clID = (idx[i] & 0x0fffffff);
949     uClus[layID][clID] |= mark;
950   }
951   //
952 }
953
954 //____________________________________________________________________
955 void AliITSMultReconstructor::FlagIfPrimary(AliESDtrack* track, const AliVertex* vtx)
956 {
957   // RS: check if the track is primary and set the flag
958   const double kPDCASPD1 = 0.1;
959   const double kPDCASPD0 = 0.3;
960   //
961   double cut = (track->HasPointOnITSLayer(0)||track->HasPointOnITSLayer(1)) ? kPDCASPD1 : kPDCASPD0;
962   // in principle, the track must already have been propagated to vertex
963   /*
964   Double_t dzRec[2]={0,0}, covdzRec[3];
965   track->PropagateToDCA(vtx, fESDEvent->GetMagneticField(), 3.0, dzRec, covdzRec);
966   */
967   double dist = track->GetD(vtx->GetX(),vtx->GetY(),fESDEvent->GetMagneticField());
968   if (TMath::Abs(dist*track->P())<cut) track->SetStatus(AliESDtrack::kMultPrimary);
969 }