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