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