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