]>
Commit | Line | Data |
---|---|---|
becfe989 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-2009, 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 | // Class for HF corrections as a function of many variables | |
18 | // | |
19 | // Example of task running on AliEn | |
20 | // which provides standard way of calculating acceptance and efficiency | |
21 | // between different steps of the procedure. | |
22 | // The ouptut of the task is a AliCFContainer from which the efficiencies | |
23 | // can be calculated | |
24 | //----------------------------------------------------------------------- | |
25 | // Author : C. Zampolli, CERN, on the basis of R. Vernet's example | |
26 | //----------------------------------------------------------------------- | |
27 | ||
28 | #include <TCanvas.h> | |
29 | #include <TParticle.h> | |
30 | #include <TH1I.h> | |
31 | #include <TStyle.h> | |
32 | ||
33 | #include "AliCFHeavyFlavourTaskMultiVar.h" | |
34 | #include "AliStack.h" | |
35 | #include "AliMCEvent.h" | |
36 | #include "AliCFManager.h" | |
37 | #include "AliCFContainer.h" | |
38 | #include "AliLog.h" | |
39 | #include "AliAODEvent.h" | |
40 | #include "AliAODRecoDecay.h" | |
41 | #include "AliAODRecoDecayHF.h" | |
42 | #include "AliAODRecoDecayHF2Prong.h" | |
43 | #include "AliAODMCParticle.h" | |
44 | ||
45 | //__________________________________________________________________________ | |
46 | AliCFHeavyFlavourTaskMultiVar::AliCFHeavyFlavourTaskMultiVar() : | |
47 | AliAnalysisTaskSE(), | |
48 | fPDG(0), | |
49 | fCFManager(0x0), | |
50 | fHistEventsProcessed(0x0), | |
51 | fCountMC(0), | |
52 | fEvents(0), | |
53 | fFillFromGenerated(kFALSE), | |
54 | fSkipped(0) | |
55 | { | |
56 | // | |
57 | //Default ctor | |
58 | // | |
59 | } | |
60 | //___________________________________________________________________________ | |
61 | AliCFHeavyFlavourTaskMultiVar::AliCFHeavyFlavourTaskMultiVar(const Char_t* name) : | |
62 | AliAnalysisTaskSE(name), | |
63 | fPDG(0), | |
64 | fCFManager(0x0), | |
65 | fHistEventsProcessed(0x0), | |
66 | fCountMC(0), | |
67 | fEvents(0), | |
68 | fFillFromGenerated(kFALSE), | |
69 | fSkipped(0) | |
70 | { | |
71 | // | |
72 | // Constructor. Initialization of Inputs and Outputs | |
73 | // | |
74 | Info("AliCFHeavyFlavourTaskMultiVar","Calling Constructor"); | |
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 | } | |
82 | ||
83 | //___________________________________________________________________________ | |
84 | AliCFHeavyFlavourTaskMultiVar& AliCFHeavyFlavourTaskMultiVar::operator=(const AliCFHeavyFlavourTaskMultiVar& c) | |
85 | { | |
86 | // | |
87 | // Assignment operator | |
88 | // | |
89 | if (this!=&c) { | |
90 | AliAnalysisTaskSE::operator=(c) ; | |
91 | fPDG = c.fPDG; | |
92 | fCFManager = c.fCFManager; | |
93 | fHistEventsProcessed = c.fHistEventsProcessed; | |
94 | } | |
95 | return *this; | |
96 | } | |
97 | ||
98 | //___________________________________________________________________________ | |
99 | AliCFHeavyFlavourTaskMultiVar::AliCFHeavyFlavourTaskMultiVar(const AliCFHeavyFlavourTaskMultiVar& c) : | |
100 | AliAnalysisTaskSE(c), | |
101 | fPDG(c.fPDG), | |
102 | fCFManager(c.fCFManager), | |
103 | fHistEventsProcessed(c.fHistEventsProcessed), | |
104 | fCountMC(c.fCountMC), | |
105 | fEvents(c.fEvents), | |
106 | fFillFromGenerated(c.fFillFromGenerated), | |
107 | fSkipped(c.fSkipped) | |
108 | { | |
109 | // | |
110 | // Copy Constructor | |
111 | // | |
112 | } | |
113 | ||
114 | //___________________________________________________________________________ | |
115 | AliCFHeavyFlavourTaskMultiVar::~AliCFHeavyFlavourTaskMultiVar() { | |
116 | // | |
117 | //destructor | |
118 | // | |
119 | Info("~AliCFHeavyFlavourTaskMultiVar","Calling Destructor"); | |
120 | if (fCFManager) delete fCFManager ; | |
121 | if (fHistEventsProcessed) delete fHistEventsProcessed ; | |
122 | } | |
123 | ||
124 | //_________________________________________________ | |
125 | void AliCFHeavyFlavourTaskMultiVar::UserExec(Option_t *) | |
126 | { | |
127 | // | |
128 | // Main loop function | |
129 | // | |
130 | ||
131 | if (!fInputEvent) { | |
132 | Error("UserExec","NO EVENT FOUND!"); | |
133 | return; | |
134 | } | |
135 | ||
136 | fEvents++; | |
137 | if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents)); | |
138 | AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent); | |
139 | fCFManager->SetEventInfo(aodEvent); | |
140 | ||
141 | // MC-event selection | |
142 | Double_t containerInput[6] ; | |
143 | ||
144 | //loop on the MC event | |
145 | ||
146 | TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName())); | |
147 | if (!mcArray) AliError("Could not find Monte-Carlo in AOD"); | |
148 | Int_t icountMC = 0; | |
149 | ||
150 | for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) { | |
151 | AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart)); | |
152 | if (!mcPart) { | |
153 | AliWarning("Particle not found in tree, skipping"); | |
154 | continue; | |
155 | } | |
156 | ||
157 | // check the MC-level cuts | |
158 | if (!fCFManager->CheckParticleCuts(0, mcPart)) continue; // 0 stands for MC level | |
159 | ||
160 | // fill the container for Gen-level selection | |
161 | Double_t vectorMC[6] = {9999.,9999.,9999.,9999.,9999.,9999.}; | |
162 | if (GetGeneratedValuesFromMCParticle(mcPart, mcArray, vectorMC)){ | |
163 | containerInput[0] = vectorMC[0]; | |
164 | containerInput[1] = vectorMC[1] ; | |
165 | containerInput[2] = vectorMC[2] ; | |
166 | containerInput[3] = vectorMC[3] ; | |
167 | containerInput[4] = vectorMC[4] ; | |
168 | containerInput[5] = vectorMC[5] ; // in micron | |
169 | fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated); | |
170 | icountMC++; | |
171 | } | |
172 | else { | |
173 | AliDebug(3,"Problems in filling the container"); | |
174 | continue; | |
175 | } | |
176 | } | |
177 | ||
178 | AliDebug(2, Form("Found %i MC particles that are D0!!",icountMC)); | |
179 | AliDebug(2, Form("aodEvent = %p ",aodEvent)); | |
180 | ||
181 | // Now go to rec level | |
182 | fCountMC += icountMC; | |
183 | // load heavy flavour vertices | |
184 | ||
185 | TClonesArray *arrayVerticesHF = (TClonesArray*)((aodEvent->GetList())->FindObject("D0toKpi")); | |
186 | if (!arrayVerticesHF) AliError("Could not find array of HF vertices"); | |
187 | AliDebug(2, Form("Found %d vertices",arrayVerticesHF->GetEntriesFast())); | |
188 | ||
189 | for (Int_t iVertex = 0; iVertex<arrayVerticesHF->GetEntriesFast(); iVertex++) { | |
190 | ||
191 | AliAODRecoDecayHF2Prong* vtx = (AliAODRecoDecayHF2Prong*)arrayVerticesHF->At(iVertex); | |
192 | ||
193 | // cuts can't be applied to RecoDecays particles | |
194 | // if (!fCFManager->CheckParticleCuts(1 , vtx)) continue; // 1 stands for AOD level | |
195 | ||
196 | // find associated MC particle | |
197 | Int_t mcLabel = IsMcVtx(vtx) ; | |
198 | if (mcLabel == -1) | |
199 | { | |
200 | AliDebug(2,"No MC particle found"); | |
201 | } | |
202 | else { | |
203 | AliAODMCParticle* mcVtxHF = (AliAODMCParticle*)mcArray->At(mcLabel); | |
204 | if (!mcVtxHF) { | |
205 | AliWarning("Could not find associated MC in AOD MC tree"); | |
206 | continue; | |
207 | } | |
208 | ||
209 | // check if associated MC v0 passes the cuts | |
210 | if (!fCFManager->CheckParticleCuts(0 ,mcVtxHF)) { // 0 stands for MC | |
211 | AliDebug(2, "Skipping the particles due to cuts"); | |
212 | continue; | |
213 | } | |
214 | ||
215 | // fill the container | |
216 | // either with reconstructed values.... | |
217 | Double_t pt = vtx->Pt(); | |
218 | Double_t rapidity = vtx->YD0(); | |
219 | ||
220 | Double_t cosThetaStar = 9999.; | |
221 | Double_t pTpi = 0.; | |
222 | Double_t pTK = 0.; | |
223 | Int_t pdgCode = mcVtxHF->GetPdgCode(); | |
224 | if (pdgCode > 0){ | |
225 | cosThetaStar = vtx->CosThetaStarD0(); | |
226 | pTpi = vtx->PtProng(0); | |
227 | pTK = vtx->PtProng(1); | |
228 | } | |
229 | else { | |
230 | cosThetaStar = vtx->CosThetaStarD0bar(); | |
231 | pTpi = vtx->PtProng(1); | |
232 | pTK = vtx->PtProng(0); | |
233 | } | |
234 | ||
235 | Double_t cT = vtx->CtD0(); | |
236 | AliDebug(2, Form("cT from reconstructed vertex = %f (micron)",cT*1E4)); | |
237 | AliDebug(2, Form("pdg code from MC = %d",TMath::Abs(mcVtxHF->GetPdgCode()))); | |
238 | ||
239 | // ... or with generated values | |
240 | ||
241 | if (!fFillFromGenerated){ | |
242 | containerInput[0] = pt; | |
243 | containerInput[1] = rapidity; | |
244 | containerInput[2] = cosThetaStar; | |
245 | containerInput[3] = pTpi; | |
246 | containerInput[4] = pTK; | |
247 | containerInput[5] = cT*1.E4; // in micron | |
248 | } | |
249 | else { | |
250 | Double_t vectorMC[6] = {9999.,9999.,9999.,9999.,9999.,9999.}; | |
251 | if (GetGeneratedValuesFromMCParticle(mcVtxHF, mcArray, vectorMC)){ | |
252 | containerInput[0] = vectorMC[0]; | |
253 | containerInput[1] = vectorMC[1] ; | |
254 | containerInput[2] = vectorMC[2] ; | |
255 | containerInput[3] = vectorMC[3] ; | |
256 | containerInput[4] = vectorMC[4] ; | |
257 | containerInput[5] = vectorMC[5] ; // in micron | |
258 | } | |
259 | else { | |
260 | AliDebug(3,"Problems in filling the container"); | |
261 | continue; | |
262 | } | |
263 | } | |
264 | AliDebug(2, Form("Filling the container with pt = %f, rapidity = %f, cosThetaStar = %f, pTpi = %f, pTK = %f, cT = %f", containerInput[0], containerInput[1], containerInput[2], containerInput[3], containerInput[4], containerInput[5])); | |
265 | fCFManager->GetParticleContainer()->Fill(containerInput,1) ; | |
266 | } | |
267 | /* | |
268 | // for debugging only, filling container with dummy values | |
269 | Double_t pt = 1.5; | |
270 | Double_t rapidity = 0.5; | |
271 | Double_t cosThetaStar = 0.7; | |
272 | Double_t pTpi = 2.5; | |
273 | Double_t pTK = 1; | |
274 | Double_t cT = 3; | |
275 | containerInput[0] = pt; | |
276 | containerInput[1] = rapidity; | |
277 | containerInput[2] = cosThetaStar; | |
278 | containerInput[3] = pTpi; | |
279 | containerInput[4] = pTK; | |
280 | containerInput[5] = cT; | |
281 | fCFManager->GetParticleContainer()->Fill(containerInput,1) ; | |
282 | */ | |
283 | ||
284 | } | |
285 | ||
286 | fHistEventsProcessed->Fill(0); | |
287 | /* PostData(0) is taken care of by AliAnalysisTaskSE */ | |
288 | PostData(1,fHistEventsProcessed) ; | |
289 | PostData(2,fCFManager->GetParticleContainer()) ; | |
290 | } | |
291 | ||
292 | ||
293 | //___________________________________________________________________________ | |
294 | void AliCFHeavyFlavourTaskMultiVar::Terminate(Option_t*) | |
295 | { | |
296 | // The Terminate() function is the last function to be called during | |
297 | // a query. It always runs on the client, it can be used to present | |
298 | // the results graphically or save the results to file. | |
299 | ||
300 | Info("Terminate",""); | |
301 | AliAnalysisTaskSE::Terminate(); | |
302 | ||
303 | AliInfo(Form("Found %i MC particles that are D0 in MC in %d events",fCountMC,fEvents)); | |
304 | AliInfo(Form("Found %i reco D0 that are not decaying in K+pi in %d events",fSkipped,fEvents)); | |
305 | ||
306 | // draw some example plots.... | |
307 | ||
308 | AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2)); | |
309 | ||
310 | // projecting the containers to obtain histograms | |
311 | // first argument = variable, second argument = step | |
312 | TH1D* h00 = cont->ShowProjection(0,0) ; // pt | |
313 | TH1D* h10 = cont->ShowProjection(1,0) ; // rapidity | |
314 | TH1D* h20 = cont->ShowProjection(2,0) ; // cosThetaStar | |
315 | TH1D* h30 = cont->ShowProjection(3,0) ; // pTpi | |
316 | TH1D* h40 = cont->ShowProjection(4,0) ; // pTK | |
317 | TH1D* h50 = cont->ShowProjection(5,0) ; // cT | |
318 | ||
319 | TH1D* h01 = cont->ShowProjection(0,1) ; // pt | |
320 | TH1D* h11 = cont->ShowProjection(1,1) ; // rapidity | |
321 | TH1D* h21 = cont->ShowProjection(2,1) ; // cosThetaStar | |
322 | TH1D* h31 = cont->ShowProjection(3,1) ; // pTpi | |
323 | TH1D* h41 = cont->ShowProjection(4,1) ; // pTK | |
324 | TH1D* h51 = cont->ShowProjection(5,1) ; // cT | |
325 | ||
326 | h00->SetTitle("pT_D0 (GeV/c)"); | |
327 | h10->SetTitle("rapidity"); | |
328 | h20->SetTitle("cosThetaStar"); | |
329 | h30->SetTitle("pT_pi (GeV/c)"); | |
330 | h40->SetTitle("pT_K (Gev/c)"); | |
331 | h50->SetTitle("cT (#mum)"); | |
332 | ||
333 | h00->GetXaxis()->SetTitle("pT_D0 (GeV/c)"); | |
334 | h10->GetXaxis()->SetTitle("rapidity"); | |
335 | h20->GetXaxis()->SetTitle("cosThetaStar"); | |
336 | h30->GetXaxis()->SetTitle("pT_pi (GeV/c)"); | |
337 | h40->GetXaxis()->SetTitle("pT_K (Gev/c)"); | |
338 | h50->GetXaxis()->SetTitle("cT (#mum)"); | |
339 | ||
340 | h01->SetTitle("pT_D0 (GeV/c)"); | |
341 | h11->SetTitle("rapidity"); | |
342 | h21->SetTitle("cosThetaStar"); | |
343 | h31->SetTitle("pT_pi (GeV/c)"); | |
344 | h41->SetTitle("pT_K (Gev/c)"); | |
345 | h51->SetTitle("cT (#mum)"); | |
346 | ||
347 | h01->GetXaxis()->SetTitle("pT_D0 (GeV/c)"); | |
348 | h11->GetXaxis()->SetTitle("rapidity"); | |
349 | h21->GetXaxis()->SetTitle("cosThetaStar"); | |
350 | h31->GetXaxis()->SetTitle("pT_pi (GeV/c)"); | |
351 | h41->GetXaxis()->SetTitle("pT_K (Gev/c)"); | |
352 | h51->GetXaxis()->SetTitle("cT (#mum)"); | |
353 | ||
354 | Double_t max0 = h00->GetMaximum(); | |
355 | Double_t max1 = h10->GetMaximum(); | |
356 | Double_t max2 = h20->GetMaximum(); | |
357 | Double_t max3 = h30->GetMaximum(); | |
358 | Double_t max4 = h40->GetMaximum(); | |
359 | Double_t max5 = h50->GetMaximum(); | |
360 | ||
361 | h00->GetYaxis()->SetRangeUser(0,max0*1.2); | |
362 | h10->GetYaxis()->SetRangeUser(0,max1*1.2); | |
363 | h20->GetYaxis()->SetRangeUser(0,max2*1.2); | |
364 | h30->GetYaxis()->SetRangeUser(0,max3*1.2); | |
365 | h40->GetYaxis()->SetRangeUser(0,max4*1.2); | |
366 | h50->GetYaxis()->SetRangeUser(0,max5*1.2); | |
367 | ||
368 | h01->GetYaxis()->SetRangeUser(0,max0*1.2); | |
369 | h11->GetYaxis()->SetRangeUser(0,max1*1.2); | |
370 | h21->GetYaxis()->SetRangeUser(0,max2*1.2); | |
371 | h31->GetYaxis()->SetRangeUser(0,max3*1.2); | |
372 | h41->GetYaxis()->SetRangeUser(0,max4*1.2); | |
373 | h51->GetYaxis()->SetRangeUser(0,max5*1.2); | |
374 | ||
375 | h00->SetMarkerStyle(20); | |
376 | h10->SetMarkerStyle(24); | |
377 | h20->SetMarkerStyle(21); | |
378 | h30->SetMarkerStyle(25); | |
379 | h40->SetMarkerStyle(27); | |
380 | h50->SetMarkerStyle(28); | |
381 | ||
382 | h00->SetMarkerColor(2); | |
383 | h10->SetMarkerColor(2); | |
384 | h20->SetMarkerColor(2); | |
385 | h30->SetMarkerColor(2); | |
386 | h40->SetMarkerColor(2); | |
387 | h50->SetMarkerColor(2); | |
388 | ||
389 | h01->SetMarkerStyle(20) ; | |
390 | h11->SetMarkerStyle(24) ; | |
391 | h21->SetMarkerStyle(21) ; | |
392 | h31->SetMarkerStyle(25) ; | |
393 | h41->SetMarkerStyle(27) ; | |
394 | h51->SetMarkerStyle(28) ; | |
395 | ||
396 | h01->SetMarkerColor(4); | |
397 | h11->SetMarkerColor(4); | |
398 | h21->SetMarkerColor(4); | |
399 | h31->SetMarkerColor(4); | |
400 | h41->SetMarkerColor(4); | |
401 | h51->SetMarkerColor(4); | |
402 | ||
403 | gStyle->SetCanvasColor(0); | |
404 | gStyle->SetFrameFillColor(0); | |
405 | gStyle->SetTitleFillColor(0); | |
406 | gStyle->SetStatColor(0); | |
407 | ||
408 | // drawing in 2 separate canvas for a matter of clearity | |
409 | TCanvas * c1 =new TCanvas("c1","pT, rapidiy, cosThetaStar",800,1600); | |
410 | TCanvas * c2 =new TCanvas("c2","pTpi, pTK, cT",800,1600); | |
411 | c1->Divide(2,3); | |
412 | c2->Divide(2,3); | |
413 | ||
414 | c1->cd(1); | |
415 | h00->Draw("p"); | |
416 | c1->cd(2); | |
417 | h01->Draw("p"); | |
418 | c1->cd(3); | |
419 | h10->Draw("p"); | |
420 | c1->cd(4); | |
421 | h11->Draw("p"); | |
422 | c1->cd(5); | |
423 | h20->Draw("p"); | |
424 | c1->cd(6); | |
425 | h21->Draw("p"); | |
426 | c1->cd(); | |
427 | ||
428 | c2->cd(1); | |
429 | h30->Draw("p"); | |
430 | c2->cd(2); | |
431 | h31->Draw("p"); | |
432 | c2->cd(3); | |
433 | h40->Draw("p"); | |
434 | c2->cd(4); | |
435 | h41->Draw("p"); | |
436 | c2->cd(5); | |
437 | h50->Draw("p"); | |
438 | c2->cd(6); | |
439 | h51->Draw("p"); | |
440 | c2->cd(); | |
441 | ||
442 | c1->SaveAs("Variables/pT_rapidity_cosThetaStar.eps"); | |
443 | c2->SaveAs("Variables/pTpi_pTK_cT.eps"); | |
444 | ||
445 | } | |
446 | ||
447 | //___________________________________________________________________________ | |
448 | void AliCFHeavyFlavourTaskMultiVar::UserCreateOutputObjects() { | |
449 | //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED | |
450 | //TO BE SET BEFORE THE EXECUTION OF THE TASK | |
451 | // | |
452 | Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName()); | |
453 | ||
454 | //slot #1 | |
455 | OpenFile(1); | |
456 | fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ; | |
457 | } | |
458 | ||
459 | //___________________________________________________________________________ | |
460 | Int_t AliCFHeavyFlavourTaskMultiVar::GetVtxLabel(Int_t* labels) { | |
461 | // | |
462 | // returns the label of the vertex, given the labels of the 2 daughter tracks | |
463 | // returns -1 if the V0 is fake | |
464 | // | |
465 | ||
466 | Int_t labD0daugh0=-1; | |
467 | Int_t labD0daugh1=-1; | |
468 | Int_t labD0=-1; | |
469 | Int_t labMother = -1; | |
470 | Int_t pdgMother = -1; | |
471 | ||
472 | AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent); | |
473 | TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName())); | |
474 | if (!mcArray) AliError("Could not find Monte-Carlo in AOD"); | |
475 | ||
476 | AliAODMCParticle *part0 = (AliAODMCParticle*)mcArray->At(labels[0]); | |
477 | if(!part0) { | |
478 | printf("no MC particle\n"); | |
479 | return -1; | |
480 | } | |
481 | ||
482 | while(part0->GetMother()>=0) { | |
483 | labMother=part0->GetMother(); | |
484 | part0 = (AliAODMCParticle*)mcArray->At(labMother); | |
485 | if(!part0) { | |
486 | printf("no MC mother particle\n"); | |
487 | } | |
488 | pdgMother = TMath::Abs(part0->GetPdgCode()); | |
489 | AliDebug(2,Form("pdgMother from particle 0 = %d",pdgMother)); | |
490 | if(pdgMother==421) { | |
491 | AliDebug(2,Form("found one particle (first daughter) coming from a D0")); | |
492 | labD0daugh0=labMother; | |
493 | } | |
494 | } | |
495 | ||
496 | AliAODMCParticle *part1 = (AliAODMCParticle*)mcArray->At(labels[1]); | |
497 | if(!part1) { | |
498 | printf("no MC particle\n"); | |
499 | return -1; | |
500 | } | |
501 | ||
502 | while(part1->GetMother()>=0) { | |
503 | labMother=part1->GetMother(); | |
504 | part1 = (AliAODMCParticle*)mcArray->At(labMother); | |
505 | if(!part1) { | |
506 | printf("no MC mother particle\n"); | |
507 | } | |
508 | pdgMother = TMath::Abs(part1->GetPdgCode()); | |
509 | AliDebug(2,Form("pdgMother from particle 1 = %d",pdgMother)); | |
510 | if(pdgMother==421) { | |
511 | AliDebug(2,Form("found one particle (second daughter) coming from a D0")); | |
512 | labD0daugh1=labMother; | |
513 | } | |
514 | } | |
515 | if ((labD0daugh0 >= 0) && (labD0daugh1 >= 0) && (labD0daugh0 == labD0daugh1)){ | |
516 | //AliDebug(2, Form("found particles from the same D0 mother")); | |
517 | labD0 = labD0daugh0; | |
518 | // check that the D0 decays in 2 prongs | |
519 | AliAODMCParticle *part = (AliAODMCParticle*)mcArray->At(labD0); | |
520 | if (TMath::Abs(part->GetDaughter(1)-part->GetDaughter(0))==1) { | |
521 | return labD0; | |
522 | } | |
523 | else { | |
524 | fSkipped+=1; | |
525 | } | |
526 | } | |
527 | ||
528 | return -1; | |
529 | ||
530 | } | |
531 | ||
532 | //___________________________________________________________________________ | |
533 | Int_t AliCFHeavyFlavourTaskMultiVar::IsMcVtx(AliAODRecoDecayHF2Prong* vtx) { | |
534 | // | |
535 | // check if the passed vertex is associated to a MC one, | |
536 | // and returns the corresponding geant label. | |
537 | // returns -1 if it is fake (i.e. label<0). | |
538 | // | |
539 | ||
540 | ||
541 | AliAODTrack *trk0 = (AliAODTrack*)vtx->GetDaughter(0); | |
542 | AliAODTrack *trk1 = (AliAODTrack*)vtx->GetDaughter(1); | |
543 | ||
544 | if (!trk0 || ! trk1) { | |
545 | AliDebug(2, "problems with the tracks while looking for mother"); | |
546 | //if (!trk0) AliDebug(2, "track 0 gives problems"); | |
547 | //if (!trk1) AliDebug(2, "track 1 gives problems"); | |
548 | return -1; | |
549 | } | |
550 | ||
551 | Int_t labels[2] = {-1,-1}; | |
552 | labels[0] = trk0->GetLabel(); | |
553 | labels[1] = trk1->GetLabel(); | |
554 | ||
555 | if (labels[0] < 0 || labels[1] < 0) { | |
556 | AliDebug(2, "problems with the labels of the tracks while looking for mother"); | |
557 | return -1; | |
558 | } | |
559 | ||
560 | return GetVtxLabel(&labels[0]) ; | |
561 | } | |
562 | //___________________________________________________________________________ | |
563 | Double_t AliCFHeavyFlavourTaskMultiVar::CosThetaStar(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const { | |
564 | ||
565 | // | |
566 | // to calculate cos(ThetaStar) for generated particle | |
567 | // using the K, since mcPartDaughter0 and mcPartDaughter1 always correspond to K and pi respectively | |
568 | // (see where the function is called) | |
569 | // | |
570 | ||
571 | Int_t pdgvtx = mcPart->GetPdgCode(); | |
572 | /* if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar) | |
573 | Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode()); | |
574 | Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode()); | |
575 | AliInfo(Form("D0, with pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1)); | |
576 | AliDebug(2,"This is a D0"); | |
577 | AliAODMCParticle* mcPartdummy = mcPartDaughter0; | |
578 | mcPartDaughter0 = mcPartDaughter1; | |
579 | mcPartDaughter1 = mcPartdummy; | |
580 | } | |
581 | else{ | |
582 | AliInfo("D0bar"); | |
583 | } | |
584 | */ | |
585 | Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode()); | |
586 | Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode()); | |
587 | if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar) | |
588 | AliDebug(2,"D0"); | |
589 | } | |
590 | else{ | |
591 | AliDebug(2,"D0bar"); | |
592 | } | |
593 | if (pdgprong0 == 211){ | |
594 | AliDebug(2,Form("pdgprong0 = %d, pdgprong1 = %d, switching...",pdgprong0,pdgprong1)); | |
595 | AliAODMCParticle* mcPartdummy = mcPartDaughter0; | |
596 | mcPartDaughter0 = mcPartDaughter1; | |
597 | mcPartDaughter1 = mcPartdummy; | |
598 | pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode()); | |
599 | pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode()); | |
600 | } | |
601 | ||
602 | AliDebug(2,Form("After checking, pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1)); | |
603 | Double_t massvtx = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdgvtx))->Mass(); | |
604 | Double_t massp[2]; | |
605 | massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass(); | |
606 | massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass(); | |
607 | ||
608 | Double_t pStar = TMath::Sqrt(TMath::Power(massvtx*massvtx-massp[0]*massp[0]-massp[1]*massp[1],2.)-4.*massp[0]*massp[0]*massp[1]*massp[1])/(2.*massvtx); | |
609 | ||
610 | Double_t px = mcPartDaughter0->Px()+mcPartDaughter1->Px(); | |
611 | Double_t py = mcPartDaughter0->Py()+mcPartDaughter1->Py(); | |
612 | Double_t pz = mcPartDaughter0->Pz()+mcPartDaughter1->Pz(); | |
613 | Double_t p = TMath::Sqrt(px*px+py*py+pz*pz); | |
614 | Double_t e = TMath::Sqrt(massvtx*massvtx+p*p); | |
615 | Double_t beta = p/e; | |
616 | Double_t gamma = e/massvtx; | |
617 | TVector3 mom(mcPartDaughter0->Px(),mcPartDaughter0->Py(),mcPartDaughter0->Pz()); | |
618 | TVector3 momTot(mcPartDaughter0->Px()+mcPartDaughter1->Px(),mcPartDaughter0->Py()+mcPartDaughter1->Py(),mcPartDaughter0->Pz()+mcPartDaughter1->Pz()); | |
619 | ||
620 | Double_t qlprong = mom.Dot(momTot)/momTot.Mag(); // analog to AliAODRecoDecay::QlProng(0) | |
621 | ||
622 | AliDebug(2,Form("pStar = %f, beta = %f, gamma = %f, qlprong = %f, massp[0] = %f", pStar, beta, gamma, qlprong, massp[0])); | |
623 | Double_t cts = (qlprong/gamma-beta*TMath::Sqrt(pStar*pStar+massp[0]*massp[0]))/pStar; | |
624 | AliDebug(2,Form("cts = %f", cts)); | |
625 | return cts; | |
626 | } | |
627 | //___________________________________________________________________________ | |
628 | Double_t AliCFHeavyFlavourTaskMultiVar::CT(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const { | |
629 | ||
630 | // | |
631 | // to calculate cT for generated particle | |
632 | // | |
633 | ||
634 | Double_t xmcPart[3] = {0,0,0}; | |
635 | Double_t xdaughter0[3] = {0,0,0}; | |
636 | Double_t xdaughter1[3] = {0,0,0}; | |
637 | mcPart->XvYvZv(xmcPart); // cm | |
638 | mcPartDaughter0->XvYvZv(xdaughter0); // cm | |
639 | mcPartDaughter1->XvYvZv(xdaughter1); //cm | |
640 | Double_t prodVtxD0 = TMath::Sqrt(xmcPart[0]*xmcPart[0]+ | |
641 | xmcPart[1]*xmcPart[1]+ | |
642 | xmcPart[2]*xmcPart[2]); | |
643 | Double_t prodVtxDaughter0 = TMath::Sqrt(xdaughter0[0]*xdaughter0[0]+ | |
644 | xdaughter0[1]*xdaughter0[1]+ | |
645 | xdaughter0[2]*xdaughter0[2]); | |
646 | Double_t prodVtxDaughter1 = TMath::Sqrt(xdaughter1[0]*xdaughter1[0]+ | |
647 | xdaughter1[1]*xdaughter1[1]+ | |
648 | xdaughter1[2]*xdaughter1[2]); | |
649 | ||
650 | AliDebug(2, Form("D0: x = %f, y = %f, z = %f, production vertex distance = %f (cm), %f (micron)", xmcPart[0], xmcPart[1], xmcPart[2], prodVtxD0, prodVtxD0*1.E4)); | |
651 | AliDebug(2, Form("Daughter0: x = %f, y = %f, z = %f, production vertex distance = %f (cm) %f (micron)", xdaughter0[0], xdaughter0[1], xdaughter0[2], prodVtxDaughter0, prodVtxDaughter0*1E4)); | |
652 | AliDebug(2, Form("Daughter1: x = %f, y = %f, z = %f, production vertex distance = %f (cm) %f (micron)", xdaughter1[0], xdaughter1[1], xdaughter1[2], prodVtxDaughter1, prodVtxDaughter1*1.E4)); | |
653 | ||
654 | Double_t cT0 = TMath::Sqrt((xdaughter0[0]-xmcPart[0])*(xdaughter0[0]-xmcPart[0])+ | |
655 | (xdaughter0[1]-xmcPart[1])*(xdaughter0[1]-xmcPart[1])+ | |
656 | (xdaughter0[2]-xmcPart[2])*(xdaughter0[2]-xmcPart[2])); | |
657 | ||
658 | Double_t cT1 = TMath::Sqrt((xdaughter1[0]-xmcPart[0])*(xdaughter1[0]-xmcPart[0])+ | |
659 | (xdaughter1[1]-xmcPart[1])*(xdaughter1[1]-xmcPart[1])+ | |
660 | (xdaughter1[2]-xmcPart[2])*(xdaughter1[2]-xmcPart[2])); | |
661 | ||
662 | if (cT0 != cT1) { | |
663 | AliWarning("cT from daughter 0 different from cT from daughter 1! Using cT from daughter 0, but PLEASE, CHECK...."); | |
664 | } | |
665 | ||
666 | // calculating cT from cT0 | |
667 | ||
668 | Double_t mass = TDatabasePDG::Instance()->GetParticle(mcPart->GetPdgCode())->Mass(); // mcPart->GetPdgCode() should return +/- 421 for the D0/D0bar | |
669 | Double_t p = mcPart-> P(); | |
670 | Double_t cT = cT0*mass/p; | |
671 | AliDebug(2, Form("cT from daughter 0 = %f (micron)", cT0*1E4)); | |
672 | AliDebug(2, Form("cT from daughter 1 = %f (micron)", cT1*1E4)); | |
673 | AliDebug(2, Form("cT (with mass = %f and p = %f) = %f (micron)", mass, p, cT*1E4)); | |
674 | return cT; | |
675 | } | |
676 | //________________________________________________________________________________ | |
677 | Bool_t AliCFHeavyFlavourTaskMultiVar::GetGeneratedValuesFromMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC) const { | |
678 | ||
679 | // | |
680 | // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle | |
681 | // | |
682 | ||
683 | Bool_t isOk = kFALSE; | |
684 | ||
685 | // check whether the D0 decays in pi+K | |
686 | // to be added!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! | |
687 | // could use a cut in the CF, but not clear how to define a TDecayChannel | |
688 | // implemented for the time being as a cut on the number of daughters - see later when | |
689 | // getting the daughters | |
690 | ||
691 | // getting the daughters | |
692 | Int_t daughter0 = mcPart->GetDaughter(0); | |
693 | Int_t daughter1 = mcPart->GetDaughter(1); | |
694 | AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1)); | |
695 | if (daughter0 == 0 || daughter1 == 0) { | |
696 | AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!"); | |
697 | return isOk; | |
698 | } | |
699 | if (TMath::Abs(daughter1 - daughter0 != 1)) { | |
700 | AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!"); | |
701 | return isOk; | |
702 | } | |
703 | AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0)); | |
704 | AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1)); | |
705 | if (!mcPartDaughter0 || !mcPartDaughter1) { | |
706 | AliWarning("At least one Daughter Particle not found in tree, skipping"); | |
707 | return isOk; | |
708 | } | |
709 | ||
710 | Double_t vtx1[3] = {0,0,0}; // primary vertex | |
711 | Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0 | |
712 | Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1 | |
713 | mcPart->XvYvZv(vtx1); // cm | |
714 | // getting vertex from daughters | |
715 | mcPartDaughter0->XvYvZv(vtx2daughter0); // cm | |
716 | mcPartDaughter1->XvYvZv(vtx2daughter1); //cm | |
717 | if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) { | |
718 | AliError("Daughters have different secondary vertex, skipping the track"); | |
719 | return isOk; | |
720 | } | |
721 | Int_t nprongs = 2; | |
722 | Short_t charge = 0; | |
723 | // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second | |
724 | AliAODMCParticle* positiveDaugh = mcPartDaughter0; | |
725 | AliAODMCParticle* negativeDaugh = mcPartDaughter1; | |
726 | if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){ | |
727 | // inverting in case the positive daughter is the second one | |
728 | positiveDaugh = mcPartDaughter1; | |
729 | negativeDaugh = mcPartDaughter0; | |
730 | } | |
731 | // getting the momentum from the daughters | |
732 | Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()}; | |
733 | Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()}; | |
734 | Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()}; | |
735 | ||
736 | Double_t d0[2] = {0.,0.}; | |
737 | ||
738 | AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0); | |
739 | ||
740 | Double_t cosThetaStar = 0.; | |
741 | Double_t cosThetaStarD0 = 0.; | |
742 | Double_t cosThetaStarD0bar = 0.; | |
743 | cosThetaStarD0 = decay->CosThetaStar(1,421,211,321); | |
744 | cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211); | |
745 | if (mcPart->GetPdgCode() == 421){ // D0 | |
746 | AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode())); | |
747 | cosThetaStar = cosThetaStarD0; | |
748 | } | |
749 | else if (mcPart->GetPdgCode() == -421){ // D0bar{ | |
750 | AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode())); | |
751 | cosThetaStar = cosThetaStarD0bar; | |
752 | } | |
753 | else{ | |
754 | AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check..."); | |
755 | return vectorMC; | |
756 | } | |
757 | if (cosThetaStar < -1 || cosThetaStar > 1) { | |
758 | AliWarning("Invalid value for cosine Theta star, returning"); | |
759 | return isOk; | |
760 | } | |
761 | ||
762 | // calculate cos(Theta*) according to the method implemented herein | |
763 | ||
764 | Double_t cts = 9999.; | |
765 | cts = CosThetaStar(mcPart, mcPartDaughter0, mcPartDaughter1); | |
766 | if (cts < -1 || cts > 1) { | |
767 | AliWarning("Invalid value for cosine Theta star from AliCFHeavyFlavourTaskMultiVar method"); | |
768 | } | |
769 | if (TMath::Abs(cts - cosThetaStar)>0.001) { | |
770 | AliError(Form("cosThetaStar automatically calculated different from that manually calculated!!! cosThetaStar = %f, cosThetaStar = %f", cosThetaStar,cts)); | |
771 | } | |
772 | ||
773 | Double_t cT = decay->Ct(421); | |
774 | ||
775 | // calculate cT from the method implemented herein | |
776 | Double_t cT1 = 0.; | |
777 | cT1 = CT(mcPart, mcPartDaughter0, mcPartDaughter1); | |
778 | ||
779 | if (TMath::Abs(cT1 - cT)>0.001) { | |
780 | AliError(Form("cT automatically calculated different from that manually calculated!!! cT = %f, cT1 = %f",cT,cT1)); | |
781 | } | |
782 | ||
783 | // get the pT of the daughters | |
784 | ||
785 | Double_t pTpi = 0.; | |
786 | Double_t pTK = 0.; | |
787 | ||
788 | if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) { | |
789 | pTpi = mcPartDaughter0->Pt(); | |
790 | pTK = mcPartDaughter1->Pt(); | |
791 | } | |
792 | else { | |
793 | pTpi = mcPartDaughter1->Pt(); | |
794 | pTK = mcPartDaughter0->Pt(); | |
795 | } | |
796 | ||
797 | vectorMC[0] = mcPart->Pt(); | |
798 | vectorMC[1] = mcPart->Y() ; | |
799 | vectorMC[2] = cosThetaStar ; | |
800 | vectorMC[3] = pTpi ; | |
801 | vectorMC[4] = pTK ; | |
802 | vectorMC[5] = cT*1.E4 ; // in micron | |
803 | isOk = kTRUE; | |
804 | return isOk; | |
805 | } | |
806 |