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