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