[PWGCF] flow: limit centrality in filter; add group NUA (#9081)

Author: Luzhiyongg

PWGCF/Flow/Tasks/flowRunbyRun.cxx

@@ -63,6 +63,7 @@ struct FlowRunbyRun {
   O2_DEFINE_CONFIGURABLE(cfgCutDCAz, float, 2.0f, "max DCA to vertex z")
   O2_DEFINE_CONFIGURABLE(cfgUseNch, bool, false, "Use Nch for flow observables")
   O2_DEFINE_CONFIGURABLE(cfgOutputNUAWeightsRefPt, bool, false, "NUA weights are filled in ref pt bins")
+  O2_DEFINE_CONFIGURABLE(cfgDynamicRunNumber, bool, false, "Add runNumber during runtime")
   Configurable<std::vector<int>> cfgRunNumbers{"cfgRunNumbers", std::vector<int>{544095, 544098, 544116, 544121, 544122, 544123, 544124}, "Preconfigured run numbers"};
   Configurable<std::vector<std::string>> cfgUserDefineGFWCorr{"cfgUserDefineGFWCorr", std::vector<std::string>{"refN10 {2} refP10 {-2}"}, "User defined GFW CorrelatorConfig"};
   Configurable<std::vector<std::string>> cfgUserDefineGFWName{"cfgUserDefineGFWName", std::vector<std::string>{"Ch10Gap22"}, "User defined GFW Name"};
@@ -83,7 +84,7 @@ struct FlowRunbyRun {
 
   // Define output
   OutputObj<FlowContainer> fFC{FlowContainer("FlowContainer")};
-  OutputObj<TList> fWeightList{"WeightList", OutputObjHandlingPolicy::AnalysisObject};
+  OutputObj<TList> fWeightList{"WeightList", OutputObjHandlingPolicy::AnalysisObject, OutputObjSourceType::OutputObjSource};
   HistogramRegistry registry{"registry"};
 
   // define global variables
@@ -248,7 +249,7 @@ struct FlowRunbyRun {
     float lRandom = fRndm->Rndm();
     if (runNumber != lastRunNumer) {
       lastRunNumer = runNumber;
-      if (std::find(runNumbers.begin(), runNumbers.end(), runNumber) == runNumbers.end()) {
+      if (cfgDynamicRunNumber && std::find(runNumbers.begin(), runNumbers.end(), runNumber) == runNumbers.end()) {
         // if run number is not in the preconfigured list, create new output histograms for this run
         createOutputObjectsForRun(runNumber);
         runNumbers.push_back(runNumber);

PWGCF/Flow/Tasks/flowTask.cxx

@@ -54,6 +54,8 @@ using namespace o2::framework::expressions;
 struct FlowTask {
 
   O2_DEFINE_CONFIGURABLE(cfgCutVertex, float, 10.0f, "Accepted z-vertex range")
+  O2_DEFINE_CONFIGURABLE(cfgCentMin, float, 0.0f, "Minimum centrality")
+  O2_DEFINE_CONFIGURABLE(cfgCentMax, float, 100.0f, "Maximum centrality")
   O2_DEFINE_CONFIGURABLE(cfgCutPtPOIMin, float, 0.2f, "Minimal pT for poi tracks")
   O2_DEFINE_CONFIGURABLE(cfgCutPtPOIMax, float, 10.0f, "Maximal pT for poi tracks")
   O2_DEFINE_CONFIGURABLE(cfgCutPtRefMin, float, 0.2f, "Minimal pT for ref tracks")
@@ -84,8 +86,11 @@ struct FlowTask {
   O2_DEFINE_CONFIGURABLE(cfgNbootstrap, int, 10, "Number of subsamples")
   O2_DEFINE_CONFIGURABLE(cfgOutputNUAWeights, bool, false, "Fill and output NUA weights")
   O2_DEFINE_CONFIGURABLE(cfgOutputNUAWeightsRefPt, bool, false, "NUA weights are filled in ref pt bins")
+  O2_DEFINE_CONFIGURABLE(cfgOutputGroupNUAWeights, bool, false, "Fill and output group NUA weights")
   O2_DEFINE_CONFIGURABLE(cfgEfficiency, std::string, "", "CCDB path to efficiency object")
   O2_DEFINE_CONFIGURABLE(cfgAcceptance, std::string, "", "CCDB path to acceptance object")
+  O2_DEFINE_CONFIGURABLE(cfgAcceptanceGroup, std::string, "", "CCDB path to group acceptance object")
+  O2_DEFINE_CONFIGURABLE(cfgAcceptanceGroupUse, bool, false, "Apply group acceptance, this option overrides cfgAcceptance")
   O2_DEFINE_CONFIGURABLE(cfgMagnetField, std::string, "GLO/Config/GRPMagField", "CCDB path to Magnet field object")
   O2_DEFINE_CONFIGURABLE(cfgEvSelOccupancy, bool, true, "Occupancy cut")
   O2_DEFINE_CONFIGURABLE(cfgCutOccupancyHigh, int, 500, "High cut on TPC occupancy")
@@ -94,6 +99,7 @@ struct FlowTask {
   Configurable<std::vector<std::string>> cfgUserDefineGFWCorr{"cfgUserDefineGFWCorr", std::vector<std::string>{"refN02 {2} refP02 {-2}", "refN12 {2} refP12 {-2}"}, "User defined GFW CorrelatorConfig"};
   Configurable<std::vector<std::string>> cfgUserDefineGFWName{"cfgUserDefineGFWName", std::vector<std::string>{"Ch02Gap22", "Ch12Gap22"}, "User defined GFW Name"};
   Configurable<std::vector<int>> cfgRunRemoveList{"cfgRunRemoveList", std::vector<int>{-1}, "excluded run numbers"};
+  Configurable<std::vector<int>> cfgGroupSplitRunNumber{"cfgGroupSplitRunNumber", std::vector<int>{544510, 544653}, "runnumbers for group splitting (suppose run numbers are increasing monotonically) "};
 
   ConfigurableAxis axisVertex{"axisVertex", {40, -20, 20}, "vertex axis for histograms"};
   ConfigurableAxis axisPhi{"axisPhi", {60, 0.0, constants::math::TwoPI}, "phi axis for histograms"};
@@ -110,12 +116,13 @@ struct FlowTask {
   ConfigurableAxis axisDCAz{"axisDCAz", {200, -2, 2}, "DCA_{z} (cm)"};
   ConfigurableAxis axisDCAxy{"axisDCAxy", {200, -1, 1}, "DCA_{xy} (cm)"};
 
-  Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex;
+  Filter collisionFilter = (nabs(aod::collision::posZ) < cfgCutVertex) && (aod::cent::centFT0C > cfgCentMin) && (aod::cent::centFT0C < cfgCentMax);
   Filter trackFilter = ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true)) && (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPtMin) && (aod::track::pt < cfgCutPtMax) && (aod::track::tpcChi2NCl < cfgCutChi2prTPCcls) && (nabs(aod::track::dcaZ) < cfgCutDCAz);
 
   // Corrections
   TH1D* mEfficiency = nullptr;
   GFWWeights* mAcceptance = nullptr;
+  TList* mGroupAcceptanceList = nullptr;
   bool correctionsLoaded = false;
 
   // Connect to ccdb
@@ -127,6 +134,7 @@ struct FlowTask {
   OutputObj<FlowContainer> fFC{FlowContainer("FlowContainer")};
   OutputObj<GFWWeights> fWeights{GFWWeights("weights")};
   HistogramRegistry registry{"registry"};
+  OutputObj<TList> fGroupNUAList{"GroupNUAList", OutputObjHandlingPolicy::AnalysisObject, OutputObjSourceType::OutputObjSource};
 
   // define global variables
   GFW* fGFW = new GFW();
@@ -150,6 +158,7 @@ struct FlowTask {
   std::unordered_map<int, TH2*> gHadronicRate;
   ctpRateFetcher mRateFetcher;
   TH2* gCurrentHadronicRate;
+  std::vector<std::shared_ptr<GFWWeights>> groupNUAWeightPtr;
 
   using AodCollisions = soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs, aod::Mults>>;
   using AodTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::TracksDCA>>;
@@ -248,6 +257,27 @@ struct FlowTask {
       fWeights->Init(true, false);
     }
 
+    TList* groupNUAWeightlist = new TList();
+    groupNUAWeightlist->SetOwner(true);
+    fGroupNUAList.setObject(groupNUAWeightlist);
+
+    if (cfgOutputGroupNUAWeights) {
+      groupNUAWeightPtr.resize(cfgGroupSplitRunNumber.value.size() + 1);
+      for (uint i = 0; i < cfgGroupSplitRunNumber.value.size() + 1; i++) {
+        GFWWeights* groupweight = nullptr;
+        if (i < cfgGroupSplitRunNumber.value.size())
+          groupweight = new GFWWeights(Form("groupweight_%d", cfgGroupSplitRunNumber.value[i]));
+        else
+          groupweight = new GFWWeights(Form("groupweight_last"));
+
+        groupweight->SetPtBins(nPtBins, ptBins);
+        groupweight->Init(true, false);
+        groupNUAWeightlist->Add(groupweight);
+        std::shared_ptr<GFWWeights> sharePtrGroupWeight(groupweight);
+        groupNUAWeightPtr[i] = sharePtrGroupWeight;
+      }
+    }
+
     // add in FlowContainer to Get boostrap sample automatically
     TObjArray* oba = new TObjArray();
     oba->Add(new TNamed("ChGap22", "ChGap22"));
@@ -506,6 +536,13 @@ struct FlowTask {
       else
         LOGF(warning, "Could not load acceptance weights from %s (%p)", cfgAcceptance.value.c_str(), (void*)mAcceptance);
     }
+    if (cfgAcceptanceGroup.value.empty() == false) {
+      mGroupAcceptanceList = ccdb->getForTimeStamp<TList>(cfgAcceptance, timestamp);
+      if (mGroupAcceptanceList == nullptr) {
+        LOGF(fatal, "Could not load grouped acceptance weights from %s", cfgAcceptanceGroup.value.c_str());
+      }
+      LOGF(info, "Loaded grouped acceptance weights from %s (%p)", cfgAcceptanceGroup.value.c_str(), (void*)mGroupAcceptanceList);
+    }
     if (cfgEfficiency.value.empty() == false) {
       mEfficiency = ccdb->getForTimeStamp<TH1D>(cfgEfficiency, timestamp);
       if (mEfficiency == nullptr) {
@@ -516,7 +553,7 @@ struct FlowTask {
     correctionsLoaded = true;
   }
 
-  bool setCurrentParticleWeights(float& weight_nue, float& weight_nua, float phi, float eta, float pt, float vtxz)
+  bool setCurrentParticleWeights(float& weight_nue, float& weight_nua, float phi, float eta, float pt, float vtxz, int groupNUAIndex = 0)
   {
     float eff = 1.;
     if (mEfficiency)
@@ -526,6 +563,14 @@ struct FlowTask {
     if (eff == 0)
       return false;
     weight_nue = 1. / eff;
+    if (cfgAcceptanceGroupUse) {
+      if (mGroupAcceptanceList && mGroupAcceptanceList->At(groupNUAIndex)) {
+        weight_nua = reinterpret_cast<GFWWeights*>(mGroupAcceptanceList->At(groupNUAIndex))->GetNUA(phi, eta, vtxz);
+      } else {
+        weight_nua = 1;
+      }
+      return true;
+    }
     if (mAcceptance)
       weight_nua = mAcceptance->GetNUA(phi, eta, vtxz);
     else
@@ -717,6 +762,17 @@ struct FlowTask {
     if (cfgUseNch)
       independent = static_cast<float>(tracks.size());
 
+    int groupNUAIndex = 0;
+    if (cfgOutputGroupNUAWeights || cfgAcceptanceGroupUse) {
+      for (uint i = 0; i < cfgGroupSplitRunNumber.value.size(); i++) {
+        if (currentRunNumber < cfgGroupSplitRunNumber.value.at(i)) {
+          break;
+        } else {
+          groupNUAIndex++;
+        }
+      }
+    }
+
     for (const auto& track : tracks) {
       if (!trackSelected(track))
         continue;
@@ -733,7 +789,16 @@ struct FlowTask {
           fWeights->Fill(track.phi(), track.eta(), vtxz, track.pt(), cent, 0);
         }
       }
-      if (!setCurrentParticleWeights(weff, wacc, track.phi(), track.eta(), track.pt(), vtxz))
+      if (cfgOutputGroupNUAWeights) {
+        if (cfgOutputNUAWeightsRefPt) {
+          if (withinPtRef) {
+            groupNUAWeightPtr[groupNUAIndex]->Fill(track.phi(), track.eta(), vtxz, track.pt(), cent, 0);
+          }
+        } else {
+          groupNUAWeightPtr[groupNUAIndex]->Fill(track.phi(), track.eta(), vtxz, track.pt(), cent, 0);
+        }
+      }
+      if (!setCurrentParticleWeights(weff, wacc, track.phi(), track.eta(), track.pt(), vtxz, groupNUAIndex))
         continue;
       registry.fill(HIST("hPt"), track.pt());
       if (withinPtRef) {