Fix Coverity 15062, 15063
[u/mrichter/AliRoot.git] / PWG1 / AliAnaFwdDetsQA.cxx
CommitLineData
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
48ClassImp(AliAnaFwdDetsQA)
49
50AliAnaFwdDetsQA::AliAnaFwdDetsQA():
51AliAnalysisTaskSE("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
78AliAnaFwdDetsQA::AliAnaFwdDetsQA(const char* name):
79AliAnalysisTaskSE(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
684b93f0 99 // Define input and output slots here
100 // Input slot #0 works with a TChain
fc6c8a4a 101 AliInfo("Constructor AliAnaFwdDetsQA");
684b93f0 102 DefineInput(0, TChain::Class());
103 // Output slot #1 TList
104 DefineOutput(1, TList::Class());
105}
106
107TH1F * 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{
fc6c8a4a 111 // helper method which can be used
112 // in order to create a histogram
684b93f0 113 TH1F* result = new TH1F(name, title, nBins, xMin, xMax);
114 result->SetOption("E");
115 if (xLabel) result->GetXaxis()->SetTitle(xLabel);
116 if (yLabel) result->GetYaxis()->SetTitle(yLabel);
117 result->SetMarkerStyle(kFullCircle);
118 return result;
119}
120
121TH1F *AliAnaFwdDetsQA::CreateEffHisto(const TH1F* hGen, const TH1F* hRec)
122{
fc6c8a4a 123 // helper method which can be used
124 // in order create an efficiency histogram
684b93f0 125 Int_t nBins = hGen->GetNbinsX();
126 TH1F* hEff = (TH1F*) hGen->Clone("hEff");
127 hEff->SetTitle("");
128 hEff->SetStats(kFALSE);
129 hEff->SetMinimum(0.);
130 hEff->SetMaximum(110.);
131 hEff->GetYaxis()->SetTitle("#epsilon [%]");
132
133 for (Int_t iBin = 0; iBin <= nBins; iBin++) {
134 Double_t nGenEff = hGen->GetBinContent(iBin);
135 Double_t nRecEff = hRec->GetBinContent(iBin);
136 if (nGenEff > 0) {
137 Double_t eff = nRecEff/nGenEff;
138 hEff->SetBinContent(iBin, 100. * eff);
139 Double_t error = sqrt(eff*(1.-eff) / nGenEff);
fc6c8a4a 140 if (error < 1e-12) error = 0.0001;
684b93f0 141 hEff->SetBinError(iBin, 100. * error);
142 }
143 else {
144 hEff->SetBinContent(iBin, -100.);
145 hEff->SetBinError(iBin, 0);
146 }
147 }
148 return hEff;
149}
150
151
152Bool_t AliAnaFwdDetsQA::FitHisto(TH1* histo, Double_t& res, Double_t& resError)
153{
fc6c8a4a 154 // fit a gaussian to
155 // a histogram
684b93f0 156 static TF1* fitFunc = new TF1("fitFunc", "gaus");
157 fitFunc->SetLineWidth(2);
158 fitFunc->SetFillStyle(0);
159 Double_t maxFitRange = 2;
160
161 if (histo->Integral() > 50) {
162 Float_t mean = histo->GetMean();
163 Float_t rms = histo->GetRMS();
164 fitFunc->SetRange(mean - maxFitRange*rms, mean + maxFitRange*rms);
165 fitFunc->SetParameters(mean, rms);
166 histo->Fit(fitFunc, "QRI0");
167 histo->GetFunction("fitFunc")->ResetBit(1<<9);
168 res = TMath::Abs(fitFunc->GetParameter(2));
169 resError = TMath::Abs(fitFunc->GetParError(2));
170 return kTRUE;
171 }
172 return kFALSE;
173}
174
175void AliAnaFwdDetsQA::UserCreateOutputObjects()
176{
177 // Create histograms
684b93f0 178 // Create output container
fc6c8a4a 179 AliInfo("AliAnaFwdDetsQA::UserCreateOutputObjects");
684b93f0 180 fListOfHistos = new TList();
181
182 fT0vtxRec = CreateHisto("hT0vtxRec", "Z vertex reconstructed with T0", 100, -25, 25, "Z_{vtx} [cm]", "");
183 fT0time = CreateHisto("hT0time", "Time0 reconstruction with T0", 5000, 10000, 20000, "t_{0} [ps]", "");
c543bd05 184 fT0time2 = CreateHisto("hT0time2", "Time0 reconstruction with T0 (mult > 10)", 5000, 10000, 20000, "t_{0} [ps]", "");
684b93f0 185 fT0mult = CreateHisto("hT0mult", "Total reconstructed multiplicity in T0", 100, -0.5, 99.5, "Multiplicity", "");
186 fT0vtxRes = CreateHisto("hT0vtxRes", "T0 Z vertex resolution", 100, -25, 25, "Delta(Z_{vtx}) [cm]", "");
c543bd05 187 fT0ampl = CreateHisto("hT0ampl","T0 multiplicity in single channel (all T0 channels)",400,-0.5,99.5);
684b93f0 188
189 fT0vtxRecGen = new TH2F("hT0vtxRecGen", "T0 reconstructed vs generated Z vertex", 100, -25, 25, 100, -25, 25);
190 fT0vtxRecGen->GetXaxis()->SetTitle("Generated Z vertex");
191 fT0vtxRecGen->GetYaxis()->SetTitle("Reconstructed Z vertex");
192 fT0vtxRecGen->SetMarkerStyle(kFullCircle);
193 fT0vtxRecGen->SetMarkerSize(0.4);
194
195 fV0a = CreateHisto("hV0a","Number of fired PMTs (V0A)",65,-0.5,64.5);
196 fV0c = CreateHisto("hV0c","Number of fired PMTs (V0C)",65,-0.5,64.5);
197 fV0multA = CreateHisto("hV0multA","Total reconstructed multiplicity (V0A)",100,0.,1000.);
198 fV0multC = CreateHisto("hV0multC","Total reconstructed multiplicity (V0C)",100,0.,1000.);
c543bd05 199 fV0ampl = CreateHisto("hV0ampl","V0 multiplicity in single channel (all V0 channels)",400,-0.5,99.5);
684b93f0 200
201 fV0multAcorr = new TH2F("hV0multAcorr", "Reconstructed vs generated (primaries only) multiplicity (V0A)",100,0.,500.,100,0.,1000.);
202 fV0multAcorr->GetXaxis()->SetTitle("# of primaries in V0A acceptance");
203 fV0multAcorr->GetYaxis()->SetTitle("Reconstructed mutiplicity in V0A");
204 fV0multCcorr = new TH2F("hV0multCcorr", "Reconstructed vs generated (primaries only) multiplicity (V0C)",100,0.,500.,100,0.,1000.);
205 fV0multCcorr->GetXaxis()->SetTitle("# of primaries in V0C acceptance");
206 fV0multCcorr->GetYaxis()->SetTitle("Reconstructed mutiplicity in V0C");
207
208 fV0Acorr = new TH2F("hV0Acorr", "Number of fired PMTs vs generated (primaries only) multiplicity (V0A)",100,0.,500.,65,-0.5,64.5);
209 fV0Acorr->GetXaxis()->SetTitle("# of primaries in V0A acceptance");
210 fV0Acorr->GetYaxis()->SetTitle("Number of fired PMTs in V0A");
211 fV0Ccorr = new TH2F("hV0Ccorr", "Number of fired PMTs vs generated (primaries only) multiplicity (V0C)",100,0.,500.,65,-0.5,64.5);
212 fV0Ccorr->GetXaxis()->SetTitle("# of primaries in V0C acceptance");
213 fV0Ccorr->GetYaxis()->SetTitle("Number of fired PMTs in V0C");
214
215 fListOfHistos->Add(fT0vtxRec);
216 fListOfHistos->Add(fT0time);
217 fListOfHistos->Add(fT0mult);
218 fListOfHistos->Add(fT0vtxRecGen);
219 fListOfHistos->Add(fT0vtxRes);
220 fListOfHistos->Add(fV0a);
221 fListOfHistos->Add(fV0c);
222 fListOfHistos->Add(fV0multA);
223 fListOfHistos->Add(fV0multC);
224 fListOfHistos->Add(fV0multAcorr);
225 fListOfHistos->Add(fV0multCcorr);
226 fListOfHistos->Add(fV0Acorr);
227 fListOfHistos->Add(fV0Ccorr);
c543bd05 228 fListOfHistos->Add(fT0ampl);
229 fListOfHistos->Add(fV0ampl);
230 fListOfHistos->Add(fT0time2);
684b93f0 231}
232
233void AliAnaFwdDetsQA::UserExec(Option_t */*option*/)
234{
fc6c8a4a 235 // The analysis code
236 // goes here
684b93f0 237 AliMCEvent* mcEvent = MCEvent();
238 if (!mcEvent) {
239 Printf("ERROR: Could not retrieve MC event");
240 return;
241 }
242
243 //Primary vertex needed
244 TArrayF mcvtx(3);
245 mcEvent->GenEventHeader()->PrimaryVertex(mcvtx);
246
247 AliStack *stack = mcEvent->Stack();
248 if (!stack) {
249 Printf("ERROR: Could not retrieve MC stack");
250 return;
251 }
252
253 Int_t nV0A = 0;
254 Int_t nV0C = 0;
255 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {//Check this loop again
256 if (!stack->IsPhysicalPrimary(iTracks)) continue;
8b04dae1 257 AliMCParticle* track = (AliMCParticle*)mcEvent->GetTrack(iTracks);
684b93f0 258 TParticle* particle = track->Particle();
259 if (!particle) continue;
c543bd05 260 if (track->Charge() == 0) continue;
684b93f0 261 Double_t eta = particle->Eta();
262 if (eta > 2.8 && eta < 5.1) {
263 nV0A++;
264 }
265 if (eta > -3.7 && eta < -1.7) {
266 nV0C++;
267 }
268 }
269
270 // ESD information
271 AliVEvent* event = InputEvent();
272 if (!event) {
273 Printf("ERROR: Could not retrieve event");
274 return;
275 }
276
277 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
278 const AliESDTZERO* esdT0 = esd->GetESDTZERO();
279 Double_t t0zvtx = esdT0->GetT0zVertex();
280 Double_t t0time = esdT0->GetT0();
281
282 fT0vtxRec->Fill(t0zvtx);
283 fT0time->Fill(t0time);
284 fT0vtxRecGen->Fill(mcvtx[2],t0zvtx);
285 t0zvtx *= -1.0;
286 fT0vtxRes->Fill(t0zvtx - mcvtx[2]);
287
288 Double_t t0mult = 0;
289 for(Int_t i = 0; i < 24; i++) {
290 t0mult += esdT0->GetT0amplitude()[i];
c543bd05 291 fT0ampl->Fill(esdT0->GetT0amplitude()[i]);
684b93f0 292 }
293
294 fT0mult->Fill(t0mult);
c543bd05 295 if (t0mult > 10)
296 fT0time2->Fill(t0time);
684b93f0 297
298 AliESDVZERO* esdV0 = esd->GetVZEROData();
299 fV0a->Fill(esdV0->GetNbPMV0A());
300 fV0c->Fill(esdV0->GetNbPMV0C());
301 fV0multA->Fill(esdV0->GetMTotV0A());
302 fV0multC->Fill(esdV0->GetMTotV0C());
303
304 fV0multAcorr->Fill((Float_t)nV0A,esdV0->GetMTotV0A());
305 fV0multCcorr->Fill((Float_t)nV0C,esdV0->GetMTotV0C());
306 fV0Acorr->Fill((Float_t)nV0A,(Float_t)esdV0->GetNbPMV0A());
307 fV0Ccorr->Fill((Float_t)nV0C,(Float_t)esdV0->GetNbPMV0C());
308
c543bd05 309 for(Int_t i = 0; i < 64; i++) {
310 fV0ampl->Fill(esdV0->GetMultiplicity(i));
311 }
684b93f0 312 // Post output data.
313 PostData(1, fListOfHistos);
314}
315
316void AliAnaFwdDetsQA::Terminate(Option_t *)
317{
fc6c8a4a 318 // Terminate
319 // Store the output histos
684b93f0 320 fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
321 if (!fListOfHistos) {
322 Printf("ERROR: fListOfHistos not available");
323 return;
324 }
325
326 fT0vtxRec = dynamic_cast<TH1F*>(fListOfHistos->At(0));
327 fT0time = dynamic_cast<TH1F*>(fListOfHistos->At(1));
328 fT0mult = dynamic_cast<TH1F*>(fListOfHistos->At(2));
329 fT0vtxRecGen = dynamic_cast<TH2F*>(fListOfHistos->At(3));
330 fT0vtxRes = dynamic_cast<TH1F*>(fListOfHistos->At(4));
331
332 fV0a = dynamic_cast<TH1F*>(fListOfHistos->At(5));
333 fV0c = dynamic_cast<TH1F*>(fListOfHistos->At(6));
334 fV0multA = dynamic_cast<TH1F*>(fListOfHistos->At(7));
335 fV0multC = dynamic_cast<TH1F*>(fListOfHistos->At(8));
336 fV0multAcorr = dynamic_cast<TH2F*>(fListOfHistos->At(9));
337 fV0multCcorr = dynamic_cast<TH2F*>(fListOfHistos->At(10));
338 fV0Acorr = dynamic_cast<TH2F*>(fListOfHistos->At(11));
339 fV0Ccorr = dynamic_cast<TH2F*>(fListOfHistos->At(12));
340
c543bd05 341 fT0ampl = dynamic_cast<TH1F*>(fListOfHistos->At(13));
342 fV0ampl = dynamic_cast<TH1F*>(fListOfHistos->At(14));
343 fT0time2 = dynamic_cast<TH1F*>(fListOfHistos->At(15));
344
345
684b93f0 346 // draw the histograms if not in batch mode
347 if (!gROOT->IsBatch()) {
348 new TCanvas;
349 fT0vtxRec->DrawCopy("E");
350 new TCanvas;
351 fT0time->DrawCopy("E");
352 new TCanvas;
353 fT0mult->DrawCopy("E");
354 new TCanvas;
355 fT0vtxRes->DrawCopy("E");
356 new TCanvas;
357 fT0vtxRecGen->DrawCopy();
358 new TCanvas;
359 fV0a->DrawCopy("E");
360 new TCanvas;
361 fV0c->DrawCopy("E");
362 new TCanvas;
363 fV0multA->DrawCopy("E");
364 new TCanvas;
365 fV0multC->DrawCopy("E");
366 new TCanvas;
367 fV0multAcorr->DrawCopy();
368 new TCanvas;
369 fV0multCcorr->DrawCopy();
370 new TCanvas;
371 fV0Acorr->DrawCopy();
372 new TCanvas;
373 fV0Ccorr->DrawCopy();
c543bd05 374 new TCanvas;
375 fT0ampl->DrawCopy("E");
376 new TCanvas;
377 fV0ampl->DrawCopy("E");
378 new TCanvas;
379 fT0time2->DrawCopy("E");
684b93f0 380 }
381
382 // write the output histograms to a file
383 TFile* outputFile = TFile::Open("FwdDetsQA.root", "recreate");
384 if (!outputFile || !outputFile->IsOpen())
385 {
386 Error("AliAnaFwdDetsQA", "opening output file FwdDetsQA.root failed");
387 return;
388 }
389 fT0vtxRec->Write();
390 fT0time->Write();
391 fT0mult->Write();
392 fT0vtxRecGen->Write();
393 fT0vtxRes->Write();
c543bd05 394 fT0ampl->Write();
395 fT0time2->Write();
684b93f0 396
397 fV0a->Write();
398 fV0c->Write();
399 fV0multA->Write();
400 fV0multC->Write();
401 fV0multAcorr->Write();
402 fV0multCcorr->Write();
403 fV0Acorr->Write();
404 fV0Ccorr->Write();
c543bd05 405 fV0ampl->Write();
684b93f0 406
407 outputFile->Close();
408 delete outputFile;
409
410 //delete esd;
411 Info("AliAnaFwdDetsQA", "Successfully finished");
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
469
470
471
472
473
474