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