Rongrong: 1. All the extrapolation methods are implemented in AliEMCALTracker; 2...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALTracker.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 //                       Class AliEMCALTracker 
17 //                      -----------------------
18 // Implementation of the track matching method between barrel tracks and
19 // EMCAL clusters.
20 // Besides algorithm implementation, some cuts are required to be set
21 // in order to define, for each track, an acceptance window where clusters
22 // are searched to find best match (if any).
23 // The class accepts as input an ESD container, and works directly on it,
24 // simply setting, for each of its tracks, the fEMCALindex flag, for each
25 // track which is matched to a cluster.
26 // In order to use method, one must launch PropagateBack().
27 //
28 // ------------------------------------------------------------------------
29 // author: A. Pulvirenti (alberto.pulvirenti@ct.infn.it)
30 // Revised by Rongrong 2010-05-31 (rongrong.ma@cern.ch)
31 //=========================================================================
32
33 #include <Riostream.h>
34 #include <iomanip>
35
36 #include <TFile.h>
37 #include <TTree.h>
38 #include <TList.h>
39 #include <TString.h>
40 #include <TVector3.h>
41 #include <TClonesArray.h>
42 #include <TGeoMatrix.h>
43
44 #include "AliLog.h"
45 #include "AliESDEvent.h"
46 #include "AliESDtrack.h"
47 #include "AliESDCaloCluster.h"
48 #include "AliEMCALRecPoint.h"
49 #include "AliRunLoader.h"
50 #include "AliEMCALTrack.h"
51 #include "AliEMCALLoader.h"
52 #include "AliEMCALGeometry.h"
53 #include "AliEMCALReconstructor.h"
54 #include "AliEMCALRecParam.h"
55 #include "AliCDBEntry.h"
56 #include "AliCDBManager.h"
57 #include "AliEMCALReconstructor.h"
58
59 #include "AliEMCALTracker.h"
60
61 ClassImp(AliEMCALTracker)
62
63 //
64 //------------------------------------------------------------------------------
65 //
66 AliEMCALTracker::AliEMCALTracker() 
67 : AliTracker(),
68   fCutPt(0),
69   fCutNITS(0),
70   fCutNTPC(50),
71   fStep(10),
72   fTrackCorrMode(kTrackCorrMMB),
73   fClusterWindow(50),
74   fCutEta(0.025),
75   fCutPhi(0.05),
76   fTracks(0),
77   fClusters(0),
78   fGeom(0)
79 {
80   //
81   // Default constructor.
82   // Initializes all simple data members to default values,
83    // and all collections to NULL.
84   // Output file name is set to a default value.
85   //
86   InitParameters();
87 }
88 //
89 //------------------------------------------------------------------------------
90 //
91 AliEMCALTracker::AliEMCALTracker(const AliEMCALTracker& copy) 
92   : AliTracker(),
93     fCutPt(copy.fCutPt),
94     fCutNITS(copy.fCutNITS),
95     fCutNTPC(copy.fCutNTPC),
96     fStep(copy.fStep),
97     fTrackCorrMode(copy.fTrackCorrMode),
98     fClusterWindow(copy.fClusterWindow),
99     fCutEta(copy.fCutEta),
100     fCutPhi(copy.fCutPhi),
101     fTracks((TObjArray*)copy.fTracks->Clone()),
102     fClusters((TObjArray*)copy.fClusters->Clone()),
103     fGeom(copy.fGeom)
104 {
105   //
106   // Copy constructor
107   // Besides copying all parameters, duplicates all collections.
108   //
109 }
110 //
111 //------------------------------------------------------------------------------
112 //
113 AliEMCALTracker& AliEMCALTracker::operator=(const AliEMCALTracker& copy)
114 {
115   //
116   // Assignment operator.
117   // Besides copying all parameters, duplicates all collections.        
118   //
119
120   fCutPt  = copy.fCutPt;
121   fClusterWindow = copy.fClusterWindow;
122   fCutEta = copy.fCutEta;
123   fCutPhi = copy.fCutPhi;       
124   fStep = copy.fStep;
125   fTrackCorrMode = copy.fTrackCorrMode;
126
127   fCutNITS = copy.fCutNITS;
128   fCutNTPC = copy.fCutNTPC;
129   
130   fTracks = (TObjArray*)copy.fTracks->Clone();
131   fClusters = (TObjArray*)copy.fClusters->Clone();
132   fGeom = copy.fGeom;
133   
134   return (*this);
135 }
136 //
137 //------------------------------------------------------------------------------
138 //
139 void AliEMCALTracker::InitParameters()
140 {
141   //
142   // Retrieve initialization parameters
143   //
144         
145   // Check if the instance of AliEMCALRecParam exists, 
146   const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
147
148   if(!recParam){
149     AliFatal("Reconstruction parameters for EMCAL not set!");
150   }
151   else{
152  
153     fCutEta  =  recParam->GetMthCutEta();
154     fCutPhi  =  recParam->GetMthCutPhi();
155     fStep    =  recParam->GetExtrapolateStep();
156     fCutPt   =  recParam->GetTrkCutPt();
157     fCutNITS =  recParam->GetTrkCutNITS();
158     fCutNTPC =  recParam->GetTrkCutNTPC();
159   }
160         
161 }
162
163 //
164 //------------------------------------------------------------------------------
165 //
166 void AliEMCALTracker::Clear(Option_t* option)
167 {
168         //
169         // Clearing method
170         // Deletes all objects in arrays and the arrays themselves
171         //
172
173         TString opt(option);
174         Bool_t clearTracks = opt.Contains("TRACKS");
175         Bool_t clearClusters = opt.Contains("CLUSTERS");
176         if (opt.Contains("ALL")) {
177                 clearTracks = kTRUE;
178                 clearClusters = kTRUE;
179         }
180         
181         //fTracks is a collection of esdTrack
182         //When clearing this array, the linked objects should not be deleted
183         if (fTracks != 0x0 && clearTracks) {
184            fTracks->Clear();
185            delete fTracks;
186            fTracks = 0;
187         }
188         if (fClusters != 0x0 && clearClusters) {
189            fClusters->Delete();
190            delete fClusters;
191            fClusters = 0;
192         }
193 }
194 //
195 //------------------------------------------------------------------------------
196 //
197 Int_t AliEMCALTracker::LoadClusters(TTree *cTree) 
198 {
199         //
200         // Load EMCAL clusters in the form of AliEMCALRecPoint,
201         // from simulation temporary files.
202         // (When included in reconstruction chain, this method is used automatically)
203         //
204         
205         Clear("CLUSTERS");
206
207         cTree->SetBranchStatus("*",0); //disable all branches
208         cTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
209
210         TBranch *branch = cTree->GetBranch("EMCALECARP");
211         if (!branch) {
212                 AliError("Can't get the branch with the EMCAL clusters");
213                 return 1;
214         }
215         
216         TClonesArray *clusters = new TClonesArray("AliEMCALRecPoint", 1000);
217         branch->SetAddress(&clusters);
218         
219         //cTree->GetEvent(0);
220         branch->GetEntry(0);
221         Int_t nClusters = (Int_t)clusters->GetEntries();
222         if(fClusters) fClusters->Delete();
223         else fClusters = new TObjArray(0);
224         for (Int_t i = 0; i < nClusters; i++) {
225                 AliEMCALRecPoint *cluster = (AliEMCALRecPoint*)clusters->At(i);
226                 if (!cluster) continue;
227                 AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
228                 fClusters->AddLast(matchCluster);
229         }
230
231         branch->SetAddress(0);
232         clusters->Delete();
233         delete clusters;
234
235         AliInfo(Form("Collected %d RecPoints from Tree", fClusters->GetEntries()));
236
237         return 0;
238 }
239 //
240 //------------------------------------------------------------------------------
241 //
242 Int_t AliEMCALTracker::LoadClusters(AliESDEvent *esd) 
243 {
244   //
245   // Load EMCAL clusters in the form of AliESDCaloClusters,
246   // from an AliESD object.
247   //
248   
249   // make sure that tracks/clusters collections are empty
250   Clear("CLUSTERS");
251   fClusters = new TObjArray(0);
252   
253   Int_t nClusters = esd->GetNumberOfCaloClusters();                     
254   for (Int_t i=0; i<nClusters; i++) 
255     {
256       AliESDCaloCluster *cluster = esd->GetCaloCluster(i);
257       if (!cluster || !cluster->IsEMCAL()) continue ; 
258       AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
259       fClusters->AddLast(matchCluster);
260     }
261   
262   AliInfo(Form("Collected %d clusters from ESD", fClusters->GetEntries()));
263   return 0;
264 }
265 //
266 //------------------------------------------------------------------------------
267 //
268 Int_t AliEMCALTracker::LoadTracks(AliESDEvent *esd)
269 {
270   //
271   // Load ESD tracks.
272   //
273         
274   Clear("TRACKS");
275   fTracks = new TObjArray(0);
276         
277   Int_t nTracks = esd->GetNumberOfTracks();
278   //Bool_t isKink=kFALSE;
279   for (Int_t i = 0; i < nTracks; i++) 
280     {
281       AliESDtrack *esdTrack = esd->GetTrack(i);
282       // set by default the value corresponding to "no match"
283       esdTrack->SetEMCALcluster(kUnmatched);
284       esdTrack->ResetStatus(AliESDtrack::kEMCALmatch);
285
286       //Select good quaulity tracks
287       if(esdTrack->Pt()<fCutPt) continue;
288       if(esdTrack->GetNcls(1)<fCutNTPC)continue;
289
290       //Loose geometric cut
291       Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
292       if(TMath::Abs(esdTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
293
294       fTracks->AddLast(esdTrack);
295     }
296
297       AliInfo(Form("Collected %d tracks", fTracks->GetEntries()));
298       return 0;
299 }
300 //
301 //------------------------------------------------------------------------------
302 //
303 void AliEMCALTracker::SetTrackCorrectionMode(Option_t *option)
304 {
305   //
306   // Set track correction mode
307   // gest the choice in string format and converts into 
308   // internal enum
309   //
310   
311   TString opt(option);
312   opt.ToUpper();
313   
314   if (!opt.CompareTo("NONE")) 
315     {
316       fTrackCorrMode = kTrackCorrNone;
317     }
318   else if (!opt.CompareTo("MMB")) 
319     {
320       fTrackCorrMode = kTrackCorrMMB;
321     }
322   else 
323     {
324       cerr << "E-AliEMCALTracker::SetTrackCorrectionMode '" << option << "': Unrecognized option" << endl;
325     }
326 }
327 //
328 //------------------------------------------------------------------------------
329 //
330 Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
331 {
332         //
333         // Main operation method.
334         // Gets external AliESD containing tracks to be matched.
335         // After executing match finding, stores in the same ESD object all infos
336         // and releases the object for further reconstruction steps.
337         //
338         //
339         // Note: should always return 0=OK, because otherwise all tracking
340         // is aborted for this event
341   
342         if (!esd) {
343                 AliError("NULL ESD passed");
344                 return 1;
345         }
346         
347         // step 1: collect clusters
348         Int_t okLoadClusters, nClusters;
349         if (!fClusters || (fClusters && fClusters->IsEmpty())) {
350                 okLoadClusters = LoadClusters(esd);
351         }
352         nClusters = fClusters->GetEntries();
353                 
354         // step 2: collect ESD tracks
355         Int_t nTracks, okLoadTracks;
356         okLoadTracks = LoadTracks(esd);
357         nTracks = fTracks->GetEntries();
358         
359         // step 3: for each track, find the closest cluster as matched within residual cuts
360         Int_t index=-1;
361         for (Int_t it = 0; it < nTracks; it++) 
362           {
363             AliESDtrack *track = (AliESDtrack*)fTracks->At(it);
364             index = FindMatchedCluster(track);
365             if (index>-1) 
366               {
367                 AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(index);
368                 track->SetEMCALcluster(cluster->Index());
369                 track->SetStatus(AliESDtrack::kEMCALmatch);
370               }
371           }
372
373         return 0;
374 }
375
376 //
377 //------------------------------------------------------------------------------
378 //
379 Int_t AliEMCALTracker::FindMatchedCluster(AliESDtrack *track)
380 {         
381   //
382   // For each track, extrapolate it to all the clusters
383   // Find the closest one as matched if the residuals (dEta, dPhi) satisfy the cuts
384   //
385
386   Double_t maxEta=fCutEta;
387   Double_t maxPhi=fCutPhi;
388   Int_t index = -1;
389   
390   // If the esdFriend is available, use the TPCOuter point as the starting point of extrapolation
391   // Otherwise use the TPCInner point
392   AliExternalTrackParam *trkParam = 0;
393   const AliESDfriendTrack*  friendTrack = track->GetFriendTrack();
394   if(friendTrack && friendTrack->GetTPCOut())
395     trkParam = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
396   else
397     trkParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
398   if(!trkParam) return index;
399
400
401   AliExternalTrackParam trkParamTmp(*trkParam);
402   Double_t eta, phi;
403   if(!ExtrapolateTrackToEMCalSurface(&trkParamTmp, 430., track->GetMass(), fStep, eta, phi)) return index;
404   if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()) return index;
405
406   //Perform extrapolation
407   Double_t trkPos[3];
408   trkParamTmp.GetXYZ(trkPos);
409   Int_t nclusters = fClusters->GetEntries();
410   for(Int_t ic=0; ic<nclusters; ic++)
411     {
412       AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
413       Float_t clsPos[3] = {cluster->X(),cluster->Y(),cluster->Z()};
414       Double_t dR = TMath::Sqrt(TMath::Power(trkPos[0]-clsPos[0],2)+TMath::Power(trkPos[1]-clsPos[1],2)+TMath::Power(trkPos[2]-clsPos[2],2));
415       if(dR > fClusterWindow) continue;
416       
417       AliExternalTrackParam trkParTmp(trkParamTmp);
418
419       Double_t tmpEta, tmpPhi;
420       if(!ExtrapolateTrackToPosition(&trkParTmp, clsPos,track->GetMass(), fStep, tmpEta, tmpPhi)) continue;
421       if(TMath::Abs(tmpPhi)<TMath::Abs(maxPhi) && TMath::Abs(tmpEta)<TMath::Abs(maxEta))
422         {
423           maxPhi=tmpPhi;
424           maxEta=tmpEta;
425           index=ic;
426         }
427       }
428   return index;
429 }
430
431
432 //
433 //------------------------------------------------------------------------------
434 //
435 Bool_t AliEMCALTracker::ExtrapolateTrackToEMCalSurface(AliExternalTrackParam *trkParam, Double_t emcalR, Double_t mass, Double_t step, Double_t &eta, Double_t &phi)
436 {
437   eta = -999, phi = -999;
438   if(!trkParam) return kFALSE;
439   if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, emcalR, mass, step, kTRUE, 0.8, -1)) return kFALSE;
440   Double_t trkPos[3] = {0.,0.,0.};
441   if(!trkParam->GetXYZ(trkPos)) return kFALSE;
442   TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
443   eta = trkPosVec.Eta();
444   phi = trkPosVec.Phi();
445   if(phi<0)
446     phi += 2*TMath::Pi();
447
448   return kTRUE;
449 }
450
451
452 //
453 //------------------------------------------------------------------------------
454 //
455 Bool_t AliEMCALTracker::ExtrapolateTrackToPosition(AliExternalTrackParam *trkParam, Float_t *clsPos, Double_t mass, Double_t step, Double_t &tmpEta, Double_t &tmpPhi)
456 {
457   //
458   //Return the residual by extrapolating a track param to a global position
459   //
460   tmpEta = -999;
461   tmpPhi = -999;
462   if(!trkParam) return kFALSE;
463   Double_t trkPos[3] = {0.,0.,0.};
464   TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
465   Double_t alpha =  ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
466   vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
467   if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), mass, step,kTRUE, 0.8, -1)) return kFALSE;
468   if(!trkParam->GetXYZ(trkPos)) return kFALSE; //Get the extrapolated global position
469
470   TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
471   TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
472
473   // track cluster matching
474   tmpPhi = clsPosVec.DeltaPhi(trkPosVec);    // tmpPhi is between -pi and pi
475   tmpEta = clsPosVec.Eta()-trkPosVec.Eta();
476
477   return kTRUE;
478 }
479
480
481 //
482 //------------------------------------------------------------------------------
483 //
484 Bool_t AliEMCALTracker::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Double_t mass, Double_t step, Double_t &tmpEta, Double_t &tmpPhi)
485 {
486   //
487   //Return the residual by extrapolating a track param to a cluster
488   //
489   tmpEta = -999;
490   tmpPhi = -999;
491   if(!cluster) return kFALSE;
492
493   Float_t clsPos[3] = {0.,0.,0.};
494   cluster->GetPosition(clsPos);
495
496   return ExtrapolateTrackToPosition(trkParam, clsPos, mass, step, tmpEta, tmpPhi);
497 }
498
499
500 //
501 //------------------------------------------------------------------------------
502 //
503 void AliEMCALTracker::UnloadClusters() 
504 {
505         //
506         // Free memory from all arrays
507         // This method is called after the local tracking step
508         // so we can safely delete everything 
509         //
510         
511         Clear();
512 }
513
514 //
515 //------------------------------------------------------------------------------
516 //
517 AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliEMCALRecPoint *recPoint)
518   : fIndex(index),
519     fX(0.),
520     fY(0.),
521     fZ(0.)
522 {
523         //
524         // Translates an AliEMCALRecPoint object into the internal format.
525         // Index of passed cluster in its native array must be specified.
526         //
527         TVector3 clpos;
528         recPoint->GetGlobalPosition(clpos);
529         
530         fX = clpos.X();
531         fY = clpos.Y();
532         fZ = clpos.Z();
533 }
534 //
535 //------------------------------------------------------------------------------
536 //
537 AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliESDCaloCluster *caloCluster)
538   : fIndex(index),
539     fX(0.),
540     fY(0.),
541     fZ(0.)
542 {
543         //
544         // Translates an AliESDCaloCluster object into the internal format.
545         // Index of passed cluster in its native array must be specified.
546         //
547         Float_t clpos[3]= {0., 0., 0.};
548         caloCluster->GetPosition(clpos);
549         
550         fX = (Double_t)clpos[0];
551         fY = (Double_t)clpos[1];
552         fZ = (Double_t)clpos[2];
553 }