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