[PWGCF] flowPtEfficiency: add pT-Nch; flowTask: remove group NUA (#10044)

Author: Luzhiyongg

PWGCF/Flow/Tasks/flowPtEfficiency.cxx

@@ -76,6 +76,7 @@ struct FlowPtEfficiency {
   ConfigurableAxis axisCentrality{"axisCentrality", {VARIABLE_WIDTH, 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90}, "X axis for histograms"};
   ConfigurableAxis axisPhi{"axisPhi", {100, 0.0f, constants::math::TwoPI}, ""};
   ConfigurableAxis axisB{"axisB", {100, 0.0f, 20.0f}, "b (fm)"};
+  ConfigurableAxis axisNch{"axisNch", {6000, 0, 6000}, "N_{ch}"};
 
   // Filter the tracks
   Filter trackFilter = (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPtMin) && (aod::track::pt < cfgCutPtMax) && (aod::track::tpcChi2NCl < cfgCutChi2prTPCcls);
@@ -93,6 +94,8 @@ struct FlowPtEfficiency {
   Filter mccollisionFilter = nabs(aod::mccollision::posZ) < cfgCutVertex;
   using MyMcCollisions = soa::Filtered<aod::McCollisions>;
 
+  Preslice<aod::Tracks> perCollision = aod::track::collisionId;
+
   // Additional filters for tracks
   TrackSelection myTrackSel;
 
@@ -145,11 +148,12 @@ struct FlowPtEfficiency {
     // create histograms
     registry.add("eventCounter", "eventCounter", kTH1F, {axisCounter});
     registry.add("hPtMCRec", "Monte Carlo Reco", {HistType::kTH1D, {axisPt}});
-    registry.add("hPtCentMCRec", "Reco production; pT (GeV/c); centrality (%)", {HistType::kTH2D, {axisPt, axisCentrality}});
+    registry.add("hPtNchMCRec", "Reco production; pT (GeV/c); multiplicity", {HistType::kTH2D, {axisPt, axisNch}});
 
     registry.add("mcEventCounter", "Monte Carlo Truth EventCounter", kTH1F, {axisCounter});
     registry.add("hPtMCGen", "Monte Carlo Truth", {HistType::kTH1D, {axisPt}});
-    registry.add("hPtCentMCGen", "Truth production; pT (GeV/c); centrality (%)", {HistType::kTH2D, {axisPt, axisCentrality}});
+    registry.add("hPtNchMCGen", "Truth production; pT (GeV/c); multiplicity", {HistType::kTH2D, {axisPt, axisNch}});
+    registry.add("numberOfRecoCollisions", "numberOfRecoCollisions", kTH1F, {{10, -0.5f, 9.5f}});
 
     if (cfgFlowEnabled) {
       registry.add("hImpactParameterReco", "hImpactParameterReco", {HistType::kTH1D, {axisB}});
@@ -393,7 +397,7 @@ struct FlowPtEfficiency {
         }
         if (isStable(mcParticle.pdgCode())) {
           registry.fill(HIST("hPtMCRec"), track.pt());
-          registry.fill(HIST("hPtCentMCRec"), track.pt(), centrality);
+          registry.fill(HIST("hPtNchMCRec"), track.pt(), tracks.size());
 
           if (cfgFlowEnabled) {
             bool withinPtPOI = (cfgFlowCutPtPOIMin < track.pt()) && (track.pt() < cfgFlowCutPtPOIMax); // within POI pT range
@@ -425,20 +429,20 @@ struct FlowPtEfficiency {
   }
   PROCESS_SWITCH(FlowPtEfficiency, processReco, "process reconstructed information", true);
 
-  void processSim(MyMcCollisions::iterator const& collision, aod::BCsWithTimestamps const&, soa::SmallGroups<soa::Join<aod::McCollisionLabels, aod::Collisions>> const& collisions, MyMcParticles const& mcParticles)
+  void processSim(MyMcCollisions::iterator const& mcCollision, aod::BCsWithTimestamps const&, soa::SmallGroups<soa::Join<aod::McCollisionLabels, aod::Collisions>> const& collisions, MyMcParticles const& mcParticles, MyTracks const& tracks)
   {
     if (cfgSelRunNumberEnabled) {
-      auto bc = collision.bc_as<aod::BCsWithTimestamps>();
+      auto bc = mcCollision.bc_as<aod::BCsWithTimestamps>();
       int runNumber = bc.runNumber();
       if (!std::count(cfgRunNumberList.value.begin(), cfgRunNumberList.value.end(), runNumber))
         return;
     }
 
-    float imp = collision.impactParameter();
+    float imp = mcCollision.impactParameter();
     float centrality = 0.;
     if (cfgFlowEnabled) {
       registry.fill(HIST("hImpactParameterTruth"), imp);
-      auto bc = collision.bc_as<aod::BCsWithTimestamps>();
+      auto bc = mcCollision.bc_as<aod::BCsWithTimestamps>();
       loadCentVsIPTruth(bc.timestamp());
       centrality = mCentVsIPTruth->GetBinContent(mCentVsIPTruth->GetXaxis()->FindBin(imp));
     }
@@ -448,10 +452,20 @@ struct FlowPtEfficiency {
 
     if (collisions.size() > -1) {
       registry.fill(HIST("mcEventCounter"), 0.5);
+
+      registry.fill(HIST("numberOfRecoCollisions"), collisions.size()); // number of times coll was reco-ed
+
+      std::vector<int> numberOfTracks;
+      for (auto const& collision : collisions) {
+        auto groupedTracks = tracks.sliceBy(perCollision, collision.globalIndex());
+        numberOfTracks.emplace_back(groupedTracks.size());
+      }
+
       for (const auto& mcParticle : mcParticles) {
         if (mcParticle.isPhysicalPrimary() && isStable(mcParticle.pdgCode())) {
           registry.fill(HIST("hPtMCGen"), mcParticle.pt());
-          registry.fill(HIST("hPtCentMCGen"), mcParticle.pt(), centrality);
+          if (collisions.size() > 0)
+            registry.fill(HIST("hPtNchMCGen"), mcParticle.pt(), numberOfTracks[0]);
 
           if (cfgFlowEnabled) {
             bool withinPtPOI = (cfgFlowCutPtPOIMin < mcParticle.pt()) && (mcParticle.pt() < cfgFlowCutPtPOIMax); // within POI pT range

PWGCF/Flow/Tasks/flowTask.cxx

@@ -89,19 +89,15 @@ struct FlowTask {
   O2_DEFINE_CONFIGURABLE(cfgNbootstrap, int, 30, "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(cfgEvSelOccupancy, bool, true, "Occupancy cut")
   O2_DEFINE_CONFIGURABLE(cfgCutOccupancyHigh, int, 500, "High cut on TPC occupancy")
   O2_DEFINE_CONFIGURABLE(cfgCutOccupancyLow, int, 0, "Low cut on TPC occupancy")
   O2_DEFINE_CONFIGURABLE(cfgUseSmallMemory, bool, false, "Use small memory mode")
   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"};
@@ -123,7 +119,6 @@ struct FlowTask {
   // Corrections
   TH1D* mEfficiency = nullptr;
   GFWWeights* mAcceptance = nullptr;
-  TList* mGroupAcceptanceList = nullptr;
   bool correctionsLoaded = false;
 
   // Connect to ccdb
@@ -135,7 +130,6 @@ 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();
@@ -156,7 +150,6 @@ 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::CentFT0CVariant1s, aod::CentFT0Ms, aod::CentFV0As, aod::Mults>>;
   using AodTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::TracksDCA>>;
@@ -261,27 +254,6 @@ 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"));
@@ -505,13 +477,6 @@ 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) {
@@ -522,7 +487,7 @@ struct FlowTask {
     correctionsLoaded = true;
   }
 
-  bool setCurrentParticleWeights(float& weight_nue, float& weight_nua, float phi, float eta, float pt, float vtxz, int groupNUAIndex = 0)
+  bool setCurrentParticleWeights(float& weight_nue, float& weight_nua, float phi, float eta, float pt, float vtxz)
   {
     float eff = 1.;
     if (mEfficiency)
@@ -532,14 +497,7 @@ 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
@@ -764,17 +722,6 @@ 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;
@@ -788,16 +735,7 @@ struct FlowTask {
           fWeights->fill(track.phi(), track.eta(), vtxz, track.pt(), cent, 0);
         }
       }
-      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))
+      if (!setCurrentParticleWeights(weff, wacc, track.phi(), track.eta(), track.pt(), vtxz))
         continue;
       registry.fill(HIST("hPt"), track.pt());
       if (withinPtRef) {