add set and getter for neutral energy fraction
[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"
c61f0e70 58
fe17d4cb 59#include "AliEMCALTracker.h"
60
61ClassImp(AliEMCALTracker)
8ba062b1 62
fe17d4cb 63//
64//------------------------------------------------------------------------------
65//
66AliEMCALTracker::AliEMCALTracker()
67 : AliTracker(),
5970dfe2 68 fCutPt(0),
69 fCutNITS(0),
70 fCutNTPC(50),
71 fStep(50),
72 fTrackCorrMode(kTrackCorrMMB),
73 fCutEta(0.025),
74 fCutPhi(0.05),
fe17d4cb 75 fTracks(0),
76 fClusters(0),
c61f0e70 77 fGeom(0)
fe17d4cb 78{
5970dfe2 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 //
0832a2bf 85 InitParameters();
fe17d4cb 86}
87//
88//------------------------------------------------------------------------------
89//
90AliEMCALTracker::AliEMCALTracker(const AliEMCALTracker& copy)
91 : AliTracker(),
5970dfe2 92 fCutPt(copy.fCutPt),
dcd86c5d 93 fCutNITS(copy.fCutNITS),
94 fCutNTPC(copy.fCutNTPC),
5970dfe2 95 fStep(copy.fStep),
96 fTrackCorrMode(copy.fTrackCorrMode),
97 fCutEta(copy.fCutEta),
98 fCutPhi(copy.fCutPhi),
fe17d4cb 99 fTracks((TObjArray*)copy.fTracks->Clone()),
100 fClusters((TObjArray*)copy.fClusters->Clone()),
c61f0e70 101 fGeom(copy.fGeom)
fe17d4cb 102{
5970dfe2 103 //
104 // Copy constructor
105 // Besides copying all parameters, duplicates all collections.
106 //
fe17d4cb 107}
108//
109//------------------------------------------------------------------------------
110//
111AliEMCALTracker& AliEMCALTracker::operator=(const AliEMCALTracker& copy)
112{
5970dfe2 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);
fe17d4cb 132}
133//
134//------------------------------------------------------------------------------
135//
8ba062b1 136void AliEMCALTracker::InitParameters()
137{
5970dfe2 138 //
139 // Retrieve initialization parameters
140 //
8ba062b1 141
142 // Check if the instance of AliEMCALRecParam exists,
3e3faf55 143 const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
3a2a23e1 144
ba6de5ea 145 if(!recParam){
146 AliFatal("Reconstruction parameters for EMCAL not set!");
8ba062b1 147 }
41f05b8c 148 else{
5970dfe2 149
150 fCutEta = recParam->GetMthCutEta();
151 fCutPhi = recParam->GetMthCutPhi();
152 fStep = recParam->GetExtrapolateStep();
153 fCutPt = recParam->GetTrkCutPt();
41f05b8c 154 fCutNITS = recParam->GetTrkCutNITS();
155 fCutNTPC = recParam->GetTrkCutNTPC();
156 }
3e3faf55 157
8ba062b1 158}
5970dfe2 159
c61f0e70 160//
161//------------------------------------------------------------------------------
162//
fe17d4cb 163void AliEMCALTracker::Clear(Option_t* option)
164{
165 //
166 // Clearing method
044225d9 167 // Deletes all objects in arrays and the arrays themselves
fe17d4cb 168 //
169
170 TString opt(option);
171 Bool_t clearTracks = opt.Contains("TRACKS");
172 Bool_t clearClusters = opt.Contains("CLUSTERS");
fe17d4cb 173 if (opt.Contains("ALL")) {
174 clearTracks = kTRUE;
175 clearClusters = kTRUE;
fe17d4cb 176 }
177
5970dfe2 178 //fTracks is a collection of esdTrack
179 //When clearing this array, the linked objects should not be deleted
fe17d4cb 180 if (fTracks != 0x0 && clearTracks) {
5970dfe2 181 fTracks->Clear();
044225d9 182 delete fTracks;
183 fTracks = 0;
fe17d4cb 184 }
185 if (fClusters != 0x0 && clearClusters) {
044225d9 186 fClusters->Delete();
187 delete fClusters;
188 fClusters = 0;
fe17d4cb 189 }
fe17d4cb 190}
191//
192//------------------------------------------------------------------------------
193//
194Int_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
bce21ea7 204 cTree->SetBranchStatus("*",0); //disable all branches
205 cTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
206
fe17d4cb 207 TBranch *branch = cTree->GetBranch("EMCALECARP");
208 if (!branch) {
c61f0e70 209 AliError("Can't get the branch with the EMCAL clusters");
fe17d4cb 210 return 1;
211 }
212
9596d957 213 TClonesArray *clusters = new TClonesArray("AliEMCALRecPoint", 1000);
fe17d4cb 214 branch->SetAddress(&clusters);
fe17d4cb 215
bce21ea7 216 //cTree->GetEvent(0);
217 branch->GetEntry(0);
9596d957 218 Int_t nClusters = (Int_t)clusters->GetEntries();
bce21ea7 219 if(fClusters) fClusters->Delete();
220 else fClusters = new TObjArray(0);
fe17d4cb 221 for (Int_t i = 0; i < nClusters; i++) {
222 AliEMCALRecPoint *cluster = (AliEMCALRecPoint*)clusters->At(i);
223 if (!cluster) continue;
fe17d4cb 224 AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
225 fClusters->AddLast(matchCluster);
226 }
3d9e8b15 227
2ad4424e 228 branch->SetAddress(0);
229 clusters->Delete();
230 delete clusters;
231
5970dfe2 232 AliInfo(Form("Collected %d RecPoints from Tree", fClusters->GetEntries()));
fe17d4cb 233
234 return 0;
235}
236//
237//------------------------------------------------------------------------------
238//
af885e0f 239Int_t AliEMCALTracker::LoadClusters(AliESDEvent *esd)
fe17d4cb 240{
5970dfe2 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;
fe17d4cb 261}
262//
263//------------------------------------------------------------------------------
264//
af885e0f 265Int_t AliEMCALTracker::LoadTracks(AliESDEvent *esd)
fe17d4cb 266{
5970dfe2 267 //
268 // Load ESD tracks.
269 //
fe17d4cb 270
5970dfe2 271 Clear("TRACKS");
272 fTracks = new TObjArray(0);
fe17d4cb 273
5970dfe2 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;
fe17d4cb 291 }
5970dfe2 292 if (isKink) continue;
fe17d4cb 293
5970dfe2 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//
307void 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 }
fe17d4cb 330}
331//
332//------------------------------------------------------------------------------
333//
af885e0f 334Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
fe17d4cb 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 //
3d9e8b15 342 //
343 // Note: should always return 0=OK, because otherwise all tracking
344 // is aborted for this event
5970dfe2 345
fe17d4cb 346 if (!esd) {
347 AliError("NULL ESD passed");
348 return 1;
349 }
350
5970dfe2 351 // step 1: collect clusters
c61f0e70 352 Int_t okLoadClusters, nClusters;
fe17d4cb 353 if (!fClusters || (fClusters && fClusters->IsEmpty())) {
fe17d4cb 354 okLoadClusters = LoadClusters(esd);
fe17d4cb 355 }
c61f0e70 356 nClusters = fClusters->GetEntries();
5970dfe2 357
358 // step 2: collect ESD tracks
359 Int_t nTracks, okLoadTracks;
360 okLoadTracks = LoadTracks(esd);
c61f0e70 361 nTracks = fTracks->GetEntries();
fe17d4cb 362
5970dfe2 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());
952d023a 373 track->SetStatus(AliESDtrack::kEMCALmatch);
5970dfe2 374 }
375 }
fe17d4cb 376
377 return 0;
378}
fe17d4cb 379
380//
381//------------------------------------------------------------------------------
382//
5970dfe2 383Int_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 //
17773e2e 389
5970dfe2 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
8108ab9e 396 AliExternalTrackParam *trkParam = 0;
5970dfe2 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 {
8108ab9e 409 //AliExternalTrackParam *trkParamTmp = new AliExternalTrackParam(*trkParam);
410 AliExternalTrackParam trkParamTmp(*trkParam);
5970dfe2 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);
8108ab9e 416 trkParamTmp.Rotate(alpha);
5970dfe2 417 //extrapolation is done here
8108ab9e 418 if(!AliTrackerBase::PropagateTrackToBxByBz(&trkParamTmp, vec.X(), track->GetMass(), fStep, kFALSE, 0.8, -1)) continue;
5970dfe2 419
420 //Calculate the residuals
8108ab9e 421 trkParamTmp.GetXYZ(trkPos);
5970dfe2 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 }
8108ab9e 439 //delete trkParamTmp;
5970dfe2 440 }
441 return index;
fe17d4cb 442}
5970dfe2 443
fe17d4cb 444//
445//------------------------------------------------------------------------------
446//
447void AliEMCALTracker::UnloadClusters()
448{
449 //
044225d9 450 // Free memory from all arrays
451 // This method is called after the local tracking step
452 // so we can safely delete everything
fe17d4cb 453 //
454
044225d9 455 Clear();
fe17d4cb 456}
04475328 457
458//
459//------------------------------------------------------------------------------
460//
fe17d4cb 461AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliEMCALRecPoint *recPoint)
462 : fIndex(index),
fe17d4cb 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//
481AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliESDCaloCluster *caloCluster)
482 : fIndex(index),
fe17d4cb 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 //
53e430a3 491 Float_t clpos[3]= {0., 0., 0.};
7592dfc4 492 caloCluster->GetPosition(clpos);
fe17d4cb 493
494 fX = (Double_t)clpos[0];
495 fY = (Double_t)clpos[1];
496 fZ = (Double_t)clpos[2];
497}