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