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