change default max distance R to 0.1
[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   fEMCalSurfaceDistance(440),
77   fClusterWindow(50),
78   fCutEta(0.025),
79   fCutPhi(0.05),
80   fITSTrackSA(kFALSE),
81   fTracks(0),
82   fClusters(0),
83   fGeom(0)
84 {
85   //
86   // Default constructor.
87   // Initializes all simple data members to default values,
88    // and all collections to NULL.
89   // Output file name is set to a default value.
90   //
91   InitParameters();
92 }
93 //
94 //------------------------------------------------------------------------------
95 //
96 AliEMCALTracker::AliEMCALTracker(const AliEMCALTracker& copy) : 
97   AliTracker(),
98   fCutPt(copy.fCutPt),
99   fCutNITS(copy.fCutNITS),
100   fCutNTPC(copy.fCutNTPC),
101   fStep(copy.fStep),
102   fTrackCorrMode(copy.fTrackCorrMode),
103   fEMCalSurfaceDistance(copy.fEMCalSurfaceDistance),
104   fClusterWindow(copy.fClusterWindow),
105   fCutEta(copy.fCutEta),
106   fCutPhi(copy.fCutPhi),
107   fITSTrackSA(copy.fITSTrackSA), 
108   fTracks((TObjArray*)copy.fTracks->Clone()),
109   fClusters((TObjArray*)copy.fClusters->Clone()),
110   fGeom(copy.fGeom)
111 {
112   //
113   // Copy constructor
114   // Besides copying all parameters, duplicates all collections.
115   //
116 }
117 //
118 //------------------------------------------------------------------------------
119 //
120 AliEMCALTracker& AliEMCALTracker::operator=(const AliEMCALTracker& source)
121 { // assignment operator; use copy ctor
122   if (&source == this) return *this;
123
124   new (this) AliEMCALTracker(source);
125   return *this;
126 }
127 //
128 //------------------------------------------------------------------------------
129 //
130 void AliEMCALTracker::InitParameters()
131 {
132   //
133   // Retrieve initialization parameters
134   //
135         
136   // Check if the instance of AliEMCALRecParam exists, 
137   const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
138
139   if (!recParam) {
140     AliFatal("Reconstruction parameters for EMCAL not set!");
141   } else {
142     fCutEta  =  recParam->GetMthCutEta();
143     fCutPhi  =  recParam->GetMthCutPhi();
144     fStep    =  recParam->GetExtrapolateStep();
145     fCutPt   =  recParam->GetTrkCutPt();
146     fCutNITS =  recParam->GetTrkCutNITS();
147     fCutNTPC =  recParam->GetTrkCutNTPC();
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     AliESDCaloCluster *cluster = esd->GetCaloCluster(i);
243     if (!cluster || !cluster->IsEMCAL()) continue ; 
244     AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
245     fClusters->AddLast(matchCluster);
246   }
247   
248   AliInfo(Form("Collected %d clusters from ESD", fClusters->GetEntries()));
249   return 0;
250 }
251 //
252 //------------------------------------------------------------------------------
253 //
254 Int_t AliEMCALTracker::LoadTracks(AliESDEvent *esd)
255 {
256   //
257   // Load ESD tracks.
258   //
259
260   UInt_t mask1 = esd->GetESDRun()->GetDetectorsInDAQ();
261   UInt_t mask2 = esd->GetESDRun()->GetDetectorsInReco();
262   Bool_t desc1 = (mask1 >> 3) & 0x1;
263   Bool_t desc2 = (mask2 >> 3) & 0x1;
264   if (desc1==0 || desc2==0) {
265 //     AliError(Form("TPC not in DAQ/RECO: %u (%u)/%u (%u)",
266 //                   mask1, esd->GetESDRun()->GetDetectorsInReco(),
267 //                   mask2, esd->GetESDRun()->GetDetectorsInDAQ()));
268     fITSTrackSA = kTRUE;
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     AliESDtrack *esdTrack = esd->GetTrack(i);
278     // set by default the value corresponding to "no match"
279     esdTrack->SetEMCALcluster(kUnmatched);
280     esdTrack->ResetStatus(AliESDtrack::kEMCALmatch);
281     
282     //Select good quaulity tracks
283     if (esdTrack->Pt()<fCutPt) continue;
284     if (!fITSTrackSA)
285       if (esdTrack->GetNcls(1)<fCutNTPC) continue;
286     
287     //Loose geometric cut
288     Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
289     if (TMath::Abs(esdTrack->Eta())>0.9 || phi <= 10 || phi >= 250) continue;
290     fTracks->AddLast(esdTrack);
291   }
292   AliInfo(Form("Collected %d tracks", fTracks->GetEntries()));
293   return 0;
294 }
295 //
296 //------------------------------------------------------------------------------
297 //
298 void AliEMCALTracker::SetTrackCorrectionMode(Option_t *option)
299 {
300   //
301   // Set track correction mode
302   // gest the choice in string format and converts into 
303   // internal enum
304   //
305   
306   TString opt(option);
307   opt.ToUpper();
308   
309   if (!opt.CompareTo("NONE")) {
310     fTrackCorrMode = kTrackCorrNone;
311   } else if (!opt.CompareTo("MMB")) {
312     fTrackCorrMode = kTrackCorrMMB;
313   } else  {
314     cerr << "E-AliEMCALTracker::SetTrackCorrectionMode '" << option << "': Unrecognized option" << endl;
315   }
316 }
317 //
318 //------------------------------------------------------------------------------
319 //
320 Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
321 {
322   //
323   // Main operation method.
324   // Gets external AliESD containing tracks to be matched.
325   // After executing match finding, stores in the same ESD object all infos
326   // and releases the object for further reconstruction steps.
327   //
328   //
329   // Note: should always return 0=OK, because otherwise all tracking
330   // is aborted for this event
331   
332   if (!esd) {
333     AliError("NULL ESD passed");
334     return 1;
335   }
336         
337   // step 1: collect clusters
338   Int_t okLoadClusters, nClusters;
339   if (!fClusters || (fClusters && fClusters->IsEmpty())) {
340     okLoadClusters = LoadClusters(esd);
341   }
342   nClusters = fClusters->GetEntries();
343                 
344   // step 2: collect ESD tracks
345   Int_t nTracks, okLoadTracks;
346   okLoadTracks = LoadTracks(esd);
347   nTracks = fTracks->GetEntries();
348   
349   // step 3: for each track, find the closest cluster as matched within residual cuts
350   Int_t index=-1;
351   for (Int_t it = 0; it < nTracks; it++) {
352     AliESDtrack *track = (AliESDtrack*)fTracks->At(it);
353     index = FindMatchedCluster(track);
354     if (index>-1) {
355       AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(index);
356       track->SetEMCALcluster(cluster->Index());
357       track->SetStatus(AliESDtrack::kEMCALmatch);
358     }
359   }
360
361   return 0;
362 }
363
364 //
365 //------------------------------------------------------------------------------
366 //
367 Int_t AliEMCALTracker::FindMatchedCluster(AliESDtrack *track)
368 {         
369   //
370   // For each track, extrapolate it to all the clusters
371   // Find the closest one as matched if the residuals (dEta, dPhi) satisfy the cuts
372   //
373
374   Float_t maxEta=fCutEta;
375   Float_t maxPhi=fCutPhi;
376   Int_t index = -1;
377   
378   // If the esdFriend is available, use the TPCOuter point as the starting point of extrapolation
379   // Otherwise use the TPCInner point
380   AliExternalTrackParam *trkParam = 0;
381   
382   if (!fITSTrackSA) { 
383     const AliESDfriendTrack*  friendTrack = track->GetFriendTrack();
384     if (friendTrack && friendTrack->GetTPCOut())
385       trkParam = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
386     else if (track->GetInnerParam())
387       trkParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
388   }
389   else
390     trkParam = new AliExternalTrackParam(*track);
391   
392   if (!trkParam) return index;
393   
394   AliExternalTrackParam trkParamTmp(*trkParam);
395   Float_t eta, phi, pt;
396   if (!AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&trkParamTmp, fEMCalSurfaceDistance, track->GetMass(kTRUE), fStep, eta, phi, pt))  {
397     if (fITSTrackSA) delete trkParam;
398     return index;
399   }
400   track->SetTrackPhiEtaPtOnEMCal(phi,eta,pt);
401   if (TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()) {
402     if (fITSTrackSA) delete trkParam;
403     return index;
404   }
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     AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
412     Float_t clsPos[3] = {cluster->X(),cluster->Y(),cluster->Z()};
413     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));
414     //printf("\n dR=%f,wind=%f\n",dR,fClusterWindow); //MARCEL
415     if (dR > fClusterWindow) continue;
416       
417     AliExternalTrackParam trkParTmp(trkParamTmp);
418
419     Float_t tmpEta, tmpPhi;
420     if (!AliEMCALRecoUtils::ExtrapolateTrackToPosition(&trkParTmp, clsPos,track->GetMass(kTRUE), 5, tmpEta, tmpPhi)) continue;
421     if (TMath::Abs(tmpPhi)<TMath::Abs(maxPhi) && TMath::Abs(tmpEta)<TMath::Abs(maxEta)) {
422       maxPhi=tmpPhi;
423       maxEta=tmpEta;
424       index=ic;
425     }
426   }
427
428   if (fITSTrackSA) delete trkParam;
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 }