Compliance with AliAnalysisTaskSE ; upgraded support for QA.
[u/mrichter/AliRoot.git] / CORRFW / test / AliCFV0Task.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 on AliEn (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 ALICFV0TASK_CXX
28 #define ALICFV0TASK_CXX
29
30 #include "AliCFV0Task.h"
31 #include "TCanvas.h"
32 #include "AliStack.h"
33 #include "TParticle.h"
34 #include "TH1I.h"
35 #include "AliMCEvent.h"
36 #include "AliCFManager.h"
37 #include "AliCFContainer.h"
38 #include "AliESDtrack.h"
39 #include "AliLog.h"
40 #include "AliESDv0.h"
41 #include "AliV0vertexer.h"
42 #include "AliCFPair.h"
43 #include "AliESDEvent.h"
44 #include "TChain.h"
45 #include "AliCFParticleGenCuts.h"
46
47 //__________________________________________________________________________
48 AliCFV0Task::AliCFV0Task() :
49   AliAnalysisTaskSE(),
50   fRebuildV0s(0),
51   fV0PDG(310),
52   fCFManager(0x0),
53   fHistEventsProcessed(0x0)
54 {
55   //
56   //Default ctor
57   //
58 }
59 //___________________________________________________________________________
60 AliCFV0Task::AliCFV0Task(const Char_t* name) :
61   AliAnalysisTaskSE(name),
62   fRebuildV0s(0),
63   fV0PDG(310),
64   fCFManager(0x0),
65   fHistEventsProcessed(0x0)
66 {
67   //
68   // Constructor. Initialization of Inputs and Outputs
69   //
70   Info("AliCFV0Task","Calling Constructor");
71   /*
72     DefineInput(0) and DefineOutput(0)
73     are taken care of by AliAnalysisTaskSE constructor
74   */
75   DefineOutput(1,TH1I::Class());
76   DefineOutput(2,AliCFContainer::Class());
77 //   DefineOutput(3,TList::Class());
78 }
79
80 //___________________________________________________________________________
81 AliCFV0Task& AliCFV0Task::operator=(const AliCFV0Task& c) 
82 {
83   //
84   // Assignment operator
85   //
86   if (this!=&c) {
87     AliAnalysisTaskSE::operator=(c) ;
88     fRebuildV0s = c.fRebuildV0s;
89     fV0PDG      = c.fV0PDG;
90     fCFManager  = c.fCFManager;
91     fHistEventsProcessed = c.fHistEventsProcessed;
92   }
93   return *this;
94 }
95
96 //___________________________________________________________________________
97 AliCFV0Task::AliCFV0Task(const AliCFV0Task& c) :
98   AliAnalysisTaskSE(c),
99   fRebuildV0s(c.fRebuildV0s),
100   fV0PDG(c.fV0PDG),
101   fCFManager(c.fCFManager),
102   fHistEventsProcessed(c.fHistEventsProcessed)
103 {
104   //
105   // Copy Constructor
106   //
107 }
108
109 //___________________________________________________________________________
110 AliCFV0Task::~AliCFV0Task() {
111   //
112   //destructor
113   //
114   Info("~AliCFV0Task","Calling Destructor");
115   if (fCFManager)           delete fCFManager ;
116   if (fHistEventsProcessed) delete fHistEventsProcessed ;
117 }
118
119 //_________________________________________________
120 void AliCFV0Task::UserExec(Option_t *)
121 {
122   //
123   // Main loop function
124   //
125   Info("UserExec","") ;
126
127   AliESDEvent* fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
128   if (!fESD) {
129     Error("UserExec","NO ESD FOUND!");
130     return;
131   }
132
133   if (!fMCEvent) Error("UserExec","NO MC INFO FOUND!");
134   fCFManager->SetEventInfo(fMCEvent);
135
136   // MC-event selection
137   Double_t containerInput[2] ;
138         
139   //loop on the MC event
140   Info("UserExec","Looping on MC event");
141   for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) { 
142     AliMCParticle *mcPart  = fMCEvent->GetTrack(ipart);
143
144     //check the MC-level cuts
145     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;
146     containerInput[0] = mcPart->Pt();
147     containerInput[1] = mcPart->Eta() ;
148     //fill the container for Gen-level selection
149     fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);
150     
151     //check the Acceptance-level cuts
152     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcPart)) continue;
153     //fill the container for Acceptance-level selection
154     fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructible);
155   }    
156
157
158   //Now go to rec level
159   Info("UserExec","Looping on ESD event");
160
161   if (fRebuildV0s) RebuildV0s(fESD) ;
162   Info("UserExec","There are %d V0s in event",fESD->GetNumberOfV0s());
163   for (Int_t iV0 = 0; iV0<fESD->GetNumberOfV0s(); iV0++) {
164
165     AliESDv0* esdV0 = fESD->GetV0(iV0);
166
167     //check if mother reconstructed V0 can be associated to a MC V0
168     Int_t labMCV0 = IsMcV0(esdV0) ;
169     if (labMCV0 <0) continue;
170
171     esdV0->ChangeMassHypothesis(fV0PDG); //important to do that before entering the cut check
172
173     AliCFPair pair(esdV0,fESD);
174     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,&pair)) continue;
175
176     //check if associated MC v0 passes the cuts
177     AliMCParticle* mcV0 = fMCEvent->GetTrack(labMCV0);
178     if (!mcV0) continue;
179
180     Info("UserExec","is v0 primary : %d",AliCFParticleGenCuts::IsPrimary(mcV0,fMCEvent->Stack()));
181
182     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcV0)) continue; 
183     
184     //fill the container
185     Double_t mom[3];
186     esdV0->GetPxPyPz(mom[0],mom[1],mom[2]);
187     Double_t pt2 = mom[0]*mom[0]+mom[1]*mom[1] ;
188     Double_t pt  = TMath::Sqrt(pt2);
189     Double_t energy  = TMath::Sqrt(pt2 + mom[2]*mom[2] + TMath::Power(TDatabasePDG::Instance()->GetParticle(fV0PDG)->Mass(),2));
190
191     containerInput[0] = pt ;
192     containerInput[1] = GetRapidity(energy,mom[2]);
193     fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;   
194     if (!fCFManager->CheckParticleCuts(AliCFManager::kPartSelCuts,&pair)) continue ;
195     fCFManager->GetParticleContainer()->Fill(containerInput,kStepSelected);
196   }
197   
198   fHistEventsProcessed->Fill(0);
199   /* PostData(0) is taken care of by AliAnalysisTaskSE */
200   PostData(1,fHistEventsProcessed) ;
201   PostData(2,fCFManager->GetParticleContainer()) ;
202 }
203
204
205 //___________________________________________________________________________
206 void AliCFV0Task::Terminate(Option_t*)
207 {
208   // The Terminate() function is the last function to be called during
209   // a query. It always runs on the client, it can be used to present
210   // the results graphically or save the results to file.
211
212   Info("Terminate","");
213   AliAnalysisTaskSE::Terminate();
214
215
216   //draw some example plots....
217
218   AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
219
220   TH1D* h00 =   cont->ShowProjection(0,0) ;
221   TH1D* h01 =   cont->ShowProjection(0,1) ;
222   TH1D* h02 =   cont->ShowProjection(0,2) ;
223   TH1D* h03 =   cont->ShowProjection(0,3) ;
224
225   TH1D* h10 =   cont->ShowProjection(1,0) ;
226   TH1D* h11 =   cont->ShowProjection(1,1) ;
227   TH1D* h12 =   cont->ShowProjection(1,2) ;
228   TH1D* h13 =   cont->ShowProjection(1,3) ;
229
230   Double_t max1 = h00->GetMaximum();
231   Double_t max2 = h10->GetMaximum();
232
233   h00->GetYaxis()->SetRangeUser(0,max1*1.2);
234   h01->GetYaxis()->SetRangeUser(0,max1*1.2);
235   h02->GetYaxis()->SetRangeUser(0,max1*1.2);
236   h03->GetYaxis()->SetRangeUser(0,max1*1.2);
237
238   h10->GetYaxis()->SetRangeUser(0,max2*1.2);
239   h11->GetYaxis()->SetRangeUser(0,max2*1.2);
240   h12->GetYaxis()->SetRangeUser(0,max2*1.2);
241   h13->GetYaxis()->SetRangeUser(0,max2*1.2);
242
243   h00->SetMarkerStyle(23) ;
244   h01->SetMarkerStyle(24) ;
245   h02->SetMarkerStyle(25) ;
246   h03->SetMarkerStyle(26) ;
247
248   h10->SetMarkerStyle(23) ;
249   h11->SetMarkerStyle(24) ;
250   h12->SetMarkerStyle(25) ;
251   h13->SetMarkerStyle(26) ;
252
253   TCanvas * c =new TCanvas("c","",1400,800);
254   c->Divide(4,2);
255
256   c->cd(1);
257   h00->Draw("p");
258   c->cd(2);
259   h01->Draw("p");
260   c->cd(3);
261   h02->Draw("p");
262   c->cd(4);
263   h03->Draw("p");
264   c->cd(5);
265   h10->Draw("p");
266   c->cd(6);
267   h11->Draw("p");
268   c->cd(7);
269   h12->Draw("p");
270   c->cd(8);
271   h13->Draw("p");
272
273   c->SaveAs("plots.eps");
274 }
275
276 //___________________________________________________________________________
277 void AliCFV0Task::UserCreateOutputObjects() {
278   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
279   //TO BE SET BEFORE THE EXECUTION OF THE TASK
280   //
281   Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
282
283   //slot #1
284   OpenFile(1);
285   fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
286 }
287
288 //___________________________________________________________________________
289 Int_t AliCFV0Task::GetV0Label(UInt_t lab1, UInt_t lab2) const {
290   //
291   // returns the label of the V0, given the labels of the 2 daughter tracks
292   // returns -1 if the V0 is fake
293   //
294   
295   AliStack* stack = fMCEvent->Stack();
296   TParticle* part1 = stack->Particle(lab1) ;
297   TParticle* part2 = stack->Particle(lab2) ;
298   
299   Int_t part1MotherLab=part1->GetFirstMother();
300   Int_t part2MotherLab=part2->GetFirstMother();
301
302   if (part1MotherLab==-1 || part2MotherLab==-1) return -1 ;
303   if (part1MotherLab != part2MotherLab )        return -1 ;
304   
305   if (stack->Particle(part1MotherLab)->GetPdgCode() != fV0PDG ) return -1 ;
306
307   switch (fV0PDG) {
308   case kK0Short : 
309     if ( (part1->GetPdgCode()==  -211 && part2->GetPdgCode()==  211)  || 
310          (part1->GetPdgCode()==   211 && part2->GetPdgCode()== -211)   )   return part1MotherLab ;
311     break ;
312   case kLambda0 :
313     if ( (part1->GetPdgCode()==  -211 && part2->GetPdgCode()== 2212)  || 
314          (part1->GetPdgCode()==  2212 && part2->GetPdgCode()== -211)   )   return part1MotherLab ;
315     break ;
316   case kLambda0Bar :
317     if ( (part1->GetPdgCode()==   211 && part2->GetPdgCode()== -2212) || 
318          (part1->GetPdgCode()== -2212 && part2->GetPdgCode()==   211)  )   return part1MotherLab ;
319     break ;
320   default :
321     return -1;
322     break ;
323   }
324   
325   return -1 ;
326 }
327
328 //___________________________________________________________________________
329 Double_t AliCFV0Task::GetRapidity(Double_t energy, Double_t pz) {
330   if (energy == pz || energy == -pz) {
331     printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
332     return 999;
333   }
334   if (energy < pz) {
335     printf("GetRapidity : ERROR : rapidity for 4-vector with E = Pz -- infinite result");
336     return 999;
337   }
338   Double_t y = 0.5 * log((energy + pz) / (energy - pz));
339   return y;
340
341
342 //___________________________________________________________________________
343 Int_t AliCFV0Task::IsMcV0(AliESDv0* v0) const {
344
345   AliESDEvent* fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
346
347   Int_t nindex = v0->GetNindex();
348   Int_t pindex = v0->GetPindex();
349   AliESDtrack *nTrack = fESD->GetTrack(nindex) ;
350   AliESDtrack *pTrack = fESD->GetTrack(pindex) ;
351   
352   if (!nTrack || !pTrack) return -1 ;
353
354   Int_t nlab  = nTrack->GetLabel() ;
355   Int_t plab  = pTrack->GetLabel() ;
356   
357   if (nlab <0 || plab <0) return -1 ;
358
359   return GetV0Label((UInt_t)nlab,(UInt_t)plab) ;
360 }
361
362
363 //___________________________________________________________________________
364 void AliCFV0Task::RebuildV0s(AliESDEvent* fESD) {
365   
366   fESD->ResetV0s();
367
368   //These are pp cuts : to change if Pb+Pb !
369   Double_t cuts[]={33,  // max. allowed chi2
370                    0.01,// min. allowed negative daughter's impact parameter 
371                    0.01,// min. allowed positive daughter's impact parameter 
372                    0.5,// max. allowed DCA between the daughter tracks
373                    0.98,// max. allowed cosine of V0's pointing angle
374                    0.2,  // min. radius of the fiducial volume
375                    100.   // max. radius of the fiducial volume
376   };
377   //   AliV0vertexer* v0Vertexer = new AliV0vertexer(cuts);
378   //   v0Vertexer->SetVertex(primVertexPos);
379   AliV0vertexer* v0Vertexer = new AliV0vertexer();
380   v0Vertexer->SetCuts(cuts) ;
381   v0Vertexer->Tracks2V0vertices(fESD);
382   delete v0Vertexer ;
383 }
384
385 #endif