fix histogram title
[u/mrichter/AliRoot.git] / PWG1 / AliAnaFwdDetsQA.cxx
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
50 AliAnaFwdDetsQA::AliAnaFwdDetsQA():
51 AliAnalysisTaskSE("AliAnaFwdDetsQA"),
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
78 AliAnaFwdDetsQA::AliAnaFwdDetsQA(const char* name):
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
152 Bool_t AliAnaFwdDetsQA::FitHisto(TH1* histo, Double_t& res, Double_t& resError)
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
175 void AliAnaFwdDetsQA::UserCreateOutputObjects()
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);
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);
228   fListOfHistos->Add(fT0ampl);
229   fListOfHistos->Add(fV0ampl);
230   fListOfHistos->Add(fT0time2);
231 }
232
233 void AliAnaFwdDetsQA::UserExec(Option_t */*option*/)
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
320 void AliAnaFwdDetsQA::Terminate(Option_t *)
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