tweak analysis and add total pT

This commit is contained in:
hiro98 2020-05-26 20:36:07 +02:00
parent 64562a5e45
commit 9163b25c8d
3 changed files with 24 additions and 12 deletions

View file

@ -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

View file

@ -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]);

View file

@ -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