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