adding shift correction for event plane

Author: jikim1290

PWGCF/JCorran/Tasks/jEPFlowAnalysis.cxx

@@ -32,7 +32,6 @@ using namespace o2::framework::expressions;
 using namespace std;
 
 using MyCollisions = soa::Join<aod::Collisions, aod::EvSels, aod::Qvectors>;
-
 using MyTracks = aod::Tracks;
 
 struct jEPFlowAnalysis {
@@ -53,6 +52,10 @@ struct jEPFlowAnalysis {
   Configurable<bool> cfgAddEvtSel{"cfgAddEvtSel", true, "Use event selection"};
   Configurable<int> cfgnTotalSystem{"cfgnTotalSystem", 7, "Total number of detectors in qVectorsTable"};
 
+  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"};
@@ -64,6 +67,12 @@ struct jEPFlowAnalysis {
   int RefBId;
   int harmInd;
 
+  int currentRunNumber = -999;
+  int lastRunNumber = -999;
+
+  std::vector<TProfile3D*> shiftprofile{};
+  std::string fullCCDBShiftCorrPath;
+
   template <typename T>
   int GetDetId(const T& name)
   {
@@ -96,20 +105,64 @@ struct jEPFlowAnalysis {
     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.};
+
+    if (cfgShiftCorr) {
+      auto bc = collision.bc_as<aod::BCsWithTimestamps>();
+      currentRunNumber = bc.runNumber();
+      if (currentRunNumber != lastRunNumber) {
+        shiftprofile.clear();
+        for (int i = 2; i < 5; 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 = 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);
 
+      auto deltapsiDet = 0.0;
+      auto deltapsiRefA = 0.0;
+      auto deltapsiRefB = 0.0;
+
+      float weight = 1.0;
+
+      if (cfgShiftCorr) {
+        for (int ishift = 1; ishift <= 10; ishift++) {
+          auto coeffshiftxDet = shiftprofile.at(i - 2)->GetBinContent(shiftprofile.at(i - 2)->FindBin(centrality, 0.5, ishift - 0.5));
+          auto coeffshiftyDet = shiftprofile.at(i - 2)->GetBinContent(shiftprofile.at(i - 2)->FindBin(centrality, 1.5, ishift - 0.5));
+          auto coeffshiftxRefA = shiftprofile.at(i - 2)->GetBinContent(shiftprofile.at(i - 2)->FindBin(centrality, 2.5, ishift - 0.5));
+          auto coeffshiftyRefA = shiftprofile.at(i - 2)->GetBinContent(shiftprofile.at(i - 2)->FindBin(centrality, 3.5, ishift - 0.5));
+          auto coeffshiftxRefB = shiftprofile.at(i - 2)->GetBinContent(shiftprofile.at(i - 2)->FindBin(centrality, 4.5, ishift - 0.5));
+          auto coeffshiftyRefB = shiftprofile.at(i - 2)->GetBinContent(shiftprofile.at(i - 2)->FindBin(centrality, 5.5, ishift - 0.5)); //currently only FT0C/TPCpos/TPCneg
+
+          deltapsiDet += ((1 / (1.0 * ishift)) * (-coeffshiftxDet * TMath::Cos(ishift * static_cast<float>(i) * EPs[0]) + coeffshiftyDet * TMath::Sin(ishift * static_cast<float>(i) * EPs[0])));
+          deltapsiRefA += ((1 / (1.0 * ishift)) * (-coeffshiftxRefA * TMath::Cos(ishift * static_cast<float>(i) * EPs[1]) + coeffshiftyRefA * TMath::Sin(ishift * static_cast<float>(i) * EPs[1])));
+          deltapsiRefB += ((1 / (1.0 * ishift)) * (-coeffshiftxRefB * TMath::Cos(ishift * static_cast<float>(i) * EPs[2]) + coeffshiftyRefB * TMath::Sin(ishift * static_cast<float>(i) * EPs[2])));
+        }
+
+        EPs[0] += deltapsiDet;
+        EPs[1] += deltapsiRefA;
+        EPs[2] += deltapsiRefB;
+      }
+
+      if (cfgSPmethod) weight *= sqrt(pow(coll.qvecRe()[DetId + harmInd], 2) + pow(coll.qvecIm()[DetId + harmInd], 2));
+
       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);
@@ -118,7 +171,7 @@ struct jEPFlowAnalysis {
         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);
+          epAnalysis.FillVnHistograms(i, cent, static_cast<float>(j + 1), track.pt(), vn * weight, vn_sin * weight);
         }
       }
     }