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