diff --git a/PWGLF/Tasks/Strangeness/CMakeLists.txt b/PWGLF/Tasks/Strangeness/CMakeLists.txt index 31a0087b3fa..5c0fecbfeb8 100644 --- a/PWGLF/Tasks/Strangeness/CMakeLists.txt +++ b/PWGLF/Tasks/Strangeness/CMakeLists.txt @@ -174,4 +174,4 @@ o2physics_add_dpl_workflow(cascadeanalysislightions o2physics_add_dpl_workflow(strangecasctrack SOURCES strangecasctrack.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2Physics::EventFilteringUtils - COMPONENT_NAME myo2) + COMPONENT_NAME Analysis) diff --git a/PWGLF/Tasks/Strangeness/strangecasctrack.cxx b/PWGLF/Tasks/Strangeness/strangecasctrack.cxx index 609537fcf36..f4f2369cee5 100644 --- a/PWGLF/Tasks/Strangeness/strangecasctrack.cxx +++ b/PWGLF/Tasks/Strangeness/strangecasctrack.cxx @@ -43,20 +43,20 @@ using namespace o2::framework; using namespace o2::framework::expressions; // tables for derived data -using DerCollisionWMult = soa::Join::iterator; -using DerCascDatas = soa::Join; +using DerCollisionWMults = soa::Join; +using DerCascDatas = soa::Join; using DerTraCascDatas = soa::Join; // tables for derived MC using DerMCGenCascades = soa::Join; -using DerMCRecCollision = soa::Join::iterator; -using DerMCRecCascDatas = soa::Join; +using DerMCRecCollisions = soa::Join; +using DerMCRecCascDatas = soa::Join; using DerMCRecTraCascDatas = soa::Join; // tables for PID selection using DauTracks = soa::Join; -struct StrangeCascTrack { +struct strangecasctrack { Service ccdb; Service pdgDB; @@ -66,9 +66,12 @@ struct StrangeCascTrack { HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; // subprocess switches: - Configurable doProcessPP{"doProcessPP", true, "true for pp, false for PbPb and OO"}; + Configurable doProcesspp{"doProcesspp", true, "true for pp"}; + Configurable doProcessPbPb{"doProcessPbPb", false, "true for PbPb"}; + Configurable doProcessOO{"doProcessOO", false, "true for OO"}; + Configurable doProcesspO{"doProcesspO", false, "true for pO"}; // selections - Configurable doApplyCuts{"doApplyCuts", true, "apply cuts"}; // dca for filtering data primaries + Configurable doApplyCuts{"doApplyCuts", true, "apply cuts"}; // general cascade cuts - dca, cosPA etc. Configurable doApplyTPCPID{"doApplyTPCPID", true, "apply tpc pid to dau tracks"}; Configurable doApplyTOFPID{"doApplyTOFPID", true, "apply tof pid to dau tracks"}; Configurable doCompetingMassRej{"doCompetingMassRej", true, "competing mass rejection for omegas"}; @@ -78,52 +81,69 @@ struct StrangeCascTrack { Configurable doApplyPurity{"doApplyPurity", false, "apply purity correction"}; Configurable doPropagatePurity{"doPropagatePurity", false, "apply purity propagation"}; Configurable ccdburl{"ccdburl", "http://alice-ccdb.cern.ch", "url of the ccdb repository to use"}; - Configurable efficiencyCCDBPath{"efficiencyCCDBPath", "GLO/Config/GeometryAligned", "Path of the efficiency corrections"}; - // mc settings - Configurable doFillTruth{"doFillTruth", false, "require MC truth for reco"}; + Configurable efficiencyCCDBPath_pp{"efficiencyCCDBPath_pp", "Users/y/yparovia/LHC24f4d/Efficiency", "Path of the efficiency corrections"}; + Configurable efficiencyCCDBPath_PbPb{"efficiencyCCDBPath_PbPb", "Users/y/yparovia/LHC25f3/Efficiency", "Path of the efficiency corrections"}; + Configurable efficiencyCCDBPath_OO{"efficiencyCCDBPath_OO", "Users/y/yparovia/LHC25h3/Efficiency", "Path of the efficiency corrections"}; + Configurable efficiencyCCDBPath_pO{"efficiencyCCDBPath_pO", "Users/y/yparovia/LHC25h2/Efficiency", "Path of the efficiency corrections"}; // event and dau track selection struct : ConfigurableGroup { + Configurable cutZVertex{"cutZVertex", 10.0f, "max Z-vertex position"}; Configurable cutDCAtoPVxy{"cutDCAtoPVxy", 0.02f, "max cascade dca to PV in xy"}; Configurable cutDCAtoPVz{"cutDCAtoPVz", 0.02f, "max cascade dca to PV in z"}; - Configurable compMassRej{"compMassRej", 0.008, "Competing mass rejection"}; + Configurable compMassRej{"compMassRej", 0.008f, "Competing mass rejection"}; + Configurable cutV0CosPA{"cutV0CosPA", 0.97f, "max V0 cosPA"}; + Configurable cutBachCosPA{"cutBachCosPA", 0.97f, "max Bachelor cosPA"}; // TPC PID selection - Configurable NSigmaTPCPion{"NSigmaTPCPion", 4, "NSigmaTPCPion"}; - Configurable NSigmaTPCKaon{"NSigmaTPCKaon", 4, "NSigmaTPCKaon"}; - Configurable NSigmaTPCProton{"NSigmaTPCProton", 4, "NSigmaTPCProton"}; + Configurable nSigmaTPCPion{"nSigmaTPCPion", 4, "NSigmaTPCPion"}; + Configurable nSigmaTPCKaon{"nSigmaTPCKaon", 4, "NSigmaTPCKaon"}; + Configurable nSigmaTPCProton{"nSigmaTPCProton", 4, "NSigmaTPCProton"}; // TOF PID selection - Configurable NSigmaTOFPion{"NSigmaTOFPion", 3, "NSigmaTOFPion"}; - Configurable NSigmaTOFKaon{"NSigmaTOFKaon", 3, "NSigmaTOFKaon"}; - Configurable NSigmaTOFProton{"NSigmaTOFProton", 3, "NSigmaTOFProton"}; + Configurable nSigmaTOFXi{"nSigmaTOFXi", 3, "nSigmaTOFXi"}; + Configurable nSigmaTOFOmega{"nSigmaTOFOmega", 3, "nSigmaTOFOmega"}; } selCuts; // axes struct : ConfigurableGroup { - ConfigurableAxis axisPhi{"Phi", {72, 0, TwoPI}, "#phi"}; - ConfigurableAxis axisEta{"Eta", {102, -2.01, 2.01}, "#eta"}; - ConfigurableAxis axisDCAxy{"DCA to xy plane", {500, 0., 0.5}, "cm"}; - ConfigurableAxis axisDCAz{"DCA to z plane", {500, 0., 0.5}, "cm"}; + ConfigurableAxis axisPhi{"axisPhi", {72, 0, TwoPI}, "#phi"}; + ConfigurableAxis axisEta{"axisEta", {102, -2.01, 2.01}, "#eta"}; + ConfigurableAxis axisDCAxy{"axisDCAxy", {500, 0., 0.5}, "cm"}; + ConfigurableAxis axisDCAz{"axisDCAz", {500, 0., 0.5}, "cm"}; ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.0, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0, 10.0}, "p_{T} (GeV/c)"}; - ConfigurableAxis axisMult{"axisMult", {VARIABLE_WIDTH, 0.0f, 0.01f, 1.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 70.0f, 100.0f}, "FT0 mult %"}; + ConfigurableAxis axisMult{"axisMult", {VARIABLE_WIDTH, 0.0f, 5.0, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 70.0f, 100.0f}, "FT0 mult %"}; ConfigurableAxis axisOmegaMass{"axisOmegaMass", {2000, 1.6, 1.8}, "#Omega M_{inv} (GeV/c^{2})"}; ConfigurableAxis axisXiMass{"axisXiMass", {2000, 1.2, 1.4}, "#Xi M_{inv} (GeV/c^{2})"}; } axesConfig; - // // Filters events - // if (doFilterEvents) { - // Filter eventFilter = (o2::aod::evsel::sel8 == true); - // Filter posZFilter = (nabs(o2::aod::collision::posZ) < cutzvertex); - // Filter posZFilterMC = (nabs(o2::aod::mccollision::posZ) < cutzvertex); - // } + // Filters events + // Filter eventFilter = (o2::aod::evsel::sel8 == true); + // Filter posZFilter = (nabs(o2::aod::collision::posZ) < selCuts.cutZVertex); + // Filter posZFilterMC = (nabs(o2::aod::mccollision::posZ) < selCuts.cutZVertex); // cascade reconstruction types static constexpr std::string_view kTypeNames[] = {"Standard", "Tracked"}; // for efficiency and purity corrections - TH2F* hEfficiency; - TH2F* hEfficiencyUncertainty; - TH2F* hPurity; - TH2F* hPurityUncertainty; + TH2F* hEfficiencyOmegaStd; + TH2F* hEfficiencyOmegaTra; + TH2F* hEfficiencyXiStd; + TH2F* hEfficiencyXiTra; + + TH2F* hEfficiencyErrOmegaStd; + TH2F* hEfficiencyErrOmegaTra; + TH2F* hEfficiencyErrXiStd; + TH2F* hEfficiencyErrXiTra; + + TH2F* hPurityOmegaStd; + TH2F* hPurityOmegaTra; + TH2F* hPurityXiStd; + TH2F* hPurityXiTra; + + TH2F* hPurityErrOmegaStd; + TH2F* hPurityErrOmegaTra; + TH2F* hPurityErrXiStd; + TH2F* hPurityErrXiTra; + int mRunNumber; // loads efficiencies and purities void initEfficiencyFromCCDB(int64_t runNumber, int64_t timestamp) @@ -135,19 +155,46 @@ struct StrangeCascTrack { LOG(info) << "Loading efficiencies and purities from CCDB for run " << mRunNumber << " now..."; auto timeStamp = timestamp; + std::string efficiencyCCDBPath = [&]() { + if (doProcesspp) { + return efficiencyCCDBPath_pp; + } else if (doProcesspO) { + return efficiencyCCDBPath_pO; + } else if (doProcessPbPb) { + return efficiencyCCDBPath_PbPb; + } + return efficiencyCCDBPath_OO; + }(); + TList* listEfficiencies = ccdb->getForTimeStamp(efficiencyCCDBPath, timeStamp); if (!listEfficiencies) { LOG(fatal) << "Problem getting TList object with efficiencies and purities!"; } - hEfficiency = static_cast(listEfficiencies->FindObject("hEfficiency")); - hPurity = static_cast(listEfficiencies->FindObject("hPurity")); - hEfficiencyUncertainty = static_cast(listEfficiencies->FindObject("hEfficiencyUncertainty")); - hPurityUncertainty = static_cast(listEfficiencies->FindObject("hPurityUncertainty")); - if (doPropagateEfficiency && !hEfficiencyUncertainty) + hEfficiencyOmegaStd = static_cast(listEfficiencies->FindObject("Eff_Omega_Standard")); + hEfficiencyOmegaTra = static_cast(listEfficiencies->FindObject("Eff_Omega_Tracked")); + hEfficiencyXiStd = static_cast(listEfficiencies->FindObject("Eff_Xi_Standard")); + hEfficiencyXiTra = static_cast(listEfficiencies->FindObject("Eff_Xi_Tracked")); + + hEfficiencyErrOmegaStd = static_cast(listEfficiencies->FindObject("EffErr_Omega_Standard")); + hEfficiencyErrOmegaTra = static_cast(listEfficiencies->FindObject("EffErr_Omega_Tracked")); + hEfficiencyErrXiStd = static_cast(listEfficiencies->FindObject("EffErr_Xi_Standard")); + hEfficiencyErrXiTra = static_cast(listEfficiencies->FindObject("EffErr_Xi_Tracked")); + + hPurityOmegaStd = static_cast(listEfficiencies->FindObject("Pur_Omega_Standard")); + hPurityOmegaTra = static_cast(listEfficiencies->FindObject("Pur_Omega_Tracked")); + hPurityXiStd = static_cast(listEfficiencies->FindObject("Pur_Xi_Standard")); + hPurityXiTra = static_cast(listEfficiencies->FindObject("Pur_Xi_Tracked")); + + hPurityErrOmegaStd = static_cast(listEfficiencies->FindObject("PurErr_Omega_Standard")); + hPurityErrOmegaTra = static_cast(listEfficiencies->FindObject("PurErr_Omega_Tracked")); + hPurityErrXiStd = static_cast(listEfficiencies->FindObject("PurErr_Xi_Standard")); + hPurityErrXiTra = static_cast(listEfficiencies->FindObject("PurErr_Xi_Tracked")); + + if (doPropagateEfficiency && (!hEfficiencyErrOmegaStd || !hEfficiencyErrOmegaTra || !hEfficiencyErrXiStd || !hEfficiencyErrXiTra)) LOG(fatal) << "Problem getting hEfficiencyUncertainty!"; - if (doPropagatePurity && !hPurityUncertainty) + if (doPropagatePurity && (!hPurityErrOmegaStd || !hPurityErrOmegaTra || !hPurityErrXiStd || !hPurityErrXiTra)) LOG(fatal) << "Problem getting hPurityUncertainty!"; LOG(info) << "Efficiencies and purities now loaded for " << mRunNumber; } @@ -156,7 +203,7 @@ struct StrangeCascTrack { void fillEvents(TEvent const& collision) { histos.fill(HIST("Events/EvCounter"), 0.5); - double mult = doProcessPP ? collision.centFT0M() : collision.centFT0C(); + double mult = (doProcesspp || doProcesspO) ? collision.centFT0M() : collision.centFT0C(); histos.fill(HIST("Events/Mult"), mult); double pvx = collision.posX(); double pvy = collision.posY(); @@ -166,13 +213,19 @@ struct StrangeCascTrack { histos.fill(HIST("Events/PVz"), pvz); } // checks general selection criteria - template - bool isValidCasc(TCascade cascade) + template + bool isValidCasc(TEvent collision, TCascade cascade) { + if (std::abs(collision.posZ()) > selCuts.cutZVertex) + return false; if (cascade.dcaXYCascToPV() > selCuts.cutDCAtoPVxy) return false; if (cascade.dcaZCascToPV() > selCuts.cutDCAtoPVz) return false; + if (cascade.v0cosPA(collision.posX(), collision.posY(), collision.posZ()) < selCuts.cutV0CosPA) + return false; + if (cascade.bachBaryonCosPA() < selCuts.cutBachCosPA) + return false; return true; } // checks TPC PID of dau tracks @@ -182,96 +235,39 @@ struct StrangeCascTrack { const auto& posDaughterTrackCasc = cascade.template posTrackExtra_as(); const auto& negDaughterTrackCasc = cascade.template negTrackExtra_as(); if (cascade.sign() < 0) { - if (std::abs(posDaughterTrackCasc.tpcNSigmaPr()) > selCuts.NSigmaTPCProton) { + if (std::abs(posDaughterTrackCasc.tpcNSigmaPr()) > selCuts.nSigmaTPCProton) { return false; } - if (std::abs(negDaughterTrackCasc.tpcNSigmaPi()) > selCuts.NSigmaTPCPion) { + if (std::abs(negDaughterTrackCasc.tpcNSigmaPi()) > selCuts.nSigmaTPCPion) { return false; } } else { - if (std::abs(negDaughterTrackCasc.tpcNSigmaPr()) > selCuts.NSigmaTPCProton) { + if (std::abs(negDaughterTrackCasc.tpcNSigmaPr()) > selCuts.nSigmaTPCProton) { return false; } - if (std::abs(posDaughterTrackCasc.tpcNSigmaPi()) > selCuts.NSigmaTPCPion) { + if (std::abs(posDaughterTrackCasc.tpcNSigmaPi()) > selCuts.nSigmaTPCPion) { return false; } } return true; } // checks TOF PID of dau tracks - // template - // bool passesTOF(TCascade cascade, TString particle) - // { - // return true; - // // const auto& bachDaughterTrackCasc = cascade.bachTrackExtra_as(); - // // const auto& posDaughterTrackCasc = cascade.posTrackExtra_as(); - // // const auto& negDaughterTrackCasc = cascade.negTrackExtra_as(); - // // bool xiPassTOFSelection = true; - // // bool omegaPassTOFSelection = true; - // // if (cascade.sign() < 0) { - // // if (posDaughterTrackCasc.hasTOF()) { - // // if (std::abs(cascade.tofNSigmaXiLaPr()) > selCuts.NSigmaTOFProton) { - // // xiPassTOFSelection &= false; - // // } - // // if (std::abs(cascade.tofNSigmaOmLaPr()) > selCuts.NSigmaTOFProton) { - // // omegaPassTOFSelection &= false; - // // } - // // } - // // if (negDaughterTrackCasc.hasTOF()) { - // // if (std::abs(cascade.tofNSigmaXiLaPi()) > selCuts.NSigmaTOFPion) { - // // xiPassTOFSelection &= false; - // // } - // // if (std::abs(cascade.tofNSigmaOmLaPi()) > selCuts.NSigmaTOFPion) { - // // omegaPassTOFSelection &= false; - // // } - // // } - // // } else { - // // if (posDaughterTrackCasc.hasTOF()) { - // // if (std::abs(cascade.tofNSigmaXiLaPi()) > selCuts.NSigmaTOFPion) { - // // xiPassTOFSelection &= false; - // // } - // // if (std::abs(cascade.tofNSigmaOmLaPi()) > selCuts.NSigmaTOFPion) { - // // omegaPassTOFSelection &= false; - // // } - // // } - // // if (negDaughterTrackCasc.hasTOF()) { - // // if (std::abs(cascade.tofNSigmaXiLaPr()) > selCuts.NSigmaTOFProton) { - // // xiPassTOFSelection &= false; - // // } - // // if (std::abs(cascade.tofNSigmaOmLaPr()) > selCuts.NSigmaTOFProton) { - // // omegaPassTOFSelection &= false; - // // } - // // } - // // } - // - // // if (bachDaughterTrackCasc.hasTOF()) { - // // if (std::abs(cascade.tofNSigmaXiPi()) > selCuts.NSigmaTOFPion) { - // // xiPassTOFSelection &= false; - // // } - // // if (std::abs(cascade.tofNSigmaOmKa()) > selCuts.NSigmaTOFKaon) { - // // omegaPassTOFSelection &= false; - // // } - // // } - // - // // if (bachDaughterTrackCasc.hasTOF()) { - // // if (std::abs(cascade.tofNSigmaXiPi()) > selCuts.NSigmaTOFPion) { - // // xiPassTOFSelection &= false; - // // } - // // if (std::abs(cascade.tofNSigmaOmKa()) > selCuts.NSigmaTOFKaon) { - // // omegaPassTOFSelection &= false; - // // } - // // } - // - // // if (particle == "xi") {return xiPassTOFSelection;} else {return omegaPassTOFSelection;} - // } + template + bool passesTOF(TCascade cascade, TString particle) + { + if (particle == "xi") + return cascade.tofXiCompatibility(selCuts.nSigmaTOFXi); + if (particle == "omega") + return cascade.tofOmegaCompatibility(selCuts.nSigmaTOFOmega); + return true; + } // checks whether gen cascade corresponds to PDG code template bool isValidPDG(TCascade cascade, TString particle) { - static constexpr int kPdgCodes[] = {3312, 3334}; // "XiMinus", "OmegaMinus" - if (particle == "xi" && std::abs(cascade.pdgCode()) == kPdgCodes[0]) + if (particle == "xi" && std::abs(cascade.pdgCode()) == PDG_t::kXiMinus) return true; - if (particle == "omega" && std::abs(cascade.pdgCode()) == kPdgCodes[1]) + if (particle == "omega" && std::abs(cascade.pdgCode()) == PDG_t::kOmegaMinus) return true; return false; } @@ -285,9 +281,9 @@ struct StrangeCascTrack { return false; int pdg = std::abs(cascmccore.pdgCode()); if (particle == "xi") - return (pdg == 3312); + return (pdg == PDG_t::kXiMinus); if (particle == "omega") - return (pdg == 3334); + return (pdg == PDG_t::kOmegaMinus); } return false; } @@ -337,110 +333,186 @@ struct StrangeCascTrack { } }(); - double mult = doProcessPP ? collision.centFT0M() : collision.centFT0C(); // ion collisions use FT0C for multiplicity, pp uses both + double mult = (doProcesspp || doProcesspO) ? collision.centFT0M() : collision.centFT0C(); // ion collisions use FT0C for multiplicity, pp uses both + + float efficiencyOmega = 1.0f; + float efficiencyXi = 1.0f; + float efficiencyOmegaErr = 0.0f; + float efficiencyXiErr = 0.0f; + float purityOmega = 1.0f; + float purityXi = 1.0f; + float purityOmegaErr = 0.0f; + float purityXiErr = 0.0f; - float efficiency = 1.0f; - float efficiencyError = 0.0f; - float purity = 1.0f; - float purityError = 0.0f; if (doApplyEfficiency) { - efficiency = hEfficiency->Interpolate(cascade.pt(), mult); - if (doPropagateEfficiency) { - efficiencyError = hEfficiencyUncertainty->Interpolate(cascade.pt(), mult); - } - if (efficiency == 0) { // check for zero efficiency, do not apply if the case - efficiency = 1.; - efficiencyError = 0.; + if constexpr (requires { cascade.topologyChi2(); }) { + efficiencyOmega = hEfficiencyOmegaTra->Interpolate(cascade.pt(), mult); + efficiencyXi = hEfficiencyXiTra->Interpolate(cascade.pt(), mult); + if (doPropagateEfficiency) { + efficiencyOmegaErr = hEfficiencyErrOmegaTra->Interpolate(cascade.pt(), mult); + efficiencyXiErr = hEfficiencyErrXiTra->Interpolate(cascade.pt(), mult); + } + if (efficiencyOmega == 0) { // check for zero efficiency, do not apply if the case + efficiencyOmega = 1.; + efficiencyOmegaErr = 0.; + } + if (efficiencyXi == 0) { // check for zero efficiency, do not apply if the case + efficiencyXi = 1.; + efficiencyXiErr = 0.; + } + } else { + efficiencyOmega = hEfficiencyOmegaStd->Interpolate(cascade.pt(), mult); + efficiencyXi = hEfficiencyXiStd->Interpolate(cascade.pt(), mult); + if (doPropagateEfficiency) { + efficiencyOmegaErr = hEfficiencyErrOmegaStd->Interpolate(cascade.pt(), mult); + efficiencyXiErr = hEfficiencyErrXiStd->Interpolate(cascade.pt(), mult); + } + if (efficiencyOmega == 0) { // check for zero efficiency, do not apply if the case + efficiencyOmega = 1.; + efficiencyOmegaErr = 0.; + } + if (efficiencyXi == 0) { // check for zero efficiency, do not apply if the case + efficiencyXi = 1.; + efficiencyXiErr = 0.; + } } } + if (doApplyPurity) { - purity = hPurity->Interpolate(cascade.pt(), mult); - if (doPropagatePurity) { - purityError = hPurityUncertainty->Interpolate(cascade.pt(), mult); - } - if (purity == 0) { // check for zero purity, do not apply if the case - purity = 1.; - purityError = 0.; + if constexpr (requires { cascade.topologyChi2(); }) { + purityOmega = hPurityOmegaTra->Interpolate(cascade.pt(), mult); + purityXi = hPurityXiTra->Interpolate(cascade.pt(), mult); + if (doPropagatePurity) { + purityOmegaErr = hPurityErrOmegaTra->Interpolate(cascade.pt(), mult); + purityXiErr = hPurityErrXiTra->Interpolate(cascade.pt(), mult); + } + if (purityOmega == 0) { // check for zero purity, do not apply if the case + purityOmega = 1.; + purityOmegaErr = 0.; + } + if (purityXi == 0) { // check for zero purity, do not apply if the case + purityXi = 1.; + purityXiErr = 0.; + } + } else { + purityOmega = hPurityOmegaStd->Interpolate(cascade.pt(), mult); + purityXi = hPurityXiStd->Interpolate(cascade.pt(), mult); + if (doPropagatePurity) { + purityOmegaErr = hPurityErrOmegaStd->Interpolate(cascade.pt(), mult); + purityXiErr = hPurityErrXiStd->Interpolate(cascade.pt(), mult); + } + if (purityOmega == 0) { // check for zero purity, do not apply if the case + purityOmega = 1.; + purityOmegaErr = 0.; + } + if (purityXi == 0) { // check for zero purity, do not apply if the case + purityXi = 1.; + purityXiErr = 0.; + } } } if (collision.index() != casccollid) { - histos.fill(HIST(kTypeNames[type]) + HIST("/EvMult"), mult); + histos.fill(HIST(kTypeNames[type]) + HIST("/NoSel/EvMult"), mult); casccollid = collision.index(); } double massXi = cascade.mXi(); double massOmega = cascade.mOmega(); double pt = cascade.pt(); - - histos.fill(HIST(kTypeNames[type]) + HIST("/DCAxy"), cascade.dcaXYCascToPV()); - histos.fill(HIST(kTypeNames[type]) + HIST("/DCAz"), cascade.dcaZCascToPV()); - histos.fill(HIST(kTypeNames[type]) + HIST("/Phi"), cascade.phi()); - histos.fill(HIST(kTypeNames[type]) + HIST("/Eta"), cascade.eta()); - histos.fill(HIST(kTypeNames[type]) + HIST("/MassXiNoSel"), massXi); - histos.fill(HIST(kTypeNames[type]) + HIST("/MassOmegaNoSel"), massOmega); + double v0cosPA = stdCasc.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); + + histos.fill(HIST(kTypeNames[type]) + HIST("/NoSel/DCAxy"), cascade.dcaXYCascToPV()); + histos.fill(HIST(kTypeNames[type]) + HIST("/NoSel/DCAz"), cascade.dcaZCascToPV()); + histos.fill(HIST(kTypeNames[type]) + HIST("/NoSel/Phi"), cascade.phi()); + histos.fill(HIST(kTypeNames[type]) + HIST("/NoSel/Eta"), cascade.eta()); + histos.fill(HIST(kTypeNames[type]) + HIST("/NoSel/BachCosPA"), stdCasc.bachBaryonCosPA()); + histos.fill(HIST(kTypeNames[type]) + HIST("/NoSel/V0CosPA"), v0cosPA); + histos.fill(HIST(kTypeNames[type]) + HIST("/NoSel/MassXi"), massXi); + histos.fill(HIST(kTypeNames[type]) + HIST("/NoSel/MassOmega"), massOmega); + + if constexpr (requires { collision.straMCCollisionId(); }) { + if (isMCTruth(stdCasc, "xi") || isMCTruth(stdCasc, "omega")) { + histos.fill(HIST(kTypeNames[type]) + HIST("/NoSel-Truth/DCAxy"), cascade.dcaXYCascToPV()); + histos.fill(HIST(kTypeNames[type]) + HIST("/NoSel-Truth/DCAz"), cascade.dcaZCascToPV()); + histos.fill(HIST(kTypeNames[type]) + HIST("/NoSel-Truth/DCAzVSpt"), pt, cascade.dcaZCascToPV()); + histos.fill(HIST(kTypeNames[type]) + HIST("/NoSel-Truth/Phi"), cascade.phi()); + histos.fill(HIST(kTypeNames[type]) + HIST("/NoSel-Truth/Eta"), cascade.eta()); + histos.fill(HIST(kTypeNames[type]) + HIST("/NoSel-Truth/EvMult"), mult); + histos.fill(HIST(kTypeNames[type]) + HIST("/NoSel-Truth/BachCosPA"), stdCasc.bachBaryonCosPA()); + histos.fill(HIST(kTypeNames[type]) + HIST("/NoSel-Truth/V0CosPA"), v0cosPA); + if (isMCTruth(stdCasc, "xi")) + histos.fill(HIST(kTypeNames[type]) + HIST("/NoSel-Truth/MassXi"), massXi); + if (isMCTruth(stdCasc, "omega")) + histos.fill(HIST(kTypeNames[type]) + HIST("/NoSel-Truth/MassOmega"), massOmega); + } + } // start checking selections bool passedAllSels = true; // apply general selection criteria if (doApplyCuts) { - if (isValidCasc(cascade)) { - histos.fill(HIST(kTypeNames[type]) + HIST("/MassXiGenSel"), massXi); - histos.fill(HIST(kTypeNames[type]) + HIST("/MassOmegaGenSel"), massOmega); - } else { + if (!isValidCasc(collision, stdCasc)) passedAllSels = false; - } } // apply tpc pid if (doApplyTPCPID) { - if (passesTPC(stdCasc)) { - histos.fill(HIST(kTypeNames[type]) + HIST("/MassXiTPCPID"), massXi); - histos.fill(HIST(kTypeNames[type]) + HIST("/MassOmegaTPCPID"), massOmega); - } else { + if (!passesTPC(stdCasc)) passedAllSels = false; - } } // apply tof pid bool passedAllSelsXi = passedAllSels; bool passedAllSelsOmega = passedAllSels; - // if (doApplyTOFPID) { - // if (passesTOF(cascade, "xi")) { - // histos.fill(HIST(kTypeNames[type]) + HIST("/MassXiTOFPID"), massXi); - // } else { - // passedAllSelsXi = false; - // } - // if (passesTOF(cascade, "omega")) { - // histos.fill(HIST(kTypeNames[type]) + HIST("/MassOmegaTOFPID"), massOmega); - // } else { - // passedAllSelsOmega = false; - // } - // } + if (doApplyTOFPID) { + if (!passesTOF(stdCasc, "xi")) + passedAllSelsXi = false; + if (!passesTOF(stdCasc, "omega")) + passedAllSelsOmega = false; + } // apply competing mass rej if (doCompetingMassRej) { - if (std::abs(massXi - pdgDB->Mass(3312)) > selCuts.compMassRej) { - histos.fill(HIST(kTypeNames[type]) + HIST("/MassOmegaMassSel"), massOmega); - } else { + if (!(std::abs(massXi - o2::constants::physics::MassXiMinus) > selCuts.compMassRej)) passedAllSelsOmega = false; - } } - // fill w/ cascs that passed all applied sels + // fill truth w/ cascs that passed all applied sels double binFillXi[3] = {massXi, pt, mult}; + + if constexpr (requires { collision.straMCCollisionId(); }) { + if (passedAllSels && (passedAllSelsXi || passedAllSelsOmega)) { // fill once for every desired cascade + if (isMCTruth(stdCasc, "xi") || isMCTruth(stdCasc, "omega")) { + histos.fill(HIST(kTypeNames[type]) + HIST("/Rec-Truth/DCAxy"), cascade.dcaXYCascToPV()); + histos.fill(HIST(kTypeNames[type]) + HIST("/Rec-Truth/DCAz"), cascade.dcaZCascToPV()); + histos.fill(HIST(kTypeNames[type]) + HIST("/Rec-Truth/DCAzVSpt"), pt, cascade.dcaZCascToPV()); + histos.fill(HIST(kTypeNames[type]) + HIST("/Rec-Truth/Phi"), cascade.phi()); + histos.fill(HIST(kTypeNames[type]) + HIST("/Rec-Truth/Eta"), cascade.eta()); + histos.fill(HIST(kTypeNames[type]) + HIST("/Rec-Truth/EvMult"), mult); + histos.fill(HIST(kTypeNames[type]) + HIST("/Rec-Truth/BachCosPA"), stdCasc.bachBaryonCosPA()); + histos.fill(HIST(kTypeNames[type]) + HIST("/Rec-Truth/V0CosPA"), v0cosPA); + } + } + } + + // fill rec if (passedAllSelsXi) { - histos.fill(HIST(kTypeNames[type]) + HIST("/MassXi"), massXi); - fillHist(histos.get(HIST(kTypeNames[type]) + HIST("/Xi")), binFillXi, efficiency, efficiencyError, purity, purityError); + histos.fill(HIST(kTypeNames[type]) + HIST("/Rec/MassXi"), massXi); + fillHist(histos.get(HIST(kTypeNames[type]) + HIST("/Rec/Xi")), binFillXi, efficiencyXi, efficiencyXiErr, purityXi, purityXiErr); if constexpr (requires { collision.straMCCollisionId(); }) { - if (doFillTruth && isMCTruth(stdCasc, "xi")) - histos.fill(HIST("MC/RecTruth/") + HIST(kTypeNames[type]) + HIST("/Xi"), massXi, pt, mult); + if (isMCTruth(stdCasc, "xi")) { + histos.fill(HIST(kTypeNames[type]) + HIST("/Rec-Truth/MassXi"), massXi); + histos.fill(HIST(kTypeNames[type]) + HIST("/Rec-Truth/Xi"), massXi, pt, mult); + } } } double binFillOmega[3] = {massOmega, pt, mult}; if (passedAllSelsOmega) { - histos.fill(HIST(kTypeNames[type]) + HIST("/MassOmega"), massOmega); - fillHist(histos.get(HIST(kTypeNames[type]) + HIST("/Omega")), binFillOmega, efficiency, efficiencyError, purity, purityError); + histos.fill(HIST(kTypeNames[type]) + HIST("/Rec/MassOmega"), massOmega); + fillHist(histos.get(HIST(kTypeNames[type]) + HIST("/Rec/Omega")), binFillOmega, efficiencyOmega, efficiencyOmegaErr, purityOmega, purityOmegaErr); if constexpr (requires { collision.straMCCollisionId(); }) { - if (doFillTruth && isMCTruth(stdCasc, "omega")) - histos.fill(HIST("MC/RecTruth/") + HIST(kTypeNames[type]) + HIST("/Omega"), massOmega, pt, mult); + if (isMCTruth(stdCasc, "omega")) { + histos.fill(HIST(kTypeNames[type]) + HIST("/Rec-Truth/MassOmega"), massOmega); + histos.fill(HIST(kTypeNames[type]) + HIST("/Rec-Truth/Omega"), massOmega, pt, mult); + } } } } @@ -456,44 +528,57 @@ struct StrangeCascTrack { histos.add("Events/Mult", "Multiplicity", kTH1F, {axesConfig.axisMult}); // for cascade processing static_for<0, 1>([&](auto type) { - histos.add(Form("%s/Phi", kTypeNames[type].data()), "Phi", kTH1F, {axesConfig.axisPhi}); - histos.add(Form("%s/Eta", kTypeNames[type].data()), "Eta", kTH1F, {axesConfig.axisEta}); - histos.add(Form("%s/DCAxy", kTypeNames[type].data()), "DCA to xy", kTH1F, {axesConfig.axisDCAxy}); - histos.add(Form("%s/DCAz", kTypeNames[type].data()), "DCA to z", kTH1F, {axesConfig.axisDCAz}); - histos.add(Form("%s/EvMult", kTypeNames[type].data()), "Multiplicity of events with >=1 cascade", kTH1F, {axesConfig.axisMult}); - // no selection applied - histos.add(Form("%s/MassOmegaNoSel", kTypeNames[type].data()), "Invariant mass hypothesis", kTH1F, {axesConfig.axisOmegaMass}); - histos.add(Form("%s/MassXiNoSel", kTypeNames[type].data()), "Invariant mass hypothesis", kTH1F, {axesConfig.axisXiMass}); - // only gen selection applied - histos.add(Form("%s/MassOmegaGenSel", kTypeNames[type].data()), "Invariant mass hypothesis", kTH1F, {axesConfig.axisOmegaMass}); - histos.add(Form("%s/MassXiGenSel", kTypeNames[type].data()), "Invariant mass hypothesis", kTH1F, {axesConfig.axisXiMass}); - // only tpc pid selection applied - histos.add(Form("%s/MassOmegaTPCPID", kTypeNames[type].data()), "Invariant mass hypothesis", kTH1F, {axesConfig.axisOmegaMass}); - histos.add(Form("%s/MassXiTPCPID", kTypeNames[type].data()), "Invariant mass hypothesis", kTH1F, {axesConfig.axisXiMass}); - // only tof pid selection applied - histos.add(Form("%s/MassOmegaTOFPID", kTypeNames[type].data()), "Invariant mass hypothesis", kTH1F, {axesConfig.axisOmegaMass}); - histos.add(Form("%s/MassXiTOFPID", kTypeNames[type].data()), "Invariant mass hypothesis", kTH1F, {axesConfig.axisXiMass}); - // only competing mass rejection selection applied - histos.add(Form("%s/MassOmegaMassSel", kTypeNames[type].data()), "Invariant mass hypothesis", kTH1F, {axesConfig.axisOmegaMass}); + // no selections applied + histos.add(Form("%s/NoSel/Phi", kTypeNames[type].data()), "Phi", kTH1F, {axesConfig.axisPhi}); + histos.add(Form("%s/NoSel/Eta", kTypeNames[type].data()), "Eta", kTH1F, {axesConfig.axisEta}); + histos.add(Form("%s/NoSel/DCAxy", kTypeNames[type].data()), "DCA to xy", kTH1F, {axesConfig.axisDCAxy}); + histos.add(Form("%s/NoSel/DCAz", kTypeNames[type].data()), "DCA to z", kTH1F, {axesConfig.axisDCAz}); + histos.add(Form("%s/NoSel/BachCosPA", kTypeNames[type].data()), "Bachelor cosPA", kTH1F, {{202, -1.1, 1.1}}); + histos.add(Form("%s/NoSel/V0CosPA", kTypeNames[type].data()), "V0 cosPA", kTH1F, {{202, -1.1, 1.1}}); + histos.add(Form("%s/NoSel/EvMult", kTypeNames[type].data()), "Multiplicity of events with >=1 cascade", kTH1F, {axesConfig.axisMult}); + histos.add(Form("%s/NoSel/MassXi", kTypeNames[type].data()), "Invariant mass hypothesis", kTH1F, {axesConfig.axisXiMass}); + histos.add(Form("%s/NoSel/MassOmega", kTypeNames[type].data()), "Invariant mass hypothesis", kTH1F, {axesConfig.axisOmegaMass}); + // mc truth for no selectrion + histos.add(Form("%s/NoSel-Truth/Phi", kTypeNames[type].data()), "Phi", kTH1F, {axesConfig.axisPhi}); + histos.add(Form("%s/NoSel-Truth/Eta", kTypeNames[type].data()), "Eta", kTH1F, {axesConfig.axisEta}); + histos.add(Form("%s/NoSel-Truth/DCAxy", kTypeNames[type].data()), "DCA to xy", kTH1F, {axesConfig.axisDCAxy}); + histos.add(Form("%s/NoSel-Truth/DCAz", kTypeNames[type].data()), "DCA to z", kTH1F, {axesConfig.axisDCAz}); + histos.add(Form("%s/NoSel-Truth/DCAzVSpt", kTypeNames[type].data()), "DCA to z vs pT", kTH2F, {axesConfig.axisPt, axesConfig.axisDCAz}); + histos.add(Form("%s/NoSel-Truth/BachCosPA", kTypeNames[type].data()), "Bachelor cosPA", kTH1F, {{202, -1.1, 1.1}}); + histos.add(Form("%s/NoSel-Truth/V0CosPA", kTypeNames[type].data()), "V0 cosPA", kTH1F, {{202, -1.1, 1.1}}); + histos.add(Form("%s/NoSel-Truth/EvMult", kTypeNames[type].data()), "Multiplicity of events with >=1 cascade", kTH1F, {axesConfig.axisMult}); + histos.add(Form("%s/NoSel-Truth/MassXi", kTypeNames[type].data()), "Invariant mass hypothesis", kTH1F, {axesConfig.axisXiMass}); + histos.add(Form("%s/NoSel-Truth/MassOmega", kTypeNames[type].data()), "Invariant mass hypothesis", kTH1F, {axesConfig.axisOmegaMass}); // passed all applied sels - histos.add(Form("%s/MassOmega", kTypeNames[type].data()), "Invariant mass hypothesis", kTH1F, {axesConfig.axisOmegaMass}); - histos.add(Form("%s/MassXi", kTypeNames[type].data()), "Invariant mass hypothesis", kTH1F, {axesConfig.axisXiMass}); - histos.add(Form("%s/Omega", kTypeNames[type].data()), "", kTHnD, {axesConfig.axisOmegaMass, axesConfig.axisPt, axesConfig.axisMult}); - histos.add(Form("%s/Xi", kTypeNames[type].data()), "", kTHnD, {axesConfig.axisXiMass, axesConfig.axisPt, axesConfig.axisMult}); + histos.add(Form("%s/Rec/MassOmega", kTypeNames[type].data()), "Invariant mass hypothesis", kTH1F, {axesConfig.axisOmegaMass}); + histos.add(Form("%s/Rec/MassXi", kTypeNames[type].data()), "Invariant mass hypothesis", kTH1F, {axesConfig.axisXiMass}); + histos.add(Form("%s/Rec/Omega", kTypeNames[type].data()), "", kTHnD, {axesConfig.axisOmegaMass, axesConfig.axisPt, axesConfig.axisMult}); + histos.add(Form("%s/Rec/Xi", kTypeNames[type].data()), "", kTHnD, {axesConfig.axisXiMass, axesConfig.axisPt, axesConfig.axisMult}); + // mc truth for all passed selections + histos.add(Form("%s/Rec-Truth/Phi", kTypeNames[type].data()), "Phi", kTH1F, {axesConfig.axisPhi}); + histos.add(Form("%s/Rec-Truth/Eta", kTypeNames[type].data()), "Eta", kTH1F, {axesConfig.axisEta}); + histos.add(Form("%s/Rec-Truth/DCAxy", kTypeNames[type].data()), "DCA to xy", kTH1F, {axesConfig.axisDCAxy}); + histos.add(Form("%s/Rec-Truth/DCAz", kTypeNames[type].data()), "DCA to z", kTH1F, {axesConfig.axisDCAz}); + histos.add(Form("%s/Rec-Truth/DCAzVSpt", kTypeNames[type].data()), "DCA to z vs pT", kTH2F, {axesConfig.axisPt, axesConfig.axisDCAz}); + histos.add(Form("%s/Rec-Truth/BachCosPA", kTypeNames[type].data()), "Bachelor cosPA", kTH1F, {{202, -1.1, 1.1}}); + histos.add(Form("%s/Rec-Truth/V0CosPA", kTypeNames[type].data()), "V0 cosPA", kTH1F, {{202, -1.1, 1.1}}); + histos.add(Form("%s/Rec-Truth/EvMult", kTypeNames[type].data()), "Multiplicity of events with >=1 cascade", kTH1F, {axesConfig.axisMult}); + histos.add(Form("%s/Rec-Truth/MassXi", kTypeNames[type].data()), "Invariant mass hypothesis", kTH1F, {axesConfig.axisXiMass}); + histos.add(Form("%s/Rec-Truth/MassOmega", kTypeNames[type].data()), "Invariant mass hypothesis", kTH1F, {axesConfig.axisOmegaMass}); + histos.add(Form("%s/Rec-Truth/Omega", kTypeNames[type].data()), "", kTHnD, {axesConfig.axisOmegaMass, axesConfig.axisPt, axesConfig.axisMult}); + histos.add(Form("%s/Rec-Truth/Xi", kTypeNames[type].data()), "", kTHnD, {axesConfig.axisXiMass, axesConfig.axisPt, axesConfig.axisMult}); }); // for MC-specific processing histos.add("MC/Gen/EvCounter", "Event Counter", kTH1F, {{1, 0, 1}}); - histos.add("MC/Gen/Xi", "Xi", kTH2F, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Xis - histos.add("MC/Gen/Omega", "Omega", kTH2F, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Omegas - histos.add("MC/Rec/EvCounter", "Event Counter", kTH1F, {{1, 0, 1}}); - histos.add("MC/Rec/EvMult", "Multiplicity", kTH1F, {axesConfig.axisMult}); - histos.add("MC/RecTruth/Standard/Omega", "", kTHnD, {axesConfig.axisOmegaMass, axesConfig.axisPt, axesConfig.axisMult}); - histos.add("MC/RecTruth/Standard/Xi", "", kTHnD, {axesConfig.axisXiMass, axesConfig.axisPt, axesConfig.axisMult}); - histos.add("MC/RecTruth/Tracked/Omega", "", kTHnD, {axesConfig.axisOmegaMass, axesConfig.axisPt, axesConfig.axisMult}); - histos.add("MC/RecTruth/Tracked/Xi", "", kTHnD, {axesConfig.axisXiMass, axesConfig.axisPt, axesConfig.axisMult}); + histos.add("MC/Gen/Xi", "Xi", kTH2F, {axesConfig.axisPt, axesConfig.axisMult}); // generated Xis + histos.add("MC/Gen/Omega", "Omega", kTH2F, {axesConfig.axisPt, axesConfig.axisMult}); // generated Omegas + histos.add("MC/Gen/PrimaryXi", "Xi primaries", kTH2F, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Xis + histos.add("MC/Gen/PrimaryOmega", "Omega primaries", kTH2F, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Omegas + histos.add("MC/Rec/EvCounter", "Event Counter", kTH1F, {{1, 0, 1}}); // counter of all recreated events + histos.add("MC/Rec/EvMult", "Multiplicity", kTH1F, {axesConfig.axisMult}); // multiplicity of all recreated events } - void processDerivedData(DerCollisionWMult const& collision, DerCascDatas const& allCascs, DerTraCascDatas const& traCascs, DauTracks const&) + void processDerivedData(DerCollisionWMults::iterator const& collision, DerCascDatas const& allCascs, DerTraCascDatas const& traCascs, DauTracks const&) { fillEvents(collision); // save info about all processed events if (doApplyEfficiency) { @@ -511,25 +596,29 @@ struct StrangeCascTrack { auto slicedRecColls = recColls.sliceBy(perMcCollision, genColl.globalIndex()); for (auto const& recColl : slicedRecColls) { histos.fill(HIST("MC/Rec/EvCounter"), 0.5); - double casc_mult = doProcessPP ? recColl.centFT0M() : recColl.centFT0C(); + double casc_mult = (doProcesspp || doProcesspO) ? recColl.centFT0M() : recColl.centFT0C(); histos.fill(HIST("MC/Rec/EvMult"), casc_mult); int64_t genCollId = recColl.straMCCollisionId(); for (auto const& casc : genCascs) { if (casc.straMCCollisionId() != genCollId) continue; // safety check - if (!casc.isPhysicalPrimary()) - continue; // skip non-primary particles double casc_pt = std::sqrt(std::pow(casc.pxMC(), 2) + std::pow(casc.pyMC(), 2)); if (isValidPDG(casc, "xi")) histos.fill(HIST("MC/Gen/Xi"), casc_pt, casc_mult); if (isValidPDG(casc, "omega")) histos.fill(HIST("MC/Gen/Omega"), casc_pt, casc_mult); + if (casc.isPhysicalPrimary()) { + if (isValidPDG(casc, "xi")) + histos.fill(HIST("MC/Gen/PrimaryXi"), casc_pt, casc_mult); + if (isValidPDG(casc, "omega")) + histos.fill(HIST("MC/Gen/PrimaryOmega"), casc_pt, casc_mult); + } } } } } - void processDerivedMCRec(DerMCRecCollision const& collision, DerMCRecCascDatas const& allCascs, DerMCRecTraCascDatas const& traCascs, DauTracks const&, DerMCGenCascades const&) + void processDerivedMCRec(DerMCRecCollisions::iterator const& collision, DerMCRecCascDatas const& allCascs, DerMCRecTraCascDatas const& traCascs, DauTracks const&, DerMCGenCascades const&) { fillEvents(collision); // save info about all processed events if (doApplyEfficiency) { @@ -539,14 +628,14 @@ struct StrangeCascTrack { analyseCascs(collision, traCascs); // process tracked cascades } - PROCESS_SWITCH(StrangeCascTrack, processDerivedData, "process derived data", true); - PROCESS_SWITCH(StrangeCascTrack, processDerivedMCGen, "process derived generated mc data", false); - PROCESS_SWITCH(StrangeCascTrack, processDerivedMCRec, "process derived reconstructed mc data", false); // mc and data are mutually exclusive! + PROCESS_SWITCH(strangecasctrack, processDerivedData, "process derived data", true); + PROCESS_SWITCH(strangecasctrack, processDerivedMCGen, "process derived generated mc data", false); + PROCESS_SWITCH(strangecasctrack, processDerivedMCRec, "process derived reconstructed mc data", false); // mc and data are mutually exclusive! }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), }; }