From 9163b25c8d5d432e0485ddf947a3fc8365451812 Mon Sep 17 00:00:00 2001
From: hiro98 <hiro@protagon.space>
Date: Tue, 26 May 2020 20:36:07 +0200
Subject: [PATCH] tweak analysis and add total pT

---
 prog/runcards/pp_phaeno/makefile              |  2 +-
 .../qqgg_proton/MC_DIPHOTON_PROTON.cc         | 27 +++++++++++--------
 .../qqgg_proton/MC_DIPHOTON_PROTON.plot       |  7 +++++
 3 files changed, 24 insertions(+), 12 deletions(-)

diff --git a/prog/runcards/pp_phaeno/makefile b/prog/runcards/pp_phaeno/makefile
index 872cacd..a077475 100644
--- a/prog/runcards/pp_phaeno/makefile
+++ b/prog/runcards/pp_phaeno/makefile
@@ -18,7 +18,7 @@ OUT_DIR = out
 EVENT_COUNT = 100000
 
 # List of Runcards, see RUNCARD_FOLDER
-RUNCARDS = basic with_jets with_jets_and_pT with_pT_and_fragmentation with_pT_and_fragmentation_and_mi
+RUNCARDS = with_pT_and_fragmentation_and_mi basic with_jets with_jets_and_pT with_pT_and_fragmentation
 
 # Runcard Names
 BASE_RUNCARD = Sherpa_Base.yaml
diff --git a/prog/runcards/pp_phaeno/qqgg_proton/MC_DIPHOTON_PROTON.cc b/prog/runcards/pp_phaeno/qqgg_proton/MC_DIPHOTON_PROTON.cc
index 1b9bdc7..9b89c53 100644
--- a/prog/runcards/pp_phaeno/qqgg_proton/MC_DIPHOTON_PROTON.cc
+++ b/prog/runcards/pp_phaeno/qqgg_proton/MC_DIPHOTON_PROTON.cc
@@ -6,8 +6,8 @@
 #include "Rivet/Projections/MissingMomentum.hh"
 #include "Rivet/Projections/PromptFinalState.hh"
 #include "Rivet/Projections/VetoedFinalState.hh"
-#include <iostream>
 #include <fstream>
+#include <iostream>
 
 namespace Rivet {
 
@@ -25,7 +25,7 @@ public:
     // Calorimeter particles for photon isolation
     VisibleFinalState visFS;
     VetoedFinalState calo_fs(visFS);
-    calo_fs.addVetoPairId(PID::MUON);  // we also want photon
+    calo_fs.addVetoPairId(PID::MUON); // we also want photon
     declare(calo_fs, "calo_fs");
 
     // we chain in a prompt final state just to be save
@@ -38,20 +38,22 @@ public:
     double min_pT = 20;
     double eta = 2.5;
 
-    _observables = {"pT", "eta", "cos_theta", "inv_m", "o_angle", "o_angle_cs"};
+    _observables = {"pT",      "eta",        "cos_theta", "inv_m",
+                    "o_angle", "o_angle_cs", "total_pT"};
 
     book(_histos["pT"], "pT", logspace(50, min_pT, energy, true));
     book(_histos["eta"], "eta", 50, -eta, eta);
     book(_histos["cos_theta"], "cos_theta", 50, -1, 1);
-    book(_histos["inv_m"], "inv_m", logspace(50, 1, 2 * energy, true));
+    book(_histos["inv_m"], "inv_m", logspace(50, 2 * min_pT, 2 * energy, true));
     book(_histos["o_angle"], "o_angle", 50, 0, 1);
     book(_histos["o_angle_cs"], "o_angle_cs", 50, 0, 1);
+    book(_histos["total_pT"], "total_pT", logspace(50, .01, 2 * energy));
   }
 
   /// Perform the per-event analysis
   void analyze(const Event &event) {
     Particles photons =
-      apply<IdentifiedFinalState>(event, "Photons").particlesByPt();
+        apply<IdentifiedFinalState>(event, "Photons").particlesByPt();
 
     // make sure that there are only two photons
     if (photons.size() < 2)
@@ -60,7 +62,8 @@ public:
     photons.resize(2);
 
     // Require the two photons to be separated in dR
-    if (deltaR(photons[0], photons[1]) < 0.45) vetoEvent;
+    if (deltaR(photons[0], photons[1]) < 0.45)
+      vetoEvent;
 
     const Particles fs = apply<VetoedFinalState>(event, "calo_fs").particles();
     // Loop over photons and require isolation
@@ -94,17 +97,19 @@ public:
 
     const auto &moms = photons.moms();
     const auto total_momentum =
-      std::accumulate(moms.begin(), moms.end(), FourMomentum(0, 0, 0, 0));
+        std::accumulate(moms.begin(), moms.end(), FourMomentum(0, 0, 0, 0));
 
     obs["inv_m"] = total_momentum.mass();
     obs["o_angle"] = std::abs(tanh((photons[1].eta() - photons[0].eta()) / 2));
 
-    //std::abs(photons[0].theta() + photons[1].theta());
+    // std::abs(photons[0].theta() + photons[1].theta());
 
     obs["o_angle_cs"] = std::abs(
-                                 sinh((photons[0].eta() - photons[1].eta())) * 2.0 * photons[0].pT() *
-                                 photons[1].pT() / sqrt(sqr(obs["inv_m"]) + sqr(total_momentum.pT())) /
-                                 obs["inv_m"]);
+        sinh((photons[0].eta() - photons[1].eta())) * 2.0 * photons[0].pT() *
+        photons[1].pT() / sqrt(sqr(obs["inv_m"]) + sqr(total_momentum.pT())) /
+        obs["inv_m"]);
+
+    obs["total_pT"] = total_momentum.pT();
 
     for (const auto &name : _observables) {
       _histos[name]->fill(obs[name]);
diff --git a/prog/runcards/pp_phaeno/qqgg_proton/MC_DIPHOTON_PROTON.plot b/prog/runcards/pp_phaeno/qqgg_proton/MC_DIPHOTON_PROTON.plot
index 40d8e2e..2e97e64 100644
--- a/prog/runcards/pp_phaeno/qqgg_proton/MC_DIPHOTON_PROTON.plot
+++ b/prog/runcards/pp_phaeno/qqgg_proton/MC_DIPHOTON_PROTON.plot
@@ -43,3 +43,10 @@ Title=scattering angle in CS frame
 XLabel=$\cos\theta^\ast_\text{CS}$
 LogY=0
 # END PLOT
+
+# BEGIN PLOT /MC_DIPHOTON_PROTON/total_pT
+Title=total $\pT$
+XLabel=$\pT$
+LogY=1
+LogX=1
+# END PLOT