]>
Commit | Line | Data |
---|---|---|
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) | |
30 | // Revised by Rongrong 2010-05-31 (rongrong.ma@cern.ch) | |
31 | //========================================================================= | |
32 | ||
33 | #include <Riostream.h> | |
34 | #include <iomanip> | |
35 | ||
36 | #include <TFile.h> | |
37 | #include <TTree.h> | |
38 | #include <TList.h> | |
39 | #include <TString.h> | |
40 | #include <TVector3.h> | |
41 | #include <TClonesArray.h> | |
42 | #include <TGeoMatrix.h> | |
43 | ||
44 | #include "AliLog.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" | |
58 | #include "AliEMCALRecoUtils.h" | |
59 | ||
60 | #include "AliEMCALTracker.h" | |
61 | ||
62 | using std::cerr; | |
63 | using std::endl; | |
64 | ClassImp(AliEMCALTracker) | |
65 | ||
66 | // | |
67 | //------------------------------------------------------------------------------ | |
68 | // | |
69 | AliEMCALTracker::AliEMCALTracker() | |
70 | : AliTracker(), | |
71 | fCutPt(0), | |
72 | fCutNITS(0), | |
73 | fCutNTPC(50), | |
74 | fStep(20), | |
75 | fTrackCorrMode(kTrackCorrMMB), | |
76 | fEMCalSurfaceDistance(440), | |
77 | fClusterWindow(50), | |
78 | fCutEta(0.025), | |
79 | fCutPhi(0.05), | |
80 | fITSTrackSA(kFALSE), | |
81 | fTracks(0), | |
82 | fClusters(0), | |
83 | fGeom(0) | |
84 | { | |
85 | // | |
86 | // Default constructor. | |
87 | // Initializes all simple data members to default values, | |
88 | // and all collections to NULL. | |
89 | // Output file name is set to a default value. | |
90 | // | |
91 | InitParameters(); | |
92 | } | |
93 | // | |
94 | //------------------------------------------------------------------------------ | |
95 | // | |
96 | AliEMCALTracker::AliEMCALTracker(const AliEMCALTracker& copy) | |
97 | : AliTracker(), | |
98 | fCutPt(copy.fCutPt), | |
99 | fCutNITS(copy.fCutNITS), | |
100 | fCutNTPC(copy.fCutNTPC), | |
101 | fStep(copy.fStep), | |
102 | fTrackCorrMode(copy.fTrackCorrMode), | |
103 | fEMCalSurfaceDistance(copy.fEMCalSurfaceDistance), | |
104 | fClusterWindow(copy.fClusterWindow), | |
105 | fCutEta(copy.fCutEta), | |
106 | fCutPhi(copy.fCutPhi), | |
107 | fITSTrackSA(copy.fITSTrackSA), | |
108 | fTracks((TObjArray*)copy.fTracks->Clone()), | |
109 | fClusters((TObjArray*)copy.fClusters->Clone()), | |
110 | fGeom(copy.fGeom) | |
111 | { | |
112 | // | |
113 | // Copy constructor | |
114 | // Besides copying all parameters, duplicates all collections. | |
115 | // | |
116 | } | |
117 | // | |
118 | //------------------------------------------------------------------------------ | |
119 | // | |
120 | AliEMCALTracker& AliEMCALTracker::operator=(const AliEMCALTracker& source) | |
121 | { // assignment operator; use copy ctor | |
122 | if (&source == this) return *this; | |
123 | ||
124 | new (this) AliEMCALTracker(source); | |
125 | return *this; | |
126 | } | |
127 | // | |
128 | //------------------------------------------------------------------------------ | |
129 | // | |
130 | void AliEMCALTracker::InitParameters() | |
131 | { | |
132 | // | |
133 | // Retrieve initialization parameters | |
134 | // | |
135 | ||
136 | // Check if the instance of AliEMCALRecParam exists, | |
137 | const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam(); | |
138 | ||
139 | if(!recParam){ | |
140 | AliFatal("Reconstruction parameters for EMCAL not set!"); | |
141 | } | |
142 | else{ | |
143 | ||
144 | fCutEta = recParam->GetMthCutEta(); | |
145 | fCutPhi = recParam->GetMthCutPhi(); | |
146 | fStep = recParam->GetExtrapolateStep(); | |
147 | fCutPt = recParam->GetTrkCutPt(); | |
148 | fCutNITS = recParam->GetTrkCutNITS(); | |
149 | fCutNTPC = recParam->GetTrkCutNTPC(); | |
150 | } | |
151 | ||
152 | } | |
153 | ||
154 | // | |
155 | //------------------------------------------------------------------------------ | |
156 | // | |
157 | void AliEMCALTracker::Clear(Option_t* option) | |
158 | { | |
159 | // | |
160 | // Clearing method | |
161 | // Deletes all objects in arrays and the arrays themselves | |
162 | // | |
163 | ||
164 | TString opt(option); | |
165 | Bool_t clearTracks = opt.Contains("TRACKS"); | |
166 | Bool_t clearClusters = opt.Contains("CLUSTERS"); | |
167 | if (opt.Contains("ALL")) { | |
168 | clearTracks = kTRUE; | |
169 | clearClusters = kTRUE; | |
170 | } | |
171 | ||
172 | //fTracks is a collection of esdTrack | |
173 | //When clearing this array, the linked objects should not be deleted | |
174 | if (fTracks != 0x0 && clearTracks) { | |
175 | fTracks->Clear(); | |
176 | delete fTracks; | |
177 | fTracks = 0; | |
178 | } | |
179 | if (fClusters != 0x0 && clearClusters) { | |
180 | fClusters->Delete(); | |
181 | delete fClusters; | |
182 | fClusters = 0; | |
183 | } | |
184 | } | |
185 | // | |
186 | //------------------------------------------------------------------------------ | |
187 | // | |
188 | Int_t AliEMCALTracker::LoadClusters(TTree *cTree) | |
189 | { | |
190 | // | |
191 | // Load EMCAL clusters in the form of AliEMCALRecPoint, | |
192 | // from simulation temporary files. | |
193 | // (When included in reconstruction chain, this method is used automatically) | |
194 | // | |
195 | ||
196 | Clear("CLUSTERS"); | |
197 | ||
198 | cTree->SetBranchStatus("*",0); //disable all branches | |
199 | cTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need | |
200 | ||
201 | TBranch *branch = cTree->GetBranch("EMCALECARP"); | |
202 | if (!branch) { | |
203 | AliError("Can't get the branch with the EMCAL clusters"); | |
204 | return 1; | |
205 | } | |
206 | ||
207 | TClonesArray *clusters = new TClonesArray("AliEMCALRecPoint", 1000); | |
208 | branch->SetAddress(&clusters); | |
209 | ||
210 | //cTree->GetEvent(0); | |
211 | branch->GetEntry(0); | |
212 | Int_t nClusters = (Int_t)clusters->GetEntries(); | |
213 | if(fClusters) fClusters->Delete(); | |
214 | else fClusters = new TObjArray(0); | |
215 | for (Int_t i = 0; i < nClusters; i++) { | |
216 | AliEMCALRecPoint *cluster = (AliEMCALRecPoint*)clusters->At(i); | |
217 | if (!cluster) continue; | |
218 | AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster); | |
219 | fClusters->AddLast(matchCluster); | |
220 | } | |
221 | ||
222 | branch->SetAddress(0); | |
223 | clusters->Delete(); | |
224 | delete clusters; | |
225 | ||
226 | AliInfo(Form("Collected %d RecPoints from Tree", fClusters->GetEntries())); | |
227 | ||
228 | return 0; | |
229 | } | |
230 | // | |
231 | //------------------------------------------------------------------------------ | |
232 | // | |
233 | Int_t AliEMCALTracker::LoadClusters(AliESDEvent *esd) | |
234 | { | |
235 | // | |
236 | // Load EMCAL clusters in the form of AliESDCaloClusters, | |
237 | // from an AliESD object. | |
238 | // | |
239 | ||
240 | // make sure that tracks/clusters collections are empty | |
241 | Clear("CLUSTERS"); | |
242 | fClusters = new TObjArray(0); | |
243 | ||
244 | Int_t nClusters = esd->GetNumberOfCaloClusters(); | |
245 | for (Int_t i=0; i<nClusters; i++) | |
246 | { | |
247 | AliESDCaloCluster *cluster = esd->GetCaloCluster(i); | |
248 | if (!cluster || !cluster->IsEMCAL()) continue ; | |
249 | AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster); | |
250 | fClusters->AddLast(matchCluster); | |
251 | } | |
252 | ||
253 | AliInfo(Form("Collected %d clusters from ESD", fClusters->GetEntries())); | |
254 | return 0; | |
255 | } | |
256 | // | |
257 | //------------------------------------------------------------------------------ | |
258 | // | |
259 | Int_t AliEMCALTracker::LoadTracks(AliESDEvent *esd) | |
260 | { | |
261 | // | |
262 | // Load ESD tracks. | |
263 | // | |
264 | ||
265 | UInt_t mask1 = esd->GetESDRun()->GetDetectorsInDAQ(); | |
266 | UInt_t mask2 = esd->GetESDRun()->GetDetectorsInReco(); | |
267 | Bool_t desc1 = (mask1 >> 3) & 0x1; | |
268 | Bool_t desc2 = (mask2 >> 3) & 0x1; | |
269 | if (desc1==0 || desc2==0) { | |
270 | // AliError(Form("TPC not in DAQ/RECO: %u (%u)/%u (%u)", | |
271 | // mask1, esd->GetESDRun()->GetDetectorsInReco(), | |
272 | // mask2, esd->GetESDRun()->GetDetectorsInDAQ())); | |
273 | fITSTrackSA = kTRUE; | |
274 | } | |
275 | ||
276 | Clear("TRACKS"); | |
277 | fTracks = new TObjArray(0); | |
278 | ||
279 | Int_t nTracks = esd->GetNumberOfTracks(); | |
280 | //Bool_t isKink=kFALSE; | |
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); | |
286 | esdTrack->ResetStatus(AliESDtrack::kEMCALmatch); | |
287 | ||
288 | //Select good quaulity tracks | |
289 | if(esdTrack->Pt()<fCutPt) continue; | |
290 | if(!fITSTrackSA) | |
291 | if(esdTrack->GetNcls(1)<fCutNTPC)continue; | |
292 | ||
293 | //Loose geometric cut | |
294 | Double_t phi = esdTrack->Phi()*TMath::RadToDeg(); | |
295 | if(TMath::Abs(esdTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue; | |
296 | ||
297 | fTracks->AddLast(esdTrack); | |
298 | } | |
299 | ||
300 | AliInfo(Form("Collected %d tracks", fTracks->GetEntries())); | |
301 | return 0; | |
302 | } | |
303 | // | |
304 | //------------------------------------------------------------------------------ | |
305 | // | |
306 | void 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 | } | |
329 | } | |
330 | // | |
331 | //------------------------------------------------------------------------------ | |
332 | // | |
333 | Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd) | |
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 | // | |
341 | // | |
342 | // Note: should always return 0=OK, because otherwise all tracking | |
343 | // is aborted for this event | |
344 | ||
345 | if (!esd) { | |
346 | AliError("NULL ESD passed"); | |
347 | return 1; | |
348 | } | |
349 | ||
350 | // step 1: collect clusters | |
351 | Int_t okLoadClusters, nClusters; | |
352 | if (!fClusters || (fClusters && fClusters->IsEmpty())) { | |
353 | okLoadClusters = LoadClusters(esd); | |
354 | } | |
355 | nClusters = fClusters->GetEntries(); | |
356 | ||
357 | // step 2: collect ESD tracks | |
358 | Int_t nTracks, okLoadTracks; | |
359 | okLoadTracks = LoadTracks(esd); | |
360 | nTracks = fTracks->GetEntries(); | |
361 | ||
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()); | |
372 | track->SetStatus(AliESDtrack::kEMCALmatch); | |
373 | } | |
374 | } | |
375 | ||
376 | return 0; | |
377 | } | |
378 | ||
379 | // | |
380 | //------------------------------------------------------------------------------ | |
381 | // | |
382 | Int_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 | // | |
388 | ||
389 | Float_t maxEta=fCutEta; | |
390 | Float_t maxPhi=fCutPhi; | |
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 | |
395 | AliExternalTrackParam *trkParam = 0; | |
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 | } | |
404 | else | |
405 | trkParam = new AliExternalTrackParam(*track); | |
406 | ||
407 | if(!trkParam) return index; | |
408 | ||
409 | AliExternalTrackParam trkParamTmp(*trkParam); | |
410 | Float_t eta, phi, pt; | |
411 | if(!AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&trkParamTmp, fEMCalSurfaceDistance, track->GetMass(kTRUE), fStep, eta, phi, pt)) { | |
412 | if(fITSTrackSA) delete trkParam; | |
413 | return index; | |
414 | } | |
415 | track->SetTrackPhiEtaPtOnEMCal(phi,eta,pt); | |
416 | if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()){ | |
417 | if(fITSTrackSA) delete trkParam; | |
418 | return index; | |
419 | } | |
420 | ||
421 | //Perform extrapolation | |
422 | Double_t trkPos[3]; | |
423 | trkParamTmp.GetXYZ(trkPos); | |
424 | Int_t nclusters = fClusters->GetEntries(); | |
425 | for(Int_t ic=0; ic<nclusters; ic++) | |
426 | { | |
427 | AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic); | |
428 | Float_t clsPos[3] = {cluster->X(),cluster->Y(),cluster->Z()}; | |
429 | 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)); | |
430 | // printf("\n dR=%f,wind=%f\n",dR,fClusterWindow); //MARCEL | |
431 | if(dR > fClusterWindow) continue; | |
432 | ||
433 | AliExternalTrackParam trkParTmp(trkParamTmp); | |
434 | ||
435 | Float_t tmpEta, tmpPhi; | |
436 | if(!AliEMCALRecoUtils::ExtrapolateTrackToPosition(&trkParTmp, clsPos,track->GetMass(kTRUE), 5, tmpEta, tmpPhi)) continue; | |
437 | if(TMath::Abs(tmpPhi)<TMath::Abs(maxPhi) && TMath::Abs(tmpEta)<TMath::Abs(maxEta)) | |
438 | { | |
439 | maxPhi=tmpPhi; | |
440 | maxEta=tmpEta; | |
441 | index=ic; | |
442 | } | |
443 | } | |
444 | ||
445 | if(fITSTrackSA) delete trkParam; | |
446 | return index; | |
447 | } | |
448 | ||
449 | // | |
450 | //------------------------------------------------------------------------------ | |
451 | // | |
452 | void AliEMCALTracker::UnloadClusters() | |
453 | { | |
454 | // | |
455 | // Free memory from all arrays | |
456 | // This method is called after the local tracking step | |
457 | // so we can safely delete everything | |
458 | // | |
459 | ||
460 | Clear(); | |
461 | } | |
462 | ||
463 | // | |
464 | //------------------------------------------------------------------------------ | |
465 | // | |
466 | AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliEMCALRecPoint *recPoint) | |
467 | : fIndex(index), | |
468 | fX(0.), | |
469 | fY(0.), | |
470 | fZ(0.) | |
471 | { | |
472 | // | |
473 | // Translates an AliEMCALRecPoint object into the internal format. | |
474 | // Index of passed cluster in its native array must be specified. | |
475 | // | |
476 | TVector3 clpos; | |
477 | recPoint->GetGlobalPosition(clpos); | |
478 | ||
479 | fX = clpos.X(); | |
480 | fY = clpos.Y(); | |
481 | fZ = clpos.Z(); | |
482 | } | |
483 | // | |
484 | //------------------------------------------------------------------------------ | |
485 | // | |
486 | AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliESDCaloCluster *caloCluster) | |
487 | : fIndex(index), | |
488 | fX(0.), | |
489 | fY(0.), | |
490 | fZ(0.) | |
491 | { | |
492 | // | |
493 | // Translates an AliESDCaloCluster object into the internal format. | |
494 | // Index of passed cluster in its native array must be specified. | |
495 | // | |
496 | Float_t clpos[3]= {0., 0., 0.}; | |
497 | caloCluster->GetPosition(clpos); | |
498 | ||
499 | fX = (Double_t)clpos[0]; | |
500 | fY = (Double_t)clpos[1]; | |
501 | fZ = (Double_t)clpos[2]; | |
502 | } |