Variable shadowing removed (B. Hippolyte, AM)
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliAnalysisElectronTask.cxx
CommitLineData
809a4336 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 * The analysis task:
17 * Filling an AliCFContainer with the quantities pt, eta and phi
18 * for tracks which survivied the particle cuts (MC resp. ESD tracks)
19 * Track selection is done using the AliHFE package
20 *
21 * Author:
22 * Markus Fasel <M.Fasel@gsi.de>
23 */
24#include <TAxis.h>
25#include <TCanvas.h>
26#include <TChain.h>
27#include <TH1F.h>
28#include <TH1I.h>
29#include <TH2F.h>
30#include <TIterator.h>
31#include <TList.h>
32#include <TLegend.h>
33#include <TMath.h>
34#include <TObjArray.h>
35#include <TParticle.h>
36#include <TProfile.h>
37#include <TTree.h>
38
39#include "AliCFContainer.h"
40#include "AliCFManager.h"
41#include "AliESDEvent.h"
42#include "AliESDInputHandler.h"
43#include "AliESDtrack.h"
44#include "AliESDtrackCuts.h"
45#include "AliLog.h"
46#include "AliAnalysisManager.h"
47#include "AliMCEvent.h"
48#include "AliMCEventHandler.h"
49#include "AliMCParticle.h"
50#include "AliPID.h"
51
52#include "AliHFEpid.h"
53#include "AliHFEcuts.h"
54#include "AliAnalysisElectronTask.h"
55
56//____________________________________________________________
57AliAnalysisElectronTask::AliAnalysisElectronTask():
58 AliAnalysisTask("PID efficiency Analysis", "")
59 , fESD(0x0)
60 , fMC(0x0)
61 , fCFM(0x0)
62 , fPID(0x0)
63 , fCuts(0x0)
64 , fNEvents(0x0)
65 , fQA(0x0)
66{
67 DefineInput(0, TChain::Class());
68 DefineOutput(0, TH1I::Class());
69 DefineOutput(1, AliCFContainer::Class());
70 DefineOutput(2, TList::Class());
71
72 // Initialize cuts
73 fCuts = new AliHFEcuts;
74 fPID = new AliHFEpid;
75}
76
77//____________________________________________________________
78AliAnalysisElectronTask::~AliAnalysisElectronTask(){
79 if(fESD) delete fESD;
80 if(fMC) delete fMC;
81 if(fPID) delete fPID;
82 if(fQA) delete fQA;
83 if(fCuts) delete fCuts;
84 if(fNEvents) delete fNEvents;
85 fQA = 0x0; fMC = 0x0; fESD = 0x0; fCFM = 0x0; fPID = 0x0; fCuts = 0x0;
86}
87
88//____________________________________________________________
89void AliAnalysisElectronTask::ConnectInputData(Option_t *){
90 TTree *esdchain = dynamic_cast<TChain *>(GetInputData(0));
91 if(!esdchain){
92 AliError("ESD chain empty");
93 return;
94 } else {
95 esdchain->SetBranchStatus("Tracks", 1);
96 }
97 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
98 if(!esdH){
99 AliError("No ESD input handler");
100 return;
101 } else {
102 fESD = esdH->GetEvent();
103 }
104 AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
105 if(!mcH){
106 AliError("No MC truth handler");
107 return;
108 } else {
109 fMC = mcH->MCEvent();
110 }
111}
112
113//____________________________________________________________
114void AliAnalysisElectronTask::CreateOutputObjects(){
115 fNEvents = new TH1I("nEvents", "Number of Events in the Analysis", 2, 0, 2); // Number of Events neccessary for the analysis and not a QA histogram
116 // First Step: TRD alone
117 if(!fQA) fQA = new TList;
118 fQA->AddAt(new TProfile("conr", "Electron PID contamination", 20, 0, 20), 0);
119 fQA->AddAt(new TH1F("alpha_rec", "Alpha from reconstructed tracks with TRD hits", 36, -TMath::Pi(), TMath::Pi()), 1);
120 fQA->AddAt(new TH1F("alpha_sim", "Alpha from simulated electron tracks", 36, -TMath::Pi(), TMath::Pi()), 2);
121 fQA->AddAt(new TH1F("nElectron", "Number of electrons", 100, 0, 100), 3);
122 fQA->AddAt(new TProfile("pidquality", "TRD PID quality as function of momentum", 20, 0, 20), 4);
123 fQA->AddAt(new TProfile("ntrdclusters", "Number of TRD clusters as function of momentum", 20, 0, 20), 5);
124 fQA->AddAt(new TH1F("chi2TRD","#chi2 per TRD cluster", 20, 0, 20), 6);
125
126 // Initialize correction Framework and Cuts
127 fCFM = new AliCFManager;
128 MakeParticleContainer();
129 // Temporary fix: Initialize particle cuts with 0x0
130 for(Int_t istep = 0; istep < fCFM->GetParticleContainer()->GetNStep(); istep++)
131 fCFM->SetParticleCutsList(istep, 0x0);
132 if(IsQAOn()){
133 fCuts->SetDebugMode();
134 fQA->AddAt(fCuts->GetQAhistograms(), 7);
135 }
136 fCuts->CreateStandardCuts();
137 fCuts->Initialize(fCFM);
138
139 // Initialize PID
140 if(IsQAOn()){
141 fPID->SetQAOn();
142 fQA->AddAt(fPID->GetQAhistograms(), 8);
143 }
144 fPID->SetHasMCData(kTRUE);
145 fPID->InitializePID("TRD");
146}
147
148//____________________________________________________________
149void AliAnalysisElectronTask::Exec(Option_t *){
150 //
151 // Run the analysis
152 //
153 if(!fESD){
154 AliError("No ESD Event");
155 return;
156 }
157 if(!fMC){
158 AliError("No MC Event");
159 return;
160 }
161 fCFM->SetEventInfo(fMC);
162 fPID->SetMCEvent(fMC);
163
164 //fCFM->CheckEventCuts(AliCFManager::kEvtGenCuts, fMC);
165
166 Double_t pt = 0;
167 Double_t container[3];
168
169 // Loop over the Monte Carlo tracks to see whether we have overlooked any track
170 AliMCParticle *mctrack = 0x0;
171 Int_t nElectrons = 0;
172 for(Int_t imc = fMC->GetNumberOfTracks(); imc--;){
173 mctrack = fMC->GetTrack(imc);
174 container[0] = mctrack->Pt();
175 container[1] = mctrack->Eta();
176 container[2] = mctrack->Phi();
177
178 if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) continue;
179 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCGenerated);
180 (dynamic_cast<TH1F *>(fQA->At(2)))->Fill(mctrack->Phi() - TMath::Pi());
181 if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, mctrack)) continue;
182 // find the label in the vector
183 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCInAcceptance);
184 nElectrons++;
185 }
186 (dynamic_cast<TH1F *>(fQA->At(3)))->Fill(nElectrons);
187
188 // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
189 AliESDtrack *track = 0x0;
190 for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
191 track = fESD->GetTrack(itrack);
192 container[0] = track->Pt();
193 container[1] = track->Eta();
194 container[2] = track->Phi();
195 if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKine, track)) continue;
196 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecKine);
197
198 // Check TRD criterions (outside the correction framework)
199 if(track->GetTRDncls()){
200 (dynamic_cast<TH1F *>(fQA->At(6)))->Fill(track->GetTRDchi2()/track->GetTRDncls());
201 (dynamic_cast<TH1F *>(fQA->At(1)))->Fill(track->GetAlpha()); // Check the acceptance without tight cuts
202 (dynamic_cast<TProfile *>(fQA->At(4)))->Fill(container[0], track->GetTRDpidQuality());
203 (dynamic_cast<TProfile *>(fQA->At(5)))->Fill(container[0], track->GetTRDncls());
204 }
205 if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) continue;
206 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecPrim);
207 if(fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcuts, track)) continue;
208 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcuts);
209 // track accepted, do PID
210 if(!fPID->IsSelected(track)) continue;
211 fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcuts + 1);
212 }
213 fNEvents->Fill(1);
214
215 // Done!!!
216 PostData(0, fNEvents);
217 PostData(1, fCFM->GetParticleContainer());
218 PostData(2, fQA);
219}
220
221//____________________________________________________________
222void AliAnalysisElectronTask::Terminate(Option_t *){
223 //
224 // Terminate not implemented at the moment
225 //
226}
227
228//____________________________________________________________
229void AliAnalysisElectronTask::MakeParticleContainer(){
230 //
231 // Create the particle container for the correction framework manager and
232 // link it
233 //
234 const Int_t nvar = 3 ; //number of variables on the grid:pt,eta, phi
235 const Double_t ptmin = 0., ptmax = 10.;
236 const Double_t etamin = -0.9, etamax = 0.9;
237 const Double_t phimin = 0., phimax = 2. * TMath::Pi();
238
239
240 //arrays for the number of bins in each dimension
241 Int_t iBin[nvar];
242 iBin[0] = 20; //bins in pt
243 iBin[1] = 8; //bins in eta
244 iBin[2] = 18; // bins in phi
245
246 //arrays for lower bounds :
247 Double_t *binLim1 = new Double_t[iBin[0] + 1];
248 Double_t *binLim2 = new Double_t[iBin[1] + 1];
249 Double_t *binLim3 = new Double_t[iBin[2] + 1];
250
251 //values for bin lower bounds
252 for(Int_t i=0; i<=iBin[0]; i++) binLim1[i]=(Double_t)ptmin + (ptmax-ptmin)/iBin[0]*(Double_t)i;
253 for(Int_t i=0; i<=iBin[1]; i++) binLim2[i]=(Double_t)etamin + (etamax-etamin)/iBin[1]*(Double_t)i;
254 for(Int_t i=0; i<=iBin[2]; i++) binLim3[i]=(Double_t)phimin + (phimax-phimin)/iBin[2]*(Double_t)i;
255
256 //one "container" for MC
257 AliCFContainer* container = new AliCFContainer("container","container for tracks", AliHFEcuts::kNcutSteps + 1, nvar, iBin);
258 //setting the bin limits
259 container -> SetBinLimits(0,binLim1);
260 container -> SetBinLimits(1,binLim2);
261 container -> SetBinLimits(2,binLim3);
262 fCFM->SetParticleContainer(container);
263}