git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New class for corrections in steps and macro to plot efficiency for different cut...
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliCFHeavyFlavourTaskMultiVarMultiStep.cxx
CommitLineData
f5ec6c30 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// 6 Steps introduced: MC, MC Acc, Reco, Reco Acc, Reco Acc + ITS Cl,
19// Reco Acc + ITS Cl + PPR cuts
20// 11 variables used: pt, y, cosThetaStar, ptPi, ptK, ct,
21// dca, d0Pi, d0K, d0Pixd0K, cosPointingAngle
22//
23//-----------------------------------------------------------------------
24// Author : C. Zampolli, CERN
25//-----------------------------------------------------------------------
26
27#include <TCanvas.h>
28#include <TParticle.h>
29#include <TH1I.h>
30#include <TStyle.h>
31
32#include "AliCFHeavyFlavourTaskMultiVarMultiStep.h"
33#include "AliStack.h"
34#include "AliMCEvent.h"
35#include "AliCFManager.h"
36#include "AliCFContainer.h"
37#include "AliLog.h"
38#include "AliAODEvent.h"
39#include "AliAODRecoDecay.h"
40#include "AliAODRecoDecayHF.h"
41#include "AliAODRecoDecayHF2Prong.h"
42#include "AliAODMCParticle.h"
43
44//__________________________________________________________________________
45AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep() :
46 AliAnalysisTaskSE(),
47 fPDG(0),
48 fCFManager(0x0),
49 fHistEventsProcessed(0x0),
50 fCountMC(0),
51 fCountAcc(0),
52 fCountReco(0),
53 fCountRecoAcc(0),
54 fCountRecoITSClusters(0),
55 fCountRecoPPR(0),
56 fEvents(0),
57 fFillFromGenerated(kFALSE),
58 fMinITSClusters(5)
59{
60 //
61 //Default ctor
62 //
63}
64//___________________________________________________________________________
65AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep(const Char_t* name) :
66 AliAnalysisTaskSE(name),
67 fPDG(0),
68 fCFManager(0x0),
69 fHistEventsProcessed(0x0),
70 fCountMC(0),
71 fCountAcc(0),
72 fCountReco(0),
73 fCountRecoAcc(0),
74 fCountRecoITSClusters(0),
75 fCountRecoPPR(0),
76 fEvents(0),
77 fFillFromGenerated(kFALSE),
78 fMinITSClusters(5)
79{
80 //
81 // Constructor. Initialization of Inputs and Outputs
82 //
83 Info("AliCFHeavyFlavourTaskMultiVarMultiStep","Calling Constructor");
84 /*
85 DefineInput(0) and DefineOutput(0)
86 are taken care of by AliAnalysisTaskSE constructor
87 */
88 DefineOutput(1,TH1I::Class());
89 DefineOutput(2,AliCFContainer::Class());
90}
91
92//___________________________________________________________________________
93AliCFHeavyFlavourTaskMultiVarMultiStep& AliCFHeavyFlavourTaskMultiVarMultiStep::operator=(const AliCFHeavyFlavourTaskMultiVarMultiStep& c)
94{
95 //
96 // Assignment operator
97 //
98 if (this!=&c) {
99 AliAnalysisTaskSE::operator=(c) ;
100 fPDG = c.fPDG;
101 fCFManager = c.fCFManager;
102 fHistEventsProcessed = c.fHistEventsProcessed;
103 }
104 return *this;
105}
106
107//___________________________________________________________________________
108AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep(const AliCFHeavyFlavourTaskMultiVarMultiStep& c) :
109 AliAnalysisTaskSE(c),
110 fPDG(c.fPDG),
111 fCFManager(c.fCFManager),
112 fHistEventsProcessed(c.fHistEventsProcessed),
113 fCountMC(c.fCountMC),
114 fCountAcc(c.fCountAcc),
115 fCountReco(c.fCountReco),
116 fCountRecoAcc(c.fCountRecoAcc),
117 fCountRecoITSClusters(c.fCountRecoITSClusters),
118 fCountRecoPPR(c.fCountRecoPPR),
119 fEvents(c.fEvents),
120 fFillFromGenerated(c.fFillFromGenerated),
121 fMinITSClusters(c.fMinITSClusters)
122{
123 //
124 // Copy Constructor
125 //
126}
127
128//___________________________________________________________________________
129AliCFHeavyFlavourTaskMultiVarMultiStep::~AliCFHeavyFlavourTaskMultiVarMultiStep() {
130 //
131 //destructor
132 //
133 if (fCFManager) delete fCFManager ;
134 if (fHistEventsProcessed) delete fHistEventsProcessed ;
135}
136
137//_________________________________________________
138void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
139{
140 //
141 // Main loop function
142 //
143
144 if (fFillFromGenerated){
145 AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
146 }
147
148 if (!fInputEvent) {
149 Error("UserExec","NO EVENT FOUND!");
150 return;
151 }
152
153 fEvents++;
154 if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
155 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
156 fCFManager->SetEventInfo(aodEvent);
157
158 // MC-event selection
159 Double_t containerInput[11] ;
160
161 //loop on the MC event
162
163 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
164 if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
165 Int_t icountMC = 0;
166 Int_t icountAcc = 0;
167 Int_t icountReco = 0;
168 Int_t icountRecoAcc = 0;
169 Int_t icountRecoITSClusters = 0;
170 Int_t icountRecoPPR = 0;
171
172 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
173 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
174 if (!mcPart) {
175 AliWarning("Particle not found in tree, skipping");
176 continue;
177 }
178
179 // check the MC-level cuts
180 if (!fCFManager->CheckParticleCuts(0, mcPart)) continue; // 0 stands for MC level
181
182 // fill the container for Gen-level selection
183 Double_t vectorMC[6] = {9999.,9999.,9999.,9999.,9999.,9999.};
184 if (GetGeneratedValuesFromMCParticle(mcPart, mcArray, vectorMC)){
185 containerInput[0] = vectorMC[0];
186 containerInput[1] = vectorMC[1] ;
187 containerInput[2] = vectorMC[2] ;
188 containerInput[3] = vectorMC[3] ;
189 containerInput[4] = vectorMC[4] ;
190 containerInput[5] = vectorMC[5] ; // in micron
191 containerInput[6] = 0.; // dummy value, meaningless in MC, in micron
192 containerInput[7] = 0.; // dummy value, meaningless in MC, in micron
193 containerInput[8] = 0.; // dummy value, meaningless in MC, in micron
194 containerInput[9] = -100000.; // dummy value, meaningless in MC, in micron^2
195 containerInput[10] = 1.01; // dummy value, meaningless in MC
196 fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);
197 icountMC++;
198
199 // check the MC-Acceptance level cuts
200 // since standard CF functions are not applicable, using Kine Cuts on daughters
201
202 Int_t daughter0 = mcPart->GetDaughter(0);
203 Int_t daughter1 = mcPart->GetDaughter(1);
204 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
205 if (daughter0 == 0 || daughter1 == 0) {
206 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!! But it should have, this check was already done...");
207 }
208 if (TMath::Abs(daughter1 - daughter0 != 1)) {
209 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, but it should be, this check was already done...");
210 }
211 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
212 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
213 if (!mcPartDaughter0 || !mcPartDaughter1) {
214 AliWarning("At least one Daughter Particle not found in tree, but it should be, this check was already done...");
215 }
216 Double_t eta0 = mcPartDaughter0->Eta();
217 Double_t eta1 = mcPartDaughter1->Eta();
218 Double_t y0 = mcPartDaughter0->Y();
219 Double_t y1 = mcPartDaughter1->Y();
220 Double_t pt0 = mcPartDaughter0->Pt();
221 Double_t pt1 = mcPartDaughter1->Pt();
222 AliDebug(2, Form("Daughter 0: eta = %f, y = %f, pt = %f", eta0, y0, pt0));
223 AliDebug(2, Form("Daughter 1: eta = %f, y = %f, pt = %f", eta1, y1, pt1));
224 Bool_t daught0inAcceptance = (TMath::Abs(eta0) < 0.9 && pt0 > 0.1);
225 Bool_t daught1inAcceptance = (TMath::Abs(eta1) < 0.9 && pt1 > 0.1);
226 if (daught0inAcceptance && daught1inAcceptance) {
227 // checking whether the cuts implemented in the CF are equivalent - simply a cross-check
228 AliDebug(2, "Daughter particles in acceptance");
229 if (!fCFManager->CheckParticleCuts(1, mcPartDaughter0)) {
230 AliInfo("Inconsistency with CF cut for daughter 0!");
231 }
232 if (!fCFManager->CheckParticleCuts(1, mcPartDaughter1)) {
233 AliInfo("Inconsistency with CF cut for daughter 1!");
234 }
235 fCFManager->GetParticleContainer()->Fill(containerInput,kStepAcceptance);
236 icountAcc++;
237 }
238 }
239 else {
240 AliDebug(3,"Problems in filling the container");
241 continue;
242 }
243 }
244
245 AliDebug(2, Form("Found %i MC particles that are D0!!",icountMC));
246 AliDebug(2, Form("Found %i MC particles that are D0 and satisfy Acc cuts!!",icountAcc));
247
248 // Now go to rec level
249 fCountMC += icountMC;
250 fCountAcc += icountAcc;
251
252 // load heavy flavour vertices
253 TClonesArray *arrayVerticesHF = (TClonesArray*)((aodEvent->GetList())->FindObject("D0toKpi"));
254 if (!arrayVerticesHF) AliError("Could not find array of HF vertices");
255 AliDebug(2, Form("Found %d vertices",arrayVerticesHF->GetEntriesFast()));
256
257 for (Int_t iVertex = 0; iVertex<arrayVerticesHF->GetEntriesFast(); iVertex++) {
258
259 AliAODRecoDecayHF2Prong* vtx = (AliAODRecoDecayHF2Prong*)arrayVerticesHF->At(iVertex);
260
261 // find associated MC particle
262 Int_t mcLabel = vtx->MatchToMC(421, mcArray) ;
263 if (mcLabel == -1)
264 {
265 AliDebug(2,"No MC particle found");
266 }
267 else {
268 AliAODMCParticle* mcVtxHF = (AliAODMCParticle*)mcArray->At(mcLabel);
269 if (!mcVtxHF) {
270 AliWarning("Could not find associated MC in AOD MC tree");
271 continue;
272 }
273
274 // check if associated MC v0 passes the cuts
275 if (!fCFManager->CheckParticleCuts(0 ,mcVtxHF)) { // 0 stands for MC
276 AliDebug(2, "Skipping the particles due to cuts");
277 continue;
278 }
279
280
281 // fill the container...
282 Double_t pt = vtx->Pt();
283 Double_t rapidity = vtx->YD0();
284
285 Double_t cosThetaStar = 9999.;
286 Double_t pTpi = 0.;
287 Double_t pTK = 0.;
288 Double_t dca = vtx->GetDCA();
289 Double_t d0pi = 0.;
290 Double_t d0K = 0.;
291 Double_t d0xd0 = vtx->Prodd0d0();
292 Double_t cosPointingAngle = vtx->CosPointingAngle();
293 Int_t pdgCode = mcVtxHF->GetPdgCode();
294 if (pdgCode > 0){
295 cosThetaStar = vtx->CosThetaStarD0();
296 pTpi = vtx->PtProng(0);
297 pTK = vtx->PtProng(1);
298 d0pi = vtx->Getd0Prong(0);
299 d0K = vtx->Getd0Prong(1);
300 }
301 else {
302 cosThetaStar = vtx->CosThetaStarD0bar();
303 pTpi = vtx->PtProng(1);
304 pTK = vtx->PtProng(0);
305 d0pi = vtx->Getd0Prong(1);
306 d0K = vtx->Getd0Prong(0);
307 }
308
309 Double_t cT = vtx->CtD0();
310
311 if (!fFillFromGenerated){
312 // ...either with reconstructed values....
313 containerInput[0] = pt;
314 containerInput[1] = rapidity;
315 containerInput[2] = cosThetaStar;
316 containerInput[3] = pTpi;
317 containerInput[4] = pTK;
318 containerInput[5] = cT*1.E4; // in micron
319 containerInput[6] = dca*1.E4; // in micron
320 containerInput[7] = d0pi*1.E4; // in micron
321 containerInput[8] = d0K*1.E4; // in micron
322 containerInput[9] = d0xd0*1.E8; // in micron^2
323 containerInput[10] = cosPointingAngle; // in micron
324 }
325 else {
326 // ... or with generated values
327 Double_t vectorMC[6] = {9999.,9999.,9999.,9999.,9999.,9999.};
328 if (GetGeneratedValuesFromMCParticle(mcVtxHF, mcArray, vectorMC)){
329 containerInput[0] = vectorMC[0];
330 containerInput[1] = vectorMC[1] ;
331 containerInput[2] = vectorMC[2] ;
332 containerInput[3] = vectorMC[3] ;
333 containerInput[4] = vectorMC[4] ;
334 containerInput[5] = vectorMC[5] ; // in micron
335 containerInput[6] = 0.; // dummy value, meaningless in MC, in micron
336 containerInput[7] = 0.; // dummy value, meaningless in MC, in micron
337 containerInput[8] = 0.; // dummy value, meaningless in MC, in micron
338 containerInput[9] = 100000.; // dummy value, meaningless in MC, in micron^2
339 containerInput[10] = 1.01; // dummy value, meaningless in MC
340 }
341 else {
342 AliDebug(3,"Problems in filling the container");
343 continue;
344 }
345 }
346 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]));
347 icountReco++;
348 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
349
350 // cut in acceptance
351 Bool_t acceptanceProng0 = (TMath::Abs(vtx->EtaProng(0))< 0.9 && vtx->PtProng(0) > 0.1);
352 Bool_t acceptanceProng1 = (TMath::Abs(vtx->EtaProng(1))< 0.9 && vtx->PtProng(1) > 0.1);
353 if (acceptanceProng0 && acceptanceProng1) {
354 AliDebug(2,"D0 reco daughters in acceptance");
355 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance) ;
356 icountRecoAcc++;
357
358 // cut on the min n. of clusters in ITS
359 Int_t ncls0=0;
360 for(Int_t l=0;l<6;l++) if(TESTBIT(vtx->GetITSClusterMap(),l)) ncls0++;
361 AliDebug(2, Form("n clusters = %d", ncls0));
362 if (ncls0 >= fMinITSClusters){
363 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters) ;
364 icountRecoITSClusters++;
365 AliDebug(2,Form("pT = %f, dca = %f, cosThetaStar = %f, pTpi = %f, pTK = %f, d0pi = %f, d0K = %f, d0xd0 = %f, cosPointingAngle = %f", pt, dca, cosThetaStar,pTpi, pTK, d0pi*1E4, d0K*1E4, d0xd0*1E8, cosPointingAngle));
366
367 // PPR cuts
368 Double_t cuts[6];
369 if (pt <= 1){
370 cuts[0] = 400;
371 cuts[1] = 0.8;
372 cuts[2] = 0.5;
373 cuts[3] = 500;
374 cuts[4] = -20000;
375 cuts[5] = 0.5;
376 }
377 else if (pt > 1 && pt <= 2){
378 cuts[0] = 300;
379 cuts[1] = 0.8;
380 cuts[2] = 0.6;
381 cuts[3] = 500;
382 cuts[4] = -20000;
383 cuts[5] = 0.6;
384 }
385 else if (pt > 2 && pt <= 3){
386 cuts[0] = 200;
387 cuts[1] = 0.8;
388 cuts[2] = 0.7;
389 cuts[3] = 500;
390 cuts[4] = -20000;
391 cuts[5] = 0.8;
392 }
393 else if (pt > 3 && pt <= 5){
394 cuts[0] = 200;
395 cuts[1] = 0.8;
396 cuts[2] = 0.7;
397 cuts[3] = 500;
398 cuts[4] = -10000;
399 cuts[5] = 0.8;
400 }
401 else if (pt > 5){
402 cuts[0] = 200;
403 cuts[1] = 0.8;
404 cuts[2] = 0.7;
405 cuts[3] = 500;
406 cuts[4] = -5000;
407 cuts[5] = 0.8;
408 }
409 if (dca*1E4 < cuts[0]
410 && TMath::Abs(cosThetaStar) < cuts[1]
411 && pTpi > cuts[2]
412 && pTK > cuts[2]
413 && TMath::Abs(d0pi*1E4) < cuts[3]
414 && TMath::Abs(d0K*1E4) < cuts[3]
415 && d0xd0*1E8 < cuts[4]
416 && cosPointingAngle > cuts[5]
417 ){
418
419 AliDebug(2,"Particle passed PPR cuts");
420 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoPPR) ;
421 icountRecoPPR++;
422 }
423 else{
424 AliDebug(2,"Particle skipped due to PPR cuts");
425 if (dca*1E4 > cuts[0]){
426 AliDebug(2,"Problems with dca");
427 }
428 if ( TMath::Abs(cosThetaStar) > cuts[1]){
429 AliDebug(2,"Problems with cosThetaStar");
430 }
431 if (pTpi < cuts[2]){
432 AliDebug(2,"Problems with pTpi");
433 }
434 if (pTK < cuts[2]){
435 AliDebug(2,"Problems with pTK");
436 }
437 if (TMath::Abs(d0pi*1E4) > cuts[3]){
438 AliDebug(2,"Problems with d0pi");
439 }
440 if (TMath::Abs(d0K*1E4) > cuts[3]){
441 AliDebug(2,"Problems with d0K");
442 }
443 if (d0xd0*1E8 > cuts[4]){
444 AliDebug(2,"Problems with d0xd0");
445 }
446 if (cosPointingAngle < cuts[5]){
447 AliDebug(2,"Problems with cosPointingAngle");
448 }
449 }
450 }
451 }
452 }
453 }
454
455 AliDebug(2, Form("Found %i Reco particles that are D0!!",icountReco));
456
457 fCountReco+= icountReco;
458 fCountRecoAcc+= icountRecoAcc;
459 fCountRecoITSClusters+= icountRecoITSClusters;
460 fCountRecoPPR+= icountRecoPPR;
461
462 fHistEventsProcessed->Fill(0);
463 /* PostData(0) is taken care of by AliAnalysisTaskSE */
464 PostData(1,fHistEventsProcessed) ;
465 PostData(2,fCFManager->GetParticleContainer()) ;
466}
467
468
469//___________________________________________________________________________
470void AliCFHeavyFlavourTaskMultiVarMultiStep::Terminate(Option_t*)
471{
472 // The Terminate() function is the last function to be called during
473 // a query. It always runs on the client, it can be used to present
474 // the results graphically or save the results to file.
475
476 Info("Terminate","");
477 AliAnalysisTaskSE::Terminate();
478
479 AliInfo(Form("Found %i MC particles that are D0 in MC, in %d events",fCountMC,fEvents));
480 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, in %d events",fCountAcc,fEvents));
481 AliInfo(Form("Found %i reco D0 that are decaying in K+pi, in %d events",fCountReco,fEvents));
482 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and are in the requested acceptance, in %d events",fCountRecoAcc,fEvents));
483 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and have at least %d clusters in ITS, in %d events",fCountRecoITSClusters,fMinITSClusters,fEvents));
484 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and satisfy PPR cuts, in %d events",fCountRecoPPR,fEvents));
485
486 // draw some example plots....
487
488 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
489
490 // projecting the containers to obtain histograms
491 // first argument = variable, second argument = step
492
493 // MC-level
494 TH1D* h00 = cont->ShowProjection(0,0) ; // pt
495 TH1D* h10 = cont->ShowProjection(1,0) ; // rapidity
496 TH1D* h20 = cont->ShowProjection(2,0) ; // cosThetaStar
497 TH1D* h30 = cont->ShowProjection(3,0) ; // pTpi
498 TH1D* h40 = cont->ShowProjection(4,0) ; // pTK
499 TH1D* h50 = cont->ShowProjection(5,0) ; // cT
500 TH1D* h60 = cont->ShowProjection(6,0) ; // dca
501 TH1D* h70 = cont->ShowProjection(7,0) ; // d0pi
502 TH1D* h80 = cont->ShowProjection(8,0) ; // d0K
503 TH1D* h90 = cont->ShowProjection(9,0) ; // d0xd0
504 TH1D* h100 = cont->ShowProjection(10,0) ; // cosPointingAngle
505
506 // MC-Acceptance level
507 TH1D* h01 = cont->ShowProjection(0,1) ; // pt
508 TH1D* h11 = cont->ShowProjection(1,1) ; // rapidity
509 TH1D* h21 = cont->ShowProjection(2,1) ; // cosThetaStar
510 TH1D* h31 = cont->ShowProjection(3,1) ; // pTpi
511 TH1D* h41 = cont->ShowProjection(4,1) ; // pTK
512 TH1D* h51 = cont->ShowProjection(5,1) ; // cT
513 TH1D* h61 = cont->ShowProjection(6,1) ; // dca
514 TH1D* h71 = cont->ShowProjection(7,1) ; // d0pi
515 TH1D* h81 = cont->ShowProjection(8,1) ; // d0K
516 TH1D* h91 = cont->ShowProjection(9,1) ; // d0xd0
517 TH1D* h101 = cont->ShowProjection(10,1) ; // cosPointingAngle
518
519 // Reco-level
520 TH1D* h02 = cont->ShowProjection(0,2) ; // pt
521 TH1D* h12 = cont->ShowProjection(1,2) ; // rapidity
522 TH1D* h22 = cont->ShowProjection(2,2) ; // cosThetaStar
523 TH1D* h32 = cont->ShowProjection(3,2) ; // pTpi
524 TH1D* h42 = cont->ShowProjection(4,2) ; // pTK
525 TH1D* h52 = cont->ShowProjection(5,2) ; // cT
526 TH1D* h62 = cont->ShowProjection(6,2) ; // dca
527 TH1D* h72 = cont->ShowProjection(7,2) ; // d0pi
528 TH1D* h82 = cont->ShowProjection(8,2) ; // d0K
529 TH1D* h92 = cont->ShowProjection(9,2) ; // d0xd0
530 TH1D* h102 = cont->ShowProjection(10,2) ; // cosPointingAngle
531
532 h00->SetTitle("pT_D0 (GeV/c)");
533 h10->SetTitle("rapidity");
534 h20->SetTitle("cosThetaStar");
535 h30->SetTitle("pT_pi (GeV/c)");
536 h40->SetTitle("pT_K (Gev/c)");
537 h50->SetTitle("cT (#mum)");
538 h60->SetTitle("dca (#mum)");
539 h70->SetTitle("d0_pi (#mum)");
540 h80->SetTitle("d0_K (#mum)");
541 h90->SetTitle("d0xd0 (#mum^2)");
542 h100->SetTitle("cosPointingAngle");
543
544 h00->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
545 h10->GetXaxis()->SetTitle("rapidity");
546 h20->GetXaxis()->SetTitle("cosThetaStar");
547 h30->GetXaxis()->SetTitle("pT_pi (GeV/c)");
548 h40->GetXaxis()->SetTitle("pT_K (Gev/c)");
549 h50->GetXaxis()->SetTitle("cT (#mum)");
550 h60->GetXaxis()->SetTitle("dca (#mum)");
551 h70->GetXaxis()->SetTitle("d0_pi (#mum)");
552 h80->GetXaxis()->SetTitle("d0_K (#mum)");
553 h90->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
554 h100->GetXaxis()->SetTitle("cosPointingAngle");
555
556 h01->SetTitle("pT_D0 (GeV/c)");
557 h11->SetTitle("rapidity");
558 h21->SetTitle("cosThetaStar");
559 h31->SetTitle("pT_pi (GeV/c)");
560 h41->SetTitle("pT_K (Gev/c)");
561 h51->SetTitle("cT (#mum)");
562 h61->SetTitle("dca (#mum)");
563 h71->SetTitle("d0_pi (#mum)");
564 h81->SetTitle("d0_K (#mum)");
565 h91->SetTitle("d0xd0 (#mum^2)");
566 h101->SetTitle("cosPointingAngle");
567
568 h01->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
569 h11->GetXaxis()->SetTitle("rapidity");
570 h21->GetXaxis()->SetTitle("cosThetaStar");
571 h31->GetXaxis()->SetTitle("pT_pi (GeV/c)");
572 h41->GetXaxis()->SetTitle("pT_K (Gev/c)");
573 h51->GetXaxis()->SetTitle("cT (#mum)");
574 h61->GetXaxis()->SetTitle("dca (#mum)");
575 h71->GetXaxis()->SetTitle("d0_pi (#mum)");
576 h81->GetXaxis()->SetTitle("d0_K (#mum)");
577 h91->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
578 h101->GetXaxis()->SetTitle("cosPointingAngle");
579
580 h02->SetTitle("pT_D0 (GeV/c)");
581 h12->SetTitle("rapidity");
582 h22->SetTitle("cosThetaStar");
583 h32->SetTitle("pT_pi (GeV/c)");
584 h42->SetTitle("pT_K (Gev/c)");
585 h52->SetTitle("cT (#mum)");
586 h62->SetTitle("dca (#mum)");
587 h72->SetTitle("d0_pi (#mum)");
588 h82->SetTitle("d0_K (#mum)");
589 h92->SetTitle("d0xd0 (#mum^2)");
590 h102->SetTitle("cosPointingAngle");
591
592 h02->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
593 h12->GetXaxis()->SetTitle("rapidity");
594 h22->GetXaxis()->SetTitle("cosThetaStar");
595 h32->GetXaxis()->SetTitle("pT_pi (GeV/c)");
596 h42->GetXaxis()->SetTitle("pT_K (Gev/c)");
597 h52->GetXaxis()->SetTitle("cT (#mum)");
598 h62->GetXaxis()->SetTitle("dca (#mum)");
599 h72->GetXaxis()->SetTitle("d0_pi (#mum)");
600 h82->GetXaxis()->SetTitle("d0_K (#mum)");
601 h92->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
602 h102->GetXaxis()->SetTitle("cosPointingAngle");
603
604 Double_t max0 = h00->GetMaximum();
605 Double_t max1 = h10->GetMaximum();
606 Double_t max2 = h20->GetMaximum();
607 Double_t max3 = h30->GetMaximum();
608 Double_t max4 = h40->GetMaximum();
609 Double_t max5 = h50->GetMaximum();
610 Double_t max6 = h60->GetMaximum();
611 Double_t max7 = h70->GetMaximum();
612 Double_t max8 = h80->GetMaximum();
613 Double_t max9 = h90->GetMaximum();
614 Double_t max10 = h100->GetMaximum();
615
616 h00->GetYaxis()->SetRangeUser(0,max0*1.2);
617 h10->GetYaxis()->SetRangeUser(0,max1*1.2);
618 h20->GetYaxis()->SetRangeUser(0,max2*1.2);
619 h30->GetYaxis()->SetRangeUser(0,max3*1.2);
620 h40->GetYaxis()->SetRangeUser(0,max4*1.2);
621 h50->GetYaxis()->SetRangeUser(0,max5*1.2);
622 h60->GetYaxis()->SetRangeUser(0,max6*1.2);
623 h70->GetYaxis()->SetRangeUser(0,max7*1.2);
624 h80->GetYaxis()->SetRangeUser(0,max8*1.2);
625 h90->GetYaxis()->SetRangeUser(0,max9*1.2);
626 h100->GetYaxis()->SetRangeUser(0,max10*1.2);
627
628 h01->GetYaxis()->SetRangeUser(0,max0*1.2);
629 h11->GetYaxis()->SetRangeUser(0,max1*1.2);
630 h21->GetYaxis()->SetRangeUser(0,max2*1.2);
631 h31->GetYaxis()->SetRangeUser(0,max3*1.2);
632 h41->GetYaxis()->SetRangeUser(0,max4*1.2);
633 h51->GetYaxis()->SetRangeUser(0,max5*1.2);
634 h61->GetYaxis()->SetRangeUser(0,max6*1.2);
635 h71->GetYaxis()->SetRangeUser(0,max7*1.2);
636 h81->GetYaxis()->SetRangeUser(0,max8*1.2);
637 h91->GetYaxis()->SetRangeUser(0,max9*1.2);
638 h101->GetYaxis()->SetRangeUser(0,max10*1.2);
639
640 h00->SetMarkerStyle(20);
641 h10->SetMarkerStyle(24);
642 h20->SetMarkerStyle(21);
643 h30->SetMarkerStyle(25);
644 h40->SetMarkerStyle(27);
645 h50->SetMarkerStyle(28);
646 h60->SetMarkerStyle(20);
647 h70->SetMarkerStyle(24);
648 h80->SetMarkerStyle(21);
649 h90->SetMarkerStyle(25);
650 h100->SetMarkerStyle(27);
651
652 h00->SetMarkerColor(2);
653 h10->SetMarkerColor(2);
654 h20->SetMarkerColor(2);
655 h30->SetMarkerColor(2);
656 h40->SetMarkerColor(2);
657 h50->SetMarkerColor(2);
658 h60->SetMarkerColor(2);
659 h70->SetMarkerColor(2);
660 h80->SetMarkerColor(2);
661 h90->SetMarkerColor(2);
662 h100->SetMarkerColor(2);
663
664 h01->SetMarkerStyle(20) ;
665 h11->SetMarkerStyle(24) ;
666 h21->SetMarkerStyle(21) ;
667 h31->SetMarkerStyle(25) ;
668 h41->SetMarkerStyle(27) ;
669 h51->SetMarkerStyle(28) ;
670 h61->SetMarkerStyle(20);
671 h71->SetMarkerStyle(24);
672 h81->SetMarkerStyle(21);
673 h91->SetMarkerStyle(25);
674 h101->SetMarkerStyle(27);
675
676 h01->SetMarkerColor(8);
677 h11->SetMarkerColor(8);
678 h21->SetMarkerColor(8);
679 h31->SetMarkerColor(8);
680 h41->SetMarkerColor(8);
681 h51->SetMarkerColor(8);
682 h61->SetMarkerColor(8);
683 h71->SetMarkerColor(8);
684 h81->SetMarkerColor(8);
685 h91->SetMarkerColor(8);
686 h101->SetMarkerColor(8);
687
688 h02->SetMarkerStyle(20) ;
689 h12->SetMarkerStyle(24) ;
690 h22->SetMarkerStyle(21) ;
691 h32->SetMarkerStyle(25) ;
692 h42->SetMarkerStyle(27) ;
693 h52->SetMarkerStyle(28) ;
694 h62->SetMarkerStyle(20);
695 h72->SetMarkerStyle(24);
696 h82->SetMarkerStyle(21);
697 h92->SetMarkerStyle(25);
698 h102->SetMarkerStyle(27);
699
700 h02->SetMarkerColor(4);
701 h12->SetMarkerColor(4);
702 h22->SetMarkerColor(4);
703 h32->SetMarkerColor(4);
704 h42->SetMarkerColor(4);
705 h52->SetMarkerColor(4);
706 h62->SetMarkerColor(4);
707 h72->SetMarkerColor(4);
708 h82->SetMarkerColor(4);
709 h92->SetMarkerColor(4);
710 h102->SetMarkerColor(4);
711
712 gStyle->SetCanvasColor(0);
713 gStyle->SetFrameFillColor(0);
714 gStyle->SetTitleFillColor(0);
715 gStyle->SetStatColor(0);
716
717 // drawing in 2 separate canvas for a matter of clearity
718 TCanvas * c1 =new TCanvas("c1","pT, rapidiy, cosThetaStar",1100,1600);
719 c1->Divide(3,3);
720
721 c1->cd(1);
722 h00->Draw("p");
723 c1->cd(1);
724 c1->cd(2);
725 h01->Draw("p");
726 c1->cd(2);
727 c1->cd(3);
728 h02->Draw("p");
729 c1->cd(3);
730 c1->cd(4);
731 h10->Draw("p");
732 c1->cd(4);
733 c1->cd(5);
734 h11->Draw("p");
735 c1->cd(5);
736 c1->cd(6);
737 h12->Draw("p");
738 c1->cd(6);
739 c1->cd(7);
740 h20->Draw("p");
741 c1->cd(7);
742 c1->cd(8);
743 h21->Draw("p");
744 c1->cd(8);
745 c1->cd(9);
746 h22->Draw("p");
747 c1->cd(9);
748 c1->cd();
749
750 TCanvas * c2 =new TCanvas("c2","pTpi, pTK, cT",1100,1600);
751 c2->Divide(3,3);
752 c2->cd(1);
753 h30->Draw("p");
754 c2->cd(1);
755 c2->cd(2);
756 h31->Draw("p");
757 c2->cd(2);
758 c2->cd(3);
759 h32->Draw("p");
760 c2->cd(3);
761 c2->cd(4);
762 h40->Draw("p");
763 c2->cd(4);
764 c2->cd(5);
765 h41->Draw("p");
766 c2->cd(5);
767 c2->cd(6);
768 h42->Draw("p");
769 c2->cd(6);
770 c2->cd(7);
771 h50->Draw("p");
772 c2->cd(7);
773 c2->cd(8);
774 h51->Draw("p");
775 c2->cd(8);
776 c2->cd(9);
777 h52->Draw("p");
778 c2->cd(9);
779 c2->cd();
780
781 TCanvas * c3 =new TCanvas("c3","dca, d0pi, d0K",1100,1600);
782 c3->Divide(3,3);
783 c3->cd(1);
784 h60->Draw("p");
785 c3->cd(1);
786 c3->cd(2);
787 h61->Draw("p");
788 c3->cd(2);
789 c3->cd(3);
790 h62->Draw("p");
791 c3->cd(3);
792 c3->cd(4);
793 h70->Draw("p");
794 c3->cd(4);
795 c3->cd(5);
796 h71->Draw("p");
797 c3->cd(5);
798 c3->cd(6);
799 h72->Draw("p");
800 c3->cd(6);
801 c3->cd(7);
802 h80->Draw("p");
803 c3->cd(7);
804 c3->cd(8);
805 h81->Draw("p");
806 c3->cd(8);
807 c3->cd(9);
808 h82->Draw("p");
809 c3->cd(9);
810 c3->cd();
811
812 TCanvas * c4 =new TCanvas("c4","d0xd0, cosPointingAngle",1100,770);
813 c4->Divide(3,2);
814 c4->cd(1);
815 h90->Draw("p");
816 c4->cd(1);
817 c4->cd(2);
818 h91->Draw("p");
819 c4->cd(2);
820 c4->cd(3);
821 h92->Draw("p");
822 c4->cd(3);
823 c4->cd(4);
824 h100->Draw("p");
825 c4->cd(4);
826 c4->cd(5);
827 h101->Draw("p");
828 c4->cd(5);
829 c4->cd(6);
830 h102->Draw("p");
831 c4->cd(6);
832 c4->cd();
833
834 /*
835 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.eps");
836 c2->SaveAs("Plots/pTpi_pTK_cT.eps");
837 c3->SaveAs("Plots/dca_d0pi_d0TK.eps");
838 c4->SaveAs("Plots/d0xd0_cosPointingAngle.eps");
839
840 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.gif");
841 c2->SaveAs("Plots/pTpi_pTK_cT.gif");
842 c3->SaveAs("Plots/dca_d0pi_d0TK.gif");
843 c4->SaveAs("Plots/d0xd0_cosPointingAngle.gif");
844 */
845}
846
847//___________________________________________________________________________
848void AliCFHeavyFlavourTaskMultiVarMultiStep::UserCreateOutputObjects() {
849 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
850 //TO BE SET BEFORE THE EXECUTION OF THE TASK
851 //
852 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
853
854 //slot #1
855 OpenFile(1);
856 fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
857}
858
859//___________________________________________________________________________
860Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CosThetaStar(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
861
862 //
863 // to calculate cos(ThetaStar) for generated particle
864 // using the K, since mcPartDaughter0 and mcPartDaughter1 always correspond to K and pi respectively
865 // (see where the function is called)
866 //
867
868 Int_t pdgvtx = mcPart->GetPdgCode();
869 /* if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
870 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
871 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
872 AliInfo(Form("D0, with pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
873 AliDebug(2,"This is a D0");
874 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
875 mcPartDaughter0 = mcPartDaughter1;
876 mcPartDaughter1 = mcPartdummy;
877 }
878 else{
879 AliInfo("D0bar");
880 }
881 */
882 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
883 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
884 if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
885 AliDebug(2,"D0");
886 }
887 else{
888 AliDebug(2,"D0bar");
889 }
890 if (pdgprong0 == 211){
891 AliDebug(2,Form("pdgprong0 = %d, pdgprong1 = %d, switching...",pdgprong0,pdgprong1));
892 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
893 mcPartDaughter0 = mcPartDaughter1;
894 mcPartDaughter1 = mcPartdummy;
895 pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
896 pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
897 }
898
899 AliDebug(2,Form("After checking, pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
900 Double_t massvtx = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdgvtx))->Mass();
901 Double_t massp[2];
902 massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
903 massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
904
905 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);
906
907 Double_t px = mcPartDaughter0->Px()+mcPartDaughter1->Px();
908 Double_t py = mcPartDaughter0->Py()+mcPartDaughter1->Py();
909 Double_t pz = mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
910 Double_t p = TMath::Sqrt(px*px+py*py+pz*pz);
911 Double_t e = TMath::Sqrt(massvtx*massvtx+p*p);
912 Double_t beta = p/e;
913 Double_t gamma = e/massvtx;
914 TVector3 mom(mcPartDaughter0->Px(),mcPartDaughter0->Py(),mcPartDaughter0->Pz());
915 TVector3 momTot(mcPartDaughter0->Px()+mcPartDaughter1->Px(),mcPartDaughter0->Py()+mcPartDaughter1->Py(),mcPartDaughter0->Pz()+mcPartDaughter1->Pz());
916
917 Double_t qlprong = mom.Dot(momTot)/momTot.Mag(); // analog to AliAODRecoDecay::QlProng(0)
918
919 AliDebug(2,Form("pStar = %f, beta = %f, gamma = %f, qlprong = %f, massp[0] = %f", pStar, beta, gamma, qlprong, massp[0]));
920 Double_t cts = (qlprong/gamma-beta*TMath::Sqrt(pStar*pStar+massp[0]*massp[0]))/pStar;
921 AliDebug(2,Form("cts = %f", cts));
922 return cts;
923}
924//___________________________________________________________________________
925Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CT(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
926
927 //
928 // to calculate cT for generated particle
929 //
930
931 Double_t xmcPart[3] = {0,0,0};
932 Double_t xdaughter0[3] = {0,0,0};
933 Double_t xdaughter1[3] = {0,0,0};
934 mcPart->XvYvZv(xmcPart); // cm
935 mcPartDaughter0->XvYvZv(xdaughter0); // cm
936 mcPartDaughter1->XvYvZv(xdaughter1); //cm
937 Double_t prodVtxD0 = TMath::Sqrt(xmcPart[0]*xmcPart[0]+
938 xmcPart[1]*xmcPart[1]+
939 xmcPart[2]*xmcPart[2]);
940 Double_t prodVtxDaughter0 = TMath::Sqrt(xdaughter0[0]*xdaughter0[0]+
941 xdaughter0[1]*xdaughter0[1]+
942 xdaughter0[2]*xdaughter0[2]);
943 Double_t prodVtxDaughter1 = TMath::Sqrt(xdaughter1[0]*xdaughter1[0]+
944 xdaughter1[1]*xdaughter1[1]+
945 xdaughter1[2]*xdaughter1[2]);
946
947 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));
948 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));
949 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));
950
951 Double_t cT0 = TMath::Sqrt((xdaughter0[0]-xmcPart[0])*(xdaughter0[0]-xmcPart[0])+
952 (xdaughter0[1]-xmcPart[1])*(xdaughter0[1]-xmcPart[1])+
953 (xdaughter0[2]-xmcPart[2])*(xdaughter0[2]-xmcPart[2]));
954
955 Double_t cT1 = TMath::Sqrt((xdaughter1[0]-xmcPart[0])*(xdaughter1[0]-xmcPart[0])+
956 (xdaughter1[1]-xmcPart[1])*(xdaughter1[1]-xmcPart[1])+
957 (xdaughter1[2]-xmcPart[2])*(xdaughter1[2]-xmcPart[2]));
958
959 if (cT0 != cT1) {
960 AliWarning("cT from daughter 0 different from cT from daughter 1! Using cT from daughter 0, but PLEASE, CHECK....");
961 }
962
963 // calculating cT from cT0
964
965 Double_t mass = TDatabasePDG::Instance()->GetParticle(mcPart->GetPdgCode())->Mass(); // mcPart->GetPdgCode() should return +/- 421 for the D0/D0bar
966 Double_t p = mcPart-> P();
967 Double_t cT = cT0*mass/p;
968 AliDebug(2, Form("cT from daughter 0 = %f (micron)", cT0*1E4));
969 AliDebug(2, Form("cT from daughter 1 = %f (micron)", cT1*1E4));
970 AliDebug(2, Form("cT (with mass = %f and p = %f) = %f (micron)", mass, p, cT*1E4));
971 return cT;
972}
973//________________________________________________________________________________
974Bool_t AliCFHeavyFlavourTaskMultiVarMultiStep::GetGeneratedValuesFromMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC) const {
975
976 //
977 // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
978 //
979
980 Bool_t isOk = kFALSE;
981
982 // check whether the D0 decays in pi+K
983 // to be added!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
984 // could use a cut in the CF, but not clear how to define a TDecayChannel
985 // implemented for the time being as a cut on the number of daughters - see later when
986 // getting the daughters
987
988 // getting the daughters
989 Int_t daughter0 = mcPart->GetDaughter(0);
990 Int_t daughter1 = mcPart->GetDaughter(1);
991 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
992 if (daughter0 == 0 || daughter1 == 0) {
993 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
994 return isOk;
995 }
996 if (TMath::Abs(daughter1 - daughter0 != 1)) {
997 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
998 return isOk;
999 }
1000 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
1001 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
1002 if (!mcPartDaughter0 || !mcPartDaughter1) {
1003 AliWarning("At least one Daughter Particle not found in tree, skipping");
1004 return isOk;
1005 }
1006
1007 Double_t vtx1[3] = {0,0,0}; // primary vertex
1008 Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0
1009 Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1
1010 mcPart->XvYvZv(vtx1); // cm
1011 // getting vertex from daughters
1012 mcPartDaughter0->XvYvZv(vtx2daughter0); // cm
1013 mcPartDaughter1->XvYvZv(vtx2daughter1); //cm
1014 if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
1015 AliError("Daughters have different secondary vertex, skipping the track");
1016 return isOk;
1017 }
1018 Int_t nprongs = 2;
1019 Short_t charge = 0;
1020 // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
1021 AliAODMCParticle* positiveDaugh = mcPartDaughter0;
1022 AliAODMCParticle* negativeDaugh = mcPartDaughter1;
1023 if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){
1024 // inverting in case the positive daughter is the second one
1025 positiveDaugh = mcPartDaughter1;
1026 negativeDaugh = mcPartDaughter0;
1027 }
1028 // getting the momentum from the daughters
1029 Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};
1030 Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};
1031 Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
1032
1033 Double_t d0[2] = {0.,0.};
1034
1035 AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
1036
1037 Double_t cosThetaStar = 0.;
1038 Double_t cosThetaStarD0 = 0.;
1039 Double_t cosThetaStarD0bar = 0.;
1040 cosThetaStarD0 = decay->CosThetaStar(1,421,211,321);
1041 cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211);
1042 if (mcPart->GetPdgCode() == 421){ // D0
1043 AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1044 cosThetaStar = cosThetaStarD0;
1045 }
1046 else if (mcPart->GetPdgCode() == -421){ // D0bar{
1047 AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1048 cosThetaStar = cosThetaStarD0bar;
1049 }
1050 else{
1051 AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check...");
1052 return vectorMC;
1053 }
1054 if (cosThetaStar < -1 || cosThetaStar > 1) {
1055 AliWarning("Invalid value for cosine Theta star, returning");
1056 return isOk;
1057 }
1058
1059 // calculate cos(Theta*) according to the method implemented herein
1060
1061 Double_t cts = 9999.;
1062 cts = CosThetaStar(mcPart, mcPartDaughter0, mcPartDaughter1);
1063 if (cts < -1 || cts > 1) {
1064 AliWarning("Invalid value for cosine Theta star from AliCFHeavyFlavourTaskMultiVarMultiStep method");
1065 }
1066 if (TMath::Abs(cts - cosThetaStar)>0.001) {
1067 AliError(Form("cosThetaStar automatically calculated different from that manually calculated!!! cosThetaStar = %f, cosThetaStar = %f", cosThetaStar,cts));
1068 }
1069
1070 Double_t cT = decay->Ct(421);
1071
1072 // calculate cT from the method implemented herein
1073 Double_t cT1 = 0.;
1074 cT1 = CT(mcPart, mcPartDaughter0, mcPartDaughter1);
1075
1076 if (TMath::Abs(cT1 - cT)>0.001) {
1077 AliError(Form("cT automatically calculated different from that manually calculated!!! cT = %f, cT1 = %f",cT,cT1));
1078 }
1079
1080 // get the pT of the daughters
1081
1082 Double_t pTpi = 0.;
1083 Double_t pTK = 0.;
1084
1085 if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
1086 pTpi = mcPartDaughter0->Pt();
1087 pTK = mcPartDaughter1->Pt();
1088 }
1089 else {
1090 pTpi = mcPartDaughter1->Pt();
1091 pTK = mcPartDaughter0->Pt();
1092 }
1093
1094 vectorMC[0] = mcPart->Pt();
1095 vectorMC[1] = mcPart->Y() ;
1096 vectorMC[2] = cosThetaStar ;
1097 vectorMC[3] = pTpi ;
1098 vectorMC[4] = pTK ;
1099 vectorMC[5] = cT*1.E4 ; // in micron
1100 isOk = kTRUE;
1101 return isOk;
1102}
1103