]>
Commit | Line | Data |
---|---|---|
684b93f0 | 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 | //------------------------------ | |
17 | // Analysis task for quality-assurance | |
18 | // of forward detectors ESD | |
19 | // | |
20 | // 12/06/2009 cvetan.cheshkov@cern.ch | |
21 | //------------------------------ | |
22 | ||
23 | #include "TChain.h" | |
24 | #include "TROOT.h" | |
25 | #include "TFile.h" | |
26 | #include "TH1F.h" | |
27 | #include "TH2F.h" | |
28 | #include "TF1.h" | |
29 | #include "TCanvas.h" | |
30 | #include "TVector3.h" | |
31 | #include "TParticle.h" | |
8b04dae1 | 32 | #include "AliVParticle.h" |
33 | #include "AliMCParticle.h" | |
684b93f0 | 34 | #include "AliESDEvent.h" |
35 | #include "AliESDv0.h" | |
36 | #include "AliESDcascade.h" | |
37 | #include "AliESDMuonTrack.h" | |
38 | #include "AliESDCaloCluster.h" | |
39 | #include "AliRun.h" | |
40 | #include "AliMCEvent.h" | |
41 | #include "AliStack.h" | |
42 | #include "AliGenEventHeader.h" | |
43 | #include "AliPID.h" | |
44 | #include "AliESDVZERO.h" | |
45 | ||
46 | #include "AliAnaFwdDetsQA.h" | |
47 | ||
48 | ClassImp(AliAnaFwdDetsQA) | |
49 | ||
50 | AliAnaFwdDetsQA::AliAnaFwdDetsQA(): | |
51 | AliAnalysisTaskSE("AliAnaFwdDetsQA"), | |
52 | fListOfHistos(0), | |
53 | fT0vtxRec(0), | |
54 | fT0vtxRecGen(0), | |
55 | fT0time(0), | |
c543bd05 | 56 | fT0time2(0), |
684b93f0 | 57 | fT0mult(0), |
58 | fT0vtxRes(0), | |
c543bd05 | 59 | fT0ampl(0), |
684b93f0 | 60 | fV0a(0), |
61 | fV0c(0), | |
62 | fV0multA(0), | |
63 | fV0multC(0), | |
64 | fV0multAcorr(0), | |
65 | fV0multCcorr(0), | |
66 | fV0Acorr(0), | |
c543bd05 | 67 | fV0Ccorr(0), |
68 | fV0ampl(0) | |
684b93f0 | 69 | { |
70 | // Default constructor | |
71 | // Define input and output slots here | |
72 | // Input slot #0 works with a TChain | |
73 | DefineInput(0, TChain::Class()); | |
74 | // Output slot #1 TList | |
75 | DefineOutput(1, TList::Class()); | |
76 | } | |
77 | ||
78 | AliAnaFwdDetsQA::AliAnaFwdDetsQA(const char* name): | |
79 | AliAnalysisTaskSE(name), | |
80 | fListOfHistos(0), | |
81 | fT0vtxRec(0), | |
82 | fT0vtxRecGen(0), | |
83 | fT0time(0), | |
c543bd05 | 84 | fT0time2(0), |
684b93f0 | 85 | fT0mult(0), |
86 | fT0vtxRes(0), | |
c543bd05 | 87 | fT0ampl(0), |
684b93f0 | 88 | fV0a(0), |
89 | fV0c(0), | |
90 | fV0multA(0), | |
91 | fV0multC(0), | |
92 | fV0multAcorr(0), | |
93 | fV0multCcorr(0), | |
94 | fV0Acorr(0), | |
c543bd05 | 95 | fV0Ccorr(0), |
96 | fV0ampl(0) | |
684b93f0 | 97 | { |
98 | // Constructor | |
99 | AliInfo("Constructor AliAnaFwdDetsQA"); | |
100 | // Define input and output slots here | |
101 | // Input slot #0 works with a TChain | |
102 | DefineInput(0, TChain::Class()); | |
103 | // Output slot #1 TList | |
104 | DefineOutput(1, TList::Class()); | |
105 | } | |
106 | ||
107 | TH1F * AliAnaFwdDetsQA::CreateHisto(const char* name, const char* title,Int_t nBins, | |
108 | Double_t xMin, Double_t xMax, | |
109 | const char* xLabel, const char* yLabel) | |
110 | { | |
111 | // create a histogram | |
112 | TH1F* result = new TH1F(name, title, nBins, xMin, xMax); | |
113 | result->SetOption("E"); | |
114 | if (xLabel) result->GetXaxis()->SetTitle(xLabel); | |
115 | if (yLabel) result->GetYaxis()->SetTitle(yLabel); | |
116 | result->SetMarkerStyle(kFullCircle); | |
117 | return result; | |
118 | } | |
119 | ||
120 | TH1F *AliAnaFwdDetsQA::CreateEffHisto(const TH1F* hGen, const TH1F* hRec) | |
121 | { | |
122 | // create an efficiency histogram | |
123 | Int_t nBins = hGen->GetNbinsX(); | |
124 | TH1F* hEff = (TH1F*) hGen->Clone("hEff"); | |
125 | hEff->SetTitle(""); | |
126 | hEff->SetStats(kFALSE); | |
127 | hEff->SetMinimum(0.); | |
128 | hEff->SetMaximum(110.); | |
129 | hEff->GetYaxis()->SetTitle("#epsilon [%]"); | |
130 | ||
131 | for (Int_t iBin = 0; iBin <= nBins; iBin++) { | |
132 | Double_t nGenEff = hGen->GetBinContent(iBin); | |
133 | Double_t nRecEff = hRec->GetBinContent(iBin); | |
134 | if (nGenEff > 0) { | |
135 | Double_t eff = nRecEff/nGenEff; | |
136 | hEff->SetBinContent(iBin, 100. * eff); | |
137 | Double_t error = sqrt(eff*(1.-eff) / nGenEff); | |
138 | if (error == 0) error = 0.0001; | |
139 | hEff->SetBinError(iBin, 100. * error); | |
140 | } | |
141 | else { | |
142 | hEff->SetBinContent(iBin, -100.); | |
143 | hEff->SetBinError(iBin, 0); | |
144 | } | |
145 | } | |
146 | return hEff; | |
147 | } | |
148 | ||
149 | ||
150 | Bool_t AliAnaFwdDetsQA::FitHisto(TH1* histo, Double_t& res, Double_t& resError) | |
151 | { | |
152 | // fit a gaussian to a histogram | |
153 | static TF1* fitFunc = new TF1("fitFunc", "gaus"); | |
154 | fitFunc->SetLineWidth(2); | |
155 | fitFunc->SetFillStyle(0); | |
156 | Double_t maxFitRange = 2; | |
157 | ||
158 | if (histo->Integral() > 50) { | |
159 | Float_t mean = histo->GetMean(); | |
160 | Float_t rms = histo->GetRMS(); | |
161 | fitFunc->SetRange(mean - maxFitRange*rms, mean + maxFitRange*rms); | |
162 | fitFunc->SetParameters(mean, rms); | |
163 | histo->Fit(fitFunc, "QRI0"); | |
164 | histo->GetFunction("fitFunc")->ResetBit(1<<9); | |
165 | res = TMath::Abs(fitFunc->GetParameter(2)); | |
166 | resError = TMath::Abs(fitFunc->GetParError(2)); | |
167 | return kTRUE; | |
168 | } | |
169 | return kFALSE; | |
170 | } | |
171 | ||
172 | void AliAnaFwdDetsQA::UserCreateOutputObjects() | |
173 | { | |
174 | // Create histograms | |
175 | AliInfo("AliAnaFwdDetsQA::UserCreateOutputObjects"); | |
176 | // Create output container | |
177 | fListOfHistos = new TList(); | |
178 | ||
179 | fT0vtxRec = CreateHisto("hT0vtxRec", "Z vertex reconstructed with T0", 100, -25, 25, "Z_{vtx} [cm]", ""); | |
180 | fT0time = CreateHisto("hT0time", "Time0 reconstruction with T0", 5000, 10000, 20000, "t_{0} [ps]", ""); | |
c543bd05 | 181 | fT0time2 = CreateHisto("hT0time2", "Time0 reconstruction with T0 (mult > 10)", 5000, 10000, 20000, "t_{0} [ps]", ""); |
684b93f0 | 182 | fT0mult = CreateHisto("hT0mult", "Total reconstructed multiplicity in T0", 100, -0.5, 99.5, "Multiplicity", ""); |
183 | fT0vtxRes = CreateHisto("hT0vtxRes", "T0 Z vertex resolution", 100, -25, 25, "Delta(Z_{vtx}) [cm]", ""); | |
c543bd05 | 184 | fT0ampl = CreateHisto("hT0ampl","T0 multiplicity in single channel (all T0 channels)",400,-0.5,99.5); |
684b93f0 | 185 | |
186 | fT0vtxRecGen = new TH2F("hT0vtxRecGen", "T0 reconstructed vs generated Z vertex", 100, -25, 25, 100, -25, 25); | |
187 | fT0vtxRecGen->GetXaxis()->SetTitle("Generated Z vertex"); | |
188 | fT0vtxRecGen->GetYaxis()->SetTitle("Reconstructed Z vertex"); | |
189 | fT0vtxRecGen->SetMarkerStyle(kFullCircle); | |
190 | fT0vtxRecGen->SetMarkerSize(0.4); | |
191 | ||
192 | fV0a = CreateHisto("hV0a","Number of fired PMTs (V0A)",65,-0.5,64.5); | |
193 | fV0c = CreateHisto("hV0c","Number of fired PMTs (V0C)",65,-0.5,64.5); | |
194 | fV0multA = CreateHisto("hV0multA","Total reconstructed multiplicity (V0A)",100,0.,1000.); | |
195 | fV0multC = CreateHisto("hV0multC","Total reconstructed multiplicity (V0C)",100,0.,1000.); | |
c543bd05 | 196 | fV0ampl = CreateHisto("hV0ampl","V0 multiplicity in single channel (all V0 channels)",400,-0.5,99.5); |
684b93f0 | 197 | |
198 | fV0multAcorr = new TH2F("hV0multAcorr", "Reconstructed vs generated (primaries only) multiplicity (V0A)",100,0.,500.,100,0.,1000.); | |
199 | fV0multAcorr->GetXaxis()->SetTitle("# of primaries in V0A acceptance"); | |
200 | fV0multAcorr->GetYaxis()->SetTitle("Reconstructed mutiplicity in V0A"); | |
201 | fV0multCcorr = new TH2F("hV0multCcorr", "Reconstructed vs generated (primaries only) multiplicity (V0C)",100,0.,500.,100,0.,1000.); | |
202 | fV0multCcorr->GetXaxis()->SetTitle("# of primaries in V0C acceptance"); | |
203 | fV0multCcorr->GetYaxis()->SetTitle("Reconstructed mutiplicity in V0C"); | |
204 | ||
205 | fV0Acorr = new TH2F("hV0Acorr", "Number of fired PMTs vs generated (primaries only) multiplicity (V0A)",100,0.,500.,65,-0.5,64.5); | |
206 | fV0Acorr->GetXaxis()->SetTitle("# of primaries in V0A acceptance"); | |
207 | fV0Acorr->GetYaxis()->SetTitle("Number of fired PMTs in V0A"); | |
208 | fV0Ccorr = new TH2F("hV0Ccorr", "Number of fired PMTs vs generated (primaries only) multiplicity (V0C)",100,0.,500.,65,-0.5,64.5); | |
209 | fV0Ccorr->GetXaxis()->SetTitle("# of primaries in V0C acceptance"); | |
210 | fV0Ccorr->GetYaxis()->SetTitle("Number of fired PMTs in V0C"); | |
211 | ||
212 | fListOfHistos->Add(fT0vtxRec); | |
213 | fListOfHistos->Add(fT0time); | |
214 | fListOfHistos->Add(fT0mult); | |
215 | fListOfHistos->Add(fT0vtxRecGen); | |
216 | fListOfHistos->Add(fT0vtxRes); | |
217 | fListOfHistos->Add(fV0a); | |
218 | fListOfHistos->Add(fV0c); | |
219 | fListOfHistos->Add(fV0multA); | |
220 | fListOfHistos->Add(fV0multC); | |
221 | fListOfHistos->Add(fV0multAcorr); | |
222 | fListOfHistos->Add(fV0multCcorr); | |
223 | fListOfHistos->Add(fV0Acorr); | |
224 | fListOfHistos->Add(fV0Ccorr); | |
c543bd05 | 225 | fListOfHistos->Add(fT0ampl); |
226 | fListOfHistos->Add(fV0ampl); | |
227 | fListOfHistos->Add(fT0time2); | |
684b93f0 | 228 | } |
229 | ||
230 | void AliAnaFwdDetsQA::UserExec(Option_t */*option*/) | |
231 | { | |
232 | ||
233 | AliMCEvent* mcEvent = MCEvent(); | |
234 | if (!mcEvent) { | |
235 | Printf("ERROR: Could not retrieve MC event"); | |
236 | return; | |
237 | } | |
238 | ||
239 | //Primary vertex needed | |
240 | TArrayF mcvtx(3); | |
241 | mcEvent->GenEventHeader()->PrimaryVertex(mcvtx); | |
242 | ||
243 | AliStack *stack = mcEvent->Stack(); | |
244 | if (!stack) { | |
245 | Printf("ERROR: Could not retrieve MC stack"); | |
246 | return; | |
247 | } | |
248 | ||
249 | Int_t nV0A = 0; | |
250 | Int_t nV0C = 0; | |
251 | for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {//Check this loop again | |
252 | if (!stack->IsPhysicalPrimary(iTracks)) continue; | |
8b04dae1 | 253 | AliMCParticle* track = (AliMCParticle*)mcEvent->GetTrack(iTracks); |
684b93f0 | 254 | TParticle* particle = track->Particle(); |
255 | if (!particle) continue; | |
c543bd05 | 256 | if (track->Charge() == 0) continue; |
684b93f0 | 257 | Double_t eta = particle->Eta(); |
258 | if (eta > 2.8 && eta < 5.1) { | |
259 | nV0A++; | |
260 | } | |
261 | if (eta > -3.7 && eta < -1.7) { | |
262 | nV0C++; | |
263 | } | |
264 | } | |
265 | ||
266 | // ESD information | |
267 | AliVEvent* event = InputEvent(); | |
268 | if (!event) { | |
269 | Printf("ERROR: Could not retrieve event"); | |
270 | return; | |
271 | } | |
272 | ||
273 | AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event); | |
274 | const AliESDTZERO* esdT0 = esd->GetESDTZERO(); | |
275 | Double_t t0zvtx = esdT0->GetT0zVertex(); | |
276 | Double_t t0time = esdT0->GetT0(); | |
277 | ||
278 | fT0vtxRec->Fill(t0zvtx); | |
279 | fT0time->Fill(t0time); | |
280 | fT0vtxRecGen->Fill(mcvtx[2],t0zvtx); | |
281 | t0zvtx *= -1.0; | |
282 | fT0vtxRes->Fill(t0zvtx - mcvtx[2]); | |
283 | ||
284 | Double_t t0mult = 0; | |
285 | for(Int_t i = 0; i < 24; i++) { | |
286 | t0mult += esdT0->GetT0amplitude()[i]; | |
c543bd05 | 287 | fT0ampl->Fill(esdT0->GetT0amplitude()[i]); |
684b93f0 | 288 | } |
289 | ||
290 | fT0mult->Fill(t0mult); | |
c543bd05 | 291 | if (t0mult > 10) |
292 | fT0time2->Fill(t0time); | |
684b93f0 | 293 | |
294 | AliESDVZERO* esdV0 = esd->GetVZEROData(); | |
295 | fV0a->Fill(esdV0->GetNbPMV0A()); | |
296 | fV0c->Fill(esdV0->GetNbPMV0C()); | |
297 | fV0multA->Fill(esdV0->GetMTotV0A()); | |
298 | fV0multC->Fill(esdV0->GetMTotV0C()); | |
299 | ||
300 | fV0multAcorr->Fill((Float_t)nV0A,esdV0->GetMTotV0A()); | |
301 | fV0multCcorr->Fill((Float_t)nV0C,esdV0->GetMTotV0C()); | |
302 | fV0Acorr->Fill((Float_t)nV0A,(Float_t)esdV0->GetNbPMV0A()); | |
303 | fV0Ccorr->Fill((Float_t)nV0C,(Float_t)esdV0->GetNbPMV0C()); | |
304 | ||
c543bd05 | 305 | for(Int_t i = 0; i < 64; i++) { |
306 | fV0ampl->Fill(esdV0->GetMultiplicity(i)); | |
307 | } | |
684b93f0 | 308 | // Post output data. |
309 | PostData(1, fListOfHistos); | |
310 | } | |
311 | ||
312 | void AliAnaFwdDetsQA::Terminate(Option_t *) | |
313 | { | |
314 | fListOfHistos = dynamic_cast<TList*>(GetOutputData(1)); | |
315 | if (!fListOfHistos) { | |
316 | Printf("ERROR: fListOfHistos not available"); | |
317 | return; | |
318 | } | |
319 | ||
320 | fT0vtxRec = dynamic_cast<TH1F*>(fListOfHistos->At(0)); | |
321 | fT0time = dynamic_cast<TH1F*>(fListOfHistos->At(1)); | |
322 | fT0mult = dynamic_cast<TH1F*>(fListOfHistos->At(2)); | |
323 | fT0vtxRecGen = dynamic_cast<TH2F*>(fListOfHistos->At(3)); | |
324 | fT0vtxRes = dynamic_cast<TH1F*>(fListOfHistos->At(4)); | |
325 | ||
326 | fV0a = dynamic_cast<TH1F*>(fListOfHistos->At(5)); | |
327 | fV0c = dynamic_cast<TH1F*>(fListOfHistos->At(6)); | |
328 | fV0multA = dynamic_cast<TH1F*>(fListOfHistos->At(7)); | |
329 | fV0multC = dynamic_cast<TH1F*>(fListOfHistos->At(8)); | |
330 | fV0multAcorr = dynamic_cast<TH2F*>(fListOfHistos->At(9)); | |
331 | fV0multCcorr = dynamic_cast<TH2F*>(fListOfHistos->At(10)); | |
332 | fV0Acorr = dynamic_cast<TH2F*>(fListOfHistos->At(11)); | |
333 | fV0Ccorr = dynamic_cast<TH2F*>(fListOfHistos->At(12)); | |
334 | ||
c543bd05 | 335 | fT0ampl = dynamic_cast<TH1F*>(fListOfHistos->At(13)); |
336 | fV0ampl = dynamic_cast<TH1F*>(fListOfHistos->At(14)); | |
337 | fT0time2 = dynamic_cast<TH1F*>(fListOfHistos->At(15)); | |
338 | ||
339 | ||
684b93f0 | 340 | // draw the histograms if not in batch mode |
341 | if (!gROOT->IsBatch()) { | |
342 | new TCanvas; | |
343 | fT0vtxRec->DrawCopy("E"); | |
344 | new TCanvas; | |
345 | fT0time->DrawCopy("E"); | |
346 | new TCanvas; | |
347 | fT0mult->DrawCopy("E"); | |
348 | new TCanvas; | |
349 | fT0vtxRes->DrawCopy("E"); | |
350 | new TCanvas; | |
351 | fT0vtxRecGen->DrawCopy(); | |
352 | new TCanvas; | |
353 | fV0a->DrawCopy("E"); | |
354 | new TCanvas; | |
355 | fV0c->DrawCopy("E"); | |
356 | new TCanvas; | |
357 | fV0multA->DrawCopy("E"); | |
358 | new TCanvas; | |
359 | fV0multC->DrawCopy("E"); | |
360 | new TCanvas; | |
361 | fV0multAcorr->DrawCopy(); | |
362 | new TCanvas; | |
363 | fV0multCcorr->DrawCopy(); | |
364 | new TCanvas; | |
365 | fV0Acorr->DrawCopy(); | |
366 | new TCanvas; | |
367 | fV0Ccorr->DrawCopy(); | |
c543bd05 | 368 | new TCanvas; |
369 | fT0ampl->DrawCopy("E"); | |
370 | new TCanvas; | |
371 | fV0ampl->DrawCopy("E"); | |
372 | new TCanvas; | |
373 | fT0time2->DrawCopy("E"); | |
684b93f0 | 374 | } |
375 | ||
376 | // write the output histograms to a file | |
377 | TFile* outputFile = TFile::Open("FwdDetsQA.root", "recreate"); | |
378 | if (!outputFile || !outputFile->IsOpen()) | |
379 | { | |
380 | Error("AliAnaFwdDetsQA", "opening output file FwdDetsQA.root failed"); | |
381 | return; | |
382 | } | |
383 | fT0vtxRec->Write(); | |
384 | fT0time->Write(); | |
385 | fT0mult->Write(); | |
386 | fT0vtxRecGen->Write(); | |
387 | fT0vtxRes->Write(); | |
c543bd05 | 388 | fT0ampl->Write(); |
389 | fT0time2->Write(); | |
684b93f0 | 390 | |
391 | fV0a->Write(); | |
392 | fV0c->Write(); | |
393 | fV0multA->Write(); | |
394 | fV0multC->Write(); | |
395 | fV0multAcorr->Write(); | |
396 | fV0multCcorr->Write(); | |
397 | fV0Acorr->Write(); | |
398 | fV0Ccorr->Write(); | |
c543bd05 | 399 | fV0ampl->Write(); |
684b93f0 | 400 | |
401 | outputFile->Close(); | |
402 | delete outputFile; | |
403 | ||
404 | //delete esd; | |
405 | Info("AliAnaFwdDetsQA", "Successfully finished"); | |
406 | } | |
407 | ||
408 | ||
409 | ||
410 | ||
411 | ||
412 | ||
413 | ||
414 | ||
415 | ||
416 | ||
417 | ||
418 | ||
419 | ||
420 | ||
421 | ||
422 | ||
423 | ||
424 | ||
425 | ||
426 | ||
427 | ||
428 | ||
429 | ||
430 | ||
431 | ||
432 | ||
433 | ||
434 | ||
435 | ||
436 | ||
437 | ||
438 | ||
439 | ||
440 | ||
441 | ||
442 | ||
443 | ||
444 | ||
445 | ||
446 | ||
447 | ||
448 | ||
449 | ||
450 | ||
451 | ||
452 | ||
453 | ||
454 | ||
455 | ||
456 | ||
457 | ||
458 | ||
459 | ||
460 | ||
461 | ||
462 | ||
463 | ||
464 | ||
465 | ||
466 | ||
467 | ||
468 |