Global coordinates from AliCluster::GetGlobalXYZ. fGeometry no longer used in AliITSM...
[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 //
50 //  fGeometry removed
51 //____________________________________________________________________
52
53 #include <TClonesArray.h>
54 #include <TH1F.h>
55 #include <TH2F.h>
56 #include <TTree.h>
57
58 #include "AliITSMultReconstructor.h"
59 #include "AliITSRecPoint.h"
60 #include "AliITSgeom.h"
61 #include "AliLog.h"
62
63 //____________________________________________________________________
64 ClassImp(AliITSMultReconstructor)
65
66
67 //____________________________________________________________________
68 AliITSMultReconstructor::AliITSMultReconstructor():
69 fClustersLay1(0),
70 fClustersLay2(0),
71 fTracklets(0),
72 fSClusters(0),
73 fAssociationFlag(0),
74 fNClustersLay1(0),
75 fNClustersLay2(0),
76 fNTracklets(0),
77 fNSingleCluster(0),
78 fPhiWindow(0),
79 fZetaWindow(0),
80 fOnlyOneTrackletPerC2(0),
81 fHistOn(0),
82 fhClustersDPhiAcc(0),
83 fhClustersDThetaAcc(0),
84 fhClustersDZetaAcc(0),
85 fhClustersDPhiAll(0),
86 fhClustersDThetaAll(0),
87 fhClustersDZetaAll(0),
88 fhDPhiVsDThetaAll(0),
89 fhDPhiVsDThetaAcc(0),
90 fhDPhiVsDZetaAll(0),
91 fhDPhiVsDZetaAcc(0),
92 fhetaTracklets(0),
93 fhphiTracklets(0),
94 fhetaClustersLay1(0),
95 fhphiClustersLay1(0){
96   // Method to reconstruct the charged particles multiplicity with the 
97   // SPD (tracklets).
98
99
100   SetHistOn();
101   SetPhiWindow();
102   SetZetaWindow();
103   SetOnlyOneTrackletPerC2();
104
105   fClustersLay1       = new Float_t*[300000];
106   fClustersLay2       = new Float_t*[300000];
107   fTracklets          = new Float_t*[300000];
108   fSClusters          = new Float_t*[300000];
109   fAssociationFlag    = new Bool_t[300000];
110
111   for(Int_t i=0; i<300000; i++) {
112     fClustersLay1[i]       = new Float_t[6];
113     fClustersLay2[i]       = new Float_t[6];
114     fTracklets[i]          = new Float_t[4];
115     fSClusters[i]           = new Float_t[2];
116     fAssociationFlag[i]    = kFALSE;
117   }
118
119   // definition of histograms
120   fhClustersDPhiAcc   = new TH1F("dphiacc",  "dphi",  100,0.,0.1);
121   fhClustersDPhiAcc->SetDirectory(0);
122   fhClustersDThetaAcc = new TH1F("dthetaacc","dtheta",100,-0.1,0.1);
123   fhClustersDThetaAcc->SetDirectory(0);
124   fhClustersDZetaAcc = new TH1F("dzetaacc","dzeta",100,-1.,1.);
125   fhClustersDZetaAcc->SetDirectory(0);
126
127   fhDPhiVsDZetaAcc = new TH2F("dphiVsDzetaacc","",100,-1.,1.,100,0.,0.1);
128   fhDPhiVsDZetaAcc->SetDirectory(0);
129   fhDPhiVsDThetaAcc = new TH2F("dphiVsDthetaAcc","",100,-0.1,0.1,100,0.,0.1);
130   fhDPhiVsDThetaAcc->SetDirectory(0);
131
132   fhClustersDPhiAll   = new TH1F("dphiall",  "dphi",  100,0.0,0.5);
133   fhClustersDPhiAll->SetDirectory(0);
134   fhClustersDThetaAll = new TH1F("dthetaall","dtheta",100,-0.5,0.5);
135   fhClustersDThetaAll->SetDirectory(0);
136   fhClustersDZetaAll = new TH1F("dzetaall","dzeta",100,-5.,5.);
137   fhClustersDZetaAll->SetDirectory(0);
138
139   fhDPhiVsDZetaAll = new TH2F("dphiVsDzetaall","",100,-5.,5.,100,0.,0.5);
140   fhDPhiVsDZetaAll->SetDirectory(0);
141   fhDPhiVsDThetaAll = new TH2F("dphiVsDthetaAll","",100,-0.5,0.5,100,0.,0.5);
142   fhDPhiVsDThetaAll->SetDirectory(0);
143
144   fhetaTracklets  = new TH1F("etaTracklets",  "eta",  100,-2.,2.);
145   fhetaTracklets->SetDirectory(0);
146   fhphiTracklets  = new TH1F("phiTracklets",  "phi",  100, 0., 2*TMath::Pi());
147   fhphiTracklets->SetDirectory(0);
148   fhetaClustersLay1  = new TH1F("etaClustersLay1",  "etaCl1",  100,-2.,2.);
149   fhetaClustersLay1->SetDirectory(0);
150   fhphiClustersLay1  = new TH1F("phiClustersLay1", "phiCl1", 100, 0., 2*TMath::Pi());
151   fhphiClustersLay1->SetDirectory(0);
152 }
153
154 //______________________________________________________________________
155 AliITSMultReconstructor::AliITSMultReconstructor(const AliITSMultReconstructor &mr) : TObject(mr),
156 fClustersLay1(mr.fClustersLay1),
157 fClustersLay2(mr.fClustersLay2),
158 fTracklets(mr.fTracklets),
159 fSClusters(mr.fSClusters),
160 fAssociationFlag(mr.fAssociationFlag),
161 fNClustersLay1(mr.fNClustersLay1),
162 fNClustersLay2(mr.fNClustersLay2),
163 fNTracklets(mr.fNTracklets),
164 fNSingleCluster(mr.fNSingleCluster),
165 fPhiWindow(mr.fPhiWindow),
166 fZetaWindow(mr.fZetaWindow),
167 fOnlyOneTrackletPerC2(mr.fOnlyOneTrackletPerC2),
168 fHistOn(mr.fHistOn),
169 fhClustersDPhiAcc(mr.fhClustersDPhiAcc),
170 fhClustersDThetaAcc(mr.fhClustersDThetaAcc),
171 fhClustersDZetaAcc(mr.fhClustersDZetaAcc),
172 fhClustersDPhiAll(mr.fhClustersDPhiAll),
173 fhClustersDThetaAll(mr.fhClustersDThetaAll),
174 fhClustersDZetaAll(mr.fhClustersDZetaAll),
175 fhDPhiVsDThetaAll(mr.fhDPhiVsDThetaAll),
176 fhDPhiVsDThetaAcc(mr.fhDPhiVsDThetaAcc),
177 fhDPhiVsDZetaAll(mr.fhDPhiVsDZetaAll),
178 fhDPhiVsDZetaAcc(mr.fhDPhiVsDZetaAcc),
179 fhetaTracklets(mr.fhetaTracklets),
180 fhphiTracklets(mr.fhphiTracklets),
181 fhetaClustersLay1(mr.fhetaClustersLay1),
182 fhphiClustersLay1(mr.fhphiClustersLay1) {
183   // Copy constructor
184
185 }
186
187 //______________________________________________________________________
188 AliITSMultReconstructor& AliITSMultReconstructor::operator=(const AliITSMultReconstructor& mr){
189   // Assignment operator
190   this->~AliITSMultReconstructor();
191   new(this) AliITSMultReconstructor(mr);
192   return *this;
193 }
194
195 //______________________________________________________________________
196 AliITSMultReconstructor::~AliITSMultReconstructor(){
197   // Destructor
198
199   // delete histograms
200   delete fhClustersDPhiAcc;
201   delete fhClustersDThetaAcc;
202   delete fhClustersDZetaAcc;
203   delete fhClustersDPhiAll;
204   delete fhClustersDThetaAll;
205   delete fhClustersDZetaAll;
206   delete fhDPhiVsDThetaAll;
207   delete fhDPhiVsDThetaAcc;
208   delete fhDPhiVsDZetaAll;
209   delete fhDPhiVsDZetaAcc;
210   delete fhetaTracklets;
211   delete fhphiTracklets;
212   delete fhetaClustersLay1;
213   delete fhphiClustersLay1;
214
215   // delete arrays
216   for(Int_t i=0; i<300000; i++) {
217     delete [] fClustersLay1[i];
218     delete [] fClustersLay2[i];
219     delete [] fTracklets[i];
220     delete [] fSClusters[i];
221   }
222   delete [] fClustersLay1;
223   delete [] fClustersLay2;
224   delete [] fTracklets;
225   delete [] fSClusters;
226
227   delete [] fAssociationFlag;
228 }
229
230 //____________________________________________________________________
231 void
232 AliITSMultReconstructor::Reconstruct(TTree* clusterTree, Float_t* vtx, Float_t* /* vtxRes*/) {
233   //
234   // - calls LoadClusterArray that finds the position of the clusters
235   //   (in global coord) 
236   // - convert the cluster coordinates to theta, phi (seen from the
237   //   interaction vertex). The third coordinate is used for ....
238   // - makes an array of tracklets 
239   //   
240   // After this method has been called, the clusters of the two layers
241   // and the tracklets can be retrieved by calling the Get'er methods.
242
243   // reset counters
244   fNClustersLay1 = 0;
245   fNClustersLay2 = 0;
246   fNTracklets = 0; 
247   fNSingleCluster = 0; 
248   // loading the clusters 
249   LoadClusterArrays(clusterTree);
250
251   // find the tracklets
252   AliDebug(1,"Looking for tracklets... ");  
253
254   //###########################################################
255   // Loop on layer 1 : finding theta, phi and z 
256   for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {    
257     Float_t x = fClustersLay1[iC1][0] - vtx[0];
258     Float_t y = fClustersLay1[iC1][1] - vtx[1];
259     Float_t z = fClustersLay1[iC1][2] - vtx[2];
260
261     Float_t r    = TMath::Sqrt(TMath::Power(x,2) +
262                                TMath::Power(y,2) +
263                                TMath::Power(z,2));
264     
265     fClustersLay1[iC1][0] = TMath::ACos(z/r);                   // Store Theta
266     fClustersLay1[iC1][1] = TMath::Pi() + TMath::ATan2(-y,-x);  // Store Phi
267     fClustersLay1[iC1][2] = z/r;                                // Store scaled z
268     if (fHistOn) {
269       Float_t eta=fClustersLay1[iC1][0];
270       eta= TMath::Tan(eta/2.);
271       eta=-TMath::Log(eta);
272       fhetaClustersLay1->Fill(eta);    
273       fhphiClustersLay1->Fill(fClustersLay1[iC1][1]);
274     }      
275   }
276   
277   // Loop on layer 2 : finding theta, phi and r   
278   for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {    
279     Float_t x = fClustersLay2[iC2][0] - vtx[0];
280     Float_t y = fClustersLay2[iC2][1] - vtx[1];
281     Float_t z = fClustersLay2[iC2][2] - vtx[2];
282    
283     Float_t r    = TMath::Sqrt(TMath::Power(x,2) +
284                                TMath::Power(y,2) +
285                                TMath::Power(z,2));
286     
287     fClustersLay2[iC2][0] = TMath::ACos(z/r);                   // Store Theta
288     fClustersLay2[iC2][1] = TMath::Pi() + TMath::ATan2(-y,-x);  // Store Phi
289     fClustersLay2[iC2][2] = z;                                  // Store z
290
291  // this only needs to be initialized for the fNClustersLay2 first associations
292     fAssociationFlag[iC2] = kFALSE;
293   }  
294   
295   //###########################################################
296   // Loop on layer 1 
297   for (Int_t iC1=0; iC1<fNClustersLay1; iC1++) {    
298
299     // reset of variables for multiple candidates
300     Int_t  iC2WithBestDist = 0;     // reset 
301     Float_t distmin        = 100.;  // just to put a huge number! 
302     Float_t dPhimin        = 0.;  // Used for histograms only! 
303     Float_t dThetamin      = 0.;  // Used for histograms only! 
304     Float_t dZetamin       = 0.;  // Used for histograms only! 
305     
306     // Loop on layer 2 
307     for (Int_t iC2=0; iC2<fNClustersLay2; iC2++) {      
308       
309       // The following excludes double associations
310       if (!fAssociationFlag[iC2]) {
311         
312         // find the difference in angles
313         Float_t dTheta = fClustersLay2[iC2][0] - fClustersLay1[iC1][0];
314         Float_t dPhi   = TMath::Abs(fClustersLay2[iC2][1] - fClustersLay1[iC1][1]);
315         // take into account boundary condition
316         if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi; 
317
318         // find the difference in z (between linear projection from layer 1
319         // and the actual point: Dzeta= z1/r1*r2 -z2)   
320         Float_t r2   = fClustersLay2[iC2][2]/TMath::Cos(fClustersLay2[iC2][0]);
321         Float_t dZeta  = fClustersLay1[iC1][2]*r2 - fClustersLay2[iC2][2];
322
323         if (fHistOn) {
324           fhClustersDPhiAll->Fill(dPhi);    
325           fhClustersDThetaAll->Fill(dTheta);    
326           fhClustersDZetaAll->Fill(dZeta);    
327           fhDPhiVsDThetaAll->Fill(dTheta, dPhi);
328           fhDPhiVsDZetaAll->Fill(dZeta, dPhi);
329         }
330         // make "elliptical" cut in Phi and Zeta! 
331         Float_t d = TMath::Sqrt(TMath::Power(dPhi/fPhiWindow,2) + TMath::Power(dZeta/fZetaWindow,2));
332
333         if (d>1) continue;      
334         
335         //look for the minimum distance: the minimum is in iC2WithBestDist
336         if (TMath::Sqrt(dZeta*dZeta+(r2*dPhi*r2*dPhi)) < distmin ) {
337           distmin=TMath::Sqrt(dZeta*dZeta + (r2*dPhi*r2*dPhi));
338           dPhimin = dPhi;
339           dThetamin = dTheta;
340           dZetamin = dZeta; 
341           iC2WithBestDist = iC2;
342         }
343       } 
344     } // end of loop over clusters in layer 2 
345     
346     if (distmin<100) { // This means that a cluster in layer 2 was found that mathes with iC1
347
348       if (fHistOn) {
349         fhClustersDPhiAcc->Fill(dPhimin);
350         fhClustersDThetaAcc->Fill(dThetamin);    
351         fhClustersDZetaAcc->Fill(dZetamin);    
352         fhDPhiVsDThetaAcc->Fill(dThetamin, dPhimin);
353         fhDPhiVsDZetaAcc->Fill(dZetamin, dPhimin);
354       }
355       
356       if (fOnlyOneTrackletPerC2) fAssociationFlag[iC2WithBestDist] = kTRUE; // flag the association
357       
358       // store the tracklet
359       
360       // use the theta from the clusters in the first layer
361       fTracklets[fNTracklets][0] = fClustersLay1[iC1][0];
362       // use the phi from the clusters in the first layer
363       fTracklets[fNTracklets][1] = fClustersLay1[iC1][1];
364       // Store the difference between phi1 and phi2
365       fTracklets[fNTracklets][2] = fClustersLay1[iC1][1] - fClustersLay2[iC2WithBestDist][1];
366
367       // find label
368       Int_t label1 = 0;
369       Int_t label2 = 0;
370       while (label2 < 3)
371       {
372         if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
373           break;
374         label1++;
375         if (label1 == 3)
376         {
377           label1 = 0;
378           label2++;
379         }
380       }
381
382       if (label2 < 3)
383       {
384         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));
385         fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
386       }
387       else
388       {
389         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));
390         fTracklets[fNTracklets][3] = -2;
391       }
392
393       if (fHistOn) {
394         Float_t eta=fTracklets[fNTracklets][0];
395         eta= TMath::Tan(eta/2.);
396         eta=-TMath::Log(eta);
397         fhetaTracklets->Fill(eta);    
398         fhphiTracklets->Fill(fTracklets[fNTracklets][1]);    
399       }
400       
401       AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
402       AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", iC1, 
403                       iC2WithBestDist));
404       fNTracklets++;
405     }
406
407     // Delete the following else if you do not want to save Clusters! 
408
409     else { // This means that the cluster has not been associated
410
411       // store the cluster
412        
413       fSClusters[fNSingleCluster][0] = fClustersLay1[iC1][0];
414       fSClusters[fNSingleCluster][1] = fClustersLay1[iC1][1];
415       AliDebug(1,Form(" Adding a single cluster %d (cluster %d  of layer 1)",
416                       fNSingleCluster, iC1));
417       fNSingleCluster++;
418     }
419
420   } // end of loop over clusters in layer 1
421   
422   AliDebug(1,Form("%d tracklets found", fNTracklets));
423 }
424
425 //____________________________________________________________________
426 void
427 AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) {
428   // This method
429   // - gets the clusters from the cluster tree 
430   // - convert them into global coordinates 
431   // - store them in the internal arrays
432   
433   AliDebug(1,"Loading clusters ...");
434   
435   fNClustersLay1 = 0;
436   fNClustersLay2 = 0;
437   
438   TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
439   TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
440
441   itsClusterBranch->SetAddress(&itsClusters);
442
443   Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();  
444   Float_t cluGlo[3]={0.,0.,0.};
445  
446   // loop over the its subdetectors
447   for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
448     
449     if (!itsClusterTree->GetEvent(iIts)) 
450       continue;
451     
452     Int_t nClusters = itsClusters->GetEntriesFast();
453     
454     // loop over clusters
455     while(nClusters--) {
456       AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
457       
458       if (cluster->GetLayer()>1) 
459         continue;            
460       
461       cluster->GetGlobalXYZ(cluGlo);
462       Float_t x = cluGlo[0];
463       Float_t y = cluGlo[1];
464       Float_t z = cluGlo[2];      
465       
466       if (cluster->GetLayer()==0) {
467         fClustersLay1[fNClustersLay1][0] = x;
468         fClustersLay1[fNClustersLay1][1] = y;
469         fClustersLay1[fNClustersLay1][2] = z;
470         for (Int_t i=0; i<3; i++)
471                 fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
472         fNClustersLay1++;
473       }
474       if (cluster->GetLayer()==1) {
475         fClustersLay2[fNClustersLay2][0] = x;
476         fClustersLay2[fNClustersLay2][1] = y;
477         fClustersLay2[fNClustersLay2][2] = z;
478         for (Int_t i=0; i<3; i++)
479                 fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
480         fNClustersLay2++;
481       }
482       
483     }// end of cluster loop
484   } // end of its "subdetector" loop  
485   
486   if (itsClusters) {
487     itsClusters->Delete();
488     delete itsClusters;
489     itsClusters = 0;
490   }
491   AliDebug(1,Form("(clusters in layer 1 : %d,  layer 2: %d)",fNClustersLay1,fNClustersLay2));
492 }
493 //____________________________________________________________________
494 void
495 AliITSMultReconstructor::SaveHists() {
496   // This method save the histograms on the output file
497   // (only if fHistOn is TRUE). 
498   
499   if (!fHistOn)
500     return;
501
502   fhClustersDPhiAll->Write();
503   fhClustersDThetaAll->Write();
504   fhClustersDZetaAll->Write();
505   fhDPhiVsDThetaAll->Write();
506   fhDPhiVsDZetaAll->Write();
507
508   fhClustersDPhiAcc->Write();
509   fhClustersDThetaAcc->Write();
510   fhClustersDZetaAcc->Write();
511   fhDPhiVsDThetaAcc->Write();
512   fhDPhiVsDZetaAcc->Write();
513
514   fhetaTracklets->Write();
515   fhphiTracklets->Write();
516   fhetaClustersLay1->Write();
517   fhphiClustersLay1->Write();
518 }
519
520