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