]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGHF/vertexingHF/AliCFHeavyFlavourTaskMultiVarMultiStep.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Transition PWG3 --> PWGHF
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliCFHeavyFlavourTaskMultiVarMultiStep.cxx
... / ...
CommitLineData
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/* $Id$ */
17
18//-----------------------------------------------------------------------
19// Class for HF corrections as a function of many variables
20// 6 Steps introduced: MC, MC Acc, Reco, Reco Acc, Reco Acc + ITS Cl,
21// Reco Acc + ITS Cl + PPR cuts
22// 13 variables used: pt, y, cosThetaStar, ptPi, ptK, ct,
23// dca, d0Pi, d0K, d0Pixd0K, cosPointingAngle, phi, z
24//
25//-----------------------------------------------------------------------
26// Author : C. Zampolli, CERN
27//-----------------------------------------------------------------------
28//-----------------------------------------------------------------------
29// Base class for HF Unfolding (pt and eta)
30// correlation matrix filled at Acceptance and PPR level
31// Author: A.Grelli , Utrecht - agrelli@uu.nl
32//-----------------------------------------------------------------------
33#include <TCanvas.h>
34#include <TParticle.h>
35#include <TDatabasePDG.h>
36#include <TH1I.h>
37#include <TStyle.h>
38#include <TFile.h>
39
40#include "AliCFHeavyFlavourTaskMultiVarMultiStep.h"
41#include "AliStack.h"
42#include "AliMCEvent.h"
43#include "AliCFManager.h"
44#include "AliCFContainer.h"
45#include "AliLog.h"
46#include "AliAnalysisManager.h"
47#include "AliAODHandler.h"
48#include "AliAODEvent.h"
49#include "AliAODRecoDecay.h"
50#include "AliAODRecoDecayHF.h"
51#include "AliAODRecoDecayHF2Prong.h"
52#include "AliAODMCParticle.h"
53#include "AliAODMCHeader.h"
54#include "AliESDtrack.h"
55#include "AliRDHFCutsD0toKpi.h"
56#include "TChain.h"
57#include "THnSparse.h"
58#include "TH2D.h"
59#include "AliAnalysisDataSlot.h"
60#include "AliAnalysisDataContainer.h"
61
62//__________________________________________________________________________
63AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep() :
64 AliAnalysisTaskSE(),
65 fPDG(0),
66 fCFManager(0x0),
67 fHistEventsProcessed(0x0),
68 fCorrelation(0x0),
69 fCountMC(0),
70 fCountAcc(0),
71 fCountVertex(0),
72 fCountRefit(0),
73 fCountReco(0),
74 fCountRecoAcc(0),
75 fCountRecoITSClusters(0),
76 fCountRecoPPR(0),
77 fCountRecoPID(0),
78 fEvents(0),
79 fFillFromGenerated(kFALSE),
80 fMinITSClusters(5),
81 fAcceptanceUnf(kTRUE),
82 fKeepD0fromB(kFALSE),
83 fKeepD0fromBOnly(kFALSE),
84 fCuts(0),
85 fUseWeight(kFALSE),
86 fWeight(1.),
87 fSign(2)
88{
89 //
90 //Default ctor
91 //
92}
93//___________________________________________________________________________
94AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep(const Char_t* name, AliRDHFCutsD0toKpi* cuts) :
95 AliAnalysisTaskSE(name),
96 fPDG(0),
97 fCFManager(0x0),
98 fHistEventsProcessed(0x0),
99 fCorrelation(0x0),
100 fCountMC(0),
101 fCountAcc(0),
102 fCountVertex(0),
103 fCountRefit(0),
104 fCountReco(0),
105 fCountRecoAcc(0),
106 fCountRecoITSClusters(0),
107 fCountRecoPPR(0),
108 fCountRecoPID(0),
109 fEvents(0),
110 fFillFromGenerated(kFALSE),
111 fMinITSClusters(5),
112 fAcceptanceUnf(kTRUE),
113 fKeepD0fromB(kFALSE),
114 fKeepD0fromBOnly(kFALSE),
115 fCuts(cuts),
116 fUseWeight(kFALSE),
117 fWeight(1.),
118 fSign(2)
119{
120 //
121 // Constructor. Initialization of Inputs and Outputs
122 //
123 Info("AliCFHeavyFlavourTaskMultiVarMultiStep","Calling Constructor");
124 /*
125 DefineInput(0) and DefineOutput(0)
126 are taken care of by AliAnalysisTaskSE constructor
127 */
128 DefineOutput(1,TH1I::Class());
129 DefineOutput(2,AliCFContainer::Class());
130 DefineOutput(3,THnSparseD::Class());
131 DefineOutput(4,AliRDHFCutsD0toKpi::Class());
132
133 fCuts->PrintAll();
134}
135
136//___________________________________________________________________________
137AliCFHeavyFlavourTaskMultiVarMultiStep& AliCFHeavyFlavourTaskMultiVarMultiStep::operator=(const AliCFHeavyFlavourTaskMultiVarMultiStep& c)
138{
139 //
140 // Assignment operator
141 //
142 if (this!=&c) {
143 AliAnalysisTaskSE::operator=(c) ;
144 fPDG = c.fPDG;
145 fCFManager = c.fCFManager;
146 fHistEventsProcessed = c.fHistEventsProcessed;
147 fCuts = c.fCuts;
148 }
149 return *this;
150}
151
152//___________________________________________________________________________
153AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep(const AliCFHeavyFlavourTaskMultiVarMultiStep& c) :
154 AliAnalysisTaskSE(c),
155 fPDG(c.fPDG),
156 fCFManager(c.fCFManager),
157 fHistEventsProcessed(c.fHistEventsProcessed),
158 fCorrelation(c.fCorrelation),
159 fCountMC(c.fCountMC),
160 fCountAcc(c.fCountAcc),
161 fCountVertex(c.fCountVertex),
162 fCountRefit(c.fCountRefit),
163 fCountReco(c.fCountReco),
164 fCountRecoAcc(c.fCountRecoAcc),
165 fCountRecoITSClusters(c.fCountRecoITSClusters),
166 fCountRecoPPR(c.fCountRecoPPR),
167 fCountRecoPID(c.fCountRecoPID),
168 fEvents(c.fEvents),
169 fFillFromGenerated(c.fFillFromGenerated),
170 fMinITSClusters(c.fMinITSClusters),
171 fAcceptanceUnf(c.fAcceptanceUnf),
172 fKeepD0fromB(c.fKeepD0fromB),
173 fKeepD0fromBOnly(c.fKeepD0fromBOnly),
174 fCuts(c.fCuts),
175 fUseWeight(c.fUseWeight),
176 fWeight(c.fWeight),
177 fSign(c.fSign)
178{
179 //
180 // Copy Constructor
181 //
182}
183
184//___________________________________________________________________________
185AliCFHeavyFlavourTaskMultiVarMultiStep::~AliCFHeavyFlavourTaskMultiVarMultiStep() {
186 //
187 //destructor
188 //
189 if (fCFManager) delete fCFManager ;
190 if (fHistEventsProcessed) delete fHistEventsProcessed ;
191 if (fCorrelation) delete fCorrelation ;
192 // if (fCuts) delete fCuts ;
193}
194
195//________________________________________________________________________
196void AliCFHeavyFlavourTaskMultiVarMultiStep::Init(){
197
198 //
199 // Initialization
200 //
201
202 if(fDebug > 1) printf("AliCFHeavyFlavourTaskMultiVarMultiStep::Init() \n");
203
204 fMinITSClusters = fCuts->GetTrackCuts()->GetMinNClustersITS();
205 AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*fCuts);
206 const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
207 copyfCuts->SetName(nameoutput);
208 // Post the data
209 PostData(4,copyfCuts);
210
211 return;
212}
213
214//_________________________________________________
215void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
216{
217 //
218 // Main loop function
219 //
220
221 PostData(1,fHistEventsProcessed) ;
222 PostData(2,fCFManager->GetParticleContainer()) ;
223 PostData(3,fCorrelation) ;
224
225 AliESDtrackCuts* trackCuts = fCuts->GetTrackCuts(); // track cuts
226
227 if (fFillFromGenerated){
228 AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
229 }
230
231 if (!fInputEvent) {
232 Error("UserExec","NO EVENT FOUND!");
233 return;
234 }
235
236 // check that the fKeepD0fromB flag is set to true when the fKeepD0fromBOnly flag is true
237 if(fKeepD0fromBOnly) {
238 fKeepD0fromB=true;
239 if(fEvents<2) AliInfo(Form("Both fKeepD0fromB and fKeepD0fromBOnly flags are true, looking _ONLY_ at D0 FROM B"));
240 }
241
242 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
243
244 TClonesArray *arrayD0toKpi=0;
245
246 if(!aodEvent && AODEvent() && IsStandardAOD()) {
247 // In case there is an AOD handler writing a standard AOD, use the AOD
248 // event in memory rather than the input (ESD) event.
249 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
250 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
251 // have to taken from the AOD event hold by the AliAODExtension
252 AliAODHandler* aodHandler = (AliAODHandler*)
253 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
254 if(aodHandler->GetExtensions()) {
255 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
256 AliAODEvent *aodFromExt = ext->GetAOD();
257 arrayD0toKpi=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
258 }
259 } else {
260 arrayD0toKpi=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
261 }
262
263
264 if (!arrayD0toKpi) {
265 AliError("Could not find array of HF vertices");
266 return;
267 }
268
269 // fix for temporary bug in ESDfilter
270 // the AODs with null vertex pointer didn't pass the PhysSel
271 if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
272
273 fEvents++;
274
275 fCFManager->SetRecEventInfo(aodEvent);
276 fCFManager->SetMCEventInfo(aodEvent);
277
278 // MC-event selection
279 Double_t containerInput[13] ;
280 Double_t containerInputMC[13] ;
281
282 //loop on the MC event
283
284 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
285 if (!mcArray) {
286 AliError("Could not find Monte-Carlo in AOD");
287 return;
288 }
289 Int_t icountMC = 0;
290 Int_t icountAcc = 0;
291 Int_t icountReco = 0;
292 Int_t icountVertex = 0;
293 Int_t icountRefit = 0;
294 Int_t icountRecoAcc = 0;
295 Int_t icountRecoITSClusters = 0;
296 Int_t icountRecoPPR = 0;
297 Int_t icountRecoPID = 0;
298
299 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
300 if (!mcHeader) {
301 AliError("Could not find MC Header in AOD");
302 return;
303 }
304
305 Int_t cquarks = 0;
306
307 // AOD primary vertex
308 AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
309 if(!vtx1) {
310 AliError("There is no primary vertex !");
311 return;
312 }
313 Double_t zPrimVertex = vtx1->GetZ();
314 Double_t zMCVertex = mcHeader->GetVtxZ();
315 Bool_t vtxFlag = kTRUE;
316 TString title=vtx1->GetTitle();
317 if(!title.Contains("VertexerTracks")) vtxFlag=kFALSE;
318
319 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
320 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
321 if (!mcPart) {
322 AliWarning("Particle not found in tree, skipping");
323 continue;
324 }
325 if (mcPart->GetPdgCode() == 4) cquarks++;
326 if (mcPart->GetPdgCode() == -4) cquarks++;
327
328 // check the MC-level cuts
329
330 if (!fCFManager->CheckParticleCuts(0, mcPart)) continue; // 0 stands for MC level
331 Int_t pdgGranma = CheckOrigin(mcPart, mcArray);
332 Int_t abspdgGranma = TMath::Abs(pdgGranma);
333 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
334 AliDebug(2,Form("Particle has a b-meson, or b-baryon mother (pdg code mother = %d )--> not coming from a c-quark, skipping...", pdgGranma));
335 if (!fKeepD0fromB) continue; // skipping particles that don't come from c quark
336 }
337 else { if(fKeepD0fromBOnly) continue; } // skipping particles that don't come from b quark
338
339 // if (TMath::Abs(pdgGranma)!=4) {
340
341 // fill the container for Gen-level selection
342 Double_t vectorMC[7] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.};
343
344 if (GetGeneratedValuesFromMCParticle(mcPart, mcArray, vectorMC)){
345 containerInputMC[0] = vectorMC[0];
346 containerInputMC[1] = vectorMC[1] ;
347 containerInputMC[2] = vectorMC[2] ;
348 containerInputMC[3] = vectorMC[3] ;
349 containerInputMC[4] = vectorMC[4] ;
350 containerInputMC[5] = vectorMC[5] ; // in micron
351 containerInputMC[6] = 0.; // dummy value, meaningless in MC, in micron
352 containerInputMC[7] = 0.; // dummy value, meaningless in MC, in micron
353 containerInputMC[8] = 0.; // dummy value, meaningless in MC, in micron
354 containerInputMC[9] = -100000.; // dummy value, meaningless in MC, in micron^2
355 containerInputMC[10] = 1.01; // dummy value, meaningless in MC
356 containerInputMC[11] = vectorMC[6]; // dummy value, meaningless in MC
357 containerInputMC[12] = zMCVertex; // z of reconstructed of primary vertex
358 if (fUseWeight) fWeight = GetWeight(vectorMC[0]); // setting the weight according to the function defined in AliCFHeavyFlavourTaskMultiVarMultiStep::GetWeight(Float_t pt)
359 AliDebug(3,Form("weight = %f",fWeight));
360 if (!fCuts->IsInFiducialAcceptance(vectorMC[0],vectorMC[1])) continue;
361 if (TMath::Abs(vectorMC[1]) < 0.5) {
362 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc,fWeight);
363 }
364 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGenerated,fWeight);
365 icountMC++;
366
367 // check the MC-Acceptance level cuts
368 // since standard CF functions are not applicable, using Kine Cuts on daughters
369
370 Int_t daughter0 = mcPart->GetDaughter(0);
371 Int_t daughter1 = mcPart->GetDaughter(1);
372 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
373 if (daughter0 == 0 || daughter1 == 0) {
374 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!! But it should have, this check was already done...");
375 }
376 if (TMath::Abs(daughter1 - daughter0) != 1) {
377 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, but it should be, this check was already done...");
378 }
379 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
380 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
381 if (!mcPartDaughter0 || !mcPartDaughter1) {
382 AliWarning("At least one Daughter Particle not found in tree, but it should be, this check was already done...");
383 continue;
384 }
385 Double_t eta0 = mcPartDaughter0->Eta();
386 Double_t eta1 = mcPartDaughter1->Eta();
387 Double_t y0 = mcPartDaughter0->Y();
388 Double_t y1 = mcPartDaughter1->Y();
389 Double_t pt0 = mcPartDaughter0->Pt();
390 Double_t pt1 = mcPartDaughter1->Pt();
391 AliDebug(2, Form("Daughter 0: eta = %f, y = %f, pt = %f", eta0, y0, pt0));
392 AliDebug(2, Form("Daughter 1: eta = %f, y = %f, pt = %f", eta1, y1, pt1));
393 Bool_t daught0inAcceptance = (TMath::Abs(eta0) < 0.9 && pt0 > 0.1);
394 Bool_t daught1inAcceptance = (TMath::Abs(eta1) < 0.9 && pt1 > 0.1);
395 if (daught0inAcceptance && daught1inAcceptance) {
396 // checking whether the cuts implemented in the CF are equivalent - simply a cross-check
397 AliDebug(2, "Daughter particles in acceptance");
398 if (!fCFManager->CheckParticleCuts(1, mcPartDaughter0)) {
399 AliDebug(2,"Inconsistency with CF cut for daughter 0!");
400 }
401 if (!fCFManager->CheckParticleCuts(1, mcPartDaughter1)) {
402 AliDebug(2,"Inconsistency with CF cut for daughter 1!");
403 }
404 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance,fWeight);
405 icountAcc++;
406
407 // check on the vertex
408 if (fCuts->IsEventSelected(aodEvent)){
409 AliDebug(2,"Vertex cut passed\n");
410 // filling the container if the vertex is ok
411 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex,fWeight) ;
412 icountVertex++;
413 // check on the kTPCrefit and kITSrefit conditions of the daughters
414 Bool_t refitFlag = kTRUE;
415 if (trackCuts->GetRequireTPCRefit() || trackCuts->GetRequireITSRefit()){
416 Int_t foundDaughters = 0;
417 for (Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
418 AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
419 if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
420 if ((track->GetLabel() == daughter0) || (track->GetLabel() == daughter1)) {
421 foundDaughters++;
422 if (trackCuts->GetRequireTPCRefit()) {
423 if(!(track->GetStatus()&AliESDtrack::kTPCrefit)){
424 refitFlag = kFALSE;
425 break;
426 }
427 }
428 if (trackCuts->GetRequireITSRefit()) {
429 if(!(track->GetStatus()&AliESDtrack::kITSrefit)){
430 refitFlag = kFALSE;
431 break;
432 }
433 }
434 }
435 if (foundDaughters == 2){ // both daughters have been checked
436 break;
437 }
438 }
439 if (foundDaughters != 2) refitFlag = kFALSE;
440 }
441 if (refitFlag){
442 AliDebug(3,"Refit cut passed\n");
443 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit,fWeight);
444 icountRefit++;
445 }
446 else{
447 AliDebug(3,"Refit cut not passed\n");
448 }
449 }
450 else{
451 AliDebug(3,"Vertex cut not passed\n");
452 }
453 }
454 else{
455 AliDebug(3,"Acceptance cut not passed\n");
456 }
457 }
458 else {
459 AliDebug(3,"Problems in filling the container");
460 continue;
461 }
462 }
463
464 if (cquarks<2) AliDebug(2,Form("Event found with %d c-quarks", cquarks));
465
466 AliDebug(2,Form("Found %i MC particles that are D0!!",icountMC));
467 AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Acc cuts!!",icountAcc));
468 AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Vertex cuts!!",icountVertex));
469 AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Refit cuts!!",icountRefit));
470
471 // Now go to rec level
472 fCountMC += icountMC;
473 fCountAcc += icountAcc;
474
475 AliDebug(2, Form("Found %d vertices",arrayD0toKpi->GetEntriesFast()));
476
477 Int_t pdgDgD0toKpi[2]={321,211};
478 Int_t isD0D0bar=1;// 1 for D0, 2 for D0bar, to be used for the PPR and PID selection steps
479
480 for (Int_t iD0toKpi = 0; iD0toKpi<arrayD0toKpi->GetEntriesFast(); iD0toKpi++) {
481
482 AliAODRecoDecayHF2Prong* d0tokpi = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->At(iD0toKpi);
483 if(!d0tokpi) continue;
484 Bool_t unsetvtx=kFALSE;
485 if(!d0tokpi->GetOwnPrimaryVtx()) {
486 d0tokpi->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
487 unsetvtx=kTRUE;
488 }
489
490 // find associated MC particle
491 Int_t mcLabel = d0tokpi->MatchToMC(421,mcArray,2,pdgDgD0toKpi) ;
492 if (mcLabel == -1)
493 {
494 AliDebug(2,"No MC particle found");
495 continue;
496 }
497 else {
498 AliAODMCParticle* mcVtxHF = (AliAODMCParticle*)mcArray->At(mcLabel);
499 if (!mcVtxHF) {
500 AliWarning("Could not find associated MC in AOD MC tree");
501 continue;
502 }
503
504 if (mcVtxHF->GetPdgCode() == 421){ // particle is D0
505 if (fSign == 1){ // I ask for D0bar only
506 AliDebug(2,"particle is D0, I ask for D0bar only");
507 continue;
508 }
509 else{
510 isD0D0bar=1;
511 }
512 }
513 else if (mcVtxHF->GetPdgCode()== -421){ // particle is D0bar
514 if (fSign == 0){ // I ask for D0 only
515 AliDebug(2,"particle is D0bar, I ask for D0 only");
516 continue;
517 }
518 else{
519 isD0D0bar=2;
520 }
521 }
522 else continue;
523
524 // check whether the daughters have kTPCrefit and kITSrefit set
525 AliAODTrack *track0 = (AliAODTrack*)d0tokpi->GetDaughter(0);
526 AliAODTrack *track1 = (AliAODTrack*)d0tokpi->GetDaughter(1);
527 if( !track0 || !track1 ) {
528 AliWarning("Could not find associated MC daughter tracks in AOD MC tree");
529 continue;
530 }
531 if ((trackCuts->GetRequireTPCRefit() && (!(track0->GetStatus()&AliESDtrack::kTPCrefit) || !(track1->GetStatus()&AliESDtrack::kTPCrefit))) ||
532 (trackCuts->GetRequireITSRefit() && (!(track0->GetStatus()&AliESDtrack::kITSrefit) || !(track1->GetStatus()&AliESDtrack::kITSrefit)))){
533 // skipping if at least one daughter does not have kTPCrefit or kITSrefit, if they were required
534 continue;
535 }
536
537 // check on the vertex -- could be moved outside the loop on the reconstructed D0...
538 if(!fCuts->IsEventSelected(aodEvent)) {
539 // skipping cause vertex is not ok
540 continue;
541 }
542
543 const Double_t d0tokpiCuts[9] = {0.3,999999.,1.1,0.,0.,999999.,999999.,999999.,0.};
544 Int_t okD0, okD0bar;
545 if (!(d0tokpi->SelectD0(&d0tokpiCuts[0],okD0,okD0bar))){
546 // skipping candidate
547 continue;
548 }
549
550 // check if associated MC v0 passes the cuts
551 if (!fCFManager->CheckParticleCuts(0 ,mcVtxHF)) { // 0 stands for MC
552 AliDebug(2, "Skipping the particles due to cuts");
553 continue;
554 }
555 Int_t pdgGranma = CheckOrigin(mcVtxHF, mcArray);
556 Int_t abspdgGranma = TMath::Abs(pdgGranma);
557 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
558 AliDebug(2,Form("At Reco level, from MC info: Particle has a b-meson, or b-baryon mother (pdg code mother = %d )--> not coming from a c-quark, skipping...", pdgGranma));
559 if (!fKeepD0fromB) continue; // skipping particles that don't come from c quark
560 }
561 else { if(fKeepD0fromBOnly) continue; } // skipping particles that don't come from b quark
562
563 // fill the container...
564 Double_t pt = d0tokpi->Pt();
565 Double_t rapidity = d0tokpi->YD0();
566 Double_t invMass=0.;
567 Double_t cosThetaStar = 9999.;
568 Double_t pTpi = 0.;
569 Double_t pTK = 0.;
570 Double_t dca = d0tokpi->GetDCA();
571 Double_t d0pi = 0.;
572 Double_t d0K = 0.;
573 Double_t d0xd0 = d0tokpi->Prodd0d0();
574 Double_t cosPointingAngle = d0tokpi->CosPointingAngle();
575 Double_t phi = d0tokpi->Phi();
576 Int_t pdgCode = mcVtxHF->GetPdgCode();
577 if (pdgCode > 0){
578 cosThetaStar = d0tokpi->CosThetaStarD0();
579 pTpi = d0tokpi->PtProng(0);
580 pTK = d0tokpi->PtProng(1);
581 d0pi = d0tokpi->Getd0Prong(0);
582 d0K = d0tokpi->Getd0Prong(1);
583 invMass=d0tokpi->InvMassD0();
584 }
585 else {
586 cosThetaStar = d0tokpi->CosThetaStarD0bar();
587 pTpi = d0tokpi->PtProng(1);
588 pTK = d0tokpi->PtProng(0);
589 d0pi = d0tokpi->Getd0Prong(1);
590 d0K = d0tokpi->Getd0Prong(0);
591 invMass=d0tokpi->InvMassD0bar();
592 }
593
594 Double_t cT = d0tokpi->CtD0();
595
596 if (!fFillFromGenerated){
597 // ...either with reconstructed values....
598 containerInput[0] = pt;
599 containerInput[1] = rapidity;
600 containerInput[2] = cosThetaStar;
601 containerInput[3] = pTpi;
602 containerInput[4] = pTK;
603 containerInput[5] = cT*1.E4; // in micron
604 containerInput[6] = dca*1.E4; // in micron
605 containerInput[7] = d0pi*1.E4; // in micron
606 containerInput[8] = d0K*1.E4; // in micron
607 containerInput[9] = d0xd0*1.E8; // in micron^2
608 containerInput[10] = cosPointingAngle; // in micron
609 containerInput[11] = phi;
610 containerInput[12] = zPrimVertex; // z of reconstructed of primary vertex
611 }
612 else {
613 // ... or with generated values
614 Double_t vectorMC[7] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.};
615 if (GetGeneratedValuesFromMCParticle(mcVtxHF, mcArray, vectorMC)){
616 containerInput[0] = vectorMC[0];
617 containerInput[1] = vectorMC[1] ;
618 containerInput[2] = vectorMC[2] ;
619 containerInput[3] = vectorMC[3] ;
620 containerInput[4] = vectorMC[4] ;
621 containerInput[5] = vectorMC[5] ; // in micron
622 containerInput[6] = 0.; // dummy value, meaningless in MC, in micron
623 containerInput[7] = 0.; // dummy value, meaningless in MC, in micron
624 containerInput[8] = 0.; // dummy value, meaningless in MC, in micron
625 containerInput[9] = 100000.; // dummy value, meaningless in MC, in micron^2
626 containerInput[10] = 1.01; // dummy value, meaningless in MC
627 containerInput[11] = vectorMC[6];
628 containerInput[12] = zMCVertex; // z of reconstructed of primary vertex
629 }
630 else {
631 AliDebug(3,"Problems in filling the container");
632 continue;
633 }
634 }
635
636 if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue; // fiducial region
637
638 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]));
639 icountReco++;
640 AliDebug(2,Form("%d: filling RECO step\n",iD0toKpi));
641 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed,fWeight) ;
642
643 // cut in acceptance
644 Float_t etaCutMin, etaCutMax, ptCutMin, ptCutMax;
645 trackCuts->GetEtaRange(etaCutMin,etaCutMax);
646 trackCuts->GetPtRange(ptCutMin,ptCutMax);
647 Bool_t acceptanceProng0 = (d0tokpi->EtaProng(0)>etaCutMin && d0tokpi->EtaProng(0)<etaCutMax && d0tokpi->PtProng(0) > ptCutMin && d0tokpi->PtProng(0) < ptCutMax);
648 Bool_t acceptanceProng1 = (d0tokpi->EtaProng(1)>etaCutMin && d0tokpi->EtaProng(1)<etaCutMax && d0tokpi->PtProng(1) > ptCutMin && d0tokpi->PtProng(1) < ptCutMax);
649 if (acceptanceProng0 && acceptanceProng1) {
650 AliDebug(2,"D0 reco daughters in acceptance");
651 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance,fWeight) ;
652 icountRecoAcc++;
653
654 if(fAcceptanceUnf){
655
656 Double_t fill[4]; //fill response matrix
657
658 // dimensions 0&1 : pt,eta (Rec)
659
660 fill[0] = pt ;
661 fill[1] = rapidity;
662
663 // dimensions 2&3 : pt,eta (MC)
664
665 fill[2] = mcVtxHF->Pt();
666 fill[3] = mcVtxHF->Y();
667
668 fCorrelation->Fill(fill);
669
670 }
671
672 // cut on the min n. of clusters in ITS
673 if (fCuts->IsSelected(d0tokpi,AliRDHFCuts::kTracks)){
674 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters) ;
675 icountRecoITSClusters++;
676 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));
677
678 // setting the use of the PID cut when applying the selection to FALSE - whatever it was. Keeping track of the original value
679 Bool_t iscutusingpid=fCuts->GetIsUsePID();
680 Int_t isselcuts=-1,isselpid=-1;
681 fCuts->SetUsePID(kFALSE);
682 //Bool_t origFlag = fCuts->GetIsPrimaryWithoutDaughters();
683 //fCuts->SetRemoveDaughtersFromPrim(kFALSE);
684 isselcuts = fCuts->IsSelected(d0tokpi,AliRDHFCuts::kCandidate,aodEvent);
685 //fCuts->SetRemoveDaughtersFromPrim(origFlag);
686 fCuts->SetUsePID(iscutusingpid); // restoring usage of the PID from the cuts object
687 if (isselcuts == 3 || isselcuts == isD0D0bar){
688 AliDebug(2,"Particle passed PPR cuts (actually cuts for D0 analysis!)");
689 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoPPR,fWeight) ;
690 icountRecoPPR++;
691
692 if(!fAcceptanceUnf){ // unfolding
693
694 Double_t fill[4]; //fill response matrix
695
696 // dimensions 0&1 : pt,eta (Rec)
697
698 fill[0] = pt ;
699 fill[1] = rapidity;
700
701 // dimensions 2&3 : pt,eta (MC)
702
703 fill[2] = mcVtxHF->Pt();
704 fill[3] = mcVtxHF->Y();
705
706 fCorrelation->Fill(fill);
707
708 }
709
710 isselpid = fCuts->IsSelected(d0tokpi,AliRDHFCuts::kPID);
711 if((fCuts->CombineSelectionLevels(3,isselcuts,isselpid)==isD0D0bar)||(fCuts->CombineSelectionLevels(3,isselcuts,isselpid)==3)){
712 AliDebug(2,"Particle passed PID cuts");
713 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoPID,fWeight) ;
714 icountRecoPID++;
715 }
716 }
717 }
718 }
719 }
720 if(unsetvtx) d0tokpi->UnsetOwnPrimaryVtx();
721 } // end loop on D0->Kpi
722
723 AliDebug(2, Form("Found %i Reco particles that are D0!!",icountReco));
724
725 fCountReco+= icountReco;
726 fCountVertex+= icountVertex;
727 fCountRefit+= icountRefit;
728 fCountRecoAcc+= icountRecoAcc;
729 fCountRecoITSClusters+= icountRecoITSClusters;
730 fCountRecoPPR+= icountRecoPPR;
731 fCountRecoPID+= icountRecoPID;
732
733 fHistEventsProcessed->Fill(0);
734 /* PostData(0) is taken care of by AliAnalysisTaskSE */
735 //PostData(1,fHistEventsProcessed) ;
736 //PostData(2,fCFManager->GetParticleContainer()) ;
737 //PostData(3,fCorrelation) ;
738}
739
740
741//___________________________________________________________________________
742void AliCFHeavyFlavourTaskMultiVarMultiStep::Terminate(Option_t*)
743{
744 // The Terminate() function is the last function to be called during
745 // a query. It always runs on the client, it can be used to present
746 // the results graphically or save the results to file.
747
748 AliAnalysisTaskSE::Terminate();
749
750 AliInfo(Form("Found %i MC particles that are D0 in MC, in %d events",fCountMC,fEvents));
751 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, in %d events",fCountAcc,fEvents));
752 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, and satisfy Vertex requirement in %d events",fCountVertex,fEvents));
753 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, and satisfy ITS+TPC refit requirementin %d events",fCountRefit,fEvents));
754 AliInfo(Form("Found %i reco D0 that are decaying in K+pi, in %d events",fCountReco,fEvents));
755 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));
756 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));
757 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and satisfy PPR cuts, in %d events",fCountRecoPPR,fEvents));
758 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and satisfy PID cuts, in %d events",fCountRecoPID,fEvents));
759
760 // draw some example plots....
761
762 // AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
763 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
764 if(!cont) {
765 printf("CONTAINER NOT FOUND\n");
766 return;
767 }
768 // projecting the containers to obtain histograms
769 // first argument = variable, second argument = step
770
771 // MC-level
772 TH1D* h00 = cont->ShowProjection(0,0) ; // pt
773 TH1D* h10 = cont->ShowProjection(1,0) ; // rapidity
774 TH1D* h20 = cont->ShowProjection(2,0) ; // cosThetaStar
775 TH1D* h30 = cont->ShowProjection(3,0) ; // pTpi
776 TH1D* h40 = cont->ShowProjection(4,0) ; // pTK
777 TH1D* h50 = cont->ShowProjection(5,0) ; // cT
778 TH1D* h60 = cont->ShowProjection(6,0) ; // dca
779 TH1D* h70 = cont->ShowProjection(7,0) ; // d0pi
780 TH1D* h80 = cont->ShowProjection(8,0) ; // d0K
781 TH1D* h90 = cont->ShowProjection(9,0) ; // d0xd0
782 TH1D* h100 = cont->ShowProjection(10,0) ; // cosPointingAngle
783 TH1D* h110 = cont->ShowProjection(11,0) ; // phi
784
785 // MC-Acceptance level
786 TH1D* h01 = cont->ShowProjection(0,1) ; // pt
787 TH1D* h11 = cont->ShowProjection(1,1) ; // rapidity
788 TH1D* h21 = cont->ShowProjection(2,1) ; // cosThetaStar
789 TH1D* h31 = cont->ShowProjection(3,1) ; // pTpi
790 TH1D* h41 = cont->ShowProjection(4,1) ; // pTK
791 TH1D* h51 = cont->ShowProjection(5,1) ; // cT
792 TH1D* h61 = cont->ShowProjection(6,1) ; // dca
793 TH1D* h71 = cont->ShowProjection(7,1) ; // d0pi
794 TH1D* h81 = cont->ShowProjection(8,1) ; // d0K
795 TH1D* h91 = cont->ShowProjection(9,1) ; // d0xd0
796 TH1D* h101 = cont->ShowProjection(10,1) ; // cosPointingAngle
797 TH1D* h111 = cont->ShowProjection(11,1) ; // phi
798
799 // Reco-level
800 TH1D* h04 = cont->ShowProjection(0,4) ; // pt
801 TH1D* h14 = cont->ShowProjection(1,4) ; // rapidity
802 TH1D* h24 = cont->ShowProjection(2,4) ; // cosThetaStar
803 TH1D* h34 = cont->ShowProjection(3,4) ; // pTpi
804 TH1D* h44 = cont->ShowProjection(4,4) ; // pTK
805 TH1D* h54 = cont->ShowProjection(5,4) ; // cT
806 TH1D* h64 = cont->ShowProjection(6,4) ; // dca
807 TH1D* h74 = cont->ShowProjection(7,4) ; // d0pi
808 TH1D* h84 = cont->ShowProjection(8,4) ; // d0K
809 TH1D* h94 = cont->ShowProjection(9,4) ; // d0xd0
810 TH1D* h104 = cont->ShowProjection(10,4) ; // cosPointingAngle
811 TH1D* h114 = cont->ShowProjection(11,4) ; // phi
812
813 h00->SetTitle("pT_D0 (GeV/c)");
814 h10->SetTitle("rapidity");
815 h20->SetTitle("cosThetaStar");
816 h30->SetTitle("pT_pi (GeV/c)");
817 h40->SetTitle("pT_K (Gev/c)");
818 h50->SetTitle("cT (#mum)");
819 h60->SetTitle("dca (#mum)");
820 h70->SetTitle("d0_pi (#mum)");
821 h80->SetTitle("d0_K (#mum)");
822 h90->SetTitle("d0xd0 (#mum^2)");
823 h100->SetTitle("cosPointingAngle");
824 h100->SetTitle("phi (rad)");
825
826 h00->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
827 h10->GetXaxis()->SetTitle("rapidity");
828 h20->GetXaxis()->SetTitle("cosThetaStar");
829 h30->GetXaxis()->SetTitle("pT_pi (GeV/c)");
830 h40->GetXaxis()->SetTitle("pT_K (Gev/c)");
831 h50->GetXaxis()->SetTitle("cT (#mum)");
832 h60->GetXaxis()->SetTitle("dca (#mum)");
833 h70->GetXaxis()->SetTitle("d0_pi (#mum)");
834 h80->GetXaxis()->SetTitle("d0_K (#mum)");
835 h90->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
836 h100->GetXaxis()->SetTitle("cosPointingAngle");
837 h110->GetXaxis()->SetTitle("phi (rad)");
838
839 h01->SetTitle("pT_D0 (GeV/c)");
840 h11->SetTitle("rapidity");
841 h21->SetTitle("cosThetaStar");
842 h31->SetTitle("pT_pi (GeV/c)");
843 h41->SetTitle("pT_K (Gev/c)");
844 h51->SetTitle("cT (#mum)");
845 h61->SetTitle("dca (#mum)");
846 h71->SetTitle("d0_pi (#mum)");
847 h81->SetTitle("d0_K (#mum)");
848 h91->SetTitle("d0xd0 (#mum^2)");
849 h101->SetTitle("cosPointingAngle");
850 h111->GetXaxis()->SetTitle("phi (rad)");
851
852 h01->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
853 h11->GetXaxis()->SetTitle("rapidity");
854 h21->GetXaxis()->SetTitle("cosThetaStar");
855 h31->GetXaxis()->SetTitle("pT_pi (GeV/c)");
856 h41->GetXaxis()->SetTitle("pT_K (Gev/c)");
857 h51->GetXaxis()->SetTitle("cT (#mum)");
858 h61->GetXaxis()->SetTitle("dca (#mum)");
859 h71->GetXaxis()->SetTitle("d0_pi (#mum)");
860 h81->GetXaxis()->SetTitle("d0_K (#mum)");
861 h91->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
862 h101->GetXaxis()->SetTitle("cosPointingAngle");
863 h111->GetXaxis()->SetTitle("phi (rad)");
864
865 h04->SetTitle("pT_D0 (GeV/c)");
866 h14->SetTitle("rapidity");
867 h24->SetTitle("cosThetaStar");
868 h34->SetTitle("pT_pi (GeV/c)");
869 h44->SetTitle("pT_K (Gev/c)");
870 h54->SetTitle("cT (#mum)");
871 h64->SetTitle("dca (#mum)");
872 h74->SetTitle("d0_pi (#mum)");
873 h84->SetTitle("d0_K (#mum)");
874 h94->SetTitle("d0xd0 (#mum^2)");
875 h104->SetTitle("cosPointingAngle");
876 h114->GetXaxis()->SetTitle("phi (rad)");
877
878 h04->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
879 h14->GetXaxis()->SetTitle("rapidity");
880 h24->GetXaxis()->SetTitle("cosThetaStar");
881 h34->GetXaxis()->SetTitle("pT_pi (GeV/c)");
882 h44->GetXaxis()->SetTitle("pT_K (Gev/c)");
883 h54->GetXaxis()->SetTitle("cT (#mum)");
884 h64->GetXaxis()->SetTitle("dca (#mum)");
885 h74->GetXaxis()->SetTitle("d0_pi (#mum)");
886 h84->GetXaxis()->SetTitle("d0_K (#mum)");
887 h94->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
888 h104->GetXaxis()->SetTitle("cosPointingAngle");
889 h114->GetXaxis()->SetTitle("phi (rad)");
890
891 Double_t max0 = h00->GetMaximum();
892 Double_t max1 = h10->GetMaximum();
893 Double_t max2 = h20->GetMaximum();
894 Double_t max3 = h30->GetMaximum();
895 Double_t max4 = h40->GetMaximum();
896 Double_t max5 = h50->GetMaximum();
897 Double_t max6 = h60->GetMaximum();
898 Double_t max7 = h70->GetMaximum();
899 Double_t max8 = h80->GetMaximum();
900 Double_t max9 = h90->GetMaximum();
901 Double_t max10 = h100->GetMaximum();
902 Double_t max11 = h110->GetMaximum();
903
904 h00->GetYaxis()->SetRangeUser(0,max0*1.2);
905 h10->GetYaxis()->SetRangeUser(0,max1*1.2);
906 h20->GetYaxis()->SetRangeUser(0,max2*1.2);
907 h30->GetYaxis()->SetRangeUser(0,max3*1.2);
908 h40->GetYaxis()->SetRangeUser(0,max4*1.2);
909 h50->GetYaxis()->SetRangeUser(0,max5*1.2);
910 h60->GetYaxis()->SetRangeUser(0,max6*1.2);
911 h70->GetYaxis()->SetRangeUser(0,max7*1.2);
912 h80->GetYaxis()->SetRangeUser(0,max8*1.2);
913 h90->GetYaxis()->SetRangeUser(0,max9*1.2);
914 h100->GetYaxis()->SetRangeUser(0,max10*1.2);
915 h110->GetYaxis()->SetRangeUser(0,max11*1.2);
916
917 h01->GetYaxis()->SetRangeUser(0,max0*1.2);
918 h11->GetYaxis()->SetRangeUser(0,max1*1.2);
919 h21->GetYaxis()->SetRangeUser(0,max2*1.2);
920 h31->GetYaxis()->SetRangeUser(0,max3*1.2);
921 h41->GetYaxis()->SetRangeUser(0,max4*1.2);
922 h51->GetYaxis()->SetRangeUser(0,max5*1.2);
923 h61->GetYaxis()->SetRangeUser(0,max6*1.2);
924 h71->GetYaxis()->SetRangeUser(0,max7*1.2);
925 h81->GetYaxis()->SetRangeUser(0,max8*1.2);
926 h91->GetYaxis()->SetRangeUser(0,max9*1.2);
927 h101->GetYaxis()->SetRangeUser(0,max10*1.2);
928 h111->GetYaxis()->SetRangeUser(0,max11*1.2);
929
930 h00->SetMarkerStyle(20);
931 h10->SetMarkerStyle(24);
932 h20->SetMarkerStyle(21);
933 h30->SetMarkerStyle(25);
934 h40->SetMarkerStyle(27);
935 h50->SetMarkerStyle(28);
936 h60->SetMarkerStyle(20);
937 h70->SetMarkerStyle(24);
938 h80->SetMarkerStyle(21);
939 h90->SetMarkerStyle(25);
940 h100->SetMarkerStyle(27);
941 h110->SetMarkerStyle(28);
942
943 h00->SetMarkerColor(2);
944 h10->SetMarkerColor(2);
945 h20->SetMarkerColor(2);
946 h30->SetMarkerColor(2);
947 h40->SetMarkerColor(2);
948 h50->SetMarkerColor(2);
949 h60->SetMarkerColor(2);
950 h70->SetMarkerColor(2);
951 h80->SetMarkerColor(2);
952 h90->SetMarkerColor(2);
953 h100->SetMarkerColor(2);
954 h110->SetMarkerColor(2);
955
956 h01->SetMarkerStyle(20) ;
957 h11->SetMarkerStyle(24) ;
958 h21->SetMarkerStyle(21) ;
959 h31->SetMarkerStyle(25) ;
960 h41->SetMarkerStyle(27) ;
961 h51->SetMarkerStyle(28) ;
962 h61->SetMarkerStyle(20);
963 h71->SetMarkerStyle(24);
964 h81->SetMarkerStyle(21);
965 h91->SetMarkerStyle(25);
966 h101->SetMarkerStyle(27);
967 h111->SetMarkerStyle(28);
968
969 h01->SetMarkerColor(8);
970 h11->SetMarkerColor(8);
971 h21->SetMarkerColor(8);
972 h31->SetMarkerColor(8);
973 h41->SetMarkerColor(8);
974 h51->SetMarkerColor(8);
975 h61->SetMarkerColor(8);
976 h71->SetMarkerColor(8);
977 h81->SetMarkerColor(8);
978 h91->SetMarkerColor(8);
979 h101->SetMarkerColor(8);
980 h111->SetMarkerColor(8);
981
982 h04->SetMarkerStyle(20) ;
983 h14->SetMarkerStyle(24) ;
984 h24->SetMarkerStyle(21) ;
985 h34->SetMarkerStyle(25) ;
986 h44->SetMarkerStyle(27) ;
987 h54->SetMarkerStyle(28) ;
988 h64->SetMarkerStyle(20);
989 h74->SetMarkerStyle(24);
990 h84->SetMarkerStyle(21);
991 h94->SetMarkerStyle(25);
992 h104->SetMarkerStyle(27);
993 h114->SetMarkerStyle(28);
994
995 h04->SetMarkerColor(4);
996 h14->SetMarkerColor(4);
997 h24->SetMarkerColor(4);
998 h34->SetMarkerColor(4);
999 h44->SetMarkerColor(4);
1000 h54->SetMarkerColor(4);
1001 h64->SetMarkerColor(4);
1002 h74->SetMarkerColor(4);
1003 h84->SetMarkerColor(4);
1004 h94->SetMarkerColor(4);
1005 h104->SetMarkerColor(4);
1006 h114->SetMarkerColor(4);
1007
1008 gStyle->SetCanvasColor(0);
1009 gStyle->SetFrameFillColor(0);
1010 gStyle->SetTitleFillColor(0);
1011 gStyle->SetStatColor(0);
1012
1013 // drawing in 2 separate canvas for a matter of clearity
1014 TCanvas * c1 =new TCanvas("c1","pT, rapidiy, cosThetaStar",1100,1600);
1015 c1->Divide(3,3);
1016
1017 c1->cd(1);
1018 h00->Draw("p");
1019 c1->cd(1);
1020 c1->cd(2);
1021 h01->Draw("p");
1022 c1->cd(2);
1023 c1->cd(3);
1024 h04->Draw("p");
1025 c1->cd(3);
1026 c1->cd(4);
1027 h10->Draw("p");
1028 c1->cd(4);
1029 c1->cd(5);
1030 h11->Draw("p");
1031 c1->cd(5);
1032 c1->cd(6);
1033 h14->Draw("p");
1034 c1->cd(6);
1035 c1->cd(7);
1036 h20->Draw("p");
1037 c1->cd(7);
1038 c1->cd(8);
1039 h21->Draw("p");
1040 c1->cd(8);
1041 c1->cd(9);
1042 h24->Draw("p");
1043 c1->cd(9);
1044 c1->cd();
1045
1046 TCanvas * c2 =new TCanvas("c2","pTpi, pTK, cT",1100,1600);
1047 c2->Divide(3,3);
1048 c2->cd(1);
1049 h30->Draw("p");
1050 c2->cd(1);
1051 c2->cd(2);
1052 h31->Draw("p");
1053 c2->cd(2);
1054 c2->cd(3);
1055 h34->Draw("p");
1056 c2->cd(3);
1057 c2->cd(4);
1058 h40->Draw("p");
1059 c2->cd(4);
1060 c2->cd(5);
1061 h41->Draw("p");
1062 c2->cd(5);
1063 c2->cd(6);
1064 h44->Draw("p");
1065 c2->cd(6);
1066 c2->cd(7);
1067 h50->Draw("p");
1068 c2->cd(7);
1069 c2->cd(8);
1070 h51->Draw("p");
1071 c2->cd(8);
1072 c2->cd(9);
1073 h54->Draw("p");
1074 c2->cd(9);
1075 c2->cd();
1076
1077 TCanvas * c3 =new TCanvas("c3","dca, d0pi, d0K",1100,1600);
1078 c3->Divide(3,3);
1079 c3->cd(1);
1080 h60->Draw("p");
1081 c3->cd(1);
1082 c3->cd(2);
1083 h61->Draw("p");
1084 c3->cd(2);
1085 c3->cd(3);
1086 h64->Draw("p");
1087 c3->cd(3);
1088 c3->cd(4);
1089 h70->Draw("p");
1090 c3->cd(4);
1091 c3->cd(5);
1092 h71->Draw("p");
1093 c3->cd(5);
1094 c3->cd(6);
1095 h74->Draw("p");
1096 c3->cd(6);
1097 c3->cd(7);
1098 h80->Draw("p");
1099 c3->cd(7);
1100 c3->cd(8);
1101 h81->Draw("p");
1102 c3->cd(8);
1103 c3->cd(9);
1104 h84->Draw("p");
1105 c3->cd(9);
1106 c3->cd();
1107
1108 TCanvas * c4 =new TCanvas("c4","d0xd0, cosPointingAngle, phi",1100,1600);
1109 c4->Divide(3,3);
1110 c4->cd(1);
1111 h90->Draw("p");
1112 c4->cd(1);
1113 c4->cd(2);
1114 h91->Draw("p");
1115 c4->cd(2);
1116 c4->cd(3);
1117 h94->Draw("p");
1118 c4->cd(3);
1119 c4->cd(4);
1120 h100->Draw("p");
1121 c4->cd(4);
1122 c4->cd(5);
1123 h101->Draw("p");
1124 c4->cd(5);
1125 c4->cd(6);
1126 h104->Draw("p");
1127 c4->cd(6);
1128 c4->cd(7);
1129 h110->Draw("p");
1130 c4->cd(7);
1131 c4->cd(8);
1132 h111->Draw("p");
1133 c4->cd(8);
1134 c4->cd(9);
1135 h114->Draw("p");
1136 c4->cd(9);
1137 c4->cd();
1138
1139 THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
1140
1141 TH2D* corr1 =hcorr->Projection(0,2);
1142 TH2D* corr2 = hcorr->Projection(1,3);
1143
1144 TCanvas * c7 =new TCanvas("c7","",800,400);
1145 c7->Divide(2,1);
1146 c7->cd(1);
1147 corr1->Draw("text");
1148 c7->cd(2);
1149 corr2->Draw("text");
1150
1151
1152 TString projectionFilename="CFtaskHFprojection";
1153 if(fKeepD0fromBOnly) projectionFilename+="_KeepD0fromBOnly";
1154 else if(fKeepD0fromB) projectionFilename+="_KeepD0fromB";
1155 projectionFilename+=".root";
1156 TFile* fileProjection = new TFile(projectionFilename.Data(),"RECREATE");
1157
1158 corr1->Write();
1159 corr2->Write();
1160 h00->Write("pT_D0_step0");
1161 h10->Write("rapidity_step0");
1162 h20->Write("cosThetaStar_step0");
1163 h30->Write("pT_pi_step0");
1164 h40->Write("pT_K_step0");
1165 h50->Write("cT_step0");
1166 h60->Write("dca_step0");
1167 h70->Write("d0_pi_step0");
1168 h80->Write("d0_K_step0");
1169 h90->Write("d0xd0_step0");
1170 h100->Write("cosPointingAngle_step0");
1171 h110->Write("phi_step0");
1172
1173 h01->Write("pT_D0_step1");
1174 h11->Write("rapidity_step1");
1175 h21->Write("cosThetaStar_step1");
1176 h31->Write("pT_pi_step1");
1177 h41->Write("pT_K_step1");
1178 h51->Write("cT_step1");
1179 h61->Write("dca_step1");
1180 h71->Write("d0_pi_step1");
1181 h81->Write("d0_K_step1");
1182 h91->Write("d0xd0_step1");
1183 h101->Write("cosPointingAngle_step1");
1184 h111->Write("phi_step1");
1185
1186 h04->Write("pT_D0_step2");
1187 h14->Write("rapidity_step2");
1188 h24->Write("cosThetaStar_step2");
1189 h34->Write("pT_pi_step2");
1190 h44->Write("pT_K_step2");
1191 h54->Write("cT_step2");
1192 h64->Write("dca_step2");
1193 h74->Write("d0_pi_step2");
1194 h80->Write("d0_K_step2");
1195 h94->Write("d0xd0_step2");
1196 h104->Write("cosPointingAngle_step2");
1197 h114->Write("phi_step2");
1198
1199 fileProjection->Close();
1200
1201
1202
1203 /*
1204 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.eps");
1205 c2->SaveAs("Plots/pTpi_pTK_cT.eps");
1206 c3->SaveAs("Plots/dca_d0pi_d0TK.eps");
1207 c4->SaveAs("Plots/d0xd0_cosPointingAngle.eps");
1208
1209 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.gif");
1210 c2->SaveAs("Plots/pTpi_pTK_cT.gif");
1211 c3->SaveAs("Plots/dca_d0pi_d0TK.gif");
1212 c4->SaveAs("Plots/d0xd0_cosPointingAngle.gif");
1213 */
1214}
1215
1216//___________________________________________________________________________
1217void AliCFHeavyFlavourTaskMultiVarMultiStep::UserCreateOutputObjects() {
1218 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1219 //TO BE SET BEFORE THE EXECUTION OF THE TASK
1220 //
1221 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
1222
1223 //slot #1
1224 OpenFile(1);
1225 fHistEventsProcessed = new TH1I("CFHFchist0","",1,0,1) ;
1226}
1227
1228//___________________________________________________________________________
1229Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CosThetaStar(const AliAODMCParticle* mcPart, const AliAODMCParticle* mcPartDaughter0, const AliAODMCParticle* mcPartDaughter1) const {
1230
1231 //
1232 // to calculate cos(ThetaStar) for generated particle
1233 // using the K, since mcPartDaughter0 and mcPartDaughter1 always correspond to K and pi respectively
1234 // (see where the function is called)
1235 //
1236
1237 Int_t pdgvtx = mcPart->GetPdgCode();
1238 /* if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
1239 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1240 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1241 AliInfo(Form("D0, with pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
1242 AliDebug(2,"This is a D0");
1243 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
1244 mcPartDaughter0 = mcPartDaughter1;
1245 mcPartDaughter1 = mcPartdummy;
1246 }
1247 else{
1248 AliInfo("D0bar");
1249 }
1250 */
1251 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1252 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1253 if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
1254 AliDebug(2,"D0");
1255 }
1256 else{
1257 AliDebug(2,"D0bar");
1258 }
1259 if (pdgprong0 == 211){
1260 AliDebug(2,Form("pdgprong0 = %d, pdgprong1 = %d, switching...",pdgprong0,pdgprong1));
1261 const AliAODMCParticle* mcPartdummy = mcPartDaughter0;
1262 mcPartDaughter0 = mcPartDaughter1;
1263 mcPartDaughter1 = mcPartdummy;
1264 pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1265 pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1266 }
1267
1268 AliDebug(2,Form("After checking, pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
1269 Double_t massvtx = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdgvtx))->Mass();
1270 Double_t massp[2];
1271 massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
1272 massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
1273
1274 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);
1275
1276 Double_t px = mcPartDaughter0->Px()+mcPartDaughter1->Px();
1277 Double_t py = mcPartDaughter0->Py()+mcPartDaughter1->Py();
1278 Double_t pz = mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
1279 Double_t p = TMath::Sqrt(px*px+py*py+pz*pz);
1280 Double_t e = TMath::Sqrt(massvtx*massvtx+p*p);
1281 Double_t beta = p/e;
1282 Double_t gamma = e/massvtx;
1283 TVector3 mom(mcPartDaughter0->Px(),mcPartDaughter0->Py(),mcPartDaughter0->Pz());
1284 TVector3 momTot(mcPartDaughter0->Px()+mcPartDaughter1->Px(),mcPartDaughter0->Py()+mcPartDaughter1->Py(),mcPartDaughter0->Pz()+mcPartDaughter1->Pz());
1285
1286 Double_t qlprong = mom.Dot(momTot)/momTot.Mag(); // analog to AliAODRecoDecay::QlProng(0)
1287
1288 AliDebug(2,Form("pStar = %f, beta = %f, gamma = %f, qlprong = %f, massp[0] = %f", pStar, beta, gamma, qlprong, massp[0]));
1289 Double_t cts = (qlprong/gamma-beta*TMath::Sqrt(pStar*pStar+massp[0]*massp[0]))/pStar;
1290 AliDebug(2,Form("cts = %f", cts));
1291 return cts;
1292}
1293//___________________________________________________________________________
1294Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CT(const AliAODMCParticle* mcPart, const AliAODMCParticle* mcPartDaughter0, const AliAODMCParticle* mcPartDaughter1) const {
1295
1296 //
1297 // to calculate cT for generated particle
1298 //
1299
1300 Double_t xmcPart[3] = {0,0,0};
1301 Double_t xdaughter0[3] = {0,0,0};
1302 Double_t xdaughter1[3] = {0,0,0};
1303 mcPart->XvYvZv(xmcPart); // cm
1304 mcPartDaughter0->XvYvZv(xdaughter0); // cm
1305 mcPartDaughter1->XvYvZv(xdaughter1); //cm
1306 Double_t prodVtxD0 = TMath::Sqrt(xmcPart[0]*xmcPart[0]+
1307 xmcPart[1]*xmcPart[1]+
1308 xmcPart[2]*xmcPart[2]);
1309 Double_t prodVtxDaughter0 = TMath::Sqrt(xdaughter0[0]*xdaughter0[0]+
1310 xdaughter0[1]*xdaughter0[1]+
1311 xdaughter0[2]*xdaughter0[2]);
1312 Double_t prodVtxDaughter1 = TMath::Sqrt(xdaughter1[0]*xdaughter1[0]+
1313 xdaughter1[1]*xdaughter1[1]+
1314 xdaughter1[2]*xdaughter1[2]);
1315
1316 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));
1317 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));
1318 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));
1319
1320 Double_t cT0 = TMath::Sqrt((xdaughter0[0]-xmcPart[0])*(xdaughter0[0]-xmcPart[0])+
1321 (xdaughter0[1]-xmcPart[1])*(xdaughter0[1]-xmcPart[1])+
1322 (xdaughter0[2]-xmcPart[2])*(xdaughter0[2]-xmcPart[2]));
1323
1324 Double_t cT1 = TMath::Sqrt((xdaughter1[0]-xmcPart[0])*(xdaughter1[0]-xmcPart[0])+
1325 (xdaughter1[1]-xmcPart[1])*(xdaughter1[1]-xmcPart[1])+
1326 (xdaughter1[2]-xmcPart[2])*(xdaughter1[2]-xmcPart[2]));
1327
1328 if ((cT0 - cT1)>1E-5) {
1329 AliWarning(Form("cT from daughter 0 (%f) different from cT from daughter 1 (%f)! Using cT from daughter 0, but PLEASE, CHECK....",cT0,cT1));
1330 }
1331
1332 // calculating cT from cT0
1333
1334 Double_t mass = TDatabasePDG::Instance()->GetParticle(mcPart->GetPdgCode())->Mass(); // mcPart->GetPdgCode() should return +/- 421 for the D0/D0bar
1335 Double_t p = mcPart-> P();
1336 Double_t cT = cT0*mass/p;
1337 AliDebug(2, Form("cT from daughter 0 = %f (micron)", cT0*1E4));
1338 AliDebug(2, Form("cT from daughter 1 = %f (micron)", cT1*1E4));
1339 AliDebug(2, Form("cT (with mass = %f and p = %f) = %f (micron)", mass, p, cT*1E4));
1340 return cT;
1341}
1342//________________________________________________________________________________
1343Bool_t AliCFHeavyFlavourTaskMultiVarMultiStep::GetGeneratedValuesFromMCParticle(AliAODMCParticle* mcPart, const TClonesArray* mcArray, Double_t* vectorMC) const {
1344
1345 //
1346 // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
1347 //
1348
1349 Bool_t isOk = kFALSE;
1350
1351 // check whether the D0 decays in pi+K
1352 // to be added!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1353 // could use a cut in the CF, but not clear how to define a TDecayChannel
1354 // implemented for the time being as a cut on the number of daughters - see later when
1355 // getting the daughters
1356
1357 // getting the daughters
1358 Int_t daughter0 = mcPart->GetDaughter(0);
1359 Int_t daughter1 = mcPart->GetDaughter(1);
1360 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
1361 if (daughter0 == 0 || daughter1 == 0) {
1362 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
1363 return isOk;
1364 }
1365 if (TMath::Abs(daughter1 - daughter0) != 1) {
1366 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
1367 return isOk;
1368 }
1369 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
1370 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
1371 if (!mcPartDaughter0 || !mcPartDaughter1) {
1372 AliWarning("At least one Daughter Particle not found in tree, skipping");
1373 return isOk;
1374 }
1375 if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==321 &&
1376 TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
1377 !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
1378 TMath::Abs(mcPartDaughter1->GetPdgCode())==321)) {
1379 AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
1380 return isOk;
1381 }
1382
1383 Double_t vtx1[3] = {0,0,0}; // primary vertex
1384 Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0
1385 Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1
1386 mcPart->XvYvZv(vtx1); // cm
1387 // getting vertex from daughters
1388 mcPartDaughter0->XvYvZv(vtx2daughter0); // cm
1389 mcPartDaughter1->XvYvZv(vtx2daughter1); //cm
1390 if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
1391 AliError("Daughters have different secondary vertex, skipping the track");
1392 return isOk;
1393 }
1394 Int_t nprongs = 2;
1395 Short_t charge = 0;
1396 // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
1397 AliAODMCParticle* positiveDaugh = mcPartDaughter0;
1398 AliAODMCParticle* negativeDaugh = mcPartDaughter1;
1399 if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){
1400 // inverting in case the positive daughter is the second one
1401 positiveDaugh = mcPartDaughter1;
1402 negativeDaugh = mcPartDaughter0;
1403 }
1404 // getting the momentum from the daughters
1405 Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};
1406 Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};
1407 Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
1408
1409 Double_t d0[2] = {0.,0.};
1410
1411 AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
1412
1413 Double_t cosThetaStar = 0.;
1414 Double_t cosThetaStarD0 = 0.;
1415 Double_t cosThetaStarD0bar = 0.;
1416 cosThetaStarD0 = decay->CosThetaStar(1,421,211,321);
1417 cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211);
1418 if (mcPart->GetPdgCode() == 421){ // D0
1419 AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1420 cosThetaStar = cosThetaStarD0;
1421 }
1422 else if (mcPart->GetPdgCode() == -421){ // D0bar{
1423 AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1424 cosThetaStar = cosThetaStarD0bar;
1425 }
1426 else{
1427 AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check...");
1428 return vectorMC;
1429 }
1430 if (cosThetaStar < -1 || cosThetaStar > 1) {
1431 AliWarning("Invalid value for cosine Theta star, returning");
1432 return isOk;
1433 }
1434
1435 // calculate cos(Theta*) according to the method implemented herein
1436
1437 Double_t cts = 9999.;
1438 cts = CosThetaStar(mcPart, mcPartDaughter0, mcPartDaughter1);
1439 if (cts < -1 || cts > 1) {
1440 AliWarning("Invalid value for cosine Theta star from AliCFHeavyFlavourTaskMultiVarMultiStep method");
1441 }
1442 if (TMath::Abs(cts - cosThetaStar)>0.001) {
1443 AliError(Form("cosThetaStar automatically calculated different from that manually calculated!!! cosThetaStar = %f, cosThetaStar = %f", cosThetaStar,cts));
1444 }
1445
1446 Double_t cT = decay->Ct(421);
1447
1448 // calculate cT from the method implemented herein
1449 Double_t cT1 = 0.;
1450 cT1 = CT(mcPart, mcPartDaughter0, mcPartDaughter1);
1451
1452 if (TMath::Abs(cT1 - cT)>0.001) {
1453 AliError(Form("cT automatically calculated different from that manually calculated!!! cT = %f, cT1 = %f",cT,cT1));
1454 }
1455
1456 // get the pT of the daughters
1457
1458 Double_t pTpi = 0.;
1459 Double_t pTK = 0.;
1460
1461 if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
1462 pTpi = mcPartDaughter0->Pt();
1463 pTK = mcPartDaughter1->Pt();
1464 }
1465 else {
1466 pTpi = mcPartDaughter1->Pt();
1467 pTK = mcPartDaughter0->Pt();
1468 }
1469
1470 vectorMC[0] = mcPart->Pt();
1471 vectorMC[1] = mcPart->Y() ;
1472 vectorMC[2] = cosThetaStar ;
1473 vectorMC[3] = pTpi ;
1474 vectorMC[4] = pTK ;
1475 vectorMC[5] = cT*1.E4 ; // in micron
1476 vectorMC[6] = mcPart->Phi() ;
1477 isOk = kTRUE;
1478 return isOk;
1479}
1480//_________________________________________________________________________________________________
1481Int_t AliCFHeavyFlavourTaskMultiVarMultiStep::CheckOrigin(const AliAODMCParticle* mcPart, const TClonesArray* mcArray)const{
1482
1483 //
1484 // checking whether the very mother of the D0 is a charm or a bottom quark
1485 //
1486
1487 Int_t pdgGranma = 0;
1488 Int_t mother = 0;
1489 mother = mcPart->GetMother();
1490 Int_t istep = 0;
1491 while (mother >0 ){
1492 istep++;
1493 AliDebug(2,Form("mother at step %d = %d", istep, mother));
1494 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
1495 if(!mcGranma) break;
1496 pdgGranma = mcGranma->GetPdgCode();
1497 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
1498 Int_t abspdgGranma = TMath::Abs(pdgGranma);
1499 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1500 break;
1501 }
1502 mother = mcGranma->GetMother();
1503 }
1504 return pdgGranma;
1505}
1506//__________________________________________________________________________________________________
1507Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::GetWeight(Float_t pt){
1508
1509 //
1510 // calculating the weight to fill the container
1511 //
1512
1513 // FNOLL central:
1514 // p0 = 1.63297e-01 --> 0.322643
1515 // p1 = 2.96275e+00
1516 // p2 = 2.30301e+00
1517 // p3 = 2.50000e+00
1518
1519 // PYTHIA
1520 // p0 = 1.85906e-01 --> 0.36609
1521 // p1 = 1.94635e+00
1522 // p2 = 1.40463e+00
1523 // p3 = 2.50000e+00
1524
1525 Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
1526 Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
1527
1528 Double_t dndptFunc1 = DodNdptFit(pt,func1);
1529 Double_t dndptFunc2 = DodNdptFit(pt,func2);
1530 AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndptFunc1,dndptFunc2,dndptFunc1/dndptFunc2));
1531 return dndptFunc1/dndptFunc2;
1532}
1533
1534//__________________________________________________________________________________________________
1535Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::DodNdptFit(Float_t pt, const Double_t* par)const{
1536
1537 //
1538 // calculating dNdpt
1539 //
1540
1541 Double_t denom = TMath::Power((pt/par[1]), par[3] );
1542 Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
1543
1544 return dNdpt;
1545}