diff --git a/PWGHF/Core/HFSelectorCuts.h b/PWGHF/Core/HFSelectorCuts.h index 5f18d65cfcc..f86b4232ef5 100644 --- a/PWGHF/Core/HFSelectorCuts.h +++ b/PWGHF/Core/HFSelectorCuts.h @@ -27,6 +27,7 @@ enum Code { kD0 = 421, kD0bar = -421, kDPlus = 411, + kDs = 431, kLambdaCPlus = 4122, kXiCPlus = 4232, kXiCCPlusPlus = 4422, @@ -34,6 +35,7 @@ enum Code { kJpsi = 443, kChic1 = 20443, kBPlus = 521, + kBs = 531, kX3872 = 9920443 }; } // namespace pdg @@ -593,6 +595,62 @@ static const std::vector pTBinLabels = { static const std::vector cutVarLabels = {"m", "CPA", "Chi2PCA", "d0 Lc+", "d0 Pi", "pT Lc+", "pT Pi", "Lb decLen", "Lb decLenXY", "Imp. Par. Product", "DeltaMLc", "Cos ThetaStar"}; } // namespace hf_cuts_lb_tolcpi +namespace hf_cuts_bs_todspi +{ +static constexpr int npTBins = 12; +static constexpr int nCutVars = 12; +// default values for the pT bin edges (can be used to configure histogram axis) +// offset by 1 from the bin numbers in cuts array +constexpr double pTBins[npTBins + 1] = { + 0, + 0.5, + 1.0, + 2.0, + 3.0, + 4.0, + 5.0, + 7.0, + 10.0, + 13.0, + 16.0, + 20.0, + 24.0}; + +auto pTBins_v = std::vector{pTBins, pTBins + npTBins + 1}; + +// default values for the cuts +// DeltaM CPA chi2PCA d0Ds d0Pi pTDs pTPi BsDecayLength BsDecayLengthXY IPProd DeltaMDs CthetaStr +constexpr double cuts[npTBins][nCutVars] = {{1., 0.8, 1., 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.8}, /* 0 < pt < 0.5 */ + {1., 0.8, 1., 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.8}, /* 0.5 < pt < 1 */ + {1., 0.8, 1., 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.8}, /* 1 < pt < 2 */ + {1., 0.8, 1., 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.8}, /* 2 < pt < 3 */ + {1., 0.8, 1., 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.8}, /* 3 < pt < 4 */ + {1., 0.8, 1., 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.8}, /* 4 < pt < 5 */ + {1., 0.8, 1., 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.8}, /* 5 < pt < 7 */ + {1., 0.8, 1., 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.8}, /* 7 < pt < 10 */ + {1., 0.8, 1., 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.8}, /* 10 < pt < 13 */ + {1., 0.8, 1., 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.8}, /* 13 < pt < 16 */ + {1., 0.8, 1., 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.8}, /* 16 < pt < 20 */ + {1., 0.8, 1., 0.01, 0.01, 1.0, 0.15, 0.05, 0.05, 0., 0.1, 0.8}}; /* 20 < pt < 24 */ +// row labels +static const std::vector pTBinLabels = { + "pT bin 0", + "pT bin 1", + "pT bin 2", + "pT bin 3", + "pT bin 4", + "pT bin 5", + "pT bin 6", + "pT bin 7", + "pT bin 8", + "pT bin 9", + "pT bin 10", + "pT bin 11"}; + +// column labels +static const std::vector cutVarLabels = {"m", "CPA", "Chi2PCA", "d0 Ds+", "d0 Pi", "pT Ds+", "pT Pi", "Bs decLen", "Bs decLenXY", "Imp. Par. Product", "DeltaMDs", "Cos ThetaStar"}; +} // namespace hf_cuts_bs_todspi + namespace hf_cuts_x_tojpsipipi { static constexpr int npTBins = 9; diff --git a/PWGHF/DataModel/HFCandidateSelectionTables.h b/PWGHF/DataModel/HFCandidateSelectionTables.h index dcb9cc372e5..31e9cc05ac2 100644 --- a/PWGHF/DataModel/HFCandidateSelectionTables.h +++ b/PWGHF/DataModel/HFCandidateSelectionTables.h @@ -93,6 +93,14 @@ DECLARE_SOA_COLUMN(IsSelDplusToPiKPi, isSelDplusToPiKPi, int); //! DECLARE_SOA_TABLE(HFSelDplusToPiKPiCandidate, "AOD", "HFSELDPLUSCAND", //! hf_selcandidate_dplus::IsSelDplusToPiKPi); +namespace hf_selcandidate_ds +{ +DECLARE_SOA_COLUMN(IsSelDsKKpi, isSelDsKKpi, int); //! +DECLARE_SOA_COLUMN(IsSelDspiKK, isSelDspiKK, int); //! +} // namespace hf_selcandidate_ds +DECLARE_SOA_TABLE(HFSelDsCandidate, "AOD", "HFSELLCCAND", //! + hf_selcandidate_ds::IsSelDsKKpi, hf_selcandidate_ds::IsSelDspiKK); + namespace hf_selcandidate_lc { DECLARE_SOA_COLUMN(IsSelLcpKpi, isSelLcpKpi, int); //! @@ -183,6 +191,13 @@ DECLARE_SOA_COLUMN(IsSelLbToLcPi, isSelLbToLcPi, int); //! DECLARE_SOA_TABLE(HFSelLbToLcPiCandidate, "AOD", "HFSELLBCAND", //! hf_selcandidate_lb::IsSelLbToLcPi); +namespace hf_selcandidate_bs +{ +DECLARE_SOA_COLUMN(IsSelBsToDsPi, isSelBsToDsPi, int); //! +} // namespace hf_selcandidate_bs +DECLARE_SOA_TABLE(HFSelBsToDsPiCandidate, "AOD", "HFSELBSCAND", //! + hf_selcandidate_bs::IsSelBsToDsPi); + namespace hf_selcandidate_x { DECLARE_SOA_COLUMN(IsSelXToJpsiToEEPiPi, isSelXToJpsiToEEPiPi, int); //! diff --git a/PWGHF/DataModel/HFSecondaryVertex.h b/PWGHF/DataModel/HFSecondaryVertex.h index 8b51151a7b6..92d0b337e8a 100644 --- a/PWGHF/DataModel/HFSecondaryVertex.h +++ b/PWGHF/DataModel/HFSecondaryVertex.h @@ -643,6 +643,18 @@ auto InvMassDPlus(const T& candidate) return candidate.m(array{RecoDecay::getMassPDG(kPiPlus), RecoDecay::getMassPDG(kKPlus), RecoDecay::getMassPDG(kPiPlus)}); } +template +auto InvMassDsKKpi(const T& candidate) +{ + return candidate.m(array{RecoDecay::getMassPDG(kKPlus), RecoDecay::getMassPDG(kKMinus), RecoDecay::getMassPDG(kPiPlus)}); +} + +template +auto InvMassDspiKK(const T& candidate) +{ + return candidate.m(array{RecoDecay::getMassPDG(kPiPlus), RecoDecay::getMassPDG(kKMinus), RecoDecay::getMassPDG(kKPlus)}); +} + // Λc± → p± K∓ π± template @@ -1205,6 +1217,95 @@ DECLARE_SOA_TABLE(HfCandLbMCRec, "AOD", "HFCANDLBMCREC", //! DECLARE_SOA_TABLE(HfCandLbMCGen, "AOD", "HFCANDLBMCGEN", //! hf_cand_lb::FlagMCMatchGen, hf_cand_lb::OriginMCGen); + +// specific Bs candidate properties +namespace hf_cand_bs +{ +DECLARE_SOA_INDEX_COLUMN_FULL(Index0, index0, int, HfCandProng3, "_0"); // Bs index +// MC matching result: +DECLARE_SOA_COLUMN(FlagMCMatchRec, flagMCMatchRec, int8_t); // reconstruction level +DECLARE_SOA_COLUMN(FlagMCMatchGen, flagMCMatchGen, int8_t); // generator level +DECLARE_SOA_COLUMN(OriginMCRec, originMCRec, int8_t); // particle origin, reconstruction level +DECLARE_SOA_COLUMN(OriginMCGen, originMCGen, int8_t); // particle origin, generator level +DECLARE_SOA_COLUMN(DebugMCRec, debugMCRec, int8_t); // debug flag for mis-association reconstruction level +// mapping of decay types +enum DecayType { BsToDsPi }; // move this to a dedicated cascade namespace in the future? + +// Bs → Ds+ π- → K+ K- π+ π- +// float massBs = RecoDecay::getMassPDG(pdg::Code::kBs); +template +auto CtBs(const T& candidate) +{ + return candidate.ct(RecoDecay::getMassPDG(pdg::Code::kBs)); +} + +template +auto YBs(const T& candidate) +{ + return candidate.y(RecoDecay::getMassPDG(pdg::Code::kBs)); +} + +template +auto EBs(const T& candidate) +{ + return candidate.e(RecoDecay::getMassPDG(pdg::Code::kBs)); +} +template +auto InvMassBsToDsPi(const T& candidate) +{ + return candidate.m(array{RecoDecay::getMassPDG(pdg::Code::kDs), RecoDecay::getMassPDG(kPiMinus)}); +} +} // namespace hf_cand_bs + +// declare dedicated Bs candidate table +DECLARE_SOA_TABLE(HfCandBsBase, "AOD", "HFCANDBSBASE", + // general columns + HFCAND_COLUMNS, + // 3-prong specific columns + hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, + hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1, + hf_cand::ImpactParameter0, hf_cand::ImpactParameter1, + hf_cand::ErrorImpactParameter0, hf_cand::ErrorImpactParameter1, + hf_cand_bs::Index0Id, hf_track_index::Index1Id, + hf_track_index::HFflag, + /* dynamic columns */ + hf_cand_prong2::M, + hf_cand_prong2::M2, + hf_cand_prong2::ImpactParameterProduct, + /* dynamic columns that use candidate momentum components */ + hf_cand::Pt, + hf_cand::Pt2, + hf_cand::P, + hf_cand::P2, + hf_cand::PVector, + hf_cand::CPA, + hf_cand::CPAXY, + hf_cand::Ct, + hf_cand::ImpactParameterXY, + hf_cand_prong2::MaxNormalisedDeltaIP, + hf_cand::Eta, + hf_cand::Phi, + hf_cand::Y, + hf_cand::E, + hf_cand::E2); + +// extended table with expression columns that can be used as arguments of dynamic columns +DECLARE_SOA_EXTENDED_TABLE_USER(HfCandBsExt, HfCandBsBase, "HFCANDBSEXT", + hf_cand_prong2::Px, hf_cand_prong2::Py, hf_cand_prong2::Pz); + +using HfCandBs = HfCandBsExt; + +// table with results of reconstruction level MC matching +DECLARE_SOA_TABLE(HfCandBsMCRec, "AOD", "HFCANDBSMCREC", //! + hf_cand_bs::FlagMCMatchRec, + hf_cand_bs::OriginMCRec, + hf_cand_bs::DebugMCRec); + +// table with results of generator level MC matching +DECLARE_SOA_TABLE(HfCandBsMCGen, "AOD", "HFCANDBSMCGEN", //! + hf_cand_bs::FlagMCMatchGen, + hf_cand_bs::OriginMCGen); + } // namespace o2::aod #endif // O2_ANALYSIS_HFSECONDARYVERTEX_H_ diff --git a/PWGHF/TableProducer/CMakeLists.txt b/PWGHF/TableProducer/CMakeLists.txt index b97060d5a3f..68cfd6a7ee7 100644 --- a/PWGHF/TableProducer/CMakeLists.txt +++ b/PWGHF/TableProducer/CMakeLists.txt @@ -49,6 +49,11 @@ o2physics_add_dpl_workflow(candidate-creator-lb PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsVertexing ROOT::EG COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(candidate-creator-bs + SOURCES HFCandidateCreatorBs.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsVertexing ROOT::EG + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(candidate-creator-x SOURCES HFCandidateCreatorX.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsVertexing ROOT::EG @@ -139,6 +144,11 @@ o2physics_add_dpl_workflow(bplus-tod0pi-candidate-selector PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsVertexing COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(bs-todspi-candidate-selector + SOURCES HFBsToDsPiCandidateSelector.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsVertexing + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(lb-tolcpi-candidate-selector SOURCES HFLbToLcPiCandidateSelector.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsVertexing diff --git a/PWGHF/TableProducer/HFBsToDsPiCandidateSelector.cxx b/PWGHF/TableProducer/HFBsToDsPiCandidateSelector.cxx new file mode 100644 index 00000000000..643f6ba52e7 --- /dev/null +++ b/PWGHF/TableProducer/HFBsToDsPiCandidateSelector.cxx @@ -0,0 +1,160 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file HFBsToDsPiCandidateSelector.cxx +/// \brief Bs → Ds+ π- candidate selector +/// +/// \author Panos Christakoglou , Nikhef + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "PWGHF/Core/HFSelectorCuts.h" +#include "PWGHF/DataModel/HFSecondaryVertex.h" +#include "PWGHF/DataModel/HFCandidateSelectionTables.h" +#include "PWGHF/Core/HFSelectorCuts.h" + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::aod::hf_cand_bs; +using namespace o2::analysis; +using namespace o2::aod::hf_cand_prong2; +using namespace o2::analysis::hf_cuts_bs_todspi; + +struct HfBsToDsPiCandidateSelector { + Produces hfSelBsToDsPiCandidate; + + Configurable pTCandMin{"pTCandMin", 0., "Lower bound of candidate pT"}; + Configurable pTCandMax{"pTCandMax", 50., "Upper bound of candidate pT"}; + + // Track quality + Configurable TPCNClsFindablePIDCut{"TPCNClsFindablePIDCut", 70., "Lower bound of TPC findable clusters for good PID"}; + + // TPC PID + Configurable pidTPCMinpT{"pidTPCMinpT", 0.15, "Lower bound of track pT for TPC PID"}; + Configurable pidTPCMaxpT{"pidTPCMaxpT", 10., "Upper bound of track pT for TPC PID"}; + Configurable nSigmaTPC{"nSigmaTPC", 5., "Nsigma cut on TPC only"}; + Configurable nSigmaTPCCombined{"nSigmaTPCCombined", 5., "Nsigma cut on TPC combined with TOF"}; + + // TOF PID + Configurable pidTOFMinpT{"pidTOFMinpT", 0.15, "Lower bound of track pT for TOF PID"}; + Configurable pidTOFMaxpT{"pidTOFMaxpT", 10., "Upper bound of track pT for TOF PID"}; + Configurable nSigmaTOF{"nSigmaTOF", 5., "Nsigma cut on TOF only"}; + Configurable nSigmaTOFCombined{"nSigmaTOFCombined", 5., "Nsigma cut on TOF combined with TPC"}; + + Configurable> pTBins{"pTBins", std::vector{hf_cuts_bs_todspi::pTBins_v}, "pT bin limits"}; + Configurable> cuts{"Bs_to_dspi_cuts", {hf_cuts_bs_todspi::cuts[0], npTBins, nCutVars, pTBinLabels, cutVarLabels}, "Bs0 candidate selection per pT bin"}; + Configurable selectionFlagDs{"selectionFlagDs", 1, "Selection Flag for Ds+"}; + + // Apply topological cuts as defined in HFSelectorCuts.h; return true if candidate passes all cuts + template + bool selectionTopol(const T1& hfCandBs, const T2& hfCandDs, const T3& trackPi) + { + auto candpT = hfCandBs.pt(); + int pTBin = findBin(pTBins, candpT); + if (pTBin == -1) { + // Printf("Bs topol selection failed at getpTBin"); + return false; + } + + // Bs mass cut + if (std::abs(InvMassBsToDsPi(hfCandBs) - RecoDecay::getMassPDG(pdg::Code::kBs)) > cuts->get(pTBin, "m")) { + // Printf("Bs topol selection failed at mass diff check"); + return false; + } + + // pion pt + if (trackPi.pt() < cuts->get(pTBin, "pT Pi")) { + return false; + } + + // Ds+ pt + if (hfCandDs.pt() < cuts->get(pTBin, "pT Ds+")) { + return false; + } + + // Ds mass + // if (trackPi.sign() < 0) { + // if (std::abs(InvMassDspKpi(hfCandDs) - RecoDecay::getMassPDG(pdg::Code::kDs)) > cuts->get(pTBin, "DeltaMDs")) { + // return false; + // } + // } + + // Bs Decay length + if (hfCandBs.decayLength() < cuts->get(pTBin, "Bs decLen")) { + return false; + } + + // Bs Decay length XY + if (hfCandBs.decayLengthXY() < cuts->get(pTBin, "Bs decLenXY")) { + return false; + } + + // Bs chi2PCA cut + if (hfCandBs.chi2PCA() > cuts->get(pTBin, "Chi2PCA")) { + // Printf("Bs selection failed at chi2PCA"); + return false; + } + + // Bs CPA cut + if (hfCandBs.cpa() < cuts->get(pTBin, "CPA")) { + return false; + } + + // d0 of pi + if (std::abs(hfCandBs.impactParameter1()) < cuts->get(pTBin, "d0 Pi")) { + return false; + } + + // d0 of Ds+ + if (std::abs(hfCandBs.impactParameter0()) < cuts->get(pTBin, "d0 Ds+")) { + return false; + } + + return true; + } + + void process(aod::HfCandBs const& hfCandBss, soa::Join, aod::BigTracksPID const&) + { + for (auto& hfCandBs : hfCandBss) { // looping over Bs candidates + + int statusBs = 0; + + // check if flagged as Bs --> Ds+ π- + if (!(hfCandBs.hfflag() & 1 << hf_cand_bs::DecayType::BsToDsPi)) { + hfSelBsToDsPiCandidate(statusBs); + // Printf("Bs candidate selection failed at hfflag check"); + continue; + } + + // Ds is always index0 and pi is index1 by default + auto candDs = hfCandBs.index0_as>(); + auto trackPi = hfCandBs.index1_as(); + + // topological cuts + if (!selectionTopol(hfCandBs, candDs, trackPi)) { + hfSelBsToDsPiCandidate(statusBs); + // Printf("Bs candidate selection failed at selection topology"); + continue; + } + + hfSelBsToDsPiCandidate(1); + // Printf("Bs candidate selection successful, candidate should be selected"); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow{}; + workflow.push_back(adaptAnalysisTask(cfgc)); + return workflow; +} diff --git a/PWGHF/TableProducer/HFCandidateCreatorBs.cxx b/PWGHF/TableProducer/HFCandidateCreatorBs.cxx new file mode 100644 index 00000000000..358c1b39e47 --- /dev/null +++ b/PWGHF/TableProducer/HFCandidateCreatorBs.cxx @@ -0,0 +1,296 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file HFCandidateCreatorBs.cxx +/// \brief Reconstruction of Bs candidates +/// \note Adapted from HFCandidateCreatorXicc +/// +/// \author Panos Christakoglou , Nikhef + +#include "Framework/AnalysisTask.h" +#include "DetectorsVertexing/DCAFitterN.h" +#include "PWGHF/DataModel/HFSecondaryVertex.h" +#include "Common/Core/trackUtilities.h" +#include "ReconstructionDataFormats/DCA.h" +#include "ReconstructionDataFormats/V0.h" +#include "PWGHF/DataModel/HFCandidateSelectionTables.h" + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::aod::hf_cand; +using namespace o2::aod::hf_cand_prong2; +using namespace o2::aod::hf_cand_prong3; +using namespace o2::aod::hf_cand_bs; +using namespace o2::framework::expressions; + +void customize(std::vector& workflowOptions) +{ + ConfigParamSpec optionDoMC{"doMC", VariantType::Bool, true, {"Perform MC matching."}}; + workflowOptions.push_back(optionDoMC); +} + +#include "Framework/runDataProcessing.h" + +/// Reconstruction of Λb candidates +struct HFCandidateCreatorBs { + Produces rowCandidateBase; + + Configurable magneticField{"magneticField", 20., "magnetic field"}; + Configurable b_propdca{"b_propdca", true, "create tracks version propagated to PCA"}; + Configurable d_maxr{"d_maxr", 200., "reject PCA's above this radius"}; + Configurable d_maxdzini{"d_maxdzini", 4., "reject (if>0) PCA candidate if tracks DZ exceeds threshold"}; + Configurable d_minparamchange{"d_minparamchange", 1.e-3, "stop iterations if largest change of any Bs is smaller than this"}; + Configurable d_minrelchi2change{"d_minrelchi2change", 0.9, "stop iterations is chi2/chi2old > this"}; + Configurable ptPionMin{"ptPionMin", 0.5, "minimum pion pT threshold (GeV/c)"}; + + OutputObj hMassDsToKKPi{TH1F("hMassDsToKKPi", "D_{s}^{#plus} candidates;inv. mass (K^{#plus}K^{#minus} #pi^{#plus}) (GeV/#it{c}^{2});entries", 500, 0., 5.)}; + OutputObj hPtDs{TH1F("hPtDs", "D_{s}^{#plus} candidates;D_{s}^{#plus} candidate #it{p}_{T} (GeV/#it{c});entries", 100, 0., 10.)}; + OutputObj hPtPion{TH1F("hPtPion", "#pi^{#minus} candidates;#pi^{#minus} candidate #it{p}_{T} (GeV/#it{c});entries", 100, 0., 10.)}; + OutputObj hCPADs{TH1F("hCPADs", "D_{s}^{#plus} candidates;D_{s}^{#plus} cosine of pointing angle;entries", 110, -1.1, 1.1)}; + OutputObj hMassBsToDsPi{TH1F("hMassBsToDsPi", "2-prong candidates;inv. mass (B_{s}^{0} #rightarrow D_{s}^{#plus}#pi^{#minus} #rightarrow K^{#plus}K^{#minus}#pi^{#plus}#pi^{#minus}) (GeV/#it{c}^{2});entries", 500, 3., 8.)}; + OutputObj hCovPVXX{TH1F("hCovPVXX", "2-prong candidates;XX element of cov. matrix of prim. vtx position (cm^{2});entries", 100, 0., 1.e-4)}; + OutputObj hCovSVXX{TH1F("hCovSVXX", "2-prong candidates;XX element of cov. matrix of sec. vtx position (cm^{2});entries", 100, 0., 0.2)}; + + double massPi = RecoDecay::getMassPDG(kPiMinus); + double massDs = RecoDecay::getMassPDG(pdg::Code::kDs); + double massDsPi = 0.; + + Configurable d_selectionFlagDs{"d_selectionFlagDs", 1, "Selection Flag for Ds"}; + Configurable cutYCandMax{"cutYCandMax", -1., "max. cand. rapidity"}; + Filter filterSelectCandidates = (aod::hf_selcandidate_ds::isSelDsKKpi >= d_selectionFlagDs || aod::hf_selcandidate_ds::isSelDspiKK >= d_selectionFlagDs); + + void process(aod::Collision const& collision, + soa::Filtered> const& dsCands, + aod::BigTracks const& tracks) + { + // 2-prong vertex fitter + o2::vertexing::DCAFitterN<2> df2; + df2.setBz(magneticField); + df2.setPropagateToPCA(b_propdca); + df2.setMaxR(d_maxr); + df2.setMaxDZIni(d_maxdzini); + df2.setMinParamChange(d_minparamchange); + df2.setMinRelChi2Change(d_minrelchi2change); + df2.setUseAbsDCA(true); + + // 3-prong vertex fitter (to rebuild Ds vertex) + o2::vertexing::DCAFitterN<3> df3; + df3.setBz(magneticField); + df3.setPropagateToPCA(b_propdca); + df3.setMaxR(d_maxr); + df3.setMaxDZIni(d_maxdzini); + df3.setMinParamChange(d_minparamchange); + df3.setMinRelChi2Change(d_minrelchi2change); + df3.setUseAbsDCA(true); + + // loop over Ds candidates + for (auto& dsCand : dsCands) { + if (!(dsCand.hfflag() & 1 << o2::aod::hf_cand_prong3::DecayType::DsToPiKK)) { + continue; + } + if (dsCand.isSelDsKKpi() >= d_selectionFlagDs) { + hMassDsToKKPi->Fill(InvMassDsKKpi(dsCand), dsCand.pt()); + } + if (dsCand.isSelDspiKK() >= d_selectionFlagDs) { + hMassDsToKKPi->Fill(InvMassDspiKK(dsCand), dsCand.pt()); + } + hPtDs->Fill(dsCand.pt()); + hCPADs->Fill(dsCand.cpa()); + + auto track0 = dsCand.index0_as(); + auto track1 = dsCand.index1_as(); + auto track2 = dsCand.index2_as(); + auto trackParVar0 = getTrackParCov(track0); + auto trackParVar1 = getTrackParCov(track1); + auto trackParVar2 = getTrackParCov(track2); + auto collision = track0.collision(); + + // reconstruct the 3-prong secondary vertex + if (df3.process(trackParVar0, trackParVar1, trackParVar2) == 0) { + continue; + } + const auto& secondaryVertex = df3.getPCACandidate(); + trackParVar0.propagateTo(secondaryVertex[0], magneticField); + trackParVar1.propagateTo(secondaryVertex[0], magneticField); + trackParVar2.propagateTo(secondaryVertex[0], magneticField); + + array pvecpK = {track0.px() + track1.px(), track0.py() + track1.py(), track0.pz() + track1.pz()}; + array pvecDs = {pvecpK[0] + track2.px(), pvecpK[1] + track2.py(), pvecpK[2] + track2.pz()}; + auto trackpK = o2::dataformats::V0(df3.getPCACandidatePos(), pvecpK, df3.calcPCACovMatrixFlat(), + trackParVar0, trackParVar1, {0, 0}, {0, 0}); + auto trackDs = o2::dataformats::V0(df3.getPCACandidatePos(), pvecDs, df3.calcPCACovMatrixFlat(), + trackpK, trackParVar2, {0, 0}, {0, 0}); + + int index0Ds = track0.globalIndex(); + int index1Ds = track1.globalIndex(); + int index2Ds = track2.globalIndex(); + // int charge = track0.sign() + track1.sign() + track2.sign(); + + for (auto& trackPion : tracks) { + if (trackPion.pt() < ptPionMin) { + continue; + } + if (trackPion.sign() > 0) { + continue; + } + if (trackPion.globalIndex() == index0Ds || trackPion.globalIndex() == index1Ds || trackPion.globalIndex() == index2Ds) { + continue; + } + hPtPion->Fill(trackPion.pt()); + array pvecPion; + auto trackParVarPi = getTrackParCov(trackPion); + + // reconstruct the 3-prong Ds vertex + if (df2.process(trackDs, trackParVarPi) == 0) { + continue; + } + + // calculate relevant properties + const auto& secondaryVertexBs = df2.getPCACandidate(); + auto chi2PCA = df2.getChi2AtPCACandidate(); + auto covMatrixPCA = df2.calcPCACovMatrixFlat(); + + df2.propagateTracksToVertex(); + df2.getTrack(0).getPxPyPzGlo(pvecDs); + df2.getTrack(1).getPxPyPzGlo(pvecPion); + + auto primaryVertex = getPrimaryVertex(collision); + auto covMatrixPV = primaryVertex.getCov(); + o2::dataformats::DCA impactParameter0; + o2::dataformats::DCA impactParameter1; + trackDs.propagateToDCA(primaryVertex, magneticField, &impactParameter0); + trackParVarPi.propagateToDCA(primaryVertex, magneticField, &impactParameter1); + + hCovSVXX->Fill(covMatrixPCA[0]); + hCovPVXX->Fill(covMatrixPV[0]); + + // get uncertainty of the decay length + double phi, theta; + getPointDirection(array{collision.posX(), collision.posY(), collision.posZ()}, secondaryVertexBs, phi, theta); + auto errorDecayLength = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta)); + auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.)); + + int hfFlag = 1 << hf_cand_bs::DecayType::BsToDsPi; + + // fill the candidate table for the Bs here: + rowCandidateBase(collision.globalIndex(), + collision.posX(), collision.posY(), collision.posZ(), + secondaryVertexBs[0], secondaryVertexBs[1], secondaryVertexBs[2], + errorDecayLength, errorDecayLengthXY, + chi2PCA, + pvecDs[0], pvecDs[1], pvecDs[2], + pvecPion[0], pvecPion[1], pvecPion[2], + impactParameter0.getY(), impactParameter1.getY(), + std::sqrt(impactParameter0.getSigmaY2()), std::sqrt(impactParameter1.getSigmaY2()), + dsCand.globalIndex(), trackPion.globalIndex(), + hfFlag); + + // calculate invariant mass + auto arrayMomenta = array{pvecDs, pvecPion}; + massDsPi = RecoDecay::M(std::move(arrayMomenta), array{massDs, massPi}); + if (dsCand.isSelDsKKpi() > 0) { + hMassBsToDsPi->Fill(massDsPi); + } + if (dsCand.isSelDspiKK() > 0) { + hMassBsToDsPi->Fill(massDsPi); + } + } // pi- loop + } // Ds loop + } // process +}; // struct + +/// Extends the base table with expression columns. +struct HFCandidateCreatorBsExpressions { + Spawns rowCandidateBs; + void init(InitContext const&) {} +}; + +/// Performs MC matching. +struct HFCandidateCreatorBsMC { + Produces rowMCMatchRec; + Produces rowMCMatchGen; + + void process(aod::HfCandBs const& candidates, + aod::HfCandProng3, + aod::BigTracksMC const& tracks, + aod::McParticles const& particlesMC) + { + int indexRec = -1; + int8_t sign = 0; + int8_t flag = 0; + int8_t origin = 0; + int8_t debug = 0; + + // Match reconstructed candidates. + for (auto& candidate : candidates) { + // Printf("New rec. candidate"); + flag = 0; + origin = 0; + debug = 0; + auto dsCand = candidate.index0(); + auto arrayDaughters = array{dsCand.index0_as(), + dsCand.index1_as(), + dsCand.index2_as(), + candidate.index1_as()}; + auto arrayDaughtersDs = array{dsCand.index0_as(), + dsCand.index1_as(), + dsCand.index2_as()}; + // Bs → Ds+ π- + // Printf("Checking Bs → Ds+ π-"); + indexRec = RecoDecay::getMatchedMCRec(particlesMC, arrayDaughters, pdg::Code::kBs, array{+kKPlus, -kKPlus, +kPiPlus, -kPiPlus}, true, &sign, 2); + if (indexRec > -1) { + // Bs → Ds+ π- + // Printf("Checking Bs → Ds+ π-"); + indexRec = RecoDecay::getMatchedMCRec(particlesMC, arrayDaughtersDs, pdg::Code::kDs, array{+kKPlus, -kKPlus, +kPiPlus}, true, &sign, 1); + if (indexRec > -1) { + flag = 1 << hf_cand_bs::DecayType::BsToDsPi; + } else { + debug = 1; + LOGF(info, "WARNING: Λb in decays in the expected final state but the condition on the intermediate state is not fulfilled"); + } + } + rowMCMatchRec(flag, origin, debug); + } + + // Match generated particles. + for (auto& particle : particlesMC) { + // Printf("New gen. candidate"); + flag = 0; + origin = 0; + // Bs → Ds+ π- + if (RecoDecay::isMatchedMCGen(particlesMC, particle, pdg::Code::kBs, array{int(pdg::Code::kDs), -kPiPlus}, true)) { + // Match Ds+ -> φπ -> K+K-π + auto DsCandMC = particlesMC.iteratorAt(particle.daughter0Id()); + + // Printf("Checking Ds+ -> φπ -> K+K-π"); + if (RecoDecay::isMatchedMCGen(particlesMC, DsCandMC, int(pdg::Code::kDs), array{+kKPlus, -kKPlus, +kPiPlus}, true, &sign)) { + flag = sign * (1 << hf_cand_bs::DecayType::BsToDsPi); + } + } + rowMCMatchGen(flag, origin); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow{ + adaptAnalysisTask(cfgc, TaskName{"hf-cand-creator-bs"}), + adaptAnalysisTask(cfgc, TaskName{"hf-cand-creator-bs-expressions"})}; + const bool doMC = cfgc.options().get("doMC"); + if (doMC) { + workflow.push_back(adaptAnalysisTask(cfgc, TaskName{"hf-cand-creator-bs-mc"})); + } + return workflow; +} diff --git a/PWGHF/Tasks/CMakeLists.txt b/PWGHF/Tasks/CMakeLists.txt index cba2510e31d..a750717732b 100644 --- a/PWGHF/Tasks/CMakeLists.txt +++ b/PWGHF/Tasks/CMakeLists.txt @@ -70,6 +70,11 @@ o2physics_add_dpl_workflow(task-bplus PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsVertexing COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(task-bs + SOURCES taskBs.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsVertexing + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(task-lb SOURCES taskLb.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsVertexing diff --git a/PWGHF/Tasks/taskBs.cxx b/PWGHF/Tasks/taskBs.cxx new file mode 100644 index 00000000000..295a19605d8 --- /dev/null +++ b/PWGHF/Tasks/taskBs.cxx @@ -0,0 +1,281 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file taskBs.cxx +/// \brief Bs analysis task +/// +/// \author Panos Christakoglou , Nikhef + +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "PWGHF/DataModel/HFSecondaryVertex.h" +#include "PWGHF/Core/HFSelectorCuts.h" +#include "PWGHF/DataModel/HFCandidateSelectionTables.h" +#include "Common/DataModel/Centrality.h" + +using namespace o2; +using namespace o2::aod; +using namespace o2::analysis; +using namespace o2::framework; +using namespace o2::aod::hf_cand_prong2; +using namespace o2::aod::hf_cand_prong3; +using namespace o2::aod::hf_cand_bs; +using namespace o2::analysis::hf_cuts_bs_todspi; +using namespace o2::framework::expressions; + +void customize(std::vector& workflowOptions) +{ + ConfigParamSpec optionDoMC{"doMC", VariantType::Bool, true, {"Fill MC histograms."}}; + workflowOptions.push_back(optionDoMC); +} + +#include "Framework/runDataProcessing.h" + +/// Λb0 analysis task +struct HfTaskBs { + HistogramRegistry registry{ + "registry", + {{"hPtProng0", "Bs candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{1000, 0., 50.}}}}, + {"hPtProng1", "Bs candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 10.}}}}, + {"hPtCand", "Bs candidates;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{1000, 0., 50.}}}}, + {"hCentrality", "centrality;centrality percentile;entries", {HistType::kTH1F, {{100, 0., 100.}}}}}}; + + Configurable selectionFlagBs{"selectionFlagBs", 1, "Selection Flag for Bs"}; + Configurable cutYCandMax{"cutYCandMax", 1.44, "max. cand. rapidity"}; + Configurable> bins{"pTBins", std::vector{hf_cuts_bs_todspi::pTBins_v}, "pT bin limits"}; + + void init(o2::framework::InitContext&) + { + auto vbins = (std::vector)bins; + registry.add("hMass", "B_{s} candidates;inv. mass D_{s}^{#plus}#pi^{#minus} (GeV/#it{c}^{2});#it{p}_{T} (GeV/#it{c}); centrality", {HistType::kTH3F, {{500, 0., 10.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}, {100, 0., 100.}}}); + registry.add("hDecLength", "B_{s} candidates;decay length (cm);entries", {HistType::kTH2F, {{200, 0., 0.4}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hDecLengthXY", "B_{s} candidates;decay length xy (cm);entries", {HistType::kTH2F, {{200, 0., 0.4}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0Prong0", "B_{s} candidates;prong 0 (D_{s}^{#plus}) DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{100, -0.05, 0.05}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0Prong1", "B_{s} candidates;prong 1 (#pi^{#minus}) DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{100, -0.05, 0.05}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hCPA", "B_{s} candidates;B_{s} candidate cosine of pointing angle;entries", {HistType::kTH2F, {{110, -1.1, 1.1}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hEta", "B_{s} candidates;B_{s} candidate #it{#eta};entries", {HistType::kTH2F, {{100, -2., 2.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hRapidity", "B_{s} candidates;B_{s} candidate #it{y};entries", {HistType::kTH2F, {{100, -2., 2.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hImpParErr", "B_{s} candidates;B_{s} candidate impact parameter error (cm);entries", {HistType::kTH2F, {{100, -1., 1.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hDecLenErr", "B_{s} candidates;B_{s} candidate decay length error (cm);entries", {HistType::kTH2F, {{100, 0., 1.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hDecLenXYErr", "B_{s} candidates;B_{s} candidate decay length xy error (cm);entries", {HistType::kTH2F, {{100, 0., 1.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hIPProd", "B_{s} candidates;B_{s} candidate impact parameter product;entries", {HistType::kTH2F, {{100, -0.5, 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hInvMassDs", "B_{s} candidates;prong0, D_{s}^{+} inv. mass (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{500, 0, 5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + Filter filterSelectCandidates = (aod::hf_selcandidate_bs::isSelBsToDsPi >= selectionFlagBs); + + void process(soa::Join::iterator const& collision, soa::Filtered> const& candidates, soa::Join, aod::BigTracks) + { + float centrality = collision.centV0M(); + registry.fill(HIST("hCentrality"), centrality); + + for (auto& candidate : candidates) { + if (!(candidate.hfflag() & 1 << hf_cand_bs::DecayType::BsToDsPi)) { + continue; + } + if (cutYCandMax >= 0. && std::abs(YBs(candidate)) > cutYCandMax) { + continue; + } + + auto candDs = candidate.index0_as>(); + auto candPi = candidate.index1_as(); + + registry.fill(HIST("hMass"), InvMassBsToDsPi(candidate), candidate.pt(), centrality); + registry.fill(HIST("hPtCand"), candidate.pt()); + registry.fill(HIST("hPtProng0"), candidate.ptProng0()); + registry.fill(HIST("hPtProng1"), candidate.ptProng1()); + registry.fill(HIST("hIPProd"), candidate.impactParameterProduct(), candidate.pt()); + registry.fill(HIST("hDecLength"), candidate.decayLength(), candidate.pt()); + registry.fill(HIST("hDecLengthXY"), candidate.decayLengthXY(), candidate.pt()); + registry.fill(HIST("hd0Prong0"), candidate.impactParameter0(), candidate.pt()); + registry.fill(HIST("hd0Prong1"), candidate.impactParameter1(), candidate.pt()); + registry.fill(HIST("hCPA"), candidate.cpa(), candidate.pt()); + registry.fill(HIST("hEta"), candidate.eta(), candidate.pt()); + registry.fill(HIST("hRapidity"), YBs(candidate), candidate.pt()); + registry.fill(HIST("hImpParErr"), candidate.errorImpactParameter0(), candidate.pt()); + registry.fill(HIST("hImpParErr"), candidate.errorImpactParameter1(), candidate.pt()); + registry.fill(HIST("hDecLenErr"), candidate.errorDecayLength(), candidate.pt()); + registry.fill(HIST("hDecLenXYErr"), candidate.errorDecayLengthXY(), candidate.pt()); + if (candPi.sign() < 0) { + registry.fill(HIST("hInvMassDs"), InvMassDsKKpi(candDs), candidate.pt()); + } + } // candidate loop + } // process +}; // struct + +/// Bs MC analysis and fill histograms +struct HfTaskBsMc { + HistogramRegistry registry{ + "registry", + {{"hPtRecSig", "Bs candidates (matched);candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{300, 0., 30.}}}}, + {"hPtRecBg", "Bs candidates (unmatched);candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{300, 0., 30.}}}}, + {"hPtGenSig", "Bs candidates (matched);candidate #it{p}_{T}^{gen.} (GeV/#it{c});entries", {HistType::kTH1F, {{300, 0., 10.}}}}, + {"hPtGen", "MC particles (matched);candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{300, 0., 30.}}}}}}; + + Configurable selectionFlagBs{"selectionFlagBs", 1, "Selection Flag for Bs"}; + Configurable cutYCandMax{"cutYCandMax", 0.8, "max. cand. rapidity"}; + Configurable> bins{"pTBins", std::vector{hf_cuts_bs_todspi::pTBins_v}, "pT bin limits"}; + + void init(o2::framework::InitContext&) + { + auto vbins = (std::vector)bins; + registry.add("hEtaGen", "MC particles (matched);B_{s} candidate #it{#eta}^{gen};entries", {HistType::kTH2F, {{100, -2., 2.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hYGen", "MC particles (matched);B_{s} candidate #it{y}^{gen};entries", {HistType::kTH2F, {{100, -2., 2.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hPtProng0Gen", "MC particles (matched);prong 0 (D_{s}^{+}) #it{p}_{T}^{gen} (GeV/#it{c});entries", {HistType::kTH2F, {{100, 0., 10.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hPtProng1Gen", "MC particles (matched);prong 1 (#pi^{-}) #it{p}_{T}^{gen} (GeV/#it{c});entries", {HistType::kTH2F, {{100, 0., 10.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hYProng0Gen", "MC particles (matched);prong 0 (D_{s}^{+}) #it{y}^{gen};entries", {HistType::kTH2F, {{100, -2, 2}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hYProng1Gen", "MC particles (matched);prong 1 (#pi^{-}) #it{y}^{gen};entries", {HistType::kTH2F, {{100, -2, 2}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hEtaProng0Gen", "MC particles (matched);prong 0 (B_{s}) #it{#eta}^{gen};entries", {HistType::kTH2F, {{100, -2, 2}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hEtaProng1Gen", "MC particles (matched);prong 1 (#pi^{-}) #it{#eta}^{gen};entries", {HistType::kTH2F, {{100, -2, 2}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hCPARecSig", "B_{s} candidates (matched);B_{s} candidate cosine of pointing angle;entries", {HistType::kTH2F, {{220, 0., 1.1}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hCPARecBg", "B_{s} candidates (unmatched);B_{s} candidate cosine of pointing angle;entries", {HistType::kTH2F, {{220, 0., 1.1}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hCPAxyRecSig", "B_{s} candidates (matched);B_{s} candidate CPAxy;entries", {HistType::kTH2F, {{220, 0., 1.1}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hCPAxyRecBg", "B_{s} candidates (unmatched);B_{s} candidate CPAxy;entries", {HistType::kTH2F, {{220, 0., 1.1}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hCPADsRecSig", "B_{s} candidates (matched);prong 0 (D_{s}^{+}) cosine of pointing angle;entries", {HistType::kTH2F, {{220, 0., 1.1}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hCPADsRecBg", "B_{s} candidates (unmatched);prong 0 (D_{s}^{+}) cosine of pointing angle;entries", {HistType::kTH2F, {{220, 0., 1.1}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hEtaRecSig", "B_{s} candidates (matched);B_{s} candidate #it{#eta};entries", {HistType::kTH2F, {{100, -2., 2.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hEtaRecBg", "B_{s} candidates (unmatched);B_{s} candidate #it{#eta};entries", {HistType::kTH2F, {{100, -2., 2.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hRapidityRecSig", "B_{s} candidates (matched);B_{s} candidate #it{y};entries", {HistType::kTH2F, {{100, -2., 2.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hRapidityRecBg", "B_{s} candidates (unmatched);B_{s} candidate #it{#y};entries", {HistType::kTH2F, {{100, -2., 2.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + + registry.add("hPtProng0RecSig", "B_{s} candidates (matched);prong 0 (D_{s}^{+}) #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH2F, {{100, 0., 10.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hPtProng1RecSig", "B_{s} candidates (matched);prong 1 (#pi^{#minus}) #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH2F, {{100, 0., 10.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hPtProng0RecBg", "B_{s} candidates (unmatched);prong 0 (D_{s}^{+}) #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH2F, {{100, 0., 10.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hPtProng1RecBg", "B_{s} candidates (unmatched);prong 1 (#pi^{#minus}) #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH2F, {{100, 0., 10.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassRecSig", "B_{s} candidates (matched);inv. mass D_{s}^{+}#pi^{+} (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{300, 4.0, 7.00}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hMassRecBg", "B_{s} candidates (unmatched);inv. mass D_{s}^{+}#pi^{+} (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{300, 4.0, 7.0}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0Prong0RecSig", "B_{s} candidates (matched);prong 0 (D_{s}^{+}) DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{200, -0.05, 0.05}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0Prong1RecSig", "B_{s} candidates (matched);prong 1 (#pi^{#minus}) DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{200, -0.05, 0.05}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0Prong0RecBg", "B_{s} candidates (unmatched);prong 0 (D_{s}^{+}) DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{200, -0.05, 0.05}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hd0Prong1RecBg", "B_{s} candidates (unmatched);prong 1 (#pi^{#minus}) DCAxy to prim. vertex (cm);entries", {HistType::kTH2F, {{200, -0.05, 0.05}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hDecLengthRecSig", "B_{s} candidates (matched);B_{s} candidate decay length (cm);entries", {HistType::kTH2F, {{100, 0., 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hDecLengthXYRecSig", "B_{s} candidates (matched);B_{s} candidate decay length xy (cm);entries", {HistType::kTH2F, {{100, 0., 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hDecLengthRecBg", "B_{s} candidates (unmatched);B_{s} candidate decay length (cm);entries", {HistType::kTH2F, {{100, 0., 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hDecLengthXYRecBg", "B_{s} candidates (unmatched);B_{s} candidate decay length xy(cm);entries", {HistType::kTH2F, {{100, 0., 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hDecLengthDsRecSig", "B_{s} candidates (matched);B_{s} candidate decay length (cm);entries", {HistType::kTH2F, {{100, 0., 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hDecLengthDsRecBg", "B_{s} candidates (unmatched);B_{s} candidate decay length (cm);entries", {HistType::kTH2F, {{100, 0., 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hDecLengthNormRecSig", "B_{s} candidates (matched);B_{s} candidate decay length (cm);entries", {HistType::kTH2F, {{100, 0., 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hDecLengthNormRecBg", "B_{s} candidates (unmatched);B_{s} candidate decay length (cm);entries", {HistType::kTH2F, {{100, 0., 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hImpParProdBsRecSig", "B_{s} candidates (matched);B_{s} candidate impact parameter product ;entries", {HistType::kTH2F, {{100, -0.5, 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hImpParProdBsRecBg", "B_{s} candidates (unmatched);B_{s} candidate impact parameter product ;entries", {HistType::kTH2F, {{100, -0.5, 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + + registry.add("hChi2PCARecSig", "B_{s} candidates (matched);sum of distances of the secondary vertex to its prongs;entries", {HistType::kTH2F, {{240, -0.01, 0.1}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hChi2PCARecBg", "B_{s} candidates (unmatched);sum of distances of the secondary vertex to its prongs;entries", {HistType::kTH2F, {{240, -0.01, 0.1}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hThetaStarRecSig", "B_{s} candidates (matched);B_{s} #cos(#theta^{*});entries", {HistType::kTH2F, {{110, -1.1, 1.1}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + registry.add("hThetaStarRecBg", "B_{s} candidates (unmatched);B_{s} #cos(#theta^{*});entries", {HistType::kTH2F, {{110, -1.1, 1.1}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); + } + + Filter filterSelectCandidates = (aod::hf_selcandidate_bs::isSelBsToDsPi >= selectionFlagBs); + + void process(soa::Filtered> const& candidates, + soa::Join const& particlesMC, aod::BigTracksMC const& tracks, aod::HfCandProng3 const&) + { + // MC rec + for (auto& candidate : candidates) { + // Printf("(Panos) MC candidate: pT: %lf",candidate.pt()); + if (!(candidate.hfflag() & 1 << hf_cand_bs::DecayType::BsToDsPi)) { + continue; + } + if (cutYCandMax >= 0. && std::abs(YBs(candidate)) > cutYCandMax) { + continue; + } + auto candDs = candidate.index0_as(); + if (std::abs(candidate.flagMCMatchRec()) == 1 << hf_cand_bs::DecayType::BsToDsPi) { + + auto indexMother = RecoDecay::getMother(particlesMC, candidate.index1_as().mcParticle_as>(), pdg::Code::kBs, true); + auto particleMother = particlesMC.iteratorAt(indexMother); + registry.fill(HIST("hPtGenSig"), particleMother.pt()); + registry.fill(HIST("hPtRecSig"), candidate.pt()); + registry.fill(HIST("hCPARecSig"), candidate.cpa(), candidate.pt()); + registry.fill(HIST("hCPAxyRecSig"), candidate.cpa(), candidate.pt()); + registry.fill(HIST("hEtaRecSig"), candidate.eta(), candidate.pt()); + registry.fill(HIST("hRapidityRecSig"), YBs(candidate), candidate.pt()); + registry.fill(HIST("hDecLengthRecSig"), candidate.decayLength(), candidate.pt()); + registry.fill(HIST("hDecLengthXYRecSig"), candidate.decayLengthXY(), candidate.pt()); + registry.fill(HIST("hMassRecSig"), InvMassBsToDsPi(candidate), candidate.pt()); + registry.fill(HIST("hd0Prong0RecSig"), candidate.impactParameter0(), candidate.pt()); + registry.fill(HIST("hd0Prong1RecSig"), candidate.impactParameter1(), candidate.pt()); + registry.fill(HIST("hPtProng0RecSig"), candidate.ptProng0(), candidate.pt()); + registry.fill(HIST("hPtProng1RecSig"), candidate.ptProng1(), candidate.pt()); + registry.fill(HIST("hImpParProdBsRecSig"), candidate.impactParameterProduct(), candidate.pt()); + registry.fill(HIST("hDecLengthNormRecSig"), candidate.decayLengthXYNormalised(), candidate.pt()); + registry.fill(HIST("hCPADsRecSig"), candDs.cpa(), candidate.pt()); + registry.fill(HIST("hDecLengthDsRecSig"), candDs.decayLength(), candidate.pt()); + registry.fill(HIST("hChi2PCARecSig"), candidate.chi2PCA(), candidate.pt()); + // registry.fill(HIST("hThetaStarRecSig"), candidate.cosThetaStar(), candidate.pt()); + } else { + registry.fill(HIST("hPtRecBg"), candidate.pt()); + registry.fill(HIST("hCPARecBg"), candidate.cpa(), candidate.pt()); + registry.fill(HIST("hCPAxyRecBg"), candidate.cpa(), candidate.pt()); + registry.fill(HIST("hEtaRecBg"), candidate.eta(), candidate.pt()); + registry.fill(HIST("hRapidityRecBg"), YBs(candidate), candidate.pt()); + registry.fill(HIST("hDecLengthRecBg"), candidate.decayLength(), candidate.pt()); + registry.fill(HIST("hDecLengthXYRecBg"), candidate.decayLengthXY(), candidate.pt()); + registry.fill(HIST("hMassRecBg"), InvMassBsToDsPi(candidate), candidate.pt()); + registry.fill(HIST("hd0Prong0RecBg"), candidate.impactParameter0(), candidate.pt()); + registry.fill(HIST("hd0Prong1RecBg"), candidate.impactParameter1(), candidate.pt()); + registry.fill(HIST("hPtProng0RecBg"), candidate.ptProng0(), candidate.pt()); + registry.fill(HIST("hPtProng1RecBg"), candidate.ptProng1(), candidate.pt()); + registry.fill(HIST("hImpParProdBsRecBg"), candidate.impactParameterProduct(), candidate.pt()); + registry.fill(HIST("hDecLengthNormRecBg"), candidate.decayLengthXYNormalised(), candidate.pt()); + registry.fill(HIST("hCPADsRecBg"), candDs.cpa(), candidate.pt()); + registry.fill(HIST("hDecLengthDsRecBg"), candDs.decayLength(), candidate.pt()); + registry.fill(HIST("hChi2PCARecBg"), candidate.chi2PCA(), candidate.pt()); + // registry.fill(HIST("hThetaStarRecBg"), candidate.cosThetaStar(), candidate.pt()); + } + } // rec + + // MC gen. level + // Printf("MC Particles: %d", particlesMC.size()); + for (auto& particle : particlesMC) { + if (std::abs(particle.flagMCMatchGen()) == 1 << hf_cand_bs::DecayType::BsToDsPi) { + + auto yParticle = RecoDecay::Y(array{particle.px(), particle.py(), particle.pz()}, RecoDecay::getMassPDG(pdg::Code::kBs)); + if (cutYCandMax >= 0. && std::abs(yParticle) > cutYCandMax) { + continue; + } + + float ptProngs[2], yProngs[2], etaProngs[2]; + for (int iD = particle.daughter0Id(), counter = 0; iD <= particle.daughter1Id(); ++iD, counter++) { + auto daught = particlesMC.iteratorAt(iD); + ptProngs[counter] = daught.pt(); + etaProngs[counter] = daught.eta(); + yProngs[counter] = RecoDecay::Y(array{daught.px(), daught.py(), daught.pz()}, RecoDecay::getMassPDG(daught.pdgCode())); + } + + registry.fill(HIST("hPtProng0Gen"), ptProngs[0], particle.pt()); + registry.fill(HIST("hPtProng1Gen"), ptProngs[1], particle.pt()); + registry.fill(HIST("hYProng0Gen"), yProngs[0], particle.pt()); + registry.fill(HIST("hYProng1Gen"), yProngs[1], particle.pt()); + registry.fill(HIST("hEtaProng0Gen"), etaProngs[0], particle.pt()); + registry.fill(HIST("hEtaProng1Gen"), etaProngs[1], particle.pt()); + + // if (cutYCandMax >= 0. && (std::abs(yProngs[0]) > cutYCandMax || std::abs(yProngs[1]) > cutYCandMax)) + // continue; + + registry.fill(HIST("hPtGen"), particle.pt()); + registry.fill(HIST("hYGen"), yParticle, particle.pt()); + registry.fill(HIST("hEtaGen"), particle.eta(), particle.pt()); + } + } // gen + } // process +}; // struct + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow{}; + const bool doMC = cfgc.options().get("doMC"); + workflow.push_back(adaptAnalysisTask(cfgc)); + if (doMC) { + workflow.push_back(adaptAnalysisTask(cfgc)); + } + return workflow; +}