1 //========================================================================
2 // Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved.
4 // Author: The ALICE Off-line Project.
5 // Contributors are mentioned in the code where appropriate.
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 //========================================================================
16 // Class AliEMCALTracker
17 // -----------------------
18 // Implementation of the track matching method between barrel tracks and
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().
28 // ------------------------------------------------------------------------
29 // author: A. Pulvirenti (alberto.pulvirenti@ct.infn.it)
30 // Revised by Rongrong 2010-05-31 (rongrong.ma@cern.ch)
31 //=========================================================================
33 #include <Riostream.h>
41 #include <TClonesArray.h>
42 #include <TGeoMatrix.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"
59 #include "AliEMCALTracker.h"
61 ClassImp(AliEMCALTracker)
64 //------------------------------------------------------------------------------
66 AliEMCALTracker::AliEMCALTracker()
72 fTrackCorrMode(kTrackCorrMMB),
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.
88 //------------------------------------------------------------------------------
90 AliEMCALTracker::AliEMCALTracker(const AliEMCALTracker& copy)
93 fCutNITS(copy.fCutNITS),
94 fCutNTPC(copy.fCutNTPC),
96 fTrackCorrMode(copy.fTrackCorrMode),
97 fCutEta(copy.fCutEta),
98 fCutPhi(copy.fCutPhi),
99 fTracks((TObjArray*)copy.fTracks->Clone()),
100 fClusters((TObjArray*)copy.fClusters->Clone()),
105 // Besides copying all parameters, duplicates all collections.
109 //------------------------------------------------------------------------------
111 AliEMCALTracker& AliEMCALTracker::operator=(const AliEMCALTracker& copy)
114 // Assignment operator.
115 // Besides copying all parameters, duplicates all collections.
118 fCutPt = copy.fCutPt;
119 fCutEta = copy.fCutEta;
120 fCutPhi = copy.fCutPhi;
122 fTrackCorrMode = copy.fTrackCorrMode;
124 fCutNITS = copy.fCutNITS;
125 fCutNTPC = copy.fCutNTPC;
127 fTracks = (TObjArray*)copy.fTracks->Clone();
128 fClusters = (TObjArray*)copy.fClusters->Clone();
134 //------------------------------------------------------------------------------
136 void AliEMCALTracker::InitParameters()
139 // Retrieve initialization parameters
142 // Check if the instance of AliEMCALRecParam exists,
143 const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
146 AliFatal("Reconstruction parameters for EMCAL not set!");
150 fCutEta = recParam->GetMthCutEta();
151 fCutPhi = recParam->GetMthCutPhi();
152 fStep = recParam->GetExtrapolateStep();
153 fCutPt = recParam->GetTrkCutPt();
154 fCutNITS = recParam->GetTrkCutNITS();
155 fCutNTPC = recParam->GetTrkCutNTPC();
161 //------------------------------------------------------------------------------
163 void AliEMCALTracker::Clear(Option_t* option)
167 // Deletes all objects in arrays and the arrays themselves
171 Bool_t clearTracks = opt.Contains("TRACKS");
172 Bool_t clearClusters = opt.Contains("CLUSTERS");
173 if (opt.Contains("ALL")) {
175 clearClusters = kTRUE;
178 //fTracks is a collection of esdTrack
179 //When clearing this array, the linked objects should not be deleted
180 if (fTracks != 0x0 && clearTracks) {
185 if (fClusters != 0x0 && clearClusters) {
192 //------------------------------------------------------------------------------
194 Int_t AliEMCALTracker::LoadClusters(TTree *cTree)
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)
204 cTree->SetBranchStatus("*",0); //disable all branches
205 cTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
207 TBranch *branch = cTree->GetBranch("EMCALECARP");
209 AliError("Can't get the branch with the EMCAL clusters");
213 TClonesArray *clusters = new TClonesArray("AliEMCALRecPoint", 1000);
214 branch->SetAddress(&clusters);
216 //cTree->GetEvent(0);
218 Int_t nClusters = (Int_t)clusters->GetEntries();
219 if(fClusters) fClusters->Delete();
220 else fClusters = new TObjArray(0);
221 for (Int_t i = 0; i < nClusters; i++) {
222 AliEMCALRecPoint *cluster = (AliEMCALRecPoint*)clusters->At(i);
223 if (!cluster) continue;
224 AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
225 fClusters->AddLast(matchCluster);
228 branch->SetAddress(0);
232 AliInfo(Form("Collected %d RecPoints from Tree", fClusters->GetEntries()));
237 //------------------------------------------------------------------------------
239 Int_t AliEMCALTracker::LoadClusters(AliESDEvent *esd)
242 // Load EMCAL clusters in the form of AliESDCaloClusters,
243 // from an AliESD object.
246 // make sure that tracks/clusters collections are empty
248 fClusters = new TObjArray(0);
250 Int_t nClusters = esd->GetNumberOfCaloClusters();
251 for (Int_t i=0; i<nClusters; i++)
253 AliESDCaloCluster *cluster = esd->GetCaloCluster(i);
254 if (!cluster || !cluster->IsEMCAL()) continue ;
255 AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
256 fClusters->AddLast(matchCluster);
259 AliInfo(Form("Collected %d clusters from ESD", fClusters->GetEntries()));
263 //------------------------------------------------------------------------------
265 Int_t AliEMCALTracker::LoadTracks(AliESDEvent *esd)
272 fTracks = new TObjArray(0);
274 Int_t nTracks = esd->GetNumberOfTracks();
275 Bool_t isKink=kFALSE;
276 for (Int_t i = 0; i < nTracks; i++)
278 AliESDtrack *esdTrack = esd->GetTrack(i);
279 // set by default the value corresponding to "no match"
280 esdTrack->SetEMCALcluster(kUnmatched);
282 //Select good quaulity tracks
283 if(esdTrack->Pt()<fCutPt) continue;
284 if(esdTrack->GetNcls(1)<fCutNTPC)continue;
286 //Reject kink daughters
288 for (Int_t j = 0; j < 3; j++)
290 if (esdTrack->GetKinkIndex(j) != 0) isKink = kTRUE;
292 if (isKink) continue;
294 //Loose geometric cut
295 Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
296 if(TMath::Abs(esdTrack->Eta())>1 || phi <= 60 || phi >= 200 ) continue;
298 fTracks->AddLast(esdTrack);
301 AliInfo(Form("Collected %d tracks", fTracks->GetEntries()));
305 //------------------------------------------------------------------------------
307 void AliEMCALTracker::SetTrackCorrectionMode(Option_t *option)
310 // Set track correction mode
311 // gest the choice in string format and converts into
318 if (!opt.CompareTo("NONE"))
320 fTrackCorrMode = kTrackCorrNone;
322 else if (!opt.CompareTo("MMB"))
324 fTrackCorrMode = kTrackCorrMMB;
328 cerr << "E-AliEMCALTracker::SetTrackCorrectionMode '" << option << "': Unrecognized option" << endl;
332 //------------------------------------------------------------------------------
334 Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
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.
343 // Note: should always return 0=OK, because otherwise all tracking
344 // is aborted for this event
347 AliError("NULL ESD passed");
351 // step 1: collect clusters
352 Int_t okLoadClusters, nClusters;
353 if (!fClusters || (fClusters && fClusters->IsEmpty())) {
354 okLoadClusters = LoadClusters(esd);
356 nClusters = fClusters->GetEntries();
358 // step 2: collect ESD tracks
359 Int_t nTracks, okLoadTracks;
360 okLoadTracks = LoadTracks(esd);
361 nTracks = fTracks->GetEntries();
363 // step 3: for each track, find the closest cluster as matched within residual cuts
365 for (Int_t it = 0; it < nTracks; it++)
367 AliESDtrack *track = (AliESDtrack*)fTracks->At(it);
368 index = FindMatchedCluster(track);
371 AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(index);
372 track->SetEMCALcluster(cluster->Index());
380 //------------------------------------------------------------------------------
382 Int_t AliEMCALTracker::FindMatchedCluster(AliESDtrack *track)
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
389 Double_t maxEta=fCutEta;
390 Double_t maxPhi=fCutPhi;
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;
396 const AliESDfriendTrack* friendTrack = track->GetFriendTrack();
397 if(friendTrack && friendTrack->GetTPCOut())
398 trkParam = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
400 trkParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
401 if(!trkParam) return index;
403 //Perform extrapolation
405 Int_t nclusters = fClusters->GetEntries();
406 for(Int_t ic=0; ic<nclusters; ic++)
408 AliExternalTrackParam *trkParamTmp = new AliExternalTrackParam(*trkParam);
409 AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
410 TVector3 vec(cluster->X(),cluster->Y(),cluster->Z());
411 Double_t alpha = ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
412 //Rotate to proper local coordinate
414 trkParamTmp->Rotate(alpha);
415 //extrapolation is done here
416 if(!AliTrackerBase::PropagateTrackToBxByBz(trkParamTmp, vec.X(), track->GetMass(), fStep, kFALSE, 0.8, -1)) continue;
418 //Calculate the residuals
419 trkParamTmp->GetXYZ(trkPos);
420 TVector3 clsPosVec(cluster->X(),cluster->Y(),cluster->Z());
421 TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
423 Double_t clsPhi = clsPosVec.Phi();
424 if(clsPhi<0) clsPhi+=2*TMath::Pi();
425 Double_t trkPhi = trkPosVec.Phi();
426 if(trkPhi<0) trkPhi+=2*TMath::Pi();
428 Double_t tmpPhi = clsPhi-trkPhi;
429 Double_t tmpEta = clsPosVec.Eta()-trkPosVec.Eta();
431 if(TMath::Abs(tmpPhi)<TMath::Abs(maxPhi) && TMath::Abs(tmpEta)<TMath::Abs(maxEta))
443 //------------------------------------------------------------------------------
445 void AliEMCALTracker::UnloadClusters()
448 // Free memory from all arrays
449 // This method is called after the local tracking step
450 // so we can safely delete everything
457 //------------------------------------------------------------------------------
459 AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliEMCALRecPoint *recPoint)
466 // Translates an AliEMCALRecPoint object into the internal format.
467 // Index of passed cluster in its native array must be specified.
470 recPoint->GetGlobalPosition(clpos);
477 //------------------------------------------------------------------------------
479 AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliESDCaloCluster *caloCluster)
486 // Translates an AliESDCaloCluster object into the internal format.
487 // Index of passed cluster in its native array must be specified.
489 Float_t clpos[3]= {0., 0., 0.};
490 caloCluster->GetPosition(clpos);
492 fX = (Double_t)clpos[0];
493 fY = (Double_t)clpos[1];
494 fZ = (Double_t)clpos[2];