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