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