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