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