PWGCF: flow: add Run3ITSall7Layers option and a new run-by-run task (#8201)

Luzhiyongg · web-flow · commit 1bc5ff98cff3 · 2024-10-30T02:43:07.000+01:00
diff --git a/PWGCF/Flow/Tasks/CMakeLists.txt b/PWGCF/Flow/Tasks/CMakeLists.txt
@@ -19,6 +19,11 @@ o2physics_add_dpl_workflow(flow-task
                     PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::GFWCore
                     COMPONENT_NAME Analysis)
 
+o2physics_add_dpl_workflow(flow-runbyrun
+                    SOURCES FlowRunbyRun.cxx
+                    PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::GFWCore
+                    COMPONENT_NAME Analysis)
+
 o2physics_add_dpl_workflow(flow-gfw-pbpb
                     SOURCES FlowGFWPbPb.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::GFWCore
diff --git a/PWGCF/Flow/Tasks/FlowRunbyRun.cxx b/PWGCF/Flow/Tasks/FlowRunbyRun.cxx
@@ -0,0 +1,204 @@
+// 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.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+//
+// code author: Zhiyong Lu (zhiyong.lu@cern.ch)
+// jira: PWGCF-254
+// Produce Run-by-Run QA plots and flow analysis for Run3
+
+#include <CCDB/BasicCCDBManager.h>
+#include <cmath>
+#include <vector>
+#include <map>
+#include "Framework/runDataProcessing.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/ASoAHelpers.h"
+#include "Framework/RunningWorkflowInfo.h"
+#include "Framework/HistogramRegistry.h"
+
+#include "Common/DataModel/EventSelection.h"
+#include "Common/Core/TrackSelection.h"
+#include "Common/DataModel/TrackSelectionTables.h"
+#include "Common/DataModel/Centrality.h"
+#include "Common/DataModel/Multiplicity.h"
+
+#include "GFWPowerArray.h"
+#include "GFW.h"
+#include "GFWCumulant.h"
+#include "GFWWeights.h"
+#include "TList.h"
+#include <TProfile.h>
+#include <TRandom3.h>
+
+using namespace o2;
+using namespace o2::framework;
+using namespace o2::framework::expressions;
+
+#define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable<TYPE> NAME{#NAME, DEFAULT, HELP};
+
+struct FlowRunbyRun {
+
+  O2_DEFINE_CONFIGURABLE(cfgCutVertex, float, 10.0f, "Accepted z-vertex range")
+  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")
+  O2_DEFINE_CONFIGURABLE(cfgCutPtRefMax, float, 3.0f, "Maximal pT for ref tracks")
+  O2_DEFINE_CONFIGURABLE(cfgCutPtMin, float, 0.2f, "Minimal pT for all tracks")
+  O2_DEFINE_CONFIGURABLE(cfgCutPtMax, float, 10.0f, "Maximal pT for all tracks")
+  O2_DEFINE_CONFIGURABLE(cfgCutEta, float, 0.8f, "Eta range for tracks")
+  O2_DEFINE_CONFIGURABLE(cfgCutChi2prTPCcls, float, 2.5, "Chi2 per TPC clusters")
+  O2_DEFINE_CONFIGURABLE(cfgCutDCAz, float, 2.0f, "max DCA to vertex z")
+  O2_DEFINE_CONFIGURABLE(cfgUseNch, bool, false, "Use Nch for flow observables")
+  Configurable<std::vector<int>> cfgRunNumbers{"cfgRunNumbers", std::vector<int>{544095, 544098, 544116, 544121, 544122, 544123, 544124}, "Preconfigured run numbers"};
+
+  ConfigurableAxis axisVertex{"axisVertex", {20, -10, 10}, "vertex axis for histograms"};
+  ConfigurableAxis axisPhi{"axisPhi", {60, 0.0, constants::math::TwoPI}, "phi axis for histograms"};
+  ConfigurableAxis axisEta{"axisEta", {40, -1., 1.}, "eta axis for histograms"};
+  ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.2, 2.4, 2.6, 2.8, 3, 3.5, 4, 5, 6, 8, 10}, "pt axis for histograms"};
+  ConfigurableAxis axisIndependent{"axisIndependent", {VARIABLE_WIDTH, 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90}, "X axis for histograms"};
+
+  Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex;
+  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);
+
+  // Connect to ccdb
+  Service<ccdb::BasicCCDBManager> ccdb;
+  Configurable<int64_t> nolaterthan{"ccdb-no-later-than", std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object"};
+  Configurable<std::string> url{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
+
+  // Define output
+  HistogramRegistry registry{"registry"};
+
+  // define global variables
+  GFW* fGFW = new GFW();
+  std::vector<GFW::CorrConfig> corrconfigs;
+  std::vector<int> RunNumbers;                                        // vector of run numbers
+  std::map<int, std::vector<std::shared_ptr<TH1>>> TH1sList;          // map of histograms for all runs
+  std::map<int, std::vector<std::shared_ptr<TProfile>>> ProfilesList; // map of profiles for all runs
+  enum OutputTH1Names {
+    // here are TProfiles for vn-pt correlations that are not implemented in GFW
+    hPhi = 0,
+    hEta,
+    hVtxZ,
+    hMult,
+    hCent,
+    kCount_TH1Names
+  };
+  enum OutputTProfileNames {
+    c22 = 0,
+    c22_gap10,
+    kCount_TProfileNames
+  };
+
+  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>>;
+
+  void init(InitContext const&)
+  {
+    ccdb->setURL(url.value);
+    ccdb->setCaching(true);
+    ccdb->setCreatedNotAfter(nolaterthan.value);
+
+    // Add output histograms to the registry
+    std::vector<int> temp = cfgRunNumbers;
+    RunNumbers = temp;
+    for (auto& runNumber : RunNumbers) {
+      CreateOutputObjectsForRun(runNumber);
+    }
+
+    fGFW->AddRegion("full", -0.8, 0.8, 1, 1);
+    fGFW->AddRegion("refN10", -0.8, -0.5, 1, 1);
+    fGFW->AddRegion("refP10", 0.5, 0.8, 1, 1);
+    corrconfigs.resize(kCount_TProfileNames);
+    corrconfigs[c22] = fGFW->GetCorrelatorConfig("full {2 -2}", "ChFull22", kFALSE);
+    corrconfigs[c22_gap10] = fGFW->GetCorrelatorConfig("refN10 {2} refP10 {-2}", "Ch10Gap22", kFALSE);
+    fGFW->CreateRegions();
+  }
+
+  template <char... chars>
+  void FillProfile(const GFW::CorrConfig& corrconf, std::shared_ptr<TProfile> profile, const double& cent)
+  {
+    double dnx, val;
+    dnx = fGFW->Calculate(corrconf, 0, kTRUE).real();
+    if (dnx == 0)
+      return;
+    if (!corrconf.pTDif) {
+      val = fGFW->Calculate(corrconf, 0, kFALSE).real() / dnx;
+      if (TMath::Abs(val) < 1)
+        profile->Fill(cent, val, dnx);
+      return;
+    }
+    return;
+  }
+
+  void CreateOutputObjectsForRun(int runNumber)
+  {
+    std::vector<std::shared_ptr<TH1>> histos(kCount_TH1Names);
+    histos[hPhi] = registry.add<TH1>(Form("%d/hPhi", runNumber), "", {HistType::kTH1D, {axisPhi}});
+    histos[hEta] = registry.add<TH1>(Form("%d/hEta", runNumber), "", {HistType::kTH1D, {axisEta}});
+    histos[hVtxZ] = registry.add<TH1>(Form("%d/hVtxZ", runNumber), "", {HistType::kTH1D, {axisVertex}});
+    histos[hMult] = registry.add<TH1>(Form("%d/hMult", runNumber), "", {HistType::kTH1D, {{3000, 0.5, 3000.5}}});
+    histos[hCent] = registry.add<TH1>(Form("%d/hCent", runNumber), "", {HistType::kTH1D, {{90, 0, 90}}});
+    TH1sList.insert(std::make_pair(runNumber, histos));
+
+    std::vector<std::shared_ptr<TProfile>> profiles(kCount_TProfileNames);
+    profiles[c22] = registry.add<TProfile>(Form("%d/c22", runNumber), "", {HistType::kTProfile, {axisIndependent}});
+    profiles[c22_gap10] = registry.add<TProfile>(Form("%d/c22_gap10", runNumber), "", {HistType::kTProfile, {axisIndependent}});
+    ProfilesList.insert(std::make_pair(runNumber, profiles));
+  }
+
+  void process(aodCollisions::iterator const& collision, aod::BCsWithTimestamps const&, aodTracks const& tracks)
+  {
+    if (!collision.sel8())
+      return;
+    if (tracks.size() < 1)
+      return;
+    // detect run number
+    auto bc = collision.bc_as<aod::BCsWithTimestamps>();
+    int runNumber = bc.runNumber();
+    if (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);
+    }
+
+    if (TH1sList.find(runNumber) == TH1sList.end()) {
+      LOGF(fatal, "RunNumber %d not found in TH1sList", runNumber);
+      return;
+    }
+
+    TH1sList[runNumber][hVtxZ]->Fill(collision.posZ());
+    TH1sList[runNumber][hMult]->Fill(tracks.size());
+    TH1sList[runNumber][hCent]->Fill(collision.centFT0C());
+
+    fGFW->Clear();
+    const auto cent = collision.centFT0C();
+    float weff = 1, wacc = 1;
+    for (auto& track : tracks) {
+      TH1sList[runNumber][hPhi]->Fill(track.phi());
+      TH1sList[runNumber][hEta]->Fill(track.eta());
+      bool WithinPtPOI = (cfgCutPtPOIMin < track.pt()) && (track.pt() < cfgCutPtPOIMax); // within POI pT range
+      bool WithinPtRef = (cfgCutPtRefMin < track.pt()) && (track.pt() < cfgCutPtRefMax); // within RF pT range
+      if (WithinPtRef) {
+        fGFW->Fill(track.eta(), 1, track.phi(), wacc * weff, 1);
+      }
+    }
+
+    // Filling TProfile
+    for (uint i = 0; i < kCount_TProfileNames; ++i) {
+      FillProfile(corrconfigs[i], ProfilesList[runNumber][i], cent);
+    }
+  }
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  return WorkflowSpec{
+    adaptAnalysisTask<FlowRunbyRun>(cfgc)};
+}
diff --git a/PWGCF/Flow/Tasks/FlowTask.cxx b/PWGCF/Flow/Tasks/FlowTask.cxx
@@ -24,6 +24,7 @@
 
 #include "Common/DataModel/EventSelection.h"
 #include "Common/Core/TrackSelection.h"
+#include "Common/Core/TrackSelectionDefaults.h"
 #include "Common/DataModel/TrackSelectionTables.h"
 #include "Common/DataModel/Centrality.h"
 #include "Common/DataModel/Multiplicity.h"
@@ -58,14 +59,16 @@ struct FlowTask {
   O2_DEFINE_CONFIGURABLE(cfgCutChi2prTPCcls, float, 2.5f, "max chi2 per TPC clusters")
   O2_DEFINE_CONFIGURABLE(cfgCutTPCclu, float, 70.0f, "minimum TPC clusters")
   O2_DEFINE_CONFIGURABLE(cfgCutDCAz, float, 2.0f, "max DCA to vertex z")
+  O2_DEFINE_CONFIGURABLE(cfgTrkSelSwitch, bool, false, "switch for self-defined track selection")
+  O2_DEFINE_CONFIGURABLE(cfgTrkSelRun3ITSMatch, bool, false, "GlobalTrackRun3ITSMatching::Run3ITSall7Layers selection")
+  O2_DEFINE_CONFIGURABLE(cfgRejectionTPCsectorOverlap, bool, true, "rejection for TPC sector overlap")
   O2_DEFINE_CONFIGURABLE(cfgUseAdditionalEventCut, bool, false, "Use additional event cut on mult correlations")
   O2_DEFINE_CONFIGURABLE(cfgTriggerkTVXinTRD, bool, true, "TRD triggered")
   O2_DEFINE_CONFIGURABLE(cfgEvSelkNoSameBunchPileup, bool, true, "rejects collisions which are associated with the same found-by-T0 bunch crossing")
   O2_DEFINE_CONFIGURABLE(cfgEvSelkIsGoodZvtxFT0vsPV, bool, true, "removes collisions with large differences between z of PV by tracks and z of PV from FT0 A-C time difference, use this cut at low multiplicities with caution")
   O2_DEFINE_CONFIGURABLE(cfgEvSelkNoCollInTimeRangeStandard, bool, true, "no collisions in specified time range")
   O2_DEFINE_CONFIGURABLE(cfgEvSelMultCorrelation, bool, true, "Multiplicity correlation cut")
   O2_DEFINE_CONFIGURABLE(cfgEvSelV0AT0ACut, bool, true, "V0A T0A 5 sigma cut")
-  O2_DEFINE_CONFIGURABLE(cfgUseAdditionalTrackCut, bool, false, "Use additional track cut on phi")
   O2_DEFINE_CONFIGURABLE(cfgGetInteractionRate, bool, false, "Get interaction rate from CCDB")
   O2_DEFINE_CONFIGURABLE(cfgUseInteractionRateCut, bool, false, "Use events with low interaction rate")
   O2_DEFINE_CONFIGURABLE(cfgCutIR, float, 50.0, "maximum interaction rate (kHz)")
@@ -141,9 +144,11 @@ struct FlowTask {
   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>>;
 
-  // Additional Event selection cuts - Copy from flowGenericFramework.cxx
+  // Track selection
+  TrackSelection myTrackSel;
   TF1* fPhiCutLow = nullptr;
   TF1* fPhiCutHigh = nullptr;
+  // Additional Event selection cuts - Copy from flowGenericFramework.cxx
   TF1* fMultPVCutLow = nullptr;
   TF1* fMultPVCutHigh = nullptr;
   TF1* fMultCutLow = nullptr;
@@ -280,7 +285,7 @@ struct FlowTask {
     std::vector<std::string> UserDefineGFWCorr = cfgUserDefineGFWCorr;
     std::vector<std::string> UserDefineGFWName = cfgUserDefineGFWName;
     if (!UserDefineGFWCorr.empty() && !UserDefineGFWName.empty()) {
-      for (auto i = 0; i < UserDefineGFWName.size(); i++) {
+      for (uint i = 0; i < UserDefineGFWName.size(); i++) {
         oba->Add(new TNamed(UserDefineGFWName.at(i).c_str(), UserDefineGFWName.at(i).c_str()));
       }
     }
@@ -354,7 +359,7 @@ struct FlowTask {
     if (!UserDefineGFWCorr.empty() && !UserDefineGFWName.empty()) {
       LOGF(info, "User adding GFW CorrelatorConfig:");
       // attentaion: here we follow the index of cfgUserDefineGFWCorr
-      for (auto i = 0; i < UserDefineGFWCorr.size(); i++) {
+      for (uint i = 0; i < UserDefineGFWCorr.size(); i++) {
         if (i >= UserDefineGFWName.size()) {
           LOGF(fatal, "The names you provided are more than configurations. UserDefineGFWName.size(): %d > UserDefineGFWCorr.size(): %d", UserDefineGFWName.size(), UserDefineGFWCorr.size());
           break;
@@ -382,10 +387,17 @@ struct FlowTask {
       fT0AV0ASigma->SetParameters(463.4144, 6.796509e-02, -9.097136e-07, 7.971088e-12, -2.600581e-17);
     }
 
-    if (cfgUseAdditionalTrackCut) {
+    if (cfgRejectionTPCsectorOverlap) {
       fPhiCutLow = new TF1("fPhiCutLow", "0.06/x+pi/18.0-0.06", 0, 100);
       fPhiCutHigh = new TF1("fPhiCutHigh", "0.1/x+pi/18.0+0.06", 0, 100);
     }
+
+    if (cfgTrkSelRun3ITSMatch) {
+      myTrackSel = getGlobalTrackSelectionRun3ITSMatch(TrackSelection::GlobalTrackRun3ITSMatching::Run3ITSall7Layers, TrackSelection::GlobalTrackRun3DCAxyCut::Default);
+    } else {
+      myTrackSel = getGlobalTrackSelectionRun3ITSMatch(TrackSelection::GlobalTrackRun3ITSMatching::Run3ITSibAny, TrackSelection::GlobalTrackRun3DCAxyCut::Default);
+    }
+    myTrackSel.SetMinNClustersTPC(cfgCutTPCclu);
   }
 
   template <char... chars>
@@ -570,11 +582,15 @@ struct FlowTask {
   template <typename TTrack>
   bool trackSelected(TTrack track)
   {
-    return (track.tpcNClsFound() >= cfgCutTPCclu);
+    if (cfgTrkSelSwitch) {
+      return myTrackSel.IsSelected(track);
+    } else {
+      return (track.tpcNClsFound() >= cfgCutTPCclu);
+    }
   }
 
   template <typename TTrack>
-  bool AdditionaltrackSelected(TTrack track, const int field)
+  bool RejectionTPCoverlap(TTrack track, const int field)
   {
     double phimodn = track.phi();
     if (field < 0) // for negative polarity field
@@ -670,7 +686,7 @@ struct FlowTask {
     double sum_ptSquare_wSquare_WithinGap08 = 0., sum_pt_wSquare_WithinGap08 = 0.;
     int Magnetfield = 0;
     double NTracksCorrected = 0;
-    if (cfgUseAdditionalTrackCut) {
+    if (cfgRejectionTPCsectorOverlap) {
       // magnet field dependence cut
       Magnetfield = getMagneticField(bc.timestamp());
     }
@@ -681,7 +697,7 @@ struct FlowTask {
     for (auto& track : tracks) {
       if (!trackSelected(track))
         continue;
-      if (cfgUseAdditionalTrackCut && !AdditionaltrackSelected(track, Magnetfield))
+      if (cfgRejectionTPCsectorOverlap && !RejectionTPCoverlap(track, Magnetfield))
         continue;
       if (cfgOutputNUAWeights)
         fWeights->Fill(track.phi(), track.eta(), vtxz, track.pt(), cent, 0);
