diff --git a/PWGHF/Core/SelectorCuts.h b/PWGHF/Core/SelectorCuts.h index 2a0917cdf97..0daaad3176d 100644 --- a/PWGHF/Core/SelectorCuts.h +++ b/PWGHF/Core/SelectorCuts.h @@ -843,6 +843,93 @@ static const std::vector labelsPt = { static const std::vector labelsCutVar = {"m", "pT p", "pT K", "pT Pi", "chi2PCA", "decay length", "cos pointing angle", "decLengthXY", "normDecLXY", "ct", "impParXY"}; } // namespace hf_cuts_xic_to_p_k_pi +namespace hf_cuts_to_xi_pi +{ +static constexpr int NBinsPt = 13; +static constexpr int NCutVars = 26; +// 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 BinsPt[NBinsPt + 1] = { + 0., + 1., + 2., + 3., + 4., + 5., + 6., + 7., + 8., + 9., + 10., + 11., + 12., + 20.}; +auto vecBinsPt = std::vector{BinsPt, BinsPt + NBinsPt + 1}; + +// default values for the cuts +constexpr double Cuts[NBinsPt][NCutVars] = { + {2.0, 3.1, 0.01, 0.01, 0.2, 0.15, 1.0, 1.0, 2.0, 0.06, 0.06, 0.04, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.97, 0.97, 0.6, 1.2, 1.0, 0.8}, + {2.0, 3.1, 0.01, 0.01, 0.2, 0.15, 1.0, 1.0, 2.0, 0.06, 0.06, 0.04, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.97, 0.97, 0.6, 1.2, 1.0, 0.8}, + {2.0, 3.1, 0.01, 0.01, 0.2, 0.15, 1.0, 1.0, 2.0, 0.06, 0.06, 0.04, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.97, 0.97, 0.6, 1.2, 1.0, 0.8}, + {2.0, 3.1, 0.01, 0.01, 0.2, 0.15, 1.0, 1.0, 2.0, 0.06, 0.06, 0.04, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.97, 0.97, 0.6, 1.2, 1.0, 0.8}, + {2.0, 3.1, 0.01, 0.01, 0.2, 0.15, 1.0, 1.0, 2.0, 0.06, 0.06, 0.04, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.97, 0.97, 0.6, 1.2, 1.0, 0.8}, + {2.0, 3.1, 0.01, 0.01, 0.2, 0.15, 1.0, 1.0, 2.0, 0.06, 0.06, 0.04, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.97, 0.97, 0.6, 1.2, 1.0, 0.8}, + {2.0, 3.1, 0.01, 0.01, 0.2, 0.15, 1.0, 1.0, 2.0, 0.06, 0.06, 0.04, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.97, 0.97, 0.6, 1.2, 1.0, 0.8}, + {2.0, 3.1, 0.01, 0.01, 0.2, 0.15, 1.0, 1.0, 2.0, 0.06, 0.06, 0.04, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.97, 0.97, 0.6, 1.2, 1.0, 0.8}, + {2.0, 3.1, 0.01, 0.01, 0.2, 0.15, 1.0, 1.0, 2.0, 0.06, 0.06, 0.04, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.97, 0.97, 0.6, 1.2, 1.0, 0.8}, + {2.0, 3.1, 0.01, 0.01, 0.2, 0.15, 1.0, 1.0, 2.0, 0.06, 0.06, 0.04, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.97, 0.97, 0.6, 1.2, 1.0, 0.8}, + {2.0, 3.1, 0.01, 0.01, 0.2, 0.15, 1.0, 1.0, 2.0, 0.06, 0.06, 0.04, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.97, 0.97, 0.6, 1.2, 1.0, 0.8}, + {2.0, 3.1, 0.01, 0.01, 0.2, 0.15, 1.0, 1.0, 2.0, 0.06, 0.06, 0.04, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.97, 0.97, 0.6, 1.2, 1.0, 0.8}, + {2.0, 3.1, 0.01, 0.01, 0.2, 0.15, 1.0, 1.0, 2.0, 0.06, 0.06, 0.04, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.0, 10.0, 0.97, 0.97, 0.6, 1.2, 1.0, 0.8}, +}; + +// row labels +static const std::vector labelsPt = { + "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", + "pT bin 12"}; + +// column labels +static const std::vector labelsCutVar = { + "mCharmBaryonMin", + "mCharmBaryonMax", + "massWindowCascade", + "massWindowV0", // Invariant mass + "ptCharmBachelorMin", + "ptBachelorMin", // pt of daughter + "dcaCascDauMax", + "dcaV0DauMax", + "dcaCharmBaryonDauMax", // dca between daughters + "dcaxyV0PosDauToPvMin", + "dcaxyV0NegDauToPvMin", + "dcaxyBachToPvMin", // DCAXY(impact parameter) + "impactParXYCharmBachelorMin", + "impactParXYCharmBachelorMax", + "impactParZCharmBachelorMin", + "impactParZCharmBachelorMax", // DCAXY(impact parameter) + "impactParXYCascMin", + "impactParXYCascMax", + "impactParZCascMin", + "impactParZCascMax", // DCAXY(impact parameter) + "cosPaCascMin", + "cosPaV0Min", // CosPa + "radiusCascMin", + "radiusV0Min", // radius + "etaTrackLFDauMax", + "etaTrackCharmBachMax" // eta cut to daughter tracks +}; +} // namespace hf_cuts_to_xi_pi + namespace hf_cuts_xic_to_xi_pi_pi { static constexpr int NBinsPt = 13; diff --git a/PWGHF/TableProducer/CMakeLists.txt b/PWGHF/TableProducer/CMakeLists.txt index 26b3400ca28..6a6dde9c3b6 100644 --- a/PWGHF/TableProducer/CMakeLists.txt +++ b/PWGHF/TableProducer/CMakeLists.txt @@ -90,6 +90,11 @@ o2physics_add_dpl_workflow(candidate-creator-xic0-omegac0 PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter KFParticle::KFParticle O2Physics::AnalysisCCDB O2Physics::EventFilteringUtils COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(candidate-creator-xic0-omegac0-qa + SOURCES candidateCreatorXic0Omegac0Qa.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter KFParticle::KFParticle O2Physics::EventFilteringUtils + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(candidate-creator-xic-to-xi-pi-pi SOURCES candidateCreatorXicToXiPiPi.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter KFParticle::KFParticle O2Physics::AnalysisCCDB O2Physics::EventFilteringUtils @@ -187,6 +192,11 @@ o2physics_add_dpl_workflow(candidate-selector-to-xi-pi PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::MLCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(candidate-selector-to-xi-pi-qa + SOURCES candidateSelectorToXiPiQa.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::MLCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(candidate-selector-xic-to-p-k-pi SOURCES candidateSelectorXicToPKPi.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::MLCore @@ -279,6 +289,11 @@ o2physics_add_dpl_workflow(tree-creator-to-xi-pi PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(tree-creator-to-xi-pi-qa + SOURCES treeCreatorToXiPiQa.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(tree-creator-xic0-to-xi-pi-kf SOURCES treeCreatorXic0ToXiPiKf.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore diff --git a/PWGHF/TableProducer/candidateCreatorXic0Omegac0Qa.cxx b/PWGHF/TableProducer/candidateCreatorXic0Omegac0Qa.cxx new file mode 100644 index 00000000000..cd4aabc2077 --- /dev/null +++ b/PWGHF/TableProducer/candidateCreatorXic0Omegac0Qa.cxx @@ -0,0 +1,2024 @@ +// 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 candidateCreatorXic0Omegac0Qa.cxx +/// \brief Reconstruction of Xic0 and Xicp candiates with hadronic decay chain +/// +/// \author Jinhyun Park , Pusan National University +/// \author Krista Smith , Pusan National University + +#ifndef HomogeneousField +#define HomogeneousField // o2-linter: disable=name/macro (required by KFParticle) +#endif + +#include "PWGHF/Core/DecayChannelsLegacy.h" +#include "PWGHF/DataModel/AliasTables.h" +#include "PWGHF/DataModel/CandidateReconstructionTables.h" +#include "PWGHF/DataModel/CandidateSelectionTables.h" +#include "PWGHF/DataModel/TrackIndexSkimmingTables.h" +#include "PWGHF/Utils/utilsBfieldCCDB.h" +#include "PWGHF/Utils/utilsEvSelHf.h" +#include "PWGHF/Utils/utilsTrkCandHf.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGLF/DataModel/mcCentrality.h" +#include "PWGLF/Utils/strangenessBuilderHelper.h" // -> Added to test removal of strangeness builder workflow + +#include "Common/Core/ZorroSummary.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/CollisionAssociationTables.h" +#include "Common/DataModel/EventSelection.h" +#include "Tools/KFparticle/KFUtilities.h" + +#include "CommonConstants/PhysicsConstants.h" +#include "DCAFitter/DCAFitterN.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/RunningWorkflowInfo.h" +#include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/DCA.h" + +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::hf_centrality; +using namespace o2::hf_occupancy; +using namespace o2::constants::physics; +using namespace o2::hf_evsel; + +struct HfCandidateCreatorXic0Omegac0Qa { + + // Cursor to fill tables + struct : ProducesGroup { + // Candidates created with DCAFitter + Produces rowCandToXiPi; + Produces rowCandToOmegaPi; + Produces rowCandToOmegaKa; + // Candidate created with KFParticle + Produces rowCandToXiPiKf; + Produces rowCandToOmegaPiKf; + Produces rowCandToOmegaKaKf; + } cursors; + + // Configurables + struct : ConfigurableGroup { + // Switch for filling histograms + // ----------------------------- + Configurable fillHistograms{"fillHistograms", true, "fill validation plots"}; + // Magnetic field setting from CCDB + // -------------------------------- + Configurable isRun2{"isRun2", false, "enable Run2 or Run3 GRP objects for magnetic field"}; + Configurable ccdbUrl{"ccdbUrl", "https://alice-ccdb.cern.ch", "url of the ccdb object"}; + Configurable ccdbPathLut{"ccdbPathLut", "GLO/Param/MatLUT", "Path for LUT parameterization"}; + Configurable ccdbPathGrp{"ccdbPathGrp", "GLO/GRP/GRP", "path of the group file (Run2)"}; + Configurable ccdbPathGrpMag{"ccdbPathGrpMag", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object (Run3)"}; + + // Options for V0 building + // ...Initial values taken from PWGLF/Utiles/strangenessBuilderModule.h + // --------------------------------------------------------------------- + Configurable minCrossedRowsForV0Building{"minCrossedRowsForV0Building", 50, "minimun TPC crossed rows for daughter tracks. Used for internal V0 Building"}; + Configurable dcanegtopvForV0Building{"dcanegtopvForV0Building", .1, "DCV Neg to PV"}; + Configurable dcapostopvForV0Building{"dcapostopvForV0Building", .1, "DCV Pos To PV"}; + Configurable v0cospaForV0Building{"v0cospaForV0Building", 0.95, "V0 CosPA"}; + Configurable dcav0dauForV0Building{"dcav0dauForV0Building", 1.0, "DCA V0 Daughters"}; + Configurable v0radiusForV0Building{"v0radiusForV0Building", 0.9, "v0radius"}; + Configurable maxDaughterEtaForV0Building{"maxDaughterEtaForV0Building", 5.0, "Maximun daughter eta (in abs value)"}; + + // Options for internal cascade building + // ...Initial values taken from PWGLF/Utiles/strangenessBuilderModule.h + // -------------------------------------------------------------------- + Configurable minCrossedRowsForCascadeBuilding{"minCrossedRowsForCascadeBuilding", 50, "minimun TPC crossed rows for daughter tracks. Used for internal Cascade Building"}; + Configurable dcabachtopvForCascadeBuilding{"dcabachtopvForCascadeBuilding", .1, "DCV Bach to PV"}; + Configurable cascradiusForCascadeBuilding{"cascradiusForCascadeBuilding", .1, "DCV Bach to PV"}; + Configurable casccospaForCascadeBuilding{"casccospaForCascadeBuilding", 0.95, "Cascade CosPA"}; + Configurable dcacascdauForCascadeBuilding{"dcacascdauForCascadeBuilding", 1.0, "DCA cascade daughters"}; + Configurable lambdaMassWindowForCascadeBuilding{"lambdaMassWindowForCascadeBuilding", 0.10, "Distance from Lambda mass(does not apply to KF path)"}; + Configurable maxDaughterEtaForCascadeBuilding{"maxDaughterEtaForCascadeBuilding", 5.0, "Maximun daughter eta (in abs value)"}; + + // Options for internal cascade building - KF Building specifics + // ...Initial values taken from PWGLF/Utiles/strangenessBuilderModule.h + // -------------------------------------------------------------------- + Configurable kfTuneForOmega{"kfTuneForOmega", false, "if enabled, take main cascade properties from omega fit instread of Xi fit(=default)"}; + Configurable kfConstructMethod{"kfConstructMethod", 2, "2 : Daughter particle masses stay fixed in construction process"}; + Configurable kfUseV0MassConstraint{"kfUseV0MassConstraint", true, "KF : Use Lambda mass constraint"}; + Configurable kfUseCascadeMassConstraint{"kfUseCascadeMassConstraint", false, "KF : Use Cascade mass constraint - WARNING : Not adequate for inv mass analysis of Xi"}; + Configurable kfDoDCAFitterPreMinimV0{"kfDoDCAFitterPreMinimV0", true, "KF : do DCAFitter pre-optimization before KF fit to include material correction for V0"}; + Configurable kfDoDCAFitterPreMinimCasc{"kfDoDCAFitterPreMinimCasc", true, "KF : do DCAFitter pre-optimization before KF fit to include material correction for Xi"}; + + // Cascade pre selection + // -------------------- + Configurable doCascadePreselection{"doCascadePreselection", false, "Use invariant mass and dcaXY cuts to preselect cascade candidates"}; + Configurable massToleranceCascade{"massToleranceCascade", 0.01, "Invariant mass tolerance for cascades"}; + Configurable dcaXYToPVCascadeMax{"dcaXYToPVCascadeMax", 3, "Max cascade DCA to PV in XY plane"}; + Configurable dcaV0DaughtersMax{"dcaV0DaughtersMax", 1.0, "Max DCA of V0 daughter"}; + Configurable dcaCascDaughtersMax{"dcaCascDaughtersMax", 1.0, "Max DCA of cascade daughter"}; + + // Options for DCAFitter + // --------------------- + Configurable propagateToPCA{"propagateToPCA", true, "Create tracks version propagated to PCA"}; + Configurable maxR{"maxR", 200., "Reject PCA's above this radius"}; + Configurable maxDZIni{"maxDZIni", 4., "Reject (if>0) PCA candidate if tracks DZ exceeds this threshold"}; + Configurable minParamChange{"minParamChange", 1.e-3, "Stop iteration if largest change of any X is smaller than this"}; + Configurable minRelChi2Change{"minRelChi2Change", 0.9, "Stop iteration if Chi2/Chi2old > this"}; + Configurable useAbsDCA{"useAbsDCA", false, "Minimise abs. distance rather than chi2"}; + Configurable useWeightedFinalPCA{"useWeightedFinalPCA", false, "Recalculate vertex position using track covariance, effective only if useAbsDCA is true"}; + + // Options for KFParticle + // ---------------------- + Configurable rejDiffCollTrack{"rejDiffCollTrack", true, "Reject tracks comming from different collisions(effective only for KFParticle w/o derived data)"}; + Configurable kfDoCascadePreselection{"kfDoCascadePreselection", false, "Use invariant mass and dcaXY cuts to preselect cascade candidates"}; + // Configurable kfConstrainTopoV0ToCasc{"kfConstrainTopoV0ToCasc", false, "KF : Use Lambda topo constraint"}; // -> Not sure if this will be used... + Configurable kfConstrainInvMassV0{"kfConstrainInvMassV0", true, "use mass constraint for cascade"}; + Configurable kfConstrainInvMassCasc{"kfConstrainInvMassCasc", true, "use mass constraint for cascade"}; + // Configurable kfConstrainTopoCascToCharmBaryon{"kfConstrainTopoCascToCharmBaryon", true, "use topo constraint for cascade"}; + Configurable kfConstrainInvMassCharmBaryon{"kfConstrainInvMassCharmBaryon", false, "constrain invariant mass of charm baryon candidate"}; + // Configurable constrainTopoCascToChramBaryon{"constrainTopoCascToChramBaryon", false, "constrain Casc to Charm Baryon"}; + // Configurable kfConstrainTopoCharmBaryon{"kfConstrainTopoCharmBaryon", false, "constrain charm baryon candidate to decay vertex"}; + + // Options for QA histogram binning + // ----------------------------- + + // For Cascade + Configurable nBinMassCasc{"nBinMassCasc", 1000, "nBinCascMass"}; + Configurable minMassCasc{"minMassCasc", 1.0, "xiMassMin"}; + Configurable maxMassCasc{"maxMassCasc", 2.0, "xiMassMax"}; + Configurable nBinPtCasc{"nBinPtCasc", 100, "nBinPtXi"}; + Configurable minPtCasc{"minPtCasc", 0.0, "minimun value of cascade"}; + Configurable maxPtCasc{"maxPtCasc", 20.0, "maximum value of cascade"}; + + // For Prong0 + Configurable nBinMassProng0{"nBinMassProng0", 100, "nBinMassPi"}; + Configurable minMassProng0{"minMassProng0", 0.0, "ptPiMin"}; + Configurable maxMassProng0{"maxMassProng0", 20.0, "ptPiMax"}; + Configurable nBinPtProng0{"nBinPtProng0", 100, "nBinPtProng0"}; + Configurable minPtProng0{"minPtProng0", 0.0, "ptPiMin"}; + Configurable maxPtProng0{"maxPtProng0", 20.0, "ptPiMax"}; + + // For Charm Baryon + Configurable nBinMassCharmBaryon{"nBinMassCharmBaryon", 3000, "nBinXic0Mass"}; + Configurable minMassCharmBaryon{"minMassCharmBaryon", 1.0, "xic0MassMin"}; + Configurable maxMassCharmBaryon{"maxMassCharmBaryon", 4.0, "xic0MassMax"}; + Configurable nBinPtCharmBaryon{"nBinPtCharmBaryon", 100, "nBinXic0Pt"}; + Configurable minPtCharmBaryon{"minPtCharmBaryon", 0.0, "xic0PtMin"}; + Configurable maxPtCharmBaryon{"maxPtCharmBaryon", 20.0, "xic0PtMax"}; + + // Etc + Configurable nBinCpa2Prong{"nBinCpa2Prong", 240, "nBinCpa2Prong"}; + Configurable cpa2ProngMin{"cpa2ProngMin", -1.2, "cpa2ProngMin"}; + Configurable cpa2ProngMax{"cpa2ProngMax", 1.2, "cpa2ProngMax"}; + Configurable nBinImpParXYXi{"nBinImpParXYXi", 30, "nBinImpParXYXi"}; + Configurable impParXYXiMin{"impParXYXiMin", -1.5, "impParXYXiMin"}; + Configurable impParXYXiMax{"impParXYXiMax", 1.5, "impParXYXiMax"}; + Configurable nBinImpParXYPi{"nBinImpParXYPi", 30, "nBinImpParXYPi"}; + Configurable impParXYPiMin{"impParXYPiMin", -1.5, "impParXYPiMin"}; + Configurable impParXYPiMax{"impParXYPiMax", 1.5, "impParXYPiMax"}; + } configs; + + // For magnetic field + Service ccdb; + o2::base::MatLayerCylSet* lut; + o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrLUT; + + // DCAFitter + o2::vertexing::DCAFitterN<2> df; + + int runNumber{0}; + double magneticField{0.}; + // float massCharmBaryonCand{0.}; + + enum CharmBaryonCandCounter { All = 0, + HfFlagPass, + CascReconstructed, + VertexFit }; + + // Table aliases + using SelectedCollisions = soa::Join; + using TracksWCovDcaExtraPidPrPiKa = soa::Join; + + HistogramRegistry registry{"hists"}; + OutputObj zorroSummary{"zorroSummary"}; + HfEventSelection hfEvSel; + + // For cascade building + o2::pwglf::strangenessBuilderHelper straHelper; + + // For candidate reconstruction in different PID hypothesis + // Each decay channel assgined by following hf_cand_casc_lf_DecayType2Prong enum + // Index 0: Xic0/Omegac0 -> Xi- Pi+ Index 1: Omegac0 -> Omega- Pi+ Index 2: Omegac0 -> Omega- K+ + // Xi- -> Lambda0 pi-, Omeag- -> Lambda0 K- + // Lambda0 -> Pr+ Pi- + std::array pdgOfCharmBaryon = {+kXiC0, +kOmegaC0, +kOmegaC0}; // -> This need to be fixed...? +kXiC0, +kOmegaC0, +kOmegaC0, +kOmeagC0 + std::array pdgOfCharmBach = {+kPiPlus, +kPiPlus, +kKPlus}; + std::array pdgOfCascade = {+kXiMinus, +kOmegaMinus, +kOmegaMinus}; + std::array pdgOfBach = {+kPiMinus, +kKMinus, +kKMinus}; + std::array pdgOfV0 = {+kLambda0, +kLambda0, +kLambda0}; + std::array pdgOfV0DauPos = {+kProton, +kProton, +kProton}; + std::array pdgOfV0DauNeg = {+kPiMinus, +kPiMinus, +kPiMinus}; + + std::array trackPidOfCascade = {o2::track::PID::XiMinus, o2::track::PID::OmegaMinus, o2::track::PID::OmegaMinus}; + std::array massOfCascades = {MassXiMinus, MassOmegaMinus, MassOmegaMinus}; + std::array massOfCharmBach = {MassPiPlus, MassPiPlus, MassKPlus}; + + // Pointer of histograms for QA + std::shared_ptr hInvMassCharmBaryonToXiPi, hInvMassCharmBaryonToOmegaPi, hInvMassCharmBaryonToOmegaKa; + std::shared_ptr hCandidateCounterToXiPi, hCandidateCounterToOmegaPi, hCandidateCounterToOmegaKa; + + void init(InitContext const&) + { + std::vector processesToXiPiDca{doprocessToXiPiWithDCAFitterNoCent, /*doprocessToXiPiWithDCAFitterNoCentWithTrackedCasc,*/ doprocessToXiPiWithDCAFitterCentFT0C, doprocessToXiPiWithDCAFitterCentFT0M}; + std::vector processesToOmegaPiDca{doprocessToOmegaPiWithDCAFitterNoCent, doprocessToOmegaPiWithDCAFitterCentFT0C, doprocessToOmegaPiWithDCAFitterCentFT0M}; + std::vector processesToOmegaKaDca{doprocessToOmegaKaWithDCAFitterNoCent, doprocessToOmegaKaWithDCAFitterCentFT0C, doprocessToOmegaKaWithDCAFitterCentFT0M}; + std::vector processesToXiPiKf{doprocessToXiPiWithKFParticleNoCent, doprocessToXiPiWithKFParticleCentFT0C, doprocessToXiPiWithKFParticleCentFT0M}; + std::vector processesToOmegaPiKf{doprocessToOmegaPiWithKFParticleNoCent, doprocessToOmegaPiWithKFParticleCentFT0C, doprocessToOmegaPiWithKFParticleCentFT0M}; + std::vector processesToOmegaKaKf{doprocessToOmegaKaWithKFParticleNoCent, doprocessToOmegaKaWithKFParticleCentFT0C, doprocessToOmegaKaWithKFParticleCentFT0M}; + std::vector processesCollMonitoring{doprocessCollisionsNoCent, doprocessCollisionsCentFT0C, doprocessCollisionsCentFT0M}; + + int xipiEnabledDca = std::accumulate(processesToXiPiDca.begin(), processesToXiPiDca.end(), 0); + int xipiEnabledKf = std::accumulate(processesToXiPiKf.begin(), processesToXiPiKf.end(), 0); + int omegapiEnabledDca = std::accumulate(processesToOmegaPiDca.begin(), processesToOmegaPiDca.end(), 0); + int omegapiEnabledKf = std::accumulate(processesToOmegaPiKf.begin(), processesToOmegaPiKf.end(), 0); + int omegakaEnabledDca = std::accumulate(processesToOmegaKaDca.begin(), processesToOmegaKaDca.end(), 0); + int omegakaEnabledKf = std::accumulate(processesToOmegaKaKf.begin(), processesToOmegaKaKf.end(), 0); + + // Exit if workflow is not configured correctly - More than one process function enabled for candidate to XiPi + if ((xipiEnabledDca > 0) && (xipiEnabledKf > 0)) { + LOGP(fatal, "More than one process function enabled for candidte decaying to xi pi"); + } + + // Exit if workflow is not configured correctly - More than one process function enabled for candidate to OmegaPi + if ((omegapiEnabledDca > 0) && (omegapiEnabledKf > 0)) { + LOGP(fatal, "More than one process function enabled for candidte decaying to omega pi"); + } + + // Exit if workflow is not configured correctly - More than one process function enabled for candidate to OmegaKa + if ((omegakaEnabledDca > 0) && (omegakaEnabledKf > 0)) { + LOGP(fatal, "More than one process function enabled for candidte decaying to omega ka"); + } + + // Exit if workflow is not configured correctly - More than one process enabled for collision monitoring + if (std::accumulate(processesCollMonitoring.begin(), processesCollMonitoring.end(), 0) > 1) { + LOGP(fatal, "More than one process fucntion for CollMonitoring was enabled. Please choose only one process function"); + } + + // Add histogram to indicate which sv method was used + registry.add("hVertexerType", "Use KF or DCAFitterN;Vertexer type;entries", {kTH1F, {{2, 0.0, 2.0}}}); + registry.get(HIST("hVertexerType"))->GetXaxis()->SetBinLabel(1 + aod::hf_cand::VertexerType::DCAFitter, "DCAFitter"); + registry.get(HIST("hVertexerType"))->GetXaxis()->SetBinLabel(1 + aod::hf_cand::VertexerType::KfParticle, "KFParticle"); + + // initialize ccdb + // --------------- + ccdb->setURL(configs.ccdbUrl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->get(configs.ccdbPathLut)); + straHelper.lut = lut; + runNumber = 0; + + // Initilization for strangeness builder helper + // -------------------------------------------- + + // Settings for internal V0 building + straHelper.v0selections.minCrossedRows = configs.minCrossedRowsForV0Building; + straHelper.v0selections.dcanegtopv = configs.dcanegtopvForV0Building; + straHelper.v0selections.dcapostopv = configs.dcapostopvForV0Building; + straHelper.v0selections.v0cospa = configs.v0cospaForV0Building; + straHelper.v0selections.dcav0dau = configs.dcav0dauForV0Building; + straHelper.v0selections.v0radius = configs.v0radiusForV0Building; + straHelper.v0selections.maxDaughterEta = configs.maxDaughterEtaForV0Building; + + // Settings for internal Cascade building + straHelper.cascadeselections.minCrossedRows = configs.minCrossedRowsForCascadeBuilding; + straHelper.cascadeselections.dcabachtopv = configs.dcabachtopvForCascadeBuilding; + straHelper.cascadeselections.cascradius = configs.cascradiusForCascadeBuilding; + straHelper.cascadeselections.casccospa = configs.casccospaForCascadeBuilding; + straHelper.cascadeselections.dcacascdau = configs.dcacascdauForCascadeBuilding; + straHelper.cascadeselections.lambdaMassWindow = configs.lambdaMassWindowForCascadeBuilding; + straHelper.cascadeselections.maxDaughterEta = configs.maxDaughterEtaForCascadeBuilding; + + // Fitter setting + straHelper.fitter.setPropagateToPCA(configs.propagateToPCA); + straHelper.fitter.setMaxR(configs.maxR); + straHelper.fitter.setMaxDZIni(configs.maxDZIni); + straHelper.fitter.setMinParamChange(configs.minParamChange); + straHelper.fitter.setUseAbsDCA(configs.useAbsDCA); + straHelper.fitter.setWeightedFinalPCA(configs.useWeightedFinalPCA); + + // Extra initialization for DCAFitter + // ---------------------------------- + if (xipiEnabledDca == 1 || omegapiEnabledDca == 1 || omegakaEnabledDca == 1) { + registry.get(HIST("hVertexerType"))->Fill(aod::hf_cand::VertexerType::DCAFitter); + df.setPropagateToPCA(configs.propagateToPCA); + df.setMaxR(configs.maxR); + df.setMaxDZIni(configs.maxDZIni); + df.setMinParamChange(configs.minParamChange); + df.setUseAbsDCA(configs.useAbsDCA); + df.setWeightedFinalPCA(configs.useWeightedFinalPCA); + } + + // Extra initialization for KFParticle + // ----------------------------------- + if (xipiEnabledKf == 1 || omegapiEnabledKf == 1 || omegakaEnabledKf == 1) { + registry.get(HIST("hVertexerType"))->Fill(aod::hf_cand::VertexerType::KfParticle); + } + + // initailize HF event selection helper + // ------------------------------------ + hfEvSel.init(registry, zorroSummary); + + // Histograms for QA + // ----------------- + registry.add("ReconstructedDecayChannel", "DecayChannel", {kTH1F, {{3, 0.0, 3.0}}}); + registry.get(HIST("ReconstructedDecayChannel"))->GetXaxis()->SetBinLabel(1 + hf_cand_casc_lf::DecayType2Prong::XiczeroOmegaczeroToXiPi, "To #Xi #pi"); + registry.get(HIST("ReconstructedDecayChannel"))->GetXaxis()->SetBinLabel(1 + hf_cand_casc_lf::DecayType2Prong::OmegaczeroToOmegaPi, "To #Omega #pi"); + registry.get(HIST("ReconstructedDecayChannel"))->GetXaxis()->SetBinLabel(1 + hf_cand_casc_lf::DecayType2Prong::OmegaczeroToOmegaK, "To #Omega K"); + + if (xipiEnabledDca != 0 || xipiEnabledKf != 0) { + hInvMassCharmBaryonToXiPi = registry.add("hInvMassCharmBaryonToXiPi", "Charm baryon invariant mass - #Xi #pi decay;inv. mass (GeV/#it{c}^{2});entries", {HistType::kTH1D, {{configs.nBinMassCharmBaryon, configs.minMassCharmBaryon, configs.maxMassCharmBaryon}}}); + hCandidateCounterToXiPi = registry.add("hCandidateCounterToXiPi", "Candidate counter wrt derived data - #Xi #pi decay;status;entries", {HistType::kTH1D, {{4, -0.5, 3.5}}}); + registry.get(HIST("hCandidateCounterToXiPi"))->GetXaxis()->SetBinLabel(1 + All, "Total"); + registry.get(HIST("hCandidateCounterToXiPi"))->GetXaxis()->SetBinLabel(1 + HfFlagPass, "HfFlagPass"); + registry.get(HIST("hCandidateCounterToXiPi"))->GetXaxis()->SetBinLabel(1 + CascReconstructed, "CascReconstructed"); + registry.get(HIST("hCandidateCounterToXiPi"))->GetXaxis()->SetBinLabel(1 + VertexFit, "VertexFit"); + registry.get(HIST("ReconstructedDecayChannel"))->Fill(hf_cand_casc_lf::DecayType2Prong::XiczeroOmegaczeroToXiPi); + } + + if (omegapiEnabledDca != 0 || omegapiEnabledKf != 0) { + hInvMassCharmBaryonToOmegaPi = registry.add("hInvMassCharmBaryonToOmegaPi", "Charm baryon invariant mass - #Omega #pi decay;inv. mass (GeV/#it{c}^{2});entries", {HistType::kTH1D, {{configs.nBinMassCharmBaryon, configs.minMassCharmBaryon, configs.maxMassCharmBaryon}}}); + hCandidateCounterToOmegaPi = registry.add("hCandidateCounterToOmegaPi", "Candidate counter wrt derived data - #Omega #pi decay;status;entries", {HistType::kTH1D, {{4, -0.5, 3.5}}}); + registry.get(HIST("hCandidateCounterToOmegaPi"))->GetXaxis()->SetBinLabel(1 + All, "Total"); + registry.get(HIST("hCandidateCounterToOmegaPi"))->GetXaxis()->SetBinLabel(1 + HfFlagPass, "HfFlagPass"); + registry.get(HIST("hCandidateCounterToOmegaPi"))->GetXaxis()->SetBinLabel(1 + CascReconstructed, "CascReconstructed"); + registry.get(HIST("hCandidateCounterToOmegaPi"))->GetXaxis()->SetBinLabel(1 + VertexFit, "VertexFit"); + registry.get(HIST("ReconstructedDecayChannel"))->Fill(hf_cand_casc_lf::DecayType2Prong::OmegaczeroToOmegaPi); + } + + if (omegakaEnabledDca != 0 || omegakaEnabledKf != 0) { + hInvMassCharmBaryonToOmegaKa = registry.add("hInvMassCharmBaryonToOmegaKa", "Charm baryon invariant mass - #Omega K decay;inv. mass (GeV/#it{c}^{2});entries", {HistType::kTH1D, {{configs.nBinMassCharmBaryon, configs.minMassCharmBaryon, configs.maxMassCharmBaryon}}}); + hCandidateCounterToOmegaKa = registry.add("hCandidateCounterToOmegaKa", "Candidate counter wrt derived data - #Omega K decay;status;entries", {HistType::kTH1D, {{4, -0.5, 3.5}}}); + registry.get(HIST("hCandidateCounterToOmegaKa"))->GetXaxis()->SetBinLabel(1 + All, "Total"); + registry.get(HIST("hCandidateCounterToOmegaKa"))->GetXaxis()->SetBinLabel(1 + HfFlagPass, "HfFlagPass"); + registry.get(HIST("hCandidateCounterToOmegaKa"))->GetXaxis()->SetBinLabel(1 + CascReconstructed, "CascReconstructed"); + registry.get(HIST("hCandidateCounterToOmegaKa"))->GetXaxis()->SetBinLabel(1 + VertexFit, "VertexFit"); + registry.get(HIST("ReconstructedDecayChannel"))->Fill(hf_cand_casc_lf::DecayType2Prong::OmegaczeroToOmegaK); + } + + registry.add("hCascMass", "Inv mass of reconstructed cascade;Inv mass;Entries", {HistType::kTH1F, {{configs.nBinMassCasc, configs.minMassCasc, configs.maxMassCasc}}}); + registry.add("hCascPt", "Pt of reconstructed cascade;pT;Entries", {HistType::kTH1F, {{configs.nBinPtCasc, configs.minPtCasc, configs.maxPtCasc}}}); + + } // end of initialization + + //////////////////////////////////////////////////////////// + // // + // Candidate reconstruction with DCAFitter // + // // + //////////////////////////////////////////////////////////// + + // template function for running charm baryon reconstruction with DCAFitter + /// \brief centEstimator is for different centrality estimators + /// \brief decayChannel is for different decay channels. 0 for XiczeroOmegaczeroToXiPi, 1 for OmegaczeroToOmegaPi, 2 for OmegaczeroToOmeagaK + /// \brief Colls is for collision tables joined with different centraltiy estimators + /// \brief Hist is for QA histograms + template + void runCreatorWithDCAFitter(Colls const&, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const&, // -> Internal cascade building + aod::V0s const&, // -> Internal v0 building + TracksWCovDcaExtraPidPrPiKa const&, + aod::BCsWithTimestamps const&, + Hist& hInvMassCharmBaryon, + Hist& hCandCounter) + { + // arrays which holds values for different decay channel reconstruction + + // Loop over candidate + for (auto const& cand : candidates) { + + // Fill cascandidates before selection + if (configs.fillHistograms) { + hCandCounter->Fill(All); + } + + // Apply hfflag selection for different candidate reconstruction + if (!TESTBIT(cand.hfflag(), decayChannel)) { + continue; + } else { + if (configs.fillHistograms) { + hCandCounter->Fill(HfFlagPass); + } + } + + // Event selection + auto collision = cand.collision_as(); + float centrality{-1.f}; + const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); + if (rejectionMask != 0) { // None of the event selection satisfied -> Reject this candidate + continue; + } + + //------------------------------Set Magnetic field------------------------------ + auto bc = collision.template bc_as(); + if (runNumber != bc.runNumber()) { + LOG(info) << ">>>>>>>>>> Current run Number : " << runNumber; + initCCDB(bc, runNumber, ccdb, configs.isRun2 ? configs.ccdbPathGrp : configs.ccdbPathGrpMag, lut, configs.isRun2); + magneticField = o2::base::Propagator::Instance()->getNominalBz(); + LOG(info) << ">>>>>>>>>> Magnetic field: " << magneticField; + } + straHelper.fitter.setBz(magneticField); // -> Magnetic field setting for internal cascade building + df.setBz(magneticField); // -> Magnetic field setting for charm baryon building + + //------------------Intenal Cascade building------------------ + auto cascAodElement = cand.cascade_as(); + auto v0AodElement = cascAodElement.v0_as(); + auto posTrack = v0AodElement.posTrack_as(); + auto negTrack = v0AodElement.negTrack_as(); + auto bachTrack = cascAodElement.bachelor_as(); + + // Make cascade starting from V0 + // If success, fill Cascade and V0 information for reconstruction + if (!straHelper.buildCascadeCandidate(collision.globalIndex(), + collision.posX(), collision.posY(), collision.posZ(), + posTrack, + negTrack, + bachTrack, + false, // calculateBachelorBaryonVariable + false, // useCascadeMomentumAtPV + true)) { + LOG(info) << "!This cascade cannot be rebuilt(cascade ID : " << cand.cascadeId() << "/ collision ID : " << collision.globalIndex(); + continue; + } else { + float storeMass = (decayChannel == 0) ? straHelper.cascade.massXi : straHelper.cascade.massOmega; + float storePt = RecoDecay::pt(straHelper.cascade.cascadeMomentum); + registry.fill(HIST("hCascMass"), storeMass); + registry.fill(HIST("hCascPt"), storePt); + } + + //------------------Cascade pre selection------------------ + float massCasc = (decayChannel == 0) ? straHelper.cascade.massXi : straHelper.cascade.massOmega; + + // Perform cascade pre selection + if (configs.doCascadePreselection) { + + // pre selection of dcaXY + // FIXME does cascadeDCAxy represents dcaXYtoPV? + if (std::abs(straHelper.cascade.cascadeDCAxy) > configs.dcaXYToPVCascadeMax) { + continue; + } + + // pre selection on invariant mass + if (std::abs(massCasc - massOfCascades[decayChannel]) > configs.massToleranceCascade) { + continue; + } + } + + if (configs.fillHistograms) { + hCandCounter->Fill(CascReconstructed); + } + + //------------------------------Info of V0 and cascade------------------------------ + // -> This quantities are used for physical properties of selected candidates + // -> Not used for candidate creation + std::array vertexV0 = {straHelper.cascade.v0Position[0], straHelper.cascade.v0Position[1], straHelper.cascade.v0Position[2]}; + std::array pVecV0 = {straHelper.cascade.v0Momentum[0], straHelper.cascade.v0Momentum[1], straHelper.cascade.v0Momentum[2]}; + std::array pVecV0DauPos = {straHelper.cascade.positiveMomentum[0], straHelper.cascade.positiveMomentum[1], straHelper.cascade.positiveMomentum[2]}; + std::array pVecV0DauNeg = {straHelper.cascade.negativeMomentum[0], straHelper.cascade.negativeMomentum[1], straHelper.cascade.negativeMomentum[2]}; + + int chargeCasc = straHelper.cascade.charge > 0 ? 1 : -1; + std::array vertexCasc = {straHelper.cascade.cascadePosition[0], straHelper.cascade.cascadePosition[1], straHelper.cascade.cascadePosition[2]}; + std::array pVecCasc = {straHelper.cascade.cascadeMomentum[0], straHelper.cascade.cascadeMomentum[1], straHelper.cascade.cascadeMomentum[2]}; + std::array covCasc = {0.}; + + //------------------------------Create cascade track------------------------------ + + constexpr std::size_t NElementsCovMatrix{6u}; + constexpr std::array MomInd = {9, 13, 14, 18, 19, 20}; // Momentum index? + for (auto i = 0u; i < NElementsCovMatrix; i++) { + covCasc[i] = straHelper.cascade.covariance[i]; + covCasc[MomInd[i]] = straHelper.cascade.covariance[MomInd[i]]; + } + + o2::track::TrackParCov trackCasc; + if (chargeCasc < 0) { // Xi- or Omega- + trackCasc = o2::track::TrackParCov(vertexCasc, pVecCasc, covCasc, -1, true); + } else if (chargeCasc > 0) { // Xi+ or Omega+ + trackCasc = o2::track::TrackParCov(vertexCasc, pVecCasc, covCasc, 1, true); + } else { + continue; + } + + trackCasc.setAbsCharge(1); + trackCasc.setPID(trackPidOfCascade[decayChannel]); + + //------------------------------Fit SV & Create Charm Baryon track------------------------------ + + // Perform secondary vertex fitting + auto trackCharmBachelor = cand.prong0_as(); + auto trackParCovCharmBachelor = getTrackParCov(trackCharmBachelor); + try { + if (df.process(trackCasc, trackParCovCharmBachelor) == 0) { + continue; + } + } catch (const std::runtime_error& e) { + LOG(info) << "Run time error found : " << e.what() << ".DCAFitter cannot work with this candidate. SKIP!"; + continue; + } + + // Propagate tracks to vertex + df.propagateTracksToVertex(); + if (!df.isPropagateTracksToVertexDone()) { + continue; + } + + if (configs.fillHistograms) { + hCandCounter->Fill(VertexFit); + } + + //------------------------------Calculate physical properties----------------------------- + + // get track momenta + std::array pVecCascAsD, pVecCharmBachAsD; + df.getTrack(0).getPxPyPzGlo(pVecCascAsD); + df.getTrack(1).getPxPyPzGlo(pVecCharmBachAsD); + + std::array pVecCharmBaryon = {pVecCascAsD[0] + pVecCharmBachAsD[0], pVecCascAsD[1] + pVecCharmBachAsD[1], pVecCascAsD[2] + pVecCharmBachAsD[2]}; + std::array pVecBach = {straHelper.cascade.bachelorMomentum[0], straHelper.cascade.bachelorMomentum[1], straHelper.cascade.bachelorMomentum[2]}; + + // get PV Properties + auto primaryVertex = getPrimaryVertex(collision); + auto covMatrixPV = primaryVertex.getCov(); + std::array pvCoord = {collision.posX(), collision.posY(), collision.posZ()}; + + // get SV Properties + auto const& secondaryVertex = df.getPCACandidate(); + auto chi2SV = df.getChi2AtPCACandidate(); + auto covMatrixSV = df.calcPCACovMatrixFlat(); + + // DCAxy and DCAz. Computed with propagatToDCABxByBz method + auto trackParCovV0DauPos = getTrackParCov(posTrack); + auto trackParCovV0DauNeg = getTrackParCov(negTrack); + auto trackParCovBach = getTrackParCov(bachTrack); + + o2::dataformats::DCA impactParameterV0DauPos, impactParameterV0DauNeg, impactParameterBach; + o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackParCovV0DauPos, 2.f, matCorr, &impactParameterV0DauPos); + o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackParCovV0DauNeg, 2.f, matCorr, &impactParameterV0DauNeg); + o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackParCovBach, 2.f, matCorr, &impactParameterBach); + + float dcaxyV0DauPos = impactParameterV0DauPos.getY(); + float dcaxyV0DauNeg = impactParameterV0DauNeg.getY(); + float dcaxyBach = impactParameterBach.getY(); + float dcazV0DauPos = impactParameterV0DauPos.getZ(); + float dcazV0DauNeg = impactParameterV0DauNeg.getZ(); + float dcazBach = impactParameterBach.getZ(); + + // get impact parameter. Compute with propagateToDCABxByBz method + o2::dataformats::DCA impactParameterCasc, impactParameterCharmBach; + o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackCasc, 2.f, matCorr, &impactParameterCasc); // trackCasc is TrackParCov object + o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackParCovCharmBachelor, 2.f, matCorr, &impactParameterCharmBach); + + // get v0 invariant mass + float massLambda = straHelper.v0.massLambda; + + // Get Charm baryon invarian mass + auto arrMomenta = std::array{pVecCascAsD, pVecCharmBachAsD}; + float massCharmBaryonCand = RecoDecay::m(arrMomenta, std::array{massOfCascades[decayChannel], massOfCharmBach[decayChannel]}); + if (configs.fillHistograms) { + hInvMassCharmBaryon->Fill(massCharmBaryonCand); + } + + // calculate cosine of pointing angle + std::array vtxCoordCharmBaryon = df.getPCACandidatePos(); + float cpaV0 = RecoDecay::cpa(pvCoord, vertexV0, pVecV0); + float cpaCasc = RecoDecay::cpa(pvCoord, vertexCasc, pVecCasc); + float cpaCharmBaryon = RecoDecay::cpa(pvCoord, vtxCoordCharmBaryon, pVecCharmBaryon); + float cpaxyV0 = RecoDecay::cpaXY(pvCoord, vertexV0, pVecV0); + float cpaxyCasc = RecoDecay::cpaXY(pvCoord, vertexCasc, pVecCasc); + float cpaxyCharmBaryon = RecoDecay::cpaXY(pvCoord, vtxCoordCharmBaryon, pVecCharmBaryon); + + // calculate decay length + float decLenV0 = RecoDecay::distance(vertexCasc, vertexV0); + float decLenCasc = RecoDecay::distance(vtxCoordCharmBaryon, vertexCasc); + float decLenCharmBaryon = RecoDecay::distance(pvCoord, vtxCoordCharmBaryon); + // get uncertainty of the decay length -> Used in previous code... + float phi, theta; + getPointDirection(std::array{primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, secondaryVertex, phi, theta); + auto errorDecayLength = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixSV, phi, theta)); + auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixSV, phi, 0.)); + + // calcuate ctau + float ctV0 = RecoDecay::ct(pVecV0, decLenV0, MassLambda0); + float ctCasc = RecoDecay::ct(pVecCasc, decLenCasc, massOfCascades[decayChannel]); + float ctOmegac0 = RecoDecay::ct(pVecCharmBaryon, decLenCharmBaryon, MassOmegaC0); + float ctXic0 = RecoDecay::ct(pVecCharmBaryon, decLenCharmBaryon, MassXiC0); + + // get eta + float etaV0DauPos = posTrack.eta(); + float etaV0DauNeg = negTrack.eta(); + float etaBach = bachTrack.eta(); + float etaCharmBach = trackCharmBachelor.eta(); + float etaV0 = RecoDecay::eta(pVecV0); + float etaCasc = RecoDecay::eta(pVecCasc); + float etaCharmBaryon = RecoDecay::eta(pVecCharmBaryon); + + // DCA between daughters + float dcaV0Dau = straHelper.cascade.v0DaughterDCA; + float dcaCascDau = straHelper.cascade.cascadeDaughterDCA; + float dcaCharmBaryonDau = std::sqrt(df.getChi2AtPCACandidate()); + + //------------------------------Fill QA histograms----------------------------- + + //------------------------------Fill the table----------------------------- + if constexpr (decayChannel == hf_cand_casc_lf::DecayType2Prong::XiczeroOmegaczeroToXiPi) { + cursors.rowCandToXiPi(collision.globalIndex(), + pvCoord[0], pvCoord[1], pvCoord[2], + secondaryVertex[0], secondaryVertex[1], secondaryVertex[2], + vertexCasc[0], vertexCasc[1], vertexCasc[2], + vertexV0[0], vertexV0[1], vertexV0[2], + bachTrack.sign(), + covMatrixSV[0], covMatrixSV[1], covMatrixSV[2], covMatrixSV[3], covMatrixSV[4], covMatrixSV[5], + pVecCharmBaryon[0], pVecCharmBaryon[1], pVecCharmBaryon[2], + pVecCascAsD[0], pVecCascAsD[1], pVecCascAsD[2], + pVecCharmBachAsD[0], pVecCharmBachAsD[1], pVecCharmBachAsD[2], + pVecV0[0], pVecV0[1], pVecV0[2], + pVecBach[0], pVecBach[1], pVecBach[2], + pVecV0DauPos[0], pVecV0DauPos[1], pVecV0DauPos[2], + pVecV0DauNeg[0], pVecV0DauNeg[1], pVecV0DauNeg[2], + impactParameterCasc.getY(), impactParameterCharmBach.getY(), + impactParameterCasc.getZ(), impactParameterCharmBach.getZ(), + std::sqrt(impactParameterCasc.getSigmaY2()), std::sqrt(impactParameterCharmBach.getSigmaY2()), + cascAodElement.v0Id(), v0AodElement.posTrackId(), v0AodElement.negTrackId(), + cand.cascadeId(), trackCharmBachelor.globalIndex(), cand.prong0Id(), + massLambda, massCasc, massCharmBaryonCand, + cpaV0, cpaCharmBaryon, cpaCasc, + cpaxyV0, cpaxyCharmBaryon, cpaxyCasc, + ctOmegac0, ctCasc, ctV0, ctXic0, + etaV0DauPos, etaV0DauNeg, etaBach, etaCharmBach, etaCharmBaryon, etaCasc, etaV0, + dcaxyV0DauPos, dcaxyV0DauNeg, dcaxyBach, + dcazV0DauPos, dcazV0DauNeg, dcazBach, + dcaCascDau, dcaV0Dau, dcaCharmBaryonDau, + decLenCharmBaryon, decLenCasc, decLenV0, errorDecayLength, errorDecayLengthXY); + } else if constexpr (decayChannel == hf_cand_casc_lf::DecayType2Prong::OmegaczeroToOmegaPi) { + cursors.rowCandToOmegaPi(collision.globalIndex(), + pvCoord[0], pvCoord[1], pvCoord[2], + secondaryVertex[0], secondaryVertex[1], secondaryVertex[2], + vertexCasc[0], vertexCasc[1], vertexCasc[2], + vertexV0[0], vertexV0[1], vertexV0[2], + bachTrack.sign(), + covMatrixSV[0], covMatrixSV[1], covMatrixSV[2], covMatrixSV[3], covMatrixSV[4], covMatrixSV[5], + pVecCharmBaryon[0], pVecCharmBaryon[1], pVecCharmBaryon[2], + pVecCascAsD[0], pVecCascAsD[1], pVecCascAsD[2], + pVecCharmBachAsD[0], pVecCharmBachAsD[1], pVecCharmBachAsD[2], + pVecV0[0], pVecV0[1], pVecV0[2], + pVecBach[0], pVecBach[1], pVecBach[2], + pVecV0DauPos[0], pVecV0DauPos[1], pVecV0DauPos[2], + pVecV0DauNeg[0], pVecV0DauNeg[1], pVecV0DauNeg[2], + impactParameterCasc.getY(), impactParameterCharmBach.getY(), + impactParameterCasc.getZ(), impactParameterCharmBach.getZ(), + std::sqrt(impactParameterCasc.getSigmaY2()), std::sqrt(impactParameterCharmBach.getSigmaY2()), + cascAodElement.v0Id(), v0AodElement.posTrackId(), v0AodElement.negTrackId(), + cand.cascadeId(), trackCharmBachelor.globalIndex(), cand.prong0Id(), + massLambda, massCasc, massCharmBaryonCand, + cpaV0, cpaCharmBaryon, cpaCasc, + cpaxyV0, cpaxyCharmBaryon, cpaxyCasc, + ctOmegac0, ctCasc, ctV0, + etaV0DauPos, etaV0DauNeg, etaBach, etaCharmBach, etaCharmBaryon, etaCasc, etaV0, + dcaxyV0DauPos, dcaxyV0DauNeg, dcaxyBach, + dcazV0DauPos, dcazV0DauNeg, dcazBach, + dcaCascDau, dcaV0Dau, dcaCharmBaryonDau, + decLenCharmBaryon, decLenCasc, decLenV0, errorDecayLength, errorDecayLengthXY, cand.hfflag()); + } else { + cursors.rowCandToOmegaKa(collision.globalIndex(), + pvCoord[0], pvCoord[1], pvCoord[2], + secondaryVertex[0], secondaryVertex[1], secondaryVertex[2], + vertexCasc[0], vertexCasc[1], vertexCasc[2], + vertexV0[0], vertexV0[1], vertexV0[2], + bachTrack.sign(), + covMatrixSV[0], covMatrixSV[1], covMatrixSV[2], covMatrixSV[3], covMatrixSV[4], covMatrixSV[5], + pVecCharmBaryon[0], pVecCharmBaryon[1], pVecCharmBaryon[2], + pVecCascAsD[0], pVecCascAsD[1], pVecCascAsD[2], + pVecCharmBachAsD[0], pVecCharmBachAsD[1], pVecCharmBachAsD[2], + pVecV0[0], pVecV0[1], pVecV0[2], + pVecBach[0], pVecBach[1], pVecBach[2], + pVecV0DauPos[0], pVecV0DauPos[1], pVecV0DauPos[2], + pVecV0DauNeg[0], pVecV0DauNeg[1], pVecV0DauNeg[2], + impactParameterCasc.getY(), impactParameterCharmBach.getY(), + impactParameterCasc.getZ(), impactParameterCharmBach.getZ(), + std::sqrt(impactParameterCasc.getSigmaY2()), std::sqrt(impactParameterCharmBach.getSigmaY2()), + cascAodElement.v0Id(), v0AodElement.posTrackId(), v0AodElement.negTrackId(), + cand.cascadeId(), trackCharmBachelor.globalIndex(), cand.prong0Id(), + massLambda, massCasc, massCharmBaryonCand, + cpaV0, cpaCharmBaryon, cpaCasc, + cpaxyV0, cpaxyCharmBaryon, cpaxyCasc, + ctOmegac0, ctCasc, ctV0, + etaV0DauPos, etaV0DauNeg, etaBach, etaCharmBach, etaCharmBaryon, etaCasc, etaV0, + dcaxyV0DauPos, dcaxyV0DauNeg, dcaxyBach, + dcazV0DauPos, dcazV0DauNeg, dcazBach, + dcaCascDau, dcaV0Dau, dcaCharmBaryonDau, + decLenCharmBaryon, decLenCasc, decLenV0, errorDecayLength, errorDecayLengthXY); + } + } // candidate loop + } // end of run function + + // template function for running Charm Baryon reconstruction via KFParticle method + /// \brief centEstimator is for different centrality estimators + /// \brief decayChannel is for different decay channels. 0 for XiczeroOmegaczeroToXiPi, 1 for OmegaczeroToOmegaPi, 2 for OmegaczeroToOmegaK + /// \brief Colls is for collision tables joined with different centrality estimators + /// \brief Hist is for QA histograms + template + void runCreatorWithKfParticle(Colls const&, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const&, // -> Implemented for internal cascade building + aod::V0s const&, // -> Implemented for internal cascade building + TracksWCovDcaExtraPidPrPiKa const&, + aod::BCsWithTimestamps const&, + Hist& hInvMassCharmBaryon, + Hist& hCandCounter) + { + // Loop over candidates + for (auto const& cand : candidates) { + + // Fill cascandidates before selection + if (configs.fillHistograms) { + hCandCounter->Fill(All); + } + + // Apply hfflag selection + if (!TESTBIT(cand.hfflag(), decayChannel)) { + continue; + } else { + if (configs.fillHistograms) { + hCandCounter->Fill(HfFlagPass); + } + } + + // Event selection + auto collision = cand.collision_as(); + float centrality{-1.f}; + const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); + if (rejectionMask != 0) { // None of the event seletion satisfied -> Reject this candidate + continue; + } + + //------------------------------Set Magnetic field------------------------------ + auto bc = collision.template bc_as(); + if (runNumber != bc.runNumber()) { + LOG(info) << ">>>>>>>>>> Current run Number : " << runNumber; + initCCDB(bc, runNumber, ccdb, configs.isRun2 ? configs.ccdbPathGrp : configs.ccdbPathGrpMag, lut, configs.isRun2); + magneticField = o2::base::Propagator::Instance()->getNominalBz(); + LOG(info) << ">>>>>>>>>> Magnetic field: " << magneticField; + } + // magnetic field setting for KFParticle + straHelper.fitter.setBz(magneticField); // -> Manetic field setting for internal cascade building + KFParticle::SetField(magneticField); // -> Magnetic field setting for CharmBaryon building + + //------------------Intenal Cascade building------------------ + auto cascAodElement = cand.cascade_as(); + auto v0AodElement = cascAodElement.v0_as(); + auto posTrack = v0AodElement.posTrack_as(); + auto negTrack = v0AodElement.negTrack_as(); + auto bachTrack = cascAodElement.bachelor_as(); + + // Make cascade starting from V0 + // If success, fill Cascade and V0 information for reconstruction + if (!straHelper.buildCascadeCandidateWithKF(collision.globalIndex(), + collision.posX(), collision.posY(), collision.posZ(), + posTrack, + negTrack, + bachTrack, + false, // calculateBachelorBaryonVariables + configs.kfConstructMethod, + configs.kfTuneForOmega, + configs.kfUseV0MassConstraint, + configs.kfUseCascadeMassConstraint, + configs.kfDoDCAFitterPreMinimV0, + configs.kfDoDCAFitterPreMinimCasc)) { + LOG(info) << "This cascade cannot be rebuilt"; + continue; + } else { + float storeMass = (decayChannel == 0) ? straHelper.cascade.massXi : straHelper.cascade.massOmega; + float storePt = RecoDecay::pt(straHelper.cascade.cascadeMomentum); + registry.fill(HIST("hCascMass"), storeMass); + registry.fill(HIST("hCascPt"), storePt); + } + + //------------------------------Cascade pre-selection------------------------------ + float massCasc = (decayChannel == 0) ? straHelper.cascade.massXi : straHelper.cascade.massOmega; + + if (configs.kfDoCascadePreselection) { + + // pre selection of dcaXY + // FIXME does cascadeDCAxy represents dcaxytoPV? + if (std::abs(straHelper.cascade.cascadeDCAxy) > configs.dcaXYToPVCascadeMax) { + continue; + } + // pre selection on dcaV0Daughters + if (std::abs(straHelper.cascade.v0DaughterDCA) > configs.dcaV0DaughtersMax) { + continue; + } + // pre selection on dcaCascDaughters + if (std::abs(straHelper.cascade.cascadeDaughterDCA) > configs.dcaCascDaughtersMax) { + continue; + } + // pre selection on invariantmass + if (std::abs(massCasc - massOfCascades[decayChannel]) > configs.massToleranceCascade) { + continue; + } + } + + if (configs.fillHistograms) { + hCandCounter->Fill(CascReconstructed); + } + + //------------------------------Info of V0 and Cascade------------------------------ + // -> This quantities are used for physical properties of selected candidates + // -> Not used for candidate creation + // auto chargeCasc = casc.sign() > 0 ? 1 : -1; + // std::array vertexV0 = {casc.xlambda(), casc.ylambda(), casc.zlambda()}; + // std::array pVecV0 = {casc.pxlambda(), casc.pylambda(), casc.pzlambda()}; + // std::array vertexCasc = {casc.x(), casc.y(), casc.z()}; + // std::array pVecCasc = {casc.px(), casc.py(), casc.pz()}; + + std::array vertexV0 = {straHelper.cascade.v0Position[0], straHelper.cascade.v0Position[1], straHelper.cascade.v0Position[2]}; + std::array pVecV0 = {straHelper.cascade.v0Momentum[0], straHelper.cascade.v0Momentum[1], straHelper.cascade.v0Momentum[2]}; + std::array pVecV0DauPos = {straHelper.cascade.positiveMomentum[0], straHelper.cascade.positiveMomentum[1], straHelper.cascade.positiveMomentum[2]}; + std::array pVecV0DauNeg = {straHelper.cascade.negativeMomentum[0], straHelper.cascade.negativeMomentum[1], straHelper.cascade.negativeMomentum[2]}; + + int chargeCasc = straHelper.cascade.charge > 0 ? 1 : -1; + std::array vertexCasc = {straHelper.cascade.cascadePosition[0], straHelper.cascade.cascadePosition[1], straHelper.cascade.cascadePosition[2]}; + std::array pVecCasc = {straHelper.cascade.cascadeMomentum[0], straHelper.cascade.cascadeMomentum[1], straHelper.cascade.cascadeMomentum[2]}; + std::array covCasc = {0.}; + + //------------------------------Create Charm Baryon as KF Particle object------------------------------ + + // initialize primary vertex + KFPVertex kfpVertex = createKFPVertexFromCollision(collision); + float covMatrixPV[6]; + kfpVertex.GetCovarianceMatrix(covMatrixPV); + KFParticle kfPv(kfpVertex); // -> For calculation of DCAs to PV + + //~~~~~create charm bachelor track into KFParticle object~~~~~// + auto trackCharmBachelor = cand.prong0_as(); + KFPTrack kfpTrackCharmBachelor = createKFPTrackFromTrack(trackCharmBachelor); + KFParticle kfCharmBachelor(kfpTrackCharmBachelor, pdgOfCharmBach[decayChannel]); + + //~~~~~create Cascade as KFParticle object~~~~~// + // -> This process is essential in building Charm baryon via KFParticle package. + // Since we don't have any information given from other tables about Charm baryon, + // only way to construct Charm baryon is to construct Charm baryon from it's daughters. + constexpr std::size_t NElementsStateVectorCasc{6}; + std::array xyzpxpypzCasc = {vertexCasc[0], vertexCasc[1], vertexCasc[2], pVecCasc[0], pVecCasc[1], pVecCasc[2]}; + float parPosMomCasc[NElementsStateVectorCasc]; + std::copy(xyzpxpypzCasc.begin(), xyzpxpypzCasc.end(), parPosMomCasc); + + KFParticle kfCasc; + // kfCasc.Create(parPosMomCasc, casc.kfTrackCovMat(), casc.sign(), massCasc); + kfCasc.Create(parPosMomCasc, straHelper.cascade.covariance, straHelper.cascade.charge, massCasc); + if (configs.kfConstrainInvMassCasc) { + kfCasc.SetNonlinearMassConstraint(massOfCascades[decayChannel]); + } + + //~~~~~~Create Charm Baryon as KFParticle object~~~~~// + KFParticle kfCharmBaryon; + const KFParticle* kfDaughterCharmBaryon[2] = {&kfCharmBachelor, &kfCasc}; + kfCharmBaryon.SetConstructMethod(configs.kfConstructMethod); + try { + kfCharmBaryon.Construct(kfDaughterCharmBaryon, 2); + } catch (std::runtime_error& e) { + LOG(debug) << "Failed to construct Charm Baryon : " << e.what(); + continue; + } + + // Get covariance matrix of Charm Baryon + auto covMatrixCharmBaryon = kfCharmBaryon.CovarianceMatrix(); + + // transport Charm Baryon daughters to Charm Baryon to decay vertex + // FIXME Is this process necessary? + float secondaryVertex[3] = {kfCharmBaryon.GetX(), kfCharmBaryon.GetY(), kfCharmBaryon.GetZ()}; + kfCasc.TransportToPoint(secondaryVertex); + kfCharmBachelor.TransportToPoint(secondaryVertex); + + // Fill histogram + if (configs.fillHistograms) { + hCandCounter->Fill(VertexFit); + hInvMassCharmBaryon->Fill(kfCharmBaryon.GetMass()); + } + + // Transport cascade and charm baryon to decay vertex(just to be sure) + kfCharmBaryon.TransportToDecayVertex(); + kfCasc.TransportToDecayVertex(); + + //!~~~~~Extra calculations for QA~~~~~~!// + /// These quantities were calculated to fill in existing output tables & QA purpose + /// In the future, these calculation can be skipped? This must be dicussed. + + /// Extra calculation for V0 daughters + KFPTrack kfpTrackV0Pos = createKFPTrackFromTrack(posTrack); + KFPTrack kfpTrackV0Neg = createKFPTrackFromTrack(negTrack); + + int pdgV0Pos = (trackCharmBachelor.sign() > 0) ? kProton : kPiPlus; + int pdgV0Neg = (trackCharmBachelor.sign() > 0) ? kPiMinus : -kProton; + + KFParticle kfV0Pos(kfpTrackV0Pos, pdgV0Pos); + KFParticle kfV0Neg(kfpTrackV0Neg, pdgV0Neg); + + auto trackParCovV0Pos = getTrackParCovFromKFP(kfV0Pos, kfV0Pos.GetPDG(), 1); + auto trackParCovV0Neg = getTrackParCovFromKFP(kfV0Neg, kfV0Neg.GetPDG(), -1); + + /// Extra calculations for V0 + KFParticle kfV0; + float xyzpxpypzV0[6] = {vertexV0[0], vertexV0[1], vertexV0[2], pVecV0[0], pVecV0[1], pVecV0[2]}; + kfV0.Create(xyzpxpypzV0, straHelper.cascade.kfTrackCovarianceV0, 0, MassLambda0); + if (configs.kfConstrainInvMassV0) { + kfV0.SetNonlinearMassConstraint(MassLambda); + } + + float kfMassV0, kfSigMassV0; + kfV0.GetMass(kfMassV0, kfSigMassV0); + + KFParticle kfV0ConstrainedToPv = kfV0; + KFParticle kfV0ConstrainedToCasc = kfV0; + kfV0ConstrainedToPv.SetProductionVertex(kfPv); + kfV0ConstrainedToCasc.SetProductionVertex(kfCasc); + + /// Extra calculations for bachelor + KFPTrack kfpTrackBach = createKFPTrackFromTrack(bachTrack); + KFParticle kfBach(kfpTrackBach, pdgOfBach[decayChannel]); + + KFParticle kfBachConstrainedToPv = kfBach; + KFParticle kfBachConstrainedToCasc = kfBach; + kfBachConstrainedToPv.SetProductionVertex(kfPv); + kfBachConstrainedToCasc.SetProductionVertex(kfCasc); + + auto trackParCovBach = getTrackParCovFromKFP(kfBachConstrainedToCasc, kfBachConstrainedToCasc.GetPDG(), bachTrack.sign()); + + /// Extra calculations for casc + float kfMassCasc, kfSigMassCasc; + kfCasc.GetMass(kfMassCasc, kfSigMassCasc); + + KFParticle kfCascConstrainedToPv = kfCasc; + KFParticle kfCascConstrainedToCharmBaryon = kfCasc; + kfCascConstrainedToPv.SetProductionVertex(kfPv); + kfCascConstrainedToCharmBaryon.SetProductionVertex(kfCharmBaryon); + + /// Extra calculation for charm bachelor + KFParticle kfCharmBachelorConstrainedToPv = kfCharmBachelor; + KFParticle kfCharmBachelorConstrainedToCharmBaryon = kfCharmBachelor; + kfCharmBachelorConstrainedToPv.SetProductionVertex(kfPv); + kfCharmBachelorConstrainedToCharmBaryon.SetProductionVertex(kfCharmBaryon); + + /// Extra calculation for charm baryon + float kfMassCharmBaryon, kfSigMassCharmBaryon; + kfCharmBaryon.GetMass(kfMassCharmBaryon, kfSigMassCharmBaryon); + + KFParticle kfCharmBaryonConstrainedToPv = kfCharmBaryon; + kfCharmBaryonConstrainedToPv.SetProductionVertex(kfPv); + + //------------------------------Calculate physical quantities and fill candidate table------------------------------ + + // Get updated daughter tracks after vertex fit + o2::track::TrackParCov trackParCovCasc = getTrackParCovFromKFP(kfCascConstrainedToCharmBaryon, kfCascConstrainedToCharmBaryon.GetPDG(), kfBach.GetQ()); + trackParCovCasc.setAbsCharge(1); + o2::track::TrackParCov trackParCovCharmBachelor = getTrackParCovFromKFP(kfCharmBachelorConstrainedToCharmBaryon, kfCharmBachelorConstrainedToCharmBaryon.GetPDG(), -kfBach.GetQ()); + trackParCovCharmBachelor.setAbsCharge(1); + + // impact parameters + std::array impactParameterV0DauPos; + std::array impactParameterV0DauNeg; + std::array impactParameterBach; + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, trackParCovV0Pos, 2.f, matCorr, &impactParameterV0DauPos); + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, trackParCovV0Neg, 2.f, matCorr, &impactParameterV0DauNeg); + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, trackParCovBach, 2.f, matCorr, &impactParameterBach); + float dcaxyV0DauPos = impactParameterV0DauPos[0]; + float dcaxyV0DauNeg = impactParameterV0DauNeg[0]; + float dcaxyBach = impactParameterBach[0]; + float dcazV0DauPos = impactParameterV0DauPos[1]; + float dcazV0DauNeg = impactParameterV0DauNeg[1]; + float dcazBach = impactParameterBach[1]; + + o2::dataformats::DCA impactParameterCasc, impactParameterCharmBachelor; + auto primaryVertex = getPrimaryVertex(collision); + o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackParCovCasc, 2.f, matCorr, &impactParameterCasc); + o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackParCovCharmBachelor, 2.f, matCorr, &impactParameterCharmBachelor); + + // get Chi2Topo/NDF + float chi2NdfTopoV0ToPv = kfV0ConstrainedToPv.GetChi2() / kfV0ConstrainedToPv.GetNDF(); + float chi2NdfTopoV0ToCasc = kfV0ConstrainedToPv.GetChi2() / kfV0ConstrainedToPv.GetNDF(); + float chi2NdfTopoBachToPv = kfBachConstrainedToPv.GetChi2() / kfBachConstrainedToPv.GetNDF(); + float chi2NdfTopoBachToCasc = kfBachConstrainedToCasc.GetChi2() / kfBachConstrainedToCasc.GetNDF(); + float chi2NdfTopoCascToPv = kfCascConstrainedToPv.GetChi2() / kfCascConstrainedToPv.GetNDF(); + float chi2NdfTopoCascToCharmBaryon = kfCascConstrainedToCharmBaryon.GetChi2() / kfCascConstrainedToCharmBaryon.GetNDF(); + float chi2NdfTopoCharmBachelorToPv = kfCharmBachelorConstrainedToPv.GetChi2() / kfCharmBachelorConstrainedToPv.GetNDF(); + float chi2NdfTopoCharmBachelorToCharmBaryon = kfCharmBachelorConstrainedToCharmBaryon.GetChi2() / kfCharmBachelorConstrainedToCharmBaryon.GetNDF(); + float chi2NdfTopoCharmBaryonToPv = kfCharmBachelorConstrainedToPv.GetChi2() / kfCharmBachelorConstrainedToPv.GetNDF(); + + float deviationCharmBachelorToPv = kfCalculateChi2ToPrimaryVertex(kfCharmBachelor, kfPv); + float chi2NdfBach = kfBach.GetChi2() / kfBach.GetNDF(); + + // get ldl + float ldlV0 = ldlFromKF(kfV0, kfPv); + float ldlCasc = ldlFromKF(kfCasc, kfPv); + float ldlCharmBaryon = ldlFromKF(kfCharmBaryon, kfPv); + + // get DCAs + float kfDcaV0Daughters = kfV0Pos.GetDistanceFromParticle(kfV0Neg); + float kfDcaCascDaughters = kfBachConstrainedToCasc.GetDistanceFromParticle(kfV0ConstrainedToCasc); + float kfDcaCharmBaryonDaughters = kfCharmBachelorConstrainedToCharmBaryon.GetDistanceFromParticle(kfCascConstrainedToCharmBaryon); + float kfDcaXYCharmBachelorToPv = kfCharmBachelorConstrainedToCharmBaryon.GetDistanceFromVertexXY(kfPv); + float kfDcaXYCascToPv = kfCascConstrainedToCharmBaryon.GetDistanceFromVertexXY(kfPv); + + // get decay length - In XY + float decayLXYV0, errDecayLXYV0; + kfV0ConstrainedToCasc.GetDecayLengthXY(decayLXYV0, errDecayLXYV0); + + float decayLXYCasc, errDecayLXYCasc; + kfCascConstrainedToCharmBaryon.GetDecayLengthXY(decayLXYCasc, decayLXYCasc); + + float decayLXYCharmBaryon, errDecayLXYCharmBaryon; + kfCharmBaryonConstrainedToPv.GetDecayLengthXY(decayLXYCharmBaryon, errDecayLXYCharmBaryon); + + // get decay length - In XYZ + float decayLV0 = RecoDecay::distance(std::array{kfCasc.GetX(), kfCasc.GetY(), kfCasc.GetZ()}, std::array{kfV0.GetX(), kfV0.GetY(), kfV0.GetZ()}); + float decayLCasc = RecoDecay::distance(std::array{kfCharmBaryon.GetX(), kfCharmBaryon.GetY(), kfCharmBaryon.GetZ()}, std::array{kfCasc.GetX(), kfCasc.GetY(), kfCasc.GetZ()}); + float decayLCharmBaryon = RecoDecay::distance(std::array{collision.posX(), collision.posY(), collision.posZ()}, std::array{kfCharmBaryon.GetX(), kfCharmBaryon.GetY(), kfCharmBaryon.GetZ()}); + + double phiCharmBaryon, thetaCharmBaryon; + getPointDirection(std::array{kfV0.GetX(), kfV0.GetY(), kfV0.GetZ()}, std::array{kfCharmBaryon.GetX(), kfCharmBaryon.GetY(), kfCharmBaryon.GetZ()}, phiCharmBaryon, thetaCharmBaryon); + float errDecayLCharmBaryon = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phiCharmBaryon, thetaCharmBaryon) + getRotatedCovMatrixXX(covMatrixCharmBaryon, phiCharmBaryon, thetaCharmBaryon)); + + // get cosine of pointing angle + std::array pvCoord = {collision.posX(), collision.posY(), collision.posZ()}; + // float cosPaV0ToPv = casc.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); + float cosPaV0ToPv = RecoDecay::cpa(pvCoord, vertexV0, pVecV0); + float cosPaCascToPv = cpaFromKF(kfCasc, kfPv); + float cosPaCharmBaryonToPv = cpaFromKF(kfCharmBaryon, kfPv); + + float cosPaXYV0ToPv = RecoDecay::cpaXY(pvCoord, vertexV0, pVecV0); + float cosPaXYCascToPv = cpaXYFromKF(kfCasc, kfPv); + float cosPaXYCharmBaryonToPv = cpaXYFromKF(kfCharmBaryon, kfPv); + + float cosPaV0ToCasc = RecoDecay::cpa(vertexCasc, vertexV0, pVecV0); + float cosPaCascToCharmBaryon = cpaFromKF(kfCasc, kfCharmBaryon); + float cosPaXYV0ToCasc = RecoDecay::cpaXY(vertexCasc, vertexV0, pVecV0); + float cosPaXYCascToCharmBaryon = cpaXYFromKF(kfCasc, kfCharmBaryon); + + // KF pT, eta + float ptCasc = kfCascConstrainedToCharmBaryon.GetPt(); + float ptCharmBachelor = kfCharmBachelorConstrainedToCharmBaryon.GetPt(); + float ptCharmBaryon = kfCharmBaryon.GetPt(); + float yCharmBaryon = kfCharmBaryon.GetRapidity(); + + // get KF cosThetaStar(?) + float cosThetaStarCharmBachelorXic0, cosThetaStarCharmBachelorOmegac0; + if constexpr (decayChannel == hf_cand_casc_lf::DecayType2Prong::XiczeroOmegaczeroToXiPi) { + cosThetaStarCharmBachelorXic0 = cosThetaStarFromKF(0, pdgOfCharmBaryon[decayChannel], pdgOfCascade[decayChannel], pdgOfCharmBach[decayChannel], kfCascConstrainedToCharmBaryon, kfCharmBachelorConstrainedToCharmBaryon); + cosThetaStarCharmBachelorOmegac0 = cosThetaStarFromKF(0, pdgOfCharmBaryon[decayChannel + 1], pdgOfCascade[decayChannel], pdgOfCharmBach[decayChannel], kfCascConstrainedToCharmBaryon, kfCharmBachelorConstrainedToCharmBaryon); + } else if constexpr (decayChannel == hf_cand_casc_lf::DecayType2Prong::OmegaczeroToOmegaPi) { + cosThetaStarCharmBachelorOmegac0 = cosThetaStarFromKF(0, pdgOfCharmBaryon[decayChannel + 1], pdgOfCascade[decayChannel], pdgOfCharmBach[decayChannel], kfCascConstrainedToCharmBaryon, kfCharmBachelorConstrainedToCharmBaryon); + } else if constexpr (decayChannel == hf_cand_casc_lf::DecayType2Prong::OmegaczeroToOmegaK) { + cosThetaStarCharmBachelorXic0 = cosThetaStarFromKF(0, +kXiC0, pdgOfCascade[decayChannel], pdgOfCharmBach[decayChannel], kfCascConstrainedToCharmBaryon, kfCharmBachelorConstrainedToCharmBaryon); + cosThetaStarCharmBachelorOmegac0 = cosThetaStarFromKF(0, +kOmegaC0, pdgOfCascade[decayChannel], pdgOfCharmBach[decayChannel], kfCascConstrainedToCharmBaryon, kfCharmBachelorConstrainedToCharmBaryon); + } + + // KF ct + float ctV0 = kfV0ConstrainedToCasc.GetLifeTime(); + float ctCasc = kfCascConstrainedToCharmBaryon.GetLifeTime(); + float ctCharmBaryon = kfCharmBachelorConstrainedToPv.GetLifeTime(); + + //------------------------------Calculate physical quantities and fill candidate table------------------------------ + + //------------------------------Fill the table------------------------------ + if constexpr (decayChannel == hf_cand_casc_lf::DecayType2Prong::XiczeroOmegaczeroToXiPi) { + cursors.rowCandToXiPiKf(collision.globalIndex(), // Global index of collision + pvCoord[0], pvCoord[1], pvCoord[2], // coordination of PV + kfCasc.GetX(), kfCasc.GetY(), kfCasc.GetZ(), // Decay position of kfCasc + vertexV0[0], vertexV0[1], vertexV0[2], + bachTrack.sign(), + covMatrixCharmBaryon[0], covMatrixCharmBaryon[1], covMatrixCharmBaryon[2], covMatrixCharmBaryon[3], covMatrixCharmBaryon[4], covMatrixCharmBaryon[5], + kfCharmBaryon.GetPx(), kfCharmBaryon.GetPy(), kfCharmBaryon.GetPz(), // x, y, z momentum of charm baryon + kfCascConstrainedToCharmBaryon.GetPx(), kfCascConstrainedToCharmBaryon.GetPy(), kfCascConstrainedToCharmBaryon.GetPz(), + kfCharmBachelorConstrainedToCharmBaryon.GetPx(), kfCharmBachelorConstrainedToCharmBaryon.GetPy(), kfCharmBachelorConstrainedToCharmBaryon.GetPz(), + pVecV0[0], pVecV0[1], pVecV0[2], + kfBachConstrainedToCasc.GetPx(), kfBachConstrainedToCasc.GetPy(), kfBachConstrainedToCasc.GetPz(), + pVecV0DauPos[0], pVecV0DauPos[1], pVecV0DauPos[2], + pVecV0DauNeg[0], pVecV0DauNeg[1], pVecV0DauNeg[2], + cascAodElement.v0Id(), v0AodElement.posTrackId(), v0AodElement.negTrackId(), + cand.cascadeId(), trackCharmBachelor.globalIndex(), cand.prong0Id(), + kfMassV0, kfMassCasc, kfMassCharmBaryon, + cosPaV0ToPv, cosPaCascToPv, + // cosPaCharmBaryonToPv, cosPaXYV0ToPv, cosPaXYCharmBaryonToPv, cosPaXYCascToPv, + ctCasc, ctV0, ctCharmBaryon, + kfV0Pos.GetEta(), kfV0Neg.GetEta(), kfBach.GetEta(), kfCharmBachelorConstrainedToCharmBaryon.GetEta(), + kfCharmBaryon.GetEta(), kfCasc.GetEta(), kfV0.GetEta(), + dcaxyV0DauPos, dcaxyV0DauNeg, dcaxyBach, + kfDcaCascDaughters, kfDcaV0Daughters, kfDcaCharmBaryonDaughters, + kfDcaXYCharmBachelorToPv, kfDcaXYCascToPv, + kfV0.GetChi2(), kfCasc.GetChi2(), kfCharmBaryon.GetChi2(), + /*FIXME chi2 of mass constrained, V0 and casc*/ kfV0.GetChi2(), kfCasc.GetChi2(), + ldlV0, ldlCasc, // ldlCharmBaryon, + chi2NdfTopoV0ToPv, chi2NdfTopoCascToPv, chi2NdfTopoCharmBachelorToPv, chi2NdfTopoCharmBaryonToPv, + chi2NdfTopoV0ToCasc, chi2NdfTopoCascToCharmBaryon, + decayLXYV0, decayLXYCasc, decayLXYCharmBaryon, + cosPaV0ToCasc, cosPaCascToCharmBaryon, // cosPaXYV0ToCasc, cosPaXYCascToCharmBaryon, + yCharmBaryon, // ptCharmBachelor, ptCharmBaryon, + cosThetaStarCharmBachelorXic0, + kfV0.GetNDF(), kfCasc.GetNDF(), kfCharmBaryon.GetNDF(), /*FIXME chi2, NDF of mass constrained V0/casc*/ kfV0.GetNDF(), kfCasc.GetNDF(), + kfV0.GetChi2() / kfV0.GetNDF(), kfCasc.GetChi2() / kfCasc.GetNDF(), kfCharmBaryon.GetChi2() / kfCharmBaryon.GetNDF(), /*FIXME chi2, NDF of mass constrained V0/casc*/ kfV0.GetChi2() / kfV0.GetNDF(), kfV0.GetChi2() / kfV0.GetNDF()); + + } else if constexpr (decayChannel == hf_cand_casc_lf::DecayType2Prong::OmegaczeroToOmegaPi) { + cursors.rowCandToOmegaPi(collision.globalIndex(), + pvCoord[0], pvCoord[1], pvCoord[2], + /*vertexCharmBaryonFromFitter. For KF, this is 0*/ 0.f, 0.f, 0.f, + kfCasc.GetX(), kfCasc.GetY(), kfCasc.GetZ(), + vertexV0[0], vertexV0[1], vertexV0[2], + trackCharmBachelor.sign(), + covMatrixCharmBaryon[0], covMatrixCharmBaryon[1], covMatrixCharmBaryon[2], covMatrixCharmBaryon[3], covMatrixCharmBaryon[4], covMatrixCharmBaryon[5], + kfCharmBaryon.GetPx(), kfCharmBaryon.GetPy(), kfCharmBaryon.GetPz(), + kfCascConstrainedToCharmBaryon.GetPx(), kfCascConstrainedToCharmBaryon.GetPy(), kfCascConstrainedToCharmBaryon.GetPz(), + kfCharmBachelorConstrainedToCharmBaryon.GetPx(), kfCharmBachelorConstrainedToCharmBaryon.GetPy(), kfCharmBachelorConstrainedToCharmBaryon.GetPz(), + pVecV0[0], pVecV0[1], pVecV0[2], + kfBachConstrainedToCasc.GetPx(), kfBachConstrainedToCasc.GetPy(), kfBachConstrainedToCasc.GetPz(), + pVecV0DauPos[0], pVecV0DauPos[1], pVecV0DauPos[2], + pVecV0DauNeg[0], pVecV0DauNeg[1], pVecV0DauNeg[2], + impactParameterCasc.getY(), impactParameterCharmBachelor.getY(), + impactParameterCasc.getZ(), impactParameterCharmBachelor.getZ(), + std::sqrt(impactParameterCasc.getSigmaY2()), std::sqrt(impactParameterCharmBachelor.getSigmaY2()), + cascAodElement.v0Id(), v0AodElement.posTrackId(), v0AodElement.negTrackId(), + cand.cascadeId(), trackCharmBachelor.globalIndex(), cand.prong0Id(), + kfMassV0, kfMassCasc, kfMassCharmBaryon, + cosPaV0ToPv, cosPaCharmBaryonToPv, cosPaCascToPv, cosPaXYV0ToPv, cosPaXYCharmBaryonToPv, cosPaXYCascToPv, + ctCharmBaryon, ctCasc, ctV0, + kfV0Pos.GetEta(), kfV0Neg.GetEta(), kfBach.GetEta(), kfCharmBachelorConstrainedToCharmBaryon.GetEta(), + kfCharmBaryon.GetEta(), kfCasc.GetEta(), kfV0.GetEta(), + dcaxyV0DauPos, dcaxyV0DauNeg, dcaxyBach, + dcazV0DauPos, dcazV0DauNeg, dcazBach, + kfDcaCascDaughters, straHelper.cascade.v0DaughterDCA, kfDcaCharmBaryonDaughters, + decayLCharmBaryon, decayLCasc, decayLV0, errDecayLCharmBaryon, errDecayLXYCharmBaryon, cand.hfflag()); + + cursors.rowCandToOmegaPiKf(kfDcaXYCharmBachelorToPv, kfDcaXYCascToPv, + /*V0 chi2. Taken from LF*/ kfV0.GetChi2(), kfCasc.GetChi2(), kfCharmBaryon.GetChi2(), /*Mass constraint only done when requested*/ kfV0.GetChi2(), /*Mass constraint only done when requested*/ kfCasc.GetChi2(), + ldlV0, ldlCasc, ldlCharmBaryon, + chi2NdfTopoV0ToPv, chi2NdfTopoCascToPv, chi2NdfTopoCharmBachelorToPv, chi2NdfTopoCharmBaryonToPv, deviationCharmBachelorToPv, + chi2NdfTopoV0ToCasc, chi2NdfTopoCascToCharmBaryon, + decayLXYV0, decayLXYCasc, decayLXYCharmBaryon, + cosPaV0ToCasc, cosPaCascToCharmBaryon, cosPaXYV0ToCasc, cosPaXYCascToCharmBaryon, + kfCharmBaryon.GetRapidity(), ptCharmBachelor, ptCharmBaryon, + cosThetaStarCharmBachelorOmegac0, + kfV0.GetNDF(), kfCasc.GetNDF(), kfCharmBaryon.GetNDF(), /*Mass constraint only done when requested*/ kfV0.GetNDF(), /*Mass constraint only done when requested*/ kfCasc.GetNDF(), + kfV0.GetChi2() / kfV0.GetNDF(), kfCasc.Chi2() / kfCasc.GetNDF(), kfCharmBaryon.GetChi2() / kfCharmBaryon.GetNDF(), + /*Mass constraint only done when requested*/ kfV0.GetChi2() / kfV0.GetNDF(), kfCasc.Chi2() / kfCasc.GetNDF(), + /*FIXME casc-rej not calculated. For now, fill in mass of KFCasc*/ kfCasc.GetMass()); + + } else { + cursors.rowCandToOmegaKaKf(collision.globalIndex(), + collision.posX(), collision.posY(), collision.posZ(), + kfPv.GetX(), kfPv.GetY(), kfPv.GetZ(), + vertexV0[0], vertexV0[1], vertexV0[2], + pVecV0[0], pVecV0[1], pVecV0[2], + vertexCasc[0], vertexCasc[1], vertexCasc[2], + pVecCasc[0], pVecCasc[1], pVecCasc[2], + vertexV0[0], vertexV0[1], vertexV0[2], + pVecV0[0], pVecV0[1], pVecV0[2], + kfCasc.GetX(), kfCasc.GetY(), kfCasc.GetZ(), + kfCasc.GetPx(), kfCasc.GetPy(), kfCasc.GetPz(), + kfCharmBaryon.GetX(), kfCharmBaryon.GetY(), kfCharmBaryon.GetZ(), + kfCharmBaryon.GetPx(), kfCharmBaryon.GetPy(), kfCharmBaryon.GetPz(), + straHelper.cascade.charge, + kfV0Pos.GetEta(), kfV0Neg.GetEta(), kfBach.GetEta(), kfCharmBachelor.GetEta(), kfV0.GetEta(), kfCasc.GetEta(), kfCharmBaryon.GetEta(), kfCharmBaryon.GetRapidity(), + impactParameterCharmBachelor.getY(), std::sqrt(impactParameterCharmBachelor.getSigmaY2()), impactParameterCasc.getY(), std::sqrt(impactParameterCasc.getSigmaY2()), + kfDcaV0Daughters, kfDcaCascDaughters, kfDcaCharmBaryonDaughters, + cosPaV0ToPv, cosPaCascToPv, cosPaCharmBaryonToPv, cosPaXYV0ToPv, cosPaXYCascToPv, cosPaXYCharmBaryonToPv, cosPaV0ToCasc, cosPaCascToCharmBaryon, cosPaXYV0ToCasc, cosPaXYCascToCharmBaryon, + kfV0.GetChi2() / kfV0.GetNDF(), kfCasc.GetChi2() / kfCasc.GetNDF(), kfCharmBaryon.GetChi2() / kfCharmBaryon.GetNDF(), + /*FIXME mass constraint is optional*/ kfV0.GetChi2() / kfV0.GetNDF(), kfCasc.GetChi2() / kfCasc.GetNDF(), + chi2NdfTopoV0ToCasc, chi2NdfTopoBachToCasc, chi2NdfTopoCharmBachelorToCharmBaryon, chi2NdfTopoCascToCharmBaryon, // Topological constraints to mother + chi2NdfTopoV0ToPv, chi2NdfTopoCascToPv, chi2NdfTopoCharmBachelorToPv, chi2NdfTopoCharmBaryonToPv, // Topological constraints to PV + ldlV0, ldlCasc, ldlCharmBaryon, + decayLXYV0, decayLXYCasc, decayLXYCharmBaryon, + kfMassV0, kfSigMassV0, kfMassCasc, kfSigMassCasc, /*FIXME no -rej has been made*/ kfMassCasc, /*FIXME no -rej has been made*/ kfSigMassCasc, kfMassCharmBaryon, kfSigMassCharmBaryon, + ptCharmBaryon, ptCharmBachelor, ptCasc, + cosThetaStarCharmBachelorOmegac0, cosThetaStarCharmBachelorXic0, ctV0, ctCasc, ctCharmBaryon, + cascAodElement.v0Id(), v0AodElement.posTrackId(), v0AodElement.negTrackId(), cand.cascadeId(), cand.prong0Id(), trackCharmBachelor.globalIndex()); + } + } // end candidate loop + } // end of runCreator + + ///////////////////////////////////////////////////// + /// /// + /// Process functions with DCAFitter /// + /// /// + ///////////////////////////////////////////////////// + + /*~~~~~~~~~~~~~~*/ + /*~~~To Xi Pi~~~*/ + /*~~~~~~~~~~~~~~*/ + void processToXiPiWithDCAFitterNoCent(SelectedCollisions const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithDCAFitter(collisions, candidates, cascades, v0s, tracks, bcsWithTimestamps, hInvMassCharmBaryonToXiPi, hCandidateCounterToXiPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToXiPiWithDCAFitterNoCent, "Charm candidte reconstruction with Xi Pi via DcaFitter method, no centrality", true); + +#if 0 + void processToXiPiWithDCAFitterNoCentWithTrackedCasc(SelectedCollisions const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::TrackedCascades const& cascades, + aod::V0s const& v0s, + TrackedCascLinked const& trackedCascLinked, + TracksWCovDcaExtraPidPrPiKa const& tracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithDCAFitter(collisions, candidates, trackedCascFull, trackedCascLinked, tracks, bcsWithTimestamps, hInvMassCharmBaryonToXiPi, hCandidateCounterToXiPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToXiPiWithDCAFitterNoCentWithTrackedCasc, "Charm candidte reconstruction with Xi Pi via DcaFitter method with tracked cascade, no centrality", false); +#endif + + void processToXiPiWithDCAFitterCentFT0C(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithDCAFitter(collisions, candidates, cascades, v0s, tracks, bcsWithTimestamps, hInvMassCharmBaryonToXiPi, hCandidateCounterToXiPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToXiPiWithDCAFitterCentFT0C, "Charm candidate reconstruction with Xi Pi via DcaFitter method, centrality selection on FT0C", false); + + void processToXiPiWithDCAFitterCentFT0M(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithDCAFitter(collisions, candidates, cascades, v0s, tracks, bcsWithTimestamps, hInvMassCharmBaryonToXiPi, hCandidateCounterToXiPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToXiPiWithDCAFitterCentFT0M, "Charm candidate reconstruction with Xi Pi via DcaFitter method, centrality selection on FT0M", false); + + /*~~~~~~~~~~~~~~~~~*/ + /*~~~To Omega Pi~~~*/ + /*~~~~~~~~~~~~~~~~~*/ + void processToOmegaPiWithDCAFitterNoCent(SelectedCollisions const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithDCAFitter(collisions, candidates, cascades, v0s, tracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaPi, hCandidateCounterToOmegaPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaPiWithDCAFitterNoCent, "Charm candidte reconstruction with Omega Pi via DcaFitter method, no centrality", false); + + void processToOmegaPiWithDCAFitterCentFT0C(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithDCAFitter(collisions, candidates, cascades, v0s, tracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaPi, hCandidateCounterToOmegaPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaPiWithDCAFitterCentFT0C, "Charm candidate reconstruction with Omega Pi via DcaFitter method, centrality selection on FT0C", false); + + void processToOmegaPiWithDCAFitterCentFT0M(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithDCAFitter(collisions, candidates, cascades, v0s, tracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaPi, hInvMassCharmBaryonToOmegaPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaPiWithDCAFitterCentFT0M, "Charm candidate reconstruction with Omega Pi via DcaFitter method, centrality selection on FT0M", false); + + /*~~~~~~~~~~~~~~~~~*/ + /*~~~To Omega Ka~~~*/ + /*~~~~~~~~~~~~~~~~~*/ + void processToOmegaKaWithDCAFitterNoCent(SelectedCollisions const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithDCAFitter(collisions, candidates, cascades, v0s, tracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaKa, hCandidateCounterToOmegaKa); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaKaWithDCAFitterNoCent, "Charm candidte reconstruction with Omega Ka via DcaFitter method, no centrality", false); + + void processToOmegaKaWithDCAFitterCentFT0C(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithDCAFitter(collisions, candidates, cascades, v0s, tracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaKa, hCandidateCounterToOmegaKa); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaKaWithDCAFitterCentFT0C, "Charm candidate reconstruction with Omega Ka via DcaFitter method, centrality selection on FT0C", false); + + void processToOmegaKaWithDCAFitterCentFT0M(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithDCAFitter(collisions, candidates, cascades, v0s, tracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaKa, hCandidateCounterToOmegaPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaKaWithDCAFitterCentFT0M, "Charm candidate reconstruction with Omega Ka via DcaFitter method, centrality selection on FT0M", false); + + ///////////////////////////////////////////////////// + /// /// + /// Process functions with KFParticle /// + /// /// + ///////////////////////////////////////////////////// + + /*~~~~~~~~~~~~~~*/ + /*~~~To Xi Pi~~~*/ + /*~~~~~~~~~~~~~~*/ + void processToXiPiWithKFParticleNoCent(SelectedCollisions const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithKfParticle(collisions, candidates, cascades, v0s, tracks, bcsWithTimestamps, hInvMassCharmBaryonToXiPi, hCandidateCounterToXiPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToXiPiWithKFParticleNoCent, "Charm Baryon decaying to Xi Pi reconstruction via KFParticle method, no centrality", false); + + void processToXiPiWithKFParticleCentFT0C(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithKfParticle(collisions, candidates, cascades, v0s, tracks, bcsWithTimestamps, hInvMassCharmBaryonToXiPi, hCandidateCounterToXiPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToXiPiWithKFParticleCentFT0C, "Charm Baryon decaying to Xi Pi reconstruction via KFParticle method, centrality on FT0C", false); + + void processToXiPiWithKFParticleCentFT0M(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithKfParticle(collisions, candidates, cascades, v0s, tracks, bcsWithTimestamps, hInvMassCharmBaryonToXiPi, hCandidateCounterToXiPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToXiPiWithKFParticleCentFT0M, "Charm Baryon decaying to Xi Pireconstruction via KFParticle method, centrality on FT0M", false); + + /*~~~~~~~~~~~~~~~~~*/ + /*~~~To Omega Pi~~~*/ + /*~~~~~~~~~~~~~~~~~*/ + void processToOmegaPiWithKFParticleNoCent(SelectedCollisions const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithKfParticle(collisions, candidates, cascades, v0s, tracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaPi, hCandidateCounterToOmegaPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaPiWithKFParticleNoCent, "Charm Baryon decaying to Omega Pi reconstruction via KFParticle method, no centrality", false); + + void processToOmegaPiWithKFParticleCentFT0C(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithKfParticle(collisions, candidates, cascades, v0s, tracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaPi, hCandidateCounterToOmegaPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaPiWithKFParticleCentFT0C, "Charm Baryon decaying to Omega Pi reconstruction via KFParticle method, centrality on FT0C", false); + + void processToOmegaPiWithKFParticleCentFT0M(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithKfParticle(collisions, candidates, cascades, v0s, tracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaPi, hCandidateCounterToOmegaPi); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaPiWithKFParticleCentFT0M, "Charm Baryong decaying to Omega Pi reconstruction via KFParticle method, centrality on FT0M", false); + + /*~~~~~~~~~~~~~~~~~*/ + /*~~~To Omega Ka~~~*/ + /*~~~~~~~~~~~~~~~~~*/ + void processToOmegaKaWithKFParticleNoCent(SelectedCollisions const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithKfParticle(collisions, candidates, cascades, v0s, tracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaKa, hCandidateCounterToOmegaKa); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaKaWithKFParticleNoCent, "Charm Baryon decaying to Omega Ka reconstruction via KFParticle method, no centrality", false); + + void processToOmegaKaWithKFParticleCentFT0C(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithKfParticle(collisions, candidates, cascades, v0s, tracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaKa, hCandidateCounterToOmegaKa); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaKaWithKFParticleCentFT0C, "Charm Baryon decaying to Omega Ka reconstruction via KFParticle method, centrality on FT0C", false); + + void processToOmegaKaWithKFParticleCentFT0M(soa::Join const& collisions, + aod::HfCascLf2Prongs const& candidates, + aod::Cascades const& cascades, + aod::V0s const& v0s, + TracksWCovDcaExtraPidPrPiKa const& tracks, + aod::BCsWithTimestamps const& bcsWithTimestamps) + { + runCreatorWithKfParticle(collisions, candidates, cascades, v0s, tracks, bcsWithTimestamps, hInvMassCharmBaryonToOmegaKa, hCandidateCounterToOmegaKa); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processToOmegaKaWithKFParticleCentFT0M, "Charm Baryong decaying to Omega Ka reconstruction via KFParticle method, centrality on FT0M", false); + + /////////////////////////////////////////////////////////////// + /// /// + /// Process functions for Collision monitoring /// + /// /// + /////////////////////////////////////////////////////////////// + + void processCollisionsNoCent(soa::Join const& collisions, + aod::BCsWithTimestamps const&) + { + for (const auto& collision : collisions) { + + // bitmask with event selection info + float centrality{-1.f}; + float occupancy = getOccupancyColl(collision, OccupancyEstimator::Its); + const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); + + // monitor the satisfied event selection + hfEvSel.fillHistograms(collision, rejectionMask, centrality, occupancy); + } + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processCollisionsNoCent, "Collision monitoring - No Centrality", true); + + void processCollisionsCentFT0C(soa::Join const& collisions, + aod::BCsWithTimestamps const&) + { + for (const auto& collision : collisions) { + + // bitmask with event selection info + float centrality{-1.f}; + float occupancy = getOccupancyColl(collision, OccupancyEstimator::Its); + const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); + + // monitor the satisfied event selection + hfEvSel.fillHistograms(collision, rejectionMask, centrality, occupancy); + } + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processCollisionsCentFT0C, "Collision monitoring - Centrality selection with FT0C", false); + + void processCollisionsCentFT0M(soa::Join const& collisions, + aod::BCsWithTimestamps const&) + { + for (const auto& collision : collisions) { + + // bitmask with event selection info + float centrality{-1.f}; + float occupancy = getOccupancyColl(collision, OccupancyEstimator::Its); + const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); + + // monitor the satisfied event selection + hfEvSel.fillHistograms(collision, rejectionMask, centrality, occupancy); + } + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0Qa, processCollisionsCentFT0M, "Collision monitoring - Centrality selection with FT0M", false); +}; + +struct HfCandidateCreatorXic0Omegac0QaMc { + + struct : ProducesGroup { + + Produces rowMcMatchRecXicToXiPi; + Produces rowMcMatchGenXicToXiPi; + Produces rowMcMatchRecOmegacToXiPi; + Produces rowMcMatchGenOmegacToXiPi; + Produces rowMcMatchRecToOmegaPi; + Produces rowMcMatchGenToOmegaPi; + Produces rowMcMatchRecToOmegaKa; + Produces rowMcMatchGenToOmegaKa; + + } cursors; + + struct : ConfigurableGroup { + + Configurable rejectBackground{"rejectBackground", true, "Reject particles from background events"}; // -> Used for only Xic0 + Configurable acceptTrackInteractionWithMaterial{"acceptTrackInteractionWithMaterial", false, "Accept candidates with final daughters interacting with materials"}; + Configurable fillMcHistograms{"fillMcHistograms", true, "Fill validation plots"}; + Configurable fillResidualTable{"fillResidualTable", true, "Fill table contaning residuals and pulls of PV and SV"}; + // Configurable matchDecayedPions{"matchedDecayedPions", true, "Match also candidates with daughter pion tracks that decay with kinked toploogy"}; + + } configs; + + enum McMatchFlag : uint8_t { + None = 0, + CharmBaryonUnmatched, + CascUnmatched, + V0Unmatched, + NumberMcMatchFlag + }; + + std::array pdgOfCharmBaryon{+kXiC0, +kOmegaC0, +kOmegaC0, +kOmegaC0}; + std::array pdgOfCascade{+kXiMinus, +kXiMinus, +kOmegaMinus, +kOmegaMinus}; + std::array pdgOfCharmBachelor{+kPiPlus, +kPiPlus, +kPiPlus, +kKPlus}; + std::array pdgOfBachelor{+kPiMinus, +kPiMinus, +kKMinus, +kKMinus}; + + // Table aliases + using McCollisionsNoCents = soa::Join; + using McCollisionsFT0Cs = soa::Join; + using McCollisionsFT0Ms = soa::Join; + using McCollisionsCentFT0Ms = soa::Join; // -> Used for subscription for process functions of centrality with FT0Ms + using BCsInfo = soa::Join; + + Preslice mcParticlesPerMcCollision = aod::mcparticle::mcCollisionId; + PresliceUnsorted colPerMcCollision = aod::mccollisionlabel::mcCollisionId; // -> Why use unsorted?? + PresliceUnsorted colPerMcCollisionFT0C = aod::mccollisionlabel::mcCollisionId; + PresliceUnsorted colPerMcCollisionFT0M = aod::mccollisionlabel::mcCollisionId; + + HistogramRegistry registry{"registry"}; + HfEventSelectionMc hfEvSelMc; + + void init(InitContext& initContext) + { + const auto& workflows = initContext.services().get(); + for (const DeviceSpec& device : workflows.devices) { + if (device.name.compare("hf-candidate-creator-xic0-omegac0-qa") == 0) { + hfEvSelMc.init(device, registry); + break; + } + } + + // Add histograms for QA + if (configs.fillMcHistograms) { + registry.add("hDebugStatusRec", "hDebugStatusRec", {HistType::kTH1D, {{NumberMcMatchFlag, 0, NumberMcMatchFlag}}}); + registry.add("hDebugStatusGen", "hDebugStatusGen", {HistType::kTH1D, {{NumberMcMatchFlag, 0, NumberMcMatchFlag}}}); + TString labels[McMatchFlag::NumberMcMatchFlag]; + labels[McMatchFlag::None] = "None"; + labels[McMatchFlag::CharmBaryonUnmatched] = "CharmBaryonUnmatched"; + labels[McMatchFlag::CascUnmatched] = "CascUnmatched"; + labels[McMatchFlag::V0Unmatched] = "V0Unmatched"; + for (int iBin = 0; iBin < McMatchFlag::NumberMcMatchFlag; iBin++) { + registry.get(HIST("hDebugStatusRec"))->GetXaxis()->SetBinLabel(iBin + 1, labels[iBin]); + registry.get(HIST("hDebugStatusGen"))->GetXaxis()->SetBinLabel(iBin + 1, labels[iBin]); + } + } + } + + // Helper function to fill in table for MC Rec + ///@brief decay Channel is the chosen reconstruction channel + ///@brief flag + ///@brief debug + ///@brief origin + ///@brief collisionMatched + ///@brief pt + ///@brief pdgcode + template + void fillRecoMcTableByDecayChannel(int8_t flag, int8_t debug, int8_t origin, bool collisionMatched, float pt, int pdgCode) + { + if constexpr (decayChannel == aod::hf_cand_xic0_omegac0::DecayType::XiczeroToXiPi) { + cursors.rowMcMatchRecXicToXiPi(flag, debug, origin, collisionMatched, pt, pdgCode); + } else if constexpr (decayChannel == aod::hf_cand_xic0_omegac0::DecayType::OmegaczeroToXiPi) { + cursors.rowMcMatchRecOmegacToXiPi(flag, debug, origin, collisionMatched, pt, pdgCode); + } else if constexpr (decayChannel == aod::hf_cand_xic0_omegac0::DecayType::OmegaczeroToOmegaPi) { + cursors.rowMcMatchRecToOmegaPi(flag, debug, origin, collisionMatched, pt, pdgCode); + } else { + cursors.rowMcMatchRecToOmegaKa(flag, debug, origin, collisionMatched, pt, pdgCode); + } + } + + // Helper function to fill in table for MC Gen + ///@brief decay Channel is the chosen reconstruction channel + ///@brief flag + ///@brief debug + ///@brief origin + ///@brief collisionMatched + ///@brief pt + ///@brief pdgcode + template + void fillGenMcTableByDecayChannel(int8_t flag, int8_t debugGenCharmBaryon, int8_t debugGenCascade, int8_t debugGenLambda, float ptCharmBaryon, float yCharmBaryon, int8_t origin, int idxMother) + { + if constexpr (decayChannel == aod::hf_cand_xic0_omegac0::DecayType::XiczeroToXiPi) { + cursors.rowMcMatchGenXicToXiPi(flag, debugGenCharmBaryon, debugGenCascade, debugGenLambda, ptCharmBaryon, yCharmBaryon, origin, idxMother); + } else if constexpr (decayChannel == aod::hf_cand_xic0_omegac0::DecayType::OmegaczeroToXiPi) { + cursors.rowMcMatchGenOmegacToXiPi(flag, debugGenCharmBaryon, debugGenCascade, debugGenLambda, ptCharmBaryon, yCharmBaryon, origin, idxMother); + } else if constexpr (decayChannel == aod::hf_cand_xic0_omegac0::DecayType::OmegaczeroToOmegaPi) { + cursors.rowMcMatchGenToOmegaPi(flag, debugGenCharmBaryon, debugGenCascade, debugGenLambda, ptCharmBaryon, yCharmBaryon, origin, idxMother); + } else { + cursors.rowMcMatchGenToOmegaKa(flag, debugGenCharmBaryon, debugGenCascade, debugGenLambda, ptCharmBaryon, yCharmBaryon, origin, idxMother); + } + } + + template + void runXic0Omegac0Mc(TRecoCand const& candidates, + TracksWMc const& tracks, + aod::McParticles const& mcParticles, + Colls const& collsWithMcLabels, + McCollisions const& mcCollisions, + BCsInfo const&) + { + int indexRec{-1}; + int indexRecCharmBaryon{-1}; + int8_t sign{-9}; + int8_t signCasc{-9}; + int8_t signV0{-9}; + int8_t flag{0}; + int8_t origin{0}; + int8_t debug{0}; + int8_t debugGenCharmBaryon{0}; + int8_t debugGenCasc{0}; + int8_t debugGenV0{0}; + bool collisionMatched = false; + float ptCharmBaryonGen = -999.; + float yCharmBaryonGen = -999.; + + //////////////////////////////////// + // Match reconstructed candidates // + //////////////////////////////////// + + for (const auto& candidate : candidates) { + + flag = 0; + origin = RecoDecay::OriginType::None; + debug = McMatchFlag::None; + collisionMatched = false; + std::vector idxBhadMothers{}; + + auto arrayDaughters = std::array{candidate.template bachelorFromCharmBaryon_as(), + candidate.template bachelor_as(), + candidate.template posTrack_as(), + candidate.template negTrack_as()}; + auto arrayDaughtersCasc = std::array{candidate.template bachelor_as(), + candidate.template posTrack_as(), + candidate.template negTrack_as()}; + auto arrayDaughtersV0 = std::array{candidate.template posTrack_as(), + candidate.template negTrack_as()}; + + // Reject particles from background events + if (configs.rejectBackground) { + bool fromBkg{false}; + for (auto const& daughter : arrayDaughters) { + if (daughter.has_mcParticle()) { + auto mcParticle = daughter.mcParticle(); + if (mcParticle.fromBackgroundEvent()) { + fromBkg = true; + break; + } + } + } + if (fromBkg) { + // fill the tables + fillRecoMcTableByDecayChannel(flag, debug, origin, collisionMatched, -1.f, 0); + continue; + } + } + + // CharmBaryon -> Charm bachelor + Cascade + indexRec = RecoDecay::getMatchedMCRec(mcParticles, arrayDaughters, pdgOfCharmBaryon[decayChannel], std::array{pdgOfCharmBachelor[decayChannel], pdgOfBachelor[decayChannel], +kProton, +kPiMinus}, true, &sign, 3); + indexRecCharmBaryon = indexRec; + if (indexRec == -1) { // Xic0 not reconstructed + debug = McMatchFlag::CharmBaryonUnmatched; + } + if (indexRec > -1) { + // Cascade -> Bachelor + V0 + indexRec = RecoDecay::getMatchedMCRec(mcParticles, arrayDaughtersCasc, pdgOfCascade[decayChannel], std::array{pdgOfBachelor[decayChannel], +kProton, +kPiMinus}, true, &signCasc, 2); + if (indexRec == -1) { // Xi- not reconstructed + debug = McMatchFlag::CascUnmatched; + } + if (indexRec > -1) { + // V0 -> Pos + Neg + indexRec = RecoDecay::getMatchedMCRec(mcParticles, arrayDaughtersV0, +kLambda0, std::array{+kProton, +kPiMinus}, true, &signV0, 1); + if (indexRec == -1) { // V0 not reconstructed + debug = McMatchFlag::V0Unmatched; + } + if (indexRec > -1) { + flag = sign * (1 << decayChannel); + collisionMatched = candidate.template collision_as().mcCollisionId() == mcParticles.iteratorAt(indexRecCharmBaryon).mcCollisionId(); + } + } + } + + // Check if Xic0 is from b-hadron decay(prompt vs non-prompt) + if (flag != 0) { + auto particle = mcParticles.rawIteratorAt(indexRecCharmBaryon); + origin = RecoDecay::getCharmHadronOrigin(mcParticles, particle, false, &idxBhadMothers); + } + if (origin == RecoDecay::OriginType::NonPrompt) { + auto bHadMother = mcParticles.rawIteratorAt(idxBhadMothers[0]); + fillRecoMcTableByDecayChannel(flag, debug, origin, collisionMatched, bHadMother.pt(), bHadMother.pdgCode()); + } else { + fillRecoMcTableByDecayChannel(flag, debug, origin, collisionMatched, -1.f, 0); + } + } // candidate loop + + /////////////////////////////// + // Match generated particles // + /////////////////////////////// + + for (auto const& mcCollision : mcCollisions) { + + auto const mcParticlesPerMcColl = mcParticles.sliceBy(mcParticlesPerMcCollision, mcCollision.globalIndex()); + + float centrality{-1.f}; + uint16_t rejectionMask{0}; + int nSplitColl{0}; + + if constexpr (centEstimator == o2::hf_centrality::CentralityEstimator::None) { + const auto collSlice = collsWithMcLabels.sliceBy(colPerMcCollision, mcCollision.globalIndex()); + rejectionMask = hfEvSelMc.getHfMcCollisionRejectionMask(mcCollision, collSlice, centrality); + } else if constexpr (centEstimator == o2::hf_centrality::CentralityEstimator::FT0C) { + const auto collSlice = collsWithMcLabels.sliceBy(colPerMcCollisionFT0C, mcCollision.globalIndex()); + rejectionMask = hfEvSelMc.getHfMcCollisionRejectionMask(mcCollision, collSlice, centrality); + } else if constexpr (centEstimator == o2::hf_centrality::CentralityEstimator::FT0M) { + const auto collSlice = collsWithMcLabels.sliceBy(colPerMcCollisionFT0M, mcCollision.globalIndex()); + rejectionMask = hfEvSelMc.getHfMcCollisionRejectionMask(mcCollision, collSlice, centrality); + } + + hfEvSelMc.fillHistograms(mcCollision, rejectionMask); + + if (rejectionMask != 0) { // none of the event selection was satisfied(?) -> Reject all particles from this event + for (unsigned int i = 0; i < mcParticlesPerMcColl.size(); ++i) { + fillGenMcTableByDecayChannel(0, 0, 0, 0, -999., -999., RecoDecay::OriginType::None, -1); + } + continue; + } + + // Match generated particles + for (auto const& particle : mcParticlesPerMcColl) { + ptCharmBaryonGen = -999.; + yCharmBaryonGen = -999.; + flag = 0; + sign = 0; + debugGenCharmBaryon = 0; + debugGenCasc = 0; + debugGenV0 = 0; + origin = RecoDecay::OriginType::None; + std::vector idxBhadMothers{}; + float kYCutTight = 0.5; + float kYCutLoose = 0.8; + + // Reject particles from background events + if (particle.fromBackgroundEvent() && configs.rejectBackground) { + fillGenMcTableByDecayChannel(flag, debugGenCharmBaryon, debugGenCasc, debugGenV0, ptCharmBaryonGen, yCharmBaryonGen, origin, -1); + continue; + } + + // Charm Baryon -> Cascade + Charm bachelor + if (RecoDecay::isMatchedMCGen(mcParticles, particle, pdgOfCharmBaryon[decayChannel], std::array{pdgOfCascade[decayChannel], pdgOfCharmBachelor[decayChannel]}, true, &sign)) { + debugGenCharmBaryon = 1; + ptCharmBaryonGen = particle.pt(); + yCharmBaryonGen = particle.y(); + debug = 1; // -> Matched Xic0 + + for (auto const& daughterCharm : particle.template daughters_as()) { + if (std::abs(daughterCharm.pdgCode()) != pdgOfCascade[decayChannel]) { + continue; + } + // Xi -> Lambda + pi + if (RecoDecay::isMatchedMCGen(mcParticles, daughterCharm, pdgOfCascade[decayChannel], std::array{+kLambda0, pdgOfBachelor[decayChannel]}, true)) { + debugGenCasc = 1; // -> Matched Xi- + for (auto const& daughterCascade : daughterCharm.template daughters_as()) { + if (std::abs(daughterCascade.pdgCode()) != +kLambda0) { + continue; + } + + // Lambda -> p + pi + if (RecoDecay::isMatchedMCGen(mcParticles, daughterCascade, +kLambda0, std::array{+kProton, +kPiMinus}, true)) { + debugGenV0 = 1; // -> Matched Lambda0 + flag = sign * (1 << decayChannel); + } + } // V0 daughter loop + } // cascade daughter loop + } + } // charm daughter loop + + // Check if charm is prompt or non-prompt + if (flag != 0) { + origin = RecoDecay::getCharmHadronOrigin(mcParticles, particle, false, &idxBhadMothers); + } + if (std::abs(yCharmBaryonGen) < kYCutTight) { + // Fill in some QA histogram. Will be implemented later + } + if (std::abs(yCharmBaryonGen) < kYCutLoose) { + // Fill in some QA histograms. Will be implemented later + } + + if (origin == RecoDecay::OriginType::NonPrompt) { + fillGenMcTableByDecayChannel(flag, debugGenCharmBaryon, debugGenCasc, debugGenV0, ptCharmBaryonGen, yCharmBaryonGen, origin, idxBhadMothers[0]); + } else { + fillGenMcTableByDecayChannel(flag, debugGenCharmBaryon, debugGenCasc, debugGenV0, ptCharmBaryonGen, yCharmBaryonGen, origin, -1); + } + + } // particle loop + + } // end of collision loop + + } // template run function + + void processMcEmpty(aod::Collisions const&) + { + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcEmpty, "Empty process function to prevent workflow from getting stuck", true); + + ///////////////////////////////////////////////////// + /// /// + /// Process functions with DCAFitter /// + /// /// + ///////////////////////////////////////////////////// + + //~~~~~~~~~~~~~~~~// + //~~~~To Xi Pi~~~~// + //~~~~~~~~~~~~~~~~// + void processMcXicToXiPiWithDCAFitterNoCent(aod::HfCandToXiPi const& candidates, + aod::TracksWMc const& tracks, + aod::McParticles const& mcParticles, + aod::McCollisions const& mcCollisions, + McCollisionsNoCents const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcXicToXiPiWithDCAFitterNoCent, "Perform MC matching of DCAFitter reconstructed Xic0 to Xi Pi. No cents", false); + + void processMcXicToXiPiWithDCAFitterCentFT0C(aod::HfCandToXiPi const& candidates, + aod::TracksWMc const& tracks, + aod::McParticles const& mcParticles, + aod::McCollisions const& mcCollisions, + McCollisionsFT0Cs const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcXicToXiPiWithDCAFitterCentFT0C, "Perform MC matching of DCAFitter reconstructed Xic0 to Xi Pi. Cents with FT0C", false); + + void processMcXicToXiPiWithDCAFitterCentFT0M(aod::HfCandToXiPi const& candidates, + aod::TracksWMc const& tracks, + aod::McParticles const& mcParticles, + McCollisionsCentFT0Ms const& mcCollisions, + McCollisionsFT0Ms const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcXicToXiPiWithDCAFitterCentFT0M, "Perform MC matching of DCAFitter reconstructed Xic0 to Xi Pi. Cents with FT0M", false); + + void processMcOmegacToXiPiWithDCAFitterNoCent(aod::HfCandToXiPi const& candidates, + aod::TracksWMc const& tracks, + aod::McParticles const& mcParticles, + aod::McCollisions const& mcCollisions, + McCollisionsNoCents const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcOmegacToXiPiWithDCAFitterNoCent, "Perform MC matching of DCAFitter reconstructed Omegac0 to Xi Pi. No cents", false); + + void processMcOmegacToXiPiWithDCAFitterCentFT0C(aod::HfCandToXiPi const& candidates, + aod::TracksWMc const& tracks, + aod::McParticles const& mcParticles, + aod::McCollisions const& mcCollisions, + McCollisionsFT0Cs const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcOmegacToXiPiWithDCAFitterCentFT0C, "Perform MC matching of DCAFitter reconstructed Omeagc0 to Xi Pi. Cents with FT0C", false); + + void processMcOmegacToXiPiWithDCAFitterCentFT0M(aod::HfCandToXiPi const& candidates, + aod::TracksWMc const& tracks, + aod::McParticles const& mcParticles, + McCollisionsCentFT0Ms const& mcCollisions, + McCollisionsFT0Ms const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcOmegacToXiPiWithDCAFitterCentFT0M, "Perform MC matching of DCAFitter reconstructed Omegac0 to Xi Pi. Cents with FT0M", false); + + //~~~~~~~~~~~~~~~~~~~// + //~~~~To Omega Pi~~~~// + //~~~~~~~~~~~~~~~~~~~// + void processMcOmegacToOmegaPiWithDCAFitterNoCent(aod::HfCandToOmegaPi const& candidates, + aod::TracksWMc const& tracks, + aod::McParticles const& mcParticles, + aod::McCollisions const& mcCollisions, + McCollisionsNoCents const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcOmegacToOmegaPiWithDCAFitterNoCent, "Perform MC matching of DCAFitter reconstructed Omegac0 to Omega Pi. No cents", false); + + void processMcOmegacToOmegaPiWithDCAFitterCentFT0C(aod::HfCandToOmegaPi const& candidates, + aod::TracksWMc const& tracks, + aod::McParticles const& mcParticles, + aod::McCollisions const& mcCollisions, + McCollisionsFT0Cs const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcOmegacToOmegaPiWithDCAFitterCentFT0C, "Perform MC matching of DCAFitter reconstructed Omegac0 to Omega Pi. Cents with FT0C", false); + + void processMcOmegacToOmegaPiWithDCAFitterCentFT0M(aod::HfCandToOmegaPi const& candidates, + aod::TracksWMc const& tracks, + aod::McParticles const& mcParticles, + McCollisionsCentFT0Ms const& mcCollisions, + McCollisionsFT0Ms const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcOmegacToOmegaPiWithDCAFitterCentFT0M, "Perform MC matching of DCAFitter reconstructed Omegac0 to Omega Pi. Cents with FT0M", false); + + //~~~~~~~~~~~~~~~~~~// + //~~~~To Omega Ka~~~~// + //~~~~~~~~~~~~~~~~~~// + void processMcOmegacToOmegaKaWithDCAFitterNoCent(aod::HfCandToOmegaK const& candidates, + aod::TracksWMc const& tracks, + aod::McParticles const& mcParticles, + aod::McCollisions const& mcCollisions, + McCollisionsNoCents const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcOmegacToOmegaKaWithDCAFitterNoCent, "Perform MC matching of DCAFitter reconstructed Omegac0 to Omega Ka. No cents", false); + + void processMcOmegacToOmegaKaWithDCAFitterCentFT0C(aod::HfCandToOmegaK const& candidates, + aod::TracksWMc const& tracks, + aod::McParticles const& mcParticles, + aod::McCollisions const& mcCollisions, + McCollisionsFT0Cs const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcOmegacToOmegaKaWithDCAFitterCentFT0C, "Perform MC matching of DCAFitter reconstructed Omegac0 to Omega Ka. Cents with FT0C", false); + + void processMcOmegacToOmegaKaWithDCAFitterCentFT0M(aod::HfCandToOmegaK const& candidates, + aod::TracksWMc const& tracks, + aod::McParticles const& mcParticles, + McCollisionsCentFT0Ms const& mcCollisions, + McCollisionsFT0Ms const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcOmegacToOmegaKaWithDCAFitterCentFT0M, "Perform MC matching of DCAFitter reconstructed Omegac0 to Omega Ka. Cents with FT0M", false); + + ///////////////////////////////////////////////////// + /// /// + /// Process functions with KFParticle /// + /// /// + ///////////////////////////////////////////////////// + void processMcXicToXiPiWithKFParticleNoCent(aod::HfCandToXiPiKf const& candidates, + aod::TracksWMc const& tracks, + aod::McParticles const& mcParticles, + aod::McCollisions const& mcCollisions, + McCollisionsNoCents const& collsWithMcLabels, + BCsInfo const& bcs) + { + runXic0Omegac0Mc(candidates, tracks, mcParticles, collsWithMcLabels, mcCollisions, bcs); + } + PROCESS_SWITCH(HfCandidateCreatorXic0Omegac0QaMc, processMcXicToXiPiWithKFParticleNoCent, "Perform MC matching of DCAFitter reconstructed Xic0 to Xi Pi. No cents", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc)}; +} diff --git a/PWGHF/TableProducer/candidateSelectorToXiPiQa.cxx b/PWGHF/TableProducer/candidateSelectorToXiPiQa.cxx new file mode 100644 index 00000000000..cfd835857e7 --- /dev/null +++ b/PWGHF/TableProducer/candidateSelectorToXiPiQa.cxx @@ -0,0 +1,640 @@ +// 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 candidateSelectorToXiPiQa.cxx +/// \brief Selection of Xic0 and Xicp candidates +/// +/// \author Jinhyun Park , Pusan National University +/// \author Krista Smith , Pusan National University + +#include "PWGHF/Core/HfMlResponseXic0ToXiPiKf.h" +#include "PWGHF/Core/SelectorCuts.h" +#include "PWGHF/DataModel/AliasTables.h" +#include "PWGHF/DataModel/CandidateReconstructionTables.h" +#include "PWGHF/DataModel/CandidateSelectionTables.h" +#include "PWGHF/Utils/utilsAnalysis.h" + +#include "Common/Core/RecoDecay.h" +#include "Common/Core/TrackSelectorPID.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include +#include +#include +#include + +// Related to ML + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::analysis; + +enum PidInfoStored { + PiFromLam = 0, + PrFromLam, + PiFromCasc, + PiFromCharm +}; + +/// Struct for applying Omegac0/Xic0 selection cuts +struct HfCandidateSelectorToXiPiQa { + + // DCAFitter + Produces hfSelToXiPi; + // KFParticle + Produces hfSelToXiPiKf; + Produces hfMlToXiPiKf; // -> Suffix "kf" of table name was deleted + + // cuts from SelectorCuts.h + Configurable> binsPt{"binsPt", std::vector{hf_cuts_to_xi_pi::vecBinsPt}, "pT bin limits"}; + Configurable> cuts{"cuts", {hf_cuts_to_xi_pi::Cuts[0], hf_cuts_to_xi_pi::NBinsPt, hf_cuts_to_xi_pi::NCutVars, hf_cuts_to_xi_pi::labelsPt, hf_cuts_to_xi_pi::labelsCutVar}, "Xic0 candidate selection per Pt Bin"}; + + // ML inference + Configurable applyMl{"applyMl", false, "Flag to apply ML selections"}; + Configurable> binsPtMl{"binsPtMl", std::vector{hf_cuts_ml::vecBinsPt}, "pT bin limits for ML application"}; + Configurable> cutDirMl{"cutDirMl", std::vector{hf_cuts_ml::vecCutDir}, "Whether to reject score values greater or smaller than the threshold"}; + Configurable> cutsMl{"cutsMl", {hf_cuts_ml::Cuts[0], hf_cuts_ml::NBinsPt, hf_cuts_ml::NCutScores, hf_cuts_ml::labelsPt, hf_cuts_ml::labelsCutScore}, "ML selections per pT bin"}; + Configurable nClassesMl{"nClassesMl", static_cast(hf_cuts_ml::NCutScores), "Number of classes in ML model"}; + Configurable> namesInputFeatures{"namesInputFeatures", std::vector{"feature1", "feature2"}, "Names of ML model input features"}; + + // CCDB configuration + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable> modelPathsCCDB{"modelPathsCCDB", std::vector{"EventFiltering/PWGHF/BDTXic0ToXipiKf"}, "Paths of models on CCDB"}; + Configurable> onnxFileNames{"onnxFileNames", std::vector{"ModelHandler_onnx_Xic0ToXipiKf.onnx"}, "ONNX file names for each pT bin (if not from CCDB full path)"}; + Configurable timestampCCDB{"timestampCCDB", -1, "timestamp of the ONNX file for ML model used to query in CCDB"}; + Configurable loadModelsFromCCDB{"loadModelsFromCCDB", false, "Flag to enable or disable the loading of models from CCDB"}; + + // LF analysis selections + Configurable applyTrkSelLf{"applyTrkSelLf", true, "Apply track selection for LF daughters"}; + + // PID options + Configurable usePidTpcOnly{"usePidTpcOnly", false, "Perform PID using only TPC"}; + Configurable usePidTpcTofCombined{"usePidTpcTofCombined", true, "Perform PID using TPC & TOF"}; + + // PID - TPC selections + Configurable ptPiPidTpcMin{"ptPiPidTpcMin", -1, "Lower bound of track pT for TPC PID for pion selection"}; + Configurable ptPiPidTpcMax{"ptPiPidTpcMax", 9999.9, "Upper bound of track pT for TPC PID for pion selection"}; + Configurable nSigmaTpcPiMax{"nSigmaTpcPiMax", 3., "Nsigma cut on TPC only for pion selection"}; + Configurable nSigmaTpcCombinedPiMax{"nSigmaTpcCombinedPiMax", 0., "Nsigma cut on TPC combined with TOF for pion selection"}; + + Configurable ptPrPidTpcMin{"ptPrPidTpcMin", -1, "Lower bound of track pT for TPC PID for proton selection"}; + Configurable ptPrPidTpcMax{"ptPrPidTpcMax", 9999.9, "Upper bound of track pT for TPC PID for proton selection"}; + Configurable nSigmaTpcPrMax{"nSigmaTpcPrMax", 3., "Nsigma cut on TPC only for proton selection"}; + Configurable nSigmaTpcCombinedPrMax{"nSigmaTpcCombinedPrMax", 0., "Nsigma cut on TPC combined with TOF for proton selection"}; + + // PID - TOF selections + Configurable ptPiPidTofMin{"ptPiPidTofMin", -1, "Lower bound of track pT for TOF PID for pion selection"}; + Configurable ptPiPidTofMax{"ptPiPidTofMax", 9999.9, "Upper bound of track pT for TOF PID for pion selection"}; + Configurable nSigmaTofPiMax{"nSigmaTofPiMax", 3., "Nsigma cut on TOF only for pion selection"}; + Configurable nSigmaTofCombinedPiMax{"nSigmaTofCombinedPiMax", 0., "Nsigma cut on TOF combined with TPC for pion selection"}; + + Configurable ptPrPidTofMin{"ptPrPidTofMin", -1, "Lower bound of track pT for TOF PID for proton selection"}; + Configurable ptPrPidTofMax{"ptPrPidTofMax", 9999.9, "Upper bound of track pT for TOF PID for proton selection"}; + Configurable nSigmaTofPrMax{"nSigmaTofPrMax", 3., "Nsigma cut on TOF only for proton selection"}; + Configurable nSigmaTofCombinedPrMax{"nSigmaTofCombinedPrMax", 0., "Nsigma cut on TOF combined with TPC for proton selection"}; + + // detector track quality selections + Configurable nClustersTpcMin{"nClustersTpcMin", 70, "Minimum number of TPC clusters requirement"}; + Configurable nTpcCrossedRowsMin{"nTpcCrossedRowsMin", 70, "Minimum number of TPC crossed rows requirement"}; + Configurable tpcCrossedRowsOverFindableClustersRatioMin{"tpcCrossedRowsOverFindableClustersRatioMin", 0.8, "Minimum ratio TPC crossed rows over findable clusters requirement"}; + Configurable tpcChi2PerClusterMax{"tpcChi2PerClusterMax", 4, "Maximum value of chi2 fit over TPC clusters"}; + Configurable nClustersItsMin{"nClustersItsMin", 3, "Minimum number of ITS clusters requirement for pi <- charm baryon"}; + Configurable nClustersItsInnBarrMin{"nClustersItsInnBarrMin", 1, "Minimum number of ITS clusters in inner barrel requirement for pi <- charm baryon"}; + Configurable itsChi2PerClusterMax{"itsChi2PerClusterMax", 36, "Maximum value of chi2 fit over ITS clusters for pi <- charm baryon"}; + + o2::analysis::HfMlResponseXic0ToXiPiKf hfMlResponse; + std::vector outputMlXic0ToXiPiKf = {}; + o2::ccdb::CcdbApi ccdbApi; + + TrackSelectorPi selectorPion; + TrackSelectorPr selectorProton; + + using TracksSel = soa::Join; + + HistogramRegistry registry{"registry"}; // for QA of selections + + void init(InitContext const&) + { + selectorPion.setRangePtTpc(ptPiPidTpcMin, ptPiPidTpcMax); + selectorPion.setRangeNSigmaTpc(-nSigmaTpcPiMax, nSigmaTpcPiMax); + selectorPion.setRangeNSigmaTpcCondTof(-nSigmaTpcCombinedPiMax, nSigmaTpcCombinedPiMax); + selectorPion.setRangePtTof(ptPiPidTofMin, ptPiPidTofMax); + selectorPion.setRangeNSigmaTof(-nSigmaTofPiMax, nSigmaTofPiMax); + selectorPion.setRangeNSigmaTofCondTpc(-nSigmaTofCombinedPiMax, nSigmaTofCombinedPiMax); + + selectorProton.setRangePtTpc(ptPrPidTpcMin, ptPrPidTpcMax); + selectorProton.setRangeNSigmaTpc(-nSigmaTpcPrMax, nSigmaTpcPrMax); + selectorProton.setRangeNSigmaTpcCondTof(-nSigmaTpcCombinedPrMax, nSigmaTpcCombinedPrMax); + selectorProton.setRangePtTof(ptPrPidTofMin, ptPrPidTofMax); + selectorProton.setRangeNSigmaTof(-nSigmaTofPrMax, nSigmaTofPrMax); + selectorProton.setRangeNSigmaTofCondTpc(-nSigmaTofCombinedPrMax, nSigmaTofCombinedPrMax); + + const AxisSpec axisSel{2, -0.5, 1.5, "status"}; + + registry.add("hSelPID", "hSelPID;status;entries", {HistType::kTH1F, {{12, 0., 12.}}}); + registry.add("hStatusCheck", "Check consecutive selections status;status;entries", {HistType::kTH1F, {{12, 0., 12.}}}); + + // for QA of the selections (bin 0 -> candidates that did not pass the selection, bin 1 -> candidates that passed the selection) + registry.add("hSelSignDec", "hSelSignDec;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelEtaPosV0Dau", "hSelEtaPosV0Dau;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelEtaNegV0Dau", "hSelEtaNegV0Dau;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelEtaPiFromCasc", "hSelEtaPiFromCasc;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelEtaPiFromCharm", "hSelEtaPiFromCharm;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelRadCasc", "hSelRadCasc;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelRadV0", "hSelRadV0;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelCosPACasc", "hSelCosPACasc;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelCosPAV0", "hSelCosPAV0;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelDCACascDau", "hSelDCACascDau;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelDCAV0Dau", "hSelDCAV0Dau;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelDCACharmDau", "hSelDCACharmDau;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelDCAXYPrimPi", "hSelDCAXYPrimPi;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelDCAZPrimPi", "hSelDCAZPrimPi;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelDCAXYCasc", "hSelDCAXYCasc;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelDCAZCasc", "hSelDCAZCasc;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelPtPiFromCasc", "hSelPtPiFromCasc;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelPtPiFromCharm", "hSelPtPiFromCharm;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelTPCQualityPiFromCharm", "hSelTPCQualityPiFromCharm;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelTPCQualityPiFromLam", "hSelTPCQualityPiFromLam;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelTPCQualityPrFromLam", "hSelTPCQualityPrFromLam;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelTPCQualityPiFromCasc", "hSelTPCQualityPiFromCasc;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelITSQualityPiFromCharm", "hSelITSQualityPiFromCharm;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelMassLam", "hSelMassLam;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelMassCasc", "hSelMassCasc;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelMassCharmBaryon", "hSelMassCharmBaryon;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelDcaXYToPvV0Daughters", "hSelDcaXYToPvV0Daughters;status;entries", {HistType::kTH1F, {axisSel}}); + registry.add("hSelDcaXYToPvPiFromCasc", "hSelDcaXYToPvPiFromCasc;status;entries", {HistType::kTH1F, {axisSel}}); + + // invarinat mass histograms + registry.add("hInvMassCharmBaryon", "Charm baryon invariant mass; int mass; entries", {HistType::kTH1F, {{1500, 1.5, 4.5}}}); + registry.add("hInvMassCharmBaryonBkg", "Charm baryon invariant mass, rejected; int mass; entries", {HistType::kTH1F, {{1500, 1.5, 4.5}}}); + + if (applyMl) { + hfMlResponse.configure(binsPtMl, cutsMl, cutDirMl, nClassesMl); + if (loadModelsFromCCDB) { + ccdbApi.init(ccdbUrl); + hfMlResponse.setModelPathsCCDB(onnxFileNames, ccdbApi, modelPathsCCDB, timestampCCDB); + } else { + hfMlResponse.setModelPathsLocal(onnxFileNames); + } + hfMlResponse.cacheInputFeaturesIndices(namesInputFeatures); + hfMlResponse.init(); + } + } + + template + void runSelection(TCandTable const& candidates, TracksSel const&) + { + double massLambdaFromPDG = o2::constants::physics::MassLambda0; + double massXiFromPDG = o2::constants::physics::MassXiMinus; + + // looping over charm baryon candidates + for (const auto& candidate : candidates) { + + bool resultSelections = true; // True if the candidate passes all the selections, False otherwise + + auto trackV0PosDau = candidate.template posTrack_as(); + auto trackV0NegDau = candidate.template negTrack_as(); + auto trackPiFromCasc = candidate.template bachelor_as(); + auto trackPiFromCharm = candidate.template bachelorFromCharmBaryon_as(); + + auto trackPiFromLam = trackV0NegDau; + auto trackPrFromLam = trackV0PosDau; + + int8_t signDecay = candidate.signDecay(); // sign of pi <- cascade + + if (signDecay > 0) { + trackPiFromLam = trackV0PosDau; + trackPrFromLam = trackV0NegDau; + registry.fill(HIST("hSelSignDec"), 1); // anti-particle decay + } else if (signDecay < 0) { + registry.fill(HIST("hSelSignDec"), 0); // particle decay + } + + // pT selection + auto ptCandXic0 = RecoDecay::pt(candidate.pxCharmBaryon(), candidate.pyCharmBaryon()); + int pTBin = findBin(binsPt, ptCandXic0); + if (pTBin == -1) { + resultSelections = false; + continue; + } + + // eta selection + double etaV0PosDau = candidate.etaV0PosDau(); + double etaV0NegDau = candidate.etaV0NegDau(); + double etaPiFromCasc = candidate.etaBachFromCasc(); + double etaPiFromCharmBaryon = candidate.etaBachFromCharmBaryon(); + if (std::abs(etaV0PosDau) > cuts->get(pTBin, "etaTrackLFDauMax")) { + resultSelections = false; + registry.fill(HIST("hSelEtaPosV0Dau"), 0); + } else { + registry.fill(HIST("hSelEtaPosV0Dau"), 1); + } + if (std::abs(etaV0NegDau) > cuts->get(pTBin, "etaTrackLFDauMax")) { + resultSelections = false; + registry.fill(HIST("hSelEtaNegV0Dau"), 0); + } else { + registry.fill(HIST("hSelEtaNegV0Dau"), 1); + } + if (std::abs(etaPiFromCasc) > cuts->get(pTBin, "etaTrackLFDauMax")) { + resultSelections = false; + registry.fill(HIST("hSelEtaPiFromCasc"), 0); + } else { + registry.fill(HIST("hSelEtaPiFromCasc"), 1); + } + if (std::abs(etaPiFromCharmBaryon) > cuts->get(pTBin, "etaTrackCharmBachMax")) { + resultSelections = false; + registry.fill(HIST("hSelEtaPiFromCharm"), 0); + } else { + registry.fill(HIST("hSelEtaPiFromCharm"), 1); + } + + // minimum radius cut (LFcut) + if (RecoDecay::sqrtSumOfSquares(candidate.xDecayVtxCascade(), candidate.yDecayVtxCascade()) < cuts->get(pTBin, "radiusCascMin")) { + resultSelections = false; + registry.fill(HIST("hSelRadCasc"), 0); + } else { + registry.fill(HIST("hSelRadCasc"), 1); + } + if (RecoDecay::sqrtSumOfSquares(candidate.xDecayVtxV0(), candidate.yDecayVtxV0()) < cuts->get(pTBin, "radiusV0Min")) { + resultSelections = false; + registry.fill(HIST("hSelRadV0"), 0); + } else { + registry.fill(HIST("hSelRadV0"), 1); + } + + // cosPA (LFcut) + if (candidate.cosPACasc() < cuts->get(pTBin, "cosPaCascMin")) { + resultSelections = false; + registry.fill(HIST("hSelCosPACasc"), 0); + } else { + registry.fill(HIST("hSelCosPACasc"), 1); + } + if (candidate.cosPAV0() < cuts->get(pTBin, "cosPaV0Min")) { + resultSelections = false; + registry.fill(HIST("hSelCosPAV0"), 0); + } else { + registry.fill(HIST("hSelCosPAV0"), 1); + } + + // cascade and v0 daughters dca cut (LF cut) + if (candidate.dcaCascDau() > cuts->get(pTBin, "dcaCascDauMax")) { + resultSelections = false; + registry.fill(HIST("hSelDCACascDau"), 0); + } else { + registry.fill(HIST("hSelDCACascDau"), 1); + } + + if (candidate.dcaV0Dau() > cuts->get(pTBin, "dcaV0DauMax")) { + resultSelections = false; + registry.fill(HIST("hSelDCAV0Dau"), 0); + } else { + registry.fill(HIST("hSelDCAV0Dau"), 1); + } + + // dca charm baryon daughters cut + if (candidate.dcaCharmBaryonDau() > cuts->get(pTBin, "dcaCharmBaryonDauMax")) { + resultSelections = false; + registry.fill(HIST("hSelDCACharmDau"), 0); + } else { + registry.fill(HIST("hSelDCACharmDau"), 1); + } + + // dcaXY v0 daughters to PV cut + if (std::abs(candidate.dcaXYToPvV0Dau0()) < cuts->get(pTBin, "dcaxyV0PosDauToPvMin") || std::abs(candidate.dcaXYToPvV0Dau1()) < cuts->get(pTBin, "dcaxyV0NegDauToPvMin")) { + resultSelections = false; + registry.fill(HIST("hSelDcaXYToPvV0Daughters"), 0); + } else { + registry.fill(HIST("hSelDcaXYToPvV0Daughters"), 1); + } + + // dcaXY pi <-- cascade to PV cut + if (std::abs(candidate.dcaXYToPvCascDau()) < cuts->get(pTBin, "dcaxyBachToPvMin")) { + resultSelections = false; + registry.fill(HIST("hSelDcaXYToPvPiFromCasc"), 0); + } else { + registry.fill(HIST("hSelDcaXYToPvPiFromCasc"), 1); + } + + // // cut on charm bachelor pion dcaXY and dcaZ + // if ((std::abs(candidate.impactParBachFromCharmBaryonXY()) < cuts->get(pTBin, "impactParXYCharmBachelorMin")) || (std::abs(candidate.impactParBachFromCharmBaryonXY()) > cuts->get(pTBin, "impactParXYCharmBachelorMax"))) { + // resultSelections = false; + // registry.fill(HIST("hSelDCAXYPrimPi"), 0); + // } else { + // registry.fill(HIST("hSelDCAXYPrimPi"), 1); + // } + // if ((std::abs(candidate.impactParBachFromCharmBaryonZ()) < cuts->get(pTBin, "impactParZCharmBachelorMin")) || (std::abs(candidate.impactParBachFromCharmBaryonZ()) > cuts->get(pTBin, "impactParZCharmBachelorMax"))) { + // resultSelections = false; + // registry.fill(HIST("hSelDCAZPrimPi"), 0); + // } else { + // registry.fill(HIST("hSelDCAZPrimPi"), 1); + // } + + // // cut on cascade dcaXY and dcaZ + // if ((std::abs(candidate.impactParCascXY()) < cuts->get(pTBin, "impactParXYCascMin")) || (std::abs(candidate.impactParCascXY()) > cuts->get(pTBin, "impactParXYCascMax"))) { + // resultSelections = false; + // registry.fill(HIST("hSelDCAXYCasc"), 0); + // } else { + // registry.fill(HIST("hSelDCAXYCasc"), 1); + // } + // if ((std::abs(candidate.impactParCascZ()) < cuts->get(pTBin, "impactParZCascMin")) || (std::abs(candidate.impactParCascZ()) > cuts->get(pTBin, "impactParZCascMax"))) { + // resultSelections = false; + // registry.fill(HIST("hSelDCAZCasc"), 0); + // } else { + // registry.fill(HIST("hSelDCAZCasc"), 1); + // } + + // pT selections + double ptPiFromCasc = RecoDecay::sqrtSumOfSquares(candidate.pxBachFromCasc(), candidate.pyBachFromCasc()); + double ptPiFromCharmBaryon = RecoDecay::sqrtSumOfSquares(candidate.pxBachFromCharmBaryon(), candidate.pyBachFromCharmBaryon()); + if (std::abs(ptPiFromCasc) < cuts->get(pTBin, "ptBachelorMin")) { + resultSelections = false; + registry.fill(HIST("hSelPtPiFromCasc"), 0); + } else { + registry.fill(HIST("hSelPtPiFromCasc"), 1); + } + if (std::abs(ptPiFromCharmBaryon) < cuts->get(pTBin, "ptCharmBachelorMin")) { + resultSelections = false; + registry.fill(HIST("hSelPtPiFromCharm"), 0); + } else { + registry.fill(HIST("hSelPtPiFromCharm"), 1); + } + + // TPC clusters selections + if (applyTrkSelLf) { + if (!isSelectedTrackTpcQuality(trackPiFromLam, nClustersTpcMin, nTpcCrossedRowsMin, tpcCrossedRowsOverFindableClustersRatioMin, tpcChi2PerClusterMax)) { + resultSelections = false; + registry.fill(HIST("hSelTPCQualityPiFromLam"), 0); + } else { + registry.fill(HIST("hSelTPCQualityPiFromLam"), 1); + } + if (!isSelectedTrackTpcQuality(trackPrFromLam, nClustersTpcMin, nTpcCrossedRowsMin, tpcCrossedRowsOverFindableClustersRatioMin, tpcChi2PerClusterMax)) { + resultSelections = false; + registry.fill(HIST("hSelTPCQualityPrFromLam"), 0); + } else { + registry.fill(HIST("hSelTPCQualityPrFromLam"), 1); + } + if (!isSelectedTrackTpcQuality(trackPiFromCasc, nClustersTpcMin, nTpcCrossedRowsMin, tpcCrossedRowsOverFindableClustersRatioMin, tpcChi2PerClusterMax)) { + resultSelections = false; + registry.fill(HIST("hSelTPCQualityPiFromCasc"), 0); + } else { + registry.fill(HIST("hSelTPCQualityPiFromCasc"), 1); + } + } + if (!isSelectedTrackTpcQuality(trackPiFromCharm, nClustersTpcMin, nTpcCrossedRowsMin, tpcCrossedRowsOverFindableClustersRatioMin, tpcChi2PerClusterMax)) { + resultSelections = false; + registry.fill(HIST("hSelTPCQualityPiFromCharm"), 0); + } else { + registry.fill(HIST("hSelTPCQualityPiFromCharm"), 1); + } + + // ITS clusters selection + if (!isSelectedTrackItsQuality(trackPiFromCharm, nClustersItsMin, itsChi2PerClusterMax) || trackPiFromCharm.itsNClsInnerBarrel() < nClustersItsInnBarrMin) { + resultSelections = false; + registry.fill(HIST("hSelITSQualityPiFromCharm"), 0); + } else { + registry.fill(HIST("hSelITSQualityPiFromCharm"), 1); + } + + // track-level PID selection + + // for TrackSelectorPID + int statusPidPrFromLam = -999; + int statusPidPiFromLam = -999; + int statusPidPiFromCasc = -999; + int statusPidPiFromCharmBaryon = -999; + + bool statusPidLambda = false; + bool statusPidCascade = false; + bool statusPidCharmBaryon = false; + + int infoTpcStored = 0; + int infoTofStored = 0; + + if (usePidTpcOnly == usePidTpcTofCombined) { + LOGF(fatal, "Check the PID configurables, usePidTpcOnly and usePidTpcTofCombined can't have the same value"); + } + + if (trackPiFromLam.hasTPC()) { + SETBIT(infoTpcStored, PiFromLam); + } + if (trackPrFromLam.hasTPC()) { + SETBIT(infoTpcStored, PrFromLam); + } + if (trackPiFromCasc.hasTPC()) { + SETBIT(infoTpcStored, PiFromCasc); + } + if (trackPiFromCharm.hasTPC()) { + SETBIT(infoTpcStored, PiFromCharm); + } + if (trackPiFromLam.hasTOF()) { + SETBIT(infoTofStored, PiFromLam); + } + if (trackPrFromLam.hasTOF()) { + SETBIT(infoTofStored, PrFromLam); + } + if (trackPiFromCasc.hasTOF()) { + SETBIT(infoTofStored, PiFromCasc); + } + if (trackPiFromCharm.hasTOF()) { + SETBIT(infoTofStored, PiFromCharm); + } + + if (usePidTpcOnly) { + statusPidPrFromLam = selectorProton.statusTpc(trackPrFromLam); + statusPidPiFromLam = selectorPion.statusTpc(trackPiFromLam); + statusPidPiFromCasc = selectorPion.statusTpc(trackPiFromCasc); + statusPidPiFromCharmBaryon = selectorPion.statusTpc(trackPiFromCharm); + } else if (usePidTpcTofCombined) { + statusPidPrFromLam = selectorProton.statusTpcOrTof(trackPrFromLam); + statusPidPiFromLam = selectorPion.statusTpcOrTof(trackPiFromLam); + statusPidPiFromCasc = selectorPion.statusTpcOrTof(trackPiFromCasc); + statusPidPiFromCharmBaryon = selectorPion.statusTpcOrTof(trackPiFromCharm); + } + + if (statusPidPrFromLam == TrackSelectorPID::Accepted && statusPidPiFromLam == TrackSelectorPID::Accepted) { + statusPidLambda = true; + if (resultSelections) { + registry.fill(HIST("hStatusCheck"), 0.5); + } + } + + if (statusPidPrFromLam == TrackSelectorPID::Accepted && statusPidPiFromLam == TrackSelectorPID::Accepted && statusPidPiFromCasc == TrackSelectorPID::Accepted) { + statusPidCascade = true; + if (resultSelections) { + registry.fill(HIST("hStatusCheck"), 1.5); + } + } + + if (statusPidPrFromLam == TrackSelectorPID::Accepted && statusPidPiFromLam == TrackSelectorPID::Accepted && statusPidPiFromCasc == TrackSelectorPID::Accepted && statusPidPiFromCharmBaryon == TrackSelectorPID::Accepted) { + statusPidCharmBaryon = true; + if (resultSelections) { + registry.fill(HIST("hStatusCheck"), 2.5); + } + } + + // invariant mass cuts + bool statusInvMassLambda = false; + bool statusInvMassCascade = false; + bool statusInvMassCharmBaryon = false; + + double invMassLambda = candidate.invMassLambda(); + double invMassCascade = candidate.invMassCascade(); + double invMassCharmBaryon = candidate.invMassCharmBaryon(); + + if (std::abs(invMassLambda - massLambdaFromPDG) < cuts->get(pTBin, "massWindowV0")) { + statusInvMassLambda = true; + registry.fill(HIST("hSelMassLam"), 1); + if (statusPidLambda && statusPidCascade && statusPidCharmBaryon && resultSelections) { + registry.fill(HIST("hStatusCheck"), 3.5); + } + } else { + registry.fill(HIST("hSelMassLam"), 0); + } + + if (std::abs(invMassCascade - massXiFromPDG) < cuts->get(pTBin, "massWindowCascade")) { + statusInvMassCascade = true; + registry.fill(HIST("hSelMassCasc"), 1); + if (statusPidLambda && statusPidCascade && statusPidCharmBaryon && statusInvMassLambda && resultSelections) { + registry.fill(HIST("hStatusCheck"), 4.5); + } + } else { + registry.fill(HIST("hSelMassCasc"), 0); + } + + if ((invMassCharmBaryon >= cuts->get(pTBin, "mCharmBaryonMin")) && (invMassCharmBaryon <= cuts->get(pTBin, "mCharmBaryonMax"))) { + statusInvMassCharmBaryon = true; + registry.fill(HIST("hSelMassCharmBaryon"), 1); + if (statusPidLambda && statusPidCascade && statusPidCharmBaryon && statusInvMassLambda && statusInvMassCascade && resultSelections) { + registry.fill(HIST("hStatusCheck"), 5.5); + } + } else { + registry.fill(HIST("hSelMassCharmBaryon"), 0); + } + + // Fill in selection result + if constexpr (dokf == false) { + hfSelToXiPi(statusPidLambda, statusPidCascade, statusPidCharmBaryon, statusInvMassLambda, statusInvMassCascade, statusInvMassCharmBaryon, resultSelections, infoTpcStored, infoTofStored, + trackPiFromCharm.tpcNSigmaPi(), trackPiFromCasc.tpcNSigmaPi(), trackPiFromLam.tpcNSigmaPi(), trackPrFromLam.tpcNSigmaPr(), + trackPiFromCharm.tofNSigmaPi(), trackPiFromCasc.tofNSigmaPi(), trackPiFromLam.tofNSigmaPi(), trackPrFromLam.tofNSigmaPr()); + } else { + // ML selections -> Currently, no Ml for DCAFitter result is implemented + if (applyMl) { + bool isSelectedMlXic0 = false; + std::vector inputFeaturesXic0 = hfMlResponse.getInputFeatures(candidate, trackPiFromLam, trackPiFromCasc, trackPiFromCharm); + isSelectedMlXic0 = hfMlResponse.isSelectedMl(inputFeaturesXic0, ptCandXic0, outputMlXic0ToXiPiKf); + if (!isSelectedMlXic0) { + continue; + } + hfMlToXiPiKf(outputMlXic0ToXiPiKf); + } + + // Set result of selection to false if one of the statusPID or statusInvMass is false + // This is required because selection table for KF does not store any information about statusPID or statusInvMass while DCAFitter does. + if (!(statusPidLambda && statusPidCascade && statusPidCharmBaryon && statusInvMassCharmBaryon && statusInvMassCascade && statusInvMassLambda)) { + resultSelections = false; + } + + hfSelToXiPiKf(resultSelections, + // statusPidCharmBaryon, statusPidCascade, statusPidLambda, statusInvMassCharmBaryon, statusInvMassCascade, statusInvMassLambda, infoTpcStored, infoTofStored, + trackPiFromCharm.tpcNSigmaPi(), trackPiFromCasc.tpcNSigmaPi(), trackPiFromLam.tpcNSigmaPi(), trackPrFromLam.tpcNSigmaPr(), + trackPiFromCharm.tofNSigmaPi(), trackPiFromCasc.tofNSigmaPi(), trackPiFromLam.tofNSigmaPi(), trackPrFromLam.tofNSigmaPr()); + } + + if (resultSelections) { + if (!statusPidLambda) { + registry.fill(HIST("hSelPID"), 0.5); + } + if (statusPidLambda) { + registry.fill(HIST("hSelPID"), 1.5); + } + if (!statusPidCascade) { + registry.fill(HIST("hSelPID"), 2.5); + } + if (statusPidCascade) { + registry.fill(HIST("hSelPID"), 3.5); + } + if (!statusPidCharmBaryon) { + registry.fill(HIST("hSelPID"), 4.5); + } + if (statusPidCharmBaryon) { + registry.fill(HIST("hSelPID"), 5.5); + } + if (!statusInvMassLambda) { + registry.fill(HIST("hSelPID"), 6.5); + } + if (statusInvMassLambda) { + registry.fill(HIST("hSelPID"), 7.5); + } + if (!statusInvMassCascade) { + registry.fill(HIST("hSelPID"), 8.5); + } + if (statusInvMassCascade) { + registry.fill(HIST("hSelPID"), 9.5); + } + if (!statusInvMassCharmBaryon) { + registry.fill(HIST("hSelPID"), 10.5); + } + if (statusInvMassCharmBaryon) { + registry.fill(HIST("hSelPID"), 11.5); + } + } + + // Fill in invariant mass histogram + if (statusPidLambda && statusPidCascade && statusPidCharmBaryon && statusInvMassLambda && statusInvMassCascade && statusInvMassCharmBaryon && resultSelections) { + registry.fill(HIST("hInvMassCharmBaryon"), invMassCharmBaryon); + } else { + registry.fill(HIST("hInvMassCharmBaryonBkg"), invMassCharmBaryon); + } + + } // end of candidate loop + } // end run fuction + + /////////////////////////////////// + /// Process with DCAFitter // + /////////////////////////////////// + void processSelectionDCAFitter(aod::HfCandToXiPi const& candidates, TracksSel const& tracks) + { + runSelection(candidates, tracks); + } + PROCESS_SWITCH(HfCandidateSelectorToXiPiQa, processSelectionDCAFitter, "Xic0 candidate selection with DCAFitter output", true); + + //////////////////////////////////// + /// Process with KFParticle // + //////////////////////////////////// + void processSelectionKFParticle(aod::HfCandToXiPiKf const& candidates, TracksSel const& tracks) + { + runSelection(candidates, tracks); + } + PROCESS_SWITCH(HfCandidateSelectorToXiPiQa, processSelectionKFParticle, "Xic0 candidate selection with KFParticle output", false); + +}; // end struct + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc)}; +} diff --git a/PWGHF/TableProducer/treeCreatorToXiPiQa.cxx b/PWGHF/TableProducer/treeCreatorToXiPiQa.cxx new file mode 100644 index 00000000000..6df52b8d01e --- /dev/null +++ b/PWGHF/TableProducer/treeCreatorToXiPiQa.cxx @@ -0,0 +1,990 @@ +// 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 treeCreatorToXiPiQa.cxx +/// \brief Writer of the omegac0 or xic0 to Xi Pi candidates in the form of flat tables to be stored in TTrees. +/// In this file are defined and filled the output tables +/// +/// \author Jinhyun Park , Pusan National University +/// \author Krista Smith , Pusan National University + +#include "PWGHF/Core/CentralityEstimation.h" +#include "PWGHF/DataModel/CandidateReconstructionTables.h" +#include "PWGHF/DataModel/CandidateSelectionTables.h" + +#include "Common/Core/RecoDecay.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +using namespace o2; +using namespace o2::framework; + +namespace o2::aod +{ +namespace full +{ +// collision info +DECLARE_SOA_COLUMN(IsEventSel8, isEventSel8, bool); +DECLARE_SOA_COLUMN(IsEventSelZ, isEventSelZ, bool); +DECLARE_SOA_COLUMN(Centrality, centrality, float); +// from creator +DECLARE_SOA_COLUMN(XPv, xPv, float); +DECLARE_SOA_COLUMN(YPv, yPv, float); +DECLARE_SOA_COLUMN(ZPv, zPv, float); +DECLARE_SOA_COLUMN(XDecayVtxCharmBaryon, xDecayVtxCharmBaryon, float); +DECLARE_SOA_COLUMN(YDecayVtxCharmBaryon, yDecayVtxCharmBaryon, float); +DECLARE_SOA_COLUMN(ZDecayVtxCharmBaryon, zDecayVtxCharmBaryon, float); +DECLARE_SOA_COLUMN(XDecayVtxCascade, xDecayVtxCascade, float); +DECLARE_SOA_COLUMN(YDecayVtxCascade, yDecayVtxCascade, float); +DECLARE_SOA_COLUMN(ZDecayVtxCascade, zDecayVtxCascade, float); +DECLARE_SOA_COLUMN(XDecayVtxV0, xDecayVtxV0, float); +DECLARE_SOA_COLUMN(YDecayVtxV0, yDecayVtxV0, float); +DECLARE_SOA_COLUMN(ZDecayVtxV0, zDecayVtxV0, float); +DECLARE_SOA_COLUMN(SignDecay, signDecay, int8_t); // sign of pi <- xi +DECLARE_SOA_COLUMN(CovVtxCharmBaryonXX, covVtxCharmBaryonXX, float); +DECLARE_SOA_COLUMN(CovVtxCharmBaryonYY, covVtxCharmBaryonYY, float); +DECLARE_SOA_COLUMN(CovVtxCharmBaryonZZ, covVtxCharmBaryonZZ, float); +DECLARE_SOA_COLUMN(PxCharmBaryon, pxCharmBaryon, float); +DECLARE_SOA_COLUMN(PyCharmBaryon, pyCharmBaryon, float); +DECLARE_SOA_COLUMN(PzCharmBaryon, pzCharmBaryon, float); +DECLARE_SOA_COLUMN(PxCasc, pxCasc, float); +DECLARE_SOA_COLUMN(PyCasc, pyCasc, float); +DECLARE_SOA_COLUMN(PzCasc, pzCasc, float); +DECLARE_SOA_COLUMN(PxPiFromCharmBaryon, pxPiFromCharmBaryon, float); +DECLARE_SOA_COLUMN(PyPiFromCharmBaryon, pyPiFromCharmBaryon, float); +DECLARE_SOA_COLUMN(PzPiFromCharmBaryon, pzPiFromCharmBaryon, float); +DECLARE_SOA_COLUMN(PxLambda, pxLambda, float); +DECLARE_SOA_COLUMN(PyLambda, pyLambda, float); +DECLARE_SOA_COLUMN(PzLambda, pzLambda, float); +DECLARE_SOA_COLUMN(PxPiFromCasc, pxPiFromCasc, float); +DECLARE_SOA_COLUMN(PyPiFromCasc, pyPiFromCasc, float); +DECLARE_SOA_COLUMN(PzPiFromCasc, pzPiFromCasc, float); +DECLARE_SOA_COLUMN(PxPosV0Dau, pxPosV0Dau, float); +DECLARE_SOA_COLUMN(PyPosV0Dau, pyPosV0Dau, float); +DECLARE_SOA_COLUMN(PzPosV0Dau, pzPosV0Dau, float); +DECLARE_SOA_COLUMN(PxNegV0Dau, pxNegV0Dau, float); +DECLARE_SOA_COLUMN(PyNegV0Dau, pyNegV0Dau, float); +DECLARE_SOA_COLUMN(PzNegV0Dau, pzNegV0Dau, float); +DECLARE_SOA_COLUMN(ImpactParCascXY, impactParCascXY, float); +DECLARE_SOA_COLUMN(ImpactParPiFromCharmBaryonXY, impactParPiFromCharmBaryonXY, float); +DECLARE_SOA_COLUMN(ImpactParCascZ, impactParCascZ, float); +DECLARE_SOA_COLUMN(ImpactParPiFromCharmBaryonZ, impactParPiFromCharmBaryonZ, float); +DECLARE_SOA_COLUMN(ErrImpactParCascXY, errImpactParCascXY, float); +DECLARE_SOA_COLUMN(ErrImpactParPiFromCharmBaryonXY, errImpactParPiFromCharmBaryonXY, float); +DECLARE_SOA_COLUMN(InvMassLambda, invMassLambda, float); +DECLARE_SOA_COLUMN(InvMassCascade, invMassCascade, float); +DECLARE_SOA_COLUMN(InvMassCharmBaryon, invMassCharmBaryon, float); +DECLARE_SOA_COLUMN(CosPAV0, cosPAV0, float); +DECLARE_SOA_COLUMN(CosPACharmBaryon, cosPACharmBaryon, float); +DECLARE_SOA_COLUMN(CosPACasc, cosPACasc, float); +DECLARE_SOA_COLUMN(CosPAXYV0, cosPAXYV0, float); +DECLARE_SOA_COLUMN(CosPAXYCharmBaryon, cosPAXYCharmBaryon, float); +DECLARE_SOA_COLUMN(CosPAXYCasc, cosPAXYCasc, float); +DECLARE_SOA_COLUMN(CTauOmegac, cTauOmegac, float); +DECLARE_SOA_COLUMN(CTauCascade, cTauCascade, float); +DECLARE_SOA_COLUMN(CTauV0, cTauV0, float); +DECLARE_SOA_COLUMN(CTauXic, cTauXic, float); +DECLARE_SOA_COLUMN(EtaV0PosDau, etaV0PosDau, float); +DECLARE_SOA_COLUMN(EtaV0NegDau, etaV0NegDau, float); +DECLARE_SOA_COLUMN(EtaPiFromCasc, etaPiFromCasc, float); +DECLARE_SOA_COLUMN(EtaPiFromCharmBaryon, etaPiFromCharmBaryon, float); +DECLARE_SOA_COLUMN(EtaCharmBaryon, etaCharmBaryon, float); +DECLARE_SOA_COLUMN(EtaCascade, etaCascade, float); +DECLARE_SOA_COLUMN(EtaV0, etaV0, float); +DECLARE_SOA_COLUMN(DcaXYToPvV0Dau0, dcaXYToPvV0Dau0, float); +DECLARE_SOA_COLUMN(DcaXYToPvV0Dau1, dcaXYToPvV0Dau1, float); +DECLARE_SOA_COLUMN(DcaXYToPvCascDau, dcaXYToPvCascDau, float); +DECLARE_SOA_COLUMN(DcaZToPvV0Dau0, dcaZToPvV0Dau0, float); +DECLARE_SOA_COLUMN(DcaZToPvV0Dau1, dcaZToPvV0Dau1, float); +DECLARE_SOA_COLUMN(DcaZToPvCascDau, dcaZToPvCascDau, float); +DECLARE_SOA_COLUMN(DcaCascDau, dcaCascDau, float); +DECLARE_SOA_COLUMN(DcaV0Dau, dcaV0Dau, float); +DECLARE_SOA_COLUMN(DcaCharmBaryonDau, dcaCharmBaryonDau, float); +DECLARE_SOA_COLUMN(DecLenCharmBaryon, decLenCharmBaryon, float); +DECLARE_SOA_COLUMN(DecLenCascade, decLenCascade, float); +DECLARE_SOA_COLUMN(DecLenV0, decLenV0, float); +DECLARE_SOA_COLUMN(ErrorDecayLengthCharmBaryon, errorDecayLengthCharmBaryon, float); +DECLARE_SOA_COLUMN(ErrorDecayLengthXYCharmBaryon, errorDecayLengthXYCharmBaryon, float); +DECLARE_SOA_COLUMN(NormImpParCascade, normImpParCascade, double); +DECLARE_SOA_COLUMN(NormImpParPiFromCharmBar, normImpParPiFromCharmBar, double); +DECLARE_SOA_COLUMN(NormDecayLenCharmBar, normDecayLenCharmBar, double); +DECLARE_SOA_COLUMN(IsPionGlbTrkWoDca, isPionGlbTrkWoDca, bool); +DECLARE_SOA_COLUMN(PionItsNCls, pionItsNCls, uint8_t); +DECLARE_SOA_COLUMN(NTpcRowsPion, nTpcRowsPion, int16_t); +DECLARE_SOA_COLUMN(NTpcRowsPiFromCasc, nTpcRowsPiFromCasc, int16_t); +DECLARE_SOA_COLUMN(NTpcRowsPosV0Dau, nTpcRowsPosV0Dau, int16_t); +DECLARE_SOA_COLUMN(NTpcRowsNegV0Dau, nTpcRowsNegV0Dau, int16_t); +// from creator KF +DECLARE_SOA_COLUMN(KfDcaXYPiFromXic, kfDcaXYPiFromXic, float); +DECLARE_SOA_COLUMN(KfDcaXYCascToPv, kfDcaXYCascToPv, float); +DECLARE_SOA_COLUMN(Chi2GeoV0, chi2GeoV0, float); +DECLARE_SOA_COLUMN(Chi2GeoCasc, chi2GeoCasc, float); +DECLARE_SOA_COLUMN(Chi2GeoXic, chi2GeoXic, float); +DECLARE_SOA_COLUMN(Chi2MassV0, chi2MassV0, float); +DECLARE_SOA_COLUMN(Chi2MassCasc, chi2MassCasc, float); +DECLARE_SOA_COLUMN(V0ldl, v0ldl, float); +DECLARE_SOA_COLUMN(Cascldl, cascldl, float); +DECLARE_SOA_COLUMN(Xicldl, xicldl, float); +DECLARE_SOA_COLUMN(Chi2TopoV0ToPv, chi2TopoV0ToPv, float); +DECLARE_SOA_COLUMN(Chi2TopoCascToPv, chi2TopoCascToPv, float); +DECLARE_SOA_COLUMN(Chi2TopoPiFromXicToPv, chi2TopoPiFromXicToPv, float); +DECLARE_SOA_COLUMN(Chi2TopoXicToPv, chi2TopoXicToPv, float); +DECLARE_SOA_COLUMN(Chi2TopoV0ToCasc, chi2TopoV0ToCasc, float); +DECLARE_SOA_COLUMN(Chi2TopoCascToXic, chi2TopoCascToXic, float); +DECLARE_SOA_COLUMN(DecayLenXYLambda, decayLenXYLambda, float); +DECLARE_SOA_COLUMN(DecayLenXYCasc, decayLenXYCasc, float); +DECLARE_SOA_COLUMN(DecayLenXYXic, decayLenXYXic, float); +DECLARE_SOA_COLUMN(CosPaV0ToCasc, cosPaV0ToCasc, float); +DECLARE_SOA_COLUMN(CosPaV0ToPv, cosPaV0ToPv, float); +DECLARE_SOA_COLUMN(CosPaCascToXic, cosPaCascToXic, float); +DECLARE_SOA_COLUMN(CosPaCascToPv, cosPaCascToPv, float); +DECLARE_SOA_COLUMN(CosPaXicToPv, cosPaXicToPv, float); +DECLARE_SOA_COLUMN(KfRapXic, kfRapXic, float); +DECLARE_SOA_COLUMN(KfptPiFromXic, kfptPiFromXic, float); +DECLARE_SOA_COLUMN(KfptXic, kfptXic, float); +DECLARE_SOA_COLUMN(CosThetaStarPiFromXic, cosThetaStarPiFromXic, float); +DECLARE_SOA_COLUMN(CtXic, ctXic, float); +DECLARE_SOA_COLUMN(EtaXic, etaXic, float); +DECLARE_SOA_COLUMN(V0Ndf, v0Ndf, float); +DECLARE_SOA_COLUMN(CascNdf, cascNdf, float); +DECLARE_SOA_COLUMN(XicNdf, xicNdf, float); +DECLARE_SOA_COLUMN(MassV0Ndf, massV0Ndf, float); +DECLARE_SOA_COLUMN(MassCascNdf, massCascNdf, float); +DECLARE_SOA_COLUMN(V0Chi2OverNdf, v0Chi2OverNdf, float); +DECLARE_SOA_COLUMN(CascChi2OverNdf, cascChi2OverNdf, float); +DECLARE_SOA_COLUMN(XicChi2OverNdf, xicChi2OverNdf, float); +DECLARE_SOA_COLUMN(MassV0Chi2OverNdf, massV0Chi2OverNdf, float); +DECLARE_SOA_COLUMN(MassCascChi2OverNdf, massCascChi2OverNdf, float); +// from creator - MC +DECLARE_SOA_COLUMN(FlagMcMatchRec, flagMcMatchRec, int8_t); // reconstruction level +DECLARE_SOA_COLUMN(DebugMcRec, debugMcRec, int8_t); // debug flag for mis-association reconstruction level +DECLARE_SOA_COLUMN(OriginRec, originRec, int8_t); +DECLARE_SOA_COLUMN(CollisionMatched, collisionMatched, bool); +// from selector +DECLARE_SOA_COLUMN(StatusPidLambda, statusPidLambda, bool); +DECLARE_SOA_COLUMN(StatusPidCascade, statusPidCascade, bool); +DECLARE_SOA_COLUMN(StatusPidCharmBaryon, statusPidCharmBaryon, bool); +DECLARE_SOA_COLUMN(StatusInvMassLambda, statusInvMassLambda, bool); +DECLARE_SOA_COLUMN(StatusInvMassCascade, statusInvMassCascade, bool); +DECLARE_SOA_COLUMN(StatusInvMassCharmBaryon, statusInvMassCharmBaryon, bool); +DECLARE_SOA_COLUMN(ResultSelections, resultSelections, bool); +DECLARE_SOA_COLUMN(PidTpcInfoStored, pidTpcInfoStored, int); +DECLARE_SOA_COLUMN(PidTofInfoStored, pidTofInfoStored, int); +DECLARE_SOA_COLUMN(TpcNSigmaPiFromCharmBaryon, tpcNSigmaPiFromCharmBaryon, float); +DECLARE_SOA_COLUMN(TpcNSigmaPiFromCasc, tpcNSigmaPiFromCasc, float); +DECLARE_SOA_COLUMN(TpcNSigmaPiFromLambda, tpcNSigmaPiFromLambda, float); +DECLARE_SOA_COLUMN(TpcNSigmaPrFromLambda, tpcNSigmaPrFromLambda, float); +DECLARE_SOA_COLUMN(TofNSigmaPiFromCharmBaryon, tofNSigmaPiFromCharmBaryon, float); +DECLARE_SOA_COLUMN(TofNSigmaPiFromCasc, tofNSigmaPiFromCasc, float); +DECLARE_SOA_COLUMN(TofNSigmaPiFromLambda, tofNSigmaPiFromLambda, float); +DECLARE_SOA_COLUMN(TofNSigmaPrFromLambda, tofNSigmaPrFromLambda, float); +} // namespace full + +DECLARE_SOA_TABLE(HfToXiPiEvs, "AOD", "HFTOXIPIEV", + full::IsEventSel8, full::IsEventSelZ); + +DECLARE_SOA_TABLE(HfToXiPiFulls, "AOD", "HFTOXIPIFULL", + full::XPv, full::YPv, full::ZPv, full::Centrality, collision::NumContrib, collision::Chi2, + full::XDecayVtxCharmBaryon, full::YDecayVtxCharmBaryon, full::ZDecayVtxCharmBaryon, + full::XDecayVtxCascade, full::YDecayVtxCascade, full::ZDecayVtxCascade, + full::XDecayVtxV0, full::YDecayVtxV0, full::ZDecayVtxV0, + full::SignDecay, + full::CovVtxCharmBaryonXX, full::CovVtxCharmBaryonYY, full::CovVtxCharmBaryonZZ, + full::PxCharmBaryon, full::PyCharmBaryon, full::PzCharmBaryon, + full::PxCasc, full::PyCasc, full::PzCasc, + full::PxPiFromCharmBaryon, full::PyPiFromCharmBaryon, full::PzPiFromCharmBaryon, + full::PxLambda, full::PyLambda, full::PzLambda, + full::PxPiFromCasc, full::PyPiFromCasc, full::PzPiFromCasc, + full::PxPosV0Dau, full::PyPosV0Dau, full::PzPosV0Dau, + full::PxNegV0Dau, full::PyNegV0Dau, full::PzNegV0Dau, + full::ImpactParCascXY, full::ImpactParPiFromCharmBaryonXY, + full::ImpactParCascZ, full::ImpactParPiFromCharmBaryonZ, + full::ErrImpactParCascXY, full::ErrImpactParPiFromCharmBaryonXY, + full::InvMassLambda, full::InvMassCascade, full::InvMassCharmBaryon, + full::CosPAV0, full::CosPACharmBaryon, full::CosPACasc, full::CosPAXYV0, full::CosPAXYCharmBaryon, full::CosPAXYCasc, + full::CTauOmegac, full::CTauCascade, full::CTauV0, full::CTauXic, + full::EtaV0PosDau, full::EtaV0NegDau, full::EtaPiFromCasc, full::EtaPiFromCharmBaryon, + full::EtaCharmBaryon, full::EtaCascade, full::EtaV0, + full::DcaXYToPvV0Dau0, full::DcaXYToPvV0Dau1, full::DcaXYToPvCascDau, + full::DcaZToPvV0Dau0, full::DcaZToPvV0Dau1, full::DcaZToPvCascDau, + full::DcaCascDau, full::DcaV0Dau, full::DcaCharmBaryonDau, + full::DecLenCharmBaryon, full::DecLenCascade, full::DecLenV0, full::ErrorDecayLengthCharmBaryon, full::ErrorDecayLengthXYCharmBaryon, + full::NormImpParCascade, full::NormImpParPiFromCharmBar, full::NormDecayLenCharmBar, full::IsPionGlbTrkWoDca, full::PionItsNCls, + full::NTpcRowsPion, full::NTpcRowsPiFromCasc, full::NTpcRowsPosV0Dau, full::NTpcRowsNegV0Dau, + full::StatusPidLambda, full::StatusPidCascade, full::StatusPidCharmBaryon, + full::StatusInvMassLambda, full::StatusInvMassCascade, full::StatusInvMassCharmBaryon, full::ResultSelections, + full::PidTpcInfoStored, full::PidTofInfoStored, + full::TpcNSigmaPiFromCharmBaryon, full::TpcNSigmaPiFromCasc, full::TpcNSigmaPiFromLambda, full::TpcNSigmaPrFromLambda, + full::TofNSigmaPiFromCharmBaryon, full::TofNSigmaPiFromCasc, full::TofNSigmaPiFromLambda, full::TofNSigmaPrFromLambda, + full::FlagMcMatchRec, full::DebugMcRec, full::OriginRec, full::CollisionMatched); + +DECLARE_SOA_TABLE(HfToXiPiLites, "AOD", "HFTOXIPILITE", + full::XPv, full::YPv, full::ZPv, full::Centrality, collision::NumContrib, collision::Chi2, + full::XDecayVtxCharmBaryon, full::YDecayVtxCharmBaryon, full::ZDecayVtxCharmBaryon, + full::XDecayVtxCascade, full::YDecayVtxCascade, full::ZDecayVtxCascade, + full::XDecayVtxV0, full::YDecayVtxV0, full::ZDecayVtxV0, + full::SignDecay, + full::PxCharmBaryon, full::PyCharmBaryon, full::PzCharmBaryon, + full::PxPiFromCharmBaryon, full::PyPiFromCharmBaryon, full::PzPiFromCharmBaryon, + full::PxPiFromCasc, full::PyPiFromCasc, full::PzPiFromCasc, + full::PxPosV0Dau, full::PyPosV0Dau, full::PzPosV0Dau, + full::PxNegV0Dau, full::PyNegV0Dau, full::PzNegV0Dau, + full::ImpactParCascXY, full::ImpactParPiFromCharmBaryonXY, + full::ErrImpactParCascXY, full::ErrImpactParPiFromCharmBaryonXY, + full::InvMassLambda, full::InvMassCascade, full::InvMassCharmBaryon, + full::EtaV0PosDau, full::EtaV0NegDau, full::EtaPiFromCasc, full::EtaPiFromCharmBaryon, + full::DcaXYToPvV0Dau0, full::DcaXYToPvV0Dau1, full::DcaXYToPvCascDau, + full::DcaCascDau, full::DcaV0Dau, full::DcaCharmBaryonDau, + full::ErrorDecayLengthCharmBaryon, full::NormImpParCascade, full::NormImpParPiFromCharmBar, + full::IsPionGlbTrkWoDca, full::PionItsNCls, + full::NTpcRowsPion, full::NTpcRowsPiFromCasc, full::NTpcRowsPosV0Dau, full::NTpcRowsNegV0Dau, + full::StatusPidLambda, full::StatusPidCascade, full::StatusPidCharmBaryon, + full::StatusInvMassLambda, full::StatusInvMassCascade, full::StatusInvMassCharmBaryon, full::ResultSelections, + full::PidTpcInfoStored, full::PidTofInfoStored, + full::TpcNSigmaPiFromCharmBaryon, full::TpcNSigmaPiFromCasc, full::TpcNSigmaPiFromLambda, full::TpcNSigmaPrFromLambda, + full::TofNSigmaPiFromCharmBaryon, full::TofNSigmaPiFromCasc, full::TofNSigmaPiFromLambda, full::TofNSigmaPrFromLambda, + full::FlagMcMatchRec, full::OriginRec, full::CollisionMatched); + +DECLARE_SOA_TABLE(HfKfXicFulls, "AOD", "HFKFXICFULL", + full::Centrality, + // full::StatusPidLambda, full::StatusPidCascade, full::StatusPidCharmBaryon, + // full::StatusInvMassLambda, full::StatusInvMassCascade, full::StatusInvMassCharmBaryon, + full::ResultSelections, + full::TpcNSigmaPiFromCharmBaryon, full::TofNSigmaPiFromCharmBaryon, full::TpcNSigmaPiFromCasc, full::TofNSigmaPiFromCasc, + full::TpcNSigmaPiFromLambda, full::TofNSigmaPiFromLambda, full::TpcNSigmaPrFromLambda, full::TofNSigmaPrFromLambda, + full::KfDcaXYPiFromXic, full::DcaCascDau, full::DcaCharmBaryonDau, full::KfDcaXYCascToPv, + full::DcaXYToPvV0Dau0, full::DcaXYToPvV0Dau1, full::DcaXYToPvCascDau, + full::Chi2GeoV0, full::Chi2GeoCasc, full::Chi2GeoXic, + full::Chi2MassV0, full::Chi2MassCasc, + full::V0ldl, full::Cascldl, // full::Xicldl, + full::Chi2TopoV0ToPv, full::Chi2TopoCascToPv, full::Chi2TopoPiFromXicToPv, full::Chi2TopoXicToPv, + full::Chi2TopoV0ToCasc, full::Chi2TopoCascToXic, + full::DecayLenXYLambda, full::DecayLenXYCasc, full::DecayLenXYXic, + full::CosPaV0ToCasc, full::CosPaV0ToPv, full::CosPaCascToXic, full::CosPaCascToPv, // full::CosPaXicToPv, + full::InvMassLambda, full::InvMassCascade, full::InvMassCharmBaryon, + full::KfRapXic, // full::KfptPiFromXic, full::KfptXic, + full::CosThetaStarPiFromXic, full::CtXic, full::EtaXic, + full::V0Ndf, full::CascNdf, full::XicNdf, + full::MassV0Ndf, full::MassCascNdf, + full::V0Chi2OverNdf, full::CascChi2OverNdf, full::XicChi2OverNdf, + full::MassV0Chi2OverNdf, full::MassCascChi2OverNdf, + full::FlagMcMatchRec, full::DebugMcRec, full::OriginRec, full::CollisionMatched); + +} // namespace o2::aod + +/// Writes the full information in an output TTree +struct HfTreeCreatorToXiPiQa { + + Produces rowCandidateFull; + Produces rowCandidateLite; + Produces rowKfCandidate; + Produces rowEv; + + Configurable zPvCut{"zPvCut", 10., "Cut on absolute value of primary vertex z coordinate"}; + + using MyTrackTable = soa::Join; + using MyEventTable = soa::Join; + using MyEventTableWithFT0C = soa::Join; + using MyEventTableWithFT0M = soa::Join; + using MyEventTableWithNTracksPV = soa::Join; + + void init(InitContext const&) + { + if ((doprocessMcLiteXic0 && doprocessMcLiteOmegac0) || (doprocessMcFullXic0 && doprocessMcFullOmegac0)) { + LOGF(fatal, "Both Xic0 and Omegac0 MC processes enabled, please choose ONLY one!"); + } + } + + ////////////////////////////////////////////////////// + // // + // Fill functions to fill in the tables // + // // + ////////////////////////////////////////////////////// + + template + void fillEvent(const T& collision, float cutZPv) + { + rowEv(collision.sel8(), std::abs(collision.posZ()) < cutZPv); + } + + template + void fillCandidate(const T& candidate, int8_t flagMc, int8_t debugMc, int8_t originMc, bool collisionMatched) + { + + // Save all candidate information + float centrality = -999.f; + if constexpr (useCentrality) { + auto const& collision = candidate.template collision_as(); + centrality = o2::hf_centrality::getCentralityColl(collision); + } + + rowCandidateFull( + candidate.xPv(), + candidate.yPv(), + candidate.zPv(), + centrality, + candidate.template collision_as().numContrib(), + candidate.template collision_as().chi2(), + candidate.xDecayVtxCharmBaryon(), + candidate.yDecayVtxCharmBaryon(), + candidate.zDecayVtxCharmBaryon(), + candidate.xDecayVtxCascade(), + candidate.yDecayVtxCascade(), + candidate.zDecayVtxCascade(), + candidate.xDecayVtxV0(), + candidate.yDecayVtxV0(), + candidate.zDecayVtxV0(), + candidate.signDecay(), + candidate.covVtxCharmBaryon0(), + candidate.covVtxCharmBaryon3(), + candidate.covVtxCharmBaryon5(), + candidate.pxCharmBaryon(), + candidate.pyCharmBaryon(), + candidate.pzCharmBaryon(), + candidate.pxCasc(), + candidate.pyCasc(), + candidate.pzCasc(), + candidate.pxBachFromCharmBaryon(), + candidate.pyBachFromCharmBaryon(), + candidate.pzBachFromCharmBaryon(), + candidate.pxLambda(), + candidate.pyLambda(), + candidate.pzLambda(), + candidate.pxBachFromCasc(), + candidate.pyBachFromCasc(), + candidate.pzBachFromCasc(), + candidate.pxPosV0Dau(), + candidate.pyPosV0Dau(), + candidate.pzPosV0Dau(), + candidate.pxNegV0Dau(), + candidate.pyNegV0Dau(), + candidate.pzNegV0Dau(), + candidate.impactParCascXY(), + candidate.impactParBachFromCharmBaryonXY(), + candidate.impactParCascZ(), + candidate.impactParBachFromCharmBaryonZ(), + candidate.errImpactParCascXY(), + candidate.errImpactParBachFromCharmBaryonXY(), + candidate.invMassLambda(), + candidate.invMassCascade(), + candidate.invMassCharmBaryon(), + candidate.cosPAV0(), + candidate.cosPACharmBaryon(), + candidate.cosPACasc(), + candidate.cosPAXYV0(), + candidate.cosPAXYCharmBaryon(), + candidate.cosPAXYCasc(), + candidate.cTauOmegac(), + candidate.cTauCascade(), + candidate.cTauV0(), + candidate.cTauXic(), + candidate.etaV0PosDau(), + candidate.etaV0NegDau(), + candidate.etaBachFromCasc(), + candidate.etaBachFromCharmBaryon(), + candidate.etaCharmBaryon(), + candidate.etaCascade(), + candidate.etaV0(), + candidate.dcaXYToPvV0Dau0(), + candidate.dcaXYToPvV0Dau1(), + candidate.dcaXYToPvCascDau(), + candidate.dcaZToPvV0Dau0(), + candidate.dcaZToPvV0Dau1(), + candidate.dcaZToPvCascDau(), + candidate.dcaCascDau(), + candidate.dcaV0Dau(), + candidate.dcaCharmBaryonDau(), + candidate.decLenCharmBaryon(), + candidate.decLenCascade(), + candidate.decLenV0(), + candidate.errorDecayLengthCharmBaryon(), + candidate.errorDecayLengthXYCharmBaryon(), + candidate.impactParCascXY() / candidate.errImpactParCascXY(), + candidate.impactParBachFromCharmBaryonXY() / candidate.errImpactParBachFromCharmBaryonXY(), + candidate.decLenCharmBaryon() / candidate.errorDecayLengthCharmBaryon(), + candidate.template bachelorFromCharmBaryon_as().isGlobalTrackWoDCA(), + candidate.template bachelorFromCharmBaryon_as().itsNCls(), + candidate.template bachelorFromCharmBaryon_as().tpcNClsCrossedRows(), + candidate.template bachelor_as().tpcNClsCrossedRows(), + candidate.template posTrack_as().tpcNClsCrossedRows(), + candidate.template negTrack_as().tpcNClsCrossedRows(), + candidate.statusPidLambda(), + candidate.statusPidCascade(), + candidate.statusPidCharmBaryon(), + candidate.statusInvMassLambda(), + candidate.statusInvMassCascade(), + candidate.statusInvMassCharmBaryon(), + candidate.resultSelections(), + candidate.pidTpcInfoStored(), + candidate.pidTofInfoStored(), + candidate.tpcNSigmaPiFromCharmBaryon(), + candidate.tpcNSigmaPiFromCasc(), + candidate.tpcNSigmaPiFromLambda(), + candidate.tpcNSigmaPrFromLambda(), + candidate.tofNSigmaPiFromCharmBaryon(), + candidate.tofNSigmaPiFromCasc(), + candidate.tofNSigmaPiFromLambda(), + candidate.tofNSigmaPrFromLambda(), + flagMc, + debugMc, + originMc, + collisionMatched); + } + + template + void fillCandidateLite(const T& candidate, int8_t flagMc, int8_t originMc, bool collisionMatched) + { + float centrality = -999.f; + if constexpr (useCentrality) { + auto const& collision = candidate.template collision_as(); + centrality = o2::hf_centrality::getCentralityColl(collision); + } + + rowCandidateLite( + candidate.xPv(), + candidate.yPv(), + candidate.zPv(), + centrality, + candidate.template collision_as().numContrib(), + candidate.template collision_as().chi2(), + candidate.xDecayVtxCharmBaryon(), + candidate.yDecayVtxCharmBaryon(), + candidate.zDecayVtxCharmBaryon(), + candidate.xDecayVtxCascade(), + candidate.yDecayVtxCascade(), + candidate.zDecayVtxCascade(), + candidate.xDecayVtxV0(), + candidate.yDecayVtxV0(), + candidate.zDecayVtxV0(), + candidate.signDecay(), + candidate.pxCharmBaryon(), + candidate.pyCharmBaryon(), + candidate.pzCharmBaryon(), + candidate.pxBachFromCharmBaryon(), + candidate.pyBachFromCharmBaryon(), + candidate.pzBachFromCharmBaryon(), + candidate.pxBachFromCasc(), + candidate.pyBachFromCasc(), + candidate.pzBachFromCasc(), + candidate.pxPosV0Dau(), + candidate.pyPosV0Dau(), + candidate.pzPosV0Dau(), + candidate.pxNegV0Dau(), + candidate.pyNegV0Dau(), + candidate.pzNegV0Dau(), + candidate.impactParCascXY(), + candidate.impactParBachFromCharmBaryonXY(), + candidate.errImpactParCascXY(), + candidate.errImpactParBachFromCharmBaryonXY(), + candidate.invMassLambda(), + candidate.invMassCascade(), + candidate.invMassCharmBaryon(), + candidate.etaV0PosDau(), + candidate.etaV0NegDau(), + candidate.etaBachFromCasc(), + candidate.etaBachFromCharmBaryon(), + candidate.dcaXYToPvV0Dau0(), + candidate.dcaXYToPvV0Dau1(), + candidate.dcaXYToPvCascDau(), + candidate.dcaCascDau(), + candidate.dcaV0Dau(), + candidate.dcaCharmBaryonDau(), + candidate.errorDecayLengthCharmBaryon(), + candidate.impactParCascXY() / candidate.errImpactParCascXY(), + candidate.impactParBachFromCharmBaryonXY() / candidate.errImpactParBachFromCharmBaryonXY(), + candidate.template bachelorFromCharmBaryon_as().isGlobalTrackWoDCA(), + candidate.template bachelorFromCharmBaryon_as().itsNCls(), + candidate.template bachelorFromCharmBaryon_as().tpcNClsCrossedRows(), + candidate.template bachelor_as().tpcNClsCrossedRows(), + candidate.template posTrack_as().tpcNClsCrossedRows(), + candidate.template negTrack_as().tpcNClsCrossedRows(), + candidate.statusPidLambda(), + candidate.statusPidCascade(), + candidate.statusPidCharmBaryon(), + candidate.statusInvMassLambda(), + candidate.statusInvMassCascade(), + candidate.statusInvMassCharmBaryon(), + candidate.resultSelections(), + candidate.pidTpcInfoStored(), + candidate.pidTofInfoStored(), + candidate.tpcNSigmaPiFromCharmBaryon(), + candidate.tpcNSigmaPiFromCasc(), + candidate.tpcNSigmaPiFromLambda(), + candidate.tpcNSigmaPrFromLambda(), + candidate.tofNSigmaPiFromCharmBaryon(), + candidate.tofNSigmaPiFromCasc(), + candidate.tofNSigmaPiFromLambda(), + candidate.tofNSigmaPrFromLambda(), + flagMc, + originMc, + collisionMatched); + } + + template + void fillKfCandidate(const T& candidate, int8_t flagMc, int8_t debugMc, int8_t originMc, bool collisionMatched) + { + float centrality = -999.f; + if constexpr (useCentrality) { + auto const& collision = candidate.template collision_as(); + centrality = o2::hf_centrality::getCentralityColl(collision); + } + + rowKfCandidate( + centrality, + // candidate.statusPidLambda(), + // candidate.statusPidCascade(), + // candidate.statusPidCharmBaryon(), + // candidate.statusInvMassLambda(), + // candidate.statusInvMassCascade(), + // candidate.statusInvMassCharmBaryon(), + candidate.resultSelections(), + candidate.tpcNSigmaPiFromCharmBaryon(), + candidate.tofNSigmaPiFromCharmBaryon(), + candidate.tpcNSigmaPiFromCasc(), + candidate.tofNSigmaPiFromCasc(), + candidate.tpcNSigmaPiFromLambda(), + candidate.tofNSigmaPiFromLambda(), + candidate.tpcNSigmaPrFromLambda(), + candidate.tofNSigmaPrFromLambda(), + candidate.kfDcaXYPiFromXic(), + candidate.dcaCascDau(), + candidate.dcaCharmBaryonDau(), + candidate.kfDcaXYCascToPv(), + candidate.dcaXYToPvV0Dau0(), + candidate.dcaXYToPvV0Dau1(), + candidate.dcaXYToPvCascDau(), + candidate.chi2GeoV0(), + candidate.chi2GeoCasc(), + candidate.chi2GeoXic(), + candidate.chi2MassV0(), + candidate.chi2MassCasc(), + candidate.v0ldl(), + candidate.cascldl(), + // candidate.xicldl(), + candidate.chi2TopoV0ToPv(), + candidate.chi2TopoCascToPv(), + candidate.chi2TopoPiFromXicToPv(), + candidate.chi2TopoXicToPv(), + candidate.chi2TopoV0ToCasc(), + candidate.chi2TopoCascToXic(), + candidate.decayLenXYLambda(), + candidate.decayLenXYCasc(), + candidate.decayLenXYXic(), + candidate.cosPaV0ToCasc(), + candidate.cosPAV0(), + candidate.cosPaCascToXic(), + candidate.cosPACasc(), + // candidate.cosPACharmBaryon(), + candidate.invMassLambda(), + candidate.invMassCascade(), + candidate.invMassCharmBaryon(), + candidate.kfRapXic(), + // candidate.kfptPiFromXic(), + // candidate.kfptXic(), + candidate.cosThetaStarPiFromXic(), + candidate.cTauXic(), + candidate.etaCharmBaryon(), + candidate.v0Ndf(), + candidate.cascNdf(), + candidate.xicNdf(), + candidate.massV0Ndf(), + candidate.massCascNdf(), + candidate.v0Chi2OverNdf(), + candidate.cascChi2OverNdf(), + candidate.xicChi2OverNdf(), + candidate.massV0Chi2OverNdf(), + candidate.massCascChi2OverNdf(), + flagMc, + debugMc, + originMc, + collisionMatched); + } + + //////////////////////////////////// + // // + // Process functions // + // // + //////////////////////////////////// + + //*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*// + //*~~~~~~~Data with DCAFitter~~~~~~~~*// + //*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*// + + void processDataFull(MyEventTable const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateFull.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, -7, -7, RecoDecay::OriginType::None, false); + } + } + + void processDataLite(MyEventTable const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateLite.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidateLite(candidate, -7, RecoDecay::OriginType::None, false); + } + } + + void processDataLiteWithFT0M(MyEventTableWithFT0M const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateLite.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidateLite(candidate, -7, RecoDecay::OriginType::None, false); + } + } + + void processDataLiteWithFT0C(MyEventTableWithFT0C const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateLite.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidateLite(candidate, -7, RecoDecay::OriginType::None, false); + } + } + + void processDataLiteWithNTracksPV(MyEventTableWithNTracksPV const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateLite.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidateLite(candidate, -7, RecoDecay::OriginType::None, false); + } + } + + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processDataFull, "Process data with full information w/o centrality", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processDataLite, "Process data and produce lite table version", true); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processDataLiteWithFT0M, "Process data and produce lite table version with FT0M", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processDataLiteWithFT0C, "Process data and produce lite table version with FT0C", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processDataLiteWithNTracksPV, "Process data and produce lite table version with NTracksPV", false); + + //*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*// + //*~~~~~~~Data with KFParticle~~~~~~~~*// + //*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*// + void processKfData(MyEventTable const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowKfCandidate.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillKfCandidate(candidate, -7, -7, RecoDecay::OriginType::None, false); + } + } + + void processKfDataWithFT0M(MyEventTableWithFT0M const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowKfCandidate.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillKfCandidate(candidate, -7, -7, RecoDecay::OriginType::None, false); + } + } + + void processKfDataWithFT0C(MyEventTableWithFT0C const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowKfCandidate.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillKfCandidate(candidate, -7, -7, RecoDecay::OriginType::None, false); + } + } + + void processKfDataWithNTracksPV(MyEventTableWithNTracksPV const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowKfCandidate.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillKfCandidate(candidate, -7, -7, RecoDecay::OriginType::None, false); + } + } + + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processKfData, "Process KF data, no cent", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processKfDataWithFT0M, "Process KF data, with FT0M", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processKfDataWithFT0C, "Process KF data, with FT0C", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processKfDataWithNTracksPV, "Process KF data, with NTracksPV", false); + //*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*// + //*~~~~~~~MC with DCAFitter~~~~~~~~*// + //*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*// + + void processMcFullXic0(MyEventTable const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateFull.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, candidate.flagMcMatchRec(), candidate.debugMcRec(), candidate.originMcRec(), candidate.collisionMatched()); + } + } + + void processMcFullOmegac0(MyEventTable const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateFull.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidate(candidate, candidate.flagMcMatchRec(), candidate.debugMcRec(), candidate.originMcRec(), candidate.collisionMatched()); + } + } + + void processMcLiteXic0(MyEventTable const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateLite.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidateLite(candidate, candidate.flagMcMatchRec(), candidate.originMcRec(), candidate.collisionMatched()); + } + } + + void processMcLiteXic0WithFT0C(MyEventTableWithFT0C const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateLite.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidateLite(candidate, candidate.flagMcMatchRec(), candidate.originMcRec(), candidate.collisionMatched()); + } + } + + void processMcLiteXic0WithFT0M(MyEventTableWithFT0M const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateLite.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidateLite(candidate, candidate.flagMcMatchRec(), candidate.originMcRec(), candidate.collisionMatched()); + } + } + + void processMcLiteXic0WithNTracksPV(MyEventTableWithNTracksPV const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateLite.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidateLite(candidate, candidate.flagMcMatchRec(), candidate.originMcRec(), candidate.collisionMatched()); + } + } + + void processMcLiteOmegac0(MyEventTable const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowCandidateLite.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillCandidateLite(candidate, candidate.flagMcMatchRec(), candidate.originMcRec(), candidate.collisionMatched()); + } + } + + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processMcFullXic0, "Process MC with full information for xic0 w/o centrality", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processMcFullOmegac0, "Process MC with full information for omegac0", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processMcLiteXic0, "Process MC and produce lite table version for xic0", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processMcLiteXic0WithFT0C, "Process MC and produce lite table version for Xic0 with FT0C", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processMcLiteXic0WithFT0M, "Process MC and produce lite table version for Xic0 with FT0M", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processMcLiteXic0WithNTracksPV, "Process MC and produce lite table version for Xic0 with NTracksPV", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processMcLiteOmegac0, "Process MC and produce lite table version for omegac0", false); + + //*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*// + //*~~~~~~~MC with KFParticle~~~~~~~~*// + //*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*// + void processKfMcXic0(MyEventTable const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowKfCandidate.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillKfCandidate(candidate, candidate.flagMcMatchRec(), candidate.debugMcRec(), candidate.originMcRec(), candidate.collisionMatched()); + } + } + + void processKfMcXic0WithFT0C(MyEventTableWithFT0C const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowKfCandidate.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillKfCandidate(candidate, candidate.flagMcMatchRec(), candidate.debugMcRec(), candidate.originMcRec(), candidate.collisionMatched()); + } + } + + void processKfMcXic0WithFT0M(MyEventTableWithFT0M const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate properties + rowKfCandidate.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillKfCandidate(candidate, candidate.flagMcMatchRec(), candidate.debugMcRec(), candidate.originMcRec(), candidate.collisionMatched()); + } + } + + void processKfMcXic0WithNTracksPV(MyEventTableWithNTracksPV const& collisions, MyTrackTable const&, + soa::Join const& candidates) + { + // Filling event properties + rowEv.reserve(collisions.size()); + for (const auto& collision : collisions) { + fillEvent(collision, zPvCut); + } + + // Filling candidate table + rowKfCandidate.reserve(candidates.size()); + for (const auto& candidate : candidates) { + fillKfCandidate(candidate, candidate.flagMcMatchRec(), candidate.debugMcRec(), candidate.originMcRec(), candidate.collisionMatched()); + } + } + + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processKfMcXic0, "Process MC with information for xic0", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processKfMcXic0WithFT0C, "Process MC with information for xic0 at FT0C", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processKfMcXic0WithFT0M, "Process MC with information for xic0 at FT0M", false); + PROCESS_SWITCH(HfTreeCreatorToXiPiQa, processKfMcXic0WithNTracksPV, "Process MC with information for xic0 at Ntrack", false); +}; // end of struct + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +}