]>
Commit | Line | Data |
---|---|---|
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 | //__________________________________________________________________________ | |
51 | AliCFV0Task::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 | //___________________________________________________________________________ | |
63 | AliCFV0Task::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 | //___________________________________________________________________________ | |
84 | AliCFV0Task& 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 | //___________________________________________________________________________ | |
100 | AliCFV0Task::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 | //___________________________________________________________________________ | |
113 | AliCFV0Task::~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 | 123 | void 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 | //___________________________________________________________________________ | |
247 | void 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 | 318 | void 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 | 330 | Int_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 | //___________________________________________________________________________ | |
370 | Double_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 | 388 | Int_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 | //___________________________________________________________________________ |
413 | Int_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 | 437 | void 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 |