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