Add differential flow histos

Author: wenyaCern

PWGCF/FemtoDream/Core/femtoDreamCollisionSelection.h

@@ -219,7 +219,7 @@ class FemtoDreamCollisionSelection
     mHistogramQn = registry;
     mHistogramQn->add("Event/centFT0CBeforeQn", "; cent", kTH1F, {{10, 0, 100}});
     mHistogramQn->add("Event/centFT0CAfterQn", "; cent", kTH1F, {{10, 0, 100}});
-    mHistogramQn->add("Event/centVsqn", "; cent; qn", kTH2F, {{10, 0, 100}, {100, 0, 1000}});
+    mHistogramQn->add("Event/centVsqn", "; cent; qn", kTH2F, {{10, 0, 100}, {200, 0, 500}});
     mHistogramQn->add("Event/centVsqnVsSpher", "; cent; qn; Sphericity", kTH3F, {{10, 0, 100}, {100, 0, 1000}, {100, 0, 1}});
     mHistogramQn->add("Event/qnBin", "; qnBin; entries", kTH1F, {{20, 0, 20}});
     mHistogramQn->add("Event/psiEP", "; #Psi_{EP} (deg); entries", kTH1F, {{100, 0, 180}});
@@ -241,6 +241,11 @@ class FemtoDreamCollisionSelection
     mHistogramQn->get<TH1>(HIST("Event/hN2allQn"))->Sumw2();
     mHistogramQn->get<TH1>(HIST("Event/hD2allQn"))->Sumw2();
 
+    mHistogramQn->add("Event/hn2pQn",  ";centrality; pT; eta; #sum Re(p_{2,A} Q_{2,B}^{*})", kTH3F, {{centBins, 0, 100}, {20, 0.2, 5.0}, {32, -0.8, 0.8}});
+    mHistogramQn->add("Event/hd2pQn",  ";centrality; pT; eta; #sum (w_{A} W_{B})",           kTH3F, {{centBins, 0, 100}, {20, 0.2, 5.0}, {32, -0.8, 0.8}});
+    mHistogramQn->add("Event/hEntries",  ";centrality; pT; eta; entries",                    kTH3F, {{centBins, 0, 100}, {20, 0.2, 5.0}, {32, -0.8, 0.8}});
+    mHistogramQn->add("Event/hEvtCount",  ";centrality; count",                              kTH1F, {{centBins, 0, 100}});
+
     if (doQnSeparation) {
       for (int iqn(0); iqn < mumQnBins; ++iqn) {
         hN2.push_back(mHistogramQn->add(("Qn/hN2_" + std::to_string(iqn)).c_str(), ";centrality; #sum Re(Q_{2,A} Q_{2,B}^{*})", kTH1F, {{centBins, 0, 100}}));
@@ -479,6 +484,45 @@ class FemtoDreamCollisionSelection
 
     mHistogramQn->fill(HIST("Event/hN2allQn"), centrality, N2_evt);
     mHistogramQn->fill(HIST("Event/hD2allQn"), centrality, D2_evt);
+
+    for (auto const& trk : tracks) {
+      const double pt  = trk.pt();
+      const double eta = trk.eta();
+      if (pt < ptMin || pt > ptMax) {
+        continue;
+      }
+      
+      const double w   = 1.0; // TODO: NUA/NUE weight if you want
+      const double phi = trk.phi();
+      const double c   = w * TMath::Cos(harmonic * phi);
+      const double s   = w * TMath::Sin(harmonic * phi);
+
+      if (eta > fEtaGap) {
+        double p2A_re = c;  
+        double p2A_im = s;  
+        double wA = w;  
+
+        double n2_trk = p2A_re * Q2B_re + p2A_im * Q2B_im;
+        double d2_trk = wA * WB;
+
+        mHistogramQn->fill(HIST("Event/hn2pQn"), centrality, pt, eta, n2_trk);
+        mHistogramQn->fill(HIST("Event/hd2pQn"), centrality, pt, eta, d2_trk);
+        mHistogramQn->fill(HIST("Event/hEntries"), centrality, pt, eta);
+      } else if (eta < -1* fEtaGap) {
+        double p2B_re = c;  
+        double p2B_im = s;  
+        double wB = w;  
+
+        double n2_trk = Q2A_re * p2B_re + Q2A_im * p2B_im;
+        double d2_trk = WA * wB;
+
+        mHistogramQn->fill(HIST("Event/hn2pQn"), centrality, pt, eta, n2_trk);
+        mHistogramQn->fill(HIST("Event/hd2pQn"), centrality, pt, eta, d2_trk);
+        mHistogramQn->fill(HIST("Event/hEntries"), centrality, pt, eta);
+      }
+    }
+    mHistogramQn->fill(HIST("Event/hEvtCount"), centrality);
+
     if (doQnSeparation && mQnBin >= 0 && mQnBin < numQnBins) {
       std::get<std::shared_ptr<TH1>>(hN2[mQnBin])->Fill(centrality, N2_evt);
       std::get<std::shared_ptr<TH1>>(hD2[mQnBin])->Fill(centrality, D2_evt);