5e00415f |
1 | #if !defined( __CINT__) || defined(__MAKECINT__) |
2 | #include <TFile.h> |
3 | #include <TError.h> |
4 | #include <TH1.h> |
5 | #include <TH2.h> |
6 | #include <TF1.h> |
7 | #include <TCanvas.h> |
8 | #include <TVector3.h> |
9 | #include <TPDGCode.h> |
10 | |
11 | #include "AliRunLoader.h" |
12 | #include "AliLoader.h" |
13 | #include "AliESD.h" |
14 | #include "AliRun.h" |
15 | #include "AliStack.h" |
16 | #include "AliHeader.h" |
17 | #include "AliGenEventHeader.h" |
18 | #endif |
19 | |
20 | |
21 | TH1F* CreateHisto(const char* name, const char* title, |
22 | Int_t nBins, Double_t xMin, Double_t xMax, |
23 | const char* xLabel = NULL, const char* yLabel = NULL) |
24 | { |
25 | // create a histogram |
26 | |
27 | TH1F* result = new TH1F(name, title, nBins, xMin, xMax); |
28 | result->SetOption("E"); |
29 | if (xLabel) result->GetXaxis()->SetTitle(xLabel); |
30 | if (yLabel) result->GetYaxis()->SetTitle(yLabel); |
31 | result->SetMarkerStyle(kFullCircle); |
32 | return result; |
33 | } |
34 | |
35 | TH1F* CreateEffHisto(TH1F* hGen, TH1F* hRec) |
36 | { |
37 | // create an efficiency histogram |
38 | |
39 | Int_t nBins = hGen->GetNbinsX(); |
40 | TH1F* hEff = (TH1F*) hGen->Clone("hEff"); |
41 | hEff->SetTitle(""); |
42 | hEff->SetStats(kFALSE); |
43 | hEff->SetMinimum(0.); |
44 | hEff->SetMaximum(110.); |
45 | hEff->GetYaxis()->SetTitle("#epsilon [%]"); |
46 | |
47 | for (Int_t iBin = 0; iBin <= nBins; iBin++) { |
48 | Double_t nGen = hGen->GetBinContent(iBin); |
49 | Double_t nRec = hRec->GetBinContent(iBin); |
50 | if (nGen > 0) { |
51 | Double_t eff = nRec/nGen; |
52 | hEff->SetBinContent(iBin, 100. * eff); |
53 | Double_t error = sqrt(eff*(1.-eff) / nGen); |
54 | if (error == 0) error = 0.0001; |
55 | hEff->SetBinError(iBin, 100. * error); |
56 | } else { |
57 | hEff->SetBinContent(iBin, -100.); |
58 | hEff->SetBinError(iBin, 0); |
59 | } |
60 | } |
61 | |
62 | return hEff; |
63 | } |
64 | |
65 | Bool_t FitHisto(TH1* histo, Double_t& res, Double_t& resError) |
66 | { |
67 | // fit a gaussian to a histogram |
68 | |
69 | static TF1* fitFunc = new TF1("fitFunc", "gaus"); |
70 | fitFunc->SetLineWidth(2); |
71 | fitFunc->SetFillStyle(0); |
72 | Double_t maxFitRange = 2; |
73 | |
74 | if (histo->Integral() > 50) { |
75 | Float_t mean = histo->GetMean(); |
76 | Float_t rms = histo->GetRMS(); |
77 | fitFunc->SetRange(mean - maxFitRange*rms, mean + maxFitRange*rms); |
78 | fitFunc->SetParameters(mean, rms); |
79 | histo->Fit(fitFunc, "QRI0"); |
80 | histo->GetFunction("fitFunc")->ResetBit(1<<9); |
81 | res = TMath::Abs(fitFunc->GetParameter(2)); |
82 | resError = TMath::Abs(fitFunc->GetParError(2)); |
83 | return kTRUE; |
84 | } |
85 | |
86 | return kFALSE; |
87 | } |
88 | |
89 | |
90 | Bool_t CheckESD(const char* gAliceFileName = "galice.root", |
91 | const char* esdFileName = "AliESDs.root") |
92 | { |
93 | // check the content of the ESD |
94 | |
95 | // check values |
96 | Int_t checkNGenLow = 1; |
97 | |
98 | Double_t checkEffLow = 0.5; |
99 | Double_t checkEffSigma = 3; |
100 | Double_t checkFakeHigh = 0.5; |
101 | Double_t checkFakeSigma = 3; |
102 | |
103 | Double_t checkResPtInvHigh = 5; |
104 | Double_t checkResPtInvSigma = 3; |
105 | Double_t checkResPhiHigh = 10; |
106 | Double_t checkResPhiSigma = 3; |
107 | Double_t checkResThetaHigh = 10; |
108 | Double_t checkResThetaSigma = 3; |
109 | |
110 | Double_t checkPIDEffLow = 0.5; |
111 | Double_t checkPIDEffSigma = 3; |
112 | Double_t checkResTOFHigh = 500; |
113 | Double_t checkResTOFSigma = 3; |
114 | |
115 | Double_t checkPHOSNLow = 5; |
116 | Double_t checkPHOSEnergyLow = 0.3; |
117 | Double_t checkPHOSEnergyHigh = 1.0; |
118 | Double_t checkEMCALNLow = 50; |
119 | Double_t checkEMCALEnergyLow = 0.05; |
120 | Double_t checkEMCALEnergyHigh = 1.0; |
121 | |
122 | Double_t checkMUONNLow = 1; |
123 | Double_t checkMUONPtLow = 0.5; |
124 | Double_t checkMUONPtHigh = 10.; |
125 | |
126 | Double_t cutPtV0 = 0.3; |
127 | Double_t checkV0EffLow = 0.02; |
128 | Double_t checkV0EffSigma = 3; |
129 | Double_t cutPtCascade = 0.5; |
130 | Double_t checkCascadeEffLow = 0.01; |
131 | Double_t checkCascadeEffSigma = 3; |
132 | |
133 | // open run loader and load gAlice, kinematics and header |
134 | AliRunLoader* runLoader = AliRunLoader::Open(gAliceFileName); |
135 | if (!runLoader) { |
136 | Error("CheckESD", "getting run loader from file %s failed", |
137 | gAliceFileName); |
138 | return kFALSE; |
139 | } |
140 | runLoader->LoadgAlice(); |
141 | gAlice = runLoader->GetAliRun(); |
142 | if (!gAlice) { |
143 | Error("CheckESD", "no galice object found"); |
144 | return kFALSE; |
145 | } |
146 | runLoader->LoadKinematics(); |
147 | runLoader->LoadHeader(); |
148 | |
149 | // open the ESD file |
150 | TFile* esdFile = TFile::Open(esdFileName); |
151 | if (!esdFile || !esdFile->IsOpen()) { |
152 | Error("CheckESD", "opening ESD file %s failed", esdFileName); |
153 | return kFALSE; |
154 | } |
155 | |
156 | // efficienc and resolution histograms |
157 | Int_t nBinsPt = 15; |
158 | Float_t minPt = 0.1; |
159 | Float_t maxPt = 3.1; |
160 | TH1F* hGen = CreateHisto("hGen", "generated tracks", |
161 | nBinsPt, minPt, maxPt, "p_{t} [GeV/c]", "N"); |
162 | TH1F* hRec = CreateHisto("hRec", "reconstructed tracks", |
163 | nBinsPt, minPt, maxPt, "p_{t} [GeV/c]", "N"); |
164 | Int_t nGen = 0; |
165 | Int_t nRec = 0; |
166 | Int_t nFake = 0; |
167 | |
168 | TH1F* hResPtInv = CreateHisto("hResPtInv", "", 100, -10, 10, |
169 | "(p_{t,rec}^{-1}-p_{t,sim}^{-1}) / p_{t,sim}^{-1} [%]", "N"); |
170 | TH1F* hResPhi = CreateHisto("hResPhi", "", 100, -20, 20, |
171 | "#phi_{rec}-#phi_{sim} [mrad]", "N"); |
172 | TH1F* hResTheta = CreateHisto("hResTheta", "", 100, -20, 20, |
173 | "#theta_{rec}-#theta_{sim} [mrad]", "N"); |
174 | |
175 | // PID |
176 | Int_t partCode[AliESDtrack::kSPECIES] = |
177 | {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; |
178 | const char* partName[AliESDtrack::kSPECIES+1] = |
179 | {"electron", "muon", "pion", "kaon", "proton", "other"}; |
180 | Double_t partFrac[AliESDtrack::kSPECIES] = |
181 | {0.01, 0.01, 0.85, 0.10, 0.05}; |
182 | Int_t identified[AliESDtrack::kSPECIES+1][AliESDtrack::kSPECIES]; |
183 | for (Int_t iGen = 0; iGen < AliESDtrack::kSPECIES+1; iGen++) { |
184 | for (Int_t iRec = 0; iRec < AliESDtrack::kSPECIES; iRec++) { |
185 | identified[iGen][iRec] = 0; |
186 | } |
187 | } |
188 | Int_t nIdentified = 0; |
189 | |
190 | // dE/dx and TOF |
191 | TH2F* hDEdxRight = new TH2F("hDEdxRight", "", 300, 0, 3, 100, 0, 400); |
192 | hDEdxRight->SetStats(kFALSE); |
193 | hDEdxRight->GetXaxis()->SetTitle("p [GeV/c]"); |
194 | hDEdxRight->GetYaxis()->SetTitle("dE/dx_{TPC}"); |
195 | hDEdxRight->SetMarkerStyle(kFullCircle); |
196 | hDEdxRight->SetMarkerSize(0.4); |
197 | TH2F* hDEdxWrong = new TH2F("hDEdxWrong", "", 300, 0, 3, 100, 0, 400); |
198 | hDEdxWrong->SetStats(kFALSE); |
199 | hDEdxWrong->GetXaxis()->SetTitle("p [GeV/c]"); |
200 | hDEdxWrong->GetYaxis()->SetTitle("dE/dx_{TPC}"); |
201 | hDEdxWrong->SetMarkerStyle(kFullCircle); |
202 | hDEdxWrong->SetMarkerSize(0.4); |
203 | hDEdxWrong->SetMarkerColor(kRed); |
204 | TH1F* hResTOFRight = CreateHisto("hResTOFRight", "", 100, -1000, 1000, |
205 | "t_{TOF}-t_{track} [ps]", "N"); |
206 | TH1F* hResTOFWrong = CreateHisto("hResTOFWrong", "", 100, -1000, 1000, |
207 | "t_{TOF}-t_{track} [ps]", "N"); |
208 | hResTOFWrong->SetLineColor(kRed); |
209 | |
210 | // calorimeters |
211 | TH1F* hEPHOS = CreateHisto("hEPHOS", "PHOS", 100, 0, 5, "E [GeV]", "N"); |
212 | TH1F* hEEMCAL = CreateHisto("hEEMCAL", "EMCAL", 100, 0, 2, "E [GeV]", "N"); |
213 | |
214 | // muons |
215 | TH1F* hPtMUON = CreateHisto("hPtMUON", "MUON", 100, 0, 20, |
216 | "p_{t} [GeV/c]", "N"); |
217 | |
218 | // V0s and cascades |
219 | TH1F* hMassK0 = CreateHisto("hMassK0", "K^{0}", 100, 0.4, 0.6, |
220 | "M(#pi^{+}#pi^{-}) [GeV/c^{2}]", "N"); |
221 | TH1F* hMassLambda = CreateHisto("hMassLambda", "#Lambda", 100, 1.0, 1.2, |
222 | "M(p#pi^{-}) [GeV/c^{2}]", "N"); |
223 | TH1F* hMassLambdaBar = CreateHisto("hMassLambdaBar", "#bar{#Lambda}", |
224 | 100, 1.0, 1.2, |
225 | "M(#bar{p}#pi^{+}) [GeV/c^{2}]", "N"); |
226 | Int_t nGenV0s = 0; |
227 | Int_t nRecV0s = 0; |
228 | TH1F* hMassXi = CreateHisto("hMassXi", "#Xi", 100, 1.2, 1.5, |
229 | "M(#Lambda#pi) [GeV/c^{2}]", "N"); |
230 | TH1F* hMassOmega = CreateHisto("hMassOmega", "#Omega", 100, 1.5, 1.8, |
231 | "M(#LambdaK) [GeV/c^{2}]", "N"); |
232 | Int_t nGenCascades = 0; |
233 | Int_t nRecCascades = 0; |
234 | |
235 | // loop over events |
236 | for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) { |
237 | runLoader->GetEvent(iEvent); |
238 | |
239 | // select simulated primary particles, V0s and cascades |
240 | AliStack* stack = gAlice->Stack(); |
241 | Int_t nParticles = stack->GetNtrack(); |
242 | TArrayF vertex(3); |
243 | runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex); |
244 | TObjArray selParticles; |
245 | TObjArray selV0s; |
246 | TObjArray selCascades; |
247 | for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) { |
248 | TParticle* particle = stack->Particle(iParticle); |
249 | if (!particle) continue; |
250 | if (particle->Pt() < 0.001) continue; |
251 | if (TMath::Abs(particle->Eta()) > 0.9) continue; |
252 | TVector3 dVertex(particle->Vx() - vertex[0], |
253 | particle->Vy() - vertex[1], |
254 | particle->Vz() - vertex[2]); |
255 | if (dVertex.Mag() > 0.0001) continue; |
256 | |
257 | switch (TMath::Abs(particle->GetPdgCode())) { |
258 | case kElectron: |
259 | case kMuonMinus: |
260 | case kPiPlus: |
261 | case kKPlus: |
262 | case kProton: { |
263 | if (particle->Pt() > minPt) { |
264 | selParticles.Add(particle); |
265 | nGen++; |
266 | hGen->Fill(particle->Pt()); |
267 | } |
268 | break; |
269 | } |
270 | case kK0Short: |
271 | case kLambda0: { |
272 | if (particle->Pt() > cutPtV0) { |
273 | nGenV0s++; |
274 | selV0s.Add(particle); |
275 | } |
276 | break; |
277 | } |
278 | case kXiMinus: |
279 | case kOmegaMinus: { |
280 | if (particle->Pt() > cutPtCascade) { |
281 | nGenCascades++; |
282 | selCascades.Add(particle); |
283 | } |
284 | break; |
285 | } |
286 | default: break; |
287 | } |
288 | } |
289 | |
290 | // get the event summary data |
291 | char esdName[256]; |
292 | sprintf(esdName, "ESD%d", iEvent); |
293 | AliESD* esd = (AliESD*) esdFile->Get(esdName); |
294 | if (!esd) { |
295 | Error("CheckESD", "no ESD object found for event %d", iEvent); |
296 | return kFALSE; |
297 | } |
298 | |
299 | // loop over tracks |
300 | for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) { |
301 | AliESDtrack* track = esd->GetTrack(iTrack); |
302 | |
303 | // select tracks of selected particles |
304 | Int_t label = TMath::Abs(track->GetLabel()); |
305 | if (label > stack->GetNtrack()) continue; // background |
306 | TParticle* particle = stack->Particle(label); |
307 | if (!selParticles.Contains(particle)) continue; |
308 | if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) continue; |
309 | if (track->GetConstrainedChi2() > 1e9) continue; |
310 | selParticles.Remove(particle); // don't count multiple tracks |
311 | |
312 | nRec++; |
313 | hRec->Fill(particle->Pt()); |
314 | if (track->GetLabel() < 0) nFake++; |
315 | |
316 | // resolutions |
317 | Double_t p[3]; |
318 | track->GetConstrainedPxPyPz(p); |
319 | TVector3 pTrack(p); |
320 | hResPtInv->Fill(100. * (1./pTrack.Pt() - 1./particle->Pt()) * |
321 | particle->Pt()); |
322 | hResPhi->Fill(1000. * (pTrack.Phi() - particle->Phi())); |
323 | hResTheta->Fill(1000. * (pTrack.Theta() - particle->Theta())); |
324 | |
325 | // PID |
326 | if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) continue; |
327 | Int_t iGen = 5; |
328 | for (Int_t i = 0; i < AliESDtrack::kSPECIES; i++) { |
329 | if (TMath::Abs(particle->GetPdgCode()) == partCode[i]) iGen = i; |
330 | } |
331 | Double_t probability[AliESDtrack::kSPECIES]; |
332 | track->GetESDpid(probability); |
333 | Double_t pMax = 0; |
334 | Int_t iRec = 0; |
335 | for (Int_t i = 0; i < AliESDtrack::kSPECIES; i++) { |
336 | probability[i] *= partFrac[i]; |
337 | if (probability[i] > pMax) { |
338 | pMax = probability[i]; |
339 | iRec = i; |
340 | } |
341 | } |
342 | identified[iGen][iRec]++; |
343 | if (iGen == iRec) nIdentified++; |
344 | |
345 | // dE/dx and TOF |
346 | Double_t time[AliESDtrack::kSPECIES]; |
347 | track->GetIntegratedTimes(time); |
348 | if (iGen == iRec) { |
349 | hDEdxRight->Fill(pTrack.Mag(), track->GetTPCsignal()); |
350 | if ((track->GetStatus() & AliESDtrack::kTOFpid) != 0) { |
351 | hResTOFRight->Fill(track->GetTOFsignal() - time[iRec]); |
352 | } |
353 | } else { |
354 | hDEdxWrong->Fill(pTrack.Mag(), track->GetTPCsignal()); |
355 | if ((track->GetStatus() & AliESDtrack::kTOFpid) != 0) { |
356 | hResTOFWrong->Fill(track->GetTOFsignal() - time[iRec]); |
357 | } |
358 | } |
359 | } |
360 | |
361 | // loop over calo tracks |
362 | for (Int_t iTrack = 0; iTrack < esd->GetNumberOfCaloTracks(); iTrack++) { |
363 | AliESDCaloTrack* track = esd->GetCaloTrack(iTrack); |
364 | TParticle* recParticle = track->GetRecParticle(); |
365 | if (recParticle->InheritsFrom("AliPHOSRecParticle")) { |
366 | hEPHOS->Fill(recParticle->Energy()); |
367 | } else if (recParticle->InheritsFrom("AliEMCALRecParticle")) { |
368 | hEEMCAL->Fill(recParticle->Energy()); |
369 | } else { |
370 | Warning("CheckESD", "unknown calo particle"); |
371 | recParticle->Dump(); |
372 | } |
373 | } |
374 | |
375 | // loop over muon tracks |
376 | for (Int_t iTrack = 0; iTrack < esd->GetNumberOfMuonTracks(); iTrack++) { |
377 | AliESDMuonTrack* track = esd->GetMuonTrack(iTrack); |
378 | Double_t ptInv = TMath::Abs(track->GetInverseBendingMomentum()); |
379 | if (ptInv > 0.001) { |
380 | hPtMUON->Fill(1./ptInv); |
381 | } |
382 | } |
383 | |
384 | // loop over V0s |
385 | for (Int_t iV0 = 0; iV0 < esd->GetNumberOfV0s(); iV0++) { |
386 | AliESDv0* v0 = esd->GetV0(iV0); |
387 | switch (v0->GetPdgCode()) { |
388 | case kK0Short : hMassK0->Fill(v0->GetEffMass()); break; |
389 | case kLambda0 : hMassLambda->Fill(v0->GetEffMass()); break; |
390 | case kLambda0Bar: hMassLambdaBar->Fill(v0->GetEffMass()); break; |
391 | default : break; |
392 | } |
393 | |
394 | Int_t negLabel = TMath::Abs(esd->GetTrack(v0->GetNindex())->GetLabel()); |
395 | if (negLabel > stack->GetNtrack()) continue; // background |
396 | Int_t negMother = stack->Particle(negLabel)->GetMother(0); |
397 | if (negMother < 0) continue; |
398 | Int_t posLabel = TMath::Abs(esd->GetTrack(v0->GetNindex())->GetLabel()); |
399 | if (posLabel > stack->GetNtrack()) continue; // background |
400 | Int_t posMother = stack->Particle(posLabel)->GetMother(0); |
401 | if (negMother != posMother) continue; |
402 | TParticle* particle = stack->Particle(negMother); |
403 | if (!selV0s.Contains(particle)) continue; |
404 | selV0s.Remove(particle); |
405 | nRecV0s++; |
406 | } |
407 | |
408 | // loop over Cascades |
409 | for (Int_t iCascade = 0; iCascade < esd->GetNumberOfCascades(); |
410 | iCascade++) { |
411 | AliESDcascade* cascade = esd->GetCascade(iCascade); |
412 | switch (TMath::Abs(cascade->GetPdgCode())) { |
413 | case kXiMinus : hMassXi->Fill(cascade->GetEffMass()); break; |
414 | case kOmegaMinus: hMassOmega->Fill(cascade->GetEffMass()); break; |
415 | default : break; |
416 | } |
417 | |
418 | Int_t negLabel = TMath::Abs(esd->GetTrack(cascade->GetNindex()) |
419 | ->GetLabel()); |
420 | if (negLabel > stack->GetNtrack()) continue; // background |
421 | Int_t negMother = stack->Particle(negLabel)->GetMother(0); |
422 | if (negMother < 0) continue; |
423 | Int_t posLabel = TMath::Abs(esd->GetTrack(cascade->GetNindex()) |
424 | ->GetLabel()); |
425 | if (posLabel > stack->GetNtrack()) continue; // background |
426 | Int_t posMother = stack->Particle(posLabel)->GetMother(0); |
427 | if (negMother != posMother) continue; |
428 | Int_t v0Mother = stack->Particle(negMother)->GetMother(0); |
429 | if (v0Mother < 0) continue; |
430 | Int_t bacLabel = TMath::Abs(esd->GetTrack(cascade->GetBindex()) |
431 | ->GetLabel()); |
432 | if (bacLabel > stack->GetNtrack()) continue; // background |
433 | Int_t bacMother = stack->Particle(bacLabel)->GetMother(0); |
434 | if (v0Mother != bacMother) continue; |
435 | TParticle* particle = stack->Particle(v0Mother); |
436 | if (!selCascades.Contains(particle)) continue; |
437 | selCascades.Remove(particle); |
438 | nRecCascades++; |
439 | } |
440 | } |
441 | |
442 | // perform checks |
443 | if (nGen < checkNGenLow) { |
444 | Warning("CheckESD", "low number of generated particles: %d", Int_t(nGen)); |
445 | } |
446 | |
447 | TH1F* hEff = CreateEffHisto(hGen, hRec); |
448 | |
449 | Info("CheckESD", "%d out of %d tracks reconstructed including %d " |
450 | "fake tracks", nRec, nGen, nFake); |
451 | if (nGen > 0) { |
452 | // efficiency |
453 | Double_t eff = nRec*1./nGen; |
454 | Double_t effError = TMath::Sqrt(eff*(1.-eff) / nGen); |
455 | Double_t fake = nFake*1./nGen; |
456 | Double_t fakeError = TMath::Sqrt(fake*(1.-fake) / nGen); |
457 | Info("CheckESD", "eff = (%.1f +- %.1f) %% fake = (%.1f +- %.1f) %%", |
458 | 100.*eff, 100.*effError, 100.*fake, 100.*fakeError); |
459 | |
460 | if (eff < checkEffLow - checkEffSigma*effError) { |
461 | Warning("CheckESD", "low efficiency: (%.1f +- %.1f) %%", |
462 | 100.*eff, 100.*effError); |
463 | } |
464 | if (fake > checkFakeHigh + checkFakeSigma*fakeError) { |
465 | Warning("CheckESD", "high fake: (%.1f +- %.1f) %%", |
466 | 100.*fake, 100.*fakeError); |
467 | } |
468 | |
469 | // resolutions |
470 | Double_t res, resError; |
471 | if (FitHisto(hResPtInv, res, resError)) { |
472 | Info("CheckESD", "relative inverse pt resolution = (%.1f +- %.1f) %%", |
473 | res, resError); |
474 | if (res > checkResPtInvHigh + checkResPtInvSigma*resError) { |
475 | Warning("CheckESD", "bad pt resolution: (%.1f +- %.1f) %%", |
476 | res, resError); |
477 | } |
478 | } |
479 | |
480 | if (FitHisto(hResPhi, res, resError)) { |
481 | Info("CheckESD", "phi resolution = (%.1f +- %.1f) mrad", res, resError); |
482 | if (res > checkResPhiHigh + checkResPhiSigma*resError) { |
483 | Warning("CheckESD", "bad phi resolution: (%.1f +- %.1f) mrad", |
484 | res, resError); |
485 | } |
486 | } |
487 | |
488 | if (FitHisto(hResTheta, res, resError)) { |
489 | Info("CheckESD", "theta resolution = (%.1f +- %.1f) mrad", |
490 | res, resError); |
491 | if (res > checkResThetaHigh + checkResThetaSigma*resError) { |
492 | Warning("CheckESD", "bad theta resolution: (%.1f +- %.1f) mrad", |
493 | res, resError); |
494 | } |
495 | } |
496 | |
497 | // PID |
498 | if (nRec > 0) { |
499 | Double_t eff = nIdentified*1./nRec; |
500 | Double_t effError = TMath::Sqrt(eff*(1.-eff) / nRec); |
501 | Info("CheckESD", "PID eff = (%.1f +- %.1f) %%", |
502 | 100.*eff, 100.*effError); |
503 | if (eff < checkPIDEffLow - checkPIDEffSigma*effError) { |
504 | Warning("CheckESD", "low PID efficiency: (%.1f +- %.1f) %%", |
505 | 100.*eff, 100.*effError); |
506 | } |
507 | } |
508 | |
509 | printf("%9s:", "gen\\rec"); |
510 | for (Int_t iRec = 0; iRec < AliESDtrack::kSPECIES; iRec++) { |
511 | printf("%9s", partName[iRec]); |
512 | } |
513 | printf("\n"); |
514 | for (Int_t iGen = 0; iGen < AliESDtrack::kSPECIES+1; iGen++) { |
515 | printf("%9s:", partName[iGen]); |
516 | for (Int_t iRec = 0; iRec < AliESDtrack::kSPECIES; iRec++) { |
517 | printf("%9d", identified[iGen][iRec]); |
518 | } |
519 | printf("\n"); |
520 | } |
521 | |
522 | if (FitHisto(hResTOFRight, res, resError)) { |
523 | Info("CheckESD", "TOF resolution = (%.1f +- %.1f) ps", res, resError); |
524 | if (res > checkResTOFHigh + checkResTOFSigma*resError) { |
525 | Warning("CheckESD", "bad TOF resolution: (%.1f +- %.1f) ps", |
526 | res, resError); |
527 | } |
528 | } |
529 | |
530 | // calorimeters |
531 | if (hEPHOS->Integral() < checkPHOSNLow) { |
532 | Warning("CheckESD", "low number of PHOS particles: %d", |
533 | Int_t(hEPHOS->Integral())); |
534 | } else { |
535 | Double_t mean = hEPHOS->GetMean(); |
536 | if (mean < checkPHOSEnergyLow) { |
537 | Warning("CheckESD", "low mean PHOS energy: %.1f GeV", mean); |
538 | } else if (mean > checkPHOSEnergyHigh) { |
539 | Warning("CheckESD", "high mean PHOS energy: %.1f GeV", mean); |
540 | } |
541 | } |
542 | |
543 | if (hEEMCAL->Integral() < checkEMCALNLow) { |
544 | Warning("CheckESD", "low number of EMCAL particles: %d", |
545 | Int_t(hEEMCAL->Integral())); |
546 | } else { |
547 | Double_t mean = hEEMCAL->GetMean(); |
548 | if (mean < checkEMCALEnergyLow) { |
549 | Warning("CheckESD", "low mean EMCAL energy: %.1f GeV", mean); |
550 | } else if (mean > checkEMCALEnergyHigh) { |
551 | Warning("CheckESD", "high mean EMCAL energy: %.1f GeV", mean); |
552 | } |
553 | } |
554 | |
555 | // muons |
556 | if (hPtMUON->Integral() < checkMUONNLow) { |
557 | Warning("CheckESD", "low number of MUON particles: %d", |
558 | Int_t(hPtMUON->Integral())); |
559 | } else { |
560 | Double_t mean = hPtMUON->GetMean(); |
561 | if (mean < checkMUONPtLow) { |
562 | Warning("CheckESD", "low mean MUON pt: %.1f GeV/c", mean); |
563 | } else if (mean > checkMUONPtHigh) { |
564 | Warning("CheckESD", "high mean MUON pt: %.1f GeV/c", mean); |
565 | } |
566 | } |
567 | |
568 | // V0s |
569 | if (nGenV0s > 0) { |
570 | Double_t eff = nRecV0s*1./nGenV0s; |
571 | Double_t effError = TMath::Sqrt(eff*(1.-eff) / nGenV0s); |
572 | if (effError == 0) effError = checkV0EffLow / TMath::Sqrt(1.*nGenV0s); |
573 | Info("CheckESD", "V0 eff = (%.1f +- %.1f) %%", |
574 | 100.*eff, 100.*effError); |
575 | if (eff < checkV0EffLow - checkV0EffSigma*effError) { |
576 | Warning("CheckESD", "low V0 efficiency: (%.1f +- %.1f) %%", |
577 | 100.*eff, 100.*effError); |
578 | } |
579 | } |
580 | |
581 | // Cascades |
582 | if (nGenCascades > 0) { |
583 | Double_t eff = nRecCascades*1./nGenCascades; |
584 | Double_t effError = TMath::Sqrt(eff*(1.-eff) / nGenCascades); |
585 | if (effError == 0) effError = checkV0EffLow / |
586 | TMath::Sqrt(1.*nGenCascades); |
587 | Info("CheckESD", "Cascade eff = (%.1f +- %.1f) %%", |
588 | 100.*eff, 100.*effError); |
589 | if (eff < checkCascadeEffLow - checkCascadeEffSigma*effError) { |
590 | Warning("CheckESD", "low Cascade efficiency: (%.1f +- %.1f) %%", |
591 | 100.*eff, 100.*effError); |
592 | } |
593 | } |
594 | } |
595 | |
596 | // draw the histograms if not in batch mode |
597 | if (!gROOT->IsBatch()) { |
598 | new TCanvas; |
599 | hEff->DrawCopy(); |
600 | new TCanvas; |
601 | hResPtInv->DrawCopy("E"); |
602 | new TCanvas; |
603 | hResPhi->DrawCopy("E"); |
604 | new TCanvas; |
605 | hResTheta->DrawCopy("E"); |
606 | new TCanvas; |
607 | hDEdxRight->DrawCopy(); |
608 | hDEdxWrong->DrawCopy("SAME"); |
609 | new TCanvas; |
610 | hResTOFRight->DrawCopy("E"); |
611 | hResTOFWrong->DrawCopy("SAME"); |
612 | new TCanvas; |
613 | hEPHOS->DrawCopy("E"); |
614 | new TCanvas; |
615 | hEEMCAL->DrawCopy("E"); |
616 | new TCanvas; |
617 | hPtMUON->DrawCopy("E"); |
618 | new TCanvas; |
619 | hMassK0->DrawCopy("E"); |
620 | new TCanvas; |
621 | hMassLambda->DrawCopy("E"); |
622 | new TCanvas; |
623 | hMassLambdaBar->DrawCopy("E"); |
624 | new TCanvas; |
625 | hMassXi->DrawCopy("E"); |
626 | new TCanvas; |
627 | hMassOmega->DrawCopy("E"); |
628 | } |
629 | |
630 | // write the output histograms to a file |
631 | TFile* outputFile = TFile::Open("check.root", "recreate"); |
632 | if (!outputFile || !outputFile->IsOpen()) { |
633 | Error("CheckESD", "opening output file check.root failed"); |
634 | return kFALSE; |
635 | } |
636 | hEff->Write(); |
637 | hResPtInv->Write(); |
638 | hResPhi->Write(); |
639 | hResTheta->Write(); |
640 | hDEdxRight->Write(); |
641 | hDEdxWrong->Write(); |
642 | hResTOFRight->Write(); |
643 | hResTOFWrong->Write(); |
644 | hEPHOS->Write(); |
645 | hEEMCAL->Write(); |
646 | hPtMUON->Write(); |
647 | hMassK0->Write(); |
648 | hMassLambda->Write(); |
649 | hMassLambdaBar->Write(); |
650 | hMassXi->Write(); |
651 | hMassOmega->Write(); |
652 | outputFile->Close(); |
653 | delete outputFile; |
654 | |
655 | // clean up |
656 | delete hGen; |
657 | delete hRec; |
658 | delete hEff; |
659 | delete hResPtInv; |
660 | delete hResPhi; |
661 | delete hResTheta; |
662 | delete hDEdxRight; |
663 | delete hDEdxWrong; |
664 | delete hResTOFRight; |
665 | delete hResTOFWrong; |
666 | delete hEPHOS; |
667 | delete hEEMCAL; |
668 | delete hPtMUON; |
669 | delete hMassK0; |
670 | delete hMassLambda; |
671 | delete hMassLambdaBar; |
672 | delete hMassXi; |
673 | delete hMassOmega; |
674 | |
675 | esdFile->Close(); |
676 | delete esdFile; |
677 | |
678 | runLoader->UnloadHeader(); |
679 | runLoader->UnloadKinematics(); |
680 | delete runLoader; |
681 | |
682 | // result of check |
683 | Info("CheckESD", "check of ESD was successfull"); |
684 | return kTRUE; |
685 | } |