From 94ae72fc226f11b0eced3233ecedcd7c167660cb Mon Sep 17 00:00:00 2001 From: Victor Date: Sat, 23 Aug 2025 19:15:48 +0200 Subject: [PATCH] [PWGCF] DptDpt -- Configurable rationalization Also - multiplicities correlations based collision exclusion --- PWGCF/TableProducer/dptDptFilter.cxx | 143 ++++++++++++++--- PWGCF/TableProducer/dptDptFilter.h | 231 +++++++++++++++++++++------ 2 files changed, 301 insertions(+), 73 deletions(-) diff --git a/PWGCF/TableProducer/dptDptFilter.cxx b/PWGCF/TableProducer/dptDptFilter.cxx index 2fe925102b8..0aab214d083 100644 --- a/PWGCF/TableProducer/dptDptFilter.cxx +++ b/PWGCF/TableProducer/dptDptFilter.cxx @@ -49,6 +49,7 @@ #include #include #include +#include #include using namespace o2; @@ -112,18 +113,32 @@ static const std::vector beforeAfterSufix = {"B", "A"}; static constexpr float MultiplicityUpperLimitBase[11][8] = { /* T0A, T0C, T0M, V0A, V0C, V0M, tracks, PV contr. tracks */ {30e3f, 10e3f, 40e3f, 50e3f, 70e3f, 100e3f, 10e3f, 400.0f}, /* no system */ - {30e3f, 10e3f, 40e3f, 50e3f, 70e3f, 100e3f, 10e3f, 400.0f}, /* pp Run2 */ + {30e3f, 10e3f, 40e3f, 50e3f, 70e3f, 50e3f, 5e3f, 400.0f}, /* pp Run2 */ {30e3f, 10e3f, 40e3f, 50e3f, 70e3f, 100e3f, 10e3f, 400.0f}, /* pPb Run2 */ {30e3f, 10e3f, 40e3f, 50e3f, 70e3f, 100e3f, 10e3f, 400.0f}, /* Pbp Run2 */ {30e3f, 10e3f, 40e3f, 50e3f, 70e3f, 100e3f, 10e3f, 400.0f}, /* PbPb Run2 */ {30e3f, 10e3f, 40e3f, 50e3f, 70e3f, 100e3f, 10e3f, 400.0f}, /* XeXe Run2 */ - {30e3f, 10e3f, 40e3f, 50e3f, 70e3f, 100e3f, 10e3f, 400.0f}, /* pp Run3 */ + {30e3f, 10e3f, 40e3f, 50e3f, 70e3f, 50e3f, 5e3f, 400.0f}, /* pp Run3 */ {30e3f, 10e3f, 40e3f, 50e3f, 70e3f, 100e3f, 10e3f, 400.0f}, /* PbPb Run3 */ - {30e3f, 10e3f, 40e3f, 50e3f, 70e3f, 100e3f, 10e3f, 400.0f}, /* NeNe Run3 */ - {30e3f, 10e3f, 40e3f, 50e3f, 70e3f, 100e3f, 10e3f, 400.0f}, /* OO Run3 */ - {30e3f, 10e3f, 40e3f, 50e3f, 70e3f, 100e3f, 10e3f, 400.0f} /* pO Run3 */ + {30e3f, 10e3f, 40e3f, 50e3f, 70e3f, 50e3f, 5e3f, 400.0f}, /* NeNe Run3 */ + {30e3f, 10e3f, 40e3f, 50e3f, 70e3f, 50e3f, 5e3f, 400.0f}, /* OO Run3 */ + {30e3f, 10e3f, 40e3f, 50e3f, 70e3f, 50e3f, 5e3f, 400.0f} /* pO Run3 */ }; +/* helpers for the multiplicity/centrality correlations exclusion formulae */ +static const std::string multiplicityCentralityCorrelationsFormulaBase[11][1] = { + /* no system */ {""}, + /* pp Run2 */ {""}, + /* pPb Run2 */ {""}, + /* Pbp Run2 */ {""}, + /* PbPb Run2 */ {""}, + /* XeXe Run2 */ {""}, + /* pp Run3 */ {""}, + /* PbPb Run3 */ {""}, + /* NeNe Run3 */ {"((221.987-0.706581*[CT0C]+0.00906828*[CT0C]*[CT0C]-23.5275*sqrt([CT0C])-2.5*(38.5823+2.00054*[CT0C]+0.0125088*[CT0C]*[CT0C]-8.47825*sqrt([CT0C])-0.273493*[CT0C]*sqrt([CT0C])))<=[MNGLTRK])&&((-0.0729502+0.79403*[MNPVC]+5.0*(0.032005+0.00372503*[MNPVC]-7.18977e-06*[MNPVC]*[MNPVC]+0.412714*sqrt([MNPVC])))>[MNGLTRK])&&((-0.0729502+0.79403*[MNPVC]-5.0*(0.032005+0.00372503*[MNPVC]-7.18977e-06*[MNPVC]*[MNPVC]+0.412714*sqrt([MNPVC])))<=[MNGLTRK])"}, + /* OO Run3 */ {"((180.584-0.211625*[CT0C]+0.00568101*[CT0C]*[CT0C]-20.9838*sqrt([CT0C])-2.5*(32.665+1.67341*[CT0C]+0.0113878*[CT0C]*[CT0C]-6.83271*sqrt([CT0C])-0.239428*[CT0C]*sqrt([CT0C])))<=[MNGLTRK])&&((-0.0533229+0.79235*[MNPVC]+5.0*(0.031512+0.0045874*[MNPVC]-1.06374e-05*[MNPVC]*[MNPVC]+0.407526*sqrt([MNPVC])))>[MNGLTRK])&&((-0.0533229+0.79235*[MNPVC]-5.0*(0.031512+0.0045874*[MNPVC]-1.06374e-05*[MNPVC]*[MNPVC]+0.407526*sqrt([MNPVC])))<=[MNGLTRK])"}, + /* pO Run3 */ {""}}; + //============================================================================================ // The DptDptFilter histogram objects // TODO: consider registering in the histogram registry @@ -393,6 +408,79 @@ struct Multiplicity { } }; +////////////////////////////////////////////////////////////////////////////////// +/// Centrality/multiplicity correlations exclusion +/// TODO: adaptation to Run 1/Run 2 +////////////////////////////////////////////////////////////////////////////////// +template +inline void storeMultiplicitiesAndCentralities(CollisionObject const& collision, TrackListObject const& tracks) +{ + /* only do it if we are tracking the centrality / multiplicity correlations */ + if (useCentralityMultiplicityCorrelationsExclusion) { + int nGlobalTracks = 0; + for (auto const& track : tracks) { + if (track.isGlobalTrack()) { + nGlobalTracks++; + } + } + for (CentMultCorrelationsParams ipar = CentMultCorrelationsMT0A; ipar < CentMultCorrelationsNOOFPARAMS; ++ipar) { + switch (ipar) { + case CentMultCorrelationsMT0A: + collisionMultiplicityCentralityObservables[ipar] = collision.multFT0A(); + break; + case CentMultCorrelationsMT0C: + collisionMultiplicityCentralityObservables[ipar] = collision.multFT0C(); + break; + case CentMultCorrelationsMT0M: + collisionMultiplicityCentralityObservables[ipar] = collision.multFT0M(); + break; + case CentMultCorrelationsMV0A: + collisionMultiplicityCentralityObservables[ipar] = collision.multFV0A(); + break; + case CentMultCorrelationsMV0C: + collisionMultiplicityCentralityObservables[ipar] = collision.multFV0C(); + break; + case CentMultCorrelationsMV0M: + collisionMultiplicityCentralityObservables[ipar] = collision.multFV0M(); + break; + case CentMultCorrelationsMNGLTRK: + collisionMultiplicityCentralityObservables[ipar] = nGlobalTracks; + break; + case CentMultCorrelationsMNPVC: + collisionMultiplicityCentralityObservables[ipar] = collision.multNTracksPV(); + break; + case CentMultCorrelationsCT0A: + if constexpr (framework::has_type_v) { + collisionMultiplicityCentralityObservables[ipar] = collision.centFT0A(); + } + break; + case CentMultCorrelationsCT0C: + if constexpr (framework::has_type_v) { + collisionMultiplicityCentralityObservables[ipar] = collision.centFT0C(); + } + break; + case CentMultCorrelationsCT0M: + if constexpr (framework::has_type_v) { + collisionMultiplicityCentralityObservables[ipar] = collision.centFT0M(); + } + break; + case CentMultCorrelationsCV0A: + if constexpr (framework::has_type_v) { + collisionMultiplicityCentralityObservables[ipar] = collision.centFV0A(); + } + break; + case CentMultCorrelationsCNTPV: + if constexpr (framework::has_type_v) { + collisionMultiplicityCentralityObservables[ipar] = collision.centNTPV(); + } + break; + default: + break; + } + } + } +} + ////////////////////////////////////////////////////////////////////////////// // The filter class ////////////////////////////////////////////////////////////////////////////// @@ -414,7 +502,10 @@ struct DptDptFilter { Configurable minOrbit{"minOrbit", -1, "Lowest orbit to track"}; Configurable maxOrbit{"maxOrbit", INT64_MAX, "Highest orbit to track"}; Configurable rctSource{"rctSource", "None", "RCT selection source: None,CBT,CBT_hadronPID,CBT_electronPID,CBT_calo,CBT_muon,CBT_muon_glo. Default: None"}; - Configurable> multiplicityUpperLimit{"multiplicityUpperLimit", {&MultiplicityUpperLimitBase[0][0], 11, 8, {systemExternalNamesMap.at(0), systemExternalNamesMap.at(1), systemExternalNamesMap.at(2), systemExternalNamesMap.at(3), systemExternalNamesMap.at(4), systemExternalNamesMap.at(5), systemExternalNamesMap.at(6), systemExternalNamesMap.at(7), systemExternalNamesMap.at(8), systemExternalNamesMap.at(9), systemExternalNamesMap.at(10)}, {multiplicitySourceConfigNamesMap.at(0), multiplicitySourceConfigNamesMap.at(1), multiplicitySourceConfigNamesMap.at(2), multiplicitySourceConfigNamesMap.at(3), multiplicitySourceConfigNamesMap.at(4), multiplicitySourceConfigNamesMap.at(5), multiplicitySourceConfigNamesMap.at(6), multiplicitySourceConfigNamesMap.at(7)}}, "Upper limits for the multiplicity observables"}; +#define SYSTEMNAME(sysid) systemExternalNamesMap.at(sysid).data() +#define MULTSRCNAME(msrcid) multiplicitySourceConfigNamesMap.at(msrcid).data() + Configurable> multiplicityUpperLimit{"multiplicityUpperLimit", {MultiplicityUpperLimitBase[0], 11, 8, {SYSTEMNAME(0), SYSTEMNAME(1), SYSTEMNAME(2), SYSTEMNAME(3), SYSTEMNAME(4), SYSTEMNAME(5), SYSTEMNAME(6), SYSTEMNAME(7), SYSTEMNAME(8), SYSTEMNAME(9), SYSTEMNAME(10)}, {MULTSRCNAME(0), MULTSRCNAME(1), MULTSRCNAME(2), MULTSRCNAME(3), MULTSRCNAME(4), MULTSRCNAME(5), MULTSRCNAME(6), MULTSRCNAME(7)}}, "Upper limits for the multiplicity observables"}; + Configurable> multiplicitiesExclusionFormula{"multiplicitiesExclusionFormula", {multiplicityCentralityCorrelationsFormulaBase[0], 11, 1, {SYSTEMNAME(0), SYSTEMNAME(1), SYSTEMNAME(2), SYSTEMNAME(3), SYSTEMNAME(4), SYSTEMNAME(5), SYSTEMNAME(6), SYSTEMNAME(7), SYSTEMNAME(8), SYSTEMNAME(9), SYSTEMNAME(10)}, {"Exclusion formula"}}, "Formula for excluding outliers of the multiplicities correlations. Use any parameter from centMultCorrelationsParamsMap"}; struct : ConfigurableGroup { std::string prefix = "cfgEventSelection.occupancySelection"; Configurable occupancyEstimation{"occupancyEstimation", "None", "Occupancy estimation: None, Tracks, FT0C. Default None"}; @@ -475,7 +566,7 @@ struct DptDptFilter { } else if (doprocessOnTheFlyGeneratorLevel) { fCentMultEstimator = CentMultFV0A; } else { - fCentMultEstimator = getCentMultEstimator(cfgCentMultEstimator); + fCentMultEstimator = getCentMultEstimator(cfgCentMultEstimator.value.c_str()); } /* RCT information usage */ if (cfgEventSelection.rctSource.value == "None") { @@ -487,19 +578,22 @@ struct DptDptFilter { } /* the occupancy selection */ - fOccupancyEstimation = getOccupancyEstimator(cfgEventSelection.occupancySelection.occupancyEstimation); + fOccupancyEstimation = getOccupancyEstimator(cfgEventSelection.occupancySelection.occupancyEstimation.value.c_str()); fMinOccupancy = cfgEventSelection.occupancySelection.minOccupancy; fMaxOccupancy = cfgEventSelection.occupancySelection.maxOccupancy; /* the trigger selection */ - triggerSelectionFlags = getTriggerSelection(cfgTriggSel); + triggerSelectionFlags = getTriggerSelection(cfgTriggSel.value.c_str()); traceCollId0 = cfgTraceCollId0; /* if the system type is not known at this time, we have to put the initialization somewhere else */ - fSystem = getSystemType(cfgSystem); + fSystem = getSystemType(cfgSystem.value.c_str()); fLhcRun = multRunForSystemMap.at(fSystem); fDataType = getDataType(cfgDataType); + /* the multiplicities outliers exclusion */ + multiplicityCentralityCorrelationsExclusion = getExclusionFormula(cfgEventSelection.multiplicitiesExclusionFormula->getData()[fSystem][0].c_str()); + /* create the output list which will own the task histograms */ TList* fOutputList = new TList(); fOutputList->SetOwner(true); @@ -509,18 +603,18 @@ struct DptDptFilter { /* create the reconstructed data histograms */ fhEventSelection = new TH1D("EventSelection", ";;counts", CollSelNOOFFLAGS, -0.5f, static_cast(CollSelNOOFFLAGS) - 0.5f); for (int ix = 0; ix < CollSelNOOFFLAGS; ++ix) { - fhEventSelection->GetXaxis()->SetBinLabel(ix + 1, collisionSelectionExternalNamesMap.at(ix).c_str()); + fhEventSelection->GetXaxis()->SetBinLabel(ix + 1, collisionSelectionExternalNamesMap.at(ix).data()); } fhTriggerSelection = new TH1D("TriggerSelection", ";;counts", TriggSelNOOFTRIGGERS, -0.5f, static_cast(TriggSelNOOFTRIGGERS) - 0.5f); for (int ix = 0; ix < TriggSelNOOFTRIGGERS; ++ix) { - fhTriggerSelection->GetXaxis()->SetBinLabel(ix + 1, TString::Format("#color[%d]{%s}", triggerSelectionFlags.test(ix) ? 2 : 1, triggerSelectionExternalNamesMap.at(ix).c_str()).Data()); + fhTriggerSelection->GetXaxis()->SetBinLabel(ix + 1, TString::Format("#color[%d]{%s}", triggerSelectionFlags.test(ix) ? 2 : 1, triggerSelectionExternalNamesMap.at(ix).data()).Data()); } fhTriggerSelectionCorrelations = new TH2D("TriggerSelectionCorrelations", ";;;counts", TriggSelNOOFTRIGGERS, -0.5f, static_cast(TriggSelNOOFTRIGGERS) - 0.5f, TriggSelNOOFTRIGGERS, -0.5f, static_cast(TriggSelNOOFTRIGGERS) - 0.5f); for (int ixy = 0; ixy < TriggSelNOOFTRIGGERS; ++ixy) { - fhTriggerSelectionCorrelations->GetXaxis()->SetBinLabel(ixy + 1, TString::Format("#color[%d]{%s}", triggerSelectionFlags.test(ixy) ? 2 : 1, triggerSelectionExternalNamesMap.at(ixy).c_str()).Data()); - fhTriggerSelectionCorrelations->GetYaxis()->SetBinLabel(ixy + 1, TString::Format("#color[%d]{%s}", triggerSelectionFlags.test(ixy) ? 2 : 1, triggerSelectionExternalNamesMap.at(ixy).c_str()).Data()); + fhTriggerSelectionCorrelations->GetXaxis()->SetBinLabel(ixy + 1, TString::Format("#color[%d]{%s}", triggerSelectionFlags.test(ixy) ? 2 : 1, triggerSelectionExternalNamesMap.at(ixy).data()).Data()); + fhTriggerSelectionCorrelations->GetYaxis()->SetBinLabel(ixy + 1, TString::Format("#color[%d]{%s}", triggerSelectionFlags.test(ixy) ? 2 : 1, triggerSelectionExternalNamesMap.at(ixy).data()).Data()); } fhVertexZB = new TH1F("VertexZB", "Vertex Z; z_{vtx}", 60, -15, 15); @@ -530,32 +624,32 @@ struct DptDptFilter { #define DPTDPTCENTRALITYAXIS 105, -0.5f, 104.5f #define DPTDPTMULTIPLICITYAXIS(est) 1001, -0.5f, cfgEventSelection.multiplicityUpperLimit->getData()[fSystem][est] - 0.5f - std::string multestimator = getCentMultEstimatorName(fCentMultEstimator); + std::string_view multestimator = getCentMultEstimatorName(fCentMultEstimator); fhCentMultB = new TH1F("CentralityB", "Centrality before cut; centrality (%)", DPTDPTCENTRALITYAXIS); fhCentMultA = new TH1F("CentralityA", "Centrality; centrality (%)", DPTDPTCENTRALITYAXIS); - fhMultB = new TH1F("MultB", TString::Format("%s Multiplicity before cut;%s Multiplicity;Collisions", multestimator.c_str(), multestimator.c_str()), DPTDPTMULTIPLICITYAXIS(estimatorMultiplicitySourceMap.at(fCentMultEstimator))); - fhMultA = new TH1F("MultA", TString::Format("%s Multiplicity;%s Multiplicity;Collisions", multestimator.c_str(), multestimator.c_str()), DPTDPTMULTIPLICITYAXIS(estimatorMultiplicitySourceMap.at(fCentMultEstimator))); + fhMultB = new TH1F("MultB", TString::Format("%s Multiplicity before cut;%s Multiplicity;Collisions", multestimator.data(), multestimator.data()), DPTDPTMULTIPLICITYAXIS(estimatorMultiplicitySourceMap.at(fCentMultEstimator))); + fhMultA = new TH1F("MultA", TString::Format("%s Multiplicity;%s Multiplicity;Collisions", multestimator.data(), multestimator.data()), DPTDPTMULTIPLICITYAXIS(estimatorMultiplicitySourceMap.at(fCentMultEstimator))); if (cfgEventSelection.fillQc) { /* the quality control histograms */ for (int i = 0; i < BeforeAfterNOOFTIMES; ++i) { - fhMultiplicityVsCentrality[i] = new TH2F(TString::Format("MultiplicityVsCentrality%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;%s centrality (%%);Number of tracks", beforeAfterName[i].c_str(), multestimator.c_str()).Data(), DPTDPTCENTRALITYAXIS, DPTDPTMULTIPLICITYAXIS(MultSourceNtracks)); + fhMultiplicityVsCentrality[i] = new TH2F(TString::Format("MultiplicityVsCentrality%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;%s centrality (%%);Number of tracks", beforeAfterName[i].c_str(), multestimator.data()).Data(), DPTDPTCENTRALITYAXIS, DPTDPTMULTIPLICITYAXIS(MultSourceNtracks)); fhMultiplicityVsT0cMultiplicity[i] = new TH2F(TString::Format("MultiplicityVsT0cMultiplicity%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;T0C Multiplicity;Number of tracks", beforeAfterName[i].c_str()).Data(), DPTDPTMULTIPLICITYAXIS(MultSourceT0C), DPTDPTMULTIPLICITYAXIS(MultSourceNtracks)); fhMultiplicityVsT0aMultiplicity[i] = new TH2F(TString::Format("MultiplicityVsT0aMultiplicity%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;T0A Multiplicity;Number of tracks", beforeAfterName[i].c_str()).Data(), DPTDPTMULTIPLICITYAXIS(MultSourceT0A), DPTDPTMULTIPLICITYAXIS(MultSourceNtracks)); fhMultiplicityVsV0aMultiplicity[i] = new TH2F(TString::Format("MultiplicityVsV0aMultiplicity%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;V0A Multiplicity;Number of tracks", beforeAfterName[i].c_str()).Data(), DPTDPTMULTIPLICITYAXIS(MultSourceV0A), DPTDPTMULTIPLICITYAXIS(MultSourceNtracks)); fhMultiplicityVsPvMultiplicity[i] = new TH2F(TString::Format("MultiplicityVsPvMultiplicity%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;PV contributors;Number of tracks", beforeAfterName[i].c_str()).Data(), DPTDPTMULTIPLICITYAXIS(MultSourcePvContributors), DPTDPTMULTIPLICITYAXIS(MultSourceNtracks)); - fhPvMultiplicityVsCentrality[i] = new TH2F(TString::Format("PvMultiplicityVsCentrality%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;%s centrality (%%);PV contributors", beforeAfterName[i].c_str(), multestimator.c_str()).Data(), DPTDPTCENTRALITYAXIS, DPTDPTMULTIPLICITYAXIS(MultSourcePvContributors)); + fhPvMultiplicityVsCentrality[i] = new TH2F(TString::Format("PvMultiplicityVsCentrality%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;%s centrality (%%);PV contributors", beforeAfterName[i].c_str(), multestimator.data()).Data(), DPTDPTCENTRALITYAXIS, DPTDPTMULTIPLICITYAXIS(MultSourcePvContributors)); fhPvMultiplicityVsT0cMultiplicity[i] = new TH2F(TString::Format("PvMultiplicityVsT0cMultiplicity%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;T0C multiplicity;PV contributors", beforeAfterName[i].c_str()).Data(), DPTDPTMULTIPLICITYAXIS(MultSourceT0C), DPTDPTMULTIPLICITYAXIS(MultSourcePvContributors)); fhPvMultiplicityVsT0aMultiplicity[i] = new TH2F(TString::Format("PvMultiplicityVsT0aMultiplicity%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;T0A multiplicity;PV contributors", beforeAfterName[i].c_str()).Data(), DPTDPTMULTIPLICITYAXIS(MultSourceT0A), DPTDPTMULTIPLICITYAXIS(MultSourcePvContributors)); fhPvMultiplicityVsV0aMultiplicity[i] = new TH2F(TString::Format("PvMultiplicityVsV0aMultiplicity%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;V0A multiplicity;PV contributors", beforeAfterName[i].c_str()).Data(), DPTDPTMULTIPLICITYAXIS(MultSourceV0A), DPTDPTMULTIPLICITYAXIS(MultSourcePvContributors)); - fhV0aMultiplicityVsCentrality[i] = new TH2F(TString::Format("V0aMultiplicityVsCentrality%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;%s centrality (%%);V0A multiplicity", beforeAfterName[i].c_str(), multestimator.c_str()).Data(), DPTDPTCENTRALITYAXIS, DPTDPTMULTIPLICITYAXIS(MultSourceV0A)); + fhV0aMultiplicityVsCentrality[i] = new TH2F(TString::Format("V0aMultiplicityVsCentrality%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;%s centrality (%%);V0A multiplicity", beforeAfterName[i].c_str(), multestimator.data()).Data(), DPTDPTCENTRALITYAXIS, DPTDPTMULTIPLICITYAXIS(MultSourceV0A)); fhV0aMultiplicityVsT0cMultiplicity[i] = new TH2F(TString::Format("V0aMultiplicityVsT0cMultiplicity%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;T0C multiplicity;V0A multiplicity", beforeAfterName[i].c_str()).Data(), DPTDPTMULTIPLICITYAXIS(MultSourceT0C), DPTDPTMULTIPLICITYAXIS(MultSourceV0A)); fhV0aMultiplicityVsT0aMultiplicity[i] = new TH2F(TString::Format("V0aMultiplicityVsT0aMultiplicity%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;T0A multiplicity;V0A multiplicity", beforeAfterName[i].c_str()).Data(), DPTDPTMULTIPLICITYAXIS(MultSourceT0A), DPTDPTMULTIPLICITYAXIS(MultSourceV0A)); - fhT0cMultiplicityVsCentrality[i] = new TH2F(TString::Format("T0cMultiplicityVsCentrality%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;%s centrality (%%);T0C multiplicity", beforeAfterName[i].c_str(), multestimator.c_str()).Data(), DPTDPTCENTRALITYAXIS, DPTDPTMULTIPLICITYAXIS(MultSourceT0C)); + fhT0cMultiplicityVsCentrality[i] = new TH2F(TString::Format("T0cMultiplicityVsCentrality%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;%s centrality (%%);T0C multiplicity", beforeAfterName[i].c_str(), multestimator.data()).Data(), DPTDPTCENTRALITYAXIS, DPTDPTMULTIPLICITYAXIS(MultSourceT0C)); fhT0cMultiplicityVsT0aMultiplicity[i] = new TH2F(TString::Format("T0cMultiplicityVsT0aMultiplicity%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;T0A multiplicity;T0C multiplicity", beforeAfterName[i].c_str()).Data(), DPTDPTMULTIPLICITYAXIS(MultSourceT0A), DPTDPTMULTIPLICITYAXIS(MultSourceT0C)); - fhT0CentralityVsCentrality[i] = new TH2F(TString::Format("T0CentralityVsCentrality%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;%s centrality (%%);T0 centrality(%%)", beforeAfterName[i].c_str(), multestimator.c_str()).Data(), DPTDPTCENTRALITYAXIS, DPTDPTCENTRALITYAXIS); - fhV0aCentralityVsCentrality[i] = new TH2F(TString::Format("V0aCentralityVsCentrality%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;%s centrality (%%);V0A centrality (%%)", beforeAfterName[i].c_str(), multestimator.c_str()).Data(), DPTDPTCENTRALITYAXIS, DPTDPTCENTRALITYAXIS); - fhNtpvCentralityVsCentrality[i] = new TH2F(TString::Format("NtpvCentralityVsCentrality%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;%s centrality (%%);NTPV centrality (%%)", beforeAfterName[i].c_str(), multestimator.c_str()).Data(), DPTDPTCENTRALITYAXIS, DPTDPTCENTRALITYAXIS); + fhT0CentralityVsCentrality[i] = new TH2F(TString::Format("T0CentralityVsCentrality%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;%s centrality (%%);T0 centrality(%%)", beforeAfterName[i].c_str(), multestimator.data()).Data(), DPTDPTCENTRALITYAXIS, DPTDPTCENTRALITYAXIS); + fhV0aCentralityVsCentrality[i] = new TH2F(TString::Format("V0aCentralityVsCentrality%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;%s centrality (%%);V0A centrality (%%)", beforeAfterName[i].c_str(), multestimator.data()).Data(), DPTDPTCENTRALITYAXIS, DPTDPTCENTRALITYAXIS); + fhNtpvCentralityVsCentrality[i] = new TH2F(TString::Format("NtpvCentralityVsCentrality%s", beforeAfterSufix[i].c_str()).Data(), TString::Format("%s;%s centrality (%%);NTPV centrality (%%)", beforeAfterName[i].c_str(), multestimator.data()).Data(), DPTDPTCENTRALITYAXIS, DPTDPTCENTRALITYAXIS); } } @@ -693,6 +787,7 @@ void DptDptFilter::processReconstructed(CollisionObject const& collision, Tracks float centormult = tentativecentmult; int64_t orbit = collision.template bc_as().globalBC() / nBCsPerOrbit; bool withinOrbitOfInterest = (cfgEventSelection.minOrbit <= orbit) && (orbit < cfgEventSelection.maxOrbit); + storeMultiplicitiesAndCentralities(collision, ftracks); if (withinOrbitOfInterest && isEventSelected(collision, centormult)) { acceptedevent = true; fhCentMultA->Fill(centormult); diff --git a/PWGCF/TableProducer/dptDptFilter.h b/PWGCF/TableProducer/dptDptFilter.h index d99ae8951d8..d79924724f2 100644 --- a/PWGCF/TableProducer/dptDptFilter.h +++ b/PWGCF/TableProducer/dptDptFilter.h @@ -34,6 +34,7 @@ #include #include +#include #include #include #include @@ -47,6 +48,7 @@ #include #include #include +#include #include namespace o2 @@ -84,7 +86,7 @@ enum SystemType { /// \std::map systemInternalCodesMap /// \brief maps system names to internal system codes -static const std::map systemInternalCodesMap{ +static const std::map systemInternalCodesMap{ {"", SystemNoSystem}, {"pp", SystemPp}, {"pPb", SystemPPb}, @@ -99,7 +101,7 @@ static const std::map systemInternalCodesMap{ /// \std::map systemExternalNamesMap /// \brief maps system internal codes to system external names -static const std::map systemExternalNamesMap{ +static const std::map systemExternalNamesMap{ {SystemNoSystem, ""}, {SystemPp, "pp"}, {SystemPPb, "pPb"}, @@ -140,7 +142,7 @@ enum CentMultEstimatorType { /// \std::map estimatorInternalCodesMap /// \brief maps centrality/multiplicity estimator names to internal estimator codes -static const std::map estimatorInternalCodesMap{ +static const std::map estimatorInternalCodesMap{ {"NOCM", CentMultNOCM}, {"V0M", CentMultV0M}, {"CL0", CentMultCL0}, @@ -153,7 +155,7 @@ static const std::map estimatorInternalCodesMap{ /// \std::map estimatorExternalNamesMap /// \brief maps internal estimator codes to centrality/multiplicity estimator external names -static const std::map estimatorExternalNamesMap{ +static const std::map estimatorExternalNamesMap{ {CentMultNOCM, "NOCM"}, {CentMultV0M, "V0M"}, {CentMultCL0, "CL0"}, @@ -216,7 +218,7 @@ static const std::map estimatorMultiplicitySourceMap{ /// \std::vector multiplicitySourceExternalNamesMap /// \brief maps internal multiplicity source to external names for the LHC runs -static const std::vector> multiplicitySourceExternalNamesMap{ +static const std::vector> multiplicitySourceExternalNamesMap{ /* Run 1 and Run 2 */ { {MultSourceT0A, "T0A multiplicity"}, @@ -225,7 +227,7 @@ static const std::vector> multiplicitySourceExternalN {MultSourceV0A, "V0A multiplicity"}, {MultSourceV0C, "V0C multiplicity"}, {MultSourceV0M, "V0M multiplicity"}, - {MultSourceNtracks, "Number of tracks"}, + {MultSourceNtracks, "Global tracks"}, {MultSourcePvContributors, "PV contributors"}}, /* Run 3 */ { @@ -235,22 +237,94 @@ static const std::vector> multiplicitySourceExternalN {MultSourceV0A, "FV0A multiplicity"}, {MultSourceV0C, "WRONG SOURCE"}, {MultSourceV0M, "FV0M multiplicity"}, - {MultSourceNtracks, "Number of tracks"}, + {MultSourceNtracks, "Global tracks"}, {MultSourcePvContributors, "PV contributors"}}}; /// \std::map multiplicitySourceConfigNamesMap /// \brief maps internal multiplicity source to external configuration names /// At configuration time neither the system nor the lhc run is known -static const std::map multiplicitySourceConfigNamesMap{ - {MultSourceT0A, "FT0A"}, - {MultSourceT0C, "FT0C"}, - {MultSourceT0M, "FT0M"}, - {MultSourceV0A, "FV0A"}, - {MultSourceV0C, "V0C"}, - {MultSourceV0M, "FV0M"}, - {MultSourceNtracks, "Number of tracks"}, +static const std::map multiplicitySourceConfigNamesMap{ + {MultSourceT0A, "FT0A (T0A)"}, + {MultSourceT0C, "FT0C (T0C)"}, + {MultSourceT0M, "FT0M (T0M)"}, + {MultSourceV0A, "FV0A (V0A)"}, + {MultSourceV0C, "WRONG (V0C)"}, + {MultSourceV0M, "FV0M (V0M)"}, + {MultSourceNtracks, "Global tracks"}, {MultSourcePvContributors, "PV contributors"}}; +/// \enum CentMultCorrelationsParams +/// \brief internal codes for the supported parameters for centrality/multiplicity correlations exclusion cuts +enum CentMultCorrelationsParams { + CentMultCorrelationsMT0A = 0, ///< multiplicity from T0A + CentMultCorrelationsMT0C, ///< multiplicity from T0C + CentMultCorrelationsMT0M, ///< multiplicity from T0M + CentMultCorrelationsMV0A, ///< multiplicity from V0A + CentMultCorrelationsMV0C, ///< multiplicity from V0C (only Run 1 & Run 2) + CentMultCorrelationsMV0M, ///< multiplicity from V0M + CentMultCorrelationsMNGLTRK, ///< multiplicity from number of global tracks + CentMultCorrelationsMNPVC, ///< multiplicity from number of PV contributors + CentMultCorrelationsCT0A, ///< centrality from T0A + CentMultCorrelationsCT0C, ///< centrality from T0C + CentMultCorrelationsCT0M, ///< centrality from T0M + CentMultCorrelationsCV0A, ///< centrality from V0A + CentMultCorrelationsCNTPV, ///< centrality from number of PV contributors + CentMultCorrelationsNOOFPARAMS ///< the number of parameters supported +}; + +/// @brief prefix increment operator +/// @param ipar value +/// @return the incremented value +inline CentMultCorrelationsParams& operator++(CentMultCorrelationsParams& ipar) +{ + return ipar = static_cast(static_cast(ipar) + 1); +} + +/// @brief postfix increment operator +/// @param ipar the value +/// @param empty +/// @return the same value +inline CentMultCorrelationsParams operator++(CentMultCorrelationsParams& ipar, int) +{ + CentMultCorrelationsParams iparTmp(ipar); + ++ipar; + return iparTmp; +} + +/// \std::map centMultCorrelationsParamsMap +/// \brief maps centrality/multiplicity correlations parameters names to internal codes +static const std::map centMultCorrelationsParamsMap{ + {"MT0A", CentMultCorrelationsMT0A}, + {"MT0C", CentMultCorrelationsMT0C}, + {"MT0M", CentMultCorrelationsMT0M}, + {"MV0A", CentMultCorrelationsMV0A}, + {"MV0C", CentMultCorrelationsMV0C}, + {"MV0M", CentMultCorrelationsMV0M}, + {"MNGLTRK", CentMultCorrelationsMNGLTRK}, + {"MNPVC", CentMultCorrelationsMNPVC}, + {"CT0A", CentMultCorrelationsCT0A}, + {"CT0C", CentMultCorrelationsCT0C}, + {"CT0M", CentMultCorrelationsCT0M}, + {"CV0A", CentMultCorrelationsCV0A}, + {"CNTPV", CentMultCorrelationsCNTPV}}; + +/// \std::map centMultCorrelationsParamsNamesMap +/// \brief maps centrality/multiplicity correlations parameters internal codes to their external names +static const std::map centMultCorrelationsParamsNamesMap{ + {CentMultCorrelationsMT0A, "MT0A"}, + {CentMultCorrelationsMT0C, "MT0C"}, + {CentMultCorrelationsMT0M, "MT0M"}, + {CentMultCorrelationsMV0A, "MV0A"}, + {CentMultCorrelationsMV0C, "MV0C"}, + {CentMultCorrelationsMV0M, "MV0M"}, + {CentMultCorrelationsMNGLTRK, "MNGLTRK"}, + {CentMultCorrelationsMNPVC, "MNPVC"}, + {CentMultCorrelationsCT0A, "CT0A"}, + {CentMultCorrelationsCT0C, "CT0C"}, + {CentMultCorrelationsCT0M, "CT0M"}, + {CentMultCorrelationsCV0A, "CV0A"}, + {CentMultCorrelationsCNTPV, "CNTPV"}}; + /// \enum TriggerSelectionTags /// \brief The potential trigger tags to apply for event selection enum TriggerSelectionTags { @@ -275,7 +349,7 @@ enum TriggerSelectionTags { /// \std::map triggerSelectionBitsMap /// \brief maps trigger selection tags to internal trigger selection bits -static const std::map triggerSelectionBitsMap{ +static const std::map triggerSelectionBitsMap{ {"none", TriggSelNONE}, {"mb", TriggSelMB}, {"nosamebunchpup", TriggSelNOSAMEBUNCHPUP}, @@ -295,7 +369,7 @@ static const std::map triggerSelectionBitsMap{ /// \std::map triggerSelectionExternalNamesMap /// \brief maps trigger selection bits to external names -static const std::map triggerSelectionExternalNamesMap{ +static const std::map triggerSelectionExternalNamesMap{ {TriggSelNONE, "none"}, {TriggSelMB, "Sel8"}, ///< Sel8 includes kIsTriggerTVX, kNoTimeFrameBorder, and kNoITSROFrameBorder {TriggSelNOSAMEBUNCHPUP, "No same bunch pileup"}, @@ -325,22 +399,23 @@ enum OccupancyEstimationType { /// \enum CollisionSelectionFlags /// \brief The different criteria for selecting/rejecting collisions enum CollisionSelectionFlags { - CollSelIN = 0, ///< new unhandled, yet, event - CollSelMBBIT, ///< minimum bias - CollSelINT7BIT, ///< INT7 Run 1/2 - CollSelSEL7BIT, ///< Sel7 Run 1/2 - CollSelTRIGGSELBIT, ///< Accepted by trigger selection - CollSelRCTBIT, ///< Accetped by the RCT information - CollSelOCCUPANCYBIT, ///< occupancy within limits - CollSelCENTRALITYBIT, ///< centrality cut passed - CollSelZVERTEXBIT, ///< zvtx cut passed - CollSelSELECTED, ///< the event has passed all selections - CollSelNOOFFLAGS ///< number of flags + CollSelIN = 0, ///< new unhandled, yet, event + CollSelMBBIT, ///< minimum bias + CollSelINT7BIT, ///< INT7 Run 1/2 + CollSelSEL7BIT, ///< Sel7 Run 1/2 + CollSelTRIGGSELBIT, ///< Accepted by trigger selection + CollSelRCTBIT, ///< Accetped by the RCT information + CollSelOCCUPANCYBIT, ///< occupancy within limits + CollSelCENTRALITYBIT, ///< centrality cut passed + CollSelZVERTEXBIT, ///< zvtx cut passed + CollSelMULTCORRELATIONS, ///< multiplicities correlations passed + CollSelSELECTED, ///< the event has passed all selections + CollSelNOOFFLAGS ///< number of flags }; /// \std::mag collisionSelectionExternalNamesMap /// \brief maps collision selection bits to external names -static const std::map collisionSelectionExternalNamesMap{ +static const std::map collisionSelectionExternalNamesMap{ {CollSelIN, "In"}, {CollSelMBBIT, "MB"}, {CollSelINT7BIT, "INT7"}, @@ -350,6 +425,7 @@ static const std::map collisionSelectionExternalNamesMap{ {CollSelOCCUPANCYBIT, "Occupancy"}, {CollSelCENTRALITYBIT, "Centrality"}, {CollSelZVERTEXBIT, "z vertex"}, + {CollSelMULTCORRELATIONS, "Multiplicities correlations"}, {CollSelSELECTED, "Selected"}}; /// \enum StrongDebugging @@ -385,6 +461,14 @@ std::bitset<32> collisionSelectionFlags; std::bitset<32> triggerFlags; std::bitset<32> triggerSelectionFlags; +//============================================================================================ +// The collision exclusion using correlations between multiplicity/centrality observables +//============================================================================================ +bool useCentralityMultiplicityCorrelationsExclusion = false; +std::vector collisionMultiplicityCentralityObservables = {}; +std::vector observableIndexForCentralityMultiplicityParameter = {}; +TFormula* multiplicityCentralityCorrelationsExclusion = nullptr; + //============================================================================================ // The input data metadata access helper //============================================================================================ @@ -736,17 +820,17 @@ float particleMaxDCAxy = 999.9f; float particleMaxDCAZ = 999.9f; bool traceCollId0 = false; -inline std::bitset<32> getTriggerSelection(std::string const& triggstr) +inline std::bitset<32> getTriggerSelection(std::string_view const& triggstr) { std::bitset<32> flags; auto split = [](const auto s) { - std::vector tokens; - std::string token; + std::vector tokens; + std::string_view token; size_t posStart = 0; size_t posEnd; - while ((posEnd = s.find("+", posStart)) != std::string::npos) { + while ((posEnd = s.find("+", posStart)) != std::string_view::npos) { token = s.substr(posStart, posEnd - posStart); posStart = posEnd + 1; tokens.push_back(token); @@ -755,13 +839,13 @@ inline std::bitset<32> getTriggerSelection(std::string const& triggstr) return tokens; }; - std::vector tags = split(triggstr); + std::vector tags = split(triggstr); for (const auto& tag : tags) { if (triggerSelectionBitsMap.contains(tag)) { flags.set(triggerSelectionBitsMap.at(tag), true); } else { - LOGF(fatal, "Wrong trigger selection tag: %s", tag.c_str()); + LOGF(fatal, "Wrong trigger selection tag: %s", tag.data()); } } return flags; @@ -770,15 +854,16 @@ inline std::bitset<32> getTriggerSelection(std::string const& triggstr) inline SystemType getSytemTypeFromMetaData() { auto period = metadataInfo.get("LPMProductionTag"); + auto anchoredPeriod = metadataInfo.get("AnchorProduction"); - if (period == "LHC25ad" || period == "LHC25g5") { - LOGF(info, "Configuring for p-O LHC25ad period"); + if (period == "LHC25ad" || anchoredPeriod == "LHC25ad") { + LOGF(info, "Configuring for p-O (anchored to) LHC25ad period"); return SystemPORun3; - } else if (period == "LHC25ae" || period == "LHC25g6") { - LOGF(info, "Configuring for O-O LHC25ae period"); + } else if (period == "LHC25ae" || anchoredPeriod == "LHC25ae") { + LOGF(info, "Configuring for O-O (anchored to) LHC25ae period"); return SystemOORun3; - } else if (period == "LHC25af" || period == "LHC25g7") { - LOGF(info, "Configuring for Ne-Ne LHC25af period"); + } else if (period == "LHC25af" || anchoredPeriod == "LHC25af") { + LOGF(info, "Configuring for Ne-Ne (anchored to) LHC25af period"); return SystemNeNeRun3; } else { LOGF(fatal, "DptDptCorrelations::getSystemTypeFromMetadata(). No automatic system type configuration for %s period", period.c_str()); @@ -786,7 +871,7 @@ inline SystemType getSytemTypeFromMetaData() return SystemPbPb; } -inline SystemType getSystemType(std::string const& sysstr) +inline SystemType getSystemType(std::string_view const& sysstr) { /* we have to figure out how extract the system type */ if (sysstr == "Auto") { @@ -797,7 +882,7 @@ inline SystemType getSystemType(std::string const& sysstr) if (systemInternalCodesMap.contains(sysstr)) { return static_cast(systemInternalCodesMap.at(sysstr)); } else { - LOGF(fatal, "DptDptCorrelations::getSystemType(). Wrong system type: %s", sysstr.c_str()); + LOGF(fatal, "DptDptCorrelations::getSystemType(). Wrong system type: %s", sysstr.data()); } } return SystemPbPb; @@ -825,17 +910,17 @@ inline DataType getDataType(std::string const& datastr) return kData; } -inline CentMultEstimatorType getCentMultEstimator(std::string const& datastr) +inline CentMultEstimatorType getCentMultEstimator(std::string_view const& datastr) { if (estimatorInternalCodesMap.contains(datastr)) { return static_cast(estimatorInternalCodesMap.at(datastr)); } else { - LOGF(fatal, "Centrality/Multiplicity estimator %s not supported yet", datastr.c_str()); + LOGF(fatal, "Centrality/Multiplicity estimator %s not supported yet", datastr.data()); } return CentMultNOCM; } -inline std::string getCentMultEstimatorName(CentMultEstimatorType est) +inline std::string_view getCentMultEstimatorName(CentMultEstimatorType est) { if (estimatorExternalNamesMap.contains(est)) { return estimatorExternalNamesMap.at(est); @@ -845,7 +930,7 @@ inline std::string getCentMultEstimatorName(CentMultEstimatorType est) return "WRONG"; } -inline OccupancyEstimationType getOccupancyEstimator(const std::string& estimator) +inline OccupancyEstimationType getOccupancyEstimator(const std::string_view& estimator) { if (estimator == "None") { return OccupancyNOOCC; @@ -854,11 +939,38 @@ inline OccupancyEstimationType getOccupancyEstimator(const std::string& estimato } else if (estimator == "FT0C") { return OccupancyFT0COCC; } else { - LOGF(fatal, "Occupancy estimator %s not supported yet", estimator.c_str()); + LOGF(fatal, "Occupancy estimator %s not supported yet", estimator.data()); return OccupancyNOOCC; } } +/// @brief gets the exclusion formula from the corresponding expression and initializes the exclusion machinery +/// @param std::string_view with the formula expression +/// @return the expression TFormula +inline TFormula* getExclusionFormula(std::string_view formula) +{ + if (formula.length() != 0) { + useCentralityMultiplicityCorrelationsExclusion = true; + TFormula* f = new TFormula("Exclussion expression", formula.data()); + int nParameters = f->GetNpar(); + collisionMultiplicityCentralityObservables.resize(CentMultCorrelationsNOOFPARAMS); + observableIndexForCentralityMultiplicityParameter.resize(nParameters); + LOGF(info, "Configuring outliers exclusion with the formula %s which has %d parameters", formula.data(), nParameters); + for (int iPar = 0; iPar < nParameters; ++iPar) { + if (centMultCorrelationsParamsMap.contains(std::string(f->GetParName(iPar)))) { + observableIndexForCentralityMultiplicityParameter[iPar] = centMultCorrelationsParamsMap.at(std::string(f->GetParName(iPar))); + LOGF(info, "\tAssigned observable %s with index %d to the parameter %s with parameter index %d", centMultCorrelationsParamsNamesMap.at(centMultCorrelationsParamsMap.at(std::string(f->GetParName(iPar)))).data(), static_cast(centMultCorrelationsParamsMap.at(std::string(f->GetParName(iPar)))), f->GetParName(iPar), iPar); + } else { + LOGF(fatal, "Exclusion expression contains parameter %s which is still not supported. Please, fix it!", f->GetParName(iPar)); + } + } + return f; + } else { + useCentralityMultiplicityCorrelationsExclusion = false; + return nullptr; + } +} + ////////////////////////////////////////////////////////////////////////////////// /// Trigger selection ////////////////////////////////////////////////////////////////////////////////// @@ -1172,6 +1284,25 @@ inline bool centralitySelection(aod::McCollision const&, float } } +/// @brief evalues the exclusion formula for the current multiplicity / centrality parameters +/// @return true if the collision is not excluded according to the formula false otherwise +/// WARNING: it has always to be called after filling the multiplicity / centrality observables +inline bool isCollisionNotExcluded() +{ + bool notExcluded = true; + if (useCentralityMultiplicityCorrelationsExclusion) { + [&]() { + /* set the formula parameter values */ + for (size_t iPar = 0; iPar < observableIndexForCentralityMultiplicityParameter.size(); ++iPar) { + multiplicityCentralityCorrelationsExclusion->SetParameter(iPar, collisionMultiplicityCentralityObservables[observableIndexForCentralityMultiplicityParameter[iPar]]); + } + }(); + notExcluded = multiplicityCentralityCorrelationsExclusion->Eval() != 0; + } + collisionFlags.set(CollSelMULTCORRELATIONS, notExcluded); + return notExcluded; +} + ////////////////////////////////////////////////////////////////////////////////// /// Occupancy selection ////////////////////////////////////////////////////////////////////////////////// @@ -1308,7 +1439,9 @@ inline bool isEventSelected(CollisionObject const& collision, float& centormult) bool centmultsel = centralitySelection(collision, centormult); - bool accepted = trigsel && rctsel && occupancysel && zvtxsel && centmultsel; + bool centmultexclusion = isCollisionNotExcluded(); + + bool accepted = trigsel && rctsel && occupancysel && zvtxsel && centmultsel && centmultexclusion; if (accepted) { collisionFlags.set(CollSelSELECTED);