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