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