Line data Source code
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 : #include <Riostream.h>
17 : #include <iomanip>
18 :
19 : #include <TFile.h>
20 : #include <TTree.h>
21 : #include <TList.h>
22 : #include <TString.h>
23 : #include <TVector3.h>
24 : #include <TClonesArray.h>
25 : #include <TGeoMatrix.h>
26 :
27 : #include "AliLog.h"
28 : #include "AliESDEvent.h"
29 : #include "AliESDtrack.h"
30 : #include "AliESDCaloCluster.h"
31 : #include "AliEMCALRecPoint.h"
32 : #include "AliRunLoader.h"
33 : #include "AliEMCALTrack.h"
34 : #include "AliEMCALLoader.h"
35 : #include "AliEMCALGeometry.h"
36 : #include "AliEMCALReconstructor.h"
37 : #include "AliEMCALRecParam.h"
38 : #include "AliCDBEntry.h"
39 : #include "AliCDBManager.h"
40 : #include "AliEMCALReconstructor.h"
41 : #include "AliEMCALRecoUtils.h"
42 :
43 : #include "AliEMCALTracker.h"
44 :
45 : /// \cond CLASSIMP
46 42 : ClassImp(AliEMCALTracker) ;
47 : /// \endcond
48 :
49 : ///
50 : /// Default constructor.
51 : /// Initializes all simple data members to default values,
52 : /// and all collections to NULL.
53 : /// Output file name is set to a default value.
54 : ///
55 : //------------------------------------------------------------------------------
56 : AliEMCALTracker::AliEMCALTracker() :
57 2 : AliTracker(),
58 2 : fCutPt(0),
59 2 : fCutNITS(0),
60 2 : fCutNTPC(50),
61 2 : fStep(20),
62 2 : fTrackCorrMode(kTrackCorrMMB),
63 2 : fEMCalSurfaceDistance(440),
64 2 : fClusterWindow(50),
65 2 : fCutEta(0.025),
66 2 : fCutPhi(0.05),
67 2 : fITSTrackSA(kFALSE),
68 2 : fTrackInITS(kFALSE),
69 2 : fTracks(0),
70 2 : fClusters(0),
71 2 : fGeom(0)
72 10 : {
73 2 : InitParameters();
74 4 : }
75 :
76 : ///
77 : /// Copy constructor.
78 : /// Besides copying all parameters, duplicates all collections.
79 : ///
80 : //------------------------------------------------------------------------------
81 : AliEMCALTracker::AliEMCALTracker(const AliEMCALTracker& copy) :
82 0 : AliTracker(),
83 0 : fCutPt(copy.fCutPt),
84 0 : fCutNITS(copy.fCutNITS),
85 0 : fCutNTPC(copy.fCutNTPC),
86 0 : fStep(copy.fStep),
87 0 : fTrackCorrMode(copy.fTrackCorrMode),
88 0 : fEMCalSurfaceDistance(copy.fEMCalSurfaceDistance),
89 0 : fClusterWindow(copy.fClusterWindow),
90 0 : fCutEta(copy.fCutEta),
91 0 : fCutPhi(copy.fCutPhi),
92 0 : fITSTrackSA(copy.fITSTrackSA),
93 0 : fTrackInITS(copy.fTrackInITS),
94 0 : fTracks((TObjArray*)copy.fTracks->Clone()),
95 0 : fClusters((TObjArray*)copy.fClusters->Clone()),
96 0 : fGeom(copy.fGeom)
97 0 : {
98 0 : }
99 :
100 : ///
101 : /// Assignment operator; use copy ctor.
102 : ///
103 : //------------------------------------------------------------------------------
104 : AliEMCALTracker& AliEMCALTracker::operator=(const AliEMCALTracker& source)
105 : {
106 0 : if (&source == this) return *this;
107 :
108 0 : new (this) AliEMCALTracker(source);
109 0 : return *this;
110 0 : }
111 :
112 : ///
113 : /// Retrieve initialization parameters.
114 : ///
115 : //------------------------------------------------------------------------------
116 : void AliEMCALTracker::InitParameters()
117 : {
118 : // Check if the instance of AliEMCALRecParam exists,
119 4 : const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
120 :
121 2 : if (!recParam)
122 : {
123 0 : AliFatal("Reconstruction parameters for EMCAL not set!");
124 0 : return; // not needed, but for coverity.
125 : }
126 :
127 2 : fCutEta = recParam->GetMthCutEta();
128 2 : fCutPhi = recParam->GetMthCutPhi();
129 2 : fStep = recParam->GetExtrapolateStep();
130 2 : fCutPt = recParam->GetTrkCutPt();
131 2 : fCutNITS = recParam->GetTrkCutNITS(); // Not used?
132 2 : fCutNTPC = recParam->GetTrkCutNTPC();
133 2 : fTrackInITS = recParam->GetTrkInITS();
134 4 : }
135 :
136 : ///
137 : /// Clearing method.
138 : /// Deletes all objects in arrays and the arrays themselves.
139 : ///
140 : //------------------------------------------------------------------------------
141 : void AliEMCALTracker::Clear(Option_t* option)
142 : {
143 52 : TString opt(option);
144 52 : Bool_t clearTracks = opt.Contains("TRACKS");
145 52 : Bool_t clearClusters = opt.Contains("CLUSTERS");
146 52 : if (opt.Contains("ALL"))
147 : {
148 : clearTracks = kTRUE;
149 : clearClusters = kTRUE;
150 10 : }
151 :
152 : // fTracks is a collection of esdTrack
153 : // When clearing this array, the linked objects should not be deleted
154 34 : if (fTracks != 0x0 && clearTracks)
155 : {
156 8 : fTracks->Clear();
157 16 : delete fTracks;
158 8 : fTracks = 0;
159 8 : }
160 :
161 42 : if (fClusters != 0x0 && clearClusters)
162 : {
163 8 : fClusters->Delete();
164 16 : delete fClusters;
165 8 : fClusters = 0;
166 8 : }
167 26 : }
168 :
169 : ///
170 : /// Load EMCAL clusters in the form of AliEMCALRecPoint,
171 : /// from simulation temporary files.
172 : /// (When included in reconstruction chain, this method is used automatically).
173 : ///
174 : //------------------------------------------------------------------------------
175 : Int_t AliEMCALTracker::LoadClusters(TTree *cTree)
176 : {
177 16 : Clear("CLUSTERS");
178 :
179 8 : cTree->SetBranchStatus("*",0); //disable all branches
180 8 : cTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
181 :
182 8 : TBranch *branch = cTree->GetBranch("EMCALECARP");
183 8 : if (!branch)
184 : {
185 0 : AliError("Can't get the branch with the EMCAL clusters");
186 0 : return 1;
187 : }
188 :
189 16 : TClonesArray *clusters = new TClonesArray("AliEMCALRecPoint", 1000);
190 8 : branch->SetAddress(&clusters);
191 :
192 : //cTree->GetEvent(0);
193 8 : branch->GetEntry(0);
194 :
195 8 : Int_t nClusters = (Int_t)clusters->GetEntries();
196 8 : if (fClusters) fClusters->Delete();
197 16 : else fClusters = new TObjArray(0);
198 :
199 82 : for (Int_t i = 0; i < nClusters; i++)
200 : {
201 33 : AliEMCALRecPoint *cluster = (AliEMCALRecPoint*)clusters->At(i);
202 33 : if (!cluster) continue;
203 :
204 33 : AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
205 33 : fClusters->AddLast(matchCluster);
206 33 : }
207 :
208 8 : branch->SetAddress(0);
209 8 : clusters->Delete();
210 16 : delete clusters;
211 :
212 24 : AliDebug(1,Form("Collected %d RecPoints from Tree", fClusters->GetEntries()));
213 :
214 : return 0;
215 16 : }
216 :
217 : ///
218 : /// Load EMCAL clusters in the form of AliESDCaloClusters,
219 : /// from an AliESD object.
220 : //
221 : //------------------------------------------------------------------------------
222 : Int_t AliEMCALTracker::LoadClusters(AliESDEvent *esd)
223 : {
224 : // make sure that tracks/clusters collections are empty
225 0 : Clear("CLUSTERS");
226 0 : fClusters = new TObjArray(0);
227 :
228 0 : Int_t nClusters = esd->GetNumberOfCaloClusters();
229 0 : for (Int_t i=0; i<nClusters; i++)
230 : {
231 0 : AliESDCaloCluster *cluster = esd->GetCaloCluster(i);
232 0 : if (!cluster || !cluster->IsEMCAL()) continue ;
233 :
234 0 : AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
235 0 : fClusters->AddLast(matchCluster);
236 0 : }
237 :
238 0 : AliDebug(1,Form("Collected %d clusters from ESD", fClusters->GetEntries()));
239 0 : return 0;
240 0 : }
241 :
242 : ///
243 : /// Load ESD tracks.
244 : /// Filter them depending on loose track and acceptance cuts.
245 : ///
246 : /// If TPC is not present, try matching with ITS tracks only.
247 : /// Set the flag fITSTrackSA.
248 : ///
249 : //------------------------------------------------------------------------------
250 : Int_t AliEMCALTracker::LoadTracks(AliESDEvent *esd)
251 : {
252 : // In case of runs without TPC but ITS on,
253 : // do matching with different cuts
254 16 : UInt_t mask1 = esd->GetESDRun()->GetDetectorsInDAQ();
255 8 : UInt_t mask2 = esd->GetESDRun()->GetDetectorsInReco();
256 8 : Bool_t desc1 = (mask1 >> 3) & 0x1;
257 8 : Bool_t desc2 = (mask2 >> 3) & 0x1;
258 :
259 16 : if (desc1==0 || desc2==0)
260 : {
261 : // AliError(Form("TPC not in DAQ/RECO: %u (%u)/%u (%u)",
262 : // mask1, esd->GetESDRun()->GetDetectorsInReco(),
263 : // mask2, esd->GetESDRun()->GetDetectorsInDAQ()));
264 0 : fITSTrackSA = kTRUE;
265 0 : }
266 :
267 8 : Clear("TRACKS");
268 16 : fTracks = new TObjArray(0);
269 :
270 8 : Int_t nTracks = esd->GetNumberOfTracks();
271 : //printf("N tracks %d\n",nTracks);
272 :
273 320 : for (Int_t i = 0; i < nTracks; i++)
274 : {
275 152 : AliESDtrack *esdTrack = esd->GetTrack(i);
276 :
277 : // Set by default the value corresponding to "no match"
278 152 : esdTrack->SetEMCALcluster(kUnmatched);
279 152 : esdTrack->ResetStatus(AliESDtrack::kEMCALmatch);
280 :
281 : //
282 : // Select tracks depending on status bit
283 : // and number of clusters.
284 : //
285 152 : if ( esdTrack->Pt() < fCutPt ) continue;
286 :
287 : // When there are TPC tracks, select them depeding on the fit bit
288 : // Take TPCout and if requested ITSout (default)
289 152 : if ( !fITSTrackSA )
290 : {
291 152 : Int_t nClustersITS = esdTrack->GetITSclusters(0);
292 152 : Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
293 :
294 152 : ULong_t status = esdTrack->GetStatus();
295 :
296 : //Bool_t tpcIn = (status&AliESDtrack::kTPCin )==AliESDtrack::kTPCin ;
297 : //Bool_t itsIn = (status&AliESDtrack::kITSin )==AliESDtrack::kITSin ;
298 152 : Bool_t tpcOut = (status&AliESDtrack::kTPCout)==AliESDtrack::kTPCout;
299 152 : Bool_t itsOut = (status&AliESDtrack::kITSout)==AliESDtrack::kITSout;
300 :
301 : // Check ITSout bit, if requested
302 152 : if ( fTrackInITS )
303 : {
304 : // if(nClustersITS > 0) // fCutNITS?
305 : // printf("ITS clusters %d, Cut %2.0f, kITSin? %d, kITSout? %d - select %d\n",
306 : // nClustersITS,fCutNITS,itsIn,itsOut,fTrackInITS);
307 :
308 0 : if ( !itsOut || nClustersITS <= 0 ) continue; // fCutNITS
309 : }
310 :
311 : // if(nClustersTPC > 0)
312 : // printf("TPC clusters %d, Cut %2.0f, kTPCin? %d; kTPCout? %d\n",
313 : // nClustersTPC,fCutNTPC,tpcIn,tpcOut);
314 :
315 304 : if ( !tpcOut || nClustersTPC < fCutNTPC ) continue;
316 :
317 134 : } // Do not match with ITS SA track
318 :
319 : //printf("\t accepted!\n");
320 :
321 : //
322 : // Loose geometric cut
323 : //
324 136 : if ( TMath::Abs(esdTrack->Eta()) > 0.9 ) continue;
325 :
326 : // Save some time and memory in case of no DCal present
327 :
328 132 : if(fGeom->GetNumberOfSuperModules() < 13)
329 : {
330 66 : Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
331 92 : if ( phi <= 10 || phi >= 250 ) continue;
332 40 : }
333 :
334 106 : fTracks->AddLast(esdTrack);
335 106 : }
336 :
337 24 : AliDebug(1,Form("Collected %d tracks", fTracks->GetEntries()));
338 8 : return 0;
339 0 : }
340 :
341 : ///
342 : /// Set track correction mode
343 : /// gets the choice in string format and converts into
344 : /// internal enum.
345 : ///
346 : //------------------------------------------------------------------------------
347 : void AliEMCALTracker::SetTrackCorrectionMode(Option_t *option)
348 : {
349 0 : TString opt(option);
350 0 : opt.ToUpper();
351 :
352 0 : if ( !opt.CompareTo("NONE") )
353 0 : fTrackCorrMode = kTrackCorrNone;
354 0 : else if ( !opt.CompareTo("MMB") )
355 0 : fTrackCorrMode = kTrackCorrMMB;
356 : else
357 0 : AliError("Unrecognized option");
358 0 : }
359 :
360 : ///
361 : /// Main operation method.
362 : /// Gets external AliESD containing tracks to be matched.
363 : /// After executing match finding, stores in the same ESD object all infos
364 : /// and releases the object for further reconstruction steps.
365 : ///
366 : /// Note: should always return 0=OK, because otherwise all tracking
367 : /// is aborted for this event.
368 : //------------------------------------------------------------------------------
369 : Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
370 : {
371 16 : if (!esd)
372 : {
373 0 : AliError("NULL ESD passed");
374 0 : return 1;
375 : }
376 :
377 10 : if( !fGeom ) fGeom = AliEMCALGeometry::GetInstance();
378 :
379 : // step 1: collect clusters
380 : Int_t okLoadClusters = 0;
381 24 : if (!fClusters || (fClusters && fClusters->IsEmpty()))
382 0 : okLoadClusters = LoadClusters(esd);
383 :
384 8 : Int_t nClusters = fClusters->GetEntries();
385 :
386 : // step 2: collect ESD tracks
387 : Int_t nTracks, okLoadTracks;
388 8 : okLoadTracks = LoadTracks(esd);
389 8 : nTracks = fTracks->GetEntries();
390 :
391 24 : AliDebug(5,Form("Propagate back %d tracks ok %d, for %d clusters ok %d",
392 : nTracks,okLoadTracks,nClusters,okLoadClusters));
393 :
394 : // step 3: for each track, find the closest cluster as matched within residual cuts
395 : Int_t index=-1;
396 228 : for (Int_t it = 0; it < nTracks; it++)
397 : {
398 106 : AliESDtrack *track = (AliESDtrack*)fTracks->At(it);
399 106 : index = FindMatchedCluster(track);
400 106 : if (index>-1)
401 : {
402 18 : AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(index);
403 18 : track->SetEMCALcluster(cluster->Index());
404 18 : track->SetStatus(AliESDtrack::kEMCALmatch);
405 18 : }
406 : }
407 :
408 : return 0;
409 8 : }
410 :
411 : ///
412 : /// For each track, extrapolate it to all the clusters.
413 : /// Find the closest one as matched if the residuals (dEta, dPhi) satisfy the cuts.
414 : ///
415 : //------------------------------------------------------------------------------
416 : Int_t AliEMCALTracker::FindMatchedCluster(AliESDtrack *track)
417 : {
418 212 : Float_t maxEta=fCutEta;
419 106 : Float_t maxPhi=fCutPhi;
420 : Int_t index = -1;
421 :
422 : // If the esdFriend is available, use the TPCOuter point as the starting point of extrapolation
423 : // Otherwise use the TPCInner point
424 : AliExternalTrackParam *trkParam = 0;
425 :
426 106 : if (!fITSTrackSA)
427 : {
428 106 : const AliESDfriendTrack* friendTrack = track->GetFriendTrack();
429 :
430 212 : if (friendTrack && friendTrack->GetTPCOut())
431 106 : trkParam = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
432 0 : else if (track->GetInnerParam())
433 0 : trkParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
434 106 : }
435 : else
436 0 : trkParam = new AliExternalTrackParam(*track);
437 :
438 106 : if (!trkParam) return index;
439 :
440 106 : AliExternalTrackParam trkParamTmp(*trkParam);
441 106 : Float_t eta, phi, pt;
442 318 : if (!AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&trkParamTmp, fEMCalSurfaceDistance, track->GetMass(kTRUE), fStep, eta, phi, pt))
443 : {
444 25 : if (fITSTrackSA) delete trkParam;
445 25 : return index;
446 : }
447 :
448 81 : track->SetTrackPhiEtaPtOnEMCal(phi,eta,pt);
449 :
450 81 : if ( TMath::Abs(eta) > 0.75 )
451 : {
452 0 : if (fITSTrackSA) delete trkParam;
453 0 : return index;
454 : }
455 :
456 : // Save some time and memory in case of no DCal present
457 177 : if ( fGeom->GetNumberOfSuperModules() < 13 &&
458 46 : ( phi < 70*TMath::DegToRad() || phi > 190*TMath::DegToRad()))
459 : {
460 19 : if (fITSTrackSA) delete trkParam;
461 19 : return index;
462 : }
463 :
464 : //Perform extrapolation
465 62 : Double_t trkPos[3];
466 62 : trkParamTmp.GetXYZ(trkPos);
467 62 : Int_t nclusters = fClusters->GetEntries();
468 640 : for (Int_t ic=0; ic<nclusters; ic++)
469 : {
470 516 : AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
471 :
472 774 : Float_t clsPos[3] = {static_cast<Float_t>(cluster->X()),
473 258 : static_cast<Float_t>(cluster->Y()),
474 258 : static_cast<Float_t>(cluster->Z())};
475 :
476 258 : 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));
477 : //printf("\n dR=%f,wind=%f\n",dR,fClusterWindow); //MARCEL
478 :
479 478 : if (dR > fClusterWindow) continue;
480 :
481 38 : AliExternalTrackParam trkParTmp(trkParamTmp);
482 :
483 38 : Float_t tmpEta, tmpPhi;
484 116 : if (!AliEMCALRecoUtils::ExtrapolateTrackToPosition(&trkParTmp, clsPos,track->GetMass(kTRUE), 5, tmpEta, tmpPhi)) continue;
485 :
486 60 : if (TMath::Abs(tmpPhi)<TMath::Abs(maxPhi) && TMath::Abs(tmpEta)<TMath::Abs(maxEta))
487 : {
488 22 : maxPhi=tmpPhi;
489 22 : maxEta=tmpEta;
490 : index=ic;
491 22 : }
492 332 : }
493 :
494 62 : if (fITSTrackSA) delete trkParam;
495 :
496 : return index;
497 274 : }
498 :
499 : ///
500 : /// Free memory from all arrays.
501 : /// This method is called after the local tracking step.
502 : /// so we can safely delete everything.
503 : ///
504 : //------------------------------------------------------------------------------
505 : void AliEMCALTracker::UnloadClusters()
506 : {
507 16 : Clear();
508 8 : }
509 :
510 : ///
511 : /// Translates an AliEMCALRecPoint object into the internal format.
512 : /// Index of passed cluster in its native array must be specified.
513 : ///
514 : //------------------------------------------------------------------------------
515 33 : AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliEMCALRecPoint *recPoint) :
516 33 : fIndex(index),
517 33 : fX(0.),
518 33 : fY(0.),
519 33 : fZ(0.)
520 165 : {
521 33 : TVector3 clpos;
522 33 : recPoint->GetGlobalPosition(clpos);
523 :
524 33 : fX = clpos.X();
525 33 : fY = clpos.Y();
526 33 : fZ = clpos.Z();
527 66 : }
528 :
529 : ///
530 : /// Translates an AliESDCaloCluster object into the internal format.
531 : /// Index of passed cluster in its native array must be specified.
532 : ///
533 : //------------------------------------------------------------------------------
534 0 : AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliESDCaloCluster *caloCluster) :
535 0 : fIndex(index),
536 0 : fX(0.),
537 0 : fY(0.),
538 0 : fZ(0.)
539 0 : {
540 0 : Float_t clpos[3]= {0., 0., 0.};
541 0 : caloCluster->GetPosition(clpos);
542 :
543 0 : fX = (Double_t)clpos[0];
544 0 : fY = (Double_t)clpos[1];
545 0 : fZ = (Double_t)clpos[2];
546 0 : }
|
