- include option to avoid multiple reconstruction of tracklets in the
[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 /* $Id$ */
17
18 //____________________________________________________________________
19 // 
20 // AliITSMultReconstructor - find clusters in the pixels (theta and
21 // phi) and tracklets.
22 // 
23 // These can be used to extract charged particles multiplcicity from the ITS.
24 //
25 // A tracklet consist of two ITS clusters, one in the first pixel
26 // layer and one in the second. The clusters are associates if the 
27 // differencies in Phi (azimuth) and Zeta (longitudinal) are inside 
28 // a fiducial volume. In case of multiple candidates it is selected the
29 // candidate with minimum distance in Phi. 
30 // The parameter AssociationChoice allows to control if two clusters
31 // in layer 2 can be associated to the same cluster in layer 1 or not.
32 // (TRUE means double associations exluded; default = TRUE)
33 //
34 // Two methods return the number of traklets and the number of clusters 
35 // in the first SPD layer (GetNTracklets GetNSingleClusters)
36 //
37 // -----------------------------------------------------------------
38 // 
39 // NOTE: The cuts on phi and zeta depend on the interacting system (p-p  
40 //  or Pb-Pb). Please, check the file AliITSMultReconstructor.h and be 
41 //  sure that SetPhiWindow and SetZetaWindow are defined accordingly.
42 // 
43 //  Author :  Tiziano Virgili 
44 //  
45 //  Recent updates (D. Elia, INFN Bari):
46 //     - multiple association forbidden (fOnlyOneTrackletPerC2 = kTRUE) 
47 //     - phi definition changed to ALICE convention (0,2*TMath::pi()) 
48 //     - cluster coordinates taken with GetGlobalXYZ()
49 //     - fGeometry removed
50 //     - number of fired chips on the two layers
51 //     - option to avoid duplicates in the overlaps (fRemoveClustersFromOverlaps)
52 //     - options and fiducial cuts via AliITSRecoParam
53 //
54 //____________________________________________________________________
55
56 #include <TClonesArray.h>
57 #include <TH1F.h>
58 #include <TH2F.h>
59 #include <TTree.h>
60
61 #include "AliITSMultReconstructor.h"
62 #include "AliITSReconstructor.h"
63 #include "AliITSsegmentationSPD.h"
64 #include "AliITSRecPoint.h"
65 #include "AliITSgeom.h"
66 #include "AliLog.h"
67
68 //____________________________________________________________________
69 ClassImp(AliITSMultReconstructor)
70
71
72 //____________________________________________________________________
73 AliITSMultReconstructor::AliITSMultReconstructor():
74 fClustersLay1(0),
75 fClustersLay2(0),
76 fDetectorIndexClustersLay1(0),
77 fDetectorIndexClustersLay2(0),
78 fOverlapFlagClustersLay1(0),
79 fOverlapFlagClustersLay2(0),
80 fTracklets(0),
81 fSClusters(0),
82 fAssociationFlag(0),
83 fNClustersLay1(0),
84 fNClustersLay2(0),
85 fNTracklets(0),
86 fNSingleCluster(0),
87 fOnlyOneTrackletPerC2(0),
88 fPhiWindow(0),
89 fZetaWindow(0),
90 fRemoveClustersFromOverlaps(0),
91 fPhiOverlapCut(0),
92 fZetaOverlapCut(0),
93 fHistOn(0),
94 fhClustersDPhiAcc(0),
95 fhClustersDThetaAcc(0),
96 fhClustersDZetaAcc(0),
97 fhClustersDPhiAll(0),
98 fhClustersDThetaAll(0),
99 fhClustersDZetaAll(0),
100 fhDPhiVsDThetaAll(0),
101 fhDPhiVsDThetaAcc(0),
102 fhDPhiVsDZetaAll(0),
103 fhDPhiVsDZetaAcc(0),
104 fhetaTracklets(0),
105 fhphiTracklets(0),
106 fhetaClustersLay1(0),
107 fhphiClustersLay1(0){
108
109   fNFiredChips[0] = 0;
110   fNFiredChips[1] = 0;
111
112   // Method to reconstruct the charged particles multiplicity with the 
113   // SPD (tracklets).
114
115
116   SetHistOn();
117
118   if(AliITSReconstructor::GetRecoParam()) { 
119     SetOnlyOneTrackletPerC2(AliITSReconstructor::GetRecoParam()->GetTrackleterOnlyOneTrackletPerC2());
120     SetPhiWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiWindow());
121     SetZetaWindow(AliITSReconstructor::GetRecoParam()->GetTrackleterZetaWindow());
122     SetRemoveClustersFromOverlaps(AliITSReconstructor::GetRecoParam()->GetTrackleterRemoveClustersFromOverlaps());
123     SetPhiOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterPhiOverlapCut());
124     SetZetaOverlapCut(AliITSReconstructor::GetRecoParam()->GetTrackleterZetaOverlapCut());
125   } else {
126     SetOnlyOneTrackletPerC2();
127     SetPhiWindow();
128     SetZetaWindow();
129     SetRemoveClustersFromOverlaps();
130     SetPhiOverlapCut();
131     SetZetaOverlapCut();
132   } 
133   
134
135   fClustersLay1              = new Float_t*[300000];
136   fClustersLay2              = new Float_t*[300000];
137   fDetectorIndexClustersLay1 = new Int_t[300000];
138   fDetectorIndexClustersLay2 = new Int_t[300000];
139   fOverlapFlagClustersLay1   = new Bool_t[300000];
140   fOverlapFlagClustersLay2   = new Bool_t[300000];
141   fTracklets                 = new Float_t*[300000];
142   fSClusters                 = new Float_t*[300000];
143   fAssociationFlag           = new Bool_t[300000];
144
145   for(Int_t i=0; i<300000; i++) {
146     fClustersLay1[i]       = new Float_t[6];
147     fClustersLay2[i]       = new Float_t[6];
148     fTracklets[i]          = new Float_t[5];
149     fSClusters[i]           = new Float_t[2];
150     fOverlapFlagClustersLay1[i]   = kFALSE;
151     fOverlapFlagClustersLay2[i]   = kFALSE;
152     fAssociationFlag[i]    = kFALSE;
153   }
154
155   // definition of histograms
156   fhClustersDPhiAcc   = new TH1F("dphiacc",  "dphi",  100,0.,0.1);
157   fhClustersDPhiAcc->SetDirectory(0);
158   fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
159   fhClustersDThetaAcc->SetDirectory(0);
160   fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
161   fhClustersDZetaAcc->SetDirectory(0);
162
163   fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,0.,0.1);
164   fhDPhiVsDZetaAcc->SetDirectory(0);
165   fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,0.,0.1);
166   fhDPhiVsDThetaAcc->SetDirectory(0);
167
168   fhClustersDPhiAll   = new TH1F("dphiall",  "dphi",  100,0.0,0.5);
169   fhClustersDPhiAll->SetDirectory(0);
170   fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
171   fhClustersDThetaAll->SetDirectory(0);
172   fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
173   fhClustersDZetaAll->SetDirectory(0);
174
175   fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,0.,0.5);
176   fhDPhiVsDZetaAll->SetDirectory(0);
177   fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,0.,0.5);
178   fhDPhiVsDThetaAll->SetDirectory(0);
179
180   fhetaTracklets  = new TH1F("etaTracklets",  "eta",  100,-2.,2.);
181   fhetaTracklets->SetDirectory(0);
182   fhphiTracklets  = new TH1F("phiTracklets",  "phi",  100, 0., 2*TMath::Pi());
183   fhphiTracklets->SetDirectory(0);
184   fhetaClustersLay1  = new TH1F("etaClustersLay1",  "etaCl1",  100,-2.,2.);
185   fhetaClustersLay1->SetDirectory(0);
186   fhphiClustersLay1  = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
187   fhphiClustersLay1->SetDirectory(0);
188 }
189
190 //______________________________________________________________________
191 AliITSMultReconstructor::AliITSMultReconstructor(const AliITSMultReconstructor &mr) : TObject(mr),
192 fClustersLay1(mr.fClustersLay1),
193 fClustersLay2(mr.fClustersLay2),
194 fDetectorIndexClustersLay1(mr.fDetectorIndexClustersLay1),
195 fDetectorIndexClustersLay2(mr.fDetectorIndexClustersLay2),
196 fOverlapFlagClustersLay1(mr.fOverlapFlagClustersLay1),
197 fOverlapFlagClustersLay2(mr.fOverlapFlagClustersLay2),
198 fTracklets(mr.fTracklets),
199 fSClusters(mr.fSClusters),
200 fAssociationFlag(mr.fAssociationFlag),
201 fNClustersLay1(mr.fNClustersLay1),
202 fNClustersLay2(mr.fNClustersLay2),
203 fNTracklets(mr.fNTracklets),
204 fNSingleCluster(mr.fNSingleCluster),
205 fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2),
206 fPhiWindow(mr.fPhiWindow),
207 fZetaWindow(mr.fZetaWindow),
208 fRemoveClustersFromOverlaps(mr.fRemoveClustersFromOverlaps),
209 fPhiOverlapCut(mr.fPhiOverlapCut),
210 fZetaOverlapCut(mr.fZetaOverlapCut),
211 fHistOn(mr.fHistOn),
212 fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
213 fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
214 fhClustersDZetaAcc(mr.fhClustersDZetaAcc),
215 fhClustersDPhiAll(mr.fhClustersDPhiAll),
216 fhClustersDThetaAll(mr.fhClustersDThetaAll),
217 fhClustersDZetaAll(mr.fhClustersDZetaAll),
218 fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
219 fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
220 fhDPhiVsDZetaAll(mr.fhDPhiVsDZetaAll),
221 fhDPhiVsDZetaAcc(mr.fhDPhiVsDZetaAcc),
222 fhetaTracklets(mr.fhetaTracklets),
223 fhphiTracklets(mr.fhphiTracklets),
224 fhetaClustersLay1(mr.fhetaClustersLay1),
225 fhphiClustersLay1(mr.fhphiClustersLay1) {
226   // Copy constructor
227
228 }
229
230 //______________________________________________________________________
231 AliITSMultReconstructor& AliITSMultReconstructor::operator=(const AliITSMultReconstructor& mr){
232   // Assignment operator
233   this->~AliITSMultReconstructor();
234   new(this) AliITSMultReconstructor(mr);
235   return *this;
236 }
237
238 //______________________________________________________________________
239 AliITSMultReconstructor::~AliITSMultReconstructor(){
240   // Destructor
241
242   // delete histograms
243   delete fhClustersDPhiAcc;
244   delete fhClustersDThetaAcc;
245   delete fhClustersDZetaAcc;
246   delete fhClustersDPhiAll;
247   delete fhClustersDThetaAll;
248   delete fhClustersDZetaAll;
249   delete fhDPhiVsDThetaAll;
250   delete fhDPhiVsDThetaAcc;
251   delete fhDPhiVsDZetaAll;
252   delete fhDPhiVsDZetaAcc;
253   delete fhetaTracklets;
254   delete fhphiTracklets;
255   delete fhetaClustersLay1;
256   delete fhphiClustersLay1;
257
258   // delete arrays
259   for(Int_t i=0; i<300000; i++) {
260     delete [] fClustersLay1[i];
261     delete [] fClustersLay2[i];
262     delete [] fTracklets[i];
263     delete [] fSClusters[i];
264   }
265   delete [] fClustersLay1;
266   delete [] fClustersLay2;
267   delete [] fDetectorIndexClustersLay1;
268   delete [] fDetectorIndexClustersLay2;
269   delete [] fOverlapFlagClustersLay1;
270   delete [] fOverlapFlagClustersLay2;
271   delete [] fTracklets;
272   delete [] fSClusters;
273
274   delete [] fAssociationFlag;
275 }
276
277 //____________________________________________________________________
278 void
279 AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) {
280   //
281   // - calls LoadClusterArray that finds the position of the clusters
282   //   (in global coord) 
283   // - convert the cluster coordinates to theta, phi (seen from the
284   //   interaction vertex). The third coordinate is used for ....
285   // - makes an array of tracklets 
286   //   
287   // After this method has been called, the clusters of the two layers
288   // and the tracklets can be retrieved by calling the Get'er methods.
289
290   // reset counters
291   fNClustersLay1 = 0;
292   fNClustersLay2 = 0;
293   fNTracklets = 0; 
294   fNSingleCluster = 0; 
295   // loading the clusters 
296   LoadClusterArrays(clusterTree);
297
298   // find the tracklets
299   AliDebug(1,"Looking for tracklets... ");  
300
301   //###########################################################
302   // Loop on layer 1 : finding theta, phi and z 
303   for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {    
304     Float_t x = fClustersLay1[iC1][0] - vtx[0];
305     Float_t y = fClustersLay1[iC1][1] - vtx[1];
306     Float_t z = fClustersLay1[iC1][2] - vtx[2];
307
308     Float_t r    = TMath::Sqrt(TMath::Power(x,2) +
309                                TMath::Power(y,2) +
310                                TMath::Power(z,2));
311     
312     fClustersLay1[iC1][0] = TMath::ACos(z/r);                   // Store Theta
313     fClustersLay1[iC1][1] = TMath::Pi() + TMath::ATan2(-y,-x);  // Store Phi
314     fClustersLay1[iC1][2] = z/r;                                // Store scaled z
315     if (fHistOn) {
316       Float_t eta=fClustersLay1[iC1][0];
317       eta= TMath::Tan(eta/2.);
318       eta=-TMath::Log(eta);
319       fhetaClustersLay1->Fill(eta);    
320       fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
321     }      
322   }
323   
324   // Loop on layer 2 : finding theta, phi and r   
325   for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {    
326     Float_t x = fClustersLay2[iC2][0] - vtx[0];
327     Float_t y = fClustersLay2[iC2][1] - vtx[1];
328     Float_t z = fClustersLay2[iC2][2] - vtx[2];
329    
330     Float_t r    = TMath::Sqrt(TMath::Power(x,2) +
331                                TMath::Power(y,2) +
332                                TMath::Power(z,2));
333     
334     fClustersLay2[iC2][0] = TMath::ACos(z/r);                   // Store Theta
335     fClustersLay2[iC2][1] = TMath::Pi() + TMath::ATan2(-y,-x);  // Store Phi
336     fClustersLay2[iC2][2] = z;                                  // Store z
337
338  // this only needs to be initialized for the fNClustersLay2 first associations
339     fAssociationFlag[iC2] = kFALSE;
340   }  
341   
342   //###########################################################
343   // Loop on layer 1 
344   for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {    
345
346     // reset of variables for multiple candidates
347     Int_t  iC2WithBestDist = 0;     // reset 
348     Float_t distmin        = 100.;  // just to put a huge number! 
349     Float_t dPhimin        = 0.;  // Used for histograms only! 
350     Float_t dThetamin      = 0.;  // Used for histograms only! 
351     Float_t dZetamin       = 0.;  // Used for histograms only! 
352      
353     if (fOverlapFlagClustersLay1[iC1]) continue;
354  
355     // Loop on layer 2 
356     for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {      
357       if (fOverlapFlagClustersLay2[iC2]) continue;
358       // The following excludes double associations
359       if (!fAssociationFlag[iC2]) {
360         
361         // find the difference in angles
362         Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
363         Float_t dPhi   = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
364         // take into account boundary condition
365         if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi; 
366
367         // find the difference in z (between linear projection from layer 1
368         // and the actual point: Dzeta= z1/r1*r2 -z2)   
369         Float_t r2   = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
370         Float_t dZeta  = fClustersLay1[iC1][2]*r2 - fClustersLay2[iC2][2];
371
372         if (fHistOn) {
373           fhClustersDPhiAll->Fill(dPhi);    
374           fhClustersDThetaAll->Fill(dTheta);    
375           fhClustersDZetaAll->Fill(dZeta);    
376           fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
377           fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
378         }
379         // make "elliptical" cut in Phi and Zeta! 
380         Float_t d = TMath::Sqrt(TMath::Power(dPhi/fPhiWindow,2) + TMath::Power(dZeta/fZetaWindow,2));
381
382         if (d>1) continue;      
383         
384         //look for the minimum distance: the minimum is in iC2WithBestDist
385         if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < distmin ) {
386           distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi));
387           dPhimin = dPhi;
388           dThetamin = dTheta;
389           dZetamin = dZeta; 
390           iC2WithBestDist = iC2;
391         }
392       } 
393     } // end of loop over clusters in layer 2 
394     
395     if (distmin<100) { // This means that a cluster in layer 2 was found that mathes with iC1
396
397       if (fHistOn) {
398         fhClustersDPhiAcc->Fill(dPhimin);
399         fhClustersDThetaAcc->Fill(dThetamin);    
400         fhClustersDZetaAcc->Fill(dZetamin);    
401         fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
402         fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
403       }
404       
405       if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE; // flag the association
406       
407       // store the tracklet
408       
409       // use the theta from the clusters in the first layer
410       fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
411       // use the phi from the clusters in the first layer
412       fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
413       // store the difference between phi1 and phi2
414       fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
415
416       // define dphi in the range [0,pi] with proper sign (track charge correlated)
417       if (fTracklets[fNTracklets][2] > TMath::Pi()) 
418           fTracklets[fNTracklets][2] = fTracklets[fNTracklets][2]-2.*TMath::Pi();
419       if (fTracklets[fNTracklets][2] < -TMath::Pi()) 
420           fTracklets[fNTracklets][2] = fTracklets[fNTracklets][2]+2.*TMath::Pi();
421
422       // find label
423       // if equal label in both clusters found this label is assigned
424       // if no equal label can be found the first labels of the L1 AND L2 cluster are assigned
425       Int_t label1 = 0;
426       Int_t label2 = 0;
427       while (label2 < 3)
428       {
429         if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
430           break;
431         label1++;
432         if (label1 == 3)
433         {
434           label1 = 0;
435           label2++;
436         }
437       }
438
439       if (label2 < 3)
440       {
441         AliDebug(AliLog::kDebug, Form("Found label %d == %d for tracklet candidate %d\n", (Int_t) fClustersLay1[iC1][3+label1], (Int_t) fClustersLay2[iC2WithBestDist][3+label2], fNTracklets));
442         fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
443         fTracklets[fNTracklets][4] = fClustersLay2[iC2WithBestDist][3+label2];
444       }
445       else
446       {
447         AliDebug(AliLog::kDebug, Form("Did not find label %d %d %d %d %d %d for tracklet candidate %d\n", (Int_t) fClustersLay1[iC1][3], (Int_t) fClustersLay1[iC1][4], (Int_t) fClustersLay1[iC1][5], (Int_t) fClustersLay2[iC2WithBestDist][3], (Int_t) fClustersLay2[iC2WithBestDist][4], (Int_t) fClustersLay2[iC2WithBestDist][5], fNTracklets));
448         fTracklets[fNTracklets][3] = fClustersLay1[iC1][3];
449         fTracklets[fNTracklets][4] = fClustersLay2[iC2WithBestDist][3];
450       }
451
452       if (fHistOn) {
453         Float_t eta=fTracklets[fNTracklets][0];
454         eta= TMath::Tan(eta/2.);
455         eta=-TMath::Log(eta);
456         fhetaTracklets->Fill(eta);    
457         fhphiTracklets->Fill(fTracklets[fNTracklets][1]);    
458       }
459       
460       AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
461       AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", iC1, 
462                       iC2WithBestDist));
463       fNTracklets++;
464
465       if (fRemoveClustersFromOverlaps) FlagClustersInOverlapRegions (iC1,iC2WithBestDist);
466
467     }
468
469     // Delete the following else if you do not want to save Clusters! 
470
471     else { // This means that the cluster has not been associated
472
473       // store the cluster
474        
475       fSClusters[fNSingleCluster][0] = fClustersLay1[iC1][0];
476       fSClusters[fNSingleCluster][1] = fClustersLay1[iC1][1];
477       AliDebug(1,Form(" Adding a single cluster %d (cluster %d  of layer 1)",
478                       fNSingleCluster, iC1));
479       fNSingleCluster++;
480     }
481
482   } // end of loop over clusters in layer 1
483   
484   AliDebug(1,Form("%d tracklets found", fNTracklets));
485 }
486
487 //____________________________________________________________________
488 void
489 AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) {
490   // This method
491   // - gets the clusters from the cluster tree 
492   // - convert them into global coordinates 
493   // - store them in the internal arrays
494   // - count the number of cluster-fired chips
495   
496   AliDebug(1,"Loading clusters and cluster-fired chips ...");
497   
498   fNClustersLay1 = 0;
499   fNClustersLay2 = 0;
500   fNFiredChips[0] = 0;
501   fNFiredChips[1] = 0;
502   
503   AliITSsegmentationSPD *seg = new AliITSsegmentationSPD();
504
505   TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
506   TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
507
508   itsClusterBranch->SetAddress(&itsClusters);
509
510   Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();  
511   Float_t cluGlo[3]={0.,0.,0.};
512  
513   // loop over the its subdetectors
514   for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
515     
516     if (!itsClusterTree->GetEvent(iIts)) 
517       continue;
518     
519     Int_t nClusters = itsClusters->GetEntriesFast();
520
521     // number of clusters in each chip of the current module
522     Int_t nClustersInChip[5] = {0,0,0,0,0};
523     Int_t layer = 0;
524     
525     // loop over clusters
526     while(nClusters--) {
527       AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
528       
529       layer = cluster->GetLayer();
530       if (layer>1) continue;            
531       
532       cluster->GetGlobalXYZ(cluGlo);
533       Float_t x = cluGlo[0];
534       Float_t y = cluGlo[1];
535       Float_t z = cluGlo[2];      
536
537       // find the chip for the current cluster
538       Float_t locz = cluster->GetDetLocalZ();
539       Int_t iChip = seg->GetChipFromLocal(0,locz);
540       nClustersInChip[iChip]++; 
541       
542       if (layer==0) {
543         fClustersLay1[fNClustersLay1][0] = x;
544         fClustersLay1[fNClustersLay1][1] = y;
545         fClustersLay1[fNClustersLay1][2] = z;
546  
547         fDetectorIndexClustersLay1[fNClustersLay1]=iIts;  
548
549         for (Int_t i=0; i<3; i++)
550                 fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
551         fNClustersLay1++;
552       }
553       if (layer==1) {
554         fClustersLay2[fNClustersLay2][0] = x;
555         fClustersLay2[fNClustersLay2][1] = y;
556         fClustersLay2[fNClustersLay2][2] = z;
557
558         fDetectorIndexClustersLay2[fNClustersLay2]=iIts;
559  
560         for (Int_t i=0; i<3; i++)
561                 fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
562         fNClustersLay2++;
563       }
564       
565     }// end of cluster loop
566
567     // get number of fired chips in the current module
568     if(layer<2)
569     for(Int_t ifChip=0; ifChip<5; ifChip++) {
570       if(nClustersInChip[ifChip] >= 1)  fNFiredChips[layer]++;
571     }
572
573   } // end of its "subdetector" loop  
574
575   if (itsClusters) {
576     itsClusters->Delete();
577     delete itsClusters;
578     delete seg;
579     itsClusters = 0;
580   }
581   AliDebug(1,Form("(clusters in layer 1 : %d,  layer 2: %d)",fNClustersLay1,fNClustersLay2));
582   AliDebug(1,Form("(cluster-fired chips in layer 1 : %d,  layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
583 }
584 //____________________________________________________________________
585 void
586 AliITSMultReconstructor::LoadClusterFiredChips(TTree* itsClusterTree) {
587   // This method
588   // - gets the clusters from the cluster tree 
589   // - counts the number of (cluster)fired chips
590   
591   AliDebug(1,"Loading cluster-fired chips ...");
592   
593   fNFiredChips[0] = 0;
594   fNFiredChips[1] = 0;
595   
596   AliITSsegmentationSPD *seg = new AliITSsegmentationSPD();
597
598   TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
599   TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
600
601   itsClusterBranch->SetAddress(&itsClusters);
602
603   Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();  
604  
605   // loop over the its subdetectors
606   for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
607     
608     if (!itsClusterTree->GetEvent(iIts)) 
609       continue;
610     
611     Int_t nClusters = itsClusters->GetEntriesFast();
612
613     // number of clusters in each chip of the current module
614     Int_t nClustersInChip[5] = {0,0,0,0,0};
615     Int_t layer = 0;
616     
617     // loop over clusters
618     while(nClusters--) {
619       AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
620       
621       layer = cluster->GetLayer();
622       if (layer>1) continue;            
623
624       // find the chip for the current cluster
625       Float_t locz = cluster->GetDetLocalZ();
626       Int_t iChip = seg->GetChipFromLocal(0,locz);
627       nClustersInChip[iChip]++; 
628       
629     }// end of cluster loop
630
631     // get number of fired chips in the current module
632     if(layer<2)
633     for(Int_t ifChip=0; ifChip<5; ifChip++) {
634       if(nClustersInChip[ifChip] >= 1)  fNFiredChips[layer]++;
635     }
636
637   } // end of its "subdetector" loop  
638   
639   if (itsClusters) {
640     itsClusters->Delete();
641     delete itsClusters;
642     delete seg;
643     itsClusters = 0;
644   }
645   AliDebug(1,Form("(cluster-fired chips in layer 1 : %d,  layer 2: %d)",fNFiredChips[0],fNFiredChips[1]));
646 }
647 //____________________________________________________________________
648 void
649 AliITSMultReconstructor::SaveHists() {
650   // This method save the histograms on the output file
651   // (only if fHistOn is TRUE). 
652   
653   if (!fHistOn)
654     return;
655
656   fhClustersDPhiAll->Write();
657   fhClustersDThetaAll->Write();
658   fhClustersDZetaAll->Write();
659   fhDPhiVsDThetaAll->Write();
660   fhDPhiVsDZetaAll->Write();
661
662   fhClustersDPhiAcc->Write();
663   fhClustersDThetaAcc->Write();
664   fhClustersDZetaAcc->Write();
665   fhDPhiVsDThetaAcc->Write();
666   fhDPhiVsDZetaAcc->Write();
667
668   fhetaTracklets->Write();
669   fhphiTracklets->Write();
670   fhetaClustersLay1->Write();
671   fhphiClustersLay1->Write();
672 }
673
674 //____________________________________________________________________
675 void 
676 AliITSMultReconstructor::FlagClustersInOverlapRegions (Int_t iC1, Int_t iC2WithBestDist) {
677
678   Float_t distClSameMod=0.;
679   Float_t distClSameModMin=0.;
680   Int_t   iClOverlap =0;
681   Float_t meanRadiusLay1 = 3.99335; // average radius inner layer
682   Float_t meanRadiusLay2 = 7.37935; // average radius outer layer;
683
684   Float_t zproj1=0.;
685   Float_t zproj2=0.;
686   Float_t deZproj=0.;
687
688   // Loop on inner layer clusters
689   for (Int_t iiC1=0; iiC1<fNClustersLay1; iiC1++) {
690     if (!fOverlapFlagClustersLay1[iiC1]) {
691       // only for adjacent modules
692       if ((TMath::Abs(fDetectorIndexClustersLay1[iC1]-fDetectorIndexClustersLay1[iiC1])==4)||
693          (TMath::Abs(fDetectorIndexClustersLay1[iC1]-fDetectorIndexClustersLay1[iiC1])==76)) {
694         Float_t dePhi=TMath::Abs(fClustersLay1[iiC1][1]-fClustersLay1[iC1][1]);
695         if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
696
697         zproj1=meanRadiusLay1/TMath::Tan(fClustersLay1[iC1][0]);
698         zproj2=meanRadiusLay1/TMath::Tan(fClustersLay1[iiC1][0]);
699
700         deZproj=TMath::Abs(zproj1-zproj2);
701
702         distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
703         if (distClSameMod<=1.) fOverlapFlagClustersLay1[iiC1]=kTRUE;
704
705 //        if (distClSameMod<=1.) {
706 //          if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
707 //            distClSameModMin=distClSameMod;
708 //            iClOverlap=iiC1;
709 //          } 
710 //        }
711
712
713       } // end adjacent modules
714     } 
715   } // end Loop on inner layer clusters
716
717 //  if (distClSameModMin!=0.) fOverlapFlagClustersLay1[iClOverlap]=kTRUE;
718
719   distClSameMod=0.;
720   distClSameModMin=0.;
721   iClOverlap =0;
722   // Loop on outer layer clusters
723   for (Int_t iiC2=0; iiC2<fNClustersLay2; iiC2++) {
724     if (!fOverlapFlagClustersLay2[iiC2]) {
725       // only for adjacent modules
726       if ((TMath::Abs(fDetectorIndexClustersLay2[iC2WithBestDist]-fDetectorIndexClustersLay2[iiC2])==4) ||
727          (TMath::Abs(fDetectorIndexClustersLay2[iC2WithBestDist]-fDetectorIndexClustersLay2[iiC2])==156)) {
728         Float_t dePhi=TMath::Abs(fClustersLay2[iiC2][1]-fClustersLay2[iC2WithBestDist][1]);
729         if (dePhi>TMath::Pi()) dePhi=2.*TMath::Pi()-dePhi;
730
731         zproj1=meanRadiusLay2/TMath::Tan(fClustersLay2[iC2WithBestDist][0]);
732         zproj2=meanRadiusLay2/TMath::Tan(fClustersLay2[iiC2][0]);
733
734         deZproj=TMath::Abs(zproj1-zproj2);
735         distClSameMod = TMath::Sqrt(TMath::Power(deZproj/fZetaOverlapCut,2)+TMath::Power(dePhi/fPhiOverlapCut,2));
736         if (distClSameMod<=1.) fOverlapFlagClustersLay2[iiC2]=kTRUE;
737
738 //        if (distClSameMod<=1.) {
739 //          if (distClSameModMin==0. || distClSameMod<distClSameModMin) {
740 //            distClSameModMin=distClSameMod;
741 //            iClOverlap=iiC2;
742 //          }
743 //        }
744
745       } // end adjacent modules
746     }
747   } // end Loop on outer layer clusters
748
749 //  if (distClSameModMin!=0.) fOverlapFlagClustersLay2[iClOverlap]=kTRUE;
750
751 }