Track labels added to tracklets (J.F. Grosse-Oetringhaus)
[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[6];
110     fClustersLay2[i]       = new Float_t[6];
111     fTracklets[i]          = new Float_t[4];
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       // find label
361       Int_t label1 = 0;
362       Int_t label2 = 0;
363       while (label2 < 3)
364       {
365         if ((Int_t) fClustersLay1[iC1][3+label1] != -2 && (Int_t) fClustersLay1[iC1][3+label1] == (Int_t) fClustersLay2[iC2WithBestDist][3+label2])
366           break;
367
368         label1++;
369         if (label1 == 3)
370         {
371           label1 = 0;
372           label2++;
373         }
374       }
375
376       if (label2 < 3)
377       {
378         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));
379         fTracklets[fNTracklets][3] = fClustersLay1[iC1][3+label1];
380       }
381       else
382       {
383         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));
384         fTracklets[fNTracklets][3] = -2;
385       }
386
387       if (fHistOn) {
388         Float_t eta=fTracklets[fNTracklets][0];
389         eta= TMath::Tan(eta/2.);
390         eta=-TMath::Log(eta);
391         fhetaTracklets->Fill(eta);    
392         fhphiTracklets->Fill(fTracklets[fNTracklets][1]);    
393       }
394       
395       AliDebug(1,Form(" Adding tracklet candidate %d ", fNTracklets));
396       AliDebug(1,Form(" Cl. %d of Layer 1 and %d of Layer 2", iC1, 
397                       iC2WithBestDist));
398       fNTracklets++;
399     }
400
401     // Delete the following else if you do not want to save Clusters! 
402
403     else { // This means that the cluster has not been associated
404
405       // store the cluster
406        
407       fSClusters[fNSingleCluster][0] = fClustersLay1[iC1][0];
408       fSClusters[fNSingleCluster][1] = fClustersLay1[iC1][1];
409       AliDebug(1,Form(" Adding a single cluster %d (cluster %d  of layer 1)",
410                       fNSingleCluster, iC1));
411       fNSingleCluster++;
412     }
413
414   } // end of loop over clusters in layer 1
415   
416   AliDebug(1,Form("%d tracklets found", fNTracklets));
417 }
418
419 //____________________________________________________________________
420 void
421 AliITSMultReconstructor::LoadClusterArrays(TTree* itsClusterTree) {
422   // This method
423   // - gets the clusters from the cluster tree 
424   // - convert them into global coordinates 
425   // - store them in the internal arrays
426   
427   AliDebug(1,"Loading clusters ...");
428   
429   fNClustersLay1 = 0;
430   fNClustersLay2 = 0;
431   
432   TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
433   TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
434
435   itsClusterBranch->SetAddress(&itsClusters);
436
437   Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();  
438  
439   // loop over the its subdetectors
440   for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
441     
442     if (!itsClusterTree->GetEvent(iIts)) 
443       continue;
444     
445     Int_t nClusters = itsClusters->GetEntriesFast();
446     
447     // stuff needed to get the global coordinates
448     Double_t rot[9];   fGeometry->GetRotMatrix(iIts,rot);
449     Int_t lay,lad,det; fGeometry->GetModuleId(iIts,lay,lad,det);
450     Float_t tx,ty,tz;  fGeometry->GetTrans(lay,lad,det,tx,ty,tz);
451     
452     // Below:
453     // "alpha" is the angle from the global X-axis to the
454     //         local GEANT X'-axis  ( rot[0]=cos(alpha) and rot[1]=sin(alpha) )
455     // "phi" is the angle from the global X-axis to the
456     //       local cluster X"-axis
457     
458     Double_t alpha   = TMath::ATan2(rot[1],rot[0])+TMath::Pi();
459     Double_t itsPhi = TMath::Pi()/2+alpha;
460     
461     if (lay==1) itsPhi+=TMath::Pi();
462     Double_t cp=TMath::Cos(itsPhi), sp=TMath::Sin(itsPhi);
463     Double_t r=tx*cp+ty*sp;
464     
465     // loop over clusters
466     while(nClusters--) {
467       AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
468       
469       if (cluster->GetLayer()>1) 
470         continue;            
471       
472       Float_t x = r*cp - cluster->GetY()*sp;
473       Float_t y = r*sp + cluster->GetY()*cp;
474       Float_t z = cluster->GetZ();      
475       
476       if (cluster->GetLayer()==0) {
477         fClustersLay1[fNClustersLay1][0] = x;
478         fClustersLay1[fNClustersLay1][1] = y;
479         fClustersLay1[fNClustersLay1][2] = z;
480         for (Int_t i=0; i<3; i++)
481                 fClustersLay1[fNClustersLay1][3+i] = cluster->GetLabel(i);
482         fNClustersLay1++;
483       }
484       if (cluster->GetLayer()==1) {
485         fClustersLay2[fNClustersLay2][0] = x;
486         fClustersLay2[fNClustersLay2][1] = y;
487         fClustersLay2[fNClustersLay2][2] = z;
488         for (Int_t i=0; i<3; i++)
489                 fClustersLay2[fNClustersLay2][3+i] = cluster->GetLabel(i);
490         fNClustersLay2++;
491       }
492       
493     }// end of cluster loop
494   } // end of its "subdetector" loop  
495   
496   AliDebug(1,Form("(clusters in layer 1 : %d,  layer 2: %d)",fNClustersLay1,fNClustersLay2));
497 }
498 //____________________________________________________________________
499 void
500 AliITSMultReconstructor::SaveHists() {
501   // This method save the histograms on the output file
502   // (only if fHistOn is TRUE). 
503   
504   if (!fHistOn)
505     return;
506
507   fhClustersDPhiAll->Write();
508   fhClustersDThetaAll->Write();
509   fhClustersDZetaAll->Write();
510   fhDPhiVsDThetaAll->Write();
511   fhDPhiVsDZetaAll->Write();
512
513   fhClustersDPhiAcc->Write();
514   fhClustersDThetaAcc->Write();
515   fhClustersDZetaAcc->Write();
516   fhDPhiVsDThetaAcc->Write();
517   fhDPhiVsDZetaAcc->Write();
518
519   fhetaTracklets->Write();
520   fhphiTracklets->Write();
521   fhetaClustersLay1->Write();
522   fhphiClustersLay1->Write();
523 }
524