From 5cd139bc06308be319b53ff2c65f557a75120330 Mon Sep 17 00:00:00 2001 From: prino Date: Fri, 28 Sep 2012 12:28:11 +0000 Subject: [PATCH] Code for computing MC efficiency for Lc->K0s+proton (Annalisa) --- PWGHF/CMakelibPWGHFvertexingHF.pkg | 1 + PWGHF/PWGHFvertexingHFLinkDef.h | 1 + PWGHF/vertexingHF/AliCFTaskVertexingHF.cxx | 2038 +++++++++-------- PWGHF/vertexingHF/AliCFTaskVertexingHF.h | 18 +- PWGHF/vertexingHF/AliCFVertexingHF.cxx | 34 + PWGHF/vertexingHF/AliCFVertexingHF.h | 2 +- .../AliCFVertexingHFLctoV0bachelor.cxx | 726 ++++++ .../AliCFVertexingHFLctoV0bachelor.h | 83 + .../AddTaskCFVertexingHFLctoV0bachelor.C | 691 ++++++ 9 files changed, 2646 insertions(+), 948 deletions(-) create mode 100644 PWGHF/vertexingHF/AliCFVertexingHFLctoV0bachelor.cxx create mode 100644 PWGHF/vertexingHF/AliCFVertexingHFLctoV0bachelor.h create mode 100644 PWGHF/vertexingHF/macros/AddTaskCFVertexingHFLctoV0bachelor.C diff --git a/PWGHF/CMakelibPWGHFvertexingHF.pkg b/PWGHF/CMakelibPWGHFvertexingHF.pkg index 8574afe9f36..7148443292c 100644 --- a/PWGHF/CMakelibPWGHFvertexingHF.pkg +++ b/PWGHF/CMakelibPWGHFvertexingHF.pkg @@ -62,6 +62,7 @@ set ( SRCS vertexingHF/AliCFVertexingHF2Prong.cxx vertexingHF/AliCFVertexingHF3Prong.cxx vertexingHF/AliCFVertexingHFCascade.cxx + vertexingHF/AliCFVertexingHFLctoV0bachelor.cxx vertexingHF/AliCFTaskVertexingHF.cxx vertexingHF/AliCFTaskForDStarAnalysis.cxx vertexingHF/AliAnalysisTaskSEDStarJets.cxx diff --git a/PWGHF/PWGHFvertexingHFLinkDef.h b/PWGHF/PWGHFvertexingHFLinkDef.h index f1c50f2e615..740682a7daf 100644 --- a/PWGHF/PWGHFvertexingHFLinkDef.h +++ b/PWGHF/PWGHFvertexingHFLinkDef.h @@ -40,6 +40,7 @@ #pragma link C++ class AliCFVertexingHF2Prong+; #pragma link C++ class AliCFVertexingHF3Prong+; #pragma link C++ class AliCFVertexingHFCascade+; +#pragma link C++ class AliCFVertexingHFLctoV0bachelor+; #pragma link C++ class AliCFTaskVertexingHF+; #pragma link C++ class AliCFTaskForDStarAnalysis+; #pragma link C++ class AliAnalysisTaskSEDStarJets+; diff --git a/PWGHF/vertexingHF/AliCFTaskVertexingHF.cxx b/PWGHF/vertexingHF/AliCFTaskVertexingHF.cxx index 6028937fba9..d9712b70449 100644 --- a/PWGHF/vertexingHF/AliCFTaskVertexingHF.cxx +++ b/PWGHF/vertexingHF/AliCFTaskVertexingHF.cxx @@ -43,6 +43,7 @@ #include "AliCFManager.h" #include "AliCFContainer.h" #include "AliLog.h" +#include "AliInputEventHandler.h" #include "AliAnalysisManager.h" #include "AliAODHandler.h" #include "AliAODEvent.h" @@ -66,840 +67,895 @@ #include "AliRDHFCutsDstoKKpi.h" #include "AliRDHFCutsLctopKpi.h" #include "AliRDHFCutsD0toKpipipi.h" +#include "AliRDHFCutsLctoV0.h" #include "AliCFVertexingHF2Prong.h" #include "AliCFVertexingHF3Prong.h" #include "AliCFVertexingHFCascade.h" +#include "AliCFVertexingHFLctoV0bachelor.h" #include "AliCFVertexingHF.h" #include "AliVertexingHFUtils.h" #include "AliAnalysisDataSlot.h" #include "AliAnalysisDataContainer.h" +#include "AliPIDResponse.h" //__________________________________________________________________________ AliCFTaskVertexingHF::AliCFTaskVertexingHF() : - AliAnalysisTaskSE(), - fCFManager(0x0), - fHistEventsProcessed(0x0), - fCorrelation(0x0), - fCountMC(0), - fCountAcc(0), - fCountVertex(0), - fCountRefit(0), - fCountReco(0), - fCountRecoAcc(0), - fCountRecoITSClusters(0), - fCountRecoPPR(0), - fCountRecoPID(0), - fEvents(0), - fDecayChannel(0), - fFillFromGenerated(kFALSE), - fOriginDselection(0), - fAcceptanceUnf(kTRUE), - fCuts(0), - fUseWeight(kFALSE), - fWeight(1.), - fUseFlatPtWeight(kFALSE), - fUseZWeight(kFALSE), - fUseNchWeight(kFALSE), - fNvar(0), - fPartName(""), - fDauNames(""), - fSign(2), - fCentralitySelection(kTRUE), - fFakeSelection(0), - fRejectIfNoQuark(kTRUE), - fUseMCVertex(kFALSE), - fDsOption(1), - fGenDsOption(3), - fConfiguration(kCheetah), // by default, setting the fast configuration - fFuncWeight(0x0), - fHistoMeasNch(0x0), - fHistoMCNch(0x0), - fResonantDecay(0) + AliAnalysisTaskSE(), + fCFManager(0x0), + fHistEventsProcessed(0x0), + fCorrelation(0x0), + fCountMC(0), + fCountAcc(0), + fCountVertex(0), + fCountRefit(0), + fCountReco(0), + fCountRecoAcc(0), + fCountRecoITSClusters(0), + fCountRecoPPR(0), + fCountRecoPID(0), + fEvents(0), + fDecayChannel(0), + fFillFromGenerated(kFALSE), + fOriginDselection(0), + fAcceptanceUnf(kTRUE), + fCuts(0), + fUseWeight(kFALSE), + fWeight(1.), + fUseFlatPtWeight(kFALSE), + fUseZWeight(kFALSE), + fUseNchWeight(kFALSE), + fNvar(0), + fPartName(""), + fDauNames(""), + fSign(2), + fCentralitySelection(kTRUE), + fFakeSelection(0), + fRejectIfNoQuark(kTRUE), + fUseMCVertex(kFALSE), + fDsOption(1), + fGenDsOption(3), + fConfiguration(kCheetah), // by default, setting the fast configuration + fFuncWeight(0x0), + fHistoMeasNch(0x0), + fHistoMCNch(0x0), + fResonantDecay(0), + fLctoV0bachelorOption(1), + fGenLctoV0bachelorOption(0) { - // - //Default ctor - // + // + //Default ctor + // } //___________________________________________________________________________ AliCFTaskVertexingHF::AliCFTaskVertexingHF(const Char_t* name, AliRDHFCuts* cuts, TF1* func) : - AliAnalysisTaskSE(name), - fCFManager(0x0), - fHistEventsProcessed(0x0), - fCorrelation(0x0), - fCountMC(0), - fCountAcc(0), - fCountVertex(0), - fCountRefit(0), - fCountReco(0), - fCountRecoAcc(0), - fCountRecoITSClusters(0), - fCountRecoPPR(0), - fCountRecoPID(0), - fEvents(0), - fDecayChannel(0), - fFillFromGenerated(kFALSE), - fOriginDselection(0), - fAcceptanceUnf(kTRUE), - fCuts(cuts), - fUseWeight(kFALSE), - fWeight(1.), - fUseFlatPtWeight(kFALSE), - fUseZWeight(kFALSE), - fUseNchWeight(kFALSE), - fNvar(0), - fPartName(""), - fDauNames(""), - fSign(2), - fCentralitySelection(kTRUE), - fFakeSelection(0), - fRejectIfNoQuark(kTRUE), - fUseMCVertex(kFALSE), - fDsOption(1), - fGenDsOption(3), - fConfiguration(kCheetah), // by default, setting the fast configuration - fFuncWeight(func), - fHistoMeasNch(0x0), - fHistoMCNch(0x0), - fResonantDecay(0) + AliAnalysisTaskSE(name), + fCFManager(0x0), + fHistEventsProcessed(0x0), + fCorrelation(0x0), + fCountMC(0), + fCountAcc(0), + fCountVertex(0), + fCountRefit(0), + fCountReco(0), + fCountRecoAcc(0), + fCountRecoITSClusters(0), + fCountRecoPPR(0), + fCountRecoPID(0), + fEvents(0), + fDecayChannel(0), + fFillFromGenerated(kFALSE), + fOriginDselection(0), + fAcceptanceUnf(kTRUE), + fCuts(cuts), + fUseWeight(kFALSE), + fWeight(1.), + fUseFlatPtWeight(kFALSE), + fUseZWeight(kFALSE), + fUseNchWeight(kFALSE), + fNvar(0), + fPartName(""), + fDauNames(""), + fSign(2), + fCentralitySelection(kTRUE), + fFakeSelection(0), + fRejectIfNoQuark(kTRUE), + fUseMCVertex(kFALSE), + fDsOption(1), + fGenDsOption(3), + fConfiguration(kCheetah), // by default, setting the fast configuration + fFuncWeight(func), + fHistoMeasNch(0x0), + fHistoMCNch(0x0), + fResonantDecay(0), + fLctoV0bachelorOption(1), + fGenLctoV0bachelorOption(0) { - // - // Constructor. Initialization of Inputs and Outputs - // - /* - DefineInput(0) and DefineOutput(0) - are taken care of by AliAnalysisTaskSE constructor - */ - DefineOutput(1,TH1I::Class()); - DefineOutput(2,AliCFContainer::Class()); - DefineOutput(3,THnSparseD::Class()); - DefineOutput(4,AliRDHFCuts::Class()); + // + // Constructor. Initialization of Inputs and Outputs + // + /* + DefineInput(0) and DefineOutput(0) + are taken care of by AliAnalysisTaskSE constructor + */ + DefineOutput(1,TH1I::Class()); + DefineOutput(2,AliCFContainer::Class()); + DefineOutput(3,THnSparseD::Class()); + DefineOutput(4,AliRDHFCuts::Class()); - fCuts->PrintAll(); + fCuts->PrintAll(); } //___________________________________________________________________________ AliCFTaskVertexingHF& AliCFTaskVertexingHF::operator=(const AliCFTaskVertexingHF& c) { - // - // Assignment operator - // - if (this!=&c) { - AliAnalysisTaskSE::operator=(c) ; - fCFManager = c.fCFManager; - fHistEventsProcessed = c.fHistEventsProcessed; - fCuts = c.fCuts; - fFuncWeight = c.fFuncWeight; - fHistoMeasNch = c.fHistoMeasNch; - fHistoMCNch = c.fHistoMCNch; - } - return *this; + // + // Assignment operator + // + if (this!=&c) { + AliAnalysisTaskSE::operator=(c) ; + fCFManager = c.fCFManager; + fHistEventsProcessed = c.fHistEventsProcessed; + fCuts = c.fCuts; + fFuncWeight = c.fFuncWeight; + fHistoMeasNch = c.fHistoMeasNch; + fHistoMCNch = c.fHistoMCNch; + } + return *this; } //___________________________________________________________________________ AliCFTaskVertexingHF::AliCFTaskVertexingHF(const AliCFTaskVertexingHF& c) : - AliAnalysisTaskSE(c), - fCFManager(c.fCFManager), - fHistEventsProcessed(c.fHistEventsProcessed), - fCorrelation(c.fCorrelation), - fCountMC(c.fCountMC), - fCountAcc(c.fCountAcc), - fCountVertex(c.fCountVertex), - fCountRefit(c.fCountRefit), - fCountReco(c.fCountReco), - fCountRecoAcc(c.fCountRecoAcc), - fCountRecoITSClusters(c.fCountRecoITSClusters), - fCountRecoPPR(c.fCountRecoPPR), - fCountRecoPID(c.fCountRecoPID), - fEvents(c.fEvents), - fDecayChannel(c.fDecayChannel), - fFillFromGenerated(c.fFillFromGenerated), - fOriginDselection(c.fOriginDselection), - fAcceptanceUnf(c.fAcceptanceUnf), - fCuts(c.fCuts), - fUseWeight(c.fUseWeight), - fWeight(c.fWeight), - fUseFlatPtWeight(c.fUseFlatPtWeight), - fUseZWeight(c.fUseZWeight), - fUseNchWeight(c.fUseNchWeight), - fNvar(c.fNvar), - fPartName(c.fPartName), - fDauNames(c.fDauNames), - fSign(c.fSign), - fCentralitySelection(c.fCentralitySelection), - fFakeSelection(c.fFakeSelection), - fRejectIfNoQuark(c.fRejectIfNoQuark), - fUseMCVertex(c.fUseMCVertex), - fDsOption(c.fDsOption), - fGenDsOption(c.fGenDsOption), - fConfiguration(c.fConfiguration), - fFuncWeight(c.fFuncWeight), - fHistoMeasNch(c.fHistoMeasNch), - fHistoMCNch(c.fHistoMCNch), - fResonantDecay(c.fResonantDecay) + AliAnalysisTaskSE(c), + fCFManager(c.fCFManager), + fHistEventsProcessed(c.fHistEventsProcessed), + fCorrelation(c.fCorrelation), + fCountMC(c.fCountMC), + fCountAcc(c.fCountAcc), + fCountVertex(c.fCountVertex), + fCountRefit(c.fCountRefit), + fCountReco(c.fCountReco), + fCountRecoAcc(c.fCountRecoAcc), + fCountRecoITSClusters(c.fCountRecoITSClusters), + fCountRecoPPR(c.fCountRecoPPR), + fCountRecoPID(c.fCountRecoPID), + fEvents(c.fEvents), + fDecayChannel(c.fDecayChannel), + fFillFromGenerated(c.fFillFromGenerated), + fOriginDselection(c.fOriginDselection), + fAcceptanceUnf(c.fAcceptanceUnf), + fCuts(c.fCuts), + fUseWeight(c.fUseWeight), + fWeight(c.fWeight), + fUseFlatPtWeight(c.fUseFlatPtWeight), + fUseZWeight(c.fUseZWeight), + fUseNchWeight(c.fUseNchWeight), + fNvar(c.fNvar), + fPartName(c.fPartName), + fDauNames(c.fDauNames), + fSign(c.fSign), + fCentralitySelection(c.fCentralitySelection), + fFakeSelection(c.fFakeSelection), + fRejectIfNoQuark(c.fRejectIfNoQuark), + fUseMCVertex(c.fUseMCVertex), + fDsOption(c.fDsOption), + fGenDsOption(c.fGenDsOption), + fConfiguration(c.fConfiguration), + fFuncWeight(c.fFuncWeight), + fHistoMeasNch(c.fHistoMeasNch), + fHistoMCNch(c.fHistoMCNch), + fResonantDecay(c.fResonantDecay), + fLctoV0bachelorOption(c.fLctoV0bachelorOption), + fGenLctoV0bachelorOption(c.fGenLctoV0bachelorOption) { - // - // Copy Constructor - // + // + // Copy Constructor + // } //___________________________________________________________________________ AliCFTaskVertexingHF::~AliCFTaskVertexingHF() { - // - //destructor - // - if (fCFManager) delete fCFManager ; - if (fHistEventsProcessed) delete fHistEventsProcessed ; - if (fCorrelation) delete fCorrelation ; - if (fCuts) delete fCuts; - if (fFuncWeight) delete fFuncWeight; - if (fHistoMeasNch) delete fHistoMeasNch; - if (fHistoMCNch) delete fHistoMCNch; + // + //destructor + // + if (fCFManager) delete fCFManager ; + if (fHistEventsProcessed) delete fHistEventsProcessed ; + if (fCorrelation) delete fCorrelation ; + if (fCuts) delete fCuts; + if (fFuncWeight) delete fFuncWeight; + if (fHistoMeasNch) delete fHistoMeasNch; + if (fHistoMCNch) delete fHistoMCNch; } //_________________________________________________________________________- void AliCFTaskVertexingHF::Init() { - // - // Initialization - // + // + // Initialization + // - if (fDebug>1) printf("AliCFTaskVertexingHF::Init()"); - if(fUseWeight && fUseZWeight) { AliFatal("Can not use at the same time pt and z-vtx weights, please choose"); return; } - if(fUseWeight && fUseNchWeight) { AliFatal("Can not use at the same time pt and Nch weights, please choose"); return; } - if(fUseNchWeight && !fHistoMCNch) { AliFatal("Need to pass the MC Nch distribution to use Nch weights"); return; } - if(fUseNchWeight) CreateMeasuredNchHisto(); - - AliRDHFCuts *copyfCuts = 0x0; - if (!fCuts){ - AliFatal("No cuts defined - Exiting..."); - return; - } + if (fDebug>1) printf("AliCFTaskVertexingHF::Init()"); + if(fUseWeight && fUseZWeight) { AliFatal("Can not use at the same time pt and z-vtx weights, please choose"); return; } + if(fUseWeight && fUseNchWeight) { AliFatal("Can not use at the same time pt and Nch weights, please choose"); return; } + if(fUseNchWeight && !fHistoMCNch) { AliFatal("Need to pass the MC Nch distribution to use Nch weights"); return; } + if(fUseNchWeight) CreateMeasuredNchHisto(); - switch (fDecayChannel){ - case 2:{ - copyfCuts = new AliRDHFCutsD0toKpi(*(static_cast(fCuts))); - switch (fConfiguration) { - case kSnail: // slow configuration: all variables in - fNvar = 16; - break; - case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled - fNvar = 8; - break; - } - fPartName="D0"; - fDauNames="K+pi"; - break; - } - case 21:{ - copyfCuts = new AliRDHFCutsDStartoKpipi(*(static_cast(fCuts))); - switch (fConfiguration) { - case kSnail: // slow configuration: all variables in - fNvar = 16; - break; - case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled - fNvar = 8; - break; - } - fPartName="Dstar"; - fDauNames="K+pi+pi"; - break; - } - case 31:{ - copyfCuts = new AliRDHFCutsDplustoKpipi(*(static_cast(fCuts))); - switch (fConfiguration) { - case kSnail: // slow configuration: all variables in - fNvar = 14; - break; - case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled - fNvar = 8; - break; - } - fPartName="Dplus"; - fDauNames="K+pi+pi"; - break; - } - case 32:{ - copyfCuts = new AliRDHFCutsLctopKpi(*(static_cast(fCuts))); - switch (fConfiguration) { - case kSnail: // slow configuration: all variables in - fNvar = 18; - break; - case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled - fNvar = 8; - break; - } - fPartName="Lambdac"; - fDauNames="p+K+pi"; - break; - } - case 33:{ - copyfCuts = new AliRDHFCutsDstoKKpi(*(static_cast(fCuts))); - switch (fConfiguration) { - case kSnail: // slow configuration: all variables in - fNvar = 14; - break; - case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled - fNvar = 8; - break; - } - fPartName="Ds"; - fDauNames="K+K+pi"; - break; - } - case 4:{ - copyfCuts = new AliRDHFCutsD0toKpipipi(*(static_cast(fCuts))); - switch (fConfiguration) { - case kSnail: // slow configuration: all variables in - fNvar = 16; - break; - case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled - fNvar = 8; - break; - } - fPartName="D0"; - fDauNames="K+pi+pi+pi"; - break; - } - default: - AliFatal("The decay channel MUST be defined according to AliCFVertexing::DecayChannel - Exiting..."); - break; - } + AliRDHFCuts *copyfCuts = 0x0; + if (!fCuts){ + AliFatal("No cuts defined - Exiting..."); + return; + } + + switch (fDecayChannel){ + case 2:{ + copyfCuts = new AliRDHFCutsD0toKpi(*(static_cast(fCuts))); + switch (fConfiguration) { + case kSnail: // slow configuration: all variables in + fNvar = 16; + break; + case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled + fNvar = 8; + break; + } + fPartName="D0"; + fDauNames="K+pi"; + break; + } + case 21:{ + copyfCuts = new AliRDHFCutsDStartoKpipi(*(static_cast(fCuts))); + switch (fConfiguration) { + case kSnail: // slow configuration: all variables in + fNvar = 16; + break; + case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled + fNvar = 8; + break; + } + fPartName="Dstar"; + fDauNames="K+pi+pi"; + break; + } + case 22:{ + copyfCuts = new AliRDHFCutsLctoV0(*(static_cast(fCuts))); + switch (fConfiguration) { + case kSnail: // slow configuration: all variables in + fNvar = 16; + break; + case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled + fNvar = 8; + break; + } + fPartName="Lambdac"; + fDauNames="V0+bachelor"; + break; + } + case 31:{ + copyfCuts = new AliRDHFCutsDplustoKpipi(*(static_cast(fCuts))); + switch (fConfiguration) { + case kSnail: // slow configuration: all variables in + fNvar = 14; + break; + case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled + fNvar = 8; + break; + } + fPartName="Dplus"; + fDauNames="K+pi+pi"; + break; + } + case 32:{ + copyfCuts = new AliRDHFCutsLctopKpi(*(static_cast(fCuts))); + switch (fConfiguration) { + case kSnail: // slow configuration: all variables in + fNvar = 18; + break; + case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled + fNvar = 8; + break; + } + fPartName="Lambdac"; + fDauNames="p+K+pi"; + break; + } + case 33:{ + copyfCuts = new AliRDHFCutsDstoKKpi(*(static_cast(fCuts))); + switch (fConfiguration) { + case kSnail: // slow configuration: all variables in + fNvar = 14; + break; + case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled + fNvar = 8; + break; + } + fPartName="Ds"; + fDauNames="K+K+pi"; + break; + } + case 4:{ + copyfCuts = new AliRDHFCutsD0toKpipipi(*(static_cast(fCuts))); + switch (fConfiguration) { + case kSnail: // slow configuration: all variables in + fNvar = 16; + break; + case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled + fNvar = 8; + break; + } + fPartName="D0"; + fDauNames="K+pi+pi+pi"; + break; + } + default: + AliFatal("The decay channel MUST be defined according to AliCFVertexing::DecayChannel - Exiting..."); + break; + } - const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName(); - if (copyfCuts){ - copyfCuts->SetName(nameoutput); + const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName(); + if (copyfCuts){ + copyfCuts->SetName(nameoutput); - //Post the data - PostData(4, copyfCuts); - } - else{ - AliFatal("Failing initializing AliRDHFCuts object - Exiting..."); - } + //Post the data + PostData(4, copyfCuts); + } + else{ + AliFatal("Failing initializing AliRDHFCuts object - Exiting..."); + } - return; + return; } //_________________________________________________ void AliCFTaskVertexingHF::UserExec(Option_t *) { - // - // Main loop function - // + // + // Main loop function + // - PostData(1,fHistEventsProcessed) ; - PostData(2,fCFManager->GetParticleContainer()) ; - PostData(3,fCorrelation) ; + PostData(1,fHistEventsProcessed) ; + PostData(2,fCFManager->GetParticleContainer()) ; + PostData(3,fCorrelation) ; - if (fFillFromGenerated){ - AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!"); - } + if (fFillFromGenerated){ + AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!"); + } - if (!fInputEvent) { - Error("UserExec","NO EVENT FOUND!"); - return; - } + if (!fInputEvent) { + Error("UserExec","NO EVENT FOUND!"); + return; + } - AliAODEvent* aodEvent = dynamic_cast(fInputEvent); + AliAODEvent* aodEvent = dynamic_cast(fInputEvent); - TClonesArray *arrayBranch=0; + TClonesArray *arrayBranch=0; - if(!aodEvent && AODEvent() && IsStandardAOD()) { - // In case there is an AOD handler writing a standard AOD, use the AOD - // event in memory rather than the input (ESD) event. - aodEvent = dynamic_cast (AODEvent()); - // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root) - // have to taken from the AOD event hold by the AliAODExtension - AliAODHandler* aodHandler = (AliAODHandler*) - ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler()); - if(aodHandler->GetExtensions()) { - AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root"); - AliAODEvent *aodFromExt = ext->GetAOD(); + if(!aodEvent && AODEvent() && IsStandardAOD()) { + // In case there is an AOD handler writing a standard AOD, use the AOD + // event in memory rather than the input (ESD) event. + aodEvent = dynamic_cast (AODEvent()); + // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root) + // have to taken from the AOD event hold by the AliAODExtension + AliAODHandler* aodHandler = (AliAODHandler*) + ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler()); + if(aodHandler->GetExtensions()) { + AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root"); + AliAODEvent *aodFromExt = ext->GetAOD(); - switch (fDecayChannel){ - case 2:{ - arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi"); - break; - } - case 21:{ - arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar"); - break; - } - case 31: - case 32: - case 33:{ - arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong"); - break; - } - case 4:{ - arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong"); - break; - } - default: - break; - } - } - } - else { - switch (fDecayChannel){ - case 2:{ - arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi"); - break; - } - case 21:{ - arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar"); - break; - } - case 31: - case 32: - case 33:{ - arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong"); - break; - } - case 4:{ - arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm4Prong"); - break; - } - default: - break; - } - } + switch (fDecayChannel){ + case 2:{ + arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi"); + break; + } + case 21:{ + arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar"); + break; + } + case 22:{ + arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("CascadesHF"); + break; + } + case 31: + case 32: + case 33:{ + arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong"); + break; + } + case 4:{ + arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong"); + break; + } + default: + break; + } + } + } + else { + switch (fDecayChannel){ + case 2:{ + arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi"); + break; + } + case 21:{ + arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar"); + break; + } + case 22:{ + arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("CascadesHF"); + break; + } + case 31: + case 32: + case 33:{ + arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong"); + break; + } + case 4:{ + arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm4Prong"); + break; + } + default: + break; + } + } - AliAODVertex *aodVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex(); - if (!aodVtx) return; + AliAODVertex *aodVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex(); + if (!aodVtx) return; - if (!arrayBranch) { - AliError("Could not find array of HF vertices"); - return; - } + if (!arrayBranch) { + AliError("Could not find array of HF vertices"); + return; + } - fEvents++; + fEvents++; - fCFManager->SetRecEventInfo(aodEvent); - fCFManager->SetMCEventInfo(aodEvent); + fCFManager->SetRecEventInfo(aodEvent); + fCFManager->SetMCEventInfo(aodEvent); - //******** DEFINE number of variables of the container***** for now set at 13, in the future in the config macro. + //******** DEFINE number of variables of the container***** for now set at 13, in the future in the config macro. - //loop on the MC event + //loop on the MC event - TClonesArray* mcArray = dynamic_cast(aodEvent->FindListObject(AliAODMCParticle::StdBranchName())); - if (!mcArray) { - AliError("Could not find Monte-Carlo in AOD"); - return; - } - Int_t icountMC = 0; - Int_t icountAcc = 0; - Int_t icountReco = 0; - Int_t icountVertex = 0; - Int_t icountRefit = 0; - Int_t icountRecoAcc = 0; - Int_t icountRecoITSClusters = 0; - Int_t icountRecoPPR = 0; - Int_t icountRecoPID = 0; - Int_t cquarks = 0; + TClonesArray* mcArray = dynamic_cast(aodEvent->FindListObject(AliAODMCParticle::StdBranchName())); + if (!mcArray) { + AliError("Could not find Monte-Carlo in AOD"); + return; + } + Int_t icountMC = 0; + Int_t icountAcc = 0; + Int_t icountReco = 0; + Int_t icountVertex = 0; + Int_t icountRefit = 0; + Int_t icountRecoAcc = 0; + Int_t icountRecoITSClusters = 0; + Int_t icountRecoPPR = 0; + Int_t icountRecoPID = 0; + Int_t cquarks = 0; - AliAODMCHeader *mcHeader = dynamic_cast(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName())); - if (!mcHeader) { - AliError("Could not find MC Header in AOD"); - return; - } + AliAODMCHeader *mcHeader = dynamic_cast(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName())); + if (!mcHeader) { + AliError("Could not find MC Header in AOD"); + return; + } - Double_t* containerInput = new Double_t[fNvar]; - Double_t* containerInputMC = new Double_t[fNvar]; + Double_t* containerInput = new Double_t[fNvar]; + Double_t* containerInputMC = new Double_t[fNvar]; - AliCFVertexingHF* cfVtxHF=0x0; - switch (fDecayChannel){ - case 2:{ - cfVtxHF = new AliCFVertexingHF2Prong(mcArray, fOriginDselection); - break; - } - case 21:{ - cfVtxHF = new AliCFVertexingHFCascade(mcArray, fOriginDselection); - break; - } - case 31: -// case 32: - case 33:{ - cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel); - if(fDecayChannel==33){ - ((AliCFVertexingHF3Prong*)cfVtxHF)->SetGeneratedDsOption(fGenDsOption); - } - break; - } - case 32:{ - cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel,fResonantDecay); - } - case 4:{ - //cfVtxHF = new AliCFVertexingHF4Prong(mcArray, originDselection); // not there yet - break; - } - default: - break; - } - if (!cfVtxHF){ - AliError("No AliCFVertexingHF initialized"); - delete[] containerInput; - delete[] containerInputMC; - return; - } + AliCFVertexingHF* cfVtxHF=0x0; + switch (fDecayChannel){ + case 2:{ + cfVtxHF = new AliCFVertexingHF2Prong(mcArray, fOriginDselection); + break; + } + case 21:{ + cfVtxHF = new AliCFVertexingHFCascade(mcArray, fOriginDselection); + break; + } + case 22:{ + cfVtxHF = new AliCFVertexingHFLctoV0bachelor(mcArray, fOriginDselection,fGenLctoV0bachelorOption); // Lc -> K0S+proton + break; + } + case 31: + // case 32: + case 33:{ + cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel); + if(fDecayChannel==33){ + ((AliCFVertexingHF3Prong*)cfVtxHF)->SetGeneratedDsOption(fGenDsOption); + } + break; + } + case 32:{ + cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel,fResonantDecay); + } + case 4:{ + //cfVtxHF = new AliCFVertexingHF4Prong(mcArray, originDselection); // not there yet + break; + } + default: + break; + } + if (!cfVtxHF){ + AliError("No AliCFVertexingHF initialized"); + delete[] containerInput; + delete[] containerInputMC; + return; + } - Double_t zPrimVertex = aodVtx ->GetZ(); - Double_t zMCVertex = mcHeader->GetVtxZ(); - Int_t runnumber = aodEvent->GetRunNumber(); - fWeight=1.; - if(fUseZWeight) fWeight *= GetZWeight(zMCVertex,runnumber); - if(fUseNchWeight){ - Int_t nChargedMCPhysicalPrimary=AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(mcArray,-1.0,1.0); - fWeight *= GetNchWeight(nChargedMCPhysicalPrimary); - AliDebug(2,Form("Using Nch weights, Mult=%d Weight=%f\n",nChargedMCPhysicalPrimary,fWeight)); - } + Double_t zPrimVertex = aodVtx ->GetZ(); + Double_t zMCVertex = mcHeader->GetVtxZ(); + Int_t runnumber = aodEvent->GetRunNumber(); + fWeight=1.; + if(fUseZWeight) fWeight *= GetZWeight(zMCVertex,runnumber); + if(fUseNchWeight){ + Int_t nChargedMCPhysicalPrimary=AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(mcArray,-1.0,1.0); + fWeight *= GetNchWeight(nChargedMCPhysicalPrimary); + AliDebug(2,Form("Using Nch weights, Mult=%d Weight=%f\n",nChargedMCPhysicalPrimary,fWeight)); + } - if (TMath::Abs(zMCVertex) > fCuts->GetMaxVtxZ()){ - AliDebug(3,Form("z coordinate of MC vertex = %f, it was required to be within [-%f, +%f], skipping event", zMCVertex, fCuts->GetMaxVtxZ(), fCuts->GetMaxVtxZ())); - delete[] containerInput; - delete[] containerInputMC; - return; - } + if (TMath::Abs(zMCVertex) > fCuts->GetMaxVtxZ()){ + AliDebug(3,Form("z coordinate of MC vertex = %f, it was required to be within [-%f, +%f], skipping event", zMCVertex, fCuts->GetMaxVtxZ(), fCuts->GetMaxVtxZ())); + delete[] containerInput; + delete[] containerInputMC; + return; + } - AliESDtrackCuts** trackCuts = new AliESDtrackCuts*[cfVtxHF->GetNProngs()]; - if (fDecayChannel == 21){ - // for the D*, setting the third element of the array of the track cuts to those for the soft pion - for (Int_t iProng = 0; iProngGetNProngs()-1; iProng++){ - trackCuts[iProng]=fCuts->GetTrackCuts(); - } - trackCuts[2] = fCuts->GetTrackCutsSoftPi(); - } - else { - for (Int_t iProng = 0; iProngGetNProngs(); iProng++){ - trackCuts[iProng]=fCuts->GetTrackCuts(); - } - } + AliESDtrackCuts** trackCuts = new AliESDtrackCuts*[cfVtxHF->GetNProngs()]; + if (fDecayChannel == 21){ + // for the D*, setting the third element of the array of the track cuts to those for the soft pion + for (Int_t iProng = 0; iProngGetNProngs()-1; iProng++){ + trackCuts[iProng]=fCuts->GetTrackCuts(); + } + trackCuts[2] = fCuts->GetTrackCutsSoftPi(); + } + else if (fDecayChannel == 22) { + // for the Lc->V0+bachelor, setting the second and third elements of the array of the track cuts to those for the V0 daughters + trackCuts[0]=fCuts->GetTrackCuts(); + trackCuts[1]=fCuts->GetTrackCutsV0daughters(); + trackCuts[2]=fCuts->GetTrackCutsV0daughters(); + } + else { + for (Int_t iProng = 0; iProngGetNProngs(); iProng++){ + trackCuts[iProng]=fCuts->GetTrackCuts(); + } + } - //General settings: vertex, feed down and fill reco container with generated values. - cfVtxHF->SetRecoPrimVertex(zPrimVertex); - cfVtxHF->SetMCPrimaryVertex(zMCVertex); - cfVtxHF->SetFillFromGenerated(fFillFromGenerated); - cfVtxHF->SetNVar(fNvar); - cfVtxHF->SetFakeSelection(fFakeSelection); - cfVtxHF->SetRejectCandidateIfNotFromQuark(fRejectIfNoQuark); - cfVtxHF->SetConfiguration(fConfiguration); - - // switch-off the trigger class selection (doesn't work for MC) - fCuts->SetTriggerClass(""); - - // MC vertex, to be used, in case, for pp - if (fUseMCVertex) fCuts->SetUseMCVertex(); - - if (fCentralitySelection){ // keep only the requested centrality - if(fCuts->IsEventSelectedInCentrality(aodEvent)!=0) { - delete[] containerInput; - delete[] containerInputMC; - delete [] trackCuts; - return; - } - } else { // keep all centralities - fCuts->SetMinCentrality(0.); - fCuts->SetMaxCentrality(100.); - } + //General settings: vertex, feed down and fill reco container with generated values. + cfVtxHF->SetRecoPrimVertex(zPrimVertex); + cfVtxHF->SetMCPrimaryVertex(zMCVertex); + cfVtxHF->SetFillFromGenerated(fFillFromGenerated); + cfVtxHF->SetNVar(fNvar); + cfVtxHF->SetFakeSelection(fFakeSelection); + cfVtxHF->SetRejectCandidateIfNotFromQuark(fRejectIfNoQuark); + cfVtxHF->SetConfiguration(fConfiguration); + + // switch-off the trigger class selection (doesn't work for MC) + fCuts->SetTriggerClass(""); + + // MC vertex, to be used, in case, for pp + if (fUseMCVertex) fCuts->SetUseMCVertex(); + + if (fCentralitySelection){ // keep only the requested centrality + if(fCuts->IsEventSelectedInCentrality(aodEvent)!=0) { + delete[] containerInput; + delete[] containerInputMC; + delete [] trackCuts; + return; + } + } else { // keep all centralities + fCuts->SetMinCentrality(0.); + fCuts->SetMaxCentrality(100.); + } - Float_t centValue = fCuts->GetCentrality(aodEvent); - cfVtxHF->SetCentralityValue(centValue); + Float_t centValue = fCuts->GetCentrality(aodEvent); + cfVtxHF->SetCentralityValue(centValue); - // number of tracklets - multiplicity estimator - Double_t multiplicity = (Double_t)(aodEvent->GetTracklets()->GetNumberOfTracklets()); // casted to double because the CF is filled with doubles - cfVtxHF->SetMultiplicity(multiplicity); + // number of tracklets - multiplicity estimator + Double_t multiplicity = (Double_t)(aodEvent->GetTracklets()->GetNumberOfTracklets()); // casted to double because the CF is filled with doubles + cfVtxHF->SetMultiplicity(multiplicity); - for (Int_t iPart=0; iPartGetEntriesFast(); iPart++) { - AliAODMCParticle* mcPart = dynamic_cast(mcArray->At(iPart)); - if (!mcPart){ - AliError("Failed casting particle from MC array!, Skipping particle"); - continue; - } - // check the MC-level cuts, must be the desidered particle - if (!fCFManager->CheckParticleCuts(0, mcPart)) { - continue; // 0 stands for MC level - } - cfVtxHF->SetMCCandidateParam(iPart); + for (Int_t iPart=0; iPartGetEntriesFast(); iPart++) { + AliAODMCParticle* mcPart = dynamic_cast(mcArray->At(iPart)); + if (!mcPart){ + AliError("Failed casting particle from MC array!, Skipping particle"); + continue; + } + // check the MC-level cuts, must be the desidered particle + if (!fCFManager->CheckParticleCuts(0, mcPart)) { + AliDebug(2,"Check the MC-level cuts - not desidered particle"); + continue; // 0 stands for MC level + } + cfVtxHF->SetMCCandidateParam(iPart); - //counting c quarks - cquarks += cfVtxHF->MCcquarkCounting(mcPart); + //counting c quarks + cquarks += cfVtxHF->MCcquarkCounting(mcPart); - if (!(cfVtxHF->SetLabelArray())){ - AliDebug(2,Form("Impossible to set the label array (decaychannel = %d)",fDecayChannel)); - continue; - } + if (!(cfVtxHF->SetLabelArray())){ + AliDebug(2,Form("Impossible to set the label array (decaychannel = %d)",fDecayChannel)); + continue; + } - //check the candiate family at MC level - if (!(cfVtxHF->CheckMCPartFamily(mcPart, mcArray))) { - AliDebug(2,Form("Check on the family wrong!!! (decaychannel = %d)",fDecayChannel)); - continue; - } - else{ - AliDebug(2,Form("Check on the family OK!!! (decaychannel = %d)",fDecayChannel)); - } + //check the candiate family at MC level + if (!(cfVtxHF->CheckMCPartFamily(mcPart, mcArray))) { + AliDebug(2,Form("Check on the family wrong!!! (decaychannel = %d)",fDecayChannel)); + continue; + } + else{ + AliInfo(Form("Check on the family OK!!! (decaychannel = %d)",fDecayChannel)); + } - //Fill the MC container - Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC); - if (mcContainerFilled) { - if (fUseWeight){ - if (fFuncWeight){ // using user-defined function - AliDebug(2,"Using function"); - fWeight = fFuncWeight->Eval(containerInputMC[0]); - } - else{ // using FONLL - AliDebug(2,"Using FONLL"); - fWeight = GetWeight(containerInputMC[0]); - } - AliDebug(2,Form("pt = %f, weight = %f",containerInputMC[0], fWeight)); - } - if (!fCuts->IsInFiducialAcceptance(containerInputMC[0],containerInputMC[1])) continue; - //MC Limited Acceptance - if (TMath::Abs(containerInputMC[1]) < 0.5) { - fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc, fWeight); - AliDebug(3,"MC Lim Acc container filled\n"); - } + //Fill the MC container + Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC); + AliDebug(2,Form("mcContainerFilled = %d)",mcContainerFilled)); + if (mcContainerFilled) { + if (fUseWeight){ + if (fFuncWeight){ // using user-defined function + AliDebug(2,"Using function"); + fWeight = fFuncWeight->Eval(containerInputMC[0]); + } + else{ // using FONLL + AliDebug(2,"Using FONLL"); + fWeight = GetWeight(containerInputMC[0]); + } + AliDebug(2,Form("pt = %f, weight = %f",containerInputMC[0], fWeight)); + } + if (!fCuts->IsInFiducialAcceptance(containerInputMC[0],containerInputMC[1])) continue; + //MC Limited Acceptance + if (TMath::Abs(containerInputMC[1]) < 0.5) { + fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc, fWeight); + AliDebug(3,"MC Lim Acc container filled\n"); + } - //MC - fCFManager->GetParticleContainer()->Fill(containerInputMC, kStepGenerated, fWeight); - icountMC++; - AliDebug(3,"MC cointainer filled \n"); + //MC + fCFManager->GetParticleContainer()->Fill(containerInputMC, kStepGenerated, fWeight); + icountMC++; + AliDebug(3,"MC cointainer filled \n"); - // MC in acceptance - // check the MC-Acceptance level cuts - // since standard CF functions are not applicable, using Kine Cuts on daughters - Bool_t mcAccepStep = cfVtxHF-> MCAcceptanceStep(); - if (mcAccepStep){ - fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance, fWeight); - AliDebug(3,"MC acceptance cut passed\n"); - icountAcc++; + // MC in acceptance + // check the MC-Acceptance level cuts + // since standard CF functions are not applicable, using Kine Cuts on daughters + Bool_t mcAccepStep = cfVtxHF-> MCAcceptanceStep(); + if (mcAccepStep){ + fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance, fWeight); + AliDebug(3,"MC acceptance cut passed\n"); + icountAcc++; - //MC Vertex step - if (fCuts->IsEventSelected(aodEvent)){ - // filling the container if the vertex is ok - fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex, fWeight) ; - AliDebug(3,"Vertex cut passed and container filled\n"); - icountVertex++; + //MC Vertex step + if (fCuts->IsEventSelected(aodEvent)){ + // filling the container if the vertex is ok + fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex, fWeight) ; + AliDebug(3,"Vertex cut passed and container filled\n"); + icountVertex++; - //mc Refit requirement - Bool_t mcRefitStep = cfVtxHF->MCRefitStep(aodEvent, &trackCuts[0]); - if (mcRefitStep){ - fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit, fWeight); - AliDebug(3,"MC Refit cut passed and container filled\n"); - icountRefit++; - } - else{ - AliDebug(3,"MC Refit cut not passed\n"); - continue; - } - } - else{ - AliDebug (3, "MC vertex step not passed\n"); - continue; - } - } - else{ - AliDebug (3, "MC in acceptance step not passed\n"); - continue; - } - } - else { - AliDebug (3, "MC container not filled\n"); - } + //mc Refit requirement + Bool_t mcRefitStep = cfVtxHF->MCRefitStep(aodEvent, &trackCuts[0]); + if (mcRefitStep){ + fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit, fWeight); + AliDebug(3,"MC Refit cut passed and container filled\n"); + icountRefit++; + } + else{ + AliDebug(3,"MC Refit cut not passed\n"); + continue; + } } + else{ + AliDebug (3, "MC vertex step not passed\n"); + continue; + } + } + else{ + AliDebug (3, "MC in acceptance step not passed\n"); + continue; + } + } + else { + AliDebug (3, "MC container not filled\n"); + } + } - if (cquarks<2) AliDebug(2,Form("Event with %d c-quarks", cquarks)); - AliDebug(2,Form("Found %i MC particles that are %s!!",icountMC,fPartName.Data())); - AliDebug(2,Form("Found %i MC particles that are %s and satisfy Acc cuts!!",icountAcc,fPartName.Data())); - AliDebug(2,Form("Found %i MC particles that are %s and satisfy Vertex cuts!!",icountVertex,fPartName.Data())); - AliDebug(2,Form("Found %i MC particles that are %s and satisfy Refit cuts!!",icountRefit,fPartName.Data())); - - // Now go to rec level - fCountMC += icountMC; - fCountAcc += icountAcc; - fCountVertex+= icountVertex; - fCountRefit+= icountRefit; - - AliDebug(2,Form("Found %d vertices for decay channel %d",arrayBranch->GetEntriesFast(),fDecayChannel)); + if (cquarks<2) AliDebug(2,Form("Event with %d c-quarks", cquarks)); + AliDebug(2,Form("Found %i MC particles that are %s!!",icountMC,fPartName.Data())); + AliDebug(2,Form("Found %i MC particles that are %s and satisfy Acc cuts!!",icountAcc,fPartName.Data())); + AliDebug(2,Form("Found %i MC particles that are %s and satisfy Vertex cuts!!",icountVertex,fPartName.Data())); + AliDebug(2,Form("Found %i MC particles that are %s and satisfy Refit cuts!!",icountRefit,fPartName.Data())); + + // Now go to rec level + fCountMC += icountMC; + fCountAcc += icountAcc; + fCountVertex+= icountVertex; + fCountRefit+= icountRefit; + + AliDebug(2,Form("Found %d vertices for decay channel %d",arrayBranch->GetEntriesFast(),fDecayChannel)); - for(Int_t iCandid = 0; iCandidGetEntriesFast();iCandid++){ - AliAODRecoDecayHF* charmCandidate=0x0; - switch (fDecayChannel){ - case 2:{ - charmCandidate = (AliAODRecoDecayHF2Prong*)arrayBranch->At(iCandid); - break; - } - case 21:{ - charmCandidate = (AliAODRecoCascadeHF*)arrayBranch->At(iCandid); - break; - } - case 31: - case 32: - case 33:{ - charmCandidate = (AliAODRecoDecayHF3Prong*)arrayBranch->At(iCandid); - break; - } - case 4:{ - charmCandidate = (AliAODRecoDecayHF4Prong*)arrayBranch->At(iCandid); - break; - } - default: - break; - } + for(Int_t iCandid = 0; iCandidGetEntriesFast();iCandid++){ + AliAODRecoDecayHF* charmCandidate=0x0; + switch (fDecayChannel){ + case 2:{ + charmCandidate = (AliAODRecoDecayHF2Prong*)arrayBranch->At(iCandid); + break; + } + case 21: + case 22:{ + charmCandidate = (AliAODRecoCascadeHF*)arrayBranch->At(iCandid); + break; + } + case 31: + case 32: + case 33:{ + charmCandidate = (AliAODRecoDecayHF3Prong*)arrayBranch->At(iCandid); + break; + } + case 4:{ + charmCandidate = (AliAODRecoDecayHF4Prong*)arrayBranch->At(iCandid); + break; + } + default: + break; + } - Bool_t unsetvtx=kFALSE; - if(!charmCandidate->GetOwnPrimaryVtx()) { - charmCandidate->SetOwnPrimaryVtx(aodVtx); // needed to compute all variables - unsetvtx=kTRUE; - } + Bool_t unsetvtx=kFALSE; + if(!charmCandidate->GetOwnPrimaryVtx()) { + charmCandidate->SetOwnPrimaryVtx(aodVtx); // needed to compute all variables + unsetvtx=kTRUE; + } - Bool_t signAssociation = cfVtxHF->SetRecoCandidateParam((AliAODRecoDecayHF*)charmCandidate); - if (!signAssociation){ - charmCandidate = 0x0; - continue; - } + Bool_t signAssociation = cfVtxHF->SetRecoCandidateParam((AliAODRecoDecayHF*)charmCandidate); + if (!signAssociation){ + charmCandidate = 0x0; + continue; + } - Int_t isPartOrAntipart = cfVtxHF->CheckReflexion(fSign); - if (isPartOrAntipart == 0){ - AliDebug(2, Form("The candidate pdg code doesn't match the requirement set in the task (fSign = %d)",fSign)); - continue; - } + Int_t isPartOrAntipart = cfVtxHF->CheckReflexion(fSign); + if (isPartOrAntipart == 0){ + AliDebug(2, Form("The candidate pdg code doesn't match the requirement set in the task (fSign = %d)",fSign)); + continue; + } + AliDebug(3,Form("iCandid=%d - signAssociation=%d, isPartOrAntipart=%d",iCandid, signAssociation, isPartOrAntipart)); - Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput); - if (recoContFilled){ + Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput); + if (recoContFilled){ - // weight according to pt - if (fUseWeight){ - if (fFuncWeight){ // using user-defined function - AliDebug(2, "Using function"); - fWeight = fFuncWeight->Eval(containerInput[0]); - } - else{ // using FONLL - AliDebug(2, "Using FONLL"); - fWeight = GetWeight(containerInput[0]); - } - AliDebug(2, Form("pt = %f, weight = %f",containerInput[0], fWeight)); - } + // weight according to pt + if (fUseWeight){ + if (fFuncWeight){ // using user-defined function + AliDebug(2, "Using function"); + fWeight = fFuncWeight->Eval(containerInput[0]); + } + else{ // using FONLL + AliDebug(2, "Using FONLL"); + fWeight = GetWeight(containerInput[0]); + } + AliDebug(2, Form("pt = %f, weight = %f",containerInput[0], fWeight)); + } - if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue; + if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue; - //Reco Step - Bool_t recoStep = cfVtxHF->RecoStep(); - Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent); + //Reco Step + Bool_t recoStep = cfVtxHF->RecoStep(); + Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent); - if (recoStep && recoContFilled && vtxCheck){ - fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ; - icountReco++; - AliDebug(3,"Reco step passed and container filled\n"); + + if (recoStep && recoContFilled && vtxCheck){ + fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ; + icountReco++; + AliDebug(3,"Reco step passed and container filled\n"); - //Reco in the acceptance -- take care of UNFOLDING!!!! - Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(&trackCuts[0]); - if (recoAcceptanceStep) { - fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ; - icountRecoAcc++; - AliDebug(3,"Reco acceptance cut passed and container filled\n"); + //Reco in the acceptance -- take care of UNFOLDING!!!! + Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(&trackCuts[0]); + if (recoAcceptanceStep) { + fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ; + icountRecoAcc++; + AliDebug(3,"Reco acceptance cut passed and container filled\n"); - if(fAcceptanceUnf){ - Double_t fill[4]; //fill response matrix - Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill); - if (bUnfolding) fCorrelation->Fill(fill); - } + if(fAcceptanceUnf){ + Double_t fill[4]; //fill response matrix + Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill); + if (bUnfolding) fCorrelation->Fill(fill); + } - //Number of ITS cluster requirements - Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks); - if (recoITSnCluster){ - fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ; - icountRecoITSClusters++; - AliDebug(3,"Reco n ITS cluster cut passed and container filled\n"); + //Number of ITS cluster requirements + Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks); + if (recoITSnCluster){ + fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ; + icountRecoITSClusters++; + AliDebug(3,"Reco n ITS cluster cut passed and container filled\n"); - Bool_t iscutsusingpid = fCuts->GetIsUsePID(); - Int_t recoAnalysisCuts = -1, recoPidSelection = -1; - fCuts->SetUsePID(kFALSE); - recoAnalysisCuts = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent); - - if (fDecayChannel==33){ // Ds case, where more possibilities are considered - Bool_t keepDs=ProcessDs(recoAnalysisCuts); - if(keepDs) recoAnalysisCuts=3; - } + Bool_t iscutsusingpid = fCuts->GetIsUsePID(); + Int_t recoAnalysisCuts = -1, recoPidSelection = -1; + fCuts->SetUsePID(kFALSE); + recoAnalysisCuts = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent); + + if (fDecayChannel==33){ // Ds case, where more possibilities are considered + Bool_t keepDs=ProcessDs(recoAnalysisCuts); + if(keepDs) recoAnalysisCuts=3; + } + else if (fDecayChannel==22){ // Lc->V0+bachelor case, where more possibilities are considered + Bool_t keepLctoV0bachelor=ProcessLctoV0Bachelor(recoAnalysisCuts); + if (keepLctoV0bachelor) recoAnalysisCuts=3; + } + - fCuts->SetUsePID(iscutsusingpid); //restoring usage of the PID from the cuts object - Bool_t tempAn=(recoAnalysisCuts == 3 || recoAnalysisCuts == isPartOrAntipart); - if (fDecayChannel == 32) tempAn=(recoAnalysisCuts >0 || recoAnalysisCuts == isPartOrAntipart); + fCuts->SetUsePID(iscutsusingpid); //restoring usage of the PID from the cuts object + Bool_t tempAn=(recoAnalysisCuts == 3 || recoAnalysisCuts == isPartOrAntipart); + if (fDecayChannel == 32) tempAn=(recoAnalysisCuts >0 || recoAnalysisCuts == isPartOrAntipart); - if (tempAn){ - fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPPR, fWeight); - icountRecoPPR++; - AliDebug(3,"Reco Analysis cuts passed and container filled \n"); - //pid selection - //recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kPID); - //if((fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==isPartOrAntipart)||(fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==3)){ - recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent); - - if (fDecayChannel==33){ // Ds case, where more possibilities are considered - Bool_t keepDs=ProcessDs(recoPidSelection); - if(keepDs) recoPidSelection=3; - } - Bool_t tempPid=(recoPidSelection == 3 || recoPidSelection == isPartOrAntipart); - if (fDecayChannel == 32) tempPid=(recoPidSelection >0 || recoPidSelection == isPartOrAntipart); - - if (tempPid){ - fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPID, fWeight); - icountRecoPID++; - AliDebug(3,"Reco PID cuts passed and container filled \n"); - if(!fAcceptanceUnf){ - Double_t fill[4]; //fill response matrix - Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill); - if (bUnfolding) fCorrelation->Fill(fill); - } - } - else { - AliDebug(3, "Analysis Cuts step not passed \n"); - continue; - } - } - else { - AliDebug(3, "PID selection not passed \n"); - continue; - } - } - else{ - AliDebug(3, "Number of ITS cluster step not passed\n"); - continue; - } - } - else{ - AliDebug(3, "Reco acceptance step not passed\n"); - continue; - } - } - else { - AliDebug(3, "Reco step not passed\n"); - continue; - } + if (tempAn){ + fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPPR, fWeight); + icountRecoPPR++; + AliDebug(3,"Reco Analysis cuts passed and container filled \n"); + //pid selection + //recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kPID); + //if((fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==isPartOrAntipart)||(fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==3)){ + recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent); + + if (fDecayChannel==33){ // Ds case, where more possibilities are considered + Bool_t keepDs=ProcessDs(recoPidSelection); + if(keepDs) recoPidSelection=3; + } else if (fDecayChannel==22){ // Lc->V0+bachelor case, where more possibilities are considered + Bool_t keepLctoV0bachelor=ProcessLctoV0Bachelor(recoAnalysisCuts); + if (keepLctoV0bachelor) recoAnalysisCuts=3; + } + + Bool_t tempPid=(recoPidSelection == 3 || recoPidSelection == isPartOrAntipart); + if (fDecayChannel == 32) tempPid=(recoPidSelection >0 || recoPidSelection == isPartOrAntipart); + + if (tempPid){ + fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPID, fWeight); + icountRecoPID++; + AliDebug(3,"Reco PID cuts passed and container filled \n"); + if(!fAcceptanceUnf){ + Double_t fill[4]; //fill response matrix + Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill); + if (bUnfolding) fCorrelation->Fill(fill); } + } + else { + AliDebug(3, "Analysis Cuts step not passed \n"); + continue; + } + } + else { + AliDebug(3, "PID selection not passed \n"); + continue; + } + } + else{ + AliDebug(3, "Number of ITS cluster step not passed\n"); + continue; + } + } + else{ + AliDebug(3, "Reco acceptance step not passed\n"); + continue; + } + } + else { + AliDebug(3, "Reco step not passed\n"); + continue; + } + } - if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx(); - } // end loop on candidate + if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx(); + } // end loop on candidate - fCountReco+= icountReco; - fCountRecoAcc+= icountRecoAcc; - fCountRecoITSClusters+= icountRecoITSClusters; - fCountRecoPPR+= icountRecoPPR; - fCountRecoPID+= icountRecoPID; + fCountReco+= icountReco; + fCountRecoAcc+= icountRecoAcc; + fCountRecoITSClusters+= icountRecoITSClusters; + fCountRecoPPR+= icountRecoPPR; + fCountRecoPID+= icountRecoPID; - fHistEventsProcessed->Fill(0); - - delete[] containerInput; - delete[] containerInputMC; - delete cfVtxHF; - if (trackCuts){ - // for (Int_t i=0; iGetNProngs(); i++){ - // delete [] trackCuts[i]; - // } - delete [] trackCuts; - } + fHistEventsProcessed->Fill(0); + + delete[] containerInput; + delete[] containerInputMC; + delete cfVtxHF; + if (trackCuts){ + // for (Int_t i=0; iGetNProngs(); i++){ + // delete [] trackCuts[i]; + // } + delete [] trackCuts; + } } @@ -907,191 +963,243 @@ void AliCFTaskVertexingHF::UserExec(Option_t *) //___________________________________________________________________________ void AliCFTaskVertexingHF::Terminate(Option_t*) { - // The Terminate() function is the last function to be called during - // a query. It always runs on the client, it can be used to present - // the results graphically or save the results to file. + // The Terminate() function is the last function to be called during + // a query. It always runs on the client, it can be used to present + // the results graphically or save the results to file. - AliAnalysisTaskSE::Terminate(); - - AliInfo(Form("Found %i MC particles that are %s in MC, in %d events",fCountMC,fPartName.Data(),fEvents)); - AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, in %d events",fCountAcc,fPartName.Data(),fEvents)); - AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, and satisfy Vertex requirement in %d events",fCountVertex,fPartName.Data(),fEvents)); - AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, and satisfy ITS+TPC refit requirementin %d events",fCountRefit,fPartName.Data(),fEvents)); - AliInfo(Form("Found %i reco %s that are decaying in %s, in %d events",fCountReco,fPartName.Data(),fDauNames.Data(),fEvents)); - AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and are in the requested acceptance, in %d events",fCountRecoAcc,fPartName.Data(),fDauNames.Data(),fEvents)); - AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and have at least 5 clusters in ITS, in %d events",fCountRecoITSClusters,fPartName.Data(),fDauNames.Data(),fEvents)); - AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and satisfy PPR cuts, in %d events",fCountRecoPPR,fPartName.Data(),fDauNames.Data(),fEvents)); - AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and satisfy PPR+PID cuts, in %d events",fCountRecoPID,fPartName.Data(),fDauNames.Data(),fEvents)); + AliAnalysisTaskSE::Terminate(); + + AliInfo(Form("Found %i MC particles that are %s in MC, in %d events",fCountMC,fPartName.Data(),fEvents)); + AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, in %d events",fCountAcc,fPartName.Data(),fEvents)); + AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, and satisfy Vertex requirement in %d events",fCountVertex,fPartName.Data(),fEvents)); + AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, and satisfy ITS+TPC refit requirementin %d events",fCountRefit,fPartName.Data(),fEvents)); + AliInfo(Form("Found %i reco %s that are decaying in %s, in %d events",fCountReco,fPartName.Data(),fDauNames.Data(),fEvents)); + AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and are in the requested acceptance, in %d events",fCountRecoAcc,fPartName.Data(),fDauNames.Data(),fEvents)); + AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and have at least 5 clusters in ITS, in %d events",fCountRecoITSClusters,fPartName.Data(),fDauNames.Data(),fEvents)); + AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and satisfy PPR cuts, in %d events",fCountRecoPPR,fPartName.Data(),fDauNames.Data(),fEvents)); + AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and satisfy PPR+PID cuts, in %d events",fCountRecoPID,fPartName.Data(),fDauNames.Data(),fEvents)); - // draw some example plots.... - AliCFContainer *cont= dynamic_cast (GetOutputData(2)); - if(!cont) { - printf("CONTAINER NOT FOUND\n"); - return; - } - // projecting the containers to obtain histograms - // first argument = variable, second argument = step - - TH1D** h = new TH1D*[3]; - if (fConfiguration == kSnail){ - //h = new TH1D[3][12]; - for (Int_t ih = 0; ih<3; ih++){ - h[ih] = new TH1D[12]; - } - for(Int_t iC=1;iC<12; iC++){ - // MC-level - h[0][iC] = *(cont->ShowProjection(iC,0)); - // MC-Acceptance level - h[1][iC] = *(cont->ShowProjection(iC,1)); - // Reco-level - h[2][iC] = *(cont->ShowProjection(iC,4)); - } - } - else{ - //h = new TH1D[3][12]; - for (Int_t ih = 0; ih<3; ih++){ - h[ih] = new TH1D[8]; - } - for(Int_t iC=0;iC<8; iC++){ - // MC-level - h[0][iC] = *(cont->ShowProjection(iC,0)); - // MC-Acceptance level - h[1][iC] = *(cont->ShowProjection(iC,1)); - // Reco-level - h[2][iC] = *(cont->ShowProjection(iC,4)); - } - } - TString* titles; - Int_t nvarToPlot = 0; - if (fConfiguration == kSnail){ - nvarToPlot = 12; - titles = new TString[nvarToPlot]; - if(fDecayChannel==31){ - titles[0]="pT_Dplus (GeV/c)"; - titles[1]="rapidity"; - titles[2]="phi (rad)"; - titles[3]="cT (#mum)"; - titles[4]="cosPointingAngle"; - titles[5]="pT_1 (GeV/c)"; - titles[6]="pT_2 (GeV/c)"; - titles[7]="pT_3 (GeV/c)"; - titles[8]="d0_1 (#mum)"; - titles[9]="d0_2 (#mum)"; - titles[10]="d0_3 (#mum)"; - titles[11]="zVertex (cm)"; - }else{ - titles[0]="pT_D0 (GeV/c)"; - titles[1]="rapidity"; - titles[2]="cosThetaStar"; - titles[3]="pT_pi (GeV/c)"; - titles[4]="pT_K (Gev/c)"; - titles[5]="cT (#mum)"; - titles[6]="dca (#mum)"; - titles[7]="d0_pi (#mum)"; - titles[8]="d0_K (#mum)"; - titles[9]="d0xd0 (#mum^2)"; - titles[10]="cosPointingAngle"; - titles[11]="phi (rad)"; - - } - } - else{ - nvarToPlot = 8; - titles = new TString[nvarToPlot]; - titles[0]="pT_candidate (GeV/c)"; - titles[1]="rapidity"; - titles[2]="cT (#mum)"; - titles[3]="phi"; - titles[4]="z_{vtx}"; - titles[5]="centrality"; - titles[6]="fake"; - titles[7]="multiplicity"; - } + // draw some example plots.... + AliCFContainer *cont= dynamic_cast (GetOutputData(2)); + if(!cont) { + printf("CONTAINER NOT FOUND\n"); + return; + } + // projecting the containers to obtain histograms + // first argument = variable, second argument = step - Int_t markers[12]={20,24,21,25,27,28, - 20,24,21,25,27,28}; - Int_t colors[3]={2,8,4}; - for(Int_t iC=0;iCSetTitle(titles[iC].Data()); - Double_t maxh=h[iStep][iC].GetMaximum(); - h[iStep][iC].GetYaxis()->SetRangeUser(0,maxh*1.2); - h[iStep][iC].SetMarkerStyle(markers[iC]); - h[iStep][iC].SetMarkerColor(colors[iStep]); - } - } + TH1D** h = new TH1D*[3]; + Int_t nvarToPlot = 0; + if (fConfiguration == kSnail){ + //h = new TH1D[3][12]; + for (Int_t ih = 0; ih<3; ih++){ + if(fDecayChannel==22){ + nvarToPlot = 16; + } else { + nvarToPlot = 12; + } + h[ih] = new TH1D[nvarToPlot]; + } + for(Int_t iC=1;iCShowProjection(iC,0)); + // MC-Acceptance level + h[1][iC] = *(cont->ShowProjection(iC,1)); + // Reco-level + h[2][iC] = *(cont->ShowProjection(iC,4)); + } + } + else { + //h = new TH1D[3][12]; + nvarToPlot = 8; + for (Int_t ih = 0; ih<3; ih++){ + h[ih] = new TH1D[nvarToPlot]; + } + for(Int_t iC=0;iCShowProjection(iC,0)); + // MC-Acceptance level + h[1][iC] = *(cont->ShowProjection(iC,1)); + // Reco-level + h[2][iC] = *(cont->ShowProjection(iC,4)); + } + } + TString* titles; + //Int_t nvarToPlot = 0; + if (fConfiguration == kSnail){ + if(fDecayChannel==31){ + //nvarToPlot = 12; + titles = new TString[nvarToPlot]; + titles[0]="pT_Dplus (GeV/c)"; + titles[1]="rapidity"; + titles[2]="phi (rad)"; + titles[3]="cT (#mum)"; + titles[4]="cosPointingAngle"; + titles[5]="pT_1 (GeV/c)"; + titles[6]="pT_2 (GeV/c)"; + titles[7]="pT_3 (GeV/c)"; + titles[8]="d0_1 (#mum)"; + titles[9]="d0_2 (#mum)"; + titles[10]="d0_3 (#mum)"; + titles[11]="zVertex (cm)"; + } else if (fDecayChannel==22) { + //nvarToPlot = 16; + titles = new TString[nvarToPlot]; + titles[0]="pT_Lc (GeV/c)"; + titles[1]="rapidity"; + titles[2]="phi (rad)"; + titles[3]="cosPAV0"; + titles[4]="onTheFlyStatusV0"; + titles[5]="centrality"; + titles[6]="fake"; + titles[7]="multiplicity"; + titles[8]="pT_bachelor (GeV/c)"; + titles[9]="pT_V0pos (GeV/c)"; + titles[10]="pT_V0neg (GeV/c)"; + titles[11]="invMassV0 (GeV/c2)"; + titles[12]="dcaV0 (nSigma)"; + titles[13]="c#tauV0 (#mum)"; + titles[14]="c#tau (#mum)"; + titles[15]="cosPA"; + } else { + //nvarToPlot = 12; + titles = new TString[nvarToPlot]; + titles[0]="pT_D0 (GeV/c)"; + titles[1]="rapidity"; + titles[2]="cosThetaStar"; + titles[3]="pT_pi (GeV/c)"; + titles[4]="pT_K (Gev/c)"; + titles[5]="cT (#mum)"; + titles[6]="dca (#mum)"; + titles[7]="d0_pi (#mum)"; + titles[8]="d0_K (#mum)"; + titles[9]="d0xd0 (#mum^2)"; + titles[10]="cosPointingAngle"; + titles[11]="phi (rad)"; + } + } + else { + //nvarToPlot = 8; + titles = new TString[nvarToPlot]; + if (fDecayChannel==22) { + titles[0]="pT_candidate (GeV/c)"; + titles[1]="rapidity"; + titles[2]="phi (rad)"; + titles[3]="cosPAV0"; + titles[4]="onTheFlyStatusV0"; + titles[5]="centrality"; + titles[6]="fake"; + titles[7]="multiplicity"; + } else { + titles[0]="pT_candidate (GeV/c)"; + titles[1]="rapidity"; + titles[2]="cT (#mum)"; + titles[3]="phi"; + titles[4]="z_{vtx}"; + titles[5]="centrality"; + titles[6]="fake"; + titles[7]="multiplicity"; + } + } + + Int_t markers[16]={20,24,21,25,27,28, + 20,24,21,25,27,28, + 20,24,21,25}; + Int_t colors[3]={2,8,4}; + for(Int_t iC=0;iCSetTitle(titles[iC].Data()); + Double_t maxh=h[iStep][iC].GetMaximum(); + h[iStep][iC].GetYaxis()->SetRangeUser(0,maxh*1.2); + h[iStep][iC].SetMarkerStyle(markers[iC]); + h[iStep][iC].SetMarkerColor(colors[iStep]); + } + } - gStyle->SetCanvasColor(0); - gStyle->SetFrameFillColor(0); - gStyle->SetTitleFillColor(0); - gStyle->SetStatColor(0); + gStyle->SetCanvasColor(0); + gStyle->SetFrameFillColor(0); + gStyle->SetTitleFillColor(0); + gStyle->SetStatColor(0); - // drawing in 2 separate canvas for a matter of clearity - TCanvas * c1 =new TCanvas(Form("c1New_%d",fDecayChannel),"Vars 0, 1, 2, 3",1100,1200); - c1->Divide(3,4); - Int_t iPad=1; - for(Int_t iVar=0; iVar<4; iVar++){ - c1->cd(iPad++); - h[0][iVar].DrawCopy("p"); - c1->cd(iPad++); - h[1][iVar].DrawCopy("p"); - c1->cd(iPad++); - h[2][iVar].DrawCopy("p"); - } + // drawing in 2 separate canvas for a matter of clearity + TCanvas * c1 =new TCanvas(Form("c1New_%d",fDecayChannel),"Vars 0, 1, 2, 3",1100,1200); + c1->Divide(3,4); + Int_t iPad=1; + for(Int_t iVar=0; iVar<4; iVar++){ + c1->cd(iPad++); + h[0][iVar].DrawCopy("p"); + c1->cd(iPad++); + h[1][iVar].DrawCopy("p"); + c1->cd(iPad++); + h[2][iVar].DrawCopy("p"); + } - TCanvas * c2 =new TCanvas(Form("c2New_%d",fDecayChannel),"Vars 4, 5, 6, 7",1100,1200); - c2->Divide(3,4); - iPad=1; - for(Int_t iVar=4; iVar<8; iVar++){ - c2->cd(iPad++); - h[0][iVar].DrawCopy("p"); - c2->cd(iPad++); - h[1][iVar].DrawCopy("p"); - c2->cd(iPad++); - h[2][iVar].DrawCopy("p"); - } + TCanvas * c2 =new TCanvas(Form("c2New_%d",fDecayChannel),"Vars 4, 5, 6, 7",1100,1200); + c2->Divide(3,4); + iPad=1; + for(Int_t iVar=4; iVar<8; iVar++){ + c2->cd(iPad++); + h[0][iVar].DrawCopy("p"); + c2->cd(iPad++); + h[1][iVar].DrawCopy("p"); + c2->cd(iPad++); + h[2][iVar].DrawCopy("p"); + } - if (fConfiguration == kSnail){ - TCanvas * c3 =new TCanvas(Form("c3New_%d",fDecayChannel),"Vars 8, 9, 10, 11",1100,1200); - c3->Divide(3,4); - iPad=1; - for(Int_t iVar=8; iVar<12; iVar++){ - c3->cd(iPad++); - h[0][iVar].DrawCopy("p"); - c3->cd(iPad++); - h[1][iVar].DrawCopy("p"); - c3->cd(iPad++); - h[2][iVar].DrawCopy("p"); - } - } + if (fConfiguration == kSnail){ + TCanvas * c3 =new TCanvas(Form("c3New_%d",fDecayChannel),"Vars 8, 9, 10, 11",1100,1200); + c3->Divide(3,4); + iPad=1; + for(Int_t iVar=8; iVar<12; iVar++){ + c3->cd(iPad++); + h[0][iVar].DrawCopy("p"); + c3->cd(iPad++); + h[1][iVar].DrawCopy("p"); + c3->cd(iPad++); + h[2][iVar].DrawCopy("p"); + } + if (fDecayChannel==22) { + TCanvas * c4 =new TCanvas(Form("c4New_%d",fDecayChannel),"Vars 12, 13, 14, 15",1100,1200); + c4->Divide(3,4); + iPad=1; + for(Int_t iVar=12; iVar<16; iVar++){ + c4->cd(iPad++); + h[0][iVar].DrawCopy("p"); + c4->cd(iPad++); + h[1][iVar].DrawCopy("p"); + c4->cd(iPad++); + h[2][iVar].DrawCopy("p"); + } + } + } - /* - THnSparseD* hcorr = dynamic_cast (GetOutputData(3)); + /* + THnSparseD* hcorr = dynamic_cast (GetOutputData(3)); - TH2D* corr1 =hcorr->Projection(0,2); - TH2D* corr2 = hcorr->Projection(1,3); + TH2D* corr1 = hcorr->Projection(0,2); + TH2D* corr2 = hcorr->Projection(1,3); - TCanvas * c7 =new TCanvas(Form("c7New_%d",fDecayChannel),"",800,400); - c7->Divide(2,1); - c7->cd(1); - corr1->DrawCopy("text"); - c7->cd(2); - corr2->DrawCopy("text"); - */ - TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE"); + TCanvas * c7 =new TCanvas(Form("c7New_%d",fDecayChannel),"",800,400); + c7->Divide(2,1); + c7->cd(1); + corr1->DrawCopy("text"); + c7->cd(2); + corr2->DrawCopy("text"); + */ + TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE"); - // corr1->Write(); - // corr2->Write(); + // corr1->Write(); + // corr2->Write(); - for(Int_t iC=0;iCClose(); - for (Int_t ih = 0; ih<3; ih++) delete [] h[ih]; - delete [] h; + for(Int_t iC=0;iCClose(); + for (Int_t ih = 0; ih<3; ih++) delete [] h[ih]; + delete [] h; } @@ -1099,19 +1207,33 @@ void AliCFTaskVertexingHF::Terminate(Option_t*) //___________________________________________________________________________ void AliCFTaskVertexingHF::UserCreateOutputObjects() { - //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED - //TO BE SET BEFORE THE EXECUTION OF THE TASK - // - Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName()); + //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED + //TO BE SET BEFORE THE EXECUTION OF THE TASK + // + Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName()); - //slot #1 - OpenFile(1); - const char* nameoutput=GetOutputSlot(1)->GetContainer()->GetName(); - fHistEventsProcessed = new TH1I(nameoutput,"",1,0,1) ; + AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager(); + AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler()); + AliPIDResponse *localPIDResponse = (AliPIDResponse*)inputHandler->GetPIDResponse(); + + if (fCuts->GetIsUsePID() && fDecayChannel==22) { + + fCuts->GetPidHF()->SetPidResponse(localPIDResponse); + (dynamic_cast(fCuts))->GetPidV0pos()->SetPidResponse(localPIDResponse); + (dynamic_cast(fCuts))->GetPidV0neg()->SetPidResponse(localPIDResponse); + fCuts->GetPidHF()->SetOldPid(kFALSE); + (dynamic_cast(fCuts))->GetPidV0pos()->SetOldPid(kFALSE); + (dynamic_cast(fCuts))->GetPidV0neg()->SetOldPid(kFALSE); + } - PostData(1,fHistEventsProcessed) ; - PostData(2,fCFManager->GetParticleContainer()) ; - PostData(3,fCorrelation) ; + //slot #1 + OpenFile(1); + const char* nameoutput=GetOutputSlot(1)->GetContainer()->GetName(); + fHistEventsProcessed = new TH1I(nameoutput,"",1,0,1) ; + + PostData(1,fHistEventsProcessed) ; + PostData(2,fCFManager->GetParticleContainer()) ; + PostData(3,fCorrelation) ; } @@ -1119,43 +1241,43 @@ void AliCFTaskVertexingHF::UserCreateOutputObjects() //_________________________________________________________________________ Double_t AliCFTaskVertexingHF::GetWeight(Float_t pt) { - // - // calculating the weight to fill the container - // + // + // calculating the weight to fill the container + // - // FNOLL central: - // p0 = 1.63297e-01 --> 0.322643 - // p1 = 2.96275e+00 - // p2 = 2.30301e+00 - // p3 = 2.50000e+00 - - // PYTHIA - // p0 = 1.85906e-01 --> 0.36609 - // p1 = 1.94635e+00 - // p2 = 1.40463e+00 - // p3 = 2.50000e+00 - - Double_t func1[4] = {0.322643,2.96275,2.30301,2.5}; - Double_t func2[4] = {0.36609,1.94635,1.40463,2.5}; - - Double_t dndpt_func1 = dNdptFit(pt,func1); - if(fUseFlatPtWeight) dndpt_func1 = 1./30.; - Double_t dndpt_func2 = dNdptFit(pt,func2); - AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2)); - return dndpt_func1/dndpt_func2; + // FNOLL central: + // p0 = 1.63297e-01 --> 0.322643 + // p1 = 2.96275e+00 + // p2 = 2.30301e+00 + // p3 = 2.50000e+00 + + // PYTHIA + // p0 = 1.85906e-01 --> 0.36609 + // p1 = 1.94635e+00 + // p2 = 1.40463e+00 + // p3 = 2.50000e+00 + + Double_t func1[4] = {0.322643,2.96275,2.30301,2.5}; + Double_t func2[4] = {0.36609,1.94635,1.40463,2.5}; + + Double_t dndpt_func1 = dNdptFit(pt,func1); + if(fUseFlatPtWeight) dndpt_func1 = 1./30.; + Double_t dndpt_func2 = dNdptFit(pt,func2); + AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2)); + return dndpt_func1/dndpt_func2; } //__________________________________________________________________________________________________ Double_t AliCFTaskVertexingHF::dNdptFit(Float_t pt, Double_t* par) { - // - // calculating dNdpt - // + // + // calculating dNdpt + // - Double_t denom = TMath::Power((pt/par[1]), par[3] ); - Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]); + Double_t denom = TMath::Power((pt/par[1]), par[3] ); + Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]); - return dNdpt; + return dNdpt; } //__________________________________________________________________________________________________ @@ -1223,6 +1345,7 @@ void AliCFTaskVertexingHF::CreateMeasuredNchHisto(){ } } +//__________________________________________________________________________________________________ Bool_t AliCFTaskVertexingHF::ProcessDs(Int_t recoAnalysisCuts) const{ // processes output of Ds is selected Bool_t keep=kFALSE; @@ -1245,3 +1368,28 @@ Bool_t AliCFTaskVertexingHF::ProcessDs(Int_t recoAnalysisCuts) const{ } return keep; } +//__________________________________________________________________________________________________ +Bool_t AliCFTaskVertexingHF::ProcessLctoV0Bachelor(Int_t recoAnalysisCuts) const{ + // processes output of Lc->V0+bnachelor is selected + Bool_t keep=kFALSE; + if(recoAnalysisCuts > 0){ + + Int_t isK0Sp = recoAnalysisCuts&1; + Int_t isLambdaBarpi = recoAnalysisCuts&2; + Int_t isLambdapi = recoAnalysisCuts&4; + + if(fLctoV0bachelorOption==1){ + if(isK0Sp) keep=kTRUE; + } + else if(fLctoV0bachelorOption==2){ + if(isLambdaBarpi) keep=kTRUE; + } + else if(fLctoV0bachelorOption==4){ + if(isLambdapi) keep=kTRUE; + } + else if(fLctoV0bachelorOption==7) { + if (isK0Sp || isLambdaBarpi || isLambdapi) keep=kTRUE; + } + } + return keep; +} diff --git a/PWGHF/vertexingHF/AliCFTaskVertexingHF.h b/PWGHF/vertexingHF/AliCFTaskVertexingHF.h index f7c1036c0d2..b5ef3d3ed82 100644 --- a/PWGHF/vertexingHF/AliCFTaskVertexingHF.h +++ b/PWGHF/vertexingHF/AliCFTaskVertexingHF.h @@ -28,6 +28,7 @@ #include "AliAnalysisTaskSE.h" #include "AliCFVertexingHF2Prong.h" #include "AliCFVertexingHF3Prong.h" +#include "AliCFVertexingHFLctoV0bachelor.h" #include "AliCFVertexingHF.h" #include @@ -160,7 +161,18 @@ public: void SetResonantDecay(UInt_t resonantDecay) {fResonantDecay = resonantDecay;} UInt_t GetResonantDecay() const {return fResonantDecay;} - + + void SetKeepLctoK0Sp() {fLctoV0bachelorOption=1;} + void SetKeepLctoLambdaBarpi() {fLctoV0bachelorOption=2;} + void SetKeepLctoLambdapi() {fLctoV0bachelorOption=4;} + void SetKeepLctoV0bachelor() {fLctoV0bachelorOption=7;} + + void SetCountAllLctoBachelor(){fGenLctoV0bachelorOption=AliCFVertexingHFLctoV0bachelor::kCountAllLctoV0;} + void SetCountLctoK0Sp(){fGenLctoV0bachelorOption=AliCFVertexingHFLctoV0bachelor::kCountK0Sp;} + void SetCountLambdaBarpi(){fGenLctoV0bachelorOption=AliCFVertexingHFLctoV0bachelor::kCountLambdapi;} + + Bool_t ProcessLctoV0Bachelor(Int_t returnCodeDs) const; + protected: AliCFManager *fCFManager; // pointer to the CF manager TH1I *fHistEventsProcessed; //! simple histo for monitoring the number of events processed @@ -200,8 +212,10 @@ protected: TH1F* fHistoMeasNch; // histogram with measured Nch distribution (pp 7 TeV) TH1F* fHistoMCNch; // histogram with Nch distribution from MC production UInt_t fResonantDecay; // resonant deacy channel to be used if the CF should be run on resonant channels only + Int_t fLctoV0bachelorOption; // Lc->V0+bachelor decay option (selection level) + Int_t fGenLctoV0bachelorOption; // Lc->V0+bachelor decay option (generation level) - ClassDef(AliCFTaskVertexingHF,12); // class for HF corrections as a function of many variables + ClassDef(AliCFTaskVertexingHF,13); // class for HF corrections as a function of many variables }; #endif diff --git a/PWGHF/vertexingHF/AliCFVertexingHF.cxx b/PWGHF/vertexingHF/AliCFVertexingHF.cxx index c30210ac73e..cf7eea08711 100644 --- a/PWGHF/vertexingHF/AliCFVertexingHF.cxx +++ b/PWGHF/vertexingHF/AliCFVertexingHF.cxx @@ -461,6 +461,7 @@ Bool_t AliCFVertexingHF::FillRecoContainer(Double_t *containerInput) delete [] vectorReco; vectorValues = 0x0; vectorReco = 0x0; + return recoContainerFilled; } @@ -813,6 +814,39 @@ Bool_t AliCFVertexingHF::SetLabelArray() return bLabelArray; } } + // added K0S case - Start + else if (pdgCode==311) { + if (part->GetNDaughters()!=1) { + delete [] fLabelArray; + fLabelArray = 0x0; + return bLabelArray; + } + Int_t labelK0Dau = part->GetDaughter(0); + AliAODMCParticle* partK0S = dynamic_cast(fmcArray->At(labelK0Dau)); + Int_t nDauRes=partK0S->GetNDaughters(); + if(nDauRes!=2 || partK0S->GetPdgCode()!=310) { + delete [] fLabelArray; + fLabelArray = 0x0; + return bLabelArray; + } + Int_t labelFirstDauRes = partK0S->GetDaughter(0); + AliDebug(2,Form(" Found K0S (%d)",labelK0Dau)); + for(Int_t iDauRes=0; iDauRes(fmcArray->At(iLabelDauRes)); + if (dauRes){ + fLabelArray[foundDaughters] = dauRes->GetLabel(); + foundDaughters++; + } + else { + AliError("Error while casting resonant daughter! returning a NULL array"); + delete [] fLabelArray; + fLabelArray = 0x0; + return bLabelArray; + } + } + } + // added K0S case - End else{ Int_t nDauRes=part->GetNDaughters(); if(nDauRes!=2) { diff --git a/PWGHF/vertexingHF/AliCFVertexingHF.h b/PWGHF/vertexingHF/AliCFVertexingHF.h index 1f8b5dde4eb..0490d8c3f08 100644 --- a/PWGHF/vertexingHF/AliCFVertexingHF.h +++ b/PWGHF/vertexingHF/AliCFVertexingHF.h @@ -43,7 +43,7 @@ class AliESDtrackCuts; class AliCFVertexingHF : public TObject { public: - enum DecayChannel{kD0toKpi = 2, kDStartoKpipi = 21, kDplustoKpipi = 31, kLctopKpi = 32, kDstoKKpi = 33, kD0toKpipipi = 4}; + enum DecayChannel{kD0toKpi = 2, kDStartoKpipi = 21, kLctoV0bachelor = 22, kDplustoKpipi = 31, kLctopKpi = 32, kDstoKKpi = 33, kD0toKpipipi = 4}; AliCFVertexingHF() ; AliCFVertexingHF(TClonesArray *mcArray, UShort_t originDselection); diff --git a/PWGHF/vertexingHF/AliCFVertexingHFLctoV0bachelor.cxx b/PWGHF/vertexingHF/AliCFVertexingHFLctoV0bachelor.cxx new file mode 100644 index 00000000000..1bbef53c65c --- /dev/null +++ b/PWGHF/vertexingHF/AliCFVertexingHFLctoV0bachelor.cxx @@ -0,0 +1,726 @@ +/************************************************************************** + * Copyright(c) 1998-2011, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id$ */ + +//----------------------------------------------------------------------- +// Class for HF corrections as a function of many variables and steps +// For Lc->V0+bachelor +// +// - Based on AliCFVertexingHFCascade - +// +// Contact : A.De Caro - decaro@sa.infn.it +// Centro 'E.Fermi' - Rome (Italy) +// INFN and University of Salerno (Italy) +// +//----------------------------------------------------------------------- + +#include "AliAODRecoDecayHF2Prong.h" +#include "AliAODMCParticle.h" +#include "AliAODEvent.h" +#include "TClonesArray.h" +#include "AliCFVertexingHF.h" +#include "AliESDtrack.h" +#include "TDatabasePDG.h" +#include "AliAODRecoCascadeHF.h" +#include "AliAODv0.h" +#include "AliCFVertexingHFLctoV0bachelor.h" +#include "AliCFContainer.h" +#include "AliCFTaskVertexingHF.h" + +#include + +using std::cout; +using std::endl; + +ClassImp(AliCFVertexingHFLctoV0bachelor) + +//_________________________________________ +AliCFVertexingHFLctoV0bachelor::AliCFVertexingHFLctoV0bachelor(): +fGenLcOption(0) +{ + // standard constructor +} + +//_____________________________________ +AliCFVertexingHFLctoV0bachelor::AliCFVertexingHFLctoV0bachelor(TClonesArray *mcArray, UShort_t originDselection, Int_t lcDecay): +AliCFVertexingHF(mcArray, originDselection), + fGenLcOption(lcDecay) +{ + // standard constructor + + SetNProngs(3); + fPtAccCut=new Float_t[fProngs]; + fEtaAccCut=new Float_t[fProngs]; + for(Int_t iP=0; iPGetPrimaryVtx()) AliDebug(3,"fReco Candidate has a pointer to PrimVtx\n"); + + AliAODRecoCascadeHF* lcV0bachelor = (AliAODRecoCascadeHF*)fRecoCandidate; + + if ( !(lcV0bachelor->Getv0()) ) return bSignAssoc; + + Int_t nentries = fmcArray->GetEntriesFast(); + AliDebug(3,Form("nentries = %d\n", nentries)); + + Int_t pdgCand = 4122; + Int_t mcLabel = -1; + Int_t mcLabelK0S = -1; + Int_t mcLabelLambda = -1; + + // Lc->K0S+p and cc + Int_t pdgDgLctoV0bachelor[2]={310,2212}; + Int_t pdgDgV0toDaughters[2]={211,211}; + mcLabelK0S = lcV0bachelor->MatchToMC(pdgCand,pdgDgLctoV0bachelor[0],pdgDgLctoV0bachelor,pdgDgV0toDaughters,fmcArray,kTRUE); + // Lc->Lambda+pi and cc + pdgDgLctoV0bachelor[0]=3122, pdgDgLctoV0bachelor[1]=211; + pdgDgV0toDaughters[0]=2212, pdgDgV0toDaughters[1]=211; + mcLabelLambda = lcV0bachelor->MatchToMC(pdgCand,pdgDgLctoV0bachelor[0],pdgDgLctoV0bachelor,pdgDgV0toDaughters,fmcArray,kTRUE); + + if (mcLabelK0S!=-1 && mcLabelLambda!=-1) + AliInfo("Strange: current Lc->V0+bachelor candidate has two MC different labels!"); + + if (fGenLcOption==kCountAllLctoV0) { + if (mcLabelK0S!=-1) mcLabel=mcLabelK0S; + if (mcLabelLambda!=-1) mcLabel=mcLabelLambda; + } + else if (fGenLcOption==kCountK0Sp) { + if (mcLabelK0S!=-1) mcLabel=mcLabelK0S; + if (mcLabelLambda!=-1) { + mcLabel=-1; + fFake = 0; // fake candidate + if (fFakeSelection==1) return bSignAssoc; + } + } + else if (fGenLcOption==kCountLambdapi) { + if (mcLabelLambda!=-1) mcLabel=mcLabelLambda; + if (mcLabelK0S!=-1) { + mcLabel=-1; + fFake = 0; // fake candidate + if (fFakeSelection==1) return bSignAssoc; + } + } + + if (mcLabel==-1) return bSignAssoc; + + if (fRecoCandidate->NumberOfFakeDaughters()>0){ + fFake = 0; // fake candidate + if (fFakeSelection==1) return bSignAssoc; + } + if (fRecoCandidate->NumberOfFakeDaughters()==0){ + fFake = 2; // non-fake candidate + if (fFakeSelection==2) return bSignAssoc; + } + + SetMCLabel(mcLabel); + fmcPartCandidate = dynamic_cast(fmcArray->At(fmcLabel)); + + if (!fmcPartCandidate){ + AliDebug(3,"No part candidate"); + return bSignAssoc; + } + + bSignAssoc = kTRUE; + return bSignAssoc; +} + +//______________________________________________ +Bool_t AliCFVertexingHFLctoV0bachelor::GetGeneratedValuesFromMCParticle(Double_t* vectorMC) +{ + // + // collecting all the necessary info (pt, y, invMassV0, cosPAwrtPVxV0, onTheFlyStatusV0) from MC particle + // (additional infos: pTbachelor, pTV0pos, pTV0neg, phi, dcaV0, cTV0, cT, cosPA) + // + + Bool_t bGenValues = kFALSE; + + if (fmcPartCandidate->GetNDaughters()!=2) return bGenValues; + + Int_t daughter0lc = fmcPartCandidate->GetDaughter(0); + Int_t daughter1lc = fmcPartCandidate->GetDaughter(1); + + //the V0 + AliAODMCParticle* mcPartDaughterV0=0; + if (fGenLcOption==kCountLambdapi) { + mcPartDaughterV0 = dynamic_cast(fmcArray->At(daughter0lc)); // strange baryon (0) + AliInfo(Form(" Case Lc->Lambda+pi: (V0) daughter0_Lc=%d (%d)",daughter0lc,mcPartDaughterV0->GetPdgCode())); + } + else if (fGenLcOption==kCountK0Sp) { + AliAODMCParticle* mcPartDaughterK0 = dynamic_cast(fmcArray->At(daughter1lc)); // strange meson (1) + if (!mcPartDaughterK0) return bGenValues; + if (TMath::Abs(mcPartDaughterK0->GetPdgCode())!=311) return bGenValues; + Int_t daughterK0 = mcPartDaughterK0->GetDaughter(0); + mcPartDaughterV0 = dynamic_cast(fmcArray->At(daughterK0)); + if (TMath::Abs(mcPartDaughterV0->GetPdgCode())!=310) return bGenValues; + AliInfo(Form(" Case Lc->K0S+p: (V0) daughter1_Lc=%d (%d to %d)",daughter1lc,mcPartDaughterK0->GetPdgCode(), + mcPartDaughterV0->GetPdgCode())); + } + + //the bachelor + AliAODMCParticle* mcPartDaughterBachelor=0; + if (fGenLcOption==kCountLambdapi) { + mcPartDaughterBachelor = dynamic_cast(fmcArray->At(daughter1lc)); // meson (1) + AliInfo(Form(" Case Lc->Lambda+pi: (bachelor) daughter1_Lc=%d (%d)",daughter1lc,mcPartDaughterBachelor->GetPdgCode())); + } + if (fGenLcOption==kCountK0Sp) { + mcPartDaughterBachelor = dynamic_cast(fmcArray->At(daughter0lc)); // baryon (0) + AliInfo(Form(" Case Lc->K0S+p: (bachelor) daughter0_Lc=%d (%d)",daughter0lc,mcPartDaughterBachelor->GetPdgCode())); + } + + if (!mcPartDaughterV0 || !mcPartDaughterBachelor) return bGenValues; + + if (fGenLcOption==kCountLambdapi) { + if ( !(mcPartDaughterV0->GetPdgCode()==3122 && + mcPartDaughterBachelor->GetPdgCode()==211) ) + AliInfo("It isn't a Lc->Lambda+pion candidate"); + } + if (fGenLcOption==kCountK0Sp) { + if ( !(mcPartDaughterV0->GetPdgCode()==310 && + mcPartDaughterBachelor->GetPdgCode()==2212) ) + AliInfo("It isn't a Lc->K0S+proton candidate"); + } + + Double_t vtx1[3] = {0,0,0}; // primary vertex + fmcPartCandidate->XvYvZv(vtx1); // cm + + // getting vertex from daughters + Double_t vtx1daughter0[3] = {0,0,0}; // secondary vertex from daughter 0 + Double_t vtx1daughter1[3] = {0,0,0}; // secondary vertex from daughter 1 + mcPartDaughterBachelor->XvYvZv(vtx1daughter0); //cm + mcPartDaughterV0->XvYvZv(vtx1daughter1); //cm + if (TMath::Abs(vtx1daughter0[0]-vtx1daughter1[0])>1E-5 || + TMath::Abs(vtx1daughter0[1]-vtx1daughter1[1])>1E-5 || + TMath::Abs(vtx1daughter0[2]-vtx1daughter1[2])>1E-5) { + AliError("Daughters have different secondary vertex, skipping the track"); + return bGenValues; + } + + // getting the momentum from the daughters + Double_t px1[2] = {mcPartDaughterBachelor->Px(), mcPartDaughterV0->Px()}; + Double_t py1[2] = {mcPartDaughterBachelor->Py(), mcPartDaughterV0->Py()}; + Double_t pz1[2] = {mcPartDaughterBachelor->Pz(), mcPartDaughterV0->Pz()}; + + Int_t nprongs = 2; + Short_t charge = mcPartDaughterBachelor->Charge(); + Double_t d0[2] = {0.,0.}; + AliAODRecoDecayHF* decayLc = new AliAODRecoDecayHF(vtx1,vtx1daughter0,nprongs,charge,px1,py1,pz1,d0); + Double_t cosPAwrtPrimVtxLc = decayLc->CosPointingAngle(); + + //V0 daughters + Int_t daughter0 = mcPartDaughterV0->GetDaughter(0); + Int_t daughter1 = mcPartDaughterV0->GetDaughter(1); + AliAODMCParticle* mcPartDaughter0 = dynamic_cast(fmcArray->At(daughter0)); //V0daughter0 + AliAODMCParticle* mcPartDaughter1 = dynamic_cast(fmcArray->At(daughter1)); //V0daughter1 + + if (!mcPartDaughter0 || !mcPartDaughter1) return kFALSE; + + // getting vertex from daughters + Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0 + Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1 + mcPartDaughter0->XvYvZv(vtx2daughter0); //cm + mcPartDaughter1->XvYvZv(vtx2daughter1); //cm + if (TMath::Abs(vtx2daughter0[0]-vtx2daughter1[0])>1E-5 || + TMath::Abs(vtx2daughter0[1]-vtx2daughter1[1])>1E-5 || + TMath::Abs(vtx2daughter0[2]-vtx2daughter1[2])>1E-5) { + AliError("Daughters have different secondary vertex, skipping the track"); + return bGenValues; + } + + // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second + AliAODMCParticle* positiveDaugh = mcPartDaughter0; + AliAODMCParticle* negativeDaugh = mcPartDaughter1; + if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){ + // inverting in case the positive daughter is the second one + positiveDaugh = mcPartDaughter1; + negativeDaugh = mcPartDaughter0; + } + // getting the momentum from the daughters + Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()}; + Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()}; + Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()}; + + nprongs = 2; + charge = 0; + AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0); + Double_t cosPAwrtPrimVtxV0 = decay->CosPointingAngle(); + + //ct + Double_t cTV0 = 0.; + if (fGenLcOption==kCountK0Sp) { + cTV0 = decay->Ct(310); // by default wrt Primary Vtx + } else if (fGenLcOption==kCountLambdapi) { + cTV0 = decay->Ct(3122); // by default wrt Primary Vtx + } + + Double_t cTLc = Ctau(fmcPartCandidate); // by default wrt Primary Vtx + + // get the bachelor pT + Double_t pTbach = 0.; + + if (TMath::Abs(fmcPartCandidate->GetPdgCode()) == 4122) + pTbach = mcPartDaughterBachelor->Pt(); + + Double_t invMass = 0.; + if (fGenLcOption==kCountK0Sp) { + invMass = decay->InvMass2Prongs(0,1,211,211); + } else if (fGenLcOption==kCountLambdapi) { + if (fmcPartCandidate->GetPdgCode() == 4122) + invMass = decay->InvMass2Prongs(0,1,2212,211); + else if (fmcPartCandidate->GetPdgCode() ==-4122) + invMass = decay->InvMass2Prongs(0,1,211,2212); + } + + switch (fConfiguration){ + case AliCFTaskVertexingHF::kSnail: + vectorMC[0] = fmcPartCandidate->Pt(); + vectorMC[1] = fmcPartCandidate->Y() ; + vectorMC[2] = fmcPartCandidate->Phi(); + vectorMC[3] = cosPAwrtPrimVtxV0; + vectorMC[4] = 0; // dummy value x MC, onTheFlyStatus + vectorMC[5] = fCentValue; // reconstructed centrality + vectorMC[6] = 1; // dummy value x MC, fFake + vectorMC[7] = fMultiplicity; // reconstructed multiplicity + + vectorMC[8] = pTbach; + vectorMC[9] = positiveDaugh->Pt(); + vectorMC[10] = negativeDaugh->Pt(); + vectorMC[11] = invMass; + vectorMC[12] = 0; // dummy value x MC, V0 DCA + vectorMC[13] = cTV0*1.E4; // in micron + vectorMC[14] = cTLc*1.E4; // in micron + vectorMC[15] = cosPAwrtPrimVtxLc; + break; + case AliCFTaskVertexingHF::kCheetah: + vectorMC[0] = fmcPartCandidate->Pt(); + vectorMC[1] = fmcPartCandidate->Y() ; + vectorMC[2] = fmcPartCandidate->Phi(); + vectorMC[3] = cosPAwrtPrimVtxV0; + vectorMC[4] = 0; // dummy value x MC, onTheFlyStatus + vectorMC[5] = fCentValue; // reconstructed centrality + vectorMC[6] = 1; // dummy value x MC, fFake + vectorMC[7] = fMultiplicity; // reconstructed multiplicity + break; + } + + delete decay; + delete decayLc; + + bGenValues = kTRUE; + return bGenValues; + +} +//____________________________________________ +Bool_t AliCFVertexingHFLctoV0bachelor::GetRecoValuesFromCandidate(Double_t *vectorReco) const +{ + // read the variables for the container + + Bool_t bFillRecoValues = kFALSE; + + //Get the Lc and the V0 from Lc + AliAODRecoCascadeHF* lcV0bachelor = (AliAODRecoCascadeHF*)fRecoCandidate; + + if ( !(lcV0bachelor->Getv0()) ) return bFillRecoValues; + + AliAODVertex *vtx0 = (AliAODVertex*)lcV0bachelor->GetPrimaryVtx(); + if (vtx0) AliDebug(1,"lcV0bachelor has primary vtx"); + + AliAODTrack* bachelor = (AliAODTrack*)lcV0bachelor->GetBachelor(); + AliAODv0* v0toDaughters = (AliAODv0*)lcV0bachelor->Getv0(); + Bool_t onTheFlyStatus = v0toDaughters->GetOnFlyStatus(); + AliAODTrack* v0positiveTrack = (AliAODTrack*)lcV0bachelor->Getv0PositiveTrack(); + AliAODTrack* v0negativeTrack = (AliAODTrack*)lcV0bachelor->Getv0NegativeTrack(); + + Double_t pt = lcV0bachelor->Pt(); + Double_t rapidity = lcV0bachelor->Y(4122); + Double_t invMassV0 = 0.; + Double_t primVtxPos[3]; vtx0->GetXYZ(primVtxPos); + Double_t cosPAwrtPrimVtxV0 = v0toDaughters->CosPointingAngle(primVtxPos); + Double_t cTV0 = 0.; + + Double_t pTbachelor = bachelor->Pt(); + Double_t pTV0pos = v0positiveTrack->Pt(); + Double_t pTV0neg = v0negativeTrack->Pt(); + Double_t phi = lcV0bachelor->Phi(); + Double_t dcaV0 = v0toDaughters->GetDCA(); + Double_t cTLc = lcV0bachelor->Ct(4122); // wrt PrimVtx + //Double_t dcaLc = lcV0bachelor->GetDCA(); + Double_t cosPointingAngleLc = lcV0bachelor->CosPointingAngle(); + + if (fGenLcOption==kCountK0Sp) { + cTV0 = v0toDaughters->Ct(310,primVtxPos); + } else if (fGenLcOption==kCountLambdapi) { + cTV0 = v0toDaughters->Ct(3122,primVtxPos); + } + + Short_t bachelorCharge = bachelor->Charge(); + if (fGenLcOption==kCountLambdapi) { + if (bachelorCharge==1) { + invMassV0 = v0toDaughters->MassLambda(); + } else if (bachelorCharge==-1) { + invMassV0 = v0toDaughters->MassAntiLambda(); + } + + } else if (fGenLcOption==kCountK0Sp) { + invMassV0 = v0toDaughters->MassK0Short(); + } + + switch (fConfiguration){ + case AliCFTaskVertexingHF::kSnail: + vectorReco[0] = pt; + vectorReco[1] = rapidity; + vectorReco[2] = phi; + vectorReco[3] = cosPAwrtPrimVtxV0; + vectorReco[4] = onTheFlyStatus; + vectorReco[5] = fCentValue; + vectorReco[6] = fFake; // whether the reconstructed candidate was a fake (fFake = 0) or not (fFake = 2) + vectorReco[7] = fMultiplicity; + + vectorReco[8] = pTbachelor; + vectorReco[9] = pTV0pos; + vectorReco[10] = pTV0neg; + vectorReco[11] = invMassV0; + vectorReco[12] = dcaV0; + vectorReco[13] = cTV0*1.E4; // in micron + vectorReco[14] = cTLc*1.E4; // in micron + vectorReco[15] = cosPointingAngleLc; + + break; + case AliCFTaskVertexingHF::kCheetah: + vectorReco[0] = pt; + vectorReco[1] = rapidity; + vectorReco[2] = phi; + vectorReco[3] = cosPAwrtPrimVtxV0; + vectorReco[4] = onTheFlyStatus; + vectorReco[5] = fCentValue; + vectorReco[6] = fFake; // whether the reconstructed candidate was a fake (fFake = 0) or not (fFake = 2) + vectorReco[7] = fMultiplicity; + break; + } + + bFillRecoValues = kTRUE; + + return bFillRecoValues; +} + +//_____________________________________________________________ +Bool_t AliCFVertexingHFLctoV0bachelor::CheckMCChannelDecay() const +{ + // check the required decay channel + + Bool_t checkCD = kFALSE; + + if (fmcPartCandidate->GetNDaughters()!=2) return checkCD; + + Int_t daughter0 = fmcPartCandidate->GetDaughter(0); + Int_t daughter1 = fmcPartCandidate->GetDaughter(1); + AliAODMCParticle* mcPartDaughter0 = dynamic_cast(fmcArray->At(daughter0)); + AliAODMCParticle* mcPartDaughter1 = dynamic_cast(fmcArray->At(daughter1)); + + if (!mcPartDaughter0 || !mcPartDaughter1) { + AliDebug (2,"Problems in the MC Daughters\n"); + return checkCD; + } + + // Lc -> Lambda + pion AND cc + if (fGenLcOption==kCountLambdapi) { + + if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==3122 && + TMath::Abs(mcPartDaughter1->GetPdgCode())==211) && + !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 && + TMath::Abs(mcPartDaughter1->GetPdgCode())==3122)) { + AliDebug(2, "The Lc MC doesn't decay in Lambda+pion (or cc), skipping!!"); + return checkCD; + } else if (TMath::Abs(mcPartDaughter0->GetPdgCode())==3122) { + if (mcPartDaughter0->GetNDaughters()!=2) return checkCD; + Int_t daughter0D0 = mcPartDaughter0->GetDaughter(0); + AliAODMCParticle* mcPartDaughter0D0 = dynamic_cast(fmcArray->At(daughter0D0)); + Int_t daughter0D1 = mcPartDaughter0->GetDaughter(1); + AliAODMCParticle* mcPartDaughter0D1 = dynamic_cast(fmcArray->At(daughter0D1)); + if (!(TMath::Abs(mcPartDaughter0D0->GetPdgCode())==211 && + TMath::Abs(mcPartDaughter0D1->GetPdgCode())==2212) && + !(TMath::Abs(mcPartDaughter0D0->GetPdgCode())==2212 && + TMath::Abs(mcPartDaughter0D1->GetPdgCode())==211)) { + AliDebug(2, "The Lambda MC doesn't decay in pi+proton (or cc), skipping!!"); + return checkCD; + } + } else if (TMath::Abs(mcPartDaughter1->GetPdgCode())==3122) { + if (mcPartDaughter1->GetNDaughters()!=2) return checkCD; + Int_t daughter1D0 = mcPartDaughter1->GetDaughter(0); + AliAODMCParticle* mcPartDaughter1D0 = dynamic_cast(fmcArray->At(daughter1D0)); + Int_t daughter1D1 = mcPartDaughter1->GetDaughter(1); + AliAODMCParticle* mcPartDaughter1D1 = dynamic_cast(fmcArray->At(daughter1D1)); + if (!(TMath::Abs(mcPartDaughter1D0->GetPdgCode())==211 && + TMath::Abs(mcPartDaughter1D1->GetPdgCode())==2212) && + !(TMath::Abs(mcPartDaughter1D0->GetPdgCode())==2212 && + TMath::Abs(mcPartDaughter1D1->GetPdgCode())==211)) { + AliDebug(2, "The Lambda MC doesn't decay in pi+proton (or cc), skipping!!"); + return checkCD; + } + } + + } else if (fGenLcOption==kCountK0Sp) { + // Lc -> K0bar + proton AND cc + + if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==311 && + TMath::Abs(mcPartDaughter1->GetPdgCode())==2212) && + !(TMath::Abs(mcPartDaughter0->GetPdgCode())==2212 && + TMath::Abs(mcPartDaughter1->GetPdgCode())==311)) { + AliDebug(2, "The Lc MC doesn't decay in K0+proton (or cc), skipping!!"); + return checkCD; + } + + if (TMath::Abs(mcPartDaughter0->GetPdgCode())==311) { + Int_t daughter = mcPartDaughter0->GetDaughter(0); + AliAODMCParticle* mcPartDaughter = dynamic_cast(fmcArray->At(daughter)); + if (!(TMath::Abs(mcPartDaughter->GetPdgCode())==310)) { + AliDebug(2, "The K0 (or K0bar) MC doesn't go in K0S, skipping!!"); + return checkCD; + } + if (mcPartDaughter->GetNDaughters()!=2) return checkCD; + Int_t daughterD0 = mcPartDaughter->GetDaughter(0); + AliAODMCParticle* mcPartDaughterD0 = dynamic_cast(fmcArray->At(daughterD0)); + Int_t daughterD1 = mcPartDaughter->GetDaughter(1); + AliAODMCParticle* mcPartDaughterD1 = dynamic_cast(fmcArray->At(daughterD1)); + if (!(TMath::Abs(mcPartDaughterD0->GetPdgCode())==211 && + TMath::Abs(mcPartDaughterD1->GetPdgCode())==211)) { + AliDebug(2, "The K0S MC doesn't decay in pi+pi, skipping!!"); + return checkCD; + } + + } + + if (TMath::Abs(mcPartDaughter1->GetPdgCode())==311) { + Int_t daughter = mcPartDaughter1->GetDaughter(0); + AliAODMCParticle* mcPartDaughter = dynamic_cast(fmcArray->At(daughter)); + if (!(TMath::Abs(mcPartDaughter->GetPdgCode())==310)) { + AliDebug(2, "The K0 (or K0bar) MC doesn't go in K0S, skipping!!"); + return checkCD; + } + if (mcPartDaughter->GetNDaughters()!=2) return checkCD; + Int_t daughterD0 = mcPartDaughter->GetDaughter(0); + AliAODMCParticle* mcPartDaughterD0 = dynamic_cast(fmcArray->At(daughterD0)); + Int_t daughterD1 = mcPartDaughter->GetDaughter(1); + AliAODMCParticle* mcPartDaughterD1 = dynamic_cast(fmcArray->At(daughterD1)); + if (! ( TMath::Abs(mcPartDaughterD0->GetPdgCode())==211 && + TMath::Abs(mcPartDaughterD1->GetPdgCode())==211 ) ) { + AliDebug(2, "The K0S MC doesn't decay in pi+pi, skipping!!"); + return checkCD; + } + + } + + } + + checkCD = kTRUE; + return checkCD; + +} + +//___________________________________________________________ + +void AliCFVertexingHFLctoV0bachelor::SetPtAccCut(Float_t* ptAccCut) +{ + // + // setting the pt cut to be used in the Acceptance steps (MC+Reco) + // + + AliInfo("The 1st element of the pt cut array will correspond to the cut applied to the bachelor - please check that it is correct"); + if (fProngs>0){ + for (Int_t iP=0; iP0){ + for (Int_t iP=0; iP0){ + for (Int_t iP=0; iP(fmcArray->At(fLabelArray[0])); // bachelor + if (!mcPartDaughter) return; + Int_t mother = mcPartDaughter->GetMother(); + AliAODMCParticle* mcMother = dynamic_cast(fmcArray->At(mother)); + if (!mcMother) return; + /* + if (fGenLcOption==kCountK0Sp) { + if ( !(TMath::Abs(mcPartDaughter->GetPdgCode())== 2212) ) + AliFatal(Form("Apparently the proton bachelor is not in the first position (%d <- %d), causing a crash!!", + mcPartDaughter->GetPdgCode(),mcMother->GetPdgCode())); + } + else if (fGenLcOption==kCountLambdapi) { + if ( !(TMath::Abs(mcPartDaughter->GetPdgCode())== 211) ) + AliFatal(Form("Apparently the pion bachelor is not in the first position (%d <- %d), causing a crash!!", + mcPartDaughter->GetPdgCode(),mcMother->GetPdgCode())); + } + */ + if (fProngs>0) { + fPtAccCut[0]=0.3; // bachelor + fEtaAccCut[0]=0.9; // bachelor + for (Int_t iP=1; iPGetBachelor()->Eta(); + else if (iProng==1) etaProng = lcV0bachelor->Getv0PositiveTrack()->Eta(); + else if (iProng==2) etaProng = lcV0bachelor->Getv0NegativeTrack()->Eta(); + + return etaProng; + + } + return 999999; +} + +//_____________________________________________________________ + +Double_t AliCFVertexingHFLctoV0bachelor::GetPtProng(Int_t iProng) const +{ + // + // getting pt of the prong + // + + Double_t ptProng=-9999.; + + if (!fRecoCandidate) return ptProng; + + AliAODRecoCascadeHF* lcV0bachelor = (AliAODRecoCascadeHF*)fRecoCandidate; + + if (iProng==0) ptProng = lcV0bachelor->GetBachelor()->Pt(); + else if (iProng==1) ptProng = lcV0bachelor->Getv0PositiveTrack()->Pt(); + else if (iProng==2) ptProng = lcV0bachelor->Getv0NegativeTrack()->Pt(); + + return ptProng; + +} + +//_____________________________________________________________ + +Double_t AliCFVertexingHFLctoV0bachelor::Ctau(AliAODMCParticle *mcPartCandidate) +{ + + Double_t cTau = 999999.; + + Double_t vtx1[3] = {0,0,0}; // primary vertex + Bool_t hasProdVertex = mcPartCandidate->XvYvZv(vtx1); // cm + + AliAODMCParticle *mcPartDaughter0 = (AliAODMCParticle*)fmcArray->At(mcPartCandidate->GetDaughter(0)); + AliAODMCParticle *mcPartDaughter1 = (AliAODMCParticle*)fmcArray->At(mcPartCandidate->GetDaughter(1)); + if (!mcPartDaughter0 && !mcPartDaughter1) return cTau; + Double_t vtx1daughter[3] = {0,0,0}; // secondary vertex + if (mcPartDaughter0) + hasProdVertex = hasProdVertex || mcPartDaughter0->XvYvZv(vtx1daughter); //cm + if (mcPartDaughter1) + hasProdVertex = hasProdVertex || mcPartDaughter1->XvYvZv(vtx1daughter); //cm + + if (!hasProdVertex) return cTau; + + Double_t decayLength = 0.; + for (Int_t ii=0; ii<3; ii++) decayLength += (vtx1daughter[ii]-vtx1[ii])*(vtx1daughter[ii]-vtx1[ii]); + decayLength = TMath::Sqrt(decayLength); + + cTau = decayLength * mcPartCandidate->M()/mcPartCandidate->P(); + AliDebug(2,Form(" cTau(4122)=%f",cTau)); + return cTau; + +} diff --git a/PWGHF/vertexingHF/AliCFVertexingHFLctoV0bachelor.h b/PWGHF/vertexingHF/AliCFVertexingHFLctoV0bachelor.h new file mode 100644 index 00000000000..326f2e72e51 --- /dev/null +++ b/PWGHF/vertexingHF/AliCFVertexingHFLctoV0bachelor.h @@ -0,0 +1,83 @@ +#ifndef ALICFVERTEXINGHFLCTOV0BACHELOR_H +#define ALICFVERTEXINGHFLCTOV0BACHELOR_H + +/************************************************************************** + * Copyright(c) 1998-2011, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id$ */ + +//----------------------------------------------------------------------- +// Class for HF corrections as a function of many variables and steps +// For Lc->V0 bachelor +// +// Author : A. De Caro +//----------------------------------------------------------------------- + + +#include "AliCFVertexingHF.h" +#include "AliAODRecoDecayHF.h" +#include "AliAODRecoCascadeHF.h" +#include "AliAODRecoDecayHF2Prong.h" + +class AliAODMCParticle; +class TClonesArray; +class AliCFVertexingHF; +class AliESDtrack; +class TDatabasePDG; + +class AliCFVertexingHFLctoV0bachelor : public AliCFVertexingHF{ + public: + + enum ELctoV0Sel { + kCountK0Sp=0, + kCountLambdapi=1, + kCountAllLctoV0=2 + }; + + AliCFVertexingHFLctoV0bachelor(); + AliCFVertexingHFLctoV0bachelor(TClonesArray *mcArray, UShort_t originDselection, Int_t lcDecay=0); + + virtual ~AliCFVertexingHFLctoV0bachelor(){}; + + Bool_t GetGeneratedValuesFromMCParticle(Double_t* /*vectorMC*/); + Bool_t GetRecoValuesFromCandidate(Double_t* /*vectorReco*/ ) const; + Bool_t CheckMCChannelDecay()const; + + Bool_t SetRecoCandidateParam(AliAODRecoDecayHF *recoCand); + + void SetPtAccCut(Float_t* ptAccCut); + void SetEtaAccCut(Float_t* etaAccCut); + void SetAccCut(Float_t* ptAccCut, Float_t* etaAccCut); + void SetAccCut(); + + Double_t GetEtaProng(Int_t iProng)const; + Double_t GetPtProng(Int_t iProng) const; + + protected: + + Double_t Ctau(AliAODMCParticle *mcPartCandidate); + + private: + + AliCFVertexingHFLctoV0bachelor(const AliCFVertexingHFLctoV0bachelor& c); + AliCFVertexingHFLctoV0bachelor& operator= (const AliCFVertexingHFLctoV0bachelor& other); + + Int_t fGenLcOption; // option for selection Lc (see enum) + + ClassDef(AliCFVertexingHFLctoV0bachelor, 1); // CF class for Lc->V0+bachelor and other cascades + +}; + +#endif diff --git a/PWGHF/vertexingHF/macros/AddTaskCFVertexingHFLctoV0bachelor.C b/PWGHF/vertexingHF/macros/AddTaskCFVertexingHFLctoV0bachelor.C new file mode 100644 index 00000000000..cb7f6fd2b3b --- /dev/null +++ b/PWGHF/vertexingHF/macros/AddTaskCFVertexingHFLctoV0bachelor.C @@ -0,0 +1,691 @@ +//DEFINITION OF A FEW CONSTANTS +const Double_t ptmin = 0.0; +const Double_t ptmax = 6.0; +const Double_t ymin = -1.2 ; +const Double_t ymax = 1.2 ; +const Double_t cosPAV0min = 0.95; +const Double_t cosPAV0max = 1.02; +const Int_t onFlymin = 0; +const Int_t onFlymax = 2; +const Float_t centmin = 0.; +//const Float_t centmin_0_10 = 0.; +//const Float_t centmax_0_10 = 10.; +//const Float_t centmin_10_60 = 10.; +//const Float_t centmax_10_60 = 60.; +//const Float_t centmin_60_100 = 60.; +//const Float_t centmax_60_100 = 100.; +const Float_t centmax = 100.; +const Int_t fakemin = 0; +const Int_t fakemax = 3; +const Float_t multmin = 0; +//const Float_t multmin_0_20 = 0; +//const Float_t multmax_0_20 = 20; +//const Float_t multmin_20_50 = 20; +//const Float_t multmax_20_50 = 50; +//const Float_t multmin_50_102 = 50; +//const Float_t multmax_50_102 = 102; +const Float_t multmax = 102; + +const Double_t ptBachmin = 0.0; +const Double_t ptBachmax = 5.0; +const Double_t ptV0posmin = 0.0; +const Double_t ptV0posmax = 5.0; +const Double_t ptV0negmin = 0.0; +const Double_t ptV0negmax = 5.0; +const Double_t dcaV0min = 0.0; // nSigma +const Double_t dcaV0max = 1.5; // nSigma +const Double_t cTV0min = 0.0; // micron +const Double_t cTV0max = 300; // micron +const Double_t cTmin = 0.0; // micron +const Double_t cTmax = 300; // micron +const Float_t cosPAmin =-1.02; +const Float_t cosPAmax = 1.02; + +const Double_t etamin = -0.9; +const Double_t etamax = 0.9; +//const Double_t zmin = -15.; +//const Double_t zmax = 15.; + + +//---------------------------------------------------- + +AliCFTaskVertexingHF *AddTaskCFVertexingHFLctoV0bachelor(const char* cutFile = "./LctoV0bachelorCuts.root", + Int_t configuration = AliCFTaskVertexingHF::kCheetah, Bool_t isKeepDfromB = kTRUE, + Bool_t isKeepDfromBOnly = kFALSE, Int_t pdgCode = 4122, Char_t isSign = 0, TString usercomment = "username") +{ + + + printf("Adding CF task using cuts from file %s\n",cutFile); + if (configuration == AliCFTaskVertexingHF::kSnail){ + printf("The configuration is set to be SLOW --> all the variables will be used to fill the CF\n"); + } + else if (configuration == AliCFTaskVertexingHF::kCheetah){ + printf("The configuration is set to be FAST --> using only pt, y, ct, phi, zvtx, centrality, fake, multiplicity to fill the CF\n"); + } + else{ + printf("The configuration is not defined! returning\n"); + return; + } + + gSystem->Sleep(2000); + + // isSign = 0 --> K0S only + // isSign = 1 --> Lambda only + // isSign = 2 --> LambdaBar only + // isSign = 3 --> K0S and Lambda and LambdaBar + + TString expected; + if (isSign == 0 && pdgCode < 0){ + AliError(Form("Error setting PDG code (%d) and sign (0 --> K0S only): they are not compatible, returning")); + return 0x0; + } + else if (isSign == 1 && pdgCode < 0){ + AliError(Form("Error setting PDG code (%d) and sign (1 --> Lambda only): they are not compatible, returning")); + return 0x0; + } + else if (isSign == 2 && pdgCode > 0){ + AliError(Form("Error setting PDG code (%d) and sign (2 --> LambdaBar only): they are not compatible, returning")); + return 0x0; + } + else if (isSign > 3 || isSign < 0){ + AliError(Form("Sign not valid (%d, possible values are 0, 1, 2, 3), returning")); + return 0x0; + } + + TFile* fileCuts = TFile::Open(cutFile); + AliRDHFCuts *cutsLctoV0 = (AliRDHFCutsLctoV0*)fileCuts->Get("LctoV0AnalysisCuts"); + + // check that the fKeepD0fromB flag is set to true when the fKeepD0fromBOnly flag is true + // for now the binning is the same than for all D's + if (isKeepDfromBOnly) isKeepDfromB = true; + + Double_t massV0min = 0.47; + Double_t massV0max = 1.14; + if (isSign==0) { + massV0min = 0.47 ; + massV0max = 0.53 ; + } else if (isSign==1 || isSign==2) { + massV0min = 1.09; + massV0max = 1.14; + } + + const Double_t phimin = 0.0; + const Double_t phimax = 2.*TMath::Pi(); + + const Int_t nbinpt = 6; //bins in pt from 0 to 6 GeV + const Int_t nbiny = 24; //bins in y + Int_t nbininvMassV0 = 60; //bins in invMassV0 + if (isSign==3) nbininvMassV0=134; //bins in invMassV0 + const Int_t nbinpointingV0 = 12; //bins in cosPointingAngleV0 + const Int_t nbinonFly = 2; //bins in onFlyStatus x V0 + + const Int_t nbincent = 18; //bins in centrality (total number) + //const Int_t nbincent_0_10 = 4; //bins in centrality between 0 and 10 + //const Int_t nbincent_10_60 = 10; //bins in centrality between 10 and 60 + //const Int_t nbincent_60_100 = 4; //bins in centrality between 60 and 100 + + const Int_t nbinfake = 3; //bins in fake + + const Int_t nbinmult = 48; //bins in multiplicity (total number) + //const Int_t nbinmult_0_20 = 20; //bins in multiplicity between 0 and 20 + //const Int_t nbinmult_20_50 = 15; //bins in multiplicity between 20 and 50 + //const Int_t nbinmult_50_102 = 13; //bins in multiplicity between 50 and 102 + + + const Int_t nbinptBach = 50; //bins in pt from 0 to 5 GeV + const Int_t nbinptV0pos = 50; //bins in pt from 0 to 5 GeV + const Int_t nbinptV0neg = 50; //bins in pt from 0 to 5 GeV + const Int_t nbinphi = 18; //bins in Phi + const Int_t nbindcaV0 = 15; //bins in dcaV0 + const Int_t nbincTV0 = 15; //bins in cTV0 + const Int_t nbincT = 15; //bins in cT + const Int_t nbinpointing = 51; //bins in cosPointingAngle + + //the sensitive variables, their indices + + // variables' indices + const UInt_t ipT = 0; + const UInt_t iy = 1; + const UInt_t iphi = 2; + const UInt_t icosPAxV0 = 3; + const UInt_t ionFly = 4; + const UInt_t imult = 5; + const UInt_t ifake = 6; + const UInt_t icent = 7; + + const UInt_t ipTbach = 8; + const UInt_t ipTposV0 = 9; + const UInt_t ipTnegV0 = 10; + const UInt_t iinvMassV0= 11; + const UInt_t idcaV0 = 12; + const UInt_t icTv0 = 13; + const UInt_t icT = 14; + const UInt_t icosPA = 15; + + //Setting the bins: pt, ptPi, and ptK are considered seprately because for them you can either define the binning by hand, or using the cuts file + + //arrays for the number of bins in each dimension + + //if ( configuration ==AliCFTaskVertexingHF::kSnail) + const Int_t nvarTot = 16 ; //number of variables on the grid:pt, y, cosThetaStar, pTpi, pTk, cT, dca, d0pi, d0K, d0xd0, cosPointingAngle, phi, z, centrality, fake, cosPointingAngleXY, normDecayLengthXY, multiplicity + //if ( configuration ==AliCFTaskVertexingHF::kCheetah) + //const Int_t nvarTot = 8 ; //number of variables on the grid:pt, y, cosThetaStar, pTpi, pTk, cT, dca, d0pi, d0K, d0xd0, cosPointingAngle, phi, z, centrality, fake, cosPointingAngleXY, normDecayLengthXY, multiplicity + + Int_t iBin[nvarTot]; + + //OPTION 1: defining the pt, ptPi, ptK bins by hand... + iBin[ipT]=nbinpt; + iBin[iy]=nbiny; + iBin[iphi]=nbinphi; + iBin[icosPAxV0]=nbinpointingV0; + iBin[ionFly]=nbinonFly; + iBin[imult]=nbinmult; + iBin[ifake]=nbinfake; + iBin[icent]=nbincent; + + iBin[ipTbach]=nbinptBach; + iBin[ipTposV0]=nbinptV0pos; + iBin[ipTnegV0]=nbinptV0neg; + iBin[iinvMassV0]=nbininvMassV0; + iBin[idcaV0]=nbindcaV0; + iBin[icTv0]=nbincTV0; + iBin[icT]=nbincT; + iBin[icosPA]=nbinpointing; + + // values for bin lower bounds + + // pt + Double_t *binLimpT=new Double_t[iBin[0]+1]; + for(Int_t i=0; i<=iBin[0]; i++) binLimpT[i]=(Double_t)ptmin + (ptmax-ptmin)/iBin[0]*(Double_t)i ; + + // y + Double_t *binLimy=new Double_t[iBin[1]+1]; + for(Int_t i=0; i<=iBin[1]; i++) binLimy[i]=(Double_t)ymin + (ymax-ymin)/iBin[1]*(Double_t)i ; + + // phi + Double_t *binLimphi=new Double_t[iBin[2]+1]; + for(Int_t i=0; i<=iBin[2]; i++) binLimphi[i]=(Double_t)phimin + (phimax-phimin)/iBin[2]*(Double_t)i ; + + // cosPointingAngleV0 + Double_t *binLimcosPAV0=new Double_t[iBin[3]+1]; + for(Int_t i=0; i<=iBin[3]; i++) binLimcosPAV0[i]=(Double_t)cosPAV0min + (cosPAV0max-cosPAV0min)/iBin[3]*(Double_t)i ; + + // onTheFlyV0 + Double_t *binLimonFlyV0=new Double_t[iBin[4]+1]; + for(Int_t i=0; i<=iBin[4]; i++) binLimonFlyV0[i]=(Double_t)onFlymin + (onFlymax-onFlymin)/iBin[4]*(Double_t)i ; + + // centrality + Double_t *binLimcent=new Double_t[iBin[5]+1]; + for(Int_t i=0; i<=iBin[5]; i++) binLimcent[i]=(Double_t)centmin + (centmax-centmin)/iBin[5]*(Double_t)i ; + /* + Double_t *binLimcent_0_10=new Double_t[iBin[icent]+1]; + for(Int_t i=0; i<=nbincent_0_10; i++) binLimcent[i]=(Double_t)centmin_0_10 + (centmax_0_10-centmin_0_10)/nbincent_0_10*(Double_t)i ; + if (binLimcent[nbincent_0_10] != centmin_10_60) { + Error("AliCFHeavyFlavourTaskMultiVarMultiStep","Calculated bin lim for cent - 1st range - differs from expected!\n"); + } + + Double_t *binLimcent_10_60=new Double_t[iBin[icent]+1]; + for(Int_t i=0; i<=nbincent_10_60; i++) binLimcent[i+nbincent_0_10]=(Double_t)centmin_10_60 + (centmax_10_60-centmin_10_60)/nbincent_10_60*(Double_t)i ; + if (binLimcent[nbincent_0_10+nbincent_10_60] != centmin_60_100) { + Error("AliCFHeavyFlavourTaskMultiVarMultiStep","Calculated bin lim for cent - 2st range - differs from expected!\n"); + } + + Double_t *binLimcent_60_100=new Double_t[iBin[icent]+1]; + for(Int_t i=0; i<=nbincent_60_100; i++) binLimcent[i+nbincent_10_60]=(Double_t)centmin_60_100 + (centmax_60_100-centmin_60_100)/nbincent_60_100*(Double_t)i ; + */ + + // fake + Double_t *binLimfake=new Double_t[iBin[6]+1]; + for(Int_t i=0; i<=iBin[6]; i++) binLimfake[i]=(Double_t)fakemin + (fakemax-fakemin)/iBin[6] * (Double_t)i; + + // multiplicity + Double_t *binLimmult=new Double_t[iBin[7]+1]; + for(Int_t i=0; i<=iBin[7]; i++) binLimmult[i]=(Double_t)multmin + (multmax-multmin)/iBin[7]*(Double_t)i ; + /* + for(Int_t i=0; i<=nbinmult_0_20; i++) binLimmult[i]=(Double_t)multmin_0_20 + (multmax_0_20-multmin_0_20)/nbinmult_0_20*(Double_t)i ; + if (binLimmult[nbinmult_0_20] != multmin_20_50) { + Error("AliCFHeavyFlavourTaskMultiVarMultiStep","Calculated bin lim for mult - 1st range - differs from expected!\n"); + } + for(Int_t i=0; i<=nbinmult_20_50; i++) binLimmult[i+nbinmult_0_20]=(Double_t)multmin_20_50 + (multmax_20_50-multmin_20_50)/nbinmult_20_50*(Double_t)i ; + if (binLimmult[nbinmult_0_20+nbinmult_20_50] != multmin_50_102) { + Error("AliCFHeavyFlavourTaskMultiVarMultiStep","Calculated bin lim for mult - 2nd range - differs from expected!\n"); + } + for(Int_t i=0; i<=nbinmult_50_102; i++) binLimmult[i+nbinmult_0_20+nbinmult_20_50]=(Double_t)multmin_50_102 + (multmax_50_102-multmin_50_102)/nbinmult_50_102*(Double_t)i ; + */ + + + // ptBach + Double_t *binLimpTbach=new Double_t[iBin[8]+1]; + for(Int_t i=0; i<=iBin[8]; i++) binLimpTbach[i]=(Double_t)ptBachmin + (ptBachmax-ptBachmin)/iBin[8]*(Double_t)i ; + + // ptV0pos + Double_t *binLimpTV0pos=new Double_t[iBin[9]+1]; + for(Int_t i=0; i<=iBin[9]; i++) binLimpTV0pos[i]=(Double_t)ptV0posmin + (ptV0posmax-ptV0posmin)/iBin[9]*(Double_t)i ; + + // ptV0neg + Double_t *binLimpTV0neg=new Double_t[iBin[10]+1]; + for(Int_t i=0; i<=iBin[10]; i++) binLimpTV0neg[i]=(Double_t)ptV0negmin + (ptV0negmax-ptV0negmin)/iBin[10]*(Double_t)i ; + + // invMassV0 + Double_t *binLimInvMassV0=new Double_t[iBin[11]+1]; + for(Int_t i=0; i<=iBin[11]; i++) binLimInvMassV0[i]=(Double_t)massV0min + (massV0max-massV0min)/iBin[11]*(Double_t)i ; + + // dcaV0 + Double_t *binLimdcaV0=new Double_t[iBin[12]+1]; + for(Int_t i=0; i<=iBin[12]; i++) binLimdcaV0[i]=(Double_t)dcaV0min + (dcaV0max-dcaV0min)/iBin[12]*(Double_t)i ; + + // cTV0 + Double_t *binLimcTV0=new Double_t[iBin[13]+1]; + for(Int_t i=0; i<=iBin[13]; i++) binLimcTV0[i]=(Double_t)cTV0min + (cTV0max-cTV0min)/iBin[13]*(Double_t)i ; + + // cT + Double_t *binLimcT=new Double_t[iBin[14]+1]; + for(Int_t i=0; i<=iBin[14]; i++) binLimcT[i]=(Double_t)cTmin + (cTmax-cTmin)/iBin[14]*(Double_t)i ; + + // cosPointingAngle + Double_t *binLimcosPA=new Double_t[iBin[15]+1]; + for(Int_t i=0; i<=iBin[15]; i++) binLimcosPA[i]=(Double_t)cosPAmin + (cosPAmax-cosPAmin)/iBin[15]*(Double_t)i ; + + + + // z Primary Vertex + //for(Int_t i=0; i<=nbinzvtx; i++) binLimzvtx[i]=(Double_t)zmin + (zmax-zmin) /nbinzvtx*(Double_t)i ; + + + //OPTION 2: ...or from the cuts file + /* + const Int_t nbinpt = cutsLctoV0->GetNPtBins(); // bins in pT + iBin[ipT]=nbinpt; + iBin[ipTpi]=nbinpt; + iBin[ipTk]=nbinpt; + Double_t *binLimpT=new Double_t[iBin[ipT]+1]; + Double_t *binLimpTpi=new Double_t[iBin[ipTpi]+1]; + Double_t *binLimpTk=new Double_t[iBin[ipTk]+1]; + // values for bin lower bounds + Float_t* floatbinLimpT = cutsLctoV0->GetPtBinLimits(); + for (Int_t ibin0 = 0 ; ibin0 SetBinLimits(0,binLimpT); + container -> SetBinLimits(1,binLimy); + container -> SetBinLimits(2,binLimphi); + container -> SetBinLimits(3,binLimcosPAV0); + container -> SetBinLimits(4,binLimonFlyV0); + container -> SetBinLimits(5,binLimcent); + container -> SetBinLimits(6,binLimfake); + container -> SetBinLimits(7,binLimmult); + + container -> SetVarTitle(0,"pt"); + container -> SetVarTitle(1,"y"); + container -> SetVarTitle(2,"phi"); + container -> SetVarTitle(3,"cosPA -V0-"); + container -> SetVarTitle(4,"onFlyV0"); + container -> SetVarTitle(5,"centrality"); + container -> SetVarTitle(6,"fake"); + container -> SetVarTitle(7,"multiplicity"); + + if (configuration == AliCFTaskVertexingHF::kSnail) { + container -> SetBinLimits(8,binLimpTbach); + container -> SetBinLimits(9,binLimpTV0pos); + container -> SetBinLimits(10,binLimpTV0neg); + container -> SetBinLimits(11,binLimInvMassV0); + container -> SetBinLimits(12,binLimdcaV0); + container -> SetBinLimits(13,binLimcTV0); + container -> SetBinLimits(14,binLimcT); + container -> SetBinLimits(15,binLimcosPA); + + container -> SetVarTitle(8,"ptBachelor"); + container -> SetVarTitle(9,"ptV0pos"); + container -> SetVarTitle(10,"ptV0neg"); + container -> SetVarTitle(11,"mV0"); + container -> SetVarTitle(12,"DCA -V0-"); + container -> SetVarTitle(13,"c#tau -V0-"); + container -> SetVarTitle(14,"c#tau"); + container -> SetVarTitle(15,"cosPA"); + } + + container -> SetStepTitle(0, "MCLimAcc"); + container -> SetStepTitle(1, "MC"); + container -> SetStepTitle(2, "MCAcc"); + container -> SetStepTitle(3, "RecoVertex"); + container -> SetStepTitle(4, "RecoRefit"); + container -> SetStepTitle(5, "Reco"); + container -> SetStepTitle(6, "RecoAcc"); + container -> SetStepTitle(7, "RecoITSCluster"); + container -> SetStepTitle(8, "RecoCuts"); + container -> SetStepTitle(9, "RecoPID"); + + //return container; + + //CREATE THE CUTS ----------------------------------------------- + + // Gen-Level kinematic cuts + AliCFTrackKineCuts *mcKineCuts = new AliCFTrackKineCuts("mcKineCuts","MC-level kinematic cuts"); + + //Particle-Level cuts: + AliCFParticleGenCuts* mcGenCuts = new AliCFParticleGenCuts("mcGenCuts","MC particle generation cuts"); + Bool_t useAbsolute = kTRUE; + if (isSign != 3 && isSign!=0) { + useAbsolute = kFALSE; + } + mcGenCuts->SetRequirePdgCode(pdgCode, useAbsolute); // kTRUE set in order to include Lc- + mcGenCuts->SetAODMC(1); //special flag for reading MC in AOD tree (important) + + // Acceptance cuts: + AliCFAcceptanceCuts* accCuts = new AliCFAcceptanceCuts("accCuts", "Acceptance cuts"); + AliCFTrackKineCuts * kineAccCuts = new AliCFTrackKineCuts("kineAccCuts","Kine-Acceptance cuts"); + kineAccCuts->SetPtRange(ptmin,ptmax); + kineAccCuts->SetEtaRange(etamin,etamax); + + // Rec-Level kinematic cuts + AliCFTrackKineCuts *recKineCuts = new AliCFTrackKineCuts("recKineCuts","rec-level kine cuts"); + + AliCFTrackQualityCuts *recQualityCuts = new AliCFTrackQualityCuts("recQualityCuts","rec-level quality cuts"); + + AliCFTrackIsPrimaryCuts *recIsPrimaryCuts = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts","rec-level isPrimary cuts"); + + printf("CREATE MC KINE CUTS\n"); + TObjArray* mcList = new TObjArray(0) ; + mcList->AddLast(mcKineCuts); + mcList->AddLast(mcGenCuts); + + printf("CREATE ACCEPTANCE CUTS\n"); + TObjArray* accList = new TObjArray(0) ; + accList->AddLast(kineAccCuts); + + printf("CREATE RECONSTRUCTION CUTS\n"); + TObjArray* recList = new TObjArray(0) ; // not used!! + recList->AddLast(recKineCuts); + recList->AddLast(recQualityCuts); + recList->AddLast(recIsPrimaryCuts); + + TObjArray* emptyList = new TObjArray(0); + + //CREATE THE INTERFACE TO CORRECTION FRAMEWORK USED IN THE TASK + printf("CREATE INTERFACE AND CUTS\n"); + AliCFManager* man = new AliCFManager() ; + man->SetParticleContainer(container); + man->SetParticleCutsList(0 , mcList); // MC, Limited Acceptance + man->SetParticleCutsList(1 , mcList); // MC + man->SetParticleCutsList(2 , accList); // Acceptance + man->SetParticleCutsList(3 , emptyList); // Vertex + man->SetParticleCutsList(4 , emptyList); // Refit + man->SetParticleCutsList(5 , emptyList); // AOD + man->SetParticleCutsList(6 , emptyList); // AOD in Acceptance + man->SetParticleCutsList(7 , emptyList); // AOD with required n. of ITS clusters + man->SetParticleCutsList(8 , emptyList); // AOD Reco (PPR cuts implemented in Task) + man->SetParticleCutsList(9 , emptyList); // AOD Reco PID + + // Get the pointer to the existing analysis manager via the static access method. + //============================================================================== + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + ::Error("AddTaskCompareHF", "No analysis manager to connect to."); + return NULL; + } + //CREATE THE TASK + printf("CREATE TASK\n"); + + // create the task + AliCFTaskVertexingHF *task = new AliCFTaskVertexingHF("AliCFTaskVertexingHF",cutsLctoV0); + task->SetConfiguration(configuration); + task->SetFillFromGenerated(kFALSE); + task->SetCFManager(man); //here is set the CF manager + task->SetDecayChannel(22);//kLctoV0bachelor + task->SetUseWeight(kFALSE); + task->SetSign(isSign); + task->SetCentralitySelection(kFALSE); + task->SetFakeSelection(0); + task->SetRejectCandidateIfNotFromQuark(kTRUE); // put to false if you want to keep HIJING D0!! + task->SetUseMCVertex(kFALSE); // put to true if you want to do studies on pp + + if (isKeepDfromB && !isKeepDfromBOnly) task->SetDselection(2); + if (isKeepDfromB && isKeepDfromBOnly) task->SetDselection(1); + + TF1* funcWeight = 0x0; + if (task->GetUseWeight()) { + funcWeight = (TF1*)fileCuts->Get("funcWeight"); + if (funcWeight == 0x0){ + Printf("FONLL Weights will be used"); + } + else { + task->SetWeightFunction(funcWeight); + Printf("User-defined Weights will be used. The function being:"); + task->GetWeightFunction(funcWeight)->Print(); + } + } + + Printf("***************** CONTAINER SETTINGS *****************"); + Printf("decay channel = %d",(Int_t)task->GetDecayChannel()); + Printf("FillFromGenerated = %d",(Int_t)task->GetFillFromGenerated()); + Printf("Dselection = %d",(Int_t)task->GetDselection()); + Printf("UseWeight = %d",(Int_t)task->GetUseWeight()); + if (task->GetUseWeight()) { + Printf("User-defined Weight function:"); + task->GetWeightFunction(funcWeight)->Print(); + } + else { + Printf("FONLL will be used for the weights"); + } + Printf("Sign = %d",(Int_t)task->GetSign()); + Printf("Centrality selection = %d",(Int_t)task->GetCentralitySelection()); + Printf("Fake selection = %d",(Int_t)task->GetFakeSelection()); + Printf("RejectCandidateIfNotFromQuark selection = %d",(Int_t)task->GetRejectCandidateIfNotFromQuark()); + Printf("UseMCVertex selection = %d",(Int_t)task->GetUseMCVertex()); + Printf("***************END CONTAINER SETTINGS *****************\n"); + + //-----------------------------------------------------------// + // create correlation matrix for unfolding - only eta-pt // + //-----------------------------------------------------------// + + Bool_t AcceptanceUnf = kTRUE; // unfold at acceptance level, otherwise PPR + + Int_t thnDim[4]; + + //first half : reconstructed + //second half : MC + + thnDim[0] = iBin[0]; + thnDim[2] = iBin[0]; + thnDim[1] = iBin[1]; + thnDim[3] = iBin[1]; + + TString nameCorr=""; + if (!isKeepDfromB) { + nameCorr="CFHFcorr0_CommonFramework_"+usercomment; + } + else if (isKeepDfromBOnly) { + nameCorr= "CFHFcorr0KeepDfromBOnly_CommonFramework_"+usercomment; + } + else { + nameCorr="CFHFcorr0allLc_CommonFramework_"+usercomment; + } + + THnSparseD* correlation = new THnSparseD(nameCorr,"THnSparse with correlations",4,thnDim); + Double_t** binEdges = new Double_t[2]; + + // set bin limits + + binEdges[0]= binLimpT; + binEdges[1]= binLimy; + + correlation->SetBinEdges(0,binEdges[0]); + correlation->SetBinEdges(2,binEdges[0]); + + correlation->SetBinEdges(1,binEdges[1]); + correlation->SetBinEdges(3,binEdges[1]); + + correlation->Sumw2(); + + // correlation matrix ready + //------------------------------------------------// + + task->SetCorrelationMatrix(correlation); // correlation matrix for unfolding + + // Create and connect containers for input/output + + // ------ input data ------ + AliAnalysisDataContainer *cinput0 = mgr->GetCommonInputContainer(); + + // ----- output data ----- + + TString outputfile = AliAnalysisManager::GetCommonFileName(); + TString output1name="", output2name="", output3name="",output4name=""; + output2name=nameContainer; + output3name=nameCorr; + if (!isKeepDfromB) { + outputfile += ":PWG3_D2H_CFtaskLctoK0Sp_CommonFramework_"+usercomment; + output1name="CFHFchist0_CommonFramework_"+usercomment; + output4name= "Cuts_CommonFramework_"+usercomment; + } + else if (isKeepDfromBOnly) { + outputfile += ":PWG3_D2H_CFtaskLctoK0SpKeepDfromBOnly_CommonFramework_"+usercomment; + output1name="CFHFchist0DfromB_CommonFramework_"+usercomment; + output4name= "Cuts_CommonFramework_DfromB_"+usercomment; + } + else { + outputfile += ":PWG3_D2H_CFtaskLctoK0SpKeepDfromB_CommonFramework_"+usercomment; + output1name="CFHFchist0allLc_CommonFramework_"+usercomment; + output4name= "Cuts_CommonFramework_allLc_"+usercomment; + } + + //now comes user's output objects : + // output TH1I for event counting + AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(output1name, TH1I::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data()); + // output Correction Framework Container (for acceptance & efficiency calculations) + AliAnalysisDataContainer *coutput2 = mgr->CreateContainer(output2name, AliCFContainer::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data()); + // Unfolding - correlation matrix + AliAnalysisDataContainer *coutput3 = mgr->CreateContainer(output3name, THnSparseD::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data()); + // cuts + AliAnalysisDataContainer *coutput4 = mgr->CreateContainer(output4name, AliRDHFCuts::Class(),AliAnalysisManager::kOutputContainer, outputfile.Data()); + + mgr->AddTask(task); + + mgr->ConnectInput(task,0,mgr->GetCommonInputContainer()); + mgr->ConnectOutput(task,1,coutput1); + mgr->ConnectOutput(task,2,coutput2); + mgr->ConnectOutput(task,3,coutput3); + mgr->ConnectOutput(task,4,coutput4); + return task; + +} -- 2.43.0