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