diff --git a/PWGLF/Tasks/Strangeness/strangeCascTrack.cxx b/PWGLF/Tasks/Strangeness/strangeCascTrack.cxx index 6bfff6476de..def45a3b637 100644 --- a/PWGLF/Tasks/Strangeness/strangeCascTrack.cxx +++ b/PWGLF/Tasks/Strangeness/strangeCascTrack.cxx @@ -60,19 +60,17 @@ using DauTracks = soa::Join; struct StrangeCascTrack { - Service ccdb; + //* Service ccdb; Service pdgDB; PresliceUnsorted> perMcCollision = aod::v0data::straMCCollisionId; + PresliceUnsorted cascsPerMcCollision = aod::cascdata::straMCCollisionId; HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; // subprocess switches: - 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"}; + Configurable doProcessIons{"doProcessIons", false, "true for PbPb and OO, false for pp and pO"}; Configurable doApplyEventCuts{"doApplyEventCuts", true, "apply general event cuts"}; // event filter - PVz, sel8, INEL>0 // Xi selections @@ -86,22 +84,33 @@ struct StrangeCascTrack { Configurable doApplyTPCPIDOmega{"doApplyTPCPIDOmega", true, "apply tpc pid to dau tracks (Omega)"}; Configurable doApplyTOFPIDOmega{"doApplyTOFPIDOmega", true, "apply tof pid to dau tracks (Omega)"}; Configurable doCompetingMassRej{"doCompetingMassRej", true, "competing mass rejection (Omega)"}; - // efficiency and purity corrections on the fly (warning: to be avoided because interpolation causes errors): - // only correct by pt - Configurable doApplyEfficiency1D{"doApplyEfficiency1D", false, "apply efficiency correction"}; - Configurable doPropagateEfficiency1D{"doPropagateEfficiency1D", false, "apply efficiency propagation"}; - Configurable doApplyPurity1D{"doApplyPurity1D", false, "apply purity correction"}; - Configurable doPropagatePurity1D{"doPropagatePurity1D", false, "apply purity propagation"}; - // correct by both pt and mult - Configurable doApplyEfficiency2D{"doApplyEfficiency2D", false, "apply efficiency correction"}; - Configurable doPropagateEfficiency2D{"doPropagateEfficiency2D", false, "apply efficiency propagation"}; - Configurable doApplyPurity2D{"doApplyPurity2D", false, "apply purity correction"}; - Configurable doPropagatePurity2D{"doPropagatePurity2D", false, "apply purity propagation"}; - Configurable ccdburl{"ccdburl", "http://alice-ccdb.cern.ch", "url of the ccdb repository to use"}; - Configurable efficiencyCCDBPathpp{"efficiencyCCDBPathpp", "Users/y/yparovia/LHC24f4d", "Path of the efficiency and purity corrections (pp)"}; - Configurable efficiencyCCDBPathPbPb{"efficiencyCCDBPathPbPb", "Users/y/yparovia/LHC25f3", "Path of the efficiency and purity corrections (PbPb)"}; - Configurable efficiencyCCDBPathOO{"efficiencyCCDBPathOO", "Users/y/yparovia/LHC25h3", "Path of the efficiency and purity corrections (OO)"}; - Configurable efficiencyCCDBPathpO{"efficiencyCCDBPathpO", "Users/y/yparovia/LHC25h2", "Path of the efficiency and purity corrections (pO)"}; + //* // efficiency and purity corrections on the fly (warning: to be avoided because interpolation causes errors): + //* // only correct by pt + //* Configurable doApplyEfficiency1D{"doApplyEfficiency1D", false, "apply efficiency correction"}; + //* Configurable doPropagateEfficiency1D{"doPropagateEfficiency1D", false, "apply efficiency propagation"}; + //* Configurable doApplyPurity1D{"doApplyPurity1D", false, "apply purity correction"}; + //* Configurable doPropagatePurity1D{"doPropagatePurity1D", false, "apply purity propagation"}; + //* // correct by both pt and mult + //* Configurable doApplyEfficiency2D{"doApplyEfficiency2D", false, "apply efficiency correction"}; + //* Configurable doPropagateEfficiency2D{"doPropagateEfficiency2D", false, "apply efficiency propagation"}; + //* Configurable doApplyPurity2D{"doApplyPurity2D", false, "apply purity correction"}; + //* Configurable doPropagatePurity2D{"doPropagatePurity2D", false, "apply purity propagation"}; + //* Configurable ccdburl{"ccdburl", "http://alice-ccdb.cern.ch", "url of the ccdb repository to use"}; + //* Configurable efficiencyCCDBPathpp{"efficiencyCCDBPathpp", "Users/y/yparovia/LHC24f4d", "Path of the efficiency and purity corrections (pp)"}; + //* Configurable efficiencyCCDBPathPbPb{"efficiencyCCDBPathPbPb", "Users/y/yparovia/LHC25f3", "Path of the efficiency and purity corrections (PbPb)"}; + //* Configurable efficiencyCCDBPathOO{"efficiencyCCDBPathOO", "Users/y/yparovia/LHC25h3", "Path of the efficiency and purity corrections (OO)"}; + //* Configurable efficiencyCCDBPathpO{"efficiencyCCDBPathpO", "Users/y/yparovia/LHC25h2", "Path of the efficiency and purity corrections (pO)"}; + + // axes + struct : ConfigurableGroup { + 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, 2.5, 3.0, 3.5, 4.0, 5.0, 6.0, 10.0}, "p_{T} (GeV/c)"}; + 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; // event and dau track selection struct : ConfigurableGroup { @@ -109,6 +118,8 @@ struct StrangeCascTrack { Configurable cutDoINEL{"cutDoINEL", true, "choose events with INEL>0"}; Configurable cutZVertex{"cutZVertex", 10.0f, "max Z-vertex position"}; Configurable cutDoSel8{"cutDoSel8", true, "choose events with sel8"}; + Configurable cutDoNoPileup{"cutDoNoPileup", true, "choose events with no bunch pileup"}; + Configurable cutDoGoodFT0{"cutDoGoodFT0", true, "choose events with IsGoodZvtxFT0VsPV"}; // cascade cuts Configurable cutDoPropagateDCA{"cutDoPropagateDCA", false, "choose events with sel8"}; Configurable cutPropDCAtoPVxy{"cutPropDCAtoPVxy", 0.02f, "max cascade dca to PV in xy - propagated"}; @@ -118,7 +129,10 @@ struct StrangeCascTrack { Configurable cutMaxV0CosPA{"cutMaxV0CosPA", 1.1f, "max V0 cosPA"}; Configurable cutMinBachCosPA{"cutMinBachCosPA", -1.1f, "min Bachelor cosPA"}; Configurable cutMaxBachCosPA{"cutMaxBachCosPA", 1.1f, "max Bachelor cosPA"}; - Configurable cutMinCascCosPA{"cutMinCascCosPA", 0.995f, "min cascade cosPA"}; + Configurable> cutMinCascCosPaVsPt{ + "cutMinCascCosPaVsPt", + {0.993, 0.994, 0.995, 0.996, 0.997, 0.997, 0.998, 0.998, 0.999, 0.999, 0.999}, + "Min Casc CosPA per pT bin (same binning as axisPt)"}; Configurable cutRapidity{"cutRapidity", 0.5f, "max rapidity"}; Configurable cutDauEta{"cutDauEta", 1.0f, "max eta of dau tracks"}; Configurable cutCompMassRej{"cutCompMassRej", 0.008f, "Competing mass rejection"}; @@ -140,141 +154,131 @@ struct StrangeCascTrack { Configurable cutNSigmaTOFOmega{"cutNSigmaTOFOmega", 3, "cutNSigmaTOFOmega"}; } selCuts; - // axes - struct : ConfigurableGroup { - 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, 2.5, 3.0, 3.5, 4.0, 5.0, 6.0, 10.0}, "p_{T} (GeV/c)"}; - 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; - // cascade reconstruction Types static constexpr std::string_view TypeNames[] = {"Standard", "Tracked"}; - // for efficiency and purity corrections - TH1F* hEfficiencyOmegaStd1D; - TH1F* hEfficiencyOmegaTra1D; - TH1F* hEfficiencyXiStd1D; - TH1F* hEfficiencyXiTra1D; - - TH1F* hEfficiencyErrOmegaStd1D; - TH1F* hEfficiencyErrOmegaTra1D; - TH1F* hEfficiencyErrXiStd1D; - TH1F* hEfficiencyErrXiTra1D; - - TH1F* hPurityOmegaStd1D; - TH1F* hPurityOmegaTra1D; - TH1F* hPurityXiStd1D; - TH1F* hPurityXiTra1D; - - TH1F* hPurityErrOmegaStd1D; - TH1F* hPurityErrOmegaTra1D; - TH1F* hPurityErrXiStd1D; - TH1F* hPurityErrXiTra1D; - - TH2F* hEfficiencyOmegaStd2D; - TH2F* hEfficiencyOmegaTra2D; - TH2F* hEfficiencyXiStd2D; - TH2F* hEfficiencyXiTra2D; - - TH2F* hEfficiencyErrOmegaStd2D; - TH2F* hEfficiencyErrOmegaTra2D; - TH2F* hEfficiencyErrXiStd2D; - TH2F* hEfficiencyErrXiTra2D; - - TH2F* hPurityOmegaStd2D; - TH2F* hPurityOmegaTra2D; - TH2F* hPurityXiStd2D; - TH2F* hPurityXiTra2D; - - TH2F* hPurityErrOmegaStd2D; - TH2F* hPurityErrOmegaTra2D; - TH2F* hPurityErrXiStd2D; - TH2F* hPurityErrXiTra2D; - - int mRunNumber; - // loads efficiencies and purities - void initEfficiencyFromCCDB(int64_t runNumber, int64_t timestamp) - { - if (mRunNumber == runNumber) { - return; - } - mRunNumber = runNumber; - LOG(info) << "Loading efficiencies and purities from CCDB for run " << mRunNumber << " now..."; - auto timeStamp = timestamp; - - std::string efficiencyCCDBPath = [&]() { - if (doProcesspp) { - return efficiencyCCDBPathpp; - } else if (doProcesspO) { - return efficiencyCCDBPathpO; - } else if (doProcessPbPb) { - return efficiencyCCDBPathPbPb; - } - return efficiencyCCDBPathOO; - }(); - - TList* listEfficiencies = ccdb->getForTimeStamp(efficiencyCCDBPath, timeStamp); - - if (!listEfficiencies) { - LOG(fatal) << "Problem getting TList object with efficiencies and purities!"; - } - - hEfficiencyOmegaStd1D = static_cast(listEfficiencies->FindObject("Eff_Omega_Standard_byPt")); - hEfficiencyOmegaTra1D = static_cast(listEfficiencies->FindObject("Eff_Omega_Tracked_byPt")); - hEfficiencyXiStd1D = static_cast(listEfficiencies->FindObject("Eff_Xi_Standard_byPt")); - hEfficiencyXiTra1D = static_cast(listEfficiencies->FindObject("Eff_Xi_Tracked_byPt")); - hEfficiencyErrOmegaStd1D = static_cast(listEfficiencies->FindObject("EffErr_Omega_Standard_byPt")); - hEfficiencyErrOmegaTra1D = static_cast(listEfficiencies->FindObject("EffErr_Omega_Tracked_byPt")); - hEfficiencyErrXiStd1D = static_cast(listEfficiencies->FindObject("EffErr_Xi_Standard_byPt")); - hEfficiencyErrXiTra1D = static_cast(listEfficiencies->FindObject("EffErr_Xi_Tracked_byPt")); - hPurityOmegaStd1D = static_cast(listEfficiencies->FindObject("Pur_Omega_Standard_byPt")); - hPurityOmegaTra1D = static_cast(listEfficiencies->FindObject("Pur_Omega_Tracked_byPt")); - hPurityXiStd1D = static_cast(listEfficiencies->FindObject("Pur_Xi_Standard_byPt")); - hPurityXiTra1D = static_cast(listEfficiencies->FindObject("Pur_Xi_Tracked_byPt")); - hPurityErrOmegaStd1D = static_cast(listEfficiencies->FindObject("PurErr_Omega_Standard_byPt")); - hPurityErrOmegaTra1D = static_cast(listEfficiencies->FindObject("PurErr_Omega_Tracked_byPt")); - hPurityErrXiStd1D = static_cast(listEfficiencies->FindObject("PurErr_Xi_Standard_byPt")); - hPurityErrXiTra1D = static_cast(listEfficiencies->FindObject("PurErr_Xi_Tracked_byPt")); - - hEfficiencyOmegaStd2D = static_cast(listEfficiencies->FindObject("Eff_Omega_Standard_byPtMult")); - hEfficiencyOmegaTra2D = static_cast(listEfficiencies->FindObject("Eff_Omega_Tracked_byPtMult")); - hEfficiencyXiStd2D = static_cast(listEfficiencies->FindObject("Eff_Xi_Standard_byPtMult")); - hEfficiencyXiTra2D = static_cast(listEfficiencies->FindObject("Eff_Xi_Tracked_byPtMult")); - hEfficiencyErrOmegaStd2D = static_cast(listEfficiencies->FindObject("EffErr_Omega_Standard_byPtMult")); - hEfficiencyErrOmegaTra2D = static_cast(listEfficiencies->FindObject("EffErr_Omega_Tracked_byPtMult")); - hEfficiencyErrXiStd2D = static_cast(listEfficiencies->FindObject("EffErr_Xi_Standard_byPtMult")); - hEfficiencyErrXiTra2D = static_cast(listEfficiencies->FindObject("EffErr_Xi_Tracked_byPtMult")); - hPurityOmegaStd2D = static_cast(listEfficiencies->FindObject("Pur_Omega_Standard_byPtMult")); - hPurityOmegaTra2D = static_cast(listEfficiencies->FindObject("Pur_Omega_Tracked_byPtMult")); - hPurityXiStd2D = static_cast(listEfficiencies->FindObject("Pur_Xi_Standard_byPtMult")); - hPurityXiTra2D = static_cast(listEfficiencies->FindObject("Pur_Xi_Tracked_byPtMult")); - hPurityErrOmegaStd2D = static_cast(listEfficiencies->FindObject("PurErr_Omega_Standard_byPtMult")); - hPurityErrOmegaTra2D = static_cast(listEfficiencies->FindObject("PurErr_Omega_Tracked_byPtMult")); - hPurityErrXiStd2D = static_cast(listEfficiencies->FindObject("PurErr_Xi_Standard_byPtMult")); - hPurityErrXiTra2D = static_cast(listEfficiencies->FindObject("PurErr_Xi_Tracked_byPtMult")); + //* // for on-the-fly efficiency and purity corrections - working but deprecated + //* TH1F* hEfficiencyOmegaStd1D; + //* TH1F* hEfficiencyOmegaTra1D; + //* TH1F* hEfficiencyXiStd1D; + //* TH1F* hEfficiencyXiTra1D; + //* + //* TH1F* hEfficiencyErrOmegaStd1D; + //* TH1F* hEfficiencyErrOmegaTra1D; + //* TH1F* hEfficiencyErrXiStd1D; + //* TH1F* hEfficiencyErrXiTra1D; + //* + //* TH1F* hPurityOmegaStd1D; + //* TH1F* hPurityOmegaTra1D; + //* TH1F* hPurityXiStd1D; + //* TH1F* hPurityXiTra1D; + //* + //* TH1F* hPurityErrOmegaStd1D; + //* TH1F* hPurityErrOmegaTra1D; + //* TH1F* hPurityErrXiStd1D; + //* TH1F* hPurityErrXiTra1D; + //* + //* TH2F* hEfficiencyOmegaStd2D; + //* TH2F* hEfficiencyOmegaTra2D; + //* TH2F* hEfficiencyXiStd2D; + //* TH2F* hEfficiencyXiTra2D; + //* + //* TH2F* hEfficiencyErrOmegaStd2D; + //* TH2F* hEfficiencyErrOmegaTra2D; + //* TH2F* hEfficiencyErrXiStd2D; + //* TH2F* hEfficiencyErrXiTra2D; + //* + //* TH2F* hPurityOmegaStd2D; + //* TH2F* hPurityOmegaTra2D; + //* TH2F* hPurityXiStd2D; + //* TH2F* hPurityXiTra2D; + //* + //* TH2F* hPurityErrOmegaStd2D; + //* TH2F* hPurityErrOmegaTra2D; + //* TH2F* hPurityErrXiStd2D; + //* TH2F* hPurityErrXiTra2D; + //* + //* int mRunNumber; + //* // loads efficiencies and purities + //* void initEfficiencyFromCCDB(int64_t runNumber, int64_t timestamp) + //* { + //* if (mRunNumber == runNumber) { + //* return; + //* } + //* mRunNumber = runNumber; + //* LOG(info) << "Loading efficiencies and purities from CCDB for run " << mRunNumber << " now..."; + //* auto timeStamp = timestamp; + //* + //* std::string efficiencyCCDBPath = [&]() { + //* if (doProcesspp) { + //* return efficiencyCCDBPathpp; + //* } else if (doProcesspO) { + //* return efficiencyCCDBPathpO; + //* } else if (doProcessPbPb) { + //* return efficiencyCCDBPathPbPb; + //* } + //* return efficiencyCCDBPathOO; + //* }(); + //* + //* TList* listEfficiencies = ccdb->getForTimeStamp(efficiencyCCDBPath, timeStamp); + //* + //* if (!listEfficiencies) { + //* LOG(fatal) << "Problem getting TList object with efficiencies and purities!"; + //* } + //* + //* hEfficiencyOmegaStd1D = static_cast(listEfficiencies->FindObject("Eff_Omega_Standard_byPt")); + //* hEfficiencyOmegaTra1D = static_cast(listEfficiencies->FindObject("Eff_Omega_Tracked_byPt")); + //* hEfficiencyXiStd1D = static_cast(listEfficiencies->FindObject("Eff_Xi_Standard_byPt")); + //* hEfficiencyXiTra1D = static_cast(listEfficiencies->FindObject("Eff_Xi_Tracked_byPt")); + //* hEfficiencyErrOmegaStd1D = static_cast(listEfficiencies->FindObject("EffErr_Omega_Standard_byPt")); + //* hEfficiencyErrOmegaTra1D = static_cast(listEfficiencies->FindObject("EffErr_Omega_Tracked_byPt")); + //* hEfficiencyErrXiStd1D = static_cast(listEfficiencies->FindObject("EffErr_Xi_Standard_byPt")); + //* hEfficiencyErrXiTra1D = static_cast(listEfficiencies->FindObject("EffErr_Xi_Tracked_byPt")); + //* hPurityOmegaStd1D = static_cast(listEfficiencies->FindObject("Pur_Omega_Standard_byPt")); + //* hPurityOmegaTra1D = static_cast(listEfficiencies->FindObject("Pur_Omega_Tracked_byPt")); + //* hPurityXiStd1D = static_cast(listEfficiencies->FindObject("Pur_Xi_Standard_byPt")); + //* hPurityXiTra1D = static_cast(listEfficiencies->FindObject("Pur_Xi_Tracked_byPt")); + //* hPurityErrOmegaStd1D = static_cast(listEfficiencies->FindObject("PurErr_Omega_Standard_byPt")); + //* hPurityErrOmegaTra1D = static_cast(listEfficiencies->FindObject("PurErr_Omega_Tracked_byPt")); + //* hPurityErrXiStd1D = static_cast(listEfficiencies->FindObject("PurErr_Xi_Standard_byPt")); + //* hPurityErrXiTra1D = static_cast(listEfficiencies->FindObject("PurErr_Xi_Tracked_byPt")); + //* + //* hEfficiencyOmegaStd2D = static_cast(listEfficiencies->FindObject("Eff_Omega_Standard_byPtMult")); + //* hEfficiencyOmegaTra2D = static_cast(listEfficiencies->FindObject("Eff_Omega_Tracked_byPtMult")); + //* hEfficiencyXiStd2D = static_cast(listEfficiencies->FindObject("Eff_Xi_Standard_byPtMult")); + //* hEfficiencyXiTra2D = static_cast(listEfficiencies->FindObject("Eff_Xi_Tracked_byPtMult")); + //* hEfficiencyErrOmegaStd2D = static_cast(listEfficiencies->FindObject("EffErr_Omega_Standard_byPtMult")); + //* hEfficiencyErrOmegaTra2D = static_cast(listEfficiencies->FindObject("EffErr_Omega_Tracked_byPtMult")); + //* hEfficiencyErrXiStd2D = static_cast(listEfficiencies->FindObject("EffErr_Xi_Standard_byPtMult")); + //* hEfficiencyErrXiTra2D = static_cast(listEfficiencies->FindObject("EffErr_Xi_Tracked_byPtMult")); + //* hPurityOmegaStd2D = static_cast(listEfficiencies->FindObject("Pur_Omega_Standard_byPtMult")); + //* hPurityOmegaTra2D = static_cast(listEfficiencies->FindObject("Pur_Omega_Tracked_byPtMult")); + //* hPurityXiStd2D = static_cast(listEfficiencies->FindObject("Pur_Xi_Standard_byPtMult")); + //* hPurityXiTra2D = static_cast(listEfficiencies->FindObject("Pur_Xi_Tracked_byPtMult")); + //* hPurityErrOmegaStd2D = static_cast(listEfficiencies->FindObject("PurErr_Omega_Standard_byPtMult")); + //* hPurityErrOmegaTra2D = static_cast(listEfficiencies->FindObject("PurErr_Omega_Tracked_byPtMult")); + //* hPurityErrXiStd2D = static_cast(listEfficiencies->FindObject("PurErr_Xi_Standard_byPtMult")); + //* hPurityErrXiTra2D = static_cast(listEfficiencies->FindObject("PurErr_Xi_Tracked_byPtMult")); + //* + //* if (doPropagateEfficiency1D && (!hEfficiencyErrOmegaStd1D || !hEfficiencyErrOmegaTra1D || !hEfficiencyErrXiStd1D || !hEfficiencyErrXiTra1D)) + //* LOG(fatal) << "Problem getting hEfficiencyUncertainty!"; + //* if (doPropagatePurity1D && (!hPurityErrOmegaStd1D || !hPurityErrOmegaTra1D || !hPurityErrXiStd1D || !hPurityErrXiTra1D)) + //* LOG(fatal) << "Problem getting hPurityUncertainty!"; + //* LOG(info) << "Efficiencies and purities now loaded for " << mRunNumber; + //* + //* if (doPropagateEfficiency2D && (!hEfficiencyErrOmegaStd2D || !hEfficiencyErrOmegaTra2D || !hEfficiencyErrXiStd2D || !hEfficiencyErrXiTra2D)) + //* LOG(fatal) << "Problem getting hEfficiencyUncertainty!"; + //* if (doPropagatePurity2D && (!hPurityErrOmegaStd2D || !hPurityErrOmegaTra2D || !hPurityErrXiStd2D || !hPurityErrXiTra2D)) + //* LOG(fatal) << "Problem getting hPurityUncertainty!"; + //* LOG(info) << "Efficiencies and purities now loaded for " << mRunNumber; + //* } - if (doPropagateEfficiency1D && (!hEfficiencyErrOmegaStd1D || !hEfficiencyErrOmegaTra1D || !hEfficiencyErrXiStd1D || !hEfficiencyErrXiTra1D)) - LOG(fatal) << "Problem getting hEfficiencyUncertainty!"; - if (doPropagatePurity1D && (!hPurityErrOmegaStd1D || !hPurityErrOmegaTra1D || !hPurityErrXiStd1D || !hPurityErrXiTra1D)) - LOG(fatal) << "Problem getting hPurityUncertainty!"; - LOG(info) << "Efficiencies and purities now loaded for " << mRunNumber; - - if (doPropagateEfficiency2D && (!hEfficiencyErrOmegaStd2D || !hEfficiencyErrOmegaTra2D || !hEfficiencyErrXiStd2D || !hEfficiencyErrXiTra2D)) - LOG(fatal) << "Problem getting hEfficiencyUncertainty!"; - if (doPropagatePurity2D && (!hPurityErrOmegaStd2D || !hPurityErrOmegaTra2D || !hPurityErrXiStd2D || !hPurityErrXiTra2D)) - LOG(fatal) << "Problem getting hPurityUncertainty!"; - LOG(info) << "Efficiencies and purities now loaded for " << mRunNumber; - } // general info about processed events template void fillEvents(TEvent const& collision) { histos.fill(HIST("NoSel-Events/EvCounter"), 0.5); - double mult = (doProcesspp || doProcesspO) ? collision.centFT0M() : collision.centFT0C(); + double mult = (doProcessIons) ? collision.centFT0C() : collision.centFT0M(); histos.fill(HIST("NoSel-Events/Mult"), mult); double pvx = collision.posX(); double pvy = collision.posY(); @@ -287,20 +291,34 @@ struct StrangeCascTrack { bool isValidEvent(TEvent collision) { bool passedAllSels = true; + //* inel>0 cut if (!selCuts.cutDoINEL || collision.multNTracksPVeta1() > 0) histos.fill(HIST("Rec-Events/EvFilter"), 0.5); else passedAllSels = false; + //* pvz cut if (std::abs(collision.posZ()) < selCuts.cutZVertex) histos.fill(HIST("Rec-Events/EvFilter"), 1.5); else passedAllSels = false; + //* sel8 cut if (!selCuts.cutDoSel8 || collision.sel8()) histos.fill(HIST("Rec-Events/EvFilter"), 2.5); else passedAllSels = false; - if (passedAllSels) + //* pileup cut + if (!selCuts.cutDoNoPileup || collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) histos.fill(HIST("Rec-Events/EvFilter"), 3.5); + else + passedAllSels = false; + //* good ft0 z-vertex vs pv cut + if (!selCuts.cutDoGoodFT0 || collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) + histos.fill(HIST("Rec-Events/EvFilter"), 4.5); + else + passedAllSels = false; + //* all cuts + if (passedAllSels) + histos.fill(HIST("Rec-Events/EvFilter"), 5.5); return passedAllSels; } // checks cascade pt @@ -328,7 +346,7 @@ struct StrangeCascTrack { } if (cascade.pt() < ptMin || cascade.pt() > ptMax) passedSel = false; - // histos.fill(HIST(Form("%s/Rec/Filters%s", TypeNames[Type].data(), particle)), 0.5); + //* histos.fill(HIST(Form("%s/Rec/Filters%s", TypeNames[Type].data(), particle)), 0.5); return passedSel; } // checks general selection criteria for cascades @@ -385,7 +403,12 @@ struct StrangeCascTrack { // Cascade cosPA bool passedCascCosPA = true; double casccospa = cascade.casccosPA(collision.posX(), collision.posY(), collision.posZ()); - if (casccospa < selCuts.cutMinCascCosPA) { + const auto& edges = axesConfig.axisPt.value; + int ptBin = std::upper_bound(edges.begin(), edges.end(), cascade.pt()) - edges.begin() - 1; + if (ptBin < 0 || ptBin >= static_cast(selCuts.cutMinCascCosPaVsPt->size())) { + ptBin = 1; + } // safety check - if pt bin not determined, default to loosest cut + if (casccospa < selCuts.cutMinCascCosPaVsPt->at(ptBin)) { passedCascCosPA = false; passedAllSels = false; } @@ -485,17 +508,19 @@ struct StrangeCascTrack { hist->SetBinContent(bin, currentContent); hist->SetBinError2(bin, currentError2); } - // calculates DCA from cosPA - template - double calculateDCA(TEvent collision, double cosPA, double decX, double decY, double decZ) - { - double pvX = collision.posX(); - double pvY = collision.posX(); - double pvZ = collision.posX(); - double sinPA = std::sqrt(1 - cosPA * cosPA); - double dca = sinPA * std::sqrt(std::pow(decX - pvX, 2) + std::pow(decY - pvY, 2) + std::pow(decZ - pvZ, 2)); - return dca; - } + + ////* calculates DCA from cosPA - deprecated due to low precision + ////* template + ////* double calculateDCA(TEvent collision, double cosPA, double decX, double decY, double decZ) + ////* { + ////* double pvX = collision.posX(); + ////* double pvY = collision.posX(); + ////* double pvZ = collision.posX(); + ////* double sinPA = std::sqrt(1 - cosPA * cosPA); + ////* double dca = sinPA * std::sqrt(std::pow(decX - pvX, 2) + std::pow(decY - pvY, 2) + std::pow(decZ - pvZ, 2)); + ////* return dca; + ////* } + // applies selections for and fills histograms template void analyseCascs(TEvent collision, TCascs cascades) @@ -539,143 +564,143 @@ struct StrangeCascTrack { float purityOmegaErr = 0.0f; float purityXiErr = 0.0f; - double mult = (doProcesspp || doProcesspO) ? collision.centFT0M() : collision.centFT0C(); // ion collisions use FT0C for multiplicity, pp uses both - - if (doApplyEfficiency1D) { - if constexpr (requires { cascade.topologyChi2(); }) { - efficiencyOmega = hEfficiencyOmegaTra1D->Interpolate(cascade.pt()); - efficiencyXi = hEfficiencyXiTra1D->Interpolate(cascade.pt()); - if (doPropagateEfficiency1D) { - efficiencyOmegaErr = hEfficiencyErrOmegaTra1D->Interpolate(cascade.pt()); - efficiencyXiErr = hEfficiencyErrXiTra1D->Interpolate(cascade.pt()); - } - 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 = hEfficiencyOmegaStd1D->Interpolate(cascade.pt()); - efficiencyXi = hEfficiencyXiStd1D->Interpolate(cascade.pt()); - if (doPropagateEfficiency1D) { - efficiencyOmegaErr = hEfficiencyErrOmegaStd1D->Interpolate(cascade.pt()); - efficiencyXiErr = hEfficiencyErrXiStd1D->Interpolate(cascade.pt()); - } - 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.; - } - } - } + double mult = (doProcessIons) ? collision.centFT0C() : collision.centFT0M(); // ion collisions use FT0C for multiplicity, pp uses both - if (doApplyEfficiency2D) { - if constexpr (requires { cascade.topologyChi2(); }) { - efficiencyOmega = hEfficiencyOmegaTra2D->Interpolate(cascade.pt(), mult); - efficiencyXi = hEfficiencyXiTra2D->Interpolate(cascade.pt(), mult); - if (doPropagateEfficiency2D) { - efficiencyOmegaErr = hEfficiencyErrOmegaTra2D->Interpolate(cascade.pt(), mult); - efficiencyXiErr = hEfficiencyErrXiTra2D->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 = hEfficiencyOmegaStd2D->Interpolate(cascade.pt(), mult); - efficiencyXi = hEfficiencyXiStd2D->Interpolate(cascade.pt(), mult); - if (doPropagateEfficiency2D) { - efficiencyOmegaErr = hEfficiencyErrOmegaStd2D->Interpolate(cascade.pt(), mult); - efficiencyXiErr = hEfficiencyErrXiStd2D->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 (doApplyPurity1D) { - if constexpr (requires { cascade.topologyChi2(); }) { - purityOmega = hPurityOmegaTra1D->Interpolate(cascade.pt()); - purityXi = hPurityXiTra1D->Interpolate(cascade.pt()); - if (doPropagatePurity1D) { - purityOmegaErr = hPurityErrOmegaTra1D->Interpolate(cascade.pt()); - purityXiErr = hPurityErrXiTra1D->Interpolate(cascade.pt()); - } - 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 = hPurityOmegaStd1D->Interpolate(cascade.pt()); - purityXi = hPurityXiStd1D->Interpolate(cascade.pt()); - if (doPropagatePurity1D) { - purityOmegaErr = hPurityErrOmegaStd1D->Interpolate(cascade.pt()); - purityXiErr = hPurityErrXiStd1D->Interpolate(cascade.pt()); - } - 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 (doApplyPurity2D) { - if constexpr (requires { cascade.topologyChi2(); }) { - purityOmega = hPurityOmegaTra2D->Interpolate(cascade.pt(), mult); - purityXi = hPurityXiTra2D->Interpolate(cascade.pt(), mult); - if (doPropagatePurity2D) { - purityOmegaErr = hPurityErrOmegaTra2D->Interpolate(cascade.pt(), mult); - purityXiErr = hPurityErrXiTra2D->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 = hPurityOmegaStd2D->Interpolate(cascade.pt(), mult); - purityXi = hPurityXiStd2D->Interpolate(cascade.pt(), mult); - if (doPropagatePurity2D) { - purityOmegaErr = hPurityErrOmegaStd2D->Interpolate(cascade.pt(), mult); - purityXiErr = hPurityErrXiStd2D->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 (doApplyEfficiency1D) { + //* if constexpr (requires { cascade.topologyChi2(); }) { + //* efficiencyOmega = hEfficiencyOmegaTra1D->Interpolate(cascade.pt()); + //* efficiencyXi = hEfficiencyXiTra1D->Interpolate(cascade.pt()); + //* if (doPropagateEfficiency1D) { + //* efficiencyOmegaErr = hEfficiencyErrOmegaTra1D->Interpolate(cascade.pt()); + //* efficiencyXiErr = hEfficiencyErrXiTra1D->Interpolate(cascade.pt()); + //* } + //* 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 = hEfficiencyOmegaStd1D->Interpolate(cascade.pt()); + //* efficiencyXi = hEfficiencyXiStd1D->Interpolate(cascade.pt()); + //* if (doPropagateEfficiency1D) { + //* efficiencyOmegaErr = hEfficiencyErrOmegaStd1D->Interpolate(cascade.pt()); + //* efficiencyXiErr = hEfficiencyErrXiStd1D->Interpolate(cascade.pt()); + //* } + //* 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 (doApplyEfficiency2D) { + //* if constexpr (requires { cascade.topologyChi2(); }) { + //* efficiencyOmega = hEfficiencyOmegaTra2D->Interpolate(cascade.pt(), mult); + //* efficiencyXi = hEfficiencyXiTra2D->Interpolate(cascade.pt(), mult); + //* if (doPropagateEfficiency2D) { + //* efficiencyOmegaErr = hEfficiencyErrOmegaTra2D->Interpolate(cascade.pt(), mult); + //* efficiencyXiErr = hEfficiencyErrXiTra2D->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 = hEfficiencyOmegaStd2D->Interpolate(cascade.pt(), mult); + //* efficiencyXi = hEfficiencyXiStd2D->Interpolate(cascade.pt(), mult); + //* if (doPropagateEfficiency2D) { + //* efficiencyOmegaErr = hEfficiencyErrOmegaStd2D->Interpolate(cascade.pt(), mult); + //* efficiencyXiErr = hEfficiencyErrXiStd2D->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 (doApplyPurity1D) { + //* if constexpr (requires { cascade.topologyChi2(); }) { + //* purityOmega = hPurityOmegaTra1D->Interpolate(cascade.pt()); + //* purityXi = hPurityXiTra1D->Interpolate(cascade.pt()); + //* if (doPropagatePurity1D) { + //* purityOmegaErr = hPurityErrOmegaTra1D->Interpolate(cascade.pt()); + //* purityXiErr = hPurityErrXiTra1D->Interpolate(cascade.pt()); + //* } + //* 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 = hPurityOmegaStd1D->Interpolate(cascade.pt()); + //* purityXi = hPurityXiStd1D->Interpolate(cascade.pt()); + //* if (doPropagatePurity1D) { + //* purityOmegaErr = hPurityErrOmegaStd1D->Interpolate(cascade.pt()); + //* purityXiErr = hPurityErrXiStd1D->Interpolate(cascade.pt()); + //* } + //* 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 (doApplyPurity2D) { + //* if constexpr (requires { cascade.topologyChi2(); }) { + //* purityOmega = hPurityOmegaTra2D->Interpolate(cascade.pt(), mult); + //* purityXi = hPurityXiTra2D->Interpolate(cascade.pt(), mult); + //* if (doPropagatePurity2D) { + //* purityOmegaErr = hPurityErrOmegaTra2D->Interpolate(cascade.pt(), mult); + //* purityXiErr = hPurityErrXiTra2D->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 = hPurityOmegaStd2D->Interpolate(cascade.pt(), mult); + //* purityXi = hPurityXiStd2D->Interpolate(cascade.pt(), mult); + //* if (doPropagatePurity2D) { + //* purityOmegaErr = hPurityErrOmegaStd2D->Interpolate(cascade.pt(), mult); + //* purityXiErr = hPurityErrXiStd2D->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.; + //* } + //* } + //* } // fill multiplicity for events with >=1 cascade if (collision.index() != casccollid) { @@ -693,19 +718,22 @@ struct StrangeCascTrack { double pt = cascade.pt(); double v0cosPA = cascade.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); double casccosPA = cascade.casccosPA(collision.posX(), collision.posY(), collision.posZ()); - double calcDCA = calculateDCA(collision, casccosPA, cascade.x(), cascade.y(), cascade.z()); + ////* double calcDCA = calculateDCA(collision, casccosPA, cascade.x(), cascade.y(), cascade.z()); double bachEta = cascade.bacheloreta(); double negEta = cascade.negativeeta(); double posEta = cascade.positiveeta(); + // fill filters for no cascade selections histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel/Filters/PropDCAxy"), cascade.dcaXYCascToPV()); histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel/Filters/PropDCAz"), cascade.dcaZCascToPV()); - histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel/Filters/CalcDCA"), calcDCA); + ////* histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel/Filters/CalcDCA"), calcDCA); histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel/Filters/BachCosPA"), stdCasc.bachBaryonCosPA()); histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel/Filters/V0CosPA"), v0cosPA); histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel/Filters/CascCosPA"), casccosPA); histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel/Filters/RapidityXi"), cascade.yXi()); + histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel/Filters/PtVsRapidityXi"), pt, cascade.yXi()); histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel/Filters/RapidityOmega"), cascade.yOmega()); + histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel/Filters/PtVsRapidityOmega"), pt, cascade.yOmega()); histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel/Filters/EtaDau"), bachEta); histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel/Filters/EtaDau"), negEta); histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel/Filters/EtaDau"), posEta); @@ -717,7 +745,7 @@ struct StrangeCascTrack { if (isMCTruth(stdCasc, "Xi") || isMCTruth(stdCasc, "Omega")) { histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel-Truth/Filters/PropDCAxy"), cascade.dcaXYCascToPV()); histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel-Truth/Filters/PropDCAz"), cascade.dcaZCascToPV()); - histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel-Truth/Filters/CalcDCA"), calcDCA); + ////* histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel-Truth/Filters/CalcDCA"), calcDCA); histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel-Truth/Filters/BachCosPA"), stdCasc.bachBaryonCosPA()); histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel-Truth/Filters/V0CosPA"), v0cosPA); histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel-Truth/Filters/CascCosPA"), casccosPA); @@ -857,22 +885,26 @@ struct StrangeCascTrack { histos.fill(HIST(TypeNames[Type]) + HIST("/Rec/EvMult"), mult); histos.fill(HIST(TypeNames[Type]) + HIST("/Rec/Filters/PropDCAxy"), cascade.dcaXYCascToPV()); histos.fill(HIST(TypeNames[Type]) + HIST("/Rec/Filters/PropDCAz"), cascade.dcaZCascToPV()); - histos.fill(HIST(TypeNames[Type]) + HIST("/Rec/Filters/CalcDCA"), calcDCA); + ////* histos.fill(HIST(TypeNames[Type]) + HIST("/Rec/Filters/CalcDCA"), calcDCA); histos.fill(HIST(TypeNames[Type]) + HIST("/Rec/Filters/BachCosPA"), stdCasc.bachBaryonCosPA()); histos.fill(HIST(TypeNames[Type]) + HIST("/Rec/Filters/V0CosPA"), v0cosPA); histos.fill(HIST(TypeNames[Type]) + HIST("/Rec/Filters/CascCosPA"), casccosPA); histos.fill(HIST(TypeNames[Type]) + HIST("/Rec/Filters/EtaDau"), bachEta); histos.fill(HIST(TypeNames[Type]) + HIST("/Rec/Filters/EtaDau"), negEta); histos.fill(HIST(TypeNames[Type]) + HIST("/Rec/Filters/EtaDau"), posEta); - if (passedAllSelsXi) + if (passedAllSelsXi) { histos.fill(HIST(TypeNames[Type]) + HIST("/Rec/Filters/RapidityXi"), cascade.yXi()); - if (passedAllSelsOmega) + histos.fill(HIST(TypeNames[Type]) + HIST("/Rec/Filters/PtVsRapidityXi"), pt, cascade.yXi()); + } + if (passedAllSelsOmega) { histos.fill(HIST(TypeNames[Type]) + HIST("/Rec/Filters/RapidityOmega"), cascade.yOmega()); + histos.fill(HIST(TypeNames[Type]) + HIST("/Rec/Filters/PtVsRapidityOmega"), pt, cascade.yOmega()); + } if (fillTruthXi || fillTruthOmega) { histos.fill(HIST(TypeNames[Type]) + HIST("/Rec-Truth/EvMult"), mult); histos.fill(HIST(TypeNames[Type]) + HIST("/Rec-Truth/Filters/PropDCAxy"), cascade.dcaXYCascToPV()); histos.fill(HIST(TypeNames[Type]) + HIST("/Rec-Truth/Filters/PropDCAz"), cascade.dcaZCascToPV()); - histos.fill(HIST(TypeNames[Type]) + HIST("/Rec-Truth/Filters/CalcDCA"), calcDCA); + ////* histos.fill(HIST(TypeNames[Type]) + HIST("/Rec-Truth/Filters/CalcDCA"), calcDCA); histos.fill(HIST(TypeNames[Type]) + HIST("/Rec-Truth/Filters/BachCosPA"), stdCasc.bachBaryonCosPA()); histos.fill(HIST(TypeNames[Type]) + HIST("/Rec-Truth/Filters/V0CosPA"), v0cosPA); histos.fill(HIST(TypeNames[Type]) + HIST("/Rec-Truth/Filters/CascCosPA"), casccosPA); @@ -927,11 +959,49 @@ struct StrangeCascTrack { histos.fill(HIST(TypeNames[Type]) + HIST("/Rec-Truth/Omega"), massOmega, pt, mult); } } + + // statistics - compare gen and reco pt and rapidity + if constexpr (requires { collision.straMCCollisionId(); }) { + if constexpr (requires { stdCasc.has_cascMCCore(); }) { + auto cascmccore = stdCasc.template cascMCCore_as(); + double genPt = cascmccore.ptMC(); + double genYXi = cascmccore.rapidityMC(0); + double genYOmega = cascmccore.rapidityMC(2); + histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel/GenRecPt"), genPt, pt); + histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel/GenRecRapidityXi"), genYXi, cascade.yXi()); + histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel/GenRecRapidityOmega"), genYOmega, cascade.yOmega()); + if (fillTruthXi || fillTruthOmega) + histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel-Truth/GenRecPt"), genPt, pt); + if (passedAllSelsOmega || passedAllSelsXi) { + histos.fill(HIST(TypeNames[Type]) + HIST("/Rec/GenRecPt"), genPt, pt); + if (fillTruthXi || fillTruthOmega) + histos.fill(HIST(TypeNames[Type]) + HIST("/Rec-Truth/GenRecPt"), genPt, pt); + } + if (fillTruthXi) + histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel-Truth/GenRecRapidityXi"), genYXi, cascade.yXi()); + if (passedAllSelsXi) { + histos.fill(HIST(TypeNames[Type]) + HIST("/Rec/GenRecRapidityXi"), genYXi, cascade.yXi()); + if (fillTruthXi) + histos.fill(HIST(TypeNames[Type]) + HIST("/Rec-Truth/GenRecRapidityXi"), genYXi, cascade.yXi()); + } + if (fillTruthOmega) + histos.fill(HIST(TypeNames[Type]) + HIST("/NoSel-Truth/GenRecRapidityOmega"), genYOmega, cascade.yOmega()); + if (passedAllSelsOmega) { + histos.fill(HIST(TypeNames[Type]) + HIST("/Rec/GenRecRapidityOmega"), genYOmega, cascade.yOmega()); + if (fillTruthOmega) + histos.fill(HIST(TypeNames[Type]) + HIST("/Rec-Truth/GenRecRapidityOmega"), genYOmega, cascade.yOmega()); + } + } + } } } void init(InitContext const&) { + // check if cut configuration is valid + if (selCuts.cutMinCascCosPaVsPt->size() != axesConfig.axisPt.value.size() - 1) { + LOGF(fatal, "cutMinCascCosPaVsPt size does not match axisPt binning"); + } // for all events processing histos.add("NoSel-Events/EvCounter", "Event Counter", kTH1D, {{1, 0, 1}}); histos.add("NoSel-Events/PVxy", "PV xy position", kTH2D, {{200, -0.1, 0.1}, {200, -0.1, 0.1}}); @@ -942,30 +1012,37 @@ struct StrangeCascTrack { histos.add("Rec-Events/PVxy", "PV xy position", kTH2D, {{200, -0.1, 0.1}, {200, -0.1, 0.1}}); histos.add("Rec-Events/PVz", "PV z position", kTH1D, {{100, -20, 20}}); histos.add("Rec-Events/Mult", "Multiplicity", kTH1D, {axesConfig.axisMult}); - histos.add("Rec-Events/EvFilter", "Event Filter", kTH1D, {{4, 0, 4}}); + histos.add("Rec-Events/EvFilter", "Event Filter", kTH1D, {{6, 0, 6}}); histos.get(HIST("Rec-Events/EvFilter"))->GetXaxis()->SetBinLabel(1, "INEL>0"); histos.get(HIST("Rec-Events/EvFilter"))->GetXaxis()->SetBinLabel(2, "PVz cut"); histos.get(HIST("Rec-Events/EvFilter"))->GetXaxis()->SetBinLabel(3, "sel8"); - histos.get(HIST("Rec-Events/EvFilter"))->GetXaxis()->SetBinLabel(4, "all"); + histos.get(HIST("Rec-Events/EvFilter"))->GetXaxis()->SetBinLabel(4, "NoSameBunchPileup"); + histos.get(HIST("Rec-Events/EvFilter"))->GetXaxis()->SetBinLabel(5, "IsGoodZvtxFT0VsPV"); + histos.get(HIST("Rec-Events/EvFilter"))->GetXaxis()->SetBinLabel(6, "all"); // for cascade processing static_for<0, 1>([&](auto Type) { // no selections applied histos.add(Form("%s/NoSel/Filters/PropDCAxy", TypeNames[Type].data()), "DCA to xy (propagated)", kTH1D, {axesConfig.axisDCAxy}); histos.add(Form("%s/NoSel/Filters/PropDCAz", TypeNames[Type].data()), "DCA to z (propagated)", kTH1D, {axesConfig.axisDCAz}); - histos.add(Form("%s/NoSel/Filters/CalcDCA", TypeNames[Type].data()), "DCA (calculated)", kTH1D, {axesConfig.axisDCAxy}); + ////* histos.add(Form("%s/NoSel/Filters/CalcDCA", TypeNames[Type].data()), "DCA (calculated)", kTH1D, {axesConfig.axisDCAxy}); histos.add(Form("%s/NoSel/Filters/BachCosPA", TypeNames[Type].data()), "Bachelor cosPA", kTH1D, {{200, -1.0, 1.0}}); histos.add(Form("%s/NoSel/Filters/V0CosPA", TypeNames[Type].data()), "V0 cosPA", kTH1D, {{200, -1.0, 1.0}}); histos.add(Form("%s/NoSel/Filters/CascCosPA", TypeNames[Type].data()), "V0 cosPA", kTH1D, {{200, -1.0, 1.0}}); histos.add(Form("%s/NoSel/Filters/RapidityXi", TypeNames[Type].data()), "y under Xi hypothesis", kTH1D, {{200, -1.0, 1.0}}); + histos.add(Form("%s/NoSel/Filters/PtVsRapidityXi", TypeNames[Type].data()), "pt vs y under Xi hypothesis", kTH2D, {axesConfig.axisPt, {200, -1.0, 1.0}}); histos.add(Form("%s/NoSel/Filters/RapidityOmega", TypeNames[Type].data()), "y under Omega hypothesis", kTH1D, {{200, -1.0, 1.0}}); + histos.add(Form("%s/NoSel/Filters/PtVsRapidityOmega", TypeNames[Type].data()), "pt vs y under Xi hypothesis", kTH2D, {axesConfig.axisPt, {200, -1.0, 1.0}}); histos.add(Form("%s/NoSel/Filters/EtaDau", TypeNames[Type].data()), "|#eta| of dau tracks", kTH1D, {axesConfig.axisEta}); histos.add(Form("%s/NoSel/EvMult", TypeNames[Type].data()), "Multiplicity of events with >=1 cascade", kTH1D, {axesConfig.axisMult}); + histos.add(Form("%s/NoSel/GenRecPt", TypeNames[Type].data()), "Generated vs reconstructed pt", kTH2D, {axesConfig.axisPt, axesConfig.axisPt}); + histos.add(Form("%s/NoSel/GenRecRapidityXi", TypeNames[Type].data()), "Generated vs reconstructed y (Xi)", kTH2D, {{200, -1.0, 1.0}, {200, -1.0, 1.0}}); + histos.add(Form("%s/NoSel/GenRecRapidityOmega", TypeNames[Type].data()), "Generated vs reconstructed y (Omega)", kTH2D, {{200, -1.0, 1.0}, {200, -1.0, 1.0}}); histos.add(Form("%s/NoSel/MassXi", TypeNames[Type].data()), "Invariant mass hypothesis", kTH1D, {axesConfig.axisXiMass}); histos.add(Form("%s/NoSel/MassOmega", TypeNames[Type].data()), "Invariant mass hypothesis", kTH1D, {axesConfig.axisOmegaMass}); // mc truth for no selectrion histos.add(Form("%s/NoSel-Truth/Filters/PropDCAxy", TypeNames[Type].data()), "DCA to xy (propagated)", kTH1D, {axesConfig.axisDCAxy}); histos.add(Form("%s/NoSel-Truth/Filters/PropDCAz", TypeNames[Type].data()), "DCA to z (propagated)", kTH1D, {axesConfig.axisDCAz}); - histos.add(Form("%s/NoSel-Truth/Filters/CalcDCA", TypeNames[Type].data()), "DCA (calculated)", kTH1D, {axesConfig.axisDCAxy}); + ////* histos.add(Form("%s/NoSel-Truth/Filters/CalcDCA", TypeNames[Type].data()), "DCA (calculated)", kTH1D, {axesConfig.axisDCAxy}); histos.add(Form("%s/NoSel-Truth/Filters/BachCosPA", TypeNames[Type].data()), "Bachelor cosPA", kTH1D, {{200, -1.0, 1.0}}); histos.add(Form("%s/NoSel-Truth/Filters/V0CosPA", TypeNames[Type].data()), "V0 cosPA", kTH1D, {{200, -1.0, 1.0}}); histos.add(Form("%s/NoSel-Truth/Filters/CascCosPA", TypeNames[Type].data()), "V0 cosPA", kTH1D, {{200, -1.0, 1.0}}); @@ -973,74 +1050,99 @@ struct StrangeCascTrack { histos.add(Form("%s/NoSel-Truth/Filters/RapidityOmega", TypeNames[Type].data()), "y under Omega hypothesis", kTH1D, {{200, -1.0, 1.0}}); histos.add(Form("%s/NoSel-Truth/Filters/EtaDau", TypeNames[Type].data()), "|#eta| of dau tracks", kTH1D, {axesConfig.axisEta}); histos.add(Form("%s/NoSel-Truth/EvMult", TypeNames[Type].data()), "Multiplicity of events with >=1 cascade", kTH1D, {axesConfig.axisMult}); + histos.add(Form("%s/NoSel-Truth/GenRecPt", TypeNames[Type].data()), "Generated vs reconstructed pt", kTH2D, {axesConfig.axisPt, axesConfig.axisPt}); + histos.add(Form("%s/NoSel-Truth/GenRecRapidityXi", TypeNames[Type].data()), "Generated vs reconstructed y (Xi)", kTH2D, {{200, -1.0, 1.0}, {200, -1.0, 1.0}}); + histos.add(Form("%s/NoSel-Truth/GenRecRapidityOmega", TypeNames[Type].data()), "Generated vs reconstructed y (Omega)", kTH2D, {{200, -1.0, 1.0}, {200, -1.0, 1.0}}); histos.add(Form("%s/NoSel-Truth/MassXi", TypeNames[Type].data()), "Invariant mass hypothesis", kTH1D, {axesConfig.axisXiMass}); histos.add(Form("%s/NoSel-Truth/MassOmega", TypeNames[Type].data()), "Invariant mass hypothesis", kTH1D, {axesConfig.axisOmegaMass}); // xi and omega selection statistics - histos.add(Form("%s/Rec/FiltersXi", TypeNames[Type].data()), "main cascade filters for Xi", kTH1D, {{5, 0, 5}}); - histos.add(Form("%s/Rec/GenFiltersXi", TypeNames[Type].data()), "general cascade filters for Xi", kTH1D, {{9, 0, 9}}); - histos.add(Form("%s/Rec/FiltersOmega", TypeNames[Type].data()), "main cascade filters for Omega", kTH1D, {{6, 0, 6}}); - histos.add(Form("%s/Rec/GenFiltersOmega", TypeNames[Type].data()), "general cascade filters for Omega", kTH1D, {{9, 0, 9}}); histos.add(Form("%s/Rec/Filters/PropDCAxy", TypeNames[Type].data()), "DCA to xy (propagated)", kTH1D, {axesConfig.axisDCAxy}); histos.add(Form("%s/Rec/Filters/PropDCAz", TypeNames[Type].data()), "DCA to z (propagated)", kTH1D, {axesConfig.axisDCAz}); - histos.add(Form("%s/Rec/Filters/CalcDCA", TypeNames[Type].data()), "DCA (calculated)", kTH1D, {axesConfig.axisDCAxy}); + ////* histos.add(Form("%s/Rec/Filters/CalcDCA", TypeNames[Type].data()), "DCA (calculated)", kTH1D, {axesConfig.axisDCAxy}); histos.add(Form("%s/Rec/Filters/BachCosPA", TypeNames[Type].data()), "Bachelor cosPA", kTH1D, {{200, -1.0, 1.0}}); histos.add(Form("%s/Rec/Filters/V0CosPA", TypeNames[Type].data()), "V0 cosPA", kTH1D, {{200, -1.0, 1.0}}); histos.add(Form("%s/Rec/Filters/CascCosPA", TypeNames[Type].data()), "V0 cosPA", kTH1D, {{200, -1.0, 1.0}}); histos.add(Form("%s/Rec/Filters/RapidityXi", TypeNames[Type].data()), "y under Xi hypothesis", kTH1D, {{200, -1.0, 1.0}}); + histos.add(Form("%s/Rec/Filters/PtVsRapidityXi", TypeNames[Type].data()), "pt vs y under Xi hypothesis", kTH2D, {axesConfig.axisPt, {200, -1.0, 1.0}}); histos.add(Form("%s/Rec/Filters/RapidityOmega", TypeNames[Type].data()), "y under Omega hypothesis", kTH1D, {{200, -1.0, 1.0}}); + histos.add(Form("%s/Rec/Filters/PtVsRapidityOmega", TypeNames[Type].data()), "pt vs y under Omega hypothesis", kTH2D, {axesConfig.axisPt, {200, -1.0, 1.0}}); histos.add(Form("%s/Rec/Filters/EtaDau", TypeNames[Type].data()), "|#eta| of dau tracks", kTH1D, {axesConfig.axisEta}); - // passed all applied sels - histos.add(Form("%s/Rec/EvMult", TypeNames[Type].data()), "Multiplicity of events with >=1 cascade", kTH1D, {axesConfig.axisMult}); + histos.add(Form("%s/Rec/EvMult", TypeNames[Type].data()), "Multiplicity of events with >=1 selected cascade", kTH1D, {axesConfig.axisMult}); + histos.add(Form("%s/Rec/GenRecPt", TypeNames[Type].data()), "Generated vs reconstructed pt", kTH2D, {axesConfig.axisPt, axesConfig.axisPt}); + histos.add(Form("%s/Rec/GenRecRapidityXi", TypeNames[Type].data()), "Generated vs reconstructed y (Xi)", kTH2D, {{200, -1.0, 1.0}, {200, -1.0, 1.0}}); + histos.add(Form("%s/Rec/GenRecRapidityOmega", TypeNames[Type].data()), "Generated vs reconstructed y (Omega)", kTH2D, {{200, -1.0, 1.0}, {200, -1.0, 1.0}}); + // selected xi + histos.add(Form("%s/Rec/FiltersXi", TypeNames[Type].data()), "main cascade filters for Xi", kTH1D, {{5, 0, 5}}); + histos.add(Form("%s/Rec/GenFiltersXi", TypeNames[Type].data()), "general cascade filters for Xi", kTH1D, {{9, 0, 9}}); histos.add(Form("%s/Rec/MassXi", TypeNames[Type].data()), "Invariant mass hypothesis", kTH1D, {axesConfig.axisXiMass}); - histos.add(Form("%s/Rec/MassOmega", TypeNames[Type].data()), "Invariant mass hypothesis", kTH1D, {axesConfig.axisOmegaMass}); - histos.add(Form("%s/Rec/Xi", TypeNames[Type].data()), "", kTHnD, {axesConfig.axisXiMass, axesConfig.axisPt, axesConfig.axisMult}); - histos.add(Form("%s/Rec/Omega", TypeNames[Type].data()), "", kTHnD, {axesConfig.axisOmegaMass, axesConfig.axisPt, axesConfig.axisMult}); histos.add(Form("%s/Rec/MassXiMinus", TypeNames[Type].data()), "Invariant mass hypothesis", kTH1D, {axesConfig.axisXiMass}); - histos.add(Form("%s/Rec/MassOmegaMinus", TypeNames[Type].data()), "Invariant mass hypothesis", kTH1D, {axesConfig.axisOmegaMass}); - histos.add(Form("%s/Rec/XiMinus", TypeNames[Type].data()), "", kTHnD, {axesConfig.axisXiMass, axesConfig.axisPt, axesConfig.axisMult}); - histos.add(Form("%s/Rec/OmegaMinus", TypeNames[Type].data()), "", kTHnD, {axesConfig.axisOmegaMass, axesConfig.axisPt, axesConfig.axisMult}); histos.add(Form("%s/Rec/MassXiPlus", TypeNames[Type].data()), "Invariant mass hypothesis", kTH1D, {axesConfig.axisXiMass}); - histos.add(Form("%s/Rec/MassOmegaPlus", TypeNames[Type].data()), "Invariant mass hypothesis", kTH1D, {axesConfig.axisOmegaMass}); + histos.add(Form("%s/Rec/Xi", TypeNames[Type].data()), "", kTHnD, {axesConfig.axisXiMass, axesConfig.axisPt, axesConfig.axisMult}); + histos.add(Form("%s/Rec/XiMinus", TypeNames[Type].data()), "", kTHnD, {axesConfig.axisXiMass, axesConfig.axisPt, axesConfig.axisMult}); histos.add(Form("%s/Rec/XiPlus", TypeNames[Type].data()), "", kTHnD, {axesConfig.axisXiMass, axesConfig.axisPt, axesConfig.axisMult}); + // selected omega + histos.add(Form("%s/Rec/FiltersOmega", TypeNames[Type].data()), "main cascade filters for Omega", kTH1D, {{6, 0, 6}}); + histos.add(Form("%s/Rec/GenFiltersOmega", TypeNames[Type].data()), "general cascade filters for Omega", kTH1D, {{9, 0, 9}}); + histos.add(Form("%s/Rec/MassOmega", TypeNames[Type].data()), "Invariant mass hypothesis", kTH1D, {axesConfig.axisOmegaMass}); + histos.add(Form("%s/Rec/MassOmegaMinus", TypeNames[Type].data()), "Invariant mass hypothesis", kTH1D, {axesConfig.axisOmegaMass}); + histos.add(Form("%s/Rec/MassOmegaPlus", TypeNames[Type].data()), "Invariant mass hypothesis", kTH1D, {axesConfig.axisOmegaMass}); + histos.add(Form("%s/Rec/Omega", TypeNames[Type].data()), "", kTHnD, {axesConfig.axisOmegaMass, axesConfig.axisPt, axesConfig.axisMult}); + histos.add(Form("%s/Rec/OmegaMinus", TypeNames[Type].data()), "", kTHnD, {axesConfig.axisOmegaMass, axesConfig.axisPt, axesConfig.axisMult}); histos.add(Form("%s/Rec/OmegaPlus", TypeNames[Type].data()), "", kTHnD, {axesConfig.axisOmegaMass, axesConfig.axisPt, axesConfig.axisMult}); // mc truth for all passed selections // xi and omega truth selection statistics - histos.add(Form("%s/Rec-Truth/FiltersXi", TypeNames[Type].data()), "main cascade filters for Xi", kTH1D, {{5, 0, 5}}); - histos.add(Form("%s/Rec-Truth/GenFiltersXi", TypeNames[Type].data()), "general cascade filters for Xi", kTH1D, {{9, 0, 9}}); - histos.add(Form("%s/Rec-Truth/FiltersOmega", TypeNames[Type].data()), "main cascade filters for Omega", kTH1D, {{6, 0, 6}}); - histos.add(Form("%s/Rec-Truth/GenFiltersOmega", TypeNames[Type].data()), "general cascade filters for Omega", kTH1D, {{9, 0, 9}}); histos.add(Form("%s/Rec-Truth/Filters/PropDCAxy", TypeNames[Type].data()), "DCA to xy (propagated)", kTH1D, {axesConfig.axisDCAxy}); histos.add(Form("%s/Rec-Truth/Filters/PropDCAz", TypeNames[Type].data()), "DCA to z (propagated)", kTH1D, {axesConfig.axisDCAz}); - histos.add(Form("%s/Rec-Truth/Filters/CalcDCA", TypeNames[Type].data()), "DCA (calculated)", kTH1D, {axesConfig.axisDCAxy}); + ////* histos.add(Form("%s/Rec-Truth/Filters/CalcDCA", TypeNames[Type].data()), "DCA (calculated)", kTH1D, {axesConfig.axisDCAxy}); histos.add(Form("%s/Rec-Truth/Filters/BachCosPA", TypeNames[Type].data()), "Bachelor cosPA", kTH1D, {{200, -1.0, 1.0}}); histos.add(Form("%s/Rec-Truth/Filters/V0CosPA", TypeNames[Type].data()), "V0 cosPA", kTH1D, {{200, -1.0, 1.0}}); histos.add(Form("%s/Rec-Truth/Filters/CascCosPA", TypeNames[Type].data()), "V0 cosPA", kTH1D, {{200, -1.0, 1.0}}); histos.add(Form("%s/Rec-Truth/Filters/RapidityXi", TypeNames[Type].data()), "y under Xi hypothesis", kTH1D, {{200, -1.0, 1.0}}); histos.add(Form("%s/Rec-Truth/Filters/RapidityOmega", TypeNames[Type].data()), "y under Omega hypothesis", kTH1D, {{200, -1.0, 1.0}}); histos.add(Form("%s/Rec-Truth/Filters/EtaDau", TypeNames[Type].data()), "|#eta| of dau tracks", kTH1D, {axesConfig.axisEta}); - // truth that passed all sels histos.add(Form("%s/Rec-Truth/EvMult", TypeNames[Type].data()), "Multiplicity of events with >=1 cascade", kTH1D, {axesConfig.axisMult}); + histos.add(Form("%s/Rec-Truth/GenRecPt", TypeNames[Type].data()), "Generated vs reconstructed pt", kTH2D, {axesConfig.axisPt, axesConfig.axisPt}); + histos.add(Form("%s/Rec-Truth/GenRecRapidityXi", TypeNames[Type].data()), "Generated vs reconstructed y (Xi)", kTH2D, {{200, -1.0, 1.0}, {200, -1.0, 1.0}}); + histos.add(Form("%s/Rec-Truth/GenRecRapidityOmega", TypeNames[Type].data()), "Generated vs reconstructed y (Omega)", kTH2D, {{200, -1.0, 1.0}, {200, -1.0, 1.0}}); + histos.add(Form("%s/Rec-Truth/FiltersXi", TypeNames[Type].data()), "main cascade filters for Xi", kTH1D, {{5, 0, 5}}); + histos.add(Form("%s/Rec-Truth/GenFiltersXi", TypeNames[Type].data()), "general cascade filters for Xi", kTH1D, {{9, 0, 9}}); histos.add(Form("%s/Rec-Truth/MassXi", TypeNames[Type].data()), "Invariant mass hypothesis", kTH1D, {axesConfig.axisXiMass}); + histos.add(Form("%s/Rec-Truth/Xi", TypeNames[Type].data()), "", kTHnD, {axesConfig.axisXiMass, axesConfig.axisPt, axesConfig.axisMult}); + histos.add(Form("%s/Rec-Truth/FiltersOmega", TypeNames[Type].data()), "main cascade filters for Omega", kTH1D, {{6, 0, 6}}); + histos.add(Form("%s/Rec-Truth/GenFiltersOmega", TypeNames[Type].data()), "general cascade filters for Omega", kTH1D, {{9, 0, 9}}); histos.add(Form("%s/Rec-Truth/MassOmega", TypeNames[Type].data()), "Invariant mass hypothesis", kTH1D, {axesConfig.axisOmegaMass}); histos.add(Form("%s/Rec-Truth/Omega", TypeNames[Type].data()), "", kTHnD, {axesConfig.axisOmegaMass, axesConfig.axisPt, axesConfig.axisMult}); - histos.add(Form("%s/Rec-Truth/Xi", TypeNames[Type].data()), "", kTHnD, {axesConfig.axisXiMass, axesConfig.axisPt, axesConfig.axisMult}); }); // for MC-specific processing + // all generated events: histos.add("MC/Gen/EvCounter", "Event Counter", kTH1D, {{1, 0, 1}}); - histos.add("MC/Gen/Xi", "Xi", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated Xis - histos.add("MC/Gen/Omega", "Omega", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated Omegas - histos.add("MC/Gen/PrimaryXi", "Xi primaries", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Xis - histos.add("MC/Gen/PrimaryOmega", "Omega primaries", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Omegas - histos.add("MC/Gen/PrimaryXiRapidity", "Xi primaries in |y|", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Xis in selected rapidity range - histos.add("MC/Gen/PrimaryOmegaRapidity", "Omega primaries in |y|", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Omegas in selected rapidity range - histos.add("MC/Gen/PrimaryXiMinus", "Xi- primaries", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Xis - histos.add("MC/Gen/PrimaryOmegaMinus", "Omega- primaries", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Omegas - histos.add("MC/Gen/PrimaryXiMinusRapidity", "Xi- primaries in |y| ", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Xis in selected rapidity range - histos.add("MC/Gen/PrimaryOmegaMinusRapidity", "Omega- primaries in |y|", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Omegas in selected rapidity range - histos.add("MC/Gen/PrimaryXiPlus", "Xi+ primaries", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Xis - histos.add("MC/Gen/PrimaryOmegaPlus", "Omega+ primaries", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Omegas - histos.add("MC/Gen/PrimaryXiPlusRapidity", "Xi+ primaries in |y|", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Xis in selected rapidity range - histos.add("MC/Gen/PrimaryOmegaPlusRapidity", "Omega+ primaries in |y|", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Omegas in selected rapidity range + histos.add("MC/Gen/Mult", "Event charged multiplicty (all gen events)", kTH1D, {{3200, 0, 3200}}); + histos.add("MC/Gen/RapidityXi", "Generated y (Xi)", kTH1D, {{200, -1.0, 1.0}}); + histos.add("MC/Gen/Xi", "Xi", kTH2D, {{3200, 0, 3200}, {100, 0.0, 10.0}}); // generated Xis + histos.add("MC/Gen/PrimaryXi", "Xi primaries", kTH2D, {{3200, 0, 3200}, {100, 0.0, 10.0}}); // generated primary Xis + histos.add("MC/Gen/PrimaryXiRapidity", "Xi primaries in |y|", kTH2D, {{3200, 0, 3200}, {100, 0.0, 10.0}}); // generated primary Xis in selected rapidity range + histos.add("MC/Gen/PrimaryXiMinusRapidity", "Xi- primaries in |y| ", kTH2D, {{3200, 0, 3200}, {100, 0.0, 10.0}}); // generated primary Xi-minus in selected rapidity range + histos.add("MC/Gen/PrimaryXiPlusRapidity", "Xi+ primaries in |y|", kTH2D, {{3200, 0, 3200}, {100, 0.0, 10.0}}); // generated primary Xi-plus in selected rapidity range + histos.add("MC/Gen/RapidityOmega", "Generated y (Omega)", kTH1D, {{200, -1.0, 1.0}}); + histos.add("MC/Gen/Omega", "Omega", kTH2D, {{3200, 0, 3200}, {100, 0.0, 10.0}}); // generated Omegas + histos.add("MC/Gen/PrimaryOmega", "Omega primaries", kTH2D, {{3200, 0, 3200}, {100, 0.0, 10.0}}); // generated primary Omegas + histos.add("MC/Gen/PrimaryOmegaRapidity", "Omega primaries in |y|", kTH2D, {{3200, 0, 3200}, {100, 0.0, 10.0}}); // generated primary Omegas in selected rapidity range + histos.add("MC/Gen/PrimaryOmegaMinusRapidity", "Omega- primaries in |y| ", kTH2D, {{3200, 0, 3200}, {100, 0.0, 10.0}}); // generated primary Omega-minus in selected rapidity range + histos.add("MC/Gen/PrimaryOmegaPlusRapidity", "Omega+ primaries in |y|", kTH2D, {{3200, 0, 3200}, {100, 0.0, 10.0}}); // generated primary Omega-plus in selected rapidity range + // generated events with reco >= 1: + histos.add("MC/EvRec/EvCounter", "Event Counter", kTH1D, {{1, 0, 1}}); + histos.add("MC/EvRec/Mult", "Event charged multiplicty (gen events with reco>=1)", kTH1D, {{3200, 0, 3200}}); + histos.add("MC/EvRec/MultCent", "Gen multiplicity vs reco centrality", kTH2D, {{100, 0, 100}, {3200, 0, 3200}}); + histos.add("MC/EvRec/Xi", "Xi", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated Xis in reco>=1 + histos.add("MC/EvRec/PrimaryXi", "Xi primaries", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Xis in reco>=1 + histos.add("MC/EvRec/PrimaryXiRapidity", "Xi primaries in |y|", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Xis in selected rapidity range in reco>=1 + histos.add("MC/EvRec/PrimaryXiMinusRapidity", "Xi- primaries in |y| ", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Xi-minus in selected rapidity range in reco>=1 + histos.add("MC/EvRec/PrimaryXiPlusRapidity", "Xi+ primaries in |y|", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Xi-plus in selected rapidity range in reco>=1 + histos.add("MC/EvRec/Omega", "Omega", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated Omegas in reco>=1 + histos.add("MC/EvRec/PrimaryOmega", "Omega primaries", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Omegas in reco>=1 + histos.add("MC/EvRec/PrimaryOmegaRapidity", "Omega primaries in |y|", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Omegas in selected rapidity range in reco>=1 + histos.add("MC/EvRec/PrimaryOmegaMinusRapidity", "Omega- primaries in |y|", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Omega-minus in selected rapidity range in reco>=1 + histos.add("MC/EvRec/PrimaryOmegaPlusRapidity", "Omega+ primaries in |y|", kTH2D, {axesConfig.axisPt, axesConfig.axisMult}); // generated primary Omega-plus in selected rapidity range in reco>=1 // label filter statistic bins for standard cascs histos.get(HIST("Standard/Rec/FiltersXi"))->GetXaxis()->SetBinLabel(1, "p_{T}"); histos.get(HIST("Standard/Rec/FiltersXi"))->GetXaxis()->SetBinLabel(2, "gen"); @@ -1164,66 +1266,107 @@ struct StrangeCascTrack { void processDerivedData(DerCollisionWMults::iterator const& collision, DerCascDatas const& allCascs, DerTraCascDatas const& traCascs, DauTracks const&) { fillEvents(collision); // save info about all processed events - if (doApplyEfficiency1D || doApplyPurity1D || doApplyEfficiency2D || doApplyPurity2D) { - initEfficiencyFromCCDB(collision.runNumber(), collision.timestamp()); - } + //* if (doApplyEfficiency1D || doApplyPurity1D || doApplyEfficiency2D || doApplyPurity2D) { + //* initEfficiencyFromCCDB(collision.runNumber(), collision.timestamp()); + //* } if (isValidEvent(collision)) { histos.fill(HIST("Rec-Events/EvCounter"), 0.5); histos.fill(HIST("Rec-Events/PVxy"), collision.posX(), collision.posY()); histos.fill(HIST("Rec-Events/PVz"), collision.posZ()); - double mult = (doProcesspp || doProcesspO) ? collision.centFT0M() : collision.centFT0C(); + double mult = (doProcessIons) ? collision.centFT0C() : collision.centFT0M(); histos.fill(HIST("Rec-Events/Mult"), mult); analyseCascs(collision, allCascs); // process all cascades analyseCascs(collision, traCascs); // process tracked cascades } } - void processDerivedMCGen(aod::StraMCCollisions const& genColls, DerMCGenCascades const& genCascs, soa::Join const& recColls) + void processDerivedMCGen(soa::Join const& genColls, DerMCGenCascades const& genCascs, soa::Join const& recColls) { for (auto const& genColl : genColls) { histos.fill(HIST("MC/Gen/EvCounter"), 0.5); // generated events statistics - // (for efficiency calculation) only take reconstructed events + double genMult = genColl.multMCNParticlesEta05(); + histos.fill(HIST("MC/Gen/Mult"), genMult); + // loop over generated cascades to obtain pt-genMult distribution + auto slicedGenCascs = genCascs.sliceBy(cascsPerMcCollision, genColl.globalIndex()); + for (auto const& casc : slicedGenCascs) { + if (casc.straMCCollisionId() != genColl.index()) + continue; // safety check for correct slicing + // fill generated Xi and Omega (all) for statistics + double cascPt = casc.ptMC(); + if (isValidPDG(casc, "Xi")) + histos.fill(HIST("MC/Gen/Xi"), genMult, cascPt); + if (isValidPDG(casc, "Omega")) + histos.fill(HIST("MC/Gen/Omega"), genMult, cascPt); + if (casc.isPhysicalPrimary()) { + // fill generated primary Xi and Omega for statistics + if (isValidPDG(casc, "Xi")) + histos.fill(HIST("MC/Gen/PrimaryXi"), genMult, cascPt); + if (isValidPDG(casc, "Omega")) + histos.fill(HIST("MC/Gen/PrimaryOmega"), genMult, cascPt); + // fill generated primary Xi within rapidity for corrections + if (isValidPDG(casc, "Xi") && std::abs(casc.rapidityMC(0)) < selCuts.cutRapidity) { + histos.fill(HIST("MC/Gen/PrimaryXiRapidity"), genMult, cascPt); + if (casc.pdgCode() == PDG_t::kXiMinus) + histos.fill(HIST("MC/Gen/PrimaryXiMinusRapidity"), genMult, cascPt); + if (casc.pdgCode() == PDG_t::kXiPlusBar) + histos.fill(HIST("MC/Gen/PrimaryXiPlusRapidity"), genMult, cascPt); + } + // fill generated primary Omega within rapidity for corrections + if (isValidPDG(casc, "Omega") && std::abs(casc.rapidityMC(2)) < selCuts.cutRapidity) { + histos.fill(HIST("MC/Gen/PrimaryOmegaRapidity"), genMult, cascPt); + if (casc.pdgCode() == PDG_t::kOmegaMinus) + histos.fill(HIST("MC/Gen/PrimaryOmegaMinusRapidity"), genMult, cascPt); + if (casc.pdgCode() == PDG_t::kOmegaPlusBar) + histos.fill(HIST("MC/Gen/PrimaryOmegaPlusRapidity"), genMult, cascPt); + } + } + // fill gen rapidity for statistics + if (isValidPDG(casc, "Xi")) + histos.fill(HIST("MC/Gen/RapidityXi"), casc.rapidityMC(0)); + if (isValidPDG(casc, "Omega")) + histos.fill(HIST("MC/Gen/RapidityOmega"), casc.rapidityMC(2)); + } + // look at rec. ev. >= 1 for corrections auto slicedRecColls = recColls.sliceBy(perMcCollision, genColl.globalIndex()); + if (slicedRecColls.size() > 0) { + histos.fill(HIST("MC/EvRec/EvCounter"), 0.5); + histos.fill(HIST("MC/EvRec/Mult"), genMult); + } for (auto const& recColl : slicedRecColls) { - double cascMult = (doProcesspp || doProcesspO) ? recColl.centFT0M() : recColl.centFT0C(); int64_t genCollId = recColl.straMCCollisionId(); - for (auto const& casc : genCascs) { + if (genColl.index() != genCollId) + continue; // safety check for correct slicing + // fill multiplicity-centrality distribution + double recoMult = (doProcessIons) ? recColl.centFT0C() : recColl.centFT0M(); + histos.fill(HIST("MC/EvRec/MultCent"), recoMult, genMult); + // look at generated cascades within reconstructed events + for (auto const& casc : slicedGenCascs) { if (casc.straMCCollisionId() != genCollId) - continue; // safety check - double cascPt = std::sqrt(std::pow(casc.pxMC(), 2) + std::pow(casc.pyMC(), 2)); + continue; // check for cascades belonging to the correct reco collision + double cascPt = casc.ptMC(); if (isValidPDG(casc, "Xi")) - histos.fill(HIST("MC/Gen/Xi"), cascPt, cascMult); + histos.fill(HIST("MC/EvRec/Xi"), cascPt, recoMult); if (isValidPDG(casc, "Omega")) - histos.fill(HIST("MC/Gen/Omega"), cascPt, cascMult); + histos.fill(HIST("MC/EvRec/Omega"), cascPt, recoMult); if (casc.isPhysicalPrimary()) { if (isValidPDG(casc, "Xi")) { - histos.fill(HIST("MC/Gen/PrimaryXi"), cascPt, cascMult); - if (std::abs(casc.rapidityMC(0)) < selCuts.cutRapidity) - histos.fill(HIST("MC/Gen/PrimaryXiRapidity"), cascPt, cascMult); - if (casc.pdgCode() == PDG_t::kXiMinus) { - histos.fill(HIST("MC/Gen/PrimaryXiMinus"), cascPt, cascMult); - if (std::abs(casc.rapidityMC(0)) < selCuts.cutRapidity) - histos.fill(HIST("MC/Gen/PrimaryXiMinusRapidity"), cascPt, cascMult); - } - if (casc.pdgCode() == PDG_t::kXiPlusBar) { - histos.fill(HIST("MC/Gen/PrimaryXiPlus"), cascPt, cascMult); - if (std::abs(casc.rapidityMC(0)) < selCuts.cutRapidity) - histos.fill(HIST("MC/Gen/PrimaryXiPlusRapidity"), cascPt, cascMult); + histos.fill(HIST("MC/EvRec/PrimaryXi"), cascPt, recoMult); + if (std::abs(casc.rapidityMC(0)) < selCuts.cutRapidity) { + histos.fill(HIST("MC/EvRec/PrimaryXiRapidity"), cascPt, recoMult); + if (casc.pdgCode() == PDG_t::kXiMinus) + histos.fill(HIST("MC/EvRec/PrimaryXiMinusRapidity"), cascPt, recoMult); + if (casc.pdgCode() == PDG_t::kXiPlusBar) + histos.fill(HIST("MC/EvRec/PrimaryXiPlusRapidity"), cascPt, recoMult); } } if (isValidPDG(casc, "Omega")) { - histos.fill(HIST("MC/Gen/PrimaryOmega"), cascPt, cascMult); - if (std::abs(casc.rapidityMC(2)) < selCuts.cutRapidity) - histos.fill(HIST("MC/Gen/PrimaryOmegaRapidity"), cascPt, cascMult); - if (casc.pdgCode() == PDG_t::kOmegaMinus) { - histos.fill(HIST("MC/Gen/PrimaryOmegaMinus"), cascPt, cascMult); - if (std::abs(casc.rapidityMC(2)) < selCuts.cutRapidity) - histos.fill(HIST("MC/Gen/PrimaryOmegaMinusRapidity"), cascPt, cascMult); - } - if (casc.pdgCode() == PDG_t::kOmegaPlusBar) { - histos.fill(HIST("MC/Gen/PrimaryOmegaPlus"), cascPt, cascMult); - if (std::abs(casc.rapidityMC(2)) < selCuts.cutRapidity) - histos.fill(HIST("MC/Gen/PrimaryOmegaPlusRapidity"), cascPt, cascMult); + histos.fill(HIST("MC/EvRec/PrimaryOmega"), cascPt, recoMult); + if (std::abs(casc.rapidityMC(2)) < selCuts.cutRapidity) { + histos.fill(HIST("MC/EvRec/PrimaryOmegaRapidity"), cascPt, recoMult); + if (casc.pdgCode() == PDG_t::kOmegaMinus) + histos.fill(HIST("MC/EvRec/PrimaryOmegaMinusRapidity"), cascPt, recoMult); + if (casc.pdgCode() == PDG_t::kOmegaPlusBar) + histos.fill(HIST("MC/EvRec/PrimaryOmegaPlusRapidity"), cascPt, recoMult); } } } @@ -1235,14 +1378,14 @@ struct StrangeCascTrack { 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 (doApplyEfficiency1D || doApplyPurity1D || doApplyEfficiency2D || doApplyPurity2D) { - initEfficiencyFromCCDB(collision.runNumber(), collision.timestamp()); - } + //* if (doApplyEfficiency1D || doApplyPurity1D || doApplyEfficiency2D || doApplyPurity2D) { + //* initEfficiencyFromCCDB(collision.runNumber(), collision.timestamp()); + //* } if (isValidEvent(collision)) { histos.fill(HIST("Rec-Events/EvCounter"), 0.5); histos.fill(HIST("Rec-Events/PVxy"), collision.posX(), collision.posY()); histos.fill(HIST("Rec-Events/PVz"), collision.posZ()); - double mult = (doProcesspp || doProcesspO) ? collision.centFT0M() : collision.centFT0C(); + double mult = (doProcessIons) ? collision.centFT0C() : collision.centFT0M(); histos.fill(HIST("Rec-Events/Mult"), mult); analyseCascs(collision, allCascs); // process all cascades analyseCascs(collision, traCascs); // process tracked cascades