Stupid bug fix in new superlight mode (from Zurich airport)
[u/mrichter/AliRoot.git] / PWGPP / 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);
ae02804c 278 if (!esd) {
279 Printf("ERROR: Could not retrieve esd");
280 return;
281 }
684b93f0 282 const AliESDTZERO* esdT0 = esd->GetESDTZERO();
283 Double_t t0zvtx = esdT0->GetT0zVertex();
284 Double_t t0time = esdT0->GetT0();
285
286 fT0vtxRec->Fill(t0zvtx);
287 fT0time->Fill(t0time);
288 fT0vtxRecGen->Fill(mcvtx[2],t0zvtx);
289 t0zvtx *= -1.0;
290 fT0vtxRes->Fill(t0zvtx - mcvtx[2]);
291
292 Double_t t0mult = 0;
293 for(Int_t i = 0; i < 24; i++) {
294 t0mult += esdT0->GetT0amplitude()[i];
c543bd05 295 fT0ampl->Fill(esdT0->GetT0amplitude()[i]);
684b93f0 296 }
297
298 fT0mult->Fill(t0mult);
c543bd05 299 if (t0mult > 10)
300 fT0time2->Fill(t0time);
684b93f0 301
302 AliESDVZERO* esdV0 = esd->GetVZEROData();
303 fV0a->Fill(esdV0->GetNbPMV0A());
304 fV0c->Fill(esdV0->GetNbPMV0C());
305 fV0multA->Fill(esdV0->GetMTotV0A());
306 fV0multC->Fill(esdV0->GetMTotV0C());
307
308 fV0multAcorr->Fill((Float_t)nV0A,esdV0->GetMTotV0A());
309 fV0multCcorr->Fill((Float_t)nV0C,esdV0->GetMTotV0C());
310 fV0Acorr->Fill((Float_t)nV0A,(Float_t)esdV0->GetNbPMV0A());
311 fV0Ccorr->Fill((Float_t)nV0C,(Float_t)esdV0->GetNbPMV0C());
312
c543bd05 313 for(Int_t i = 0; i < 64; i++) {
314 fV0ampl->Fill(esdV0->GetMultiplicity(i));
315 }
684b93f0 316 // Post output data.
317 PostData(1, fListOfHistos);
318}
319
320void AliAnaFwdDetsQA::Terminate(Option_t *)
321{
fc6c8a4a 322 // Terminate
323 // Store the output histos
684b93f0 324 fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
325 if (!fListOfHistos) {
326 Printf("ERROR: fListOfHistos not available");
327 return;
328 }
329
684b93f0 330 //delete esd;
331 Info("AliAnaFwdDetsQA", "Successfully finished");
332}
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394