]>
Commit | Line | Data |
---|---|---|
4c904d93 | 1 | /************************************************************************** |
2 | * This file is property of and copyright by the ALICE HLT Project * | |
3 | * All rights reserved. * | |
4 | * * | |
5 | * Primary Authors: * | |
6 | * Artur Szostak <artursz@iafrica.com> * | |
7 | * * | |
8 | * Permission to use, copy, modify and distribute this software and its * | |
9 | * documentation strictly for non-commercial purposes is hereby granted * | |
10 | * without fee, provided that the above copyright notice appears in all * | |
11 | * copies and that both the copyright notice and this permission notice * | |
12 | * appear in the supporting documentation. The authors make no claims * | |
13 | * about the suitability of this software for any purpose. It is * | |
14 | * provided "as is" without express or implied warranty. * | |
15 | **************************************************************************/ | |
16 | ||
17 | /// | |
18 | /// @file AliAnalysisTaskLinkToMC.cxx | |
19 | /// @author Artur Szostak <artursz@iafrica.com> | |
20 | /// @date 27 Oct 2008 | |
21 | /// @brief Implementation of the AliAnalysisTaskLinkToMC task. | |
22 | /// | |
23 | /// The AliAnalysisTaskLinkToMC connects reconstructed tracks found in the ESD | |
24 | /// to their corresponding Monte Carlo (MC) tracks via the proximity of their | |
25 | /// hit coordinates. Each ESD track is copied over to an AOD track and put into | |
26 | /// the AOD event container, but the AOD track label is filled such that the | |
27 | /// AOD track label is the same as the MC track label. | |
28 | /// If a MC track could not be found for a given ESD track then -1 is filled | |
29 | /// into the AOD track label. | |
30 | /// The corresponding MC track is found by correlating the reconstructed track | |
31 | /// clusters with the MC track references. For each cluster the closest track | |
32 | /// reference coordinate is found and a chi-squared value is computed for each | |
33 | /// X and Y dimension independently. The chi-square's are summed up and the MC | |
34 | /// track with the smallest total chi-squared is chosen as the best match to | |
35 | /// the ESD track. | |
36 | /// | |
37 | /// The following cuts control the conditions under which we consider a MC track | |
38 | /// to match a given ESD track: | |
39 | /// 1) SigmaCut(Double_t value) - This is maximum distance expressed in standard | |
40 | /// deviations that is allowed between a reconstructed cluster coordinate | |
41 | /// and the MC track reference coordinate. The condition is expressed as: | |
42 | /// (x_{cluster} - x_{trackref})^2 / (sigma_{x})^2 < #sigma_{cut} | |
43 | /// This conditions is applied to each X and Y dimension respectively and | |
44 | /// cluster / track reference pairs that do not pass this cut are considered | |
45 | /// not to match. | |
46 | /// 2) MinClusters(Int_t value) - Indicates the minimum number of clusters that | |
47 | /// must be matched to a track reference from a MC track for the ESD and | |
48 | /// MC track to be considered to match. | |
49 | /// The following are absolute cuts which are used mainly as a speed optimisation | |
50 | /// so as to skip over cluster / track reference pairs that are far away from | |
51 | /// each other early. | |
52 | /// 3) HardCutLimitX(Double_t value) - The absolute maximum distance in centimetres | |
53 | /// allowed between a cluster's X coordinate and a track reference's X coordinate. | |
54 | /// If a cluster / track reference pair does not pass this condition then the | |
55 | /// pair is considered not to match. | |
56 | /// 4) HardCutLimitY(Double_t value) - Similar to HardCutLimitX but applied to the | |
57 | /// Y dimension of the cluster and track reference coordinate. | |
58 | /// 5) HardCutLimitZ(Double_t value) - Similar to HardCutLimitX but applied to the | |
59 | /// Z dimension. | |
60 | /// | |
61 | /// The analysis task also optionally generates histograms which show can be used | |
62 | /// to check the quality of the correlation between ESD and MC tracks. | |
63 | /// The histograms are generated in 2D for the pT and rapidity variable, | |
64 | /// and are called: | |
65 | /// 1) findableTracksHist - Filled with the MC tracks which should in principle | |
66 | /// be findable by the reconstruction software. A findable track is defined | |
67 | /// by the IsFindable method. | |
68 | /// 2) foundTracksHistMC - Filled with the MC information for each muon track that | |
69 | /// was found in the ESD and for which we could match its corresponding MC track. | |
70 | /// 3) foundTracksHist - Filled just like the foundTracksHistMC histogram, but with | |
71 | /// the momentum information taken from the ESD track and not the corresponding | |
72 | /// MC track. | |
73 | /// 4) fakeTracksHist - Filled with all the ESD tracks for which we could not find | |
74 | /// its matching MC track and assume it is a fake track. | |
75 | /// These histograms are filled by default, but can be disabled by using the | |
76 | /// GenerateHistograms(false) method call. | |
77 | /// | |
78 | ||
79 | #include "AliAnalysisTaskLinkToMC.h" | |
80 | #include "AliESDEvent.h" | |
81 | #include "AliESDMuonTrack.h" | |
82 | #include "AliESDMuonCluster.h" | |
83 | #include "AliMCEvent.h" | |
84 | #include "AliMCParticle.h" | |
85 | #include "AliAODVertex.h" | |
86 | #include "AliAODEvent.h" | |
87 | #include "AliAODTrack.h" | |
88 | #include "AliTrackReference.h" | |
89 | #include "AliStack.h" | |
90 | #include "AliLog.h" | |
91 | #include "TObjArray.h" | |
92 | #include "TCanvas.h" | |
93 | #include "TMath.h" | |
94 | #include "TMap.h" | |
95 | #include "TList.h" | |
96 | #include "TH2D.h" | |
97 | #include "TStyle.h" | |
98 | ||
805e576c | 99 | ClassImp(AliAnalysisTaskLinkToMC) |
4c904d93 | 100 | |
101 | ||
102 | void AliAnalysisTaskLinkToMC::HardCutLimitX(Double_t value) | |
103 | { | |
104 | /// Sets the absolute maximum distance in centimetres between the cluster | |
105 | /// and track reference that is allowed in the X direction. | |
106 | ||
107 | if (value >= 0) | |
108 | { | |
109 | fHardCutLimitX = value; | |
110 | } | |
111 | else | |
112 | { | |
113 | AliWarning("Setting value for fHardCutLimitX which is negative." | |
114 | " Will set it to a positive value." | |
115 | ); | |
116 | fHardCutLimitX = TMath::Abs(value); | |
117 | } | |
118 | } | |
119 | ||
120 | ||
121 | void AliAnalysisTaskLinkToMC::HardCutLimitY(Double_t value) | |
122 | { | |
123 | /// Sets the absolute maximum distance in centimetres between the cluster | |
124 | /// and track reference that is allowed in the Y direction. | |
125 | ||
126 | if (value >= 0) | |
127 | { | |
128 | fHardCutLimitY = value; | |
129 | } | |
130 | else | |
131 | { | |
132 | AliWarning("Setting value for fHardCutLimitY which is negative." | |
133 | " Will set it to a positive value." | |
134 | ); | |
135 | fHardCutLimitY = TMath::Abs(value); | |
136 | } | |
137 | } | |
138 | ||
139 | ||
140 | void AliAnalysisTaskLinkToMC::HardCutLimitZ(Double_t value) | |
141 | { | |
142 | /// Sets the absolute maximum distance in centimetres between the cluster | |
143 | /// and track reference that is allowed in the Z direction. | |
144 | ||
145 | if (value >= 0) | |
146 | { | |
147 | fHardCutLimitZ = value; | |
148 | } | |
149 | else | |
150 | { | |
151 | AliWarning("Setting value for fHardCutLimitZ which is negative." | |
152 | " Will set it to a positive value." | |
153 | ); | |
154 | fHardCutLimitZ = TMath::Abs(value); | |
155 | } | |
156 | } | |
157 | ||
158 | ||
159 | Double_t AliAnalysisTaskLinkToMC::SigmaCut() const | |
160 | { | |
161 | /// Returns the maximum sigma value to allow for each X and Y dimension | |
162 | /// between clusters and track references. | |
163 | ||
164 | return TMath::Sqrt(fSigmaCut2); | |
165 | } | |
166 | ||
167 | ||
168 | Bool_t AliAnalysisTaskLinkToMC::ChamberMustMatch(Int_t chamber) const | |
169 | { | |
170 | /// Returns the flag indicating if the ESD track cluster on the given | |
171 | /// chamber must match a corresponding track reference to be valid. | |
172 | /// \param chamber The chamber number must be in the range 1 to 14. | |
173 | ||
174 | if (chamber >= 1 and chamber <= 14) | |
175 | { | |
176 | return fChamberMustMatch[chamber-1]; | |
177 | } | |
178 | else | |
179 | { | |
180 | AliError(Form("The chamber number %d is invalid. Expected a value" | |
181 | " in the range [1..14]. Will return a value of 'false'.", | |
182 | chamber | |
183 | )); | |
184 | return false; | |
185 | } | |
186 | } | |
187 | ||
188 | ||
189 | void AliAnalysisTaskLinkToMC::ChamberMustMatch(Int_t chamber, Bool_t value) | |
190 | { | |
191 | /// Set the flag indicating if the ESD track cluster on the given chamber | |
192 | /// must match a corresponding track reference to be valid. | |
193 | /// \param chamber The chamber number must be in the range 1 to 14. | |
194 | /// \param value The new value to set to. | |
195 | ||
196 | if (chamber >= 1 and chamber <= 14) | |
197 | { | |
198 | fChamberMustMatch[chamber-1] = value; | |
199 | } | |
200 | else | |
201 | { | |
202 | AliError(Form("The chamber number %d is invalid. Expected a value" | |
203 | " in the range [1..14]. Will not set any fChamberMustMatch flag.", | |
204 | chamber | |
205 | )); | |
206 | } | |
207 | } | |
208 | ||
209 | ||
210 | Bool_t AliAnalysisTaskLinkToMC::StationMustMatch(Int_t station) const | |
211 | { | |
212 | /// Returns the flag indicating if at least one ESD track cluster on the | |
213 | /// given station must match a corresponding track reference to be valid. | |
214 | /// \param station The station number must be in the range 1..7. Station | |
215 | /// 1 to 5 are the tracking stations while 6 is MT1 and 7 is MT2. | |
216 | ||
217 | if (station >= 1 and station <= 7) | |
218 | { | |
219 | return fStationMustMatch[station-1]; | |
220 | } | |
221 | else | |
222 | { | |
223 | AliError(Form("The station number %d is invalid. Expected a value" | |
224 | " in the range [1..7]. Will return a value of 'false'.", | |
225 | station | |
226 | )); | |
227 | return false; | |
228 | } | |
229 | } | |
230 | ||
231 | ||
232 | void AliAnalysisTaskLinkToMC::StationMustMatch(Int_t station, Bool_t value) | |
233 | { | |
234 | /// Sets the flag indicating if at least one ESD track cluster on the | |
235 | /// given station must match a corresponding track reference to be valid. | |
236 | /// \param station The station number must be in the range 1..7. Station | |
237 | /// 1 to 5 are the tracking stations while 6 is MT1 and 7 is MT2. | |
238 | /// \param value The new value to set to. | |
239 | ||
240 | if (station >= 1 and station <= 7) | |
241 | { | |
242 | fStationMustMatch[station-1] = value; | |
243 | } | |
244 | else | |
245 | { | |
246 | AliError(Form("The station number %d is invalid. Expected a value" | |
247 | " in the range [1..7]. Will not set any fStationMustMatch flag.", | |
248 | station | |
249 | )); | |
250 | } | |
251 | } | |
252 | ||
253 | ||
254 | AliAnalysisTaskLinkToMC::AliAnalysisTaskLinkToMC() : | |
255 | AliAnalysisTaskSE(), | |
256 | fHistos(NULL), | |
257 | fFindableHist(NULL), | |
258 | fFoundHistMC(NULL), | |
259 | fFoundHist(NULL), | |
260 | fFakeHist(NULL), | |
261 | fHardCutLimitX(4.), // cm | |
262 | fHardCutLimitY(4.), // cm | |
263 | fHardCutLimitZ(4.), // cm | |
264 | fVarX(0.144*0.144), // cm^2 | |
265 | fVarY(0.01*0.01), // cm^2 | |
266 | fSigmaCut2(5*5), // sigma^2 | |
267 | fMinClusters(6), | |
268 | fMinClustersInSt12(0), // no minimum requirement for stations 1 and 2. | |
269 | fMinClustersInSt45(3), // 3 hits required on station 4 and 5 by default. | |
270 | fMakeHists(true), // Generate histograms for monitoring by default | |
271 | fUseZCoordinate(false) // Just use detElemId to match in the Z direction. | |
272 | { | |
273 | /// Default constructor. | |
274 | /// Sets input slot 0 to TChain input for ESDs, output slot 0 to TTree | |
275 | /// for AODs and output slot 1 to TList for the list of output histograms. | |
276 | /// The output histograms generated are to cross check the quality of the | |
277 | /// correlation and their generation can be turned off with | |
278 | /// AliAnalysisTaskLinkToMC::GenerateHistograms(false). | |
279 | ||
280 | // Do not force any matching for particular chambers by default. | |
281 | for (Int_t i = 0; i < 14; i++) | |
282 | { | |
283 | fChamberMustMatch[i] = false; | |
284 | } | |
285 | // Only force a minimum of 1 hit to be matched on stations 1, 2 and 3 by default. | |
286 | for (Int_t i = 0; i < 3; i++) | |
287 | { | |
288 | fStationMustMatch[i] = true; | |
289 | } | |
290 | for (Int_t i = 3; i < 7; i++) | |
291 | { | |
292 | fStationMustMatch[i] = false; | |
293 | } | |
294 | ||
295 | DefineOutput(1, TList::Class()); | |
296 | } | |
297 | ||
298 | ||
299 | AliAnalysisTaskLinkToMC::AliAnalysisTaskLinkToMC(const char* name) : | |
300 | AliAnalysisTaskSE(name), | |
301 | fHistos(NULL), | |
302 | fFindableHist(NULL), | |
303 | fFoundHistMC(NULL), | |
304 | fFoundHist(NULL), | |
305 | fFakeHist(NULL), | |
306 | fHardCutLimitX(4.), // cm | |
307 | fHardCutLimitY(4.), // cm | |
308 | fHardCutLimitZ(4.), // cm | |
309 | fVarX(0.144*0.144), // cm^2 | |
310 | fVarY(0.01*0.01), // cm^2 | |
311 | fSigmaCut2(5*5), // sigma^2 | |
312 | fMinClusters(6), | |
313 | fMinClustersInSt12(0), // no minimum requirement for stations 1 and 2. | |
314 | fMinClustersInSt45(3), // 3 hits required on station 4 and 5 by default. | |
315 | fMakeHists(true), // Generate histograms for monitoring by default | |
316 | fUseZCoordinate(false) // Just use detElemId to match in the Z direction. | |
317 | { | |
318 | /// Constructor with a task name. | |
319 | /// Sets input slot 0 to TChain input for ESDs, output slot 0 to TTree | |
320 | /// for AODs and output slot 1 to TList for the list of output histograms | |
321 | /// storing histograms generated to cross check the quility of the | |
322 | /// correlation. | |
323 | /// \param name The name of the task. | |
324 | ||
325 | // Do not force any matching for particular chambers by default. | |
326 | for (Int_t i = 0; i < 14; i++) | |
327 | { | |
328 | fChamberMustMatch[i] = false; | |
329 | } | |
330 | // Only force a minimum of 1 hit to be matched on stations 1, 2 and 3 by default. | |
331 | for (Int_t i = 0; i < 3; i++) | |
332 | { | |
333 | fStationMustMatch[i] = true; | |
334 | } | |
335 | for (Int_t i = 3; i < 7; i++) | |
336 | { | |
337 | fStationMustMatch[i] = false; | |
338 | } | |
339 | ||
340 | DefineOutput(1, TList::Class()); | |
341 | } | |
342 | ||
343 | ||
344 | void AliAnalysisTaskLinkToMC::UserCreateOutputObjects() | |
345 | { | |
346 | /// Creates the output histograms containing findable tracks, found tracks | |
347 | /// and fake tracks. The binning range for all histograms is 30 bins along | |
348 | /// the pT dimension (X direction) from 0 to 20 GeV/c and 30 bins in the | |
349 | /// rapidity dimension (Y direction) form -4 to -2.4 units. | |
350 | ||
351 | fHistos = new TList; | |
352 | ||
353 | char name[1024]; | |
354 | char title[1024]; | |
355 | ||
356 | Int_t nBinsX = 30; | |
357 | Double_t minX = 0.; | |
358 | Double_t maxX = 20.; | |
359 | Int_t nBinsY = 30; | |
360 | Double_t minY = -4.; | |
361 | Double_t maxY = -2.4; | |
362 | ||
363 | sprintf(name, "findableTracksHist"); | |
364 | sprintf(title, "Findable tracks"); | |
365 | fFindableHist = new TH2D(name, title, nBinsX, minX, maxX, nBinsY, minY, maxY); | |
366 | fFindableHist->SetXTitle("p_{T} [GeV/c]"); | |
367 | fFindableHist->SetYTitle("Rapidity (Y)"); | |
368 | fHistos->Add(fFindableHist); | |
369 | ||
370 | sprintf(name, "foundTracksHistMC"); | |
371 | sprintf(title, "Found tracks (filled with Monte Carlo kinematic information)"); | |
372 | fFoundHistMC = new TH2D(name, title, nBinsX, minX, maxX, nBinsY, minY, maxY); | |
373 | fFoundHistMC->SetXTitle("p_{T} [GeV/c]"); | |
374 | fFoundHistMC->SetYTitle("Rapidity (Y)"); | |
375 | fHistos->Add(fFoundHistMC); | |
376 | ||
377 | sprintf(name, "foundTracksHist"); | |
378 | sprintf(title, "Found tracks (filled with reconstructed kinematic information)"); | |
379 | fFoundHist = new TH2D(name, title, nBinsX, minX, maxX, nBinsY, minY, maxY); | |
380 | fFoundHist->SetXTitle("p_{T} [GeV/c]"); | |
381 | fFoundHist->SetYTitle("Rapidity (Y)"); | |
382 | fHistos->Add(fFoundHist); | |
383 | ||
384 | sprintf(name, "fakeTracksHist"); | |
385 | sprintf(title, "Fake tracks"); | |
386 | fFakeHist = new TH2D(name, title, nBinsX, minX, maxX, nBinsY, minY, maxY); | |
387 | fFakeHist->SetXTitle("p_{T} [GeV/c]"); | |
388 | fFakeHist->SetYTitle("Rapidity (Y)"); | |
389 | fHistos->Add(fFakeHist); | |
390 | } | |
391 | ||
392 | ||
393 | void AliAnalysisTaskLinkToMC::UserExec(Option_t* /*option*/) | |
394 | { | |
395 | /// This is the top level method that performs the work of the task. | |
396 | /// It will fill the histograms for quality checking with the Monte | |
397 | /// Carlo (MC) information if they are to be generated. | |
398 | /// The input event is then processed which must be of type AliESDEvent. | |
399 | /// For all the ESD muon tracks we try to find the corresponding MC | |
400 | /// track and build a TMap for the mapping. | |
401 | /// If the AOD output event is available then it is filled with AOD | |
402 | /// tracks which copy the information form the ESD track and the AOD | |
403 | /// track label is filled with the corresponding MC track label. | |
404 | /// If a corresponding MC track could not be found then the AOD track | |
405 | /// label is filled with -1. | |
406 | /// Finally the checking histograms are further filled with the ESD | |
407 | /// information based on the ESD to MC track map. | |
408 | ||
409 | if (fMakeHists and | |
410 | (fFindableHist == NULL or fFoundHistMC == NULL or fFoundHist == NULL or fFakeHist == NULL) | |
411 | ) | |
412 | { | |
413 | AliError("Histograms have not been created."); | |
414 | return; | |
415 | } | |
416 | ||
417 | if (fMakeHists and MCEvent() != NULL) FillHistsFromMC(); | |
418 | ||
419 | if (InputEvent() != NULL) | |
420 | { | |
421 | if (InputEvent()->IsA() == AliESDEvent::Class()) | |
422 | { | |
423 | if (MCEvent() != NULL) | |
424 | { | |
425 | TMap links; | |
426 | FindLinksToMC(links); | |
427 | if (AODEvent() != NULL) CreateAODTracks(links); | |
428 | if (fMakeHists) FillHistsFromLinks(links); | |
429 | } | |
430 | else | |
431 | { | |
432 | AliError("To process the ESD event we must have the Monte Carlo data also."); | |
433 | } | |
434 | } | |
435 | else | |
436 | { | |
437 | AliError(Form("The input event is of type \"%s\". Do not know how to handle it so it will be skipped.", | |
438 | InputEvent()->ClassName() | |
439 | )); | |
440 | } | |
441 | } | |
442 | ||
443 | if (fMakeHists) PostData(1, fHistos); | |
444 | } | |
445 | ||
446 | ||
447 | void AliAnalysisTaskLinkToMC::Terminate(Option_t* /*option*/) | |
448 | { | |
449 | /// The terminate method will draw the histograms filled by this task if | |
450 | /// they were filled by setting appropriate flag with the | |
451 | /// GenerateHistograms(true) method, which is the default. | |
452 | ||
453 | if (fMakeHists) | |
454 | { | |
455 | gStyle->SetPalette(1); | |
456 | new TCanvas; | |
457 | fFindableHist->DrawCopy("colz"); | |
458 | new TCanvas; | |
459 | fFoundHistMC->DrawCopy("colz"); | |
460 | new TCanvas; | |
461 | fFoundHist->DrawCopy("colz"); | |
462 | new TCanvas; | |
463 | fFakeHist->DrawCopy("colz"); | |
464 | ||
465 | AliInfo(Form("Findable tracks = %d", Int_t(fFindableHist->GetEntries()) )); | |
466 | AliInfo(Form("Found tracks = %d", Int_t(fFoundHistMC->GetEntries()) )); | |
467 | AliInfo(Form("Fake tracks = %d", Int_t(fFakeHist->GetEntries()) )); | |
468 | } | |
469 | } | |
470 | ||
471 | ||
472 | bool AliAnalysisTaskLinkToMC::IsFindable(AliMCParticle* track) const | |
473 | { | |
474 | /// This method is used to identify if a Monte Carlo (MC) track is in principle | |
475 | /// findable in the muon spectrometer. | |
476 | /// This means the MC track must be a muon particle and leave at least one track | |
477 | /// reference point on either chamber in each tracking station. | |
478 | /// \param track The MC track to check. | |
479 | /// \returns true if the track is findable by the muon spectrometer | |
480 | /// and false otherwise. | |
481 | /// \note The definition of track find-ability does not affect the correlation | |
482 | /// between ESD and MC tracks nor the generation of the AOD output tracks, | |
483 | /// only the histogram of findable tracks is effected if it is generated at all. | |
484 | ||
485 | // Select only muons. | |
486 | if (TMath::Abs(track->Particle()->GetPdgCode()) != 13) return false; | |
487 | ||
488 | // Build hit mask for chambers. | |
489 | bool hitOnCh[14]; | |
490 | for (Int_t i = 0; i < 14; i++) | |
491 | { | |
492 | hitOnCh[i] = false; | |
493 | } | |
494 | for (Int_t i = 0; i < track->GetNumberOfTrackReferences(); i++) | |
495 | { | |
496 | AliTrackReference* ref = track->GetTrackReference(i); | |
497 | if (ref->DetectorId() != AliTrackReference::kMUON) continue; | |
498 | Int_t chamberId = ref->UserId() / 100 - 1; | |
499 | if (0 <= chamberId and chamberId < 14) | |
500 | { | |
501 | hitOnCh[chamberId] = true; | |
502 | } | |
503 | } | |
504 | ||
505 | // Check that enough hits are available on all tracking chambers. | |
506 | bool hitOnSt1 = hitOnCh[0] or hitOnCh[1]; | |
507 | bool hitOnSt2 = hitOnCh[2] or hitOnCh[3]; | |
508 | bool hitOnSt3 = hitOnCh[4] or hitOnCh[5]; | |
509 | bool hitOnSt4 = hitOnCh[6] or hitOnCh[7]; | |
510 | bool hitOnSt5 = hitOnCh[8] or hitOnCh[9]; | |
511 | ||
512 | return (hitOnSt1 and hitOnSt2 and hitOnSt3 and hitOnSt4 and hitOnSt5); | |
513 | } | |
514 | ||
515 | ||
516 | void AliAnalysisTaskLinkToMC::FillHistsFromMC() | |
517 | { | |
518 | /// Fills the histogram containing findable tracks from the MC event information. | |
519 | /// Only MC tracks for which the IsFindable method returns true are added | |
520 | /// to the histogram. | |
521 | ||
522 | for (Int_t i = 0; i < MCEvent()->GetNumberOfTracks(); i++) | |
523 | { | |
7aad0c47 | 524 | AliMCParticle* mcTrack = (AliMCParticle*) MCEvent()->GetTrack(i); |
4c904d93 | 525 | |
526 | // Select only reconstructible tracks to fill into findable hist. | |
527 | if (IsFindable(mcTrack)) | |
528 | { | |
529 | fFindableHist->Fill(mcTrack->Pt(), mcTrack->Y()); | |
530 | } | |
531 | } | |
532 | } | |
533 | ||
534 | ||
535 | void AliAnalysisTaskLinkToMC::FillHistsFromLinks(TMap& links) | |
536 | { | |
537 | /// Fills the histograms containing found tracks and fake tracks. | |
538 | /// Found tracks are all those that have a ESD muon track in the ESD event | |
539 | /// and for which a corresponding MC track could be found. | |
540 | /// Fake tracks are ESD muon tracks for which no corresponding MC track | |
541 | /// could be matched so it is considered a fake track. | |
542 | /// \param links This contains the mapping between ESD muon tracks and | |
543 | /// MC tracks. ESD tracks for which no MC track could be matched should | |
544 | /// have the "value" of the key/value pair in the mapping set to NULL. | |
545 | ||
546 | TIter itmap(links.GetTable()); | |
547 | TPair* pair = NULL; | |
548 | while ( (pair = static_cast<TPair*>(itmap())) != NULL ) | |
549 | { | |
550 | if (pair->Value() != NULL) | |
551 | { | |
552 | AliESDMuonTrack* esdTrack = static_cast<AliESDMuonTrack*>(pair->Key()); | |
553 | fFoundHist->Fill(esdTrack->Pt(), esdTrack->Y()); | |
554 | ||
555 | AliMCParticle* mcTrack = static_cast<AliMCParticle*>(pair->Value()); | |
556 | fFoundHistMC->Fill(mcTrack->Pt(), mcTrack->Y()); | |
557 | } | |
558 | else | |
559 | { | |
560 | AliESDMuonTrack* esdTrack = static_cast<AliESDMuonTrack*>(pair->Key()); | |
561 | fFakeHist->Fill(esdTrack->Pt(), esdTrack->Y()); | |
562 | } | |
563 | } | |
564 | } | |
565 | ||
566 | ||
567 | void AliAnalysisTaskLinkToMC::CreateAODTracks(TMap& links) | |
568 | { | |
569 | /// This code is copied from AliAnalysisTaskESDMuonFilter::ConvertESDtoAOD with | |
570 | /// modifications to fill the MC track label using the ESD to MC track mapping. | |
571 | /// \param links This is the mapping between ESD tracks and MC tracks. | |
572 | ||
573 | // ESD Muon Filter analysis task executed for each event | |
574 | AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()); | |
575 | ||
576 | // Define arrays for muons | |
577 | Double_t pos[3]; | |
578 | Double_t p[3]; | |
579 | Double_t pid[10]; | |
580 | ||
581 | // has to be changed once the muon pid is provided by the ESD | |
805be8a6 | 582 | for (Int_t i = 0; i < 10; pid[i++] = 0.) {} |
4c904d93 | 583 | pid[AliAODTrack::kMuon]=1.; |
584 | ||
585 | AliAODHeader* header = AODEvent()->GetHeader(); | |
586 | AliAODTrack *aodTrack = 0x0; | |
587 | AliESDMuonTrack *esdMuTrack = 0x0; | |
588 | ||
589 | // Access to the AOD container of tracks | |
590 | TClonesArray &tracks = *(AODEvent()->GetTracks()); | |
591 | Int_t jTracks = tracks.GetEntriesFast(); | |
592 | ||
593 | // Read primary vertex from AOD event | |
594 | AliAODVertex *primary = AODEvent()->GetPrimaryVertex(); | |
595 | if (primary != NULL) primary->Print(); | |
596 | ||
597 | // Loop on muon tracks to fill the AOD track branch | |
598 | Int_t nMuTracks = esd->GetNumberOfMuonTracks(); | |
599 | printf("Number of Muon Tracks=%d\n",nMuTracks); | |
600 | ||
601 | // Update number of positive and negative tracks from AOD event (M.G.) | |
602 | Int_t nPosTracks = header->GetRefMultiplicityPos(); | |
603 | Int_t nNegTracks = header->GetRefMultiplicityNeg(); | |
604 | ||
605 | for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) | |
606 | { | |
607 | esdMuTrack = esd->GetMuonTrack(nMuTrack); | |
608 | ||
609 | //if (!esdMuTrack->ContainTrackerData()) continue; | |
610 | ||
611 | UInt_t selectInfo = 0; | |
612 | ||
613 | p[0] = esdMuTrack->Px(); | |
614 | p[1] = esdMuTrack->Py(); | |
615 | p[2] = esdMuTrack->Pz(); | |
616 | ||
617 | pos[0] = esdMuTrack->GetNonBendingCoor(); | |
618 | pos[1] = esdMuTrack->GetBendingCoor(); | |
619 | pos[2] = esdMuTrack->GetZ(); | |
620 | ||
621 | Int_t label = -1; | |
622 | TPair* pair = static_cast<TPair*>( links.FindObject(esdMuTrack) ); | |
623 | if (pair != NULL and pair->Value() != NULL) | |
624 | { | |
625 | AliMCParticle* mcTrack = static_cast<AliMCParticle*>(pair->Value()); | |
626 | label = mcTrack->GetLabel(); | |
627 | if (label == -1) | |
628 | { | |
629 | AliWarning("The MC track was set to -1."); | |
630 | } | |
631 | } | |
632 | ||
633 | aodTrack = new(tracks[jTracks++]) AliAODTrack( | |
634 | esdMuTrack->GetUniqueID(), // ID | |
635 | label, // label | |
636 | p, // momentum | |
637 | kTRUE, // cartesian coordinate system | |
638 | pos, // position | |
639 | kFALSE, // isDCA | |
640 | 0x0, // covariance matrix | |
641 | esdMuTrack->Charge(), // charge | |
642 | 0, // ITSClusterMap | |
643 | pid, // pid | |
644 | primary, // primary vertex | |
645 | kFALSE, // used for vertex fit? | |
646 | kFALSE, // used for primary vertex fit? | |
647 | AliAODTrack::kPrimary,// track type | |
648 | selectInfo | |
649 | ); | |
650 | ||
651 | aodTrack->SetXYAtDCA(esdMuTrack->GetNonBendingCoorAtDCA(), esdMuTrack->GetBendingCoorAtDCA()); | |
652 | aodTrack->SetPxPyPzAtDCA(esdMuTrack->PxAtDCA(), esdMuTrack->PyAtDCA(), esdMuTrack->PzAtDCA()); | |
653 | aodTrack->ConvertAliPIDtoAODPID(); | |
654 | aodTrack->SetChi2perNDF(esdMuTrack->GetChi2() / (2.*esdMuTrack->GetNHit() - 5.)); | |
655 | aodTrack->SetChi2MatchTrigger(esdMuTrack->GetChi2MatchTrigger()); | |
656 | aodTrack->SetHitsPatternInTrigCh(esdMuTrack->GetHitsPatternInTrigCh()); | |
657 | aodTrack->SetMuonClusterMap(esdMuTrack->GetMuonClusterMap()); | |
658 | aodTrack->SetMatchTrigger(esdMuTrack->GetMatchTrigger()); | |
659 | ||
660 | if (primary != NULL) primary->AddDaughter(aodTrack); | |
661 | ||
662 | if (esdMuTrack->Charge() > 0) nPosTracks++; | |
663 | else nNegTracks++; | |
664 | } | |
665 | ||
666 | header->SetRefMultiplicity(jTracks); | |
667 | header->SetRefMultiplicityPos(nPosTracks); | |
668 | header->SetRefMultiplicityNeg(nNegTracks); | |
669 | ||
670 | // Note: we do not have to call PostData for the AOD because the AliAnalysisTaskSE | |
671 | // base class already does so just after the UserExec. | |
672 | } | |
673 | ||
674 | ||
675 | void AliAnalysisTaskLinkToMC::FindLinksToMC(TMap& links) | |
676 | { | |
677 | /// Finds all MC tracks that correspond to their ESD tracks and creates a | |
678 | /// mapping between the ESD tracks and MC tracks which is filled into the | |
679 | /// 'links' TMap object. | |
680 | /// \param links This will be filled with the ESD and MC track pairs which | |
681 | /// have been matched. Thus, for each ESD track which is set as the TMap key, | |
682 | /// the corresponding MC track is identified by the TMap value. If no | |
683 | /// MC track could be matched then the MC track value is set to NULL. | |
684 | ||
685 | // Algorithm description for connecting reconstructed ESD muon tracks to MC: | |
686 | // | |
687 | // Build ESD track list | |
688 | // Build MC track list | |
689 | // while ESD track list not empty and MC track list not empty: | |
690 | // Find ESD/MC track pair with best compatibility. | |
691 | // if no pair is compatible then break loop. | |
692 | // Store found pair in return list. | |
693 | // Remove pair from ESD and MC track list. | |
694 | // Mark remaining tracks in ESD track list as fakes. | |
695 | // | |
696 | // To find ESD/MC track pair with best compatibility: | |
697 | // bestQuality := 1e100 | |
698 | // bestESDTrack := nil | |
699 | // bestMCTrack := nil | |
700 | // for each ESD track: | |
701 | // for each MC track: | |
702 | // chi^2 / numOfClustersMatched := calculate compatibility | |
703 | // if tracks are compatible and chi^2 / numOfClustersMatched < bestQuality then: | |
704 | // bestQuality := chi^2 / noOfClustersMatched | |
705 | // bestESDTrack := current ESD track | |
706 | // bestMCTrack := current MC track | |
707 | // | |
708 | // To calculate compatibility (ESD track, MC track): | |
709 | // chi^2 := 1e100 | |
710 | // numberOfClustersAssociated := 0 | |
711 | // for each cluster in ESD track: | |
712 | // find nearest 2 track refs from MC track | |
713 | // if 2 track refs found then find average of them into MChit | |
714 | // else MChit := nearest ref track. | |
715 | // k := calculate chi^2 distance between ESD cluster and MChit | |
716 | // if k <= maxValue then: | |
717 | // chi^2 := chi^2 + k | |
718 | // numberOfClustersAssociated := numberOfClustersAssociated + 1 | |
719 | // else: | |
720 | // chi^2 := chi^2 + maxAllowedChi2 | |
721 | // if numberOfClustersAssociated >= minGoodClusters then: | |
722 | // return (compatible, chi^2) | |
723 | // else: | |
724 | // return (not compatible, 1e100) | |
725 | ||
726 | AliESDEvent* esd = static_cast<AliESDEvent*>(InputEvent()); | |
727 | ||
728 | TObjArray esdTracks; | |
729 | TObjArray mcTracks; | |
730 | for (Int_t i = 0; i < esd->GetNumberOfMuonTracks(); i++) | |
731 | { | |
732 | esdTracks.Add(esd->GetMuonTrack(i)); | |
733 | } | |
734 | for (Int_t i = 0; i < MCEvent()->GetNumberOfTracks(); i++) | |
735 | { | |
736 | mcTracks.Add(MCEvent()->GetTrack(i)); | |
737 | } | |
738 | ||
739 | while (not esdTracks.IsEmpty() and not mcTracks.IsEmpty()) | |
740 | { | |
741 | Double_t bestQuality = 1e100; | |
742 | Int_t bestNoOfClusters = 0; | |
743 | AliESDMuonTrack* bestEsdTrack = NULL; | |
744 | AliMCParticle* bestMcTrack = NULL; | |
745 | TIter itesd(&esdTracks); | |
746 | TIter itmc(&mcTracks); | |
747 | AliESDMuonTrack* esdTrack = NULL; | |
748 | AliMCParticle* mcTrack = NULL; | |
749 | ||
750 | // Find the ESD and MC track pair that matches the best. | |
751 | while ( (esdTrack = static_cast<AliESDMuonTrack*>(itesd())) != NULL ) | |
752 | { | |
753 | while ( (mcTrack = static_cast<AliMCParticle*>(itmc())) != NULL ) | |
754 | { | |
755 | Double_t chi2 = 1e100; | |
756 | bool tracksCompatible = false; | |
757 | Int_t noOfClusters = 0; | |
758 | CalculateTrackCompatibility(esdTrack, mcTrack, tracksCompatible, chi2, noOfClusters); | |
759 | Double_t quality = noOfClusters > 0 ? chi2 / Double_t(noOfClusters) : 1e100; | |
760 | if (tracksCompatible and quality < bestQuality and noOfClusters >= bestNoOfClusters) | |
761 | { | |
762 | bestQuality = quality; | |
763 | bestNoOfClusters = noOfClusters; | |
764 | bestEsdTrack = esdTrack; | |
765 | bestMcTrack = mcTrack; | |
766 | } | |
767 | } | |
768 | itmc.Reset(); | |
769 | } | |
770 | ||
771 | // If no track pair could be matched then we are done because we | |
772 | // will not be able to match anything more next time. | |
773 | if (bestMcTrack == NULL) break; | |
774 | ||
775 | links.Add(bestEsdTrack, bestMcTrack); | |
776 | ||
777 | esdTracks.Remove(bestEsdTrack); | |
778 | mcTracks.Remove(bestMcTrack); | |
779 | } | |
780 | ||
781 | // Add all remaining ESD tracks which could not be matched to any MC track | |
782 | // to the mapping but with NULL for the MC track. | |
783 | TIter itesd(&esdTracks); | |
784 | AliESDMuonTrack* esdTrack = NULL; | |
785 | while ( (esdTrack = static_cast<AliESDMuonTrack*>(itesd())) != NULL ) | |
786 | { | |
787 | links.Add(esdTrack, NULL); | |
788 | } | |
789 | } | |
790 | ||
791 | ||
792 | void AliAnalysisTaskLinkToMC::CalculateTrackCompatibility( | |
793 | AliESDMuonTrack* esdTrack, AliMCParticle* mcTrack, | |
794 | bool& tracksCompatible, Double_t& chi2, Int_t& noOfClusters | |
795 | ) | |
796 | { | |
797 | /// Calculates the compatibility between a ESD and MC track pair. | |
798 | /// For each cluster of the ESD track, two track references are found | |
799 | /// which are spatially closest to it. If the distance between the cluster | |
800 | /// and track reference is larger than the HardCutLimit*() parameters then | |
801 | /// the track reference is not used. The average coordinate is calculated | |
802 | /// from the two track references found and the chi-square between the cluster | |
803 | /// and average coordinate is calculated. If only one track reference is | |
804 | /// found then it is used for the calculation without taking an average. | |
805 | /// If no track references are found then the cluster is not matched and | |
806 | /// we continue to the next cluster. | |
807 | /// If the chi-squared value for the cluster / track reference pair is | |
808 | /// larger than that allowed by the SigmaCut() parameter then the cluster | |
809 | /// is also considered not to be matched. | |
810 | /// The tracks are compatible if the ESD track has a minimum number of | |
811 | /// clusters matched with MC track references from the MC track. | |
812 | /// [in] \param esdTrack The ESD track we are trying to match. | |
813 | /// [in] \param mcTrack The MC track we are trying to match the ESD track to. | |
814 | /// [out] \param tracksCompatible Filled with true if the given esdTrack | |
815 | /// and mcTrack are compatible. | |
816 | /// [out] \param chi2 Filled with the total chi-squared calculated for all | |
817 | /// matched clusters. | |
818 | /// [out] \param noOfClusters Filled with the number of clusters that were | |
819 | /// matched to at least one track reference. | |
820 | ||
821 | chi2 = 0; | |
822 | noOfClusters = 0; | |
823 | tracksCompatible = false; | |
824 | Int_t noOfClustersSt12 = 0; | |
825 | Int_t noOfClustersSt45 = 0; | |
826 | bool chamberMatched[14] = { | |
827 | false, false, false, false, false, false, false, | |
828 | false, false, false, false, false, false, false | |
829 | }; | |
830 | bool stationMatched[7] = { | |
831 | false, false, false, false, false, false, false | |
832 | }; | |
833 | ||
834 | for (Int_t i = 0; i < esdTrack->GetNClusters(); i++) | |
835 | { | |
836 | AliESDMuonCluster* cluster = static_cast<AliESDMuonCluster*>( esdTrack->GetClusters()[i] ); | |
837 | Double_t varX = cluster->GetErrX2(); | |
838 | Double_t varY = cluster->GetErrY2(); | |
839 | // If the variance is zero then use the default one or just 1 | |
840 | // if the default is also not set. | |
841 | if (varX == 0) varX = fVarX != 0 ? fVarX : 1.; | |
842 | if (varY == 0) varY = fVarY != 0 ? fVarY : 1.; | |
843 | ||
844 | Double_t chi2A = 1e100; | |
845 | AliTrackReference* refA = NULL; | |
846 | Double_t chi2B = 1e100; | |
847 | AliTrackReference* refB = NULL; | |
848 | ||
849 | for (Int_t j = 0; j < mcTrack->GetNumberOfTrackReferences(); j++) | |
850 | { | |
851 | AliTrackReference* ref = mcTrack->GetTrackReference(j); | |
852 | ||
853 | if (ref->DetectorId() != AliTrackReference::kMUON) continue; | |
854 | if (ref->UserId() != cluster->GetDetElemId()) continue; | |
855 | ||
856 | Double_t dx = ref->X() - cluster->GetX(); | |
857 | if (TMath::Abs(dx) > fHardCutLimitX) continue; | |
858 | Double_t dy = ref->Y() - cluster->GetY(); | |
859 | if (TMath::Abs(dy) > fHardCutLimitY) continue; | |
860 | if (fUseZCoordinate) | |
861 | { | |
862 | Double_t dz = ref->Z() - cluster->GetZ(); | |
863 | if (TMath::Abs(dz) > fHardCutLimitZ) continue; | |
864 | } | |
865 | ||
866 | Double_t chi2sum = dx*dx / varX + dy*dy / varY; | |
867 | if (chi2sum < chi2A) | |
868 | { | |
869 | // Copy track reference A to B and set new value for A | |
870 | chi2B = chi2A; | |
871 | refB = refA; | |
872 | chi2A = chi2sum; | |
873 | refA = ref; | |
874 | } | |
875 | else if (chi2sum < chi2B) | |
876 | { | |
877 | // Track reference A still has smaller chi2 so just replace B. | |
878 | chi2B = chi2sum; | |
879 | refB = ref; | |
880 | } | |
881 | } | |
882 | ||
883 | Double_t x = 0, y = 0; | |
884 | if (refA != NULL and refB != NULL) | |
885 | { | |
886 | // Average the track references if 2 were found. | |
887 | x = (refA->X() + refB->X()) * 0.5; | |
888 | y = (refA->Y() + refB->Y()) * 0.5; | |
889 | } | |
890 | else if (refA != NULL) | |
891 | { | |
892 | x = refA->X(); | |
893 | y = refA->Y(); | |
894 | } | |
895 | else | |
896 | { | |
897 | continue; // no track reference hit found. | |
898 | } | |
899 | ||
900 | Double_t dx = x - cluster->GetX(); | |
901 | Double_t dy = y - cluster->GetY(); | |
902 | ||
903 | Double_t chi2x = dx*dx / varX; | |
904 | if (chi2x > fSigmaCut2) continue; | |
905 | Double_t chi2y = dy*dy / varY; | |
906 | if (chi2y > fSigmaCut2) continue; | |
907 | ||
908 | chi2 += chi2x + chi2y; | |
909 | noOfClusters++; | |
910 | if (cluster->GetChamberId() >= 0 and cluster->GetChamberId() < 4) | |
911 | { | |
912 | noOfClustersSt12++; | |
913 | } | |
914 | if (cluster->GetChamberId() >= 6 and cluster->GetChamberId() < 10) | |
915 | { | |
916 | noOfClustersSt45++; | |
917 | } | |
918 | ||
919 | if (0 <= cluster->GetChamberId() and cluster->GetChamberId() < 14) | |
920 | { | |
921 | chamberMatched[cluster->GetChamberId()] = true; | |
922 | stationMatched[cluster->GetChamberId()/2] = true; | |
923 | } | |
924 | } | |
925 | ||
926 | // Need to check that all the chambers and stations that had to match we | |
927 | // actually found matching track references. | |
928 | ||
929 | bool forcedChambersMatched = true; | |
930 | for (Int_t i = 0; i < 14; i++) | |
931 | { | |
932 | if (fChamberMustMatch[i] and not chamberMatched[i]) | |
933 | { | |
934 | forcedChambersMatched = false; | |
935 | } | |
936 | } | |
937 | bool forcedStationsMatched = true; | |
938 | for (Int_t i = 0; i < 7; i++) | |
939 | { | |
940 | if (fStationMustMatch[i] and not stationMatched[i]) | |
941 | { | |
942 | forcedStationsMatched = false; | |
943 | } | |
944 | } | |
945 | ||
946 | if (noOfClusters >= fMinClusters and | |
947 | noOfClustersSt12 >= fMinClustersInSt12 and | |
948 | noOfClustersSt45 >= fMinClustersInSt45 and | |
949 | forcedChambersMatched and forcedStationsMatched | |
950 | ) | |
951 | { | |
952 | tracksCompatible = true; | |
953 | } | |
954 | } | |
955 |