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