Update master to aliroot
[u/mrichter/AliRoot.git] / ANALYSIS / AliPriorsTask.cxx
CommitLineData
9ab172e7 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// Example of task (running locally) to calculate
18// the a priori concentration within the new class
19// AliEventPoolLoop
20//-----------------------------------------------------------------------
21
22
23#include <TParticle.h>
24#include <TLine.h>
25#include <TLegend.h>
26#include <TH1I.h>
27#include <TChain.h>
28#include <TH2D.h>
29
30#include "AliPriorsTask.h"
31#include "AliStack.h"
32#include "AliMCEvent.h"
33#include "AliAnalysisManager.h"
34#include "AliESDEvent.h"
35#include "AliESDtrack.h"
36#include "AliLog.h"
37
38ClassImp(AliPriorsTask)
39
40//__________________________________________________________________________
41AliPriorsTask::AliPriorsTask() :
42 fHistEventsProcessed(0x0),
43 fCalcPriors(0x0),
44 fNiterations(0x0)
45{
46 //
47 //Default ctor
48 //
49 for(Int_t i=0; i< AliPID::kSPECIES; i++){
50 fPriors[i]=0;
51 fRecId[i]=0;
52 fMCId[i]=0;
53 }
54}
55//___________________________________________________________________________
56AliPriorsTask::AliPriorsTask(const Char_t* name) :
57 AliAnalysisTaskSE(name),
58 fHistEventsProcessed(0x0),
59 fCalcPriors(0x0),
60 fNiterations(0)
61{
62 //
63 // Constructor. Initialization of Inputs and Outputs
64 //
65
66 for(Int_t i=0; i< AliPID::kSPECIES; i++){
67 fPriors[i]=0;
68 fRecId[i]=0;
69 fMCId[i]=0;
70 }
71
72
73 Info("AliPriorsTask","Calling Constructor");
74
75 /*
76 DefineInput(0) and DefineOutput(0)
77 are taken care of by AliAnalysisTaskSE constructor
78 */
79 DefineOutput(1,TH1I::Class());
80 DefineOutput(2,TH2D::Class());
81}
82
83//___________________________________________________________________________
84AliPriorsTask& AliPriorsTask::operator=(const AliPriorsTask& c)
85{
86 //
87 // Assignment operator
88 //
89 if (this!=&c) {
90 AliAnalysisTaskSE::operator=(c) ;
91 fHistEventsProcessed = c.fHistEventsProcessed;
92 fCalcPriors = c.fCalcPriors;
93 for(Int_t i=0; i< AliPID::kSPECIES; i++){
94 fPriors[i]=c.fPriors[i];
95 fRecId[i]=c.fRecId[i];
96 fMCId[i]=c.fMCId[i];
97 }
98 fNiterations=c.fNiterations;
99 }
100 return *this;
101}
102
103//___________________________________________________________________________
104AliPriorsTask::AliPriorsTask(const AliPriorsTask& c) :
105 AliAnalysisTaskSE(c),
106 fHistEventsProcessed(c.fHistEventsProcessed),
107 fCalcPriors(c.fCalcPriors),
108 fNiterations(c.fNiterations)
109{
110 //
111 // Copy Constructor
112 //
113 for(Int_t i=0; i< AliPID::kSPECIES; i++){
114 fPriors[i]=c.fPriors[i];
115 fRecId[i]=c.fRecId[i];
116 fMCId[i]=c.fMCId[i];
117 }
118}
119
120//___________________________________________________________________________
121AliPriorsTask::~AliPriorsTask() {
122 //
123 // Destructor
124 //
125 if (fHistEventsProcessed) delete fHistEventsProcessed ;
126 if (fCalcPriors) delete fCalcPriors ;
127}
128
129//_________________________________________________
130void AliPriorsTask::UserExec(Option_t *)
131{
132 //
133 // Main loop function
134 //
135 AliDebug(1,"UserExec") ;
136
137 if(!fInputEvent) {
138 Printf("No input event!");
139 return;
140 }
141
142 AliESDEvent* fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
143 if (!fESD) {
144 Error("UserExec","NO ESD FOUND!");
145 return;
146 }
147
148 if (!fMCEvent) {
149 Error("UserExec","NO MC INFO FOUND");
150 return;
151 }
152
153 //loop on tracks
154 for (Int_t iTrack = 0; iTrack<fESD->GetNumberOfTracks(); iTrack++) {
155
156 AliESDtrack* track = fESD->GetTrack(iTrack);
157 if (!track) continue;
158
159 //some quality/pid requirements
160 if (!track->IsOn(AliESDtrack::kITSrefit)) continue;
161 if (!track->IsOn(AliESDtrack::kESDpid)) continue;
162 if (!track->IsOn(AliESDtrack::kTPCpid)) continue;
163 if (!track->IsOn(AliESDtrack::kTOFpid)) continue;
164
165 if(track->GetLabel()<0) continue;
166
167 Double_t esdpid[AliPID::kSPECIES]; track->GetESDpid(esdpid);
168
169 AliPID pid(esdpid,kTRUE);
170 pid.SetPriors(fPriors);
171 AliPID::EParticleType id = pid.GetMostProbable();
172
173 if ((Int_t)id < (Int_t)AliPID::kSPECIES) fRecId[id]++;
174
175 if(fNiterations == 0 ){
b866961d 176 Int_t pdgcode = TMath::Abs(((fMCEvent->GetTrack(track->GetLabel()))->Particle())->GetPdgCode());
9ab172e7 177 switch(pdgcode){
178 case 11 : fMCId[AliPID::kElectron]++; break;
179 case 13 : fMCId[AliPID::kMuon]++; break;
180 case 211 : fMCId[AliPID::kPion]++; break;
181 case 321 : fMCId[AliPID::kKaon]++; break;
182 case 2212: fMCId[AliPID::kProton]++; break;
183 default : break;
184 }
185 }
186 }
187
188
189 fHistEventsProcessed->Fill(0);
190
191 PostData(1,fHistEventsProcessed) ;
192 PostData(2,fCalcPriors) ;
193}
194
195//___________________________________________________________________________
196void AliPriorsTask::Terminate(Option_t*)
197{
198 // The Terminate() function is the last function to be called during
199 // a query. It always runs on the client, it can be used to present
200 // the results graphically or save the results to file.
201
202 AliAnalysisTaskSE::Terminate();
203
204 Double_t MCsum =0; for(Int_t i=0; i< AliPID::kSPECIES; i++ ) MCsum+=fMCId[i];
205 if(MCsum>0) Printf(" MC Priors : %f %f %f %f %f",fMCId[0]/MCsum,fMCId[1]/MCsum,fMCId[2]/MCsum,fMCId[3]/MCsum,fMCId[4]/MCsum);
206 Printf(" Priors : %f %f %f %f %f \n",fPriors[0],fPriors[1],fPriors[2],fPriors[3],fPriors[4]);
207
208 TH2D *histo= dynamic_cast<TH2D*> (GetOutputData(2));
209 for(Int_t ip = 0; ip < AliPID::kSPECIES; ip++) histo->Fill(fNiterMax-1,fPriors[ip]);
210 histo->SetMarkerStyle(8);
211 histo->DrawCopy();
212
213 Int_t ltype[AliPID::kSPECIES] = {1,2,3,2,1};
214 Int_t lcolor[AliPID::kSPECIES] = {1,2,8,38,46};
215
216 TLine *l[AliPID::kSPECIES];
217 for(Int_t isp =0; isp < AliPID::kSPECIES; isp++) {
218 if(MCsum>0)l[isp]= new TLine(0.,fMCId[isp]/MCsum,fNiterMax,fMCId[isp]/MCsum);
219 l[isp]->SetLineWidth(isp+1);
220 l[isp]->SetLineColor(lcolor[isp]);
221 l[isp]->SetLineStyle(ltype[isp]);
222 l[isp]->Draw("same");
223 }
224
225 const char *part[AliPID::kSPECIES]={"e","#mu","#pi","K","p"};
226 TLegend *mcpriors = new TLegend(0.4,0.6,0.6,0.8);
227 for(Int_t ip=0; ip < AliPID::kSPECIES; ip++) mcpriors->AddEntry(l[ip],part[ip],"l");
228 mcpriors->Draw("same");
229
230}
231//___________________________________________________________________________
232void AliPriorsTask::UserCreateOutputObjects()
233{
234 //
235 //HERE ONE CAN CREATE OUTPUT OBJECTS
236 //BEFORE THE EXECUTION OF THE TASK
237 //
238 OpenFile(1);
239 fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
240
241 OpenFile(2);
242 fCalcPriors = new TH2D("priors",Form("priors in %i iterations",fNiterMax),fNiterMax*10,-0.5,fNiterMax+0.5,100,0,1);
243 fCalcPriors->SetXTitle("Iteration step");
244 fCalcPriors->SetYTitle("concentrations");
245}
246//________________________________________________________________________
247
248Bool_t AliPriorsTask::NotifyBinChange()
249{
250 //
251 //METHOD CALLED AT THE END OF THE EVENT LOOP.
252 //HERE CONCENTRATIONS ARE CALCULATED, STORED
253 //AND USED IN THE NEXT EVENT LOOP
254 //
255
256 Double_t sum =0;
257 for(Int_t i=0; i< AliPID::kSPECIES; i++) {
258 sum+=fRecId[i];
259 }
260 if(sum == 0) return kFALSE;
261
262 Double_t conc[AliPID::kSPECIES];
263 for(Int_t i=0; i< AliPID::kSPECIES; i++) {
264 conc[i]=fRecId[i]/sum;
265 fCalcPriors->Fill(fNiterations,conc[i]);
266 fRecId[i]=0;
267 }
268
269 SetPriors(conc);
270 fNiterations++;
271
272 return kTRUE;
273}