55549e1ee7d3770f88d8a1333760e2b1ac9d9f41
[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 #include "AliEMCALRecoUtils.h"
59
60 #include "AliEMCALTracker.h"
61
62 using std::cerr;
63 using std::endl;
64 ClassImp(AliEMCALTracker)
65
66 //
67 //------------------------------------------------------------------------------
68 //
69 AliEMCALTracker::AliEMCALTracker() 
70 : AliTracker(),
71   fCutPt(0),
72   fCutNITS(0),
73   fCutNTPC(50),
74   fStep(20),
75   fTrackCorrMode(kTrackCorrMMB),
76   fClusterWindow(50),
77   fCutEta(0.025),
78   fCutPhi(0.05),
79   fTracks(0),
80   fClusters(0),
81   fGeom(0)
82 {
83   //
84   // Default constructor.
85   // Initializes all simple data members to default values,
86    // and all collections to NULL.
87   // Output file name is set to a default value.
88   //
89   InitParameters();
90 }
91 //
92 //------------------------------------------------------------------------------
93 //
94 AliEMCALTracker::AliEMCALTracker(const AliEMCALTracker& copy) 
95   : AliTracker(),
96     fCutPt(copy.fCutPt),
97     fCutNITS(copy.fCutNITS),
98     fCutNTPC(copy.fCutNTPC),
99     fStep(copy.fStep),
100     fTrackCorrMode(copy.fTrackCorrMode),
101     fClusterWindow(copy.fClusterWindow),
102     fCutEta(copy.fCutEta),
103     fCutPhi(copy.fCutPhi),
104     fTracks((TObjArray*)copy.fTracks->Clone()),
105     fClusters((TObjArray*)copy.fClusters->Clone()),
106     fGeom(copy.fGeom)
107 {
108   //
109   // Copy constructor
110   // Besides copying all parameters, duplicates all collections.
111   //
112 }
113 //
114 //------------------------------------------------------------------------------
115 //
116 AliEMCALTracker& AliEMCALTracker::operator=(const AliEMCALTracker& source)
117 { // assignment operator; use copy ctor
118   if (&source == this) return *this;
119
120   new (this) AliEMCALTracker(source);
121   return *this;
122 }
123 //
124 //------------------------------------------------------------------------------
125 //
126 void AliEMCALTracker::InitParameters()
127 {
128   //
129   // Retrieve initialization parameters
130   //
131         
132   // Check if the instance of AliEMCALRecParam exists, 
133   const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
134
135   if(!recParam){
136     AliFatal("Reconstruction parameters for EMCAL not set!");
137   }
138   else{
139  
140     fCutEta  =  recParam->GetMthCutEta();
141     fCutPhi  =  recParam->GetMthCutPhi();
142     fStep    =  recParam->GetExtrapolateStep();
143     fCutPt   =  recParam->GetTrkCutPt();
144     fCutNITS =  recParam->GetTrkCutNITS();
145     fCutNTPC =  recParam->GetTrkCutNTPC();
146   }
147         
148 }
149
150 //
151 //------------------------------------------------------------------------------
152 //
153 void AliEMCALTracker::Clear(Option_t* option)
154 {
155         //
156         // Clearing method
157         // Deletes all objects in arrays and the arrays themselves
158         //
159
160         TString opt(option);
161         Bool_t clearTracks = opt.Contains("TRACKS");
162         Bool_t clearClusters = opt.Contains("CLUSTERS");
163         if (opt.Contains("ALL")) {
164                 clearTracks = kTRUE;
165                 clearClusters = kTRUE;
166         }
167         
168         //fTracks is a collection of esdTrack
169         //When clearing this array, the linked objects should not be deleted
170         if (fTracks != 0x0 && clearTracks) {
171            fTracks->Clear();
172            delete fTracks;
173            fTracks = 0;
174         }
175         if (fClusters != 0x0 && clearClusters) {
176            fClusters->Delete();
177            delete fClusters;
178            fClusters = 0;
179         }
180 }
181 //
182 //------------------------------------------------------------------------------
183 //
184 Int_t AliEMCALTracker::LoadClusters(TTree *cTree) 
185 {
186         //
187         // Load EMCAL clusters in the form of AliEMCALRecPoint,
188         // from simulation temporary files.
189         // (When included in reconstruction chain, this method is used automatically)
190         //
191         
192         Clear("CLUSTERS");
193
194         cTree->SetBranchStatus("*",0); //disable all branches
195         cTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
196
197         TBranch *branch = cTree->GetBranch("EMCALECARP");
198         if (!branch) {
199                 AliError("Can't get the branch with the EMCAL clusters");
200                 return 1;
201         }
202         
203         TClonesArray *clusters = new TClonesArray("AliEMCALRecPoint", 1000);
204         branch->SetAddress(&clusters);
205         
206         //cTree->GetEvent(0);
207         branch->GetEntry(0);
208         Int_t nClusters = (Int_t)clusters->GetEntries();
209         if(fClusters) fClusters->Delete();
210         else fClusters = new TObjArray(0);
211         for (Int_t i = 0; i < nClusters; i++) {
212                 AliEMCALRecPoint *cluster = (AliEMCALRecPoint*)clusters->At(i);
213                 if (!cluster) continue;
214                 AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
215                 fClusters->AddLast(matchCluster);
216         }
217
218         branch->SetAddress(0);
219         clusters->Delete();
220         delete clusters;
221
222         AliInfo(Form("Collected %d RecPoints from Tree", fClusters->GetEntries()));
223
224         return 0;
225 }
226 //
227 //------------------------------------------------------------------------------
228 //
229 Int_t AliEMCALTracker::LoadClusters(AliESDEvent *esd) 
230 {
231   //
232   // Load EMCAL clusters in the form of AliESDCaloClusters,
233   // from an AliESD object.
234   //
235   
236   // make sure that tracks/clusters collections are empty
237   Clear("CLUSTERS");
238   fClusters = new TObjArray(0);
239   
240   Int_t nClusters = esd->GetNumberOfCaloClusters();                     
241   for (Int_t i=0; i<nClusters; i++) 
242     {
243       AliESDCaloCluster *cluster = esd->GetCaloCluster(i);
244       if (!cluster || !cluster->IsEMCAL()) continue ; 
245       AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
246       fClusters->AddLast(matchCluster);
247     }
248   
249   AliInfo(Form("Collected %d clusters from ESD", fClusters->GetEntries()));
250   return 0;
251 }
252 //
253 //------------------------------------------------------------------------------
254 //
255 Int_t AliEMCALTracker::LoadTracks(AliESDEvent *esd)
256 {
257   //
258   // Load ESD tracks.
259   //
260         
261   Clear("TRACKS");
262   fTracks = new TObjArray(0);
263         
264   Int_t nTracks = esd->GetNumberOfTracks();
265   //Bool_t isKink=kFALSE;
266   for (Int_t i = 0; i < nTracks; i++) 
267     {
268       AliESDtrack *esdTrack = esd->GetTrack(i);
269       // set by default the value corresponding to "no match"
270       esdTrack->SetEMCALcluster(kUnmatched);
271       esdTrack->ResetStatus(AliESDtrack::kEMCALmatch);
272
273       //Select good quaulity tracks
274       if(esdTrack->Pt()<fCutPt) continue;
275       if(esdTrack->GetNcls(1)<fCutNTPC)continue;
276
277       //Loose geometric cut
278       Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
279       if(TMath::Abs(esdTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
280
281       fTracks->AddLast(esdTrack);
282     }
283
284       AliInfo(Form("Collected %d tracks", fTracks->GetEntries()));
285       return 0;
286 }
287 //
288 //------------------------------------------------------------------------------
289 //
290 void AliEMCALTracker::SetTrackCorrectionMode(Option_t *option)
291 {
292   //
293   // Set track correction mode
294   // gest the choice in string format and converts into 
295   // internal enum
296   //
297   
298   TString opt(option);
299   opt.ToUpper();
300   
301   if (!opt.CompareTo("NONE")) 
302     {
303       fTrackCorrMode = kTrackCorrNone;
304     }
305   else if (!opt.CompareTo("MMB")) 
306     {
307       fTrackCorrMode = kTrackCorrMMB;
308     }
309   else 
310     {
311       cerr << "E-AliEMCALTracker::SetTrackCorrectionMode '" << option << "': Unrecognized option" << endl;
312     }
313 }
314 //
315 //------------------------------------------------------------------------------
316 //
317 Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
318 {
319         //
320         // Main operation method.
321         // Gets external AliESD containing tracks to be matched.
322         // After executing match finding, stores in the same ESD object all infos
323         // and releases the object for further reconstruction steps.
324         //
325         //
326         // Note: should always return 0=OK, because otherwise all tracking
327         // is aborted for this event
328   
329         if (!esd) {
330                 AliError("NULL ESD passed");
331                 return 1;
332         }
333         
334         // step 1: collect clusters
335         Int_t okLoadClusters, nClusters;
336         if (!fClusters || (fClusters && fClusters->IsEmpty())) {
337                 okLoadClusters = LoadClusters(esd);
338         }
339         nClusters = fClusters->GetEntries();
340                 
341         // step 2: collect ESD tracks
342         Int_t nTracks, okLoadTracks;
343         okLoadTracks = LoadTracks(esd);
344         nTracks = fTracks->GetEntries();
345         
346         // step 3: for each track, find the closest cluster as matched within residual cuts
347         Int_t index=-1;
348         for (Int_t it = 0; it < nTracks; it++) 
349           {
350             AliESDtrack *track = (AliESDtrack*)fTracks->At(it);
351             index = FindMatchedCluster(track);
352             if (index>-1) 
353               {
354                 AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(index);
355                 track->SetEMCALcluster(cluster->Index());
356                 track->SetStatus(AliESDtrack::kEMCALmatch);
357               }
358           }
359
360         return 0;
361 }
362
363 //
364 //------------------------------------------------------------------------------
365 //
366 Int_t AliEMCALTracker::FindMatchedCluster(AliESDtrack *track)
367 {         
368   //
369   // For each track, extrapolate it to all the clusters
370   // Find the closest one as matched if the residuals (dEta, dPhi) satisfy the cuts
371   //
372
373   Float_t maxEta=fCutEta;
374   Float_t maxPhi=fCutPhi;
375   Int_t index = -1;
376   
377   // If the esdFriend is available, use the TPCOuter point as the starting point of extrapolation
378   // Otherwise use the TPCInner point
379   AliExternalTrackParam *trkParam = 0;
380   const AliESDfriendTrack*  friendTrack = track->GetFriendTrack();
381   if(friendTrack && friendTrack->GetTPCOut())
382     trkParam = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
383   else
384     trkParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
385   if(!trkParam) return index;
386
387
388   AliExternalTrackParam trkParamTmp(*trkParam);
389   Float_t eta, phi;
390   if(!AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&trkParamTmp, 430., track->GetMass(kTRUE), fStep, eta, phi))  return index;
391   track->SetTrackPhiEtaOnEMCal(phi,eta);
392   if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()) return index;
393
394   //Perform extrapolation
395   Double_t trkPos[3];
396   trkParamTmp.GetXYZ(trkPos);
397   Int_t nclusters = fClusters->GetEntries();
398   for(Int_t ic=0; ic<nclusters; ic++)
399     {
400       AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
401       Float_t clsPos[3] = {cluster->X(),cluster->Y(),cluster->Z()};
402       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));
403       if(dR > fClusterWindow) continue;
404       
405       AliExternalTrackParam trkParTmp(trkParamTmp);
406
407       Float_t tmpEta, tmpPhi;
408       if(!AliEMCALRecoUtils::ExtrapolateTrackToPosition(&trkParTmp, clsPos,track->GetMass(kTRUE), 5, tmpEta, tmpPhi)) continue;
409       if(TMath::Abs(tmpPhi)<TMath::Abs(maxPhi) && TMath::Abs(tmpEta)<TMath::Abs(maxEta))
410         {
411           maxPhi=tmpPhi;
412           maxEta=tmpEta;
413           index=ic;
414         }
415       }
416   return index;
417 }
418
419 //
420 //------------------------------------------------------------------------------
421 //
422 void AliEMCALTracker::UnloadClusters() 
423 {
424         //
425         // Free memory from all arrays
426         // This method is called after the local tracking step
427         // so we can safely delete everything 
428         //
429         
430         Clear();
431 }
432
433 //
434 //------------------------------------------------------------------------------
435 //
436 AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliEMCALRecPoint *recPoint)
437   : fIndex(index),
438     fX(0.),
439     fY(0.),
440     fZ(0.)
441 {
442         //
443         // Translates an AliEMCALRecPoint object into the internal format.
444         // Index of passed cluster in its native array must be specified.
445         //
446         TVector3 clpos;
447         recPoint->GetGlobalPosition(clpos);
448         
449         fX = clpos.X();
450         fY = clpos.Y();
451         fZ = clpos.Z();
452 }
453 //
454 //------------------------------------------------------------------------------
455 //
456 AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliESDCaloCluster *caloCluster)
457   : fIndex(index),
458     fX(0.),
459     fY(0.),
460     fZ(0.)
461 {
462         //
463         // Translates an AliESDCaloCluster object into the internal format.
464         // Index of passed cluster in its native array must be specified.
465         //
466         Float_t clpos[3]= {0., 0., 0.};
467         caloCluster->GetPosition(clpos);
468         
469         fX = (Double_t)clpos[0];
470         fY = (Double_t)clpos[1];
471         fZ = (Double_t)clpos[2];
472 }