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