Page MenuHomeHEPForge
Files F8310096

No OneTemporary
Actions

View Options

diff --git a/examples/WRapidityAnalysisROOT.cc b/examples/WRapidityAnalysisROOT.cc
index b57f245..5668bbd 100644
--- a/examples/WRapidityAnalysisROOT.cc
+++ b/examples/WRapidityAnalysisROOT.cc
@@ -1,131 +1,131 @@
//! Simple HEJ analysis to plot y_W with ROOT
//!
//! See README.md for installation and running instructions.
//!
#include <cassert>
#include <cmath>
#include <iterator>
#include <memory> // for std::unique_ptr
#include <iostream> // for std::cout
#include "HEJ/PDG_codes.hh"
#include "TFile.h"
#include "TH1.h"
#include "HEJ/Analysis.hh"
#include "HEJ/Event.hh"
#include "HEJ/Version.hh"
namespace YAML {
class Node;
}
namespace LHEF {
class HEPRUP;
}
namespace {
void Fill2DHistogram(TH1D & h,double x,double wt)
{
const int binnumber=h.GetXaxis()->FindBin(x);
const double binwidth=h.GetXaxis()->GetBinWidth(binnumber);
h.Fill(x,wt/binwidth);
}
class WRapidityAnalysisROOT: public HEJ::Analysis {
private:
TFile hfile_;
TH1D htotal_,hydif_,hyf_,hyb_,hyW_,hyWsampling_;
public:
WRapidityAnalysisROOT(
YAML::Node const & /* config */, LHEF::HEPRUP const & /* heprup */
):hfile_{"HEJ.root","RECREATE",
("Output from "+HEJ::Version::package_name_full()).c_str()},
htotal_{"total","total",1,-0.5,.5},
hydif_{"ydif","ydif",40,0.,10.},
hyf_{"yf","yf",40,-5.,5.},
hyb_{"yb","yb",40,-5.,5.},
hyW_{"yW","yW",40,-5.,5.},
hyWsampling_{"yWsampling","yWsampling",40,-5.,5.}
{
// Setup histograms
htotal_.Sumw2();
hydif_.Sumw2();
hyf_.Sumw2();
hyb_.Sumw2();
hyW_.Sumw2();
}
void fill(
HEJ::Event const & event,
HEJ::Event const & /* FO_event */
) override {
const double weight = event.central().weight;
htotal_.Fill(0.0,weight);
auto const & jets = event.jets(); // jets are sorted in rapidity
// find the W momentum
auto const & outgoing = event.outgoing();
auto the_boson = std::find_if(
outgoing.begin(), outgoing.end(),
- [](auto const & particle) { return std::abs(particle.pid) == HEJ::pid::Wp; }
+ [](auto const & particle) { return std::abs(particle.type) == HEJ::pid::Wp; }
);
// we already checked in `pass_cuts` that there is a boson,
// so the following assert cannot fail
assert(the_boson != outgoing.end());
Fill2DHistogram(hydif_,
std::abs(jets.back().rapidity()-jets.front().rapidity()), weight);
Fill2DHistogram(hyf_, jets.front().rapidity(), weight);
Fill2DHistogram(hyb_, jets.back().rapidity(), weight);
Fill2DHistogram(hyW_, the_boson->rapidity(), weight);
Fill2DHistogram(hyWsampling_, the_boson->rapidity(),1);
}
bool pass_cuts(
HEJ::Event const & event,
HEJ::Event const & /* FO_event */
) override {
if(event.jets().size() < 2)
return false;
auto const & outgoing = event.outgoing();
// Find the W boson
auto the_W_boson = std::find_if(
outgoing.begin(), outgoing.end(),
- [](auto const & particle) { return std::abs(particle.pid) == HEJ::pid::Wp; }
+ [](auto const & particle) { return std::abs(particle.type) == HEJ::pid::Wp; }
);
if(the_W_boson == outgoing.end())
return false;
const auto W_pos = std::distance(outgoing.begin(), the_W_boson);
// Find the decay products
const auto the_decay = event.decays().find(W_pos);
if(the_decay == event.decays().end())
return false;
// Find the charged lepton
const auto the_charged_lepton = std::find_if(
the_decay->second.begin(), the_decay->second.end(),
[](HEJ::Particle const & p) {
return HEJ::is_anylepton(p) && !HEJ::is_anyneutrino(p);
}
);
if(the_charged_lepton == the_decay->second.end())
return false;
return std::abs(the_charged_lepton->rapidity()) > 2.5;
}
void finalise() override {
hfile_.Write();
hfile_.Close();
std::cout << "WRapidityAnalysisROOT bids you farewell " << std::endl; // be polite
}
};
}
extern "C"
__attribute__((visibility("default")))
std::unique_ptr<HEJ::Analysis> make_analysis(
YAML::Node const & config, LHEF::HEPRUP const & heprup
){
return std::make_unique<WRapidityAnalysisROOT>(config, heprup);
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 5:21 PM (12 h, 17 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023615
Default Alt Text
(4 KB)

Event Timeline