Variable shadowing removed (B. Hippolyte, AM)
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliAnalysisElectronTask.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  * 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 //____________________________________________________________
57 AliAnalysisElectronTask::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 //____________________________________________________________
78 AliAnalysisElectronTask::~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 //____________________________________________________________
89 void 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 //____________________________________________________________
114 void 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 //____________________________________________________________
149 void 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 //____________________________________________________________
222 void AliAnalysisElectronTask::Terminate(Option_t *){
223         //
224         // Terminate not implemented at the moment
225         //
226 }
227
228 //____________________________________________________________
229 void 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 }