AliPhysics  8d00e07 (8d00e07)
AliAnaFwdDetsQA.cxx
Go to the documentation of this file.
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 "AliVParticle.h"
33 #include "AliMCParticle.h"
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 
52  fListOfHistos(0),
53  fT0vtxRec(0),
54  fT0vtxRecGen(0),
55  fT0time(0),
56  fT0time2(0),
57  fT0mult(0),
58  fT0vtxRes(0),
59  fT0ampl(0),
60  fV0a(0),
61  fV0c(0),
62  fV0multA(0),
63  fV0multC(0),
64  fV0multAcorr(0),
65  fV0multCcorr(0),
66  fV0Acorr(0),
67  fV0Ccorr(0),
68  fV0ampl(0)
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 
79 AliAnalysisTaskSE(name),
80  fListOfHistos(0),
81  fT0vtxRec(0),
82  fT0vtxRecGen(0),
83  fT0time(0),
84  fT0time2(0),
85  fT0mult(0),
86  fT0vtxRes(0),
87  fT0ampl(0),
88  fV0a(0),
89  fV0c(0),
90  fV0multA(0),
91  fV0multC(0),
92  fV0multAcorr(0),
93  fV0multCcorr(0),
94  fV0Acorr(0),
95  fV0Ccorr(0),
96  fV0ampl(0)
97 {
98  // Constructor
99  // Define input and output slots here
100  // Input slot #0 works with a TChain
101  AliInfo("Constructor AliAnaFwdDetsQA");
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  // helper method which can be used
112  // in order to create a histogram
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 
121 TH1F *AliAnaFwdDetsQA::CreateEffHisto(const TH1F* hGen, const TH1F* hRec)
122 {
123  // helper method which can be used
124  // in order create an efficiency histogram
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);
140  if (error < 1e-12) error = 0.0001;
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 
153 {
154  // fit a gaussian to
155  // a histogram
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 
176 {
177  // Create histograms
178  // Create output container
179  AliInfo("AliAnaFwdDetsQA::UserCreateOutputObjects");
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]", "");
184  fT0time2 = CreateHisto("hT0time2", "Time0 reconstruction with T0 (mult > 10)", 5000, 10000, 20000, "t_{0} [ps]", "");
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]", "");
187  fT0ampl = CreateHisto("hT0ampl","T0 multiplicity in single channel (all T0 channels)",400,-0.5,99.5);
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.);
199  fV0ampl = CreateHisto("hV0ampl","V0 multiplicity in single channel (all V0 channels)",400,-0.5,99.5);
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);
219  fListOfHistos->Add(fT0vtxRes);
220  fListOfHistos->Add(fV0a);
221  fListOfHistos->Add(fV0c);
222  fListOfHistos->Add(fV0multA);
223  fListOfHistos->Add(fV0multC);
226  fListOfHistos->Add(fV0Acorr);
227  fListOfHistos->Add(fV0Ccorr);
228  fListOfHistos->Add(fT0ampl);
229  fListOfHistos->Add(fV0ampl);
230  fListOfHistos->Add(fT0time2);
231 }
232 
234 {
235  // The analysis code
236  // goes here
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;
257  AliMCParticle* track = (AliMCParticle*)mcEvent->GetTrack(iTracks);
258  TParticle* particle = track->Particle();
259  if (!particle) continue;
260  if (track->Charge() == 0) continue;
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  if (!esd) {
279  Printf("ERROR: Could not retrieve esd");
280  return;
281  }
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];
295  fT0ampl->Fill(esdT0->GetT0amplitude()[i]);
296  }
297 
298  fT0mult->Fill(t0mult);
299  if (t0mult > 10)
300  fT0time2->Fill(t0time);
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 
313  for(Int_t i = 0; i < 64; i++) {
314  fV0ampl->Fill(esdV0->GetMultiplicity(i));
315  }
316  // Post output data.
317  PostData(1, fListOfHistos);
318 }
319 
321 {
322  // Terminate
323  // Store the output histos
324  fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
325  if (!fListOfHistos) {
326  Printf("ERROR: fListOfHistos not available");
327  return;
328  }
329 
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 
double Double_t
Definition: External.C:58
Definition: External.C:236
const char * title
Definition: MakeQAPdf.C:27
TH1F * CreateHisto(const char *name, const char *title, Int_t nBins, Double_t xMin, Double_t xMax, const char *xLabel=NULL, const char *yLabel=NULL)
TH1F * CreateEffHisto(const TH1F *hGen, const TH1F *hRec)
AliStack * stack
virtual void UserExec(Option_t *option)
virtual void UserCreateOutputObjects()
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
virtual void Terminate(Option_t *)
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
Bool_t FitHisto(TH1 *histo, Double_t &res, Double_t &resError)
Definition: External.C:196