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