[PWGCF] adding shift correction for event plane

PWGCF/JCorran/Core/JEPFlowAnalysis.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.
 //
@@ -8,14 +8,14 @@
 // 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.
 
 #include "JEPFlowAnalysis.h"
 
 using namespace o2;
 using namespace o2::framework;
 using namespace std;
 
-void JEPFlowAnalysis::FillVnHistograms(const Int_t harmN, Float_t cent, Float_t det, Float_t pT, Float_t vn, Float_t vn_sin)
+void JEPFlowAnalysis::fillVnHistograms(const Int_t harmN, Float_t cent, Float_t det, Float_t pT, Float_t vn, Float_t vn_sin)
 {
   switch (harmN) {
     case 2: {
@@ -35,7 +35,7 @@
   }
 }
 
-void JEPFlowAnalysis::FillResolutionHistograms(Float_t cent, Float_t harmN, Float_t ResNumA, Float_t ResNumB, Float_t ResDenom)
+void JEPFlowAnalysis::fillResolutionHistograms(Float_t cent, Float_t harmN, Float_t ResNumA, Float_t ResNumB, Float_t ResDenom)
 {
   mHistRegistry->fill(HIST("fResNumA"), ResNumA, harmN, cent, 1.);
   mHistRegistry->fill(HIST("fResNumB"), ResNumB, harmN, cent, 1.);

PWGCF/JCorran/Core/JEPFlowAnalysis.h

@@ -19,22 +19,22 @@
 // O2 headers. //
 #include "Framework/HistogramRegistry.h"
 
 using namespace o2;
 using namespace o2::framework;
 using namespace std;
 
 class JEPFlowAnalysis
 {
  public:
   JEPFlowAnalysis() = default;
-  void SetHistRegistry(HistogramRegistry* histReg) { mHistRegistry = histReg; }
+  void setHistRegistry(HistogramRegistry* histReg) { mHistRegistry = histReg; }
 
-  void FillHistograms(const Int_t fCentBin, Float_t det, Float_t v2, Float_t v3, Float_t v4);
-  void FillVnHistograms(const Int_t harmN, Float_t fCent, Float_t det, Float_t pT, Float_t vn, Float_t vn_sin);
-  void FillResolutionHistograms(Float_t fCent, Float_t harmN, Float_t ResNumA, Float_t ResNumB, Float_t ResDenom);
+  void fillHistograms(const Int_t fCentBin, Float_t det, Float_t v2, Float_t v3, Float_t v4);
+  void fillVnHistograms(const Int_t harmN, Float_t fCent, Float_t det, Float_t pT, Float_t vn, Float_t vn_sin);
+  void fillResolutionHistograms(Float_t fCent, Float_t harmN, Float_t ResNumA, Float_t ResNumB, Float_t ResDenom);
   TComplex Q(const Int_t harmN, const Int_t p);
 
-  void CreateHistograms()
+  void createHistograms()
   {
     if (!mHistRegistry) {
       LOGF(error, "Histogram registry missing. Quitting...");

PWGCF/JCorran/Tasks/CMakeLists.txt

@@ -30,7 +30,7 @@ o2physics_add_dpl_workflow(flow-nuacreation
                     COMPONENT_NAME Analysis)
 
 o2physics_add_dpl_workflow(epflow-analysis
-                    SOURCES jEPFlowAnalysis.cxx
+                    SOURCES JepFlowAnalysis.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::JCorran
                     COMPONENT_NAME Analysis)
 

PWGCF/JCorran/Tasks/jEPFlowAnalysis.cxx

@@ -8,40 +8,49 @@
 // 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.
+
 /// \author Maxim Virta (maxim.virta@cern.ch)
 /// \since Jul 2024
+/// \brief task for flow measurement
+/// \file jEPFlowAnalysis.cxx
 
-#include "Framework/AnalysisTask.h"
-#include "Framework/RunningWorkflowInfo.h"
-#include "Framework/HistogramRegistry.h"
+#include "JEPFlowAnalysis.h"
 
-#include "Common/DataModel/EventSelection.h"
+#include "FlowJHistManager.h"
+
+#include "Common/Core/EventPlaneHelper.h"
 #include "Common/Core/TrackSelection.h"
-#include "Framework/runDataProcessing.h"
+#include "Common/DataModel/EventSelection.h"
+#include "Common/DataModel/Qvectors.h"
 #include "Common/DataModel/TrackSelectionTables.h"
 
-#include "Common/DataModel/Qvectors.h"
-#include "Common/Core/EventPlaneHelper.h"
+#include "CCDB/BasicCCDBManager.h"
+#include "CCDB/CcdbApi.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/HistogramRegistry.h"
+#include "Framework/RunningWorkflowInfo.h"
+#include "Framework/runDataProcessing.h"
 
-#include "FlowJHistManager.h"
-#include "JEPFlowAnalysis.h"
+#include <string>
+#include <vector>
 
 using namespace o2;
 using namespace o2::framework;
 using namespace o2::framework::expressions;
 using namespace std;
 
 using MyCollisions = soa::Join<aod::Collisions, aod::EvSels, aod::Qvectors>;
-
 using MyTracks = aod::Tracks;
 
-struct jEPFlowAnalysis {
+struct JepFlowAnalysis {
 
-  HistogramRegistry EPFlowHistograms{"EPFlow", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
+  HistogramRegistry epFlowHistograms{"EPFlow", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
   JEPFlowAnalysis epAnalysis;
   EventPlaneHelper helperEP;
   FlowJHistManager histManager;
-  Bool_t debug = kFALSE;
+  bool debug = kFALSE;
+  Service<o2::ccdb::BasicCCDBManager> ccdb;
+  o2::ccdb::CcdbApi ccdbApi;
 
   // Set Configurables here
   struct : ConfigurableGroup {
@@ -53,19 +62,37 @@ struct jEPFlowAnalysis {
   Configurable<bool> cfgAddEvtSel{"cfgAddEvtSel", true, "Use event selection"};
   Configurable<int> cfgnTotalSystem{"cfgnTotalSystem", 7, "Total number of detectors in qVectorsTable"};
 
+  Configurable<int> cfgNmode{"cfgNmode", 1, "the number of modulations to be calculated"};
+
+  Configurable<bool> cfgShiftCorr{"cfgShiftCorr", false, "additional shift correction"};
+  Configurable<std::string> cfgShiftPath{"cfgShiftPath", "Users/j/junlee/Qvector/QvecCalib/Shift", "Path for Shift"};
+  Configurable<bool> cfgSPmethod{"cfgSPmethod", false, "flag for scalar product"};
+
   Configurable<std::string> cfgDetName{"cfgDetName", "FT0C", "The name of detector to be analyzed"};
   Configurable<std::string> cfgRefAName{"cfgRefAName", "TPCPos", "The name of detector for reference A"};
   Configurable<std::string> cfgRefBName{"cfgRefBName", "TPCNeg", "The name of detector for reference B"};
 
   Filter trackFilter = (aod::track::pt > cfgTrackCuts.cfgPtMin) && (aod::track::pt < cfgTrackCuts.cfgPtMax) && (nabs(aod::track::eta) < cfgTrackCuts.cfgEtaMax);
 
-  int DetId;
-  int RefAId;
-  int RefBId;
+  int detId;
+  int refAId;
+  int refBId;
   int harmInd;
 
+  int currentRunNumber = -999;
+  int lastRunNumber = -999;
+
+  int initflow = 2;
+  int nshift = 10;
+  int nsystem = 3;
+  int nqvec = 4;
+  int nstep = 3;
+
+  std::vector<TProfile3D*> shiftprofile{};
+  std::string fullCCDBShiftCorrPath;
+
   template <typename T>
-  int GetDetId(const T& name)
+  int getDetId(const T& name)
   {
     if (name.value == "FT0C") {
       return 0;
@@ -88,37 +115,83 @@ struct jEPFlowAnalysis {
 
   void init(InitContext const&)
   {
-    DetId = GetDetId(cfgDetName);
-    RefAId = GetDetId(cfgRefAName);
-    RefBId = GetDetId(cfgRefBName);
+    detId = getDetId(cfgDetName);
+    refAId = getDetId(cfgRefAName);
+    refBId = getDetId(cfgRefBName);
 
-    epAnalysis.SetHistRegistry(&EPFlowHistograms);
-    epAnalysis.CreateHistograms();
+    epAnalysis.setHistRegistry(&epFlowHistograms);
+    epAnalysis.createHistograms();
   }
 
-  void process(MyCollisions::iterator const& coll, soa::Filtered<MyTracks> const& tracks)
+  void process(MyCollisions::iterator const& coll, soa::Filtered<MyTracks> const& tracks, aod::BCsWithTimestamps const&)
   {
     if (cfgAddEvtSel && (!coll.sel8() || !coll.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV) || !coll.selection_bit(aod::evsel::kNoSameBunchPileup)))
       return;
 
-    Float_t cent = coll.cent();
-    EPFlowHistograms.fill(HIST("FullCentrality"), cent);
-    Float_t EPs[3] = {0.};
-    for (int i = 2; i < 5; i++) {                  // loop over different harmonic orders
-      harmInd = cfgnTotalSystem * 4 * (i - 2) + 3; // harmonic index to access corresponding Q-vector as all Q-vectors are in same vector
-      EPs[0] = helperEP.GetEventPlane(coll.qvecRe()[DetId + harmInd], coll.qvecIm()[DetId + harmInd], i);
-      EPs[1] = helperEP.GetEventPlane(coll.qvecRe()[RefAId + harmInd], coll.qvecIm()[RefAId + harmInd], i);
-      EPs[2] = helperEP.GetEventPlane(coll.qvecRe()[RefBId + harmInd], coll.qvecIm()[RefBId + harmInd], i);
-
-      Float_t resNumA = helperEP.GetResolution(EPs[0], EPs[1], i);
-      Float_t resNumB = helperEP.GetResolution(EPs[0], EPs[2], i);
-      Float_t resDenom = helperEP.GetResolution(EPs[1], EPs[2], i);
-      epAnalysis.FillResolutionHistograms(cent, static_cast<float>(i), resNumA, resNumB, resDenom);
-      for (uint j = 0; j < 3; j++) { // loop over detectors used
-        for (auto& track : tracks) {
-          Float_t vn = TMath::Cos((i) * (track.phi() - EPs[j]));
-          Float_t vn_sin = TMath::Sin((i) * (track.phi() - EPs[j]));
-          epAnalysis.FillVnHistograms(i, cent, static_cast<float>(j + 1), track.pt(), vn, vn_sin);
+    float cent = coll.cent();
+    epFlowHistograms.fill(HIST("FullCentrality"), cent);
+    float eps[3] = {0.};
+
+    if (cfgShiftCorr) {
+      auto bc = coll.bc_as<aod::BCsWithTimestamps>();
+      currentRunNumber = bc.runNumber();
+      if (currentRunNumber != lastRunNumber) {
+        shiftprofile.clear();
+        for (int i = initflow; i < initflow + cfgNmode; i++) {
+          fullCCDBShiftCorrPath = cfgShiftPath;
+          fullCCDBShiftCorrPath += "/v";
+          fullCCDBShiftCorrPath += std::to_string(i);
+          auto objshift = ccdb->getForTimeStamp<TProfile3D>(fullCCDBShiftCorrPath, bc.timestamp());
+          shiftprofile.push_back(objshift);
+        }
+        lastRunNumber = currentRunNumber;
+      }
+    }
+
+    for (int i = initflow; i < initflow + cfgNmode; i++) {                  // loop over different harmonic orders
+      harmInd = cfgnTotalSystem * nqvec * (i - initflow) + nstep; // harmonic index to access corresponding Q-vector as all Q-vectors are in same vector
+      eps[0] = helperEP.GetEventPlane(coll.qvecRe()[detId + harmInd], coll.qvecIm()[detId + harmInd], i);
+      eps[1] = helperEP.GetEventPlane(coll.qvecRe()[refAId + harmInd], coll.qvecIm()[refAId + harmInd], i);
+      eps[2] = helperEP.GetEventPlane(coll.qvecRe()[refBId + harmInd], coll.qvecIm()[refBId + harmInd], i);
+
+      auto deltapsiDet = 0.0;
+      auto deltapsiRefA = 0.0;
+      auto deltapsiRefB = 0.0;
+
+      float weight = 1.0;
+
+      if (cfgShiftCorr) {
+        for (int ishift = 1; ishift <= nshift; ishift++) {
+          auto coeffshiftxDet = shiftprofile.at(i - 2)->GetBinContent(shiftprofile.at(i - 2)->FindBin(cent, 0.5, ishift - 0.5));
+          auto coeffshiftyDet = shiftprofile.at(i - 2)->GetBinContent(shiftprofile.at(i - 2)->FindBin(cent, 1.5, ishift - 0.5));
+          auto coeffshiftxRefA = shiftprofile.at(i - 2)->GetBinContent(shiftprofile.at(i - 2)->FindBin(cent, 2.5, ishift - 0.5));
+          auto coeffshiftyRefA = shiftprofile.at(i - 2)->GetBinContent(shiftprofile.at(i - 2)->FindBin(cent, 3.5, ishift - 0.5));
+          auto coeffshiftxRefB = shiftprofile.at(i - 2)->GetBinContent(shiftprofile.at(i - 2)->FindBin(cent, 4.5, ishift - 0.5));
+          auto coeffshiftyRefB = shiftprofile.at(i - 2)->GetBinContent(shiftprofile.at(i - 2)->FindBin(cent, 5.5, ishift - 0.5)); // currently only FT0C/TPCpos/TPCneg
+
+          deltapsiDet += ((1 / (1.0 * ishift)) * (-coeffshiftxDet * std::cos(ishift * static_cast<float>(i) * eps[0]) + coeffshiftyDet * std::sin(ishift * static_cast<float>(i) * eps[0])));
+          deltapsiRefA += ((1 / (1.0 * ishift)) * (-coeffshiftxRefA * std::cos(ishift * static_cast<float>(i) * eps[1]) + coeffshiftyRefA * std::sin(ishift * static_cast<float>(i) * eps[1])));
+          deltapsiRefB += ((1 / (1.0 * ishift)) * (-coeffshiftxRefB * std::cos(ishift * static_cast<float>(i) * eps[2]) + coeffshiftyRefB * std::sin(ishift * static_cast<float>(i) * eps[2])));
+        }
+
+        eps[0] += deltapsiDet;
+        eps[1] += deltapsiRefA;
+        eps[2] += deltapsiRefB;
+      }
+
+      if (cfgSPmethod)
+        weight *= std::sqrt(std::pow(coll.qvecRe()[detId + harmInd], 2) + std::pow(coll.qvecIm()[detId + harmInd], 2));
+
+      float resNumA = helperEP.GetResolution(eps[0], eps[1], i);
+      float resNumB = helperEP.GetResolution(eps[0], eps[2], i);
+      float resDenom = helperEP.GetResolution(eps[1], eps[2], i);
+      epAnalysis.fillResolutionHistograms(cent, static_cast<float>(i), resNumA, resNumB, resDenom);
+
+      for (uint j = 0; j < nsystem; j++) { // loop over detectors used
+        for (const auto& track : tracks) {
+          float vn = std::cos((i) * (track.phi() - eps[j]));
+          float vn_sin = std::sin((i) * (track.phi() - eps[j]));
+          epAnalysis.fillVnHistograms(i, cent, static_cast<float>(j + 1), track.pt(), vn * weight, vn_sin * weight);
         }
       }
     }
@@ -128,5 +201,5 @@ struct jEPFlowAnalysis {
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
 {
   return WorkflowSpec{
-    adaptAnalysisTask<jEPFlowAnalysis>(cfgc)};
+    adaptAnalysisTask<JepFlowAnalysis>(cfgc)};
 }