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