[PWGCF] try EbE EP correction (#10333)

Author: Luzhiyongg

PWGCF/Flow/Tasks/flowQa.cxx

@@ -95,6 +95,7 @@ struct FlowQa {
   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")
+  O2_DEFINE_CONFIGURABLE(cfgUseEPcorrection, bool, false, "Use event plane efficiency correction")
   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"};
@@ -157,6 +158,12 @@ struct FlowQa {
   ctpRateFetcher mRateFetcher;
   TH2* gCurrentHadronicRate;
 
+  std::vector<TF1*> funcEff;
+  TH1D* hFindPtBin;
+  TF1* funcV2;
+  TF1* funcV3;
+  TF1* funcV4;
+
   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>>;
 
@@ -355,6 +362,31 @@ struct FlowQa {
       }
     }
     fGFW->CreateRegions();
+
+    if (cfgUseEPcorrection) {
+      hFindPtBin = new TH1D("hFindPtBin", "hFindPtBin", 7, 0.2, 3.0);
+      funcEff.resize(7);
+      funcEff[0] = new TF1("funcEff0", "[0]+[1]*x+[2]*(2*x*x - 1)+[3]*(4*x*x*x - 3*x)+[4]*(8*x*x*x*x - 8*x*x + 1)", 0, 3000);
+      funcEff[0]->SetParameters(0.7366162408, -2.11923e-05, 1.5258e-09, -2.23824e-12, 4.53824e-16);
+      funcEff[1] = new TF1("funcEff1", "[0]+[1]*x+[2]*(2*x*x - 1)+[3]*(4*x*x*x - 3*x)+[4]*(8*x*x*x*x - 8*x*x + 1)", 0, 3000);
+      funcEff[1]->SetParameters(0.7742102264, -2.50337e-05, -1.12602e-09, -1.38765e-12, 3.86733e-16);
+      funcEff[2] = new TF1("funcEff2", "[0]+[1]*x+[2]*(2*x*x - 1)+[3]*(4*x*x*x - 3*x)+[4]*(8*x*x*x*x - 8*x*x + 1)", 0, 3000);
+      funcEff[2]->SetParameters(0.7933082148, -2.24092e-05, -2.55079e-09, -8.59327e-13, 3.1966e-16);
+      funcEff[3] = new TF1("funcEff3", "[0]+[1]*x+[2]*(2*x*x - 1)+[3]*(4*x*x*x - 3*x)+[4]*(8*x*x*x*x - 8*x*x + 1)", 0, 3000);
+      funcEff[3]->SetParameters(0.8084143879, -1.88185e-05, -3.26378e-09, -6.57021e-13, 2.92597e-16);
+      funcEff[4] = new TF1("funcEff4", "[0]+[1]*x+[2]*(2*x*x - 1)+[3]*(4*x*x*x - 3*x)+[4]*(8*x*x*x*x - 8*x*x + 1)", 0, 3000);
+      funcEff[4]->SetParameters(0.8160584679, -4.0401e-06, -1.5118e-08, 3.22423e-12, -1.21729e-16);
+      funcEff[5] = new TF1("funcEff5", "[0]+[1]*x+[2]*(2*x*x - 1)+[3]*(4*x*x*x - 3*x)+[4]*(8*x*x*x*x - 8*x*x + 1)", 0, 3000);
+      funcEff[5]->SetParameters(0.8245747952, -1.88827e-05, 9.88837e-10, -2.23894e-12, 4.62594e-16);
+      funcEff[6] = new TF1("funcEff6", "[0]+[1]*x+[2]*(2*x*x - 1)+[3]*(4*x*x*x - 3*x)+[4]*(8*x*x*x*x - 8*x*x + 1)", 0, 3000);
+      funcEff[6]->SetParameters(0.8277255051, -1.07687e-05, -3.9408e-09, -8.7277e-13, 3.3104e-16);
+      funcV2 = new TF1("funcV2", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x", 0, 100);
+      funcV2->SetParameters(0.0186111, 0.00351907, -4.38264e-05, 1.35383e-07, -3.96266e-10);
+      funcV3 = new TF1("funcV3", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x", 0, 100);
+      funcV3->SetParameters(0.0174056, 0.000703329, -1.45044e-05, 1.91991e-07, -1.62137e-09);
+      funcV4 = new TF1("funcV4", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x", 0, 100);
+      funcV4->SetParameters(0.008845, 0.000259668, -3.24435e-06, 4.54837e-08, -6.01825e-10);
+    }
   }
 
   template <char... chars>
@@ -670,6 +702,32 @@ struct FlowQa {
     if (cfgUseNch)
       independent = static_cast<float>(tracks.size());
 
+    double psi2Est = 0, psi3Est = 0, psi4Est = 0;
+    float wEPeff = 1;
+    double v2 = 0, v3 = 0, v4 = 0;
+    if (cfgUseEPcorrection) {
+      double q2x = 0, q2y = 0;
+      double q3x = 0, q3y = 0;
+      double q4x = 0, q4y = 0;
+      for (const auto& track : tracks) {
+        bool withinPtRef = (cfgCutPtRefMin < track.pt()) && (track.pt() < cfgCutPtRefMax); // within RF pT rang
+        if (withinPtRef) {
+          q2x += std::cos(2 * track.phi());
+          q2y += std::sin(2 * track.phi());
+          q3x += std::cos(3 * track.phi());
+          q3y += std::sin(3 * track.phi());
+          q4x += std::cos(4 * track.phi());
+          q4y += std::sin(4 * track.phi());
+        }
+      }
+      psi2Est = std::atan2(q2y, q2x) / 2.;
+      psi3Est = std::atan2(q3y, q3x) / 3.;
+      psi4Est = std::atan2(q4y, q4x) / 4.;
+      v2 = funcV2->Eval(cent);
+      v3 = funcV3->Eval(cent);
+      v4 = funcV4->Eval(cent);
+    }
+
     for (const auto& track : tracks) {
       if (!trackSelected(track))
         continue;
@@ -685,6 +743,15 @@ struct FlowQa {
       }
       if (!setCurrentParticleWeights(weff, wacc, track.phi(), track.eta(), track.pt(), vtxz))
         continue;
+      if (cfgUseEPcorrection && withinPtRef) {
+        double fphi = v2 * std::cos(2 * (track.phi() - psi2Est)) + v3 * std::cos(3 * (track.phi() - psi3Est)) + v4 * std::cos(4 * (track.phi() - psi4Est));
+        fphi = (1 + 2 * fphi);
+        int pTBinForEff = hFindPtBin->FindBin(track.pt());
+        if (pTBinForEff >= 1 && pTBinForEff <= 7) {
+          wEPeff = funcEff[pTBinForEff - 1]->Eval(fphi * tracks.size());
+          weff *= wEPeff;
+        }
+      }
       registry.fill(HIST("hPt"), track.pt());
       if (withinPtRef) {
         registry.fill(HIST("hPhi"), track.phi());