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