cbfd09c9d61b38cf259cef0c115b59fa51784b71
[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 "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
46 ClassImp(AliAnaFwdDetsQA)
47
48 AliAnaFwdDetsQA::AliAnaFwdDetsQA():
49 AliAnalysisTaskSE("AliAnaFwdDetsQA"),
50   fListOfHistos(0),
51   fT0vtxRec(0),
52   fT0vtxRecGen(0),
53   fT0time(0),
54   fT0time2(0),
55   fT0mult(0),
56   fT0vtxRes(0),
57   fT0ampl(0),
58   fV0a(0),
59   fV0c(0),
60   fV0multA(0),
61   fV0multC(0),
62   fV0multAcorr(0),
63   fV0multCcorr(0),
64   fV0Acorr(0),
65   fV0Ccorr(0),
66   fV0ampl(0)
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
76 AliAnaFwdDetsQA::AliAnaFwdDetsQA(const char* name):
77 AliAnalysisTaskSE(name),
78   fListOfHistos(0),
79   fT0vtxRec(0),
80   fT0vtxRecGen(0),
81   fT0time(0),
82   fT0time2(0),
83   fT0mult(0),
84   fT0vtxRes(0),
85   fT0ampl(0),
86   fV0a(0),
87   fV0c(0),
88   fV0multA(0),
89   fV0multC(0),
90   fV0multAcorr(0),
91   fV0multCcorr(0),
92   fV0Acorr(0),
93   fV0Ccorr(0),
94   fV0ampl(0)
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
105 TH1F * 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
118 TH1F *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
148 Bool_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
170 void 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]", "");
179   fT0time2  = CreateHisto("hT0time2", "Time0 reconstruction with T0 (mult > 10)", 5000, 10000, 20000, "t_{0} [ps]", "");
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]", "");
182   fT0ampl   = CreateHisto("hT0ampl","T0 multiplicity in single channel (all T0 channels)",400,-0.5,99.5);
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.);
194   fV0ampl  = CreateHisto("hV0ampl","V0 multiplicity in single channel (all V0 channels)",400,-0.5,99.5);
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);
223   fListOfHistos->Add(fT0ampl);
224   fListOfHistos->Add(fV0ampl);
225   fListOfHistos->Add(fT0time2);
226 }
227
228 void 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;
254     if (track->Charge() == 0) continue;
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];
285     fT0ampl->Fill(esdT0->GetT0amplitude()[i]);
286   }
287
288   fT0mult->Fill(t0mult);
289   if (t0mult > 10)
290     fT0time2->Fill(t0time);
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
303   for(Int_t i = 0; i < 64; i++) {
304     fV0ampl->Fill(esdV0->GetMultiplicity(i));
305   }
306   // Post output data.
307   PostData(1, fListOfHistos);
308 }
309
310 void 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
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
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();
366     new TCanvas;
367     fT0ampl->DrawCopy("E");
368     new TCanvas;
369     fV0ampl->DrawCopy("E");
370     new TCanvas;
371     fT0time2->DrawCopy("E");
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();
386   fT0ampl->Write();
387   fT0time2->Write();
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();
397   fV0ampl->Write();
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