bug fix in the vertex selection
[u/mrichter/AliRoot.git] / PWG2 / AliProtonCorrectionTask.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 //-----------------------------------------------------------------------
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
47 ClassImp(AliProtonCorrectionTask)
48
49 //__________________________________________________________________________
50 AliProtonCorrectionTask::AliProtonCorrectionTask() :
51   fReadTPCTracks(0),
52   fReadAODData(0),
53   fCFManagerProtons(0x0),
54   fCFManagerAntiProtons(0x0),
55   fQAHistList(0x0),
56   fHistEventsProcessed(0x0)
57 {
58   //
59   //Default ctor
60   //
61 }
62 //___________________________________________________________________________
63 AliProtonCorrectionTask::AliProtonCorrectionTask(const Char_t* name) :
64   AliAnalysisTaskSE(name),
65   fReadTPCTracks(0),
66   fReadAODData(0),
67   fCFManagerProtons(0x0),
68   fCFManagerAntiProtons(0x0),
69   fQAHistList(0x0),
70   fHistEventsProcessed(0x0) {
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());
82   DefineOutput(3,AliCFContainer::Class());
83   DefineOutput(4,TList::Class());
84 }
85
86 //___________________________________________________________________________
87 AliProtonCorrectionTask& AliProtonCorrectionTask::operator=(const AliProtonCorrectionTask& c) {
88   //
89   // Assignment operator
90   //
91   if (this!=&c) {
92     AliAnalysisTaskSE::operator=(c) ;
93     fReadTPCTracks = c.fReadTPCTracks ;
94     fReadAODData = c.fReadAODData ;
95     fCFManagerProtons  = c.fCFManagerProtons;
96     fCFManagerAntiProtons  = c.fCFManagerAntiProtons;
97     fQAHistList = c.fQAHistList ;
98     fHistEventsProcessed = c.fHistEventsProcessed;
99   }
100   return *this;
101 }
102
103 //___________________________________________________________________________
104 AliProtonCorrectionTask::AliProtonCorrectionTask(const AliProtonCorrectionTask& c) :
105   AliAnalysisTaskSE(c),
106   fReadTPCTracks(c.fReadTPCTracks),
107   fReadAODData(c.fReadAODData),
108   fCFManagerProtons(c.fCFManagerProtons),
109   fCFManagerAntiProtons(c.fCFManagerAntiProtons),
110   fQAHistList(c.fQAHistList),
111   fHistEventsProcessed(c.fHistEventsProcessed) {
112   //
113   // Copy Constructor
114   //
115 }
116
117 //___________________________________________________________________________
118 AliProtonCorrectionTask::~AliProtonCorrectionTask() {
119   //
120   //destructor
121   //
122   Info("~AliProtonCorrectionTask","Calling Destructor");
123   if (fCFManagerProtons) delete fCFManagerProtons ;
124   if (fCFManagerAntiProtons) delete fCFManagerAntiProtons ;
125   if (fHistEventsProcessed) delete fHistEventsProcessed ;
126   if (fQAHistList) {fQAHistList->Clear(); delete fQAHistList;}
127 }
128
129 //_________________________________________________
130 void AliProtonCorrectionTask::UserExec(Option_t *) {
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 
147   fCFManagerProtons->SetEventInfo(fMCEvent);
148   fCFManagerAntiProtons->SetEventInfo(fMCEvent);
149
150   // MC-event selection
151   Double_t containerInput[2] ;
152         
153   //__________________________________________________________//
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
159     if (!fCFManagerProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
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
166     fCFManagerProtons->GetParticleContainer()->Fill(containerInput,kStepGenerated);
167
168     //check the Acceptance-level cuts
169     if (!fCFManagerProtons->CheckParticleCuts(AliCFManager::kPartAccCuts,mcPart)) continue;
170     //fill the container for Acceptance-level selection
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     }
212
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   //__________________________________________________________//
236   //ESD track loop
237   for (Int_t iTrack = 0; iTrack<fEvent->GetNumberOfTracks(); iTrack++) {
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
254     if (!fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartRecCuts,track)) continue;
255     
256     // is track associated to particle ?
257     Int_t label = track->GetLabel();
258     if (label<0) continue;
259     AliMCParticle *mcPart  = fMCEvent->GetTrack(label);
260     
261     // check if this track was part of the signal
262     if (!fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue; 
263     
264     //fill the container
265     containerInput[0] = Rapidity(track->Px(),
266                                  track->Py(),
267                                  track->Pz());
268     containerInput[1] = track->Pt();
269     fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;   
270
271     if (!fCFManagerAntiProtons->CheckParticleCuts(AliCFManager::kPartSelCuts,track)) continue ;
272     fCFManagerAntiProtons->GetParticleContainer()->Fill(containerInput,kStepSelected);
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) ;
281   PostData(2,fCFManagerProtons->GetParticleContainer()) ;
282   PostData(3,fCFManagerAntiProtons->GetParticleContainer()) ;
283   PostData(4,fQAHistList) ;
284 }
285
286
287 //___________________________________________________________________________
288 void AliProtonCorrectionTask::Terminate(Option_t*) {
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
310   TCanvas * c =new TCanvas("c","",1400,800);
311   c->Divide(4,2);
312
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");
321
322   c->SaveAs("plots.eps");
323 }
324
325
326 //___________________________________________________________________________
327 void 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 //___________________________________________________________________________ 
342 Double_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