correct previous non intended commit, but leaving some new things
[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
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
5970dfe2 276 Clear("TRACKS");
277 fTracks = new TObjArray(0);
fe17d4cb 278
5970dfe2 279 Int_t nTracks = esd->GetNumberOfTracks();
8fc351e3 280 //Bool_t isKink=kFALSE;
5970dfe2 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);
8fc351e3 286 esdTrack->ResetStatus(AliESDtrack::kEMCALmatch);
5970dfe2 287
288 //Select good quaulity tracks
289 if(esdTrack->Pt()<fCutPt) continue;
42ceff04 290 if(!fITSTrackSA)
291 if(esdTrack->GetNcls(1)<fCutNTPC)continue;
5970dfe2 292
5970dfe2 293 //Loose geometric cut
294 Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
d874f457 295 if(TMath::Abs(esdTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
5970dfe2 296
297 fTracks->AddLast(esdTrack);
298 }
299
300 AliInfo(Form("Collected %d tracks", fTracks->GetEntries()));
301 return 0;
302}
303//
304//------------------------------------------------------------------------------
305//
306void 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 }
fe17d4cb 329}
330//
331//------------------------------------------------------------------------------
332//
af885e0f 333Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
fe17d4cb 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 //
3d9e8b15 341 //
342 // Note: should always return 0=OK, because otherwise all tracking
343 // is aborted for this event
5970dfe2 344
fe17d4cb 345 if (!esd) {
346 AliError("NULL ESD passed");
347 return 1;
348 }
349
5970dfe2 350 // step 1: collect clusters
c61f0e70 351 Int_t okLoadClusters, nClusters;
fe17d4cb 352 if (!fClusters || (fClusters && fClusters->IsEmpty())) {
fe17d4cb 353 okLoadClusters = LoadClusters(esd);
fe17d4cb 354 }
c61f0e70 355 nClusters = fClusters->GetEntries();
5970dfe2 356
357 // step 2: collect ESD tracks
358 Int_t nTracks, okLoadTracks;
359 okLoadTracks = LoadTracks(esd);
c61f0e70 360 nTracks = fTracks->GetEntries();
fe17d4cb 361
5970dfe2 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());
952d023a 372 track->SetStatus(AliESDtrack::kEMCALmatch);
5970dfe2 373 }
374 }
fe17d4cb 375
376 return 0;
377}
fe17d4cb 378
379//
380//------------------------------------------------------------------------------
381//
5970dfe2 382Int_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 //
17773e2e 388
ee602376 389 Float_t maxEta=fCutEta;
390 Float_t maxPhi=fCutPhi;
5970dfe2 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
8108ab9e 395 AliExternalTrackParam *trkParam = 0;
42ceff04 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 }
5970dfe2 404 else
42ceff04 405 trkParam = new AliExternalTrackParam(*track);
406
5970dfe2 407 if(!trkParam) return index;
42ceff04 408
8fc351e3 409 AliExternalTrackParam trkParamTmp(*trkParam);
ee602376 410 Float_t eta, phi;
da25bc37 411 if(!AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&trkParamTmp, 430., track->GetMass(kTRUE), fStep, eta, phi)) return index;
ed15417e 412 track->SetTrackPhiEtaOnEMCal(phi,eta);
8fc351e3 413 if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()) return index;
414
5970dfe2 415 //Perform extrapolation
416 Double_t trkPos[3];
8fc351e3 417 trkParamTmp.GetXYZ(trkPos);
5970dfe2 418 Int_t nclusters = fClusters->GetEntries();
419 for(Int_t ic=0; ic<nclusters; ic++)
420 {
5970dfe2 421 AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
8fc351e3 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));
42ceff04 424// printf("\n dR=%f,wind=%f\n",dR,fClusterWindow); //MARCEL
8fc351e3 425 if(dR > fClusterWindow) continue;
426
427 AliExternalTrackParam trkParTmp(trkParamTmp);
428
ee602376 429 Float_t tmpEta, tmpPhi;
da25bc37 430 if(!AliEMCALRecoUtils::ExtrapolateTrackToPosition(&trkParTmp, clsPos,track->GetMass(kTRUE), 5, tmpEta, tmpPhi)) continue;
5970dfe2 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 }
5970dfe2 437 }
438 return index;
fe17d4cb 439}
5970dfe2 440
fe17d4cb 441//
442//------------------------------------------------------------------------------
443//
444void AliEMCALTracker::UnloadClusters()
445{
446 //
044225d9 447 // Free memory from all arrays
448 // This method is called after the local tracking step
449 // so we can safely delete everything
fe17d4cb 450 //
451
044225d9 452 Clear();
fe17d4cb 453}
04475328 454
455//
456//------------------------------------------------------------------------------
457//
fe17d4cb 458AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliEMCALRecPoint *recPoint)
459 : fIndex(index),
fe17d4cb 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//
478AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliESDCaloCluster *caloCluster)
479 : fIndex(index),
fe17d4cb 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 //
53e430a3 488 Float_t clpos[3]= {0., 0., 0.};
7592dfc4 489 caloCluster->GetPosition(clpos);
fe17d4cb 490
491 fX = (Double_t)clpos[0];
492 fY = (Double_t)clpos[1];
493 fZ = (Double_t)clpos[2];
494}