add external setter for MC parents list
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALTracker.cxx
CommitLineData
fe17d4cb 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)
5970dfe2 30// Revised by Rongrong 2010-05-31 (rongrong.ma@cern.ch)
fe17d4cb 31//=========================================================================
32
33#include <Riostream.h>
34#include <iomanip>
35
36#include <TFile.h>
37#include <TTree.h>
fe17d4cb 38#include <TList.h>
39#include <TString.h>
40#include <TVector3.h>
41#include <TClonesArray.h>
c61f0e70 42#include <TGeoMatrix.h>
fe17d4cb 43
44#include "AliLog.h"
af885e0f 45#include "AliESDEvent.h"
fe17d4cb 46#include "AliESDtrack.h"
89ffc0b0 47#include "AliESDCaloCluster.h"
fe17d4cb 48#include "AliEMCALRecPoint.h"
49#include "AliRunLoader.h"
50#include "AliEMCALTrack.h"
51#include "AliEMCALLoader.h"
c61f0e70 52#include "AliEMCALGeometry.h"
8ba062b1 53#include "AliEMCALReconstructor.h"
54#include "AliEMCALRecParam.h"
55#include "AliCDBEntry.h"
56#include "AliCDBManager.h"
3e3faf55 57#include "AliEMCALReconstructor.h"
ee602376 58#include "AliEMCALRecoUtils.h"
c61f0e70 59
fe17d4cb 60#include "AliEMCALTracker.h"
61
0267cfa6 62using std::cerr;
63using std::endl;
fe17d4cb 64ClassImp(AliEMCALTracker)
8ba062b1 65
fe17d4cb 66//
67//------------------------------------------------------------------------------
68//
69AliEMCALTracker::AliEMCALTracker()
8fc351e3 70: AliTracker(),
71 fCutPt(0),
72 fCutNITS(0),
73 fCutNTPC(50),
da34fafe 74 fStep(20),
8fc351e3 75 fTrackCorrMode(kTrackCorrMMB),
76 fClusterWindow(50),
77 fCutEta(0.025),
78 fCutPhi(0.05),
42ceff04 79 fITSTrackSA(kFALSE),
8fc351e3 80 fTracks(0),
81 fClusters(0),
82 fGeom(0)
fe17d4cb 83{
5970dfe2 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 //
0832a2bf 90 InitParameters();
fe17d4cb 91}
92//
93//------------------------------------------------------------------------------
94//
95AliEMCALTracker::AliEMCALTracker(const AliEMCALTracker& copy)
96 : AliTracker(),
5970dfe2 97 fCutPt(copy.fCutPt),
dcd86c5d 98 fCutNITS(copy.fCutNITS),
99 fCutNTPC(copy.fCutNTPC),
5970dfe2 100 fStep(copy.fStep),
101 fTrackCorrMode(copy.fTrackCorrMode),
8fc351e3 102 fClusterWindow(copy.fClusterWindow),
5970dfe2 103 fCutEta(copy.fCutEta),
104 fCutPhi(copy.fCutPhi),
42ceff04 105 fITSTrackSA(copy.fITSTrackSA),
fe17d4cb 106 fTracks((TObjArray*)copy.fTracks->Clone()),
107 fClusters((TObjArray*)copy.fClusters->Clone()),
c61f0e70 108 fGeom(copy.fGeom)
fe17d4cb 109{
5970dfe2 110 //
111 // Copy constructor
112 // Besides copying all parameters, duplicates all collections.
113 //
fe17d4cb 114}
115//
116//------------------------------------------------------------------------------
117//
f1d9131f 118AliEMCALTracker& AliEMCALTracker::operator=(const AliEMCALTracker& source)
119{ // assignment operator; use copy ctor
120 if (&source == this) return *this;
5970dfe2 121
f1d9131f 122 new (this) AliEMCALTracker(source);
123 return *this;
fe17d4cb 124}
125//
126//------------------------------------------------------------------------------
127//
8ba062b1 128void AliEMCALTracker::InitParameters()
129{
5970dfe2 130 //
131 // Retrieve initialization parameters
132 //
8ba062b1 133
134 // Check if the instance of AliEMCALRecParam exists,
3e3faf55 135 const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
3a2a23e1 136
ba6de5ea 137 if(!recParam){
138 AliFatal("Reconstruction parameters for EMCAL not set!");
8ba062b1 139 }
41f05b8c 140 else{
5970dfe2 141
142 fCutEta = recParam->GetMthCutEta();
143 fCutPhi = recParam->GetMthCutPhi();
8fc351e3 144 fStep = recParam->GetExtrapolateStep();
5970dfe2 145 fCutPt = recParam->GetTrkCutPt();
8fc351e3 146 fCutNITS = recParam->GetTrkCutNITS();
147 fCutNTPC = recParam->GetTrkCutNTPC();
41f05b8c 148 }
3e3faf55 149
8ba062b1 150}
5970dfe2 151
c61f0e70 152//
153//------------------------------------------------------------------------------
154//
fe17d4cb 155void AliEMCALTracker::Clear(Option_t* option)
156{
157 //
158 // Clearing method
044225d9 159 // Deletes all objects in arrays and the arrays themselves
fe17d4cb 160 //
161
162 TString opt(option);
163 Bool_t clearTracks = opt.Contains("TRACKS");
164 Bool_t clearClusters = opt.Contains("CLUSTERS");
fe17d4cb 165 if (opt.Contains("ALL")) {
166 clearTracks = kTRUE;
167 clearClusters = kTRUE;
fe17d4cb 168 }
169
5970dfe2 170 //fTracks is a collection of esdTrack
171 //When clearing this array, the linked objects should not be deleted
fe17d4cb 172 if (fTracks != 0x0 && clearTracks) {
5970dfe2 173 fTracks->Clear();
044225d9 174 delete fTracks;
175 fTracks = 0;
fe17d4cb 176 }
177 if (fClusters != 0x0 && clearClusters) {
044225d9 178 fClusters->Delete();
179 delete fClusters;
180 fClusters = 0;
fe17d4cb 181 }
fe17d4cb 182}
183//
184//------------------------------------------------------------------------------
185//
186Int_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
bce21ea7 196 cTree->SetBranchStatus("*",0); //disable all branches
197 cTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
198
fe17d4cb 199 TBranch *branch = cTree->GetBranch("EMCALECARP");
200 if (!branch) {
c61f0e70 201 AliError("Can't get the branch with the EMCAL clusters");
fe17d4cb 202 return 1;
203 }
204
9596d957 205 TClonesArray *clusters = new TClonesArray("AliEMCALRecPoint", 1000);
fe17d4cb 206 branch->SetAddress(&clusters);
fe17d4cb 207
bce21ea7 208 //cTree->GetEvent(0);
209 branch->GetEntry(0);
9596d957 210 Int_t nClusters = (Int_t)clusters->GetEntries();
bce21ea7 211 if(fClusters) fClusters->Delete();
212 else fClusters = new TObjArray(0);
fe17d4cb 213 for (Int_t i = 0; i < nClusters; i++) {
214 AliEMCALRecPoint *cluster = (AliEMCALRecPoint*)clusters->At(i);
215 if (!cluster) continue;
fe17d4cb 216 AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
217 fClusters->AddLast(matchCluster);
218 }
3d9e8b15 219
2ad4424e 220 branch->SetAddress(0);
221 clusters->Delete();
222 delete clusters;
223
5970dfe2 224 AliInfo(Form("Collected %d RecPoints from Tree", fClusters->GetEntries()));
fe17d4cb 225
226 return 0;
227}
228//
229//------------------------------------------------------------------------------
230//
af885e0f 231Int_t AliEMCALTracker::LoadClusters(AliESDEvent *esd)
fe17d4cb 232{
5970dfe2 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;
fe17d4cb 253}
254//
255//------------------------------------------------------------------------------
256//
af885e0f 257Int_t AliEMCALTracker::LoadTracks(AliESDEvent *esd)
fe17d4cb 258{
5970dfe2 259 //
260 // Load ESD tracks.
261 //
42ceff04 262
26050574 263 UInt_t mask1 = esd->GetESDRun()->GetDetectorsInDAQ();
264 UInt_t mask2 = esd->GetESDRun()->GetDetectorsInReco();
265 Bool_t desc1 = (mask1 >> 3) & 0x1;
266 Bool_t desc2 = (mask2 >> 3) & 0x1;
267 if (desc1==0 || desc2==0) {
30e29b2a 268// AliError(Form("TPC not in DAQ/RECO: %u (%u)/%u (%u)",
269// mask1, esd->GetESDRun()->GetDetectorsInReco(),
270// mask2, esd->GetESDRun()->GetDetectorsInDAQ()));
26050574 271 fITSTrackSA = kTRUE;
42ceff04 272 }
273
5970dfe2 274 Clear("TRACKS");
275 fTracks = new TObjArray(0);
fe17d4cb 276
5970dfe2 277 Int_t nTracks = esd->GetNumberOfTracks();
8fc351e3 278 //Bool_t isKink=kFALSE;
5970dfe2 279 for (Int_t i = 0; i < nTracks; i++)
280 {
281 AliESDtrack *esdTrack = esd->GetTrack(i);
282 // set by default the value corresponding to "no match"
283 esdTrack->SetEMCALcluster(kUnmatched);
8fc351e3 284 esdTrack->ResetStatus(AliESDtrack::kEMCALmatch);
5970dfe2 285
286 //Select good quaulity tracks
287 if(esdTrack->Pt()<fCutPt) continue;
42ceff04 288 if(!fITSTrackSA)
289 if(esdTrack->GetNcls(1)<fCutNTPC)continue;
5970dfe2 290
5970dfe2 291 //Loose geometric cut
292 Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
d874f457 293 if(TMath::Abs(esdTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
5970dfe2 294
295 fTracks->AddLast(esdTrack);
296 }
297
298 AliInfo(Form("Collected %d tracks", fTracks->GetEntries()));
299 return 0;
300}
301//
302//------------------------------------------------------------------------------
303//
304void 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 }
fe17d4cb 327}
328//
329//------------------------------------------------------------------------------
330//
af885e0f 331Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
fe17d4cb 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 //
3d9e8b15 339 //
340 // Note: should always return 0=OK, because otherwise all tracking
341 // is aborted for this event
5970dfe2 342
fe17d4cb 343 if (!esd) {
344 AliError("NULL ESD passed");
345 return 1;
346 }
347
5970dfe2 348 // step 1: collect clusters
c61f0e70 349 Int_t okLoadClusters, nClusters;
fe17d4cb 350 if (!fClusters || (fClusters && fClusters->IsEmpty())) {
fe17d4cb 351 okLoadClusters = LoadClusters(esd);
fe17d4cb 352 }
c61f0e70 353 nClusters = fClusters->GetEntries();
5970dfe2 354
355 // step 2: collect ESD tracks
356 Int_t nTracks, okLoadTracks;
357 okLoadTracks = LoadTracks(esd);
c61f0e70 358 nTracks = fTracks->GetEntries();
fe17d4cb 359
5970dfe2 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());
952d023a 370 track->SetStatus(AliESDtrack::kEMCALmatch);
5970dfe2 371 }
372 }
fe17d4cb 373
374 return 0;
375}
fe17d4cb 376
377//
378//------------------------------------------------------------------------------
379//
5970dfe2 380Int_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 //
17773e2e 386
ee602376 387 Float_t maxEta=fCutEta;
388 Float_t maxPhi=fCutPhi;
5970dfe2 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
8108ab9e 393 AliExternalTrackParam *trkParam = 0;
42ceff04 394
395 if(!fITSTrackSA){
396 const AliESDfriendTrack* friendTrack = track->GetFriendTrack();
397 if(friendTrack && friendTrack->GetTPCOut())
398 trkParam = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
399 else if(track->GetInnerParam())
400 trkParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
401 }
5970dfe2 402 else
42ceff04 403 trkParam = new AliExternalTrackParam(*track);
404
5970dfe2 405 if(!trkParam) return index;
42ceff04 406
8fc351e3 407 AliExternalTrackParam trkParamTmp(*trkParam);
ee602376 408 Float_t eta, phi;
30e29b2a 409 if(!AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&trkParamTmp, 430., track->GetMass(kTRUE), fStep, eta, phi)) {
410 if(fITSTrackSA) delete trkParam;
411 return index;
412 }
ed15417e 413 track->SetTrackPhiEtaOnEMCal(phi,eta);
30e29b2a 414 if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()){
415 if(fITSTrackSA) delete trkParam;
416 return index;
417 }
8fc351e3 418
5970dfe2 419 //Perform extrapolation
420 Double_t trkPos[3];
8fc351e3 421 trkParamTmp.GetXYZ(trkPos);
5970dfe2 422 Int_t nclusters = fClusters->GetEntries();
423 for(Int_t ic=0; ic<nclusters; ic++)
424 {
5970dfe2 425 AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
8fc351e3 426 Float_t clsPos[3] = {cluster->X(),cluster->Y(),cluster->Z()};
427 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));
42ceff04 428// printf("\n dR=%f,wind=%f\n",dR,fClusterWindow); //MARCEL
8fc351e3 429 if(dR > fClusterWindow) continue;
430
431 AliExternalTrackParam trkParTmp(trkParamTmp);
432
ee602376 433 Float_t tmpEta, tmpPhi;
da25bc37 434 if(!AliEMCALRecoUtils::ExtrapolateTrackToPosition(&trkParTmp, clsPos,track->GetMass(kTRUE), 5, tmpEta, tmpPhi)) continue;
5970dfe2 435 if(TMath::Abs(tmpPhi)<TMath::Abs(maxPhi) && TMath::Abs(tmpEta)<TMath::Abs(maxEta))
436 {
437 maxPhi=tmpPhi;
438 maxEta=tmpEta;
439 index=ic;
440 }
5970dfe2 441 }
30e29b2a 442
443 if(fITSTrackSA) delete trkParam;
5970dfe2 444 return index;
fe17d4cb 445}
5970dfe2 446
fe17d4cb 447//
448//------------------------------------------------------------------------------
449//
450void AliEMCALTracker::UnloadClusters()
451{
452 //
044225d9 453 // Free memory from all arrays
454 // This method is called after the local tracking step
455 // so we can safely delete everything
fe17d4cb 456 //
457
044225d9 458 Clear();
fe17d4cb 459}
04475328 460
461//
462//------------------------------------------------------------------------------
463//
fe17d4cb 464AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliEMCALRecPoint *recPoint)
465 : fIndex(index),
fe17d4cb 466 fX(0.),
467 fY(0.),
468 fZ(0.)
469{
470 //
471 // Translates an AliEMCALRecPoint object into the internal format.
472 // Index of passed cluster in its native array must be specified.
473 //
474 TVector3 clpos;
475 recPoint->GetGlobalPosition(clpos);
476
477 fX = clpos.X();
478 fY = clpos.Y();
479 fZ = clpos.Z();
480}
481//
482//------------------------------------------------------------------------------
483//
484AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliESDCaloCluster *caloCluster)
485 : fIndex(index),
fe17d4cb 486 fX(0.),
487 fY(0.),
488 fZ(0.)
489{
490 //
491 // Translates an AliESDCaloCluster object into the internal format.
492 // Index of passed cluster in its native array must be specified.
493 //
53e430a3 494 Float_t clpos[3]= {0., 0., 0.};
7592dfc4 495 caloCluster->GetPosition(clpos);
fe17d4cb 496
497 fX = (Double_t)clpos[0];
498 fY = (Double_t)clpos[1];
499 fZ = (Double_t)clpos[2];
500}