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