Outlier cut support for JCorran task

Author: jaelpark

PWGCF/JCorran/Tasks/jflucAnalysisTask.cxx

@@ -23,7 +23,11 @@
 #include "Framework/RunningWorkflowInfo.h"
 #include "ReconstructionDataFormats/V0.h"
 
+#include <TFormula.h>
+
+#include <array>
 #include <deque>
+#include <memory>
 
 // #include "CCDB/BasicCCDBManager.h"
 
@@ -56,6 +60,7 @@ struct jflucAnalysisTask {
   O2_DEFINE_CONFIGURABLE(ptmax, float, 5.0, "Maximum pt for tracks");
   O2_DEFINE_CONFIGURABLE(cfgCentBinsForMC, int, 0, "0 = OFF and 1 = ON for data like multiplicity/centrality bins for MC process");
   O2_DEFINE_CONFIGURABLE(cfgMultCorrelationsMask, uint16_t, 0, "Selection bitmask for the multiplicity correlations. This should match the filter selection cfgEstimatorBitMask.")
+  O2_DEFINE_CONFIGURABLE(cfgMultCutFormula, std::string, "", "Multiplicity correlations cut formula. A result greater than zero results in accepted event. Parameters: [cFT0C] FT0C centrality, [mFV0A] V0A multiplicity, [mGlob] global track multiplicity, [mPV] PV track multiplicity")
 
   ConfigurableAxis axisMultiplicity{"axisMultiplicity", {VARIABLE_WIDTH, 0, 5, 10, 20, 30, 40, 50, 100.1}, "multiplicity / centrality axis for histograms"};
   ConfigurableAxis phiAxis{"axisPhi", {50, 0.0, o2::constants::math::TwoPI}, "phi axis for histograms"};
@@ -71,6 +76,9 @@ struct jflucAnalysisTask {
 
   HistogramRegistry registry{"registry"};
 
+  std::unique_ptr<TFormula> multCutFormula;
+  std::array<uint, 4> multCutFormulaParamIndex;
+
   void init(InitContext const&)
   {
     auto axisSpecMult = AxisSpec(axisMultiplicity);
@@ -98,6 +106,28 @@ struct jflucAnalysisTask {
     }
     if ((doprocessCFDerivedMultSet || doprocessCFDerivedMultSetCorrected) && cfgMultCorrelationsMask == 0)
       LOGF(fatal, "cfgMultCorrelationsMask can not be 0 when MultSet process functions are in use.");
+
+    if (!cfgMultCutFormula.value.empty()) {
+      if (cfgMultCorrelationsMask == 0)
+        LOGF(fatal, "cfgMultCorrelationsMask can not be 0 when outlier cuts are enabled.");
+      multCutFormula = std::make_unique<TFormula>("multCutFormula", cfgMultCutFormula.value.c_str());
+      std::fill_n(multCutFormulaParamIndex.begin(), std::size(multCutFormulaParamIndex), ~0u);
+      std::array<std::string, 4> pars = {"cFT0C", "mFV0A", "mPV", "mGlob"}; // must correspond the order of MultiplicityEstimators
+      for (uint i = 0, n = multCutFormula->GetNpar(); i < n; ++i) {
+        auto m = std::find(pars.begin(), pars.end(), multCutFormula->GetParName(i));
+        if (m == pars.end()) {
+
+          LOGF(warning, "Unknown parameter in cfgMultCutFormula: %s", multCutFormula->GetParName(i));
+          continue;
+        }
+        if ((cfgMultCorrelationsMask.value & (1u << i)) == 0) {
+          LOGF(warning, "The centrality/multiplicity estimator %s is not available to be used in cfgMultCutFormula. Ensure cfgMultCorrelationsMask is correct and matches the CFMultSets in derived data.");
+        } else {
+          multCutFormulaParamIndex[std::distance(pars.begin(), m)] = i;
+          LOGF(info, "Multiplicity cut parameter %s in use.", m->c_str());
+        }
+      }
+    }
   }
 
   template <class T>
@@ -138,6 +168,20 @@ struct jflucAnalysisTask {
     }
   }
 
+  template <class CollType>
+  bool passOutlier(CollType const& collision)
+  {
+    if (cfgMultCutFormula.value.empty())
+      return true;
+    for (uint i = 0; i < 4; ++i) {
+      if ((cfgMultCorrelationsMask.value & (1u << i)) == 0 || multCutFormulaParamIndex[i] == ~0u)
+        continue;
+      auto estIndex = std::popcount(cfgMultCorrelationsMask.value & ((1u << i) - 1));
+      multCutFormula->SetParameter(multCutFormulaParamIndex[i], collision.multiplicities()[estIndex]);
+    }
+    return multCutFormula->Eval() > 0.0f;
+  }
+
   void processJDerived(aod::JCollision const& collision, soa::Filtered<aod::JTracks> const& tracks)
   {
     analyze(collision, tracks);
@@ -166,6 +210,8 @@ struct jflucAnalysisTask {
   {
     if (std::popcount(cfgMultCorrelationsMask.value) != static_cast<int>(collision.multiplicities().size()))
       LOGF(fatal, "Multiplicity selections (cfgMultCorrelationsMask = 0x%x) do not match the size of the table column (%ld). The histogram filling relies on the preservation of order.", cfgMultCorrelationsMask.value, collision.multiplicities().size());
+    if (!passOutlier(collision))
+      return;
     analyze(collision, tracks);
   }
   PROCESS_SWITCH(jflucAnalysisTask, processCFDerivedMultSet, "Process CF derived data with multiplicity sets", false);
@@ -174,6 +220,8 @@ struct jflucAnalysisTask {
   {
     if (std::popcount(cfgMultCorrelationsMask.value) != static_cast<int>(collision.multiplicities().size()))
       LOGF(fatal, "Multiplicity selections (cfgMultCorrelationsMask = 0x%x) do not match the size of the table column (%ld). The histogram filling relies on the preservation of order.", cfgMultCorrelationsMask.value, collision.multiplicities().size());
+    if (!passOutlier(collision))
+      return;
     analyze(collision, tracks);
   }
   PROCESS_SWITCH(jflucAnalysisTask, processCFDerivedMultSetCorrected, "Process CF derived data with corrections and multiplicity sets", false);