First version of the task for the corrections - proton analysis
[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),
53 fCFManager(0x0),
54 fQAHistList(0x0),
55 fHistEventsProcessed(0x0)
56{
57 //
58 //Default ctor
59 //
60}
61//___________________________________________________________________________
62AliProtonCorrectionTask::AliProtonCorrectionTask(const Char_t* name) :
63 AliAnalysisTaskSE(name),
64 fReadTPCTracks(0),
65 fReadAODData(0),
66 fCFManager(0x0),
67 fQAHistList(0x0),
68 fHistEventsProcessed(0x0)
69{
70 //
71 // Constructor. Initialization of Inputs and Outputs
72 //
73 Info("AliProtonCorrectionTask","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,AliCFContainer::Class());
81 DefineOutput(3,TList::Class());
82}
83
84//___________________________________________________________________________
85AliProtonCorrectionTask& AliProtonCorrectionTask::operator=(const AliProtonCorrectionTask& c)
86{
87 //
88 // Assignment operator
89 //
90 if (this!=&c) {
91 AliAnalysisTaskSE::operator=(c) ;
92 fReadTPCTracks = c.fReadTPCTracks ;
93 fReadAODData = c.fReadAODData ;
94 fCFManager = c.fCFManager;
95 fQAHistList = c.fQAHistList ;
96 fHistEventsProcessed = c.fHistEventsProcessed;
97 }
98 return *this;
99}
100
101//___________________________________________________________________________
102AliProtonCorrectionTask::AliProtonCorrectionTask(const AliProtonCorrectionTask& c) :
103 AliAnalysisTaskSE(c),
104 fReadTPCTracks(c.fReadTPCTracks),
105 fReadAODData(c.fReadAODData),
106 fCFManager(c.fCFManager),
107 fQAHistList(c.fQAHistList),
108 fHistEventsProcessed(c.fHistEventsProcessed)
109{
110 //
111 // Copy Constructor
112 //
113}
114
115//___________________________________________________________________________
116AliProtonCorrectionTask::~AliProtonCorrectionTask() {
117 //
118 //destructor
119 //
120 Info("~AliProtonCorrectionTask","Calling Destructor");
121 if (fCFManager) delete fCFManager ;
122 if (fHistEventsProcessed) delete fHistEventsProcessed ;
123 if (fQAHistList) {fQAHistList->Clear(); delete fQAHistList;}
124}
125
126//_________________________________________________
127void AliProtonCorrectionTask::UserExec(Option_t *)
128{
129 //
130 // Main loop function
131 //
132 Info("UserExec","") ;
133
134 AliVEvent* fEvent = fInputEvent ;
135 AliVParticle* track ;
136
137 if (!fEvent) {
138 Error("UserExec","NO EVENT FOUND!");
139 return;
140 }
141
142 if (!fMCEvent) Error("UserExec","NO MC INFO FOUND");
143
144 //pass the MC evt handler to the cuts that need it
145 fCFManager->SetEventInfo(fMCEvent);
146
147 // MC-event selection
148 Double_t containerInput[2] ;
149
150 //loop on the MC event
151 for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) {
152 AliMCParticle *mcPart = fMCEvent->GetTrack(ipart);
153
154 //check the MC-level cuts
155 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
156
157 containerInput[0] = Rapidity(mcPart->Px(),
158 mcPart->Py(),
159 mcPart->Pz());
160 containerInput[1] = (Float_t)mcPart->Pt();
161 //fill the container for Gen-level selection
162 fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);
163
164 //check the Acceptance-level cuts
165 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcPart)) continue;
166 //fill the container for Acceptance-level selection
167 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructible);
168 }
169
170 //Now go to rec level
171 for (Int_t iTrack = 0; iTrack<fEvent->GetNumberOfTracks(); iTrack++) {
172
173 track = fEvent->GetTrack(iTrack);
174
175 if (fReadTPCTracks) {
176 if (fReadAODData) {
177 Error("UserExec","TPC-only tracks are not supported with AOD");
178 return ;
179 }
180 AliESDtrack* esdTrack = (AliESDtrack*) track;
181 AliESDtrack* esdTrackTPC = new AliESDtrack();
182 if (!esdTrack->FillTPCOnlyTrack(*esdTrackTPC)) {
183 Error("UserExec","Could not retrieve TPC info");
184 continue;
185 }
186 track = esdTrackTPC ;
187 }
188
189 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,track)) continue;
190
191 // is track associated to particle ?
192
193 Int_t label = track->GetLabel();
194
195 if (label<0) continue;
196 AliMCParticle *mcPart = fMCEvent->GetTrack(label);
197
198 // check if this track was part of the signal
199 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
200
201 //fill the container
202 Double_t mom[3];
203 track->PxPyPz(mom);
204 Double_t pt=TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
205 containerInput[0] = Rapidity(track->Px(),
206 track->Py(),
207 track->Pz());
208 containerInput[1] = track->Pt();
209 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
210
211 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,track)) continue ;
212 fCFManager->GetParticleContainer()->Fill(containerInput,kStepSelected);
213
214 if (fReadTPCTracks) delete track;
215 }
216
217 fHistEventsProcessed->Fill(0);
218
219 /* PostData(0) is taken care of by AliAnalysisTaskSE */
220 PostData(1,fHistEventsProcessed) ;
221 PostData(2,fCFManager->GetParticleContainer()) ;
222 PostData(3,fQAHistList) ;
223}
224
225
226//___________________________________________________________________________
227void AliProtonCorrectionTask::Terminate(Option_t*)
228{
229 // The Terminate() function is the last function to be called during
230 // a query. It always runs on the client, it can be used to present
231 // the results graphically or save the results to file.
232
233 Info("Terminate","");
234 AliAnalysisTaskSE::Terminate();
235
236 //draw some example plots....
237
238 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
239
240 TH1D* h00 = cont->ShowProjection(0,0) ;
241 TH1D* h01 = cont->ShowProjection(0,1) ;
242 TH1D* h02 = cont->ShowProjection(0,2) ;
243 TH1D* h03 = cont->ShowProjection(0,3) ;
244
245 TH1D* h10 = cont->ShowProjection(1,0) ;
246 TH1D* h11 = cont->ShowProjection(1,1) ;
247 TH1D* h12 = cont->ShowProjection(1,2) ;
248 TH1D* h13 = cont->ShowProjection(1,3) ;
249
250 Double_t max1 = h00->GetMaximum();
251 Double_t max2 = h10->GetMaximum();
252
253 /*h00->GetYaxis()->SetRangeUser(0,max1*1.2);
254 h01->GetYaxis()->SetRangeUser(0,max1*1.2);
255 h02->GetYaxis()->SetRangeUser(0,max1*1.2);
256 h03->GetYaxis()->SetRangeUser(0,max1*1.2);
257
258 h10->GetYaxis()->SetRangeUser(0,max2*1.2);
259 h11->GetYaxis()->SetRangeUser(0,max2*1.2);
260 h12->GetYaxis()->SetRangeUser(0,max2*1.2);
261 h13->GetYaxis()->SetRangeUser(0,max2*1.2);
262
263 h00->SetMarkerStyle(23) ;
264 h01->SetMarkerStyle(24) ;
265 h02->SetMarkerStyle(25) ;
266 h03->SetMarkerStyle(26) ;
267
268 h10->SetMarkerStyle(23) ;
269 h11->SetMarkerStyle(24) ;
270 h12->SetMarkerStyle(25) ;
271 h13->SetMarkerStyle(26) ;*/
272
273 TCanvas * c =new TCanvas("c","",1400,800);
274 c->Divide(4,2);
275
276 c->cd(1);
277 h00->Draw("p");
278 c->cd(2);
279 h01->Draw("p");
280 c->cd(3);
281 h02->Draw("p");
282 c->cd(4);
283 h03->Draw("p");
284 c->cd(5);
285 h10->Draw("p");
286 c->cd(6);
287 h11->Draw("p");
288 c->cd(7);
289 h12->Draw("p");
290 c->cd(8);
291 h13->Draw("p");
292
293 c->SaveAs("plots.eps");
294}
295
296
297//___________________________________________________________________________
298void AliProtonCorrectionTask::UserCreateOutputObjects() {
299 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
300 //TO BE SET BEFORE THE EXECUTION OF THE TASK
301 //
302 Info("CreateOutputObjects","CreateOutputObjects of task %s", GetName());
303
304 //slot #1
305 OpenFile(1);
306 fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
307
308// OpenFile(2);
309// OpenFile(3);
310}
311
312//___________________________________________________________________________
313Double_t AliProtonCorrectionTask::Rapidity(Double_t Px, Double_t Py, Double_t Pz) {
314 //returns the rapidity of the (anti)proton
315 Double_t fMass = 9.38270000000000048e-01;
316
317 Double_t P = TMath::Sqrt(TMath::Power(Px,2) +
318 TMath::Power(Py,2) +
319 TMath::Power(Pz,2));
320 Double_t energy = TMath::Sqrt(P*P + fMass*fMass);
321 Double_t y = -999;
322 if(energy != Pz)
323 y = 0.5*TMath::Log((energy + Pz)/(energy - Pz));
324
325 return y;
326}
327
328#endif