changed the order of call of endofcycle so that images are produced
[u/mrichter/AliRoot.git] / PWG2 / AliProtonCorrectionTask.cxx
CommitLineData
b5fc8a3e 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, on AliEn and CAF),
18// which provides standard way of calculating acceptance and efficiency
19// between different steps of the procedure.
20// The ouptut of the task is a AliCFContainer from which the efficiencies
21// can be calculated
22//-----------------------------------------------------------------------
23// Author : R. Vernet, Consorzio Cometa - Catania (it)
24//-----------------------------------------------------------------------
25
26
27#ifndef ALIPROTONCORRECTIONTASK_CXX
28#define ALIPROTONCORRECTIONTASK_CXX
29
30#include "AliProtonCorrectionTask.h"
31#include "TCanvas.h"
32#include "AliStack.h"
33#include "TParticle.h"
34#include "TH1I.h"
35#include "AliMCEvent.h"
36#include "AliAnalysisManager.h"
37#include "AliESDEvent.h"
38#include "AliAODEvent.h"
39#include "AliCFManager.h"
40#include "AliCFCutBase.h"
41#include "AliCFContainer.h"
42#include "TChain.h"
43#include "AliESDtrack.h"
44#include "AliLog.h"
45//#include "AliTPCtrack.h"
46
47ClassImp(AliProtonCorrectionTask)
48
49//__________________________________________________________________________
50AliProtonCorrectionTask::AliProtonCorrectionTask() :
51 fReadTPCTracks(0),
52 fReadAODData(0),
3079041e 53 fCFManagerProtons(0x0),
54 fCFManagerAntiProtons(0x0),
b5fc8a3e 55 fQAHistList(0x0),
56 fHistEventsProcessed(0x0)
57{
58 //
59 //Default ctor
60 //
61}
62//___________________________________________________________________________
63AliProtonCorrectionTask::AliProtonCorrectionTask(const Char_t* name) :
64 AliAnalysisTaskSE(name),
65 fReadTPCTracks(0),
66 fReadAODData(0),
3079041e 67 fCFManagerProtons(0x0),
68 fCFManagerAntiProtons(0x0),
b5fc8a3e 69 fQAHistList(0x0),
340705b3 70 fHistEventsProcessed(0x0) {
b5fc8a3e 71 //
72 // Constructor. Initialization of Inputs and Outputs
73 //
74 Info("AliProtonCorrectionTask","Calling Constructor");
75
76 /*
77 DefineInput(0) and DefineOutput(0)
78 are taken care of by AliAnalysisTaskSE constructor
79 */
80 DefineOutput(1,TH1I::Class());
81 DefineOutput(2,AliCFContainer::Class());
3079041e 82 DefineOutput(3,AliCFContainer::Class());
83 DefineOutput(4,TList::Class());
b5fc8a3e 84}
85
86//___________________________________________________________________________
340705b3 87AliProtonCorrectionTask& AliProtonCorrectionTask::operator=(const AliProtonCorrectionTask& c) {
b5fc8a3e 88 //
89 // Assignment operator
90 //
91 if (this!=&c) {
92 AliAnalysisTaskSE::operator=(c) ;
93 fReadTPCTracks = c.fReadTPCTracks ;
94 fReadAODData = c.fReadAODData ;
3079041e 95 fCFManagerProtons = c.fCFManagerProtons;
96 fCFManagerAntiProtons = c.fCFManagerAntiProtons;
b5fc8a3e 97 fQAHistList = c.fQAHistList ;
98 fHistEventsProcessed = c.fHistEventsProcessed;
99 }
100 return *this;
101}
102
103//___________________________________________________________________________
104AliProtonCorrectionTask::AliProtonCorrectionTask(const AliProtonCorrectionTask& c) :
105 AliAnalysisTaskSE(c),
106 fReadTPCTracks(c.fReadTPCTracks),
107 fReadAODData(c.fReadAODData),
3079041e 108 fCFManagerProtons(c.fCFManagerProtons),
109 fCFManagerAntiProtons(c.fCFManagerAntiProtons),
b5fc8a3e 110 fQAHistList(c.fQAHistList),
340705b3 111 fHistEventsProcessed(c.fHistEventsProcessed) {
b5fc8a3e 112 //
113 // Copy Constructor
114 //
115}
116
117//___________________________________________________________________________
118AliProtonCorrectionTask::~AliProtonCorrectionTask() {
119 //
120 //destructor
121 //
122 Info("~AliProtonCorrectionTask","Calling Destructor");
3079041e 123 if (fCFManagerProtons) delete fCFManagerProtons ;
124 if (fCFManagerAntiProtons) delete fCFManagerAntiProtons ;
b5fc8a3e 125 if (fHistEventsProcessed) delete fHistEventsProcessed ;
126 if (fQAHistList) {fQAHistList->Clear(); delete fQAHistList;}
127}
128
129//_________________________________________________
340705b3 130void AliProtonCorrectionTask::UserExec(Option_t *) {
b5fc8a3e 131 //
132 // Main loop function
133 //
134 Info("UserExec","") ;
135
136 AliVEvent* fEvent = fInputEvent ;
137 AliVParticle* track ;
138
139 if (!fEvent) {
140 Error("UserExec","NO EVENT FOUND!");
141 return;
142 }
143
144 if (!fMCEvent) Error("UserExec","NO MC INFO FOUND");
145
146 //pass the MC evt handler to the cuts that need it
3079041e 147 fCFManagerProtons->SetEventInfo(fMCEvent);
148 fCFManagerAntiProtons->SetEventInfo(fMCEvent);
b5fc8a3e 149
150 // MC-event selection
151 Double_t containerInput[2] ;
152
3079041e 153 //__________________________________________________________//
b5fc8a3e 154 //loop on the MC event
155 for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) {
156 AliMCParticle *mcPart = fMCEvent->GetTrack(ipart);
157
158 //check the MC-level cuts
3079041e 159 if (!fCFManagerProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
b5fc8a3e 160
161 containerInput[0] = Rapidity(mcPart->Px(),
162 mcPart->Py(),
163 mcPart->Pz());
164 containerInput[1] = (Float_t)mcPart->Pt();
165 //fill the container for Gen-level selection
3079041e 166 fCFManagerProtons->GetParticleContainer()->Fill(containerInput,kStepGenerated);
b5fc8a3e 167
168 //check the Acceptance-level cuts
3079041e 169 if (!fCFManagerProtons->CheckParticleCuts(AliCFManager::kPartAccCuts,mcPart)) continue;
b5fc8a3e 170 //fill the container for Acceptance-level selection
3079041e 171 fCFManagerProtons->GetParticleContainer()->Fill(containerInput,kStepReconstructible);
172 }//loop over MC particles (protons)
173 //__________________________________________________________//
174 //loop on the MC event
175 for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) {
176 AliMCParticle *mcPart = fMCEvent->GetTrack(ipart);
177
178 //check the MC-level cuts
179 if (!fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
180
181 containerInput[0] = Rapidity(mcPart->Px(),
182 mcPart->Py(),
183 mcPart->Pz());
184 containerInput[1] = (Float_t)mcPart->Pt();
185 //fill the container for Gen-level selection
186 fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,kStepGenerated);
187
188 //check the Acceptance-level cuts
189 if (!fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartAccCuts,mcPart)) continue;
190 //fill the container for Acceptance-level selection
191 fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,kStepReconstructible);
192 }//loop over MC particles (antiprotons)
193
194 //__________________________________________________________//
195 //ESD track loop
196 for (Int_t iTrack = 0; iTrack<fEvent->GetNumberOfTracks(); iTrack++) {
197 track = fEvent->GetTrack(iTrack);
198
199 if (fReadTPCTracks) {
200 if (fReadAODData) {
201 Error("UserExec","TPC-only tracks are not supported with AOD");
202 return ;
203 }
204 AliESDtrack* esdTrack = (AliESDtrack*) track;
205 AliESDtrack* esdTrackTPC = new AliESDtrack();
206 if (!esdTrack->FillTPCOnlyTrack(*esdTrackTPC)) {
207 Error("UserExec","Could not retrieve TPC info");
208 continue;
209 }
210 track = esdTrackTPC ;
211 }
b5fc8a3e 212
3079041e 213 if (!fCFManagerProtons->CheckParticleCuts(AliCFManager::kPartRecCuts,track)) continue;
214
215 // is track associated to particle ?
216 Int_t label = track->GetLabel();
217 if (label<0) continue;
218 AliMCParticle *mcPart = fMCEvent->GetTrack(label);
219
220 // check if this track was part of the signal
221 if (!fCFManagerProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
222
223 //fill the container
224 containerInput[0] = Rapidity(track->Px(),
225 track->Py(),
226 track->Pz());
227 containerInput[1] = track->Pt();
228 fCFManagerProtons->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
229
230 if (!fCFManagerProtons->CheckParticleCuts(AliCFManager::kPartSelCuts,track)) continue ;
231 fCFManagerProtons->GetParticleContainer()->Fill(containerInput,kStepSelected);
232
233 if (fReadTPCTracks) delete track;
234 }
235 //__________________________________________________________//
340705b3 236 //ESD track loop
b5fc8a3e 237 for (Int_t iTrack = 0; iTrack<fEvent->GetNumberOfTracks(); iTrack++) {
b5fc8a3e 238 track = fEvent->GetTrack(iTrack);
239
240 if (fReadTPCTracks) {
241 if (fReadAODData) {
242 Error("UserExec","TPC-only tracks are not supported with AOD");
243 return ;
244 }
245 AliESDtrack* esdTrack = (AliESDtrack*) track;
246 AliESDtrack* esdTrackTPC = new AliESDtrack();
247 if (!esdTrack->FillTPCOnlyTrack(*esdTrackTPC)) {
248 Error("UserExec","Could not retrieve TPC info");
249 continue;
250 }
251 track = esdTrackTPC ;
252 }
253
3079041e 254 if (!fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartRecCuts,track)) continue;
b5fc8a3e 255
256 // is track associated to particle ?
b5fc8a3e 257 Int_t label = track->GetLabel();
b5fc8a3e 258 if (label<0) continue;
259 AliMCParticle *mcPart = fMCEvent->GetTrack(label);
260
261 // check if this track was part of the signal
3079041e 262 if (!fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
b5fc8a3e 263
264 //fill the container
b5fc8a3e 265 containerInput[0] = Rapidity(track->Px(),
266 track->Py(),
267 track->Pz());
268 containerInput[1] = track->Pt();
3079041e 269 fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
b5fc8a3e 270
3079041e 271 if (!fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartSelCuts,track)) continue ;
272 fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,kStepSelected);
b5fc8a3e 273
274 if (fReadTPCTracks) delete track;
275 }
276
277 fHistEventsProcessed->Fill(0);
278
279 /* PostData(0) is taken care of by AliAnalysisTaskSE */
280 PostData(1,fHistEventsProcessed) ;
3079041e 281 PostData(2,fCFManagerProtons->GetParticleContainer()) ;
282 PostData(3,fCFManagerAntiProtons->GetParticleContainer()) ;
283 PostData(4,fQAHistList) ;
b5fc8a3e 284}
285
286
287//___________________________________________________________________________
340705b3 288void AliProtonCorrectionTask::Terminate(Option_t*) {
b5fc8a3e 289 // The Terminate() function is the last function to be called during
290 // a query. It always runs on the client, it can be used to present
291 // the results graphically or save the results to file.
292
293 Info("Terminate","");
294 AliAnalysisTaskSE::Terminate();
295
296 //draw some example plots....
297
298 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
299
300 TH1D* h00 = cont->ShowProjection(0,0) ;
301 TH1D* h01 = cont->ShowProjection(0,1) ;
302 TH1D* h02 = cont->ShowProjection(0,2) ;
303 TH1D* h03 = cont->ShowProjection(0,3) ;
304
305 TH1D* h10 = cont->ShowProjection(1,0) ;
306 TH1D* h11 = cont->ShowProjection(1,1) ;
307 TH1D* h12 = cont->ShowProjection(1,2) ;
308 TH1D* h13 = cont->ShowProjection(1,3) ;
309
b5fc8a3e 310 TCanvas * c =new TCanvas("c","",1400,800);
311 c->Divide(4,2);
312
3079041e 313 c->cd(1); h00->Draw("p");
314 c->cd(2); h01->Draw("p");
315 c->cd(3); h02->Draw("p");
316 c->cd(4); h03->Draw("p");
317 c->cd(5); h10->Draw("p");
318 c->cd(6); h11->Draw("p");
319 c->cd(7); h12->Draw("p");
320 c->cd(8); h13->Draw("p");
b5fc8a3e 321
322 c->SaveAs("plots.eps");
323}
324
325
326//___________________________________________________________________________
327void AliProtonCorrectionTask::UserCreateOutputObjects() {
328 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
329 //TO BE SET BEFORE THE EXECUTION OF THE TASK
330 //
331 Info("CreateOutputObjects","CreateOutputObjects of task %s", GetName());
332
333 //slot #1
334 OpenFile(1);
335 fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
336
337// OpenFile(2);
338// OpenFile(3);
339}
340
341//___________________________________________________________________________
342Double_t AliProtonCorrectionTask::Rapidity(Double_t Px, Double_t Py, Double_t Pz) {
343 //returns the rapidity of the (anti)proton
344 Double_t fMass = 9.38270000000000048e-01;
345
346 Double_t P = TMath::Sqrt(TMath::Power(Px,2) +
347 TMath::Power(Py,2) +
348 TMath::Power(Pz,2));
349 Double_t energy = TMath::Sqrt(P*P + fMass*fMass);
350 Double_t y = -999;
351 if(energy != Pz)
352 y = 0.5*TMath::Log((energy + Pz)/(energy - Pz));
353
354 return y;
355}
356
357#endif