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