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