1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //------------------------------
17 // Analysis task for quality-assurance
18 // of forward detectors ESD
20 // 12/06/2009 cvetan.cheshkov@cern.ch
21 //------------------------------
31 #include "TParticle.h"
32 #include "AliVParticle.h"
33 #include "AliMCParticle.h"
34 #include "AliESDEvent.h"
36 #include "AliESDcascade.h"
37 #include "AliESDMuonTrack.h"
38 #include "AliESDCaloCluster.h"
40 #include "AliMCEvent.h"
42 #include "AliGenEventHeader.h"
44 #include "AliESDVZERO.h"
46 #include "AliAnaFwdDetsQA.h"
48 ClassImp(AliAnaFwdDetsQA)
50 AliAnaFwdDetsQA::AliAnaFwdDetsQA():
51 AliAnalysisTaskSE("AliAnaFwdDetsQA"),
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());
78 AliAnaFwdDetsQA::AliAnaFwdDetsQA(const char* name):
79 AliAnalysisTaskSE(name),
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());
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)
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);
121 TH1F *AliAnaFwdDetsQA::CreateEffHisto(const TH1F* hGen, const TH1F* hRec)
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");
128 hEff->SetStats(kFALSE);
129 hEff->SetMinimum(0.);
130 hEff->SetMaximum(110.);
131 hEff->GetYaxis()->SetTitle("#epsilon [%]");
133 for (Int_t iBin = 0; iBin <= nBins; iBin++) {
134 Double_t nGenEff = hGen->GetBinContent(iBin);
135 Double_t nRecEff = hRec->GetBinContent(iBin);
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);
144 hEff->SetBinContent(iBin, -100.);
145 hEff->SetBinError(iBin, 0);
152 Bool_t AliAnaFwdDetsQA::FitHisto(TH1* histo, Double_t& res, Double_t& resError)
156 static TF1* fitFunc = new TF1("fitFunc", "gaus");
157 fitFunc->SetLineWidth(2);
158 fitFunc->SetFillStyle(0);
159 Double_t maxFitRange = 2;
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));
175 void AliAnaFwdDetsQA::UserCreateOutputObjects()
178 // Create output container
179 AliInfo("AliAnaFwdDetsQA::UserCreateOutputObjects");
180 fListOfHistos = new TList();
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);
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);
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);
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");
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");
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);
233 void AliAnaFwdDetsQA::UserExec(Option_t */*option*/)
237 AliMCEvent* mcEvent = MCEvent();
239 Printf("ERROR: Could not retrieve MC event");
243 //Primary vertex needed
245 mcEvent->GenEventHeader()->PrimaryVertex(mcvtx);
247 AliStack *stack = mcEvent->Stack();
249 Printf("ERROR: Could not retrieve MC stack");
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) {
265 if (eta > -3.7 && eta < -1.7) {
271 AliVEvent* event = InputEvent();
273 Printf("ERROR: Could not retrieve event");
277 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
278 const AliESDTZERO* esdT0 = esd->GetESDTZERO();
279 Double_t t0zvtx = esdT0->GetT0zVertex();
280 Double_t t0time = esdT0->GetT0();
282 fT0vtxRec->Fill(t0zvtx);
283 fT0time->Fill(t0time);
284 fT0vtxRecGen->Fill(mcvtx[2],t0zvtx);
286 fT0vtxRes->Fill(t0zvtx - mcvtx[2]);
289 for(Int_t i = 0; i < 24; i++) {
290 t0mult += esdT0->GetT0amplitude()[i];
291 fT0ampl->Fill(esdT0->GetT0amplitude()[i]);
294 fT0mult->Fill(t0mult);
296 fT0time2->Fill(t0time);
298 AliESDVZERO* esdV0 = esd->GetVZEROData();
299 fV0a->Fill(esdV0->GetNbPMV0A());
300 fV0c->Fill(esdV0->GetNbPMV0C());
301 fV0multA->Fill(esdV0->GetMTotV0A());
302 fV0multC->Fill(esdV0->GetMTotV0C());
304 fV0multAcorr->Fill((Float_t)nV0A,esdV0->GetMTotV0A());
305 fV0multCcorr->Fill((Float_t)nV0C,esdV0->GetMTotV0C());
306 fV0Acorr->Fill((Float_t)nV0A,(Float_t)esdV0->GetNbPMV0A());
307 fV0Ccorr->Fill((Float_t)nV0C,(Float_t)esdV0->GetNbPMV0C());
309 for(Int_t i = 0; i < 64; i++) {
310 fV0ampl->Fill(esdV0->GetMultiplicity(i));
313 PostData(1, fListOfHistos);
316 void AliAnaFwdDetsQA::Terminate(Option_t *)
319 // Store the output histos
320 fListOfHistos = dynamic_cast<TList*>(GetOutputData(1));
321 if (!fListOfHistos) {
322 Printf("ERROR: fListOfHistos not available");
326 fT0vtxRec = dynamic_cast<TH1F*>(fListOfHistos->At(0));
327 fT0time = dynamic_cast<TH1F*>(fListOfHistos->At(1));
328 fT0mult = dynamic_cast<TH1F*>(fListOfHistos->At(2));
329 fT0vtxRecGen = dynamic_cast<TH2F*>(fListOfHistos->At(3));
330 fT0vtxRes = dynamic_cast<TH1F*>(fListOfHistos->At(4));
332 fV0a = dynamic_cast<TH1F*>(fListOfHistos->At(5));
333 fV0c = dynamic_cast<TH1F*>(fListOfHistos->At(6));
334 fV0multA = dynamic_cast<TH1F*>(fListOfHistos->At(7));
335 fV0multC = dynamic_cast<TH1F*>(fListOfHistos->At(8));
336 fV0multAcorr = dynamic_cast<TH2F*>(fListOfHistos->At(9));
337 fV0multCcorr = dynamic_cast<TH2F*>(fListOfHistos->At(10));
338 fV0Acorr = dynamic_cast<TH2F*>(fListOfHistos->At(11));
339 fV0Ccorr = dynamic_cast<TH2F*>(fListOfHistos->At(12));
341 fT0ampl = dynamic_cast<TH1F*>(fListOfHistos->At(13));
342 fV0ampl = dynamic_cast<TH1F*>(fListOfHistos->At(14));
343 fT0time2 = dynamic_cast<TH1F*>(fListOfHistos->At(15));
346 // draw the histograms if not in batch mode
347 if (!gROOT->IsBatch()) {
349 fT0vtxRec->DrawCopy("E");
351 fT0time->DrawCopy("E");
353 fT0mult->DrawCopy("E");
355 fT0vtxRes->DrawCopy("E");
357 fT0vtxRecGen->DrawCopy();
363 fV0multA->DrawCopy("E");
365 fV0multC->DrawCopy("E");
367 fV0multAcorr->DrawCopy();
369 fV0multCcorr->DrawCopy();
371 fV0Acorr->DrawCopy();
373 fV0Ccorr->DrawCopy();
375 fT0ampl->DrawCopy("E");
377 fV0ampl->DrawCopy("E");
379 fT0time2->DrawCopy("E");
382 // write the output histograms to a file
383 TFile* outputFile = TFile::Open("FwdDetsQA.root", "recreate");
384 if (!outputFile || !outputFile->IsOpen())
386 Error("AliAnaFwdDetsQA", "opening output file FwdDetsQA.root failed");
392 fT0vtxRecGen->Write();
401 fV0multAcorr->Write();
402 fV0multCcorr->Write();
411 Info("AliAnaFwdDetsQA", "Successfully finished");