Merge branch 'master' of https://github.com/AliceO2Group/O2Physics into updateworkflow

Author: Jinhyun Park

Common/CCDB/ctpRateFetcher.cxx

@@ -95,9 +95,7 @@ double ctpRateFetcher::pileUpCorrection(double triggerRate)
   if (mLHCIFdata == nullptr) {
     LOG(fatal) << "No filling" << std::endl;
   }
-  auto bfilling = mLHCIFdata->getBunchFilling();
-  std::vector<int> bcs = bfilling.getFilledBCs();
-  double nbc = bcs.size();
+  double nbc = mLHCIFdata->getBunchFilling().getPattern().count();
   double nTriggersPerFilledBC = triggerRate / nbc / constants::lhc::LHCRevFreq;
   double mu = -std::log(1 - nTriggersPerFilledBC);
   return mu * nbc * constants::lhc::LHCRevFreq;

PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx

@@ -21,19 +21,30 @@
 #include "Common/DataModel/PIDResponseTPC.h"
 #include "Common/DataModel/TrackSelectionTables.h"
 
-#include "CCDB/BasicCCDBManager.h"
-#include "Framework/ASoAHelpers.h"
-#include "Framework/AnalysisDataModel.h"
-#include "Framework/AnalysisTask.h"
-#include "Framework/HistogramRegistry.h"
-#include "Framework/HistogramSpec.h"
-#include "Framework/O2DatabasePDGPlugin.h"
-#include "Framework/runDataProcessing.h"
-
+#include <CCDB/BasicCCDBManager.h>
+#include <CommonConstants/PhysicsConstants.h>
+#include <Framework/ASoAHelpers.h>
+#include <Framework/AnalysisDataModel.h>
+#include <Framework/AnalysisHelpers.h>
+#include <Framework/AnalysisTask.h>
+#include <Framework/Configurable.h>
+#include <Framework/HistogramRegistry.h>
+#include <Framework/HistogramSpec.h>
+#include <Framework/InitContext.h>
+#include <Framework/O2DatabasePDGPlugin.h>
+#include <Framework/OutputObjHeader.h>
+#include <Framework/runDataProcessing.h>
+
+#include <TH2.h>
 #include <TPDGCode.h>
 
 #include <algorithm>
+#include <array>
+#include <cstddef>
+#include <cstdint>
+#include <cstdlib>
 #include <string>
+#include <string_view>
 #include <vector>
 
 using namespace o2;
@@ -112,6 +123,17 @@ struct NchCumulantsId {
   Configurable<float> cfgCutPtMax{"cfgCutPtMax", 3.0, "max cut for pT"};
   Configurable<float> cfgCutPtMin{"cfgCutPtMin", 0.15, "min cut for pT"};
 
+  ConfigurableAxis axisNch{"axisNch", {1200, -60.0, 60.0}, "Net_charge_dN"};
+  ConfigurableAxis axisPosCh{"axisPosCh", {3010, -0.5, 300.5}, "Pos_charge"};
+  ConfigurableAxis axisNegCh{"axisNegCh", {3010, -0.5, 300.5}, "Neg_charge"};
+  ConfigurableAxis axisNt{"axisNt", {8010, -0.5, 800.5}, "Mult_midRap_Nch"};
+  ConfigurableAxis axisPrCh{"axisPrCh", {3010, -0.5, 300.5}, "Pr_charge"};
+  ConfigurableAxis axisAPrCh{"axisAPrCh", {3010, -0.5, 300.5}, "APr_charge"};
+  ConfigurableAxis axisKaCh{"axisKaCh", {3010, -0.5, 300.5}, "Ka_charge"};
+  ConfigurableAxis axisAKaCh{"axisAKaCh", {3010, -0.5, 300.5}, "AKa_charge"};
+  ConfigurableAxis axisPiCh{"axisPiCh", {3010, -0.5, 300.5}, "Pion_Positive"};
+  ConfigurableAxis axisAPiCh{"axisAPiCh", {3010, -0.5, 300.5}, "Pion_Negative"};
+
   Configurable<bool> checkCollPosZMc{"checkCollPosZMc", false, "checkCollPosZMc"};
   Configurable<bool> flagUnusedVariableError{"flagUnusedVariableError", false, "flagUnusedVariableError"};
   Configurable<bool> cfgDoRejectionForId{"cfgDoRejectionForId", false, "Apply rejection cut before PID selection (selTrackForId)"};
@@ -121,7 +143,7 @@ struct NchCumulantsId {
   Configurable<bool> cfgEvSel03doIsGoodITSLayersAll{"cfgEvSel03doIsGoodITSLayersAll", true, "apply kIsGoodITSLayersAll"};
 
   // Configurables for particle Identification
-  Configurable<bool> cfgId01CheckVetoCut{"cfgId01CheckVetoCut", false, "cfgId01CheckVetoCut"};
+  Configurable<bool> cfgId01CheckVetoCut{"cfgId01CheckVetoCut", true, "cfgId01CheckVetoCut"};
   Configurable<bool> cfgId02DoElRejection{"cfgId02DoElRejection", true, "cfgId02DoElRejection"};
   Configurable<bool> cfgId03DoDeRejection{"cfgId03DoDeRejection", false, "cfgId03DoDeRejection"};
   Configurable<bool> cfgId04DoPdependentId{"cfgId04DoPdependentId", true, "cfgId04DoPdependentId"};
@@ -246,17 +268,6 @@ struct NchCumulantsId {
     const AxisSpec axisTOFNSigma = {200, -10.0, 10.0, "n#sigma_{TOF}"};
     const AxisSpec axisTOFExpMom = {200, 0.0f, 10.0f, "#it{p}_{tofExpMom} (GeV/#it{c})"};
 
-    const AxisSpec axisNch(100, -50, 50, "Net_charge_dN");
-    const AxisSpec axisPosCh(101, -1, 100, "Pos_charge");
-    const AxisSpec axisNegCh(101, -1, 100, "Neg_charge");
-    const AxisSpec axisNt(201, -1, 200, "Mult_midRap_Nch");
-    const AxisSpec axisPrCh(101, -1, 100, "Pr_charge");
-    const AxisSpec axisAPrCh(101, -1, 100, "APr_charge");
-    const AxisSpec axisKaCh(101, -1, 100, "Ka_charge");
-    const AxisSpec axisAKaCh(101, -1, 100, "AKa_charge");
-    const AxisSpec axisPiCh(101, -1, 100, "Pion_Positive");
-    const AxisSpec axisAPiCh(101, -1, 100, "Pion_Negative");
-
     const AxisSpec axisIdTag = {32, -0.5f, 31.5f, "idTag"};
     const AxisSpec axisMcTag = {32, -0.5f, 31.5f, "mcTag"};
 

PWGCF/EbyEFluctuations/Tasks/radialFlowDecorr.cxx

@@ -76,7 +76,7 @@ struct RadialFlowDecorr {
   static constexpr int KnFt0cCell = 96;
   static constexpr int KIntM = 3;
   static constexpr int KIntK = 3;
-  static constexpr int KNEta = 17;
+  static constexpr int KNEta = 9;
   static constexpr float KFloatEpsilon = 1e-6f;
   static constexpr int KPiPlus = 211;
   static constexpr int KKPlus = 321;
@@ -151,10 +151,10 @@ struct RadialFlowDecorr {
   static constexpr float KinvalidCentrality = -1.0f;
   inline static const std::vector<float> etaLw = {
     -0.8,
-    -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7};
+    -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6};
   inline static const std::vector<float> etaUp = {
     0.8,
-    -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8};
+    -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8};
 
   Configurable<float> cfgVtxZCut{"cfgVtxZCut", 10.f, "z-vertex range"};
   Configurable<float> cfgPtMin{"cfgPtMin", 0.2f, "min pT"};
@@ -187,6 +187,7 @@ struct RadialFlowDecorr {
   Configurable<float> cfgCutPtUpper{"cfgCutPtUpper", 10.0f, "Higher pT cut for inclusive hadron analysis"};
   Configurable<float> cfgCutPtUpperPID{"cfgCutPtUpperPID", 6.0f, "Higher pT cut for identified particle analysis"};
   Configurable<float> cfgCutEta{"cfgCutEta", 0.8f, "absolute Eta cut"};
+  Configurable<int> cfgMinTracksPerEtaBin{"cfgMinTracksPerEtaBin", 0, "Min weighted-track sum required in every narrow eta bin for inclusive species (0 = disabled)"};
   Configurable<int> cfgNsubsample{"cfgNsubsample", 10, "Number of subsamples"};
   Configurable<int> cfgCentralityChoice{"cfgCentralityChoice", 1, "Which centrality estimator? 1-->FT0C, 2-->FT0M, 3-->FDDM, 4-->FV0A"};
   Configurable<bool> cfgEvSelNoSameBunchPileup{"cfgEvSelNoSameBunchPileup", true, "Pileup removal"};
@@ -201,7 +202,7 @@ struct RadialFlowDecorr {
   Configurable<float> cfgLinPupParam2{"cfgLinPupParam2", 3.0f, "(Lower) Linear Pileup Cut Const"};
   Configurable<float> cfgLinPupParam3{"cfgLinPupParam3", 3.0f, "(Lower) Linear Pileup Slope"};
 
-  Configurable<int> cfgNchPbMax{"cfgNchPbMax", 4000, "Max Nch range for PbPb collisions"};
+  Configurable<int> cfgNchPbMax{"cfgNchPbMax", 5000, "Max Nch range for PbPb collisions"};
   Configurable<int> cfgNchOMax{"cfgNchOMax", 800, "Max Nch range for OO collisions"};
 
   Configurable<int> cfgSys{"cfgSys", 1, "Efficiency to be used for which system? 1-->PbPb, 2-->NeNe, 3-->OO, 4-->pp"};
@@ -222,21 +223,17 @@ struct RadialFlowDecorr {
   const AxisSpec vzAxis{5, -12.5, 12.5, "Vz"};
   const AxisSpec chgAxis{3, -1.5, 1.5};
   const AxisSpec pTAxis{{0.0, 0.2, 0.4, 0.6, 0.8, 1, 3, 5, 7, 10}, "pT Axis"};
-  const AxisSpec etaAxis{{-0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9}, "Eta"};
+  const AxisSpec etaAxis{{-0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8}, "Eta"};
   const AxisSpec phiAxis{KNbinsPhi, KPhiMin, TwoPI, "#phi"};
   const AxisSpec etaBinAxis{KNEta + 1, -0.5, KNEta + 0.5, "#eta bin Number"};
   const AxisSpec spBinAxis{KNsp + 1, -KBinOffset, static_cast<float>(KNsp) + KBinOffset, "species index Number"};
 
-  const AxisSpec gapAxis{{-1.55, -1.45, -1.35, -1.25, -1.15, -1.05, -0.95, -0.85,
-                          -0.75, -0.65, -0.55, -0.45, -0.35, -0.25, -0.15, -0.05,
-                          0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75,
-                          0.85, 0.95, 1.05, 1.15, 1.25, 1.35, 1.45, 1.55},
+  const AxisSpec gapAxis{{-1.5, -1.3, -1.1, -0.9, -0.7, -0.5, -0.3, -0.1,
+                          0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5},
                          "Gap"};
 
-  const AxisSpec sumAxis{{-1.55, -1.45, -1.35, -1.25, -1.15, -1.05, -0.95, -0.85,
-                          -0.75, -0.65, -0.55, -0.45, -0.35, -0.25, -0.15, -0.05,
-                          0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75,
-                          0.85, 0.95, 1.05, 1.15, 1.25, 1.35, 1.45, 1.55},
+  const AxisSpec sumAxis{{-1.5, -1.3, -1.1, -0.9, -0.7, -0.5, -0.3, -0.1,
+                          0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5},
                          "Sum"};
 
   Configurable<bool> cfgRunMCGetNSig{"cfgRunMCGetNSig", false, "Run MC pass to get mean of Nsig Plots"};
@@ -552,6 +549,32 @@ struct RadialFlowDecorr {
     return true;
   }
 
+  template <std::size_t NspT, std::size_t NetaT, std::size_t NkT>
+  bool hasMinTracksInAllEtaBins(const double (&sw)[NspT][NetaT][NkT])
+  {
+    const int minTracks = cfgMinTracksPerEtaBin;
+    if (minTracks <= 0)
+      return true;
+    for (std::size_t ieta = 1; ieta < NetaT; ++ieta) {
+      if (sw[kInclusiveIdx][ieta][1] < static_cast<double>(minTracks))
+        return false;
+    }
+    return true;
+  }
+
+  template <std::size_t NspT, std::size_t NetaT>
+  bool hasMinTracksInAllEtaBins(const double (&sw)[NspT][NetaT])
+  {
+    const int minTracks = cfgMinTracksPerEtaBin;
+    if (minTracks <= 0)
+      return true;
+    for (std::size_t ieta = 1; ieta < NetaT; ++ieta) {
+      if (sw[kInclusiveIdx][ieta] < static_cast<double>(minTracks))
+        return false;
+    }
+    return true;
+  }
+
   template <typename T>
   bool isTrackSelected(const T& trk)
   {
@@ -1996,6 +2019,9 @@ struct RadialFlowDecorr {
       }
     }
 
+    if (!hasMinTracksInAllEtaBins(sumWiTruth) || !hasMinTracksInAllEtaBins(sumWiReco))
+      return;
+
     for (int isp = 0; isp < KNsp; ++isp) {
       histos.fill(HIST("MCReco/Prof_Cent_Nsp_Nchrec"), cent, isp, sumWiReco[isp][0]);
       histos.fill(HIST("MCReco/Prof_Mult_Nsp_Nchrec"), multPV, isp, sumWiReco[isp][0]);
@@ -2352,6 +2378,9 @@ struct RadialFlowDecorr {
       }
     } // trkslice
 
+    if (!hasMinTracksInAllEtaBins(sumWkTru) || !hasMinTracksInAllEtaBins(sumWkReco))
+      return;
+
     for (int ieta = 0; ieta < KNEta; ++ieta) {
       const int ibx = state.pmeanTruNchEtabinSpbinStep2->GetXaxis()->FindBin(mcCollision.multNTracksPV());
       const int iby = ieta + 1;
@@ -2546,9 +2575,12 @@ struct RadialFlowDecorr {
         float sum = (etaValA + etaValB);
         for (int isp = 0; isp < KNsp; ++isp) {
 
-          float c2SubTru = p1kBarTru[isp][ietaA] * p1kBarTru[isp][ietaC];
-          float c2SubReco = p1kBarReco[isp][ietaA] * p1kBarReco[isp][ietaC];
-          float c2SubRecoEffCor = p1kBarRecoEffCor[isp][ietaA] * p1kBarRecoEffCor[isp][ietaC];
+          float c2SubTru = (ietaA == ietaC) ? static_cast<float>(c2Tru[isp][ietaA])
+                                            : p1kBarTru[isp][ietaA] * p1kBarTru[isp][ietaC];
+          float c2SubReco = (ietaA == ietaC) ? static_cast<float>(c2Reco[isp][ietaA])
+                                             : p1kBarReco[isp][ietaA] * p1kBarReco[isp][ietaC];
+          float c2SubRecoEffCor = (ietaA == ietaC) ? static_cast<float>(c2RecoEffCor[isp][ietaA])
+                                                   : p1kBarRecoEffCor[isp][ietaA] * p1kBarRecoEffCor[isp][ietaC];
 
           float covTru = p1kBarTruMult[isp][ietaA] * p1kBarTru[isp][ietaC];
           float covReco = p1kBarRecoMult[isp][ietaA] * p1kBarReco[isp][ietaC];
@@ -3200,6 +3232,9 @@ struct RadialFlowDecorr {
       }
     }
 
+    if (!hasMinTracksInAllEtaBins(sumWi))
+      return;
+
     for (int isp = 0; isp < KNsp; ++isp) {
       if (sumWi[isp][0] < 1.0f)
         continue;
@@ -3373,6 +3408,9 @@ struct RadialFlowDecorr {
       }
     }
 
+    if (!hasMinTracksInAllEtaBins(sumwk))
+      return;
+
     double amplFT0A = 0, amplFT0C = 0;
     if (coll.has_foundFT0()) {
       const auto& ft0 = coll.foundFT0();
@@ -3440,8 +3478,8 @@ struct RadialFlowDecorr {
           histos.fill(HIST("Prof_Cov_Mult_etabin_spbin"), coll.multNTracksPV(), ietaA, isp, covAC);
         }
         if (std::isfinite(covCA)) {
-          histos.fill(HIST("Prof_Cov_Cent_etabin_spbin"), cent, ietaA, isp, covCA);
-          histos.fill(HIST("Prof_Cov_Mult_etabin_spbin"), coll.multNTracksPV(), ietaA, isp, covCA);
+          histos.fill(HIST("Prof_Cov_Cent_etabin_spbin"), cent, ietaC, isp, covCA);
+          histos.fill(HIST("Prof_Cov_Mult_etabin_spbin"), coll.multNTracksPV(), ietaC, isp, covCA);
         }
         if (std::isfinite(covFT0A)) {
           histos.fill(HIST("Prof_CovFT0A_Cent_etabin_spbin"), cent, ietaA, isp, covFT0A);
@@ -3464,7 +3502,8 @@ struct RadialFlowDecorr {
 
         for (int isp = 0; isp < KNsp; ++isp) {
 
-          float c2Sub = p1kBar[isp][ietaA] * p1kBar[isp][ietaC];
+          float c2Sub = (ietaA == ietaC) ? static_cast<float>(c2[isp][ietaA])
+                                         : p1kBar[isp][ietaA] * p1kBar[isp][ietaC];
           float cov = p1kBarMult[isp][ietaA] * p1kBar[isp][ietaC];
           float covFT0A = p1kBarFt0A * p1kBar[isp][ietaC];
           float covFT0C = p1kBarFt0C * p1kBar[isp][ietaA];