add meanPt and deltaPt

Author: Luzhiyongg

PWGCF/Flow/Tasks/flowTask.cxx

@@ -304,6 +304,13 @@ struct FlowTask {
     registry.add("hDCAz", "DCAz after cuts; DCAz (cm); Pt", {HistType::kTH2D, {{200, -0.5, 0.5}, {200, 0, 5}}});
     registry.add("hDCAxy", "DCAxy after cuts; DCAxy (cm); Pt", {HistType::kTH2D, {{200, -0.5, 0.5}, {200, 0, 5}}});
     registry.add("hTrackCorrection2d", "Correlation table for number of tracks table; uncorrected track; corrected track", {HistType::kTH2D, {axisNch, axisNch}});
+    registry.add("hMeanPt", "", {HistType::kTProfile, {axisIndependent}});
+    registry.add("hDeltaPt", "", {HistType::kTProfile, {axisIndependent}});
+    registry.add("hMeanPtWithinGap08", "", {HistType::kTProfile, {axisIndependent}});
+    registry.add("c22_gap08_Weff", "", {HistType::kTProfile, {axisIndependent}});
+    registry.add("c22_gap08_trackMeanPt", "", {HistType::kTProfile, {axisIndependent}});
+    registry.add("PtVariance_partA_WithinGap08", "", {HistType::kTProfile, {axisIndependent}});
+    registry.add("PtVariance_partB_WithinGap08", "", {HistType::kTProfile, {axisIndependent}});
     if (doprocessMCGen) {
       registry.add("MCGen/MChPhi", "#phi distribution", {HistType::kTH1D, {axisPhi}});
       registry.add("MCGen/MChEta", "#eta distribution", {HistType::kTH1D, {axisEta}});
@@ -555,6 +562,25 @@ struct FlowTask {
     return;
   }
 
+  template <char... chars, char... chars2>
+  void fillpTvnProfile(const GFW::CorrConfig& corrconf, const double& sum_pt, const double& WeffEvent, const ConstStr<chars...>& vnWeff, const ConstStr<chars2...>& vnpT, const double& cent)
+  {
+    double meanPt = sum_pt / WeffEvent;
+    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 (std::fabs(val) < 1) {
+        registry.fill(vnWeff, cent, val, dnx * WeffEvent);
+        registry.fill(vnpT, cent, val * meanPt, dnx * WeffEvent);
+      }
+      return;
+    }
+    return;
+  }
+
   template <DataType dt>
   void fillFC(const GFW::CorrConfig& corrconf, const double& cent, const double& rndm)
   {
@@ -898,6 +924,10 @@ struct FlowTask {
 
     // track weights
     float weff = 1, wacc = 1;
+    double weffEvent = 0;
+    double ptSum = 0., ptSum_Gap08 = 0.;
+    double weffEventWithinGap08 = 0., weffEventSquareWithinGap08 = 0.;
+    double sumPtsquareWsquareWithinGap08 = 0., sumPtWsquareWithinGap08 = 0.;
     double nTracksCorrected = 0;
     int magnetfield = 0;
     float independent = cent;
@@ -943,6 +973,7 @@ struct FlowTask {
         continue;
       bool withinPtPOI = (cfgCutPtPOIMin < track.pt()) && (track.pt() < cfgCutPtPOIMax); // within POI pT range
       bool withinPtRef = (cfgCutPtRefMin < track.pt()) && (track.pt() < cfgCutPtRefMax); // within RF pT range
+      bool withinEtaGap08 = (std::abs(track.eta()) < cfgEtaPtPt);
       if (cfgOutputNUAWeights) {
         if (cfgOutputNUAWeightsRefPt) {
           if (withinPtRef)
@@ -978,7 +1009,16 @@ struct FlowTask {
         registry.fill(HIST("hnTPCCrossedRow"), track.tpcNClsCrossedRows());
         registry.fill(HIST("hDCAz"), track.dcaZ(), track.pt());
         registry.fill(HIST("hDCAxy"), track.dcaXY(), track.pt());
+        weffEvent += weff;
+        ptSum += weff * track.pt();
         nTracksCorrected += weff;
+        if (withinEtaGap08) {
+          ptSum_Gap08 += weff * track.pt();
+          sumPtWsquareWithinGap08 += weff * weff * track.pt();
+          sumPtsquareWsquareWithinGap08 += weff * weff * track.pt() * track.pt();
+          weffEventWithinGap08 += weff;
+          weffEventSquareWithinGap08 += weff * weff;
+        }
       }
       if (withinPtRef)
         fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), wacc * weff, 1);
@@ -993,6 +1033,35 @@ struct FlowTask {
     }
     registry.fill(HIST("hTrackCorrection2d"), tracks.size(), nTracksCorrected);
 
+    // additional loop to calculate standard error of pT
+    double deltaPtSum = 0.;
+    double meanPt = ptSum / weffEvent;
+    for (const auto& track : tracks) {
+      if (!setCurrentParticleWeights(weff, wacc, track.phi(), track.eta(), track.pt(), vtxz))
+        continue;
+      deltaPtSum += weff * (track.pt() - meanPt) * (track.pt() - meanPt);
+    }
+    double weffEventDiffWithGap08 = weffEventWithinGap08 * weffEventWithinGap08 - weffEventSquareWithinGap08;
+    // MeanPt
+    if (weffEvent) {
+      registry.fill(HIST("hMeanPt"), independent, ptSum / weffEvent, weffEvent);
+      registry.fill(HIST("hDeltaPt"), independent, std::sqrt(deltaPtSum / (weffEvent - 1)), weffEvent);
+    }
+    if (weffEventWithinGap08)
+      registry.fill(HIST("hMeanPtWithinGap08"), independent, ptSum_Gap08 / weffEventWithinGap08, weffEventWithinGap08);
+    // c22_gap8 * pt_withGap8
+    if (weffEventWithinGap08)
+      fillpTvnProfile(corrconfigs.at(7), ptSum_Gap08, weffEventWithinGap08, HIST("c22_gap08_Weff"), HIST("c22_gap08_trackMeanPt"), independent);
+    // PtVariance
+    if (weffEventDiffWithGap08) {
+      registry.fill(HIST("PtVariance_partA_WithinGap08"), independent,
+                    (ptSum_Gap08 * ptSum_Gap08 - sumPtsquareWsquareWithinGap08) / weffEventDiffWithGap08,
+                    weffEventDiffWithGap08);
+      registry.fill(HIST("PtVariance_partB_WithinGap08"), independent,
+                    (weffEventWithinGap08 * ptSum_Gap08 - sumPtWsquareWithinGap08) / weffEventDiffWithGap08,
+                    weffEventDiffWithGap08);
+    }
+
     // Filling Flow Container
     for (uint l_ind = 0; l_ind < corrconfigs.size(); l_ind++) {
       fillFC<kReco>(corrconfigs.at(l_ind), independent, lRandom);