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