]>
Commit | Line | Data |
---|---|---|
8858ee12 | 1 | #include "TChain.h" |
2 | #include "TROOT.h" | |
3 | #include "TFile.h" | |
4 | #include "TError.h" | |
5 | #include "TH1.h" | |
6 | #include "TH2.h" | |
7 | #include "TF1.h" | |
8 | #include "TArrayI.h" | |
9 | #include "TCanvas.h" | |
10 | #include "TVector3.h" | |
11 | #include "TPDGCode.h" | |
12 | #include "TParticle.h" | |
13 | #include "AliESDEvent.h" | |
14 | #include "AliESDv0.h" | |
15 | #include "AliESDcascade.h" | |
16 | #include "AliESDMuonTrack.h" | |
17 | #include "AliESDCaloCluster.h" | |
18 | #include "AliRun.h" | |
19 | #include "AliHeader.h" | |
20 | #include "AliMCEvent.h" | |
21 | #include "AliGenEventHeader.h" | |
22 | #include "AliPID.h" | |
23 | ||
24 | #include "AliAnalysisTaskCheckESD.h" | |
25 | ||
26 | ClassImp(AliAnalysisTaskCheckESD) | |
27 | ||
28 | AliAnalysisTaskCheckESD::AliAnalysisTaskCheckESD(): | |
29 | AliAnalysisTaskSE("AliAnalysisTaskCheckESD"), | |
30 | fListOfHistos(0), | |
31 | hGen(0), | |
32 | hRec(0), | |
33 | hResPtInv(0), | |
34 | hResPhi(0), | |
35 | hResTheta(0), | |
36 | hDEdxRight(0), | |
37 | hDEdxWrong(0), | |
38 | hResTOFRight(0), | |
39 | hResTOFWrong(0), | |
40 | hEPHOS(0), | |
41 | hEEMCAL(0), | |
42 | hPtMUON(0), | |
43 | hMassK0(0), | |
44 | hMassLambda(0), | |
45 | hMassLambdaBar(0), | |
46 | hMassXi(0), | |
47 | hMassOmega(0), | |
48 | hScalars(0), | |
49 | hArrayHist(0) | |
50 | { | |
51 | // Default constructor | |
52 | // Define input and output slots here | |
53 | // Input slot #0 works with a TChain | |
54 | DefineInput(0, TChain::Class()); | |
55 | // Output slot #1 TList | |
56 | DefineOutput(1, TList::Class()); | |
57 | } | |
58 | ||
59 | AliAnalysisTaskCheckESD::AliAnalysisTaskCheckESD(const char* name): | |
60 | AliAnalysisTaskSE(name), | |
61 | fListOfHistos(0), | |
62 | hGen(0), | |
63 | hRec(0), | |
64 | hResPtInv(0), | |
65 | hResPhi(0), | |
66 | hResTheta(0), | |
67 | hDEdxRight(0), | |
68 | hDEdxWrong(0), | |
69 | hResTOFRight(0), | |
70 | hResTOFWrong(0), | |
71 | hEPHOS(0), | |
72 | hEEMCAL(0), | |
73 | hPtMUON(0), | |
74 | hMassK0(0), | |
75 | hMassLambda(0), | |
76 | hMassLambdaBar(0), | |
77 | hMassXi(0), | |
78 | hMassOmega(0), | |
79 | hScalars(0), | |
80 | hArrayHist(0) | |
81 | { | |
82 | // Constructor | |
83 | AliInfo("Constructor AliAnalysisTaskCheckESD"); | |
84 | // Define input and output slots here | |
85 | // Input slot #0 works with a TChain | |
86 | DefineInput(0, TChain::Class()); | |
87 | // Output slot #1 TList | |
88 | DefineOutput(1, TList::Class()); | |
89 | } | |
90 | ||
91 | TH1F * AliAnalysisTaskCheckESD::CreateHisto(const char* name, const char* title,Int_t nBins, | |
92 | Double_t xMin, Double_t xMax, | |
93 | const char* xLabel, const char* yLabel) | |
94 | { | |
95 | // create a histogram | |
96 | TH1F* result = new TH1F(name, title, nBins, xMin, xMax); | |
97 | result->SetOption("E"); | |
98 | if (xLabel) result->GetXaxis()->SetTitle(xLabel); | |
99 | if (yLabel) result->GetYaxis()->SetTitle(yLabel); | |
100 | result->SetMarkerStyle(kFullCircle); | |
101 | return result; | |
102 | } | |
103 | ||
104 | TH1F * AliAnalysisTaskCheckESD::CreateEffHisto(TH1F* hGen, TH1F* hRec) | |
105 | { | |
106 | // create an efficiency histogram | |
107 | Int_t nBins = hGen->GetNbinsX(); | |
108 | TH1F* hEff = (TH1F*) hGen->Clone("hEff"); | |
109 | hEff->SetTitle(""); | |
110 | hEff->SetStats(kFALSE); | |
111 | hEff->SetMinimum(0.); | |
112 | hEff->SetMaximum(110.); | |
113 | hEff->GetYaxis()->SetTitle("#epsilon [%]"); | |
114 | ||
115 | for (Int_t iBin = 0; iBin <= nBins; iBin++) { | |
116 | Double_t nGen_eff = hGen->GetBinContent(iBin); | |
117 | Double_t nRec_eff = hRec->GetBinContent(iBin); | |
118 | if (nGen_eff > 0) { | |
119 | Double_t eff = nRec_eff/nGen_eff; | |
120 | hEff->SetBinContent(iBin, 100. * eff); | |
121 | Double_t error = sqrt(eff*(1.-eff) / nGen_eff); | |
122 | if (error == 0) error = 0.0001; | |
123 | hEff->SetBinError(iBin, 100. * error); | |
124 | } | |
125 | else { | |
126 | hEff->SetBinContent(iBin, -100.); | |
127 | hEff->SetBinError(iBin, 0); | |
128 | } | |
129 | } | |
130 | return hEff; | |
131 | } | |
132 | ||
133 | ||
134 | Bool_t AliAnalysisTaskCheckESD::FitHisto(TH1* histo, Double_t& res, Double_t& resError) | |
135 | { | |
136 | // fit a gaussian to a histogram | |
137 | static TF1* fitFunc = new TF1("fitFunc", "gaus"); | |
138 | fitFunc->SetLineWidth(2); | |
139 | fitFunc->SetFillStyle(0); | |
140 | Double_t maxFitRange = 2; | |
141 | ||
142 | if (histo->Integral() > 50) { | |
143 | Float_t mean = histo->GetMean(); | |
144 | Float_t rms = histo->GetRMS(); | |
145 | fitFunc->SetRange(mean - maxFitRange*rms, mean + maxFitRange*rms); | |
146 | fitFunc->SetParameters(mean, rms); | |
147 | histo->Fit(fitFunc, "QRI0"); | |
148 | histo->GetFunction("fitFunc")->ResetBit(1<<9); | |
149 | res = TMath::Abs(fitFunc->GetParameter(2)); | |
150 | resError = TMath::Abs(fitFunc->GetParError(2)); | |
151 | return kTRUE; | |
152 | } | |
153 | return kFALSE; | |
154 | } | |
155 | ||
156 | void AliAnalysisTaskCheckESD::UserCreateOutputObjects() | |
157 | { | |
158 | // Create histograms | |
159 | AliInfo("AliAnalysisTaskCheckESD::UserCreateOutputObjects"); | |
160 | // Create output container | |
161 | fListOfHistos = new TList(); | |
162 | ||
163 | // efficiency and resolution histog | |
164 | Int_t nBinsPt = 15; | |
165 | Float_t minPt = 0.1; | |
166 | Float_t maxPt = 3.1; | |
167 | ||
168 | hGen = CreateHisto("hGen", "generated tracks", nBinsPt, minPt, maxPt, "p_{t} [GeV/c]", "N"); | |
169 | hRec = CreateHisto("hRec", "reconstructed tracks", nBinsPt, minPt, maxPt, "p_{t} [GeV/c]", "N"); | |
170 | hResPtInv = CreateHisto("hResPtInv", "", 100, -10, 10, "(p_{t,rec}^{-1}-p_{t,sim}^{-1}) / p_{t,sim}^{-1} [%]", "N"); | |
171 | hResPhi = CreateHisto("hResPhi", "", 100, -20, 20, "#phi_{rec}-#phi_{sim} [mrad]", "N"); | |
172 | hResTheta = CreateHisto("hResTheta", "", 100, -20, 20, "#theta_{rec}-#theta_{sim} [mrad]", "N"); | |
173 | ||
174 | // dE/dx and TOF | |
175 | hDEdxRight = new TH2F("hDEdxRight", "", 300, 0, 3, 100, 0, 400); | |
176 | hDEdxRight->SetStats(kFALSE); | |
177 | hDEdxRight->GetXaxis()->SetTitle("p [GeV/c]"); | |
178 | hDEdxRight->GetYaxis()->SetTitle("dE/dx_{TPC}"); | |
179 | hDEdxRight->SetMarkerStyle(kFullCircle); | |
180 | hDEdxRight->SetMarkerSize(0.4); | |
181 | hDEdxWrong = new TH2F("hDEdxWrong", "", 300, 0, 3, 100, 0, 400); | |
182 | hDEdxWrong->SetStats(kFALSE); | |
183 | hDEdxWrong->GetXaxis()->SetTitle("p [GeV/c]"); | |
184 | hDEdxWrong->GetYaxis()->SetTitle("dE/dx_{TPC}"); | |
185 | hDEdxWrong->SetMarkerStyle(kFullCircle); | |
186 | hDEdxWrong->SetMarkerSize(0.4); | |
187 | hDEdxWrong->SetMarkerColor(kRed); | |
188 | ||
189 | hResTOFRight = CreateHisto("hResTOFRight", "", 100, -1000, 1000, "t_{TOF}-t_{track} [ps]", "N"); | |
190 | hResTOFWrong = CreateHisto("hResTOFWrong", "", 100, -1000, 1000, "t_{TOF}-t_{track} [ps]", "N"); | |
191 | hResTOFWrong->SetLineColor(kRed); | |
192 | ||
193 | // calorimeters | |
8d071f77 | 194 | hEPHOS = CreateHisto("hEPHOS", "PHOS", 100, 0, 50, "E [GeV]", "N"); |
8858ee12 | 195 | hEEMCAL = CreateHisto("hEEMCAL", "EMCAL", 100, 0, 50, "E [GeV]", "N"); |
196 | ||
197 | // muons | |
198 | hPtMUON = CreateHisto("hPtMUON", "MUON", 100, 0, 20, "p_{t} [GeV/c]", "N"); | |
199 | ||
200 | // V0s and cascades | |
201 | hMassK0 = CreateHisto("hMassK0", "K^{0}", 100, 0.4, 0.6, "M(#pi^{+}#pi^{-}) [GeV/c^{2}]", "N"); | |
202 | hMassLambda = CreateHisto("hMassLambda", "#Lambda", 100, 1.0, 1.2, "M(p#pi^{-}) [GeV/c^{2}]", "N"); | |
203 | ||
204 | hMassLambdaBar = CreateHisto("hMassLambdaBar", "#bar{#Lambda}", 100, 1.0, 1.2, "M(#bar{p}#pi^{+}) [GeV/c^{2}]", "N"); | |
205 | hMassXi = CreateHisto("hMassXi", "#Xi", 100, 1.2, 1.5, "M(#Lambda#pi) [GeV/c^{2}]", "N"); | |
206 | hMassOmega = CreateHisto("hMassOmega", "#Omega", 100, 1.5, 1.8, "M(#LambdaK) [GeV/c^{2}]", "N"); | |
207 | hScalars = new TH1F("hScalars","Container of scalars",8,0,8); | |
208 | hArrayHist = new TH1F("hArrayHist","Container for Array", | |
209 | (AliPID::kSPECIES+1)*AliPID::kSPECIES,0,(AliPID::kSPECIES+1)*AliPID::kSPECIES); | |
210 | ||
211 | fListOfHistos->Add(hGen); | |
212 | fListOfHistos->Add(hRec); | |
213 | fListOfHistos->Add(hResPtInv); | |
214 | fListOfHistos->Add(hResPhi); | |
215 | fListOfHistos->Add(hResTheta); | |
216 | fListOfHistos->Add(hDEdxRight); | |
217 | fListOfHistos->Add(hDEdxWrong); | |
218 | fListOfHistos->Add(hResTOFRight); | |
219 | fListOfHistos->Add(hResTOFWrong); | |
220 | fListOfHistos->Add(hEPHOS); | |
221 | fListOfHistos->Add(hEEMCAL); | |
222 | fListOfHistos->Add(hPtMUON); | |
223 | fListOfHistos->Add(hMassK0); | |
224 | fListOfHistos->Add(hMassLambda); | |
225 | fListOfHistos->Add(hMassLambdaBar); | |
226 | fListOfHistos->Add(hMassXi); | |
227 | fListOfHistos->Add(hMassOmega); | |
228 | fListOfHistos->Add(hScalars); | |
229 | fListOfHistos->Add(hArrayHist); | |
230 | } | |
231 | ||
232 | void AliAnalysisTaskCheckESD::UserExec(Option_t *option) | |
233 | { | |
234 | const Int_t kXiMinus = 3312; | |
235 | const Int_t kOmegaMinus = 3334; | |
236 | // check the content of the ESD | |
237 | Double_t cutPtV0 = 0.3; | |
238 | Double_t cutPtCascade = 0.5; | |
239 | Float_t minPt = 0.1; | |
240 | ||
241 | // PID | |
242 | Int_t partCode[AliPID::kSPECIES] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; | |
243 | Double_t partFrac[AliPID::kSPECIES] = {0.01, 0.01, 0.85, 0.10, 0.05}; | |
244 | ||
245 | AliMCEvent* mcEvent = MCEvent(); | |
246 | if (!mcEvent) { | |
247 | Printf("ERROR: Could not retrieve MC event"); | |
248 | return; | |
249 | } | |
250 | ||
251 | //Primary vertex needed | |
252 | TArrayF vertex(3); | |
253 | mcEvent->GenEventHeader()->PrimaryVertex(vertex); | |
254 | ||
255 | TObjArray selParticles; //Use TClonesArray? | |
256 | TObjArray selV0s; | |
257 | TObjArray selCascades; | |
258 | for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {//Check this loop again | |
259 | AliMCParticle* track = mcEvent->GetTrack(iTracks); | |
260 | TParticle* particle = track->Particle(); | |
261 | if (!particle) continue; | |
262 | if (particle->Pt() < 0.001) continue; | |
263 | if (TMath::Abs(particle->Eta()) > 0.9) continue; | |
264 | TVector3 dVertex(particle->Vx() - vertex[0], particle->Vy() - vertex[1], particle->Vz() - vertex[2]); | |
265 | if (dVertex.Mag() > 0.0001) continue; | |
266 | ||
267 | switch (TMath::Abs(particle->GetPdgCode())) { | |
268 | case kElectron: | |
269 | case kMuonMinus: | |
270 | case kPiPlus: | |
271 | case kKPlus: | |
272 | case kProton: | |
273 | { | |
274 | if (particle->Pt() > minPt) { | |
275 | selParticles.Add(particle); | |
276 | hScalars->Fill(0); | |
277 | hGen->Fill(particle->Pt()); | |
278 | } | |
279 | break; | |
280 | } | |
281 | case kK0Short: | |
282 | case kLambda0: | |
283 | { | |
284 | if (particle->Pt() > cutPtV0) { | |
285 | hScalars->Fill(3); | |
286 | selV0s.Add(particle); | |
287 | } | |
288 | break; | |
289 | } | |
290 | case kXiMinus: | |
291 | case kOmegaMinus: | |
292 | { | |
293 | if (particle->Pt() > cutPtCascade) { | |
294 | hScalars->Fill(6); | |
295 | selCascades.Add(particle); | |
296 | } | |
297 | break; | |
298 | } | |
299 | default: break; | |
300 | } | |
301 | } | |
302 | ||
303 | // ESD information | |
304 | AliVEvent* event = InputEvent(); | |
305 | if (!event) { | |
306 | Printf("ERROR: Could not retrieve event"); | |
307 | return; | |
308 | } | |
309 | ||
310 | AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event); | |
311 | // loop over tracks | |
312 | for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) { | |
313 | AliESDtrack* track = esd->GetTrack(iTrack); | |
314 | ||
315 | // select tracks of selected particles | |
316 | Int_t label = TMath::Abs(track->GetLabel()); | |
317 | if (label > mcEvent->GetNumberOfTracks()) continue; // background | |
318 | AliMCParticle* mctrack = mcEvent->GetTrack(label); | |
319 | TParticle* particle = mctrack->Particle(); | |
320 | if (!selParticles.Contains(particle)) continue; | |
321 | if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) continue; | |
322 | if (track->GetConstrainedChi2() > 1e9) continue; | |
323 | selParticles.Remove(particle); // don't count multiple tracks | |
324 | ||
325 | hScalars->Fill(1); | |
326 | hRec->Fill(particle->Pt()); | |
327 | if (track->GetLabel() < 0) { | |
328 | hScalars->Fill(2); | |
329 | } | |
330 | ||
331 | // resolutions | |
332 | hResPtInv->Fill(100. * (TMath::Abs(track->GetSigned1Pt()) - 1./particle->Pt()) *particle->Pt()); | |
333 | hResPhi->Fill(1000. * (track->Phi() - particle->Phi())); | |
334 | hResTheta->Fill(1000. * (track->Theta() - particle->Theta())); | |
335 | ||
336 | // PID | |
337 | if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) continue; | |
338 | Int_t iGen = 5; | |
339 | ||
340 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
341 | if (TMath::Abs(particle->GetPdgCode()) == partCode[i]) iGen = i; | |
342 | } | |
343 | Double_t probability[AliPID::kSPECIES]; | |
344 | track->GetESDpid(probability); | |
345 | Double_t pMax = 0; | |
346 | Int_t iRec = 0; | |
347 | for (Int_t i = 0; i < AliPID::kSPECIES; i++) { | |
348 | probability[i] *= partFrac[i]; | |
349 | if (probability[i] > pMax) { | |
350 | pMax = probability[i]; | |
351 | iRec = i; | |
352 | } | |
353 | ||
354 | } | |
355 | hArrayHist->Fill(AliPID::kSPECIES*iGen + iRec); | |
356 | if (iGen == iRec) { | |
357 | hScalars->Fill(5); | |
358 | } | |
359 | ||
360 | // dE/dx and TOF | |
361 | Double_t time[AliPID::kSPECIES]; | |
362 | track->GetIntegratedTimes(time); | |
363 | if (iGen == iRec) { | |
364 | hDEdxRight->Fill(particle->P(), track->GetTPCsignal()); | |
365 | if ((track->GetStatus() & AliESDtrack::kTOFpid) != 0) { | |
366 | hResTOFRight->Fill(track->GetTOFsignal() - time[iRec]); | |
367 | } | |
368 | } | |
369 | else { | |
370 | hDEdxWrong->Fill(particle->P(), track->GetTPCsignal()); | |
371 | if ((track->GetStatus() & AliESDtrack::kTOFpid) != 0) { | |
372 | hResTOFWrong->Fill(track->GetTOFsignal() - time[iRec]); | |
373 | } | |
374 | } | |
375 | } | |
376 | ||
377 | // loop over muon tracks | |
378 | for (Int_t iTrack = 0; iTrack < esd->GetNumberOfMuonTracks(); iTrack++) { | |
379 | AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack); | |
380 | Double_t ptInv = TMath::Abs(muonTrack->GetInverseBendingMomentum()); | |
381 | if (ptInv > 0.001) { | |
382 | hPtMUON->Fill(1./ptInv); | |
383 | } | |
384 | } | |
385 | ||
386 | // loop over V0s | |
387 | for (Int_t iV0 = 0; iV0 < esd->GetNumberOfV0s(); iV0++) { | |
388 | AliESDv0* v0 = esd->GetV0(iV0); | |
389 | if (v0->GetOnFlyStatus()) continue; | |
390 | v0->ChangeMassHypothesis(kK0Short); | |
391 | hMassK0->Fill(v0->GetEffMass()); | |
392 | v0->ChangeMassHypothesis(kLambda0); | |
393 | hMassLambda->Fill(v0->GetEffMass()); | |
394 | v0->ChangeMassHypothesis(kLambda0Bar); | |
395 | hMassLambdaBar->Fill(v0->GetEffMass()); | |
396 | ||
397 | Int_t negLabel = TMath::Abs(esd->GetTrack(v0->GetNindex())->GetLabel()); | |
398 | if (negLabel > mcEvent->GetNumberOfTracks()) continue; // background | |
399 | AliMCParticle* track = mcEvent->GetTrack(negLabel); | |
400 | TParticle* particle = track->Particle(); | |
401 | Int_t negMother = particle->GetMother(0); | |
402 | if (negMother < 0) continue; | |
403 | Int_t posLabel = TMath::Abs(esd->GetTrack(v0->GetPindex())->GetLabel()); | |
404 | if (posLabel > mcEvent->GetNumberOfTracks()) continue; // background | |
405 | track = mcEvent->GetTrack(posLabel); | |
406 | particle = track->Particle(); | |
407 | Int_t posMother = particle->GetMother(0); | |
408 | if (negMother != posMother) continue; | |
409 | track = mcEvent->GetTrack(negMother); | |
410 | particle = track->Particle(); | |
411 | if (!selV0s.Contains(particle)) continue; | |
412 | selV0s.Remove(particle); | |
413 | hScalars->Fill(4); | |
414 | } | |
415 | ||
416 | // loop over Cascades | |
417 | for (Int_t iCascade = 0; iCascade < esd->GetNumberOfCascades(); iCascade++) { | |
418 | AliESDcascade* cascade = esd->GetCascade(iCascade); | |
419 | Double_t v0q; | |
420 | cascade->ChangeMassHypothesis(v0q,kXiMinus); | |
421 | hMassXi->Fill(cascade->GetEffMass()); | |
422 | cascade->ChangeMassHypothesis(v0q,kOmegaMinus); | |
423 | hMassOmega->Fill(cascade->GetEffMass()); | |
424 | ||
425 | Int_t negLabel = TMath::Abs(esd->GetTrack(cascade->GetNindex())->GetLabel()); | |
426 | if (negLabel > mcEvent->GetNumberOfTracks()) continue; // background | |
427 | AliMCParticle* track = mcEvent->GetTrack(negLabel); | |
428 | TParticle* particle = track->Particle(); | |
429 | Int_t negMother = particle->GetMother(0); | |
430 | if (negMother < 0) continue; | |
431 | Int_t posLabel = TMath::Abs(esd->GetTrack(cascade->GetPindex())->GetLabel()); | |
432 | if (posLabel > mcEvent->GetNumberOfTracks()) continue; // background | |
433 | track = mcEvent->GetTrack(posLabel); | |
434 | particle = track->Particle(); | |
435 | Int_t posMother = particle->GetMother(0); | |
436 | if (negMother != posMother) continue; | |
437 | track = mcEvent->GetTrack(negMother); | |
438 | particle = track->Particle(); | |
439 | Int_t v0Mother = particle->GetMother(0); | |
440 | if (v0Mother < 0) continue; | |
441 | Int_t bacLabel = TMath::Abs(esd->GetTrack(cascade->GetBindex())->GetLabel()); | |
442 | if (bacLabel > mcEvent->GetNumberOfTracks()) continue; // background | |
443 | track = mcEvent->GetTrack(bacLabel); | |
444 | particle = track->Particle(); | |
445 | Int_t bacMother = particle->GetMother(0); | |
446 | if (v0Mother != bacMother) continue; | |
447 | track = mcEvent->GetTrack(v0Mother); | |
448 | particle = track->Particle(); | |
449 | if (!selCascades.Contains(particle)) continue; | |
450 | selCascades.Remove(particle); | |
451 | hScalars->Fill(7); | |
452 | } | |
453 | ||
454 | // loop over the clusters | |
455 | for (Int_t iCluster=0; iCluster<esd->GetNumberOfCaloClusters(); iCluster++) { | |
456 | AliESDCaloCluster * clust = esd->GetCaloCluster(iCluster); | |
457 | if (clust->IsPHOS()) hEPHOS->Fill(clust->E()); | |
458 | if (clust->IsEMCAL()) hEEMCAL->Fill(clust->E()); | |
459 | } | |
460 | ||
461 | // Post output data. | |
462 | PostData(1, fListOfHistos); | |
463 | } | |
464 | ||
465 | void AliAnalysisTaskCheckESD::Terminate(Option_t *) | |
466 | { | |
467 | // check values | |
468 | Int_t checkNGenLow = 1; | |
469 | Double_t checkEffLow = 0.5; | |
470 | Double_t checkEffSigma = 3; | |
471 | Double_t checkFakeHigh = 0.5; | |
472 | Double_t checkFakeSigma = 3; | |
473 | ||
474 | Double_t checkResPtInvHigh = 5; | |
475 | Double_t checkResPtInvSigma = 3; | |
476 | Double_t checkResPhiHigh = 10; | |
477 | Double_t checkResPhiSigma = 3; | |
478 | Double_t checkResThetaHigh = 10; | |
479 | Double_t checkResThetaSigma = 3; | |
480 | ||
481 | Double_t checkPIDEffLow = 0.5; | |
482 | Double_t checkPIDEffSigma = 3; | |
483 | Double_t checkResTOFHigh = 500; | |
484 | Double_t checkResTOFSigma = 3; | |
485 | ||
486 | Double_t checkPHOSNLow = 5; | |
487 | Double_t checkPHOSEnergyLow = 0.3; | |
488 | Double_t checkPHOSEnergyHigh = 1.0; | |
489 | Double_t checkEMCALNLow = 50; | |
490 | Double_t checkEMCALEnergyLow = 0.05; | |
491 | Double_t checkEMCALEnergyHigh = 1.0; | |
492 | ||
493 | Double_t checkMUONNLow = 1; | |
494 | Double_t checkMUONPtLow = 0.5; | |
495 | Double_t checkMUONPtHigh = 10.; | |
496 | ||
497 | Double_t checkV0EffLow = 0.02; | |
498 | Double_t checkV0EffSigma = 3; | |
499 | Double_t checkCascadeEffLow = 0.01; | |
500 | Double_t checkCascadeEffSigma = 3; | |
501 | ||
502 | const char* partName[AliPID::kSPECIES+1] = {"electron", "muon", "pion", "kaon", "proton", "other"}; | |
503 | fListOfHistos = dynamic_cast<TList*>(GetOutputData(1)); | |
504 | if (!fListOfHistos) { | |
505 | Printf("ERROR: fListOfHistos not available"); | |
506 | return; | |
507 | } | |
508 | ||
509 | hGen = dynamic_cast<TH1F*>(fListOfHistos->At(0)); | |
510 | hRec = dynamic_cast<TH1F*>(fListOfHistos->At(1)); | |
511 | hResPtInv = dynamic_cast<TH1F*>(fListOfHistos->At(2)); | |
512 | hResPhi = dynamic_cast<TH1F*>(fListOfHistos->At(3)); | |
513 | hResTheta = dynamic_cast<TH1F*>(fListOfHistos->At(4)); | |
514 | hDEdxRight = dynamic_cast<TH2F*>(fListOfHistos->At(5)); | |
515 | hDEdxWrong = dynamic_cast<TH2F*>(fListOfHistos->At(6)); | |
516 | hResTOFRight = dynamic_cast<TH1F*>(fListOfHistos->At(7)); | |
517 | hResTOFWrong = dynamic_cast<TH1F*>(fListOfHistos->At(8)); | |
518 | hEPHOS = dynamic_cast<TH1F*>(fListOfHistos->At(9)); | |
519 | hEEMCAL = dynamic_cast<TH1F*>(fListOfHistos->At(10)); | |
520 | hPtMUON = dynamic_cast<TH1F*>(fListOfHistos->At(11)); | |
521 | hMassK0 = dynamic_cast<TH1F*>(fListOfHistos->At(12)); | |
522 | hMassLambda = dynamic_cast<TH1F*>(fListOfHistos->At(13)); | |
523 | hMassLambdaBar = dynamic_cast<TH1F*>(fListOfHistos->At(14)); | |
524 | hMassXi = dynamic_cast<TH1F*>(fListOfHistos->At(15)); | |
525 | hMassOmega = dynamic_cast<TH1F*>(fListOfHistos->At(16)); | |
526 | hScalars = dynamic_cast<TH1F*>(fListOfHistos->At(17)); | |
527 | hArrayHist = dynamic_cast<TH1F*>(fListOfHistos->At(18)); | |
528 | ||
529 | Int_t nGen = Int_t(hScalars->GetBinContent(1)); | |
530 | Int_t nRec = Int_t(hScalars->GetBinContent(2)); | |
531 | Int_t nFake = Int_t(hScalars->GetBinContent(3)); | |
532 | Int_t nGenV0s = Int_t(hScalars->GetBinContent(4)); | |
533 | Int_t nRecV0s = Int_t(hScalars->GetBinContent(5)); | |
534 | Int_t nIdentified = Int_t(hScalars->GetBinContent(6)); | |
535 | Int_t nGenCascades = Int_t(hScalars->GetBinContent(7)); | |
536 | Int_t nRecCascades = Int_t(hScalars->GetBinContent(8)); | |
537 | ||
538 | Int_t k = 1; | |
539 | ||
540 | Int_t identified[AliPID::kSPECIES+1][AliPID::kSPECIES]; | |
541 | for(Int_t i = 0; i < (AliPID::kSPECIES+1); i++) | |
542 | for(Int_t j = 0; j < AliPID::kSPECIES; j++) { | |
543 | identified[i][j] = Int_t(hArrayHist->GetBinContent(k)); | |
544 | k++; | |
545 | } | |
546 | ||
547 | ||
548 | // perform checks | |
549 | if (nGen < checkNGenLow) { | |
550 | Warning("CheckESD", "low number of generated particles: %d", Int_t(nGen)); | |
551 | } | |
552 | ||
553 | TH1F* hEff = CreateEffHisto(hGen, hRec); | |
554 | ||
555 | Info("CheckESD", "%d out of %d tracks reconstructed including %d " | |
556 | "fake tracks", nRec, nGen, nFake); | |
557 | if (nGen > 0) { | |
558 | // efficiency | |
559 | Double_t eff = nRec*1./nGen; | |
560 | Double_t effError = TMath::Sqrt(eff*(1.-eff) / nGen); | |
561 | Double_t fake = nFake*1./nGen; | |
562 | Double_t fakeError = TMath::Sqrt(fake*(1.-fake) / nGen); | |
563 | Info("CheckESD", "eff = (%.1f +- %.1f) %% fake = (%.1f +- %.1f) %%", | |
564 | 100.*eff, 100.*effError, 100.*fake, 100.*fakeError); | |
565 | if (eff < checkEffLow - checkEffSigma*effError) { | |
566 | Warning("CheckESD", "low efficiency: (%.1f +- %.1f) %%",100.*eff, 100.*effError); | |
567 | } | |
568 | if (fake > checkFakeHigh + checkFakeSigma*fakeError) { | |
569 | Warning("CheckESD", "high fake: (%.1f +- %.1f) %%",100.*fake, 100.*fakeError); | |
570 | } | |
571 | // resolutions | |
572 | Double_t res, resError; | |
573 | if (FitHisto(hResPtInv, res, resError)) { | |
574 | Info("CheckESD", "relative inverse pt resolution = (%.1f +- %.1f) %%",res, resError); | |
575 | if (res > checkResPtInvHigh + checkResPtInvSigma*resError) { | |
576 | Warning("CheckESD", "bad pt resolution: (%.1f +- %.1f) %%",res, resError); | |
577 | } | |
578 | } | |
579 | if (FitHisto(hResPhi, res, resError)) { | |
580 | Info("CheckESD", "phi resolution = (%.1f +- %.1f) mrad", res, resError); | |
581 | if (res > checkResPhiHigh + checkResPhiSigma*resError) { | |
582 | Warning("CheckESD", "bad phi resolution: (%.1f +- %.1f) mrad", | |
583 | res, resError); | |
584 | } | |
585 | } | |
586 | ||
587 | if (FitHisto(hResTheta, res, resError)) { | |
588 | Info("CheckESD", "theta resolution = (%.1f +- %.1f) mrad", | |
589 | res, resError); | |
590 | if (res > checkResThetaHigh + checkResThetaSigma*resError) { | |
591 | Warning("CheckESD", "bad theta resolution: (%.1f +- %.1f) mrad", | |
592 | res, resError); | |
593 | } | |
594 | } | |
595 | // PID | |
596 | if (nRec > 0) { | |
597 | Double_t eff = nIdentified*1./nRec; | |
598 | Double_t effError = TMath::Sqrt(eff*(1.-eff) / nRec); | |
599 | Info("CheckESD", "PID eff = (%.1f +- %.1f) %%", | |
600 | 100.*eff, 100.*effError); | |
601 | if (eff < checkPIDEffLow - checkPIDEffSigma*effError) { | |
602 | Warning("CheckESD", "low PID efficiency: (%.1f +- %.1f) %%", | |
603 | 100.*eff, 100.*effError); | |
604 | } | |
605 | } | |
606 | ||
607 | printf("%9s:", "gen\\rec"); | |
608 | for (Int_t iRec = 0; iRec < AliPID::kSPECIES; iRec++) { | |
609 | printf("%9s", partName[iRec]); | |
610 | } | |
611 | printf("\n"); | |
612 | for (Int_t iGen = 0; iGen < AliPID::kSPECIES+1; iGen++) { | |
613 | printf("%9s:", partName[iGen]); | |
614 | for (Int_t iRec = 0; iRec < AliPID::kSPECIES; iRec++) { | |
615 | printf("%9d", identified[iGen][iRec]); | |
616 | } | |
617 | printf("\n"); | |
618 | } | |
619 | ||
620 | if (FitHisto(hResTOFRight, res, resError)) { | |
621 | Info("CheckESD", "TOF resolution = (%.1f +- %.1f) ps", res, resError); | |
622 | if (res > checkResTOFHigh + checkResTOFSigma*resError) { | |
623 | Warning("CheckESD", "bad TOF resolution: (%.1f +- %.1f) ps", | |
624 | res, resError); | |
625 | } | |
626 | } | |
627 | ||
628 | // calorimeters | |
629 | if (hEPHOS->Integral() < checkPHOSNLow) { | |
630 | Warning("CheckESD", "low number of PHOS particles: %d", | |
631 | Int_t(hEPHOS->Integral())); | |
632 | } | |
633 | else { | |
634 | Double_t mean = hEPHOS->GetMean(); | |
635 | if (mean < checkPHOSEnergyLow) { | |
636 | Warning("CheckESD", "low mean PHOS energy: %.1f GeV", mean); | |
637 | } else if (mean > checkPHOSEnergyHigh) { | |
638 | Warning("CheckESD", "high mean PHOS energy: %.1f GeV", mean); | |
639 | } | |
640 | } | |
641 | ||
642 | if (hEEMCAL->Integral() < checkEMCALNLow) { | |
643 | Warning("CheckESD", "low number of EMCAL particles: %d", | |
644 | Int_t(hEEMCAL->Integral())); | |
645 | } | |
646 | else { | |
647 | Double_t mean = hEEMCAL->GetMean(); | |
648 | if (mean < checkEMCALEnergyLow) { | |
649 | Warning("CheckESD", "low mean EMCAL energy: %.1f GeV", mean); | |
650 | } | |
651 | else if (mean > checkEMCALEnergyHigh) { | |
652 | Warning("CheckESD", "high mean EMCAL energy: %.1f GeV", mean); | |
653 | } | |
654 | } | |
655 | ||
656 | // muons | |
657 | if (hPtMUON->Integral() < checkMUONNLow) { | |
658 | Warning("CheckESD", "low number of MUON particles: %d", | |
659 | Int_t(hPtMUON->Integral())); | |
660 | } | |
661 | else { | |
662 | Double_t mean = hPtMUON->GetMean(); | |
663 | if (mean < checkMUONPtLow) { | |
664 | Warning("CheckESD", "low mean MUON pt: %.1f GeV/c", mean); | |
665 | } | |
666 | else if (mean > checkMUONPtHigh) { | |
667 | Warning("CheckESD", "high mean MUON pt: %.1f GeV/c", mean); | |
668 | } | |
669 | } | |
670 | ||
671 | // V0s | |
672 | if (nGenV0s > 0) { | |
673 | Double_t eff = nRecV0s*1./nGenV0s; | |
674 | Double_t effError = TMath::Sqrt(eff*(1.-eff) / nGenV0s); | |
675 | if (effError == 0) effError = checkV0EffLow / TMath::Sqrt(1.*nGenV0s); | |
676 | Info("CheckESD", "V0 eff = (%.1f +- %.1f) %%", | |
677 | 100.*eff, 100.*effError); | |
678 | if (eff < checkV0EffLow - checkV0EffSigma*effError) { | |
679 | Warning("CheckESD", "low V0 efficiency: (%.1f +- %.1f) %%", | |
680 | 100.*eff, 100.*effError); | |
681 | } | |
682 | } | |
683 | ||
684 | // Cascades | |
685 | if (nGenCascades > 0) { | |
686 | Double_t eff = nRecCascades*1./nGenCascades; | |
687 | Double_t effError = TMath::Sqrt(eff*(1.-eff) / nGenCascades); | |
688 | if (effError == 0) effError = checkV0EffLow / | |
689 | TMath::Sqrt(1.*nGenCascades); | |
690 | Info("CheckESD", "Cascade eff = (%.1f +- %.1f) %%", | |
691 | 100.*eff, 100.*effError); | |
692 | if (eff < checkCascadeEffLow - checkCascadeEffSigma*effError) { | |
693 | Warning("CheckESD", "low Cascade efficiency: (%.1f +- %.1f) %%", | |
694 | 100.*eff, 100.*effError); | |
695 | } | |
696 | } | |
697 | ||
698 | } | |
699 | ||
700 | // draw the histograms if not in batch mode | |
701 | if (!gROOT->IsBatch()) { | |
702 | new TCanvas; | |
703 | hEff->DrawCopy(); | |
704 | new TCanvas; | |
705 | hResPtInv->DrawCopy("E"); | |
706 | new TCanvas; | |
707 | hResPhi->DrawCopy("E"); | |
708 | new TCanvas; | |
709 | hResTheta->DrawCopy("E"); | |
710 | new TCanvas; | |
711 | hDEdxRight->DrawCopy(); | |
712 | hDEdxWrong->DrawCopy("SAME"); | |
713 | new TCanvas; | |
714 | hResTOFRight->DrawCopy("E"); | |
715 | hResTOFWrong->DrawCopy("SAME"); | |
716 | new TCanvas; | |
717 | hEPHOS->DrawCopy("E"); | |
718 | new TCanvas; | |
719 | hEEMCAL->DrawCopy("E"); | |
720 | new TCanvas; | |
721 | hPtMUON->DrawCopy("E"); | |
722 | new TCanvas; | |
723 | hMassK0->DrawCopy("E"); | |
724 | new TCanvas; | |
725 | hMassLambda->DrawCopy("E"); | |
726 | new TCanvas; | |
727 | hMassLambdaBar->DrawCopy("E"); | |
728 | new TCanvas; | |
729 | hMassXi->DrawCopy("E"); | |
730 | new TCanvas; | |
731 | hMassOmega->DrawCopy("E"); | |
732 | } | |
733 | ||
734 | // write the output histograms to a file | |
735 | TFile* outputFile = TFile::Open("check.root", "recreate"); | |
736 | if (!outputFile || !outputFile->IsOpen()) | |
737 | { | |
738 | Error("CheckESD", "opening output file check.root failed"); | |
739 | return; | |
740 | } | |
741 | hEff->Write(); | |
742 | hResPtInv->Write(); | |
743 | hResPhi->Write(); | |
744 | hResTheta->Write(); | |
745 | hDEdxRight->Write(); | |
746 | hDEdxWrong->Write(); | |
747 | hResTOFRight->Write(); | |
748 | hResTOFWrong->Write(); | |
749 | hEPHOS->Write(); | |
750 | hEEMCAL->Write(); | |
751 | hPtMUON->Write(); | |
752 | hMassK0->Write(); | |
753 | hMassLambda->Write(); | |
754 | hMassLambdaBar->Write(); | |
755 | hMassXi->Write(); | |
756 | hMassOmega->Write(); | |
757 | outputFile->Close(); | |
758 | delete outputFile; | |
759 | ||
760 | // clean up | |
761 | /* delete hGen; | |
762 | delete hRec; | |
763 | delete hEff; | |
764 | delete hResPtInv; | |
765 | delete hResPhi; | |
766 | delete hResTheta; | |
767 | delete hDEdxRight; | |
768 | delete hDEdxWrong; | |
769 | delete hResTOFRight; | |
770 | delete hResTOFWrong; | |
771 | delete hEPHOS; | |
772 | delete hEEMCAL; | |
773 | delete hPtMUON; | |
774 | delete hMassK0; | |
775 | delete hMassLambda; | |
776 | delete hMassLambdaBar; | |
777 | delete hMassXi; | |
778 | delete hMassOmega; | |
779 | delete hScalars; | |
780 | delete hArrayHist; */ | |
781 | ||
782 | //delete esd; | |
783 | // result of check | |
784 | Info("CheckESD", "check of ESD was successfull"); | |
785 | } | |
786 | ||
787 | ||
788 | ||
789 | ||
790 | ||
791 | ||
792 | ||
793 | ||
794 | ||
795 | ||
796 | ||
797 | ||
798 | ||
799 | ||
800 | ||
801 | ||
802 | ||
803 | ||
804 | ||
805 | ||
806 | ||
807 | ||
808 | ||
809 | ||
810 | ||
811 | ||
812 | ||
813 | ||
814 | ||
815 | ||
816 | ||
817 | ||
818 | ||
819 | ||
820 | ||
821 | ||
822 | ||
823 | ||
824 | ||
825 | ||
826 | ||
827 | ||
828 | ||
829 | ||
830 | ||
831 | ||
832 | ||
833 | ||
834 | ||
835 | ||
836 | ||
837 | ||
838 | ||
839 | ||
840 | ||
841 | ||
842 | ||
843 | ||
844 | ||
845 | ||
846 | ||
847 |