[PWGCF] Outlier cut support for JCorran task

PWGCF/JCorran/Tasks/jflucAnalysisTask.cxx

@@ -1,4 +1,4 @@
 // Copyright 2019-2020 CERN and copyright holders of ALICE O2.
 // See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
 // All rights not expressly granted are reserved.
 //
@@ -10,7 +10,7 @@
 // or submit itself to any jurisdiction.
 /// \author Jasper Parkkila (jparkkil@cern.ch)
 /// \author Dong Jo Kim (djkim@jyu.fi)
 /// \since Sep 2022
 
 #include "Common/Core/TrackSelection.h"
 #include "Common/DataModel/Centrality.h"
@@ -23,7 +23,12 @@
 #include "Framework/RunningWorkflowInfo.h"
 #include "ReconstructionDataFormats/V0.h"
 
+#include <TFormula.h>
+
+#include <array>
 #include <deque>
+#include <memory>
+#include <string>
 
 // #include "CCDB/BasicCCDBManager.h"
 
@@ -41,7 +46,7 @@
 
 #define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable<TYPE> NAME{#NAME, DEFAULT, HELP};
 
 struct jflucAnalysisTask {
   ~jflucAnalysisTask()
   {
     if (pcf)
@@ -56,10 +61,11 @@
   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"};
   ConfigurableAxis etaAxis{"axisEta", {40, -2.0, 2.0}, "eta axis for histograms"};
   ConfigurableAxis zvtAxis{"axisZvt", {20, -10.0, 10.0}, "zvertex axis for histograms"};
   ConfigurableAxis ptAxis{"axisPt", {60, 0.0, 300.0}, "pt axis for histograms"};
   ConfigurableAxis massAxis{"axisMass", {1, 0.0, 10.0}, "mass axis for histograms"};
@@ -71,6 +77,9 @@
 
   HistogramRegistry registry{"registry"};
 
+  std::unique_ptr<TFormula> multCutFormula;
+  std::array<uint, 4> multCutFormulaParamIndex;
+
   void init(InitContext const&)
   {
     auto axisSpecMult = AxisSpec(axisMultiplicity);
@@ -98,10 +107,32 @@
     }
     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>
   using hasInvMass = decltype(std::declval<T&>().invMass());
 
   template <class CollisionT, class TrackT>
   void analyze(CollisionT const& collision, TrackT const& tracks)
@@ -120,7 +151,7 @@
   void analyze(CollisionT const& collision, POITrackT const& poiTracks, REFTrackT const& refTracks)
   {
     if constexpr (std::experimental::is_detected<hasInvMass, typename POITrackT::iterator>::value) {
       for (auto& track : poiTracks)
         registry.fill(HIST("invMass"), track.invMass(), track.pt(), collision.multiplicity());
     }
     pcf2Prong->Init();
@@ -138,6 +169,20 @@
     }
   }
 
+  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 +211,8 @@
   {
     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 +221,8 @@
   {
     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);