Page MenuHomeHEPForge
Files F11222397

No OneTemporary
Actions

View Options

diff --git a/analyses/pluginATLAS/ATLAS_2017_I1591327.cc b/analyses/pluginATLAS/ATLAS_2017_I1591327.cc
--- a/analyses/pluginATLAS/ATLAS_2017_I1591327.cc
+++ b/analyses/pluginATLAS/ATLAS_2017_I1591327.cc
@@ -1,186 +1,183 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/IdentifiedFinalState.hh"
#include "Rivet/Projections/FastJets.hh"
namespace Rivet {
/// Isolated diphoton + X differential cross-sections
class ATLAS_2017_I1591327 : public Analysis {
public:
// Constructor
- ATLAS_2017_I1591327() : Analysis("ATLAS_2017_I1591327") {
-
- }
+ DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2017_I1591327);
// Book histograms and initialise projections before the run
void init() {
const FinalState fs;
declare(fs, "FS");
FastJets fj(fs, FastJets::KT, 0.5);
_area_def = new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec());
fj.useJetArea(_area_def);
declare(fj, "KtJetsD05");
IdentifiedFinalState photonfs(Cuts::abseta < 2.37 && Cuts::pT > 30*GeV);
photonfs.acceptId(PID::PHOTON);
declare(photonfs, "Photon");
// Histograms
book(_h_M , 2, 1, 1);
book(_h_pT , 3, 1, 1);
book(_h_at , 4, 1, 1);
book(_h_phistar, 5, 1, 1);
book(_h_costh , 6, 1, 1);
book(_h_dPhi , 7, 1, 1);
}
// Perform the per-event analysis
void analyze(const Event& event) {
// Require at least 2 photons in final state
const Particles photons = apply<IdentifiedFinalState>(event, "Photon").particlesByPt();
if (photons.size() < 2) vetoEvent;
// Compute the median energy density
_ptDensity.clear();
_sigma.clear();
_Njets.clear();
vector<vector<double> > ptDensities;
vector<double> emptyVec;
ptDensities.assign(ETA_BINS.size()-1, emptyVec);
// Get jets, and corresponding jet areas
const shared_ptr<fastjet::ClusterSequenceArea> clust_seq_area = applyProjection<FastJets>(event, "KtJetsD05").clusterSeqArea();
for (const fastjet::PseudoJet& jet : apply<FastJets>(event, "KtJetsD05").pseudoJets(0.0*GeV)) {
const double aeta = fabs(jet.eta());
const double pt = jet.perp();
const double area = clust_seq_area->area(jet);
if (area < 1e-3) continue;
const int ieta = binIndex(aeta, ETA_BINS);
if (ieta != -1) ptDensities[ieta].push_back(pt/area);
}
// Compute median jet properties over the jets in the event
for (size_t b = 0; b < ETA_BINS.size()-1; ++b) {
double median = 0.0, sigma = 0.0;
int Njets = 0;
if (ptDensities[b].size() > 0) {
std::sort(ptDensities[b].begin(), ptDensities[b].end());
int nDens = ptDensities[b].size();
median = (nDens % 2 == 0) ? (ptDensities[b][nDens/2]+ptDensities[b][(nDens-2)/2])/2 : ptDensities[b][(nDens-1)/2];
sigma = ptDensities[b][(int)(.15865*nDens)];
Njets = nDens;
}
_ptDensity.push_back(median);
_sigma.push_back(sigma);
_Njets.push_back(Njets);
}
// Loop over photons and fill vector of isolated ones
Particles isolated_photons;
for (const Particle& photon : photons) {
// Check if it's a prompt photon (needed for SHERPA 2->5 sample, otherwise I also get photons from hadron decays in jets)
- if (photon.fromDecay()) continue;
+ if (!photon.isPrompt()) continue;
// Remove photons in ECAL crack region
if (inRange(photon.abseta(), 1.37, 1.56)) continue;
const double eta_P = photon.eta();
const double phi_P = photon.phi();
// Compute isolation via particles within an R=0.4 cone of the photon
const Particles fs = apply<FinalState>(event, "FS").particles();
FourMomentum mom_in_EtCone;
for (const Particle& p : fs) {
// Reject if not in cone
if (deltaR(photon.momentum(), p.momentum()) > 0.4) continue;
// Reject if in the 5x7 cell central core
if (fabs(eta_P - p.eta()) < 0.025 * 5 * 0.5 &&
fabs(phi_P - p.phi()) < PI/128. * 7 * 0.5) continue;
// Sum momentum
mom_in_EtCone += p.momentum();
}
// Now figure out the correction (area*density)
const double EtCone_area = M_PI*sqr(0.4) - (7*.025)*(5*M_PI/128.); // cone area - central core rectangle
const double correction = _ptDensity[binIndex(fabs(eta_P), ETA_BINS)] * EtCone_area;
// Discard the photon if there is more than 11 GeV of cone activity
// NOTE: Shouldn't need to subtract photon itself (it's in the central core)
if (mom_in_EtCone.Et() - correction > 11*GeV) continue;
// Add isolated photon to list
isolated_photons.push_back(photon);
}
// Require at least two isolated photons
if (isolated_photons.size() < 2) vetoEvent;
// Select leading pT pair
sortByPt(isolated_photons);
const FourMomentum y1 = isolated_photons[0];
const FourMomentum y2 = isolated_photons[1];
// Leading photon should have pT > 40 GeV, subleading > 30 GeV
if (y1.pT() < 40.*GeV) vetoEvent;
if (y2.pT() < 30.*GeV) vetoEvent;
// Require the two photons to be separated (dR>0.4)
if (deltaR(y1,y2) < 0.4) vetoEvent;
const FourMomentum yy = y1 + y2;
const double Myy = yy.mass();
const double pTyy = yy.pT();
const double dPhiyy = mapAngle0ToPi(y1.phi() - y2.phi());
// phi*
const double costhetastar_ = fabs(tanh(( y1.eta() - y2.eta() ) / 2.));
const double sinthetastar_ = sqrt(1. - pow(costhetastar_, 2));
const double phistar = tan(0.5 * (PI - dPhiyy)) * sinthetastar_;
// a_t
const Vector3 t_hat(y1.x()-y2.x(), y1.y()-y2.y(), 0.);
const double factor = t_hat.mod();
const Vector3 t_hatx(t_hat.x()/factor, t_hat.y()/factor, t_hat.z()/factor);
const Vector3 At(y1.x()+y2.x(), y1.y()+y2.y(), 0.);
// Compute a_t transverse component with respect to t_hat
const double at = At.cross(t_hatx).mod();
// Fill histograms
_h_M->fill(Myy);
_h_pT->fill(pTyy);
_h_dPhi->fill(dPhiyy);
_h_costh->fill(costhetastar_);
_h_phistar->fill(phistar);
_h_at->fill(at);
}
// Normalise histograms etc., after the run
void finalize() {
- const double sf = crossSection() / (femtobarn * sumOfWeights());
+ const double sf = crossSection()/femtobarn / sumOfWeights();
scale({_h_M, _h_pT, _h_dPhi, _h_costh, _h_phistar, _h_at}, sf);
}
private:
Histo1DPtr _h_M, _h_pT, _h_dPhi, _h_costh, _h_phistar, _h_at;
fastjet::AreaDefinition* _area_def;
const vector<double> ETA_BINS = {0.0, 1.5, 3.0};
vector<double> _ptDensity, _sigma, _Njets;
};
- // The hook for the plugin system
DECLARE_RIVET_PLUGIN(ATLAS_2017_I1591327);
}
diff --git a/analyses/pluginATLAS/ATLAS_2017_I1604029.cc b/analyses/pluginATLAS/ATLAS_2017_I1604029.cc
--- a/analyses/pluginATLAS/ATLAS_2017_I1604029.cc
+++ b/analyses/pluginATLAS/ATLAS_2017_I1604029.cc
@@ -1,149 +1,153 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/DressedLeptons.hh"
#include "Rivet/Projections/PromptFinalState.hh"
#include "Rivet/Projections/UnstableFinalState.hh"
namespace Rivet {
- ///@brief: ttbar + gamma at 8 TeV
+
+ /// @brief ttbar + gamma at 8 TeV
class ATLAS_2017_I1604029 : public Analysis {
public:
// Constructor
DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2017_I1604029);
// Book histograms and initialise projections before the run
void init() {
const FinalState fs;
// signal photons
PromptFinalState prompt_ph(Cuts::abspid == PID::PHOTON && Cuts::pT > 15*GeV && Cuts::abseta < 2.37);
declare(prompt_ph, "photons");
- // bare leptons
+ // bare leptons
Cut base_cuts = (Cuts::abseta < 2.7) && (Cuts::pT > 10*GeV);
IdentifiedFinalState bare_leps(base_cuts);
bare_leps.acceptIdPair(PID::MUON);
bare_leps.acceptIdPair(PID::ELECTRON);
declare(bare_leps, "bare_leptons");
// dressed leptons
Cut dressed_cuts = (Cuts::abseta < 2.5) && (Cuts::pT > 25*GeV);
PromptFinalState prompt_mu(base_cuts && Cuts::abspid == PID::MUON);
PromptFinalState prompt_el(base_cuts && Cuts::abspid == PID::ELECTRON);
IdentifiedFinalState all_photons(fs, PID::PHOTON);
DressedLeptons elecs(all_photons, prompt_el, 0.1, dressed_cuts);
declare(elecs, "elecs");
DressedLeptons muons(all_photons, prompt_mu, 0.1, dressed_cuts);
declare(muons, "muons");
// auxiliary projections for 'single-lepton ttbar filter'
PromptFinalState prompt_lep(Cuts::abspid == PID::MUON || Cuts::abspid == PID::ELECTRON);
declare(prompt_lep, "prompt_leps");
declare(UnstableFinalState(), "ufs");
// jets
- FastJets jets(fs, FastJets::ANTIKT, 0.4, JetAlg::NO_MUONS, JetAlg::NO_INVISIBLES);
+ FastJets jets(fs, FastJets::ANTIKT, 0.4, JetAlg::Muons::NONE, JetAlg::Invisibles::NONE);
declare(jets, "jets");
+
// BOOK HISTOGRAMS
- _h["pt"] = bookHisto1D(2,1,1);
- _h["eta"] = bookHisto1D(3,1,1);
+ book(_h["pt"], 2,1,1);
+ book(_h["eta"], 3,1,1);
}
+
// Perform the per-event analysis
void analyze(const Event& event) {
-
- const double weight = event.weight();
// analysis extrapolated to 1-lepton-plus-jets channel, where "lepton" cannot be a tau
- // (i.e. contribution from dileptonic ttbar where one of the leptons is outside
+ // (i.e. contribution from dileptonic ttbar where one of the leptons is outside
// the detector acceptance has been subtracted as a background)
if (applyProjection<PromptFinalState>(event, "prompt_leps").particles().size() != 1) vetoEvent;
for (const auto& p : apply<UnstableFinalState>(event, "ufs").particles()) {
if (p.fromPromptTau()) vetoEvent;
}
// photon selection
Particles photons = applyProjection<PromptFinalState>(event, "photons").particlesByPt();
Particles bare_leps = apply<IdentifiedFinalState>(event, "bare_leptons").particles();
for (const Particle& lep : bare_leps)
ifilter_discard(photons, deltaRLess(lep, 0.1));
if (photons.size() != 1) vetoEvent;
const Particle& photon = photons[0];
// jet selection
Jets jets = apply<JetAlg>(event, "jets").jetsByPt(Cuts::abseta < 2.5 && Cuts::pT > 25*GeV);
// lepton selection
const vector<DressedLepton>& elecs = apply<DressedLeptons>(event, "elecs").dressedLeptons();
const vector<DressedLepton>& all_muons = apply<DressedLeptons>(event, "muons").dressedLeptons();
// jet photon/electron overlap removal
for (const DressedLepton& e : elecs)
ifilter_discard(jets, deltaRLess(e, 0.2, RAPIDITY));
for (const Particle& ph : photons)
ifilter_discard(jets, deltaRLess(ph, 0.1, RAPIDITY));
if (jets.size() < 4) vetoEvent;
// photon-jet minimum deltaR
double mindR_phjet = 999.;
for (Jet jet : jets) {
const double dR_phjet = deltaR(photon, jet);
if (dR_phjet < mindR_phjet) mindR_phjet = dR_phjet;
}
if (mindR_phjet < 0.5) vetoEvent;
// muon jet overlap removal
vector<DressedLepton> muons;
- foreach (DressedLepton mu, all_muons) {
+ for (const DressedLepton& mu : all_muons) {
bool overlaps = false;
- foreach (Jet jet, jets) {
+ for (const Jet& jet : jets) {
if (deltaR(mu, jet) < 0.4) {
overlaps = true;
break;
}
}
if (overlaps) continue;
muons.push_back(mu);
}
- // one electron XOR one muon
+ // one electron XOR one muon
bool isEl = elecs.size() == 1 && muons.size() == 0;
bool isMu = muons.size() == 1 && elecs.size() == 0;
if (!isEl && !isMu) vetoEvent;
// photon-lepton deltaR
double mindR_phlep = deltaR(photon, isEl? elecs[0] : muons[0]);
if (mindR_phlep < 0.7) vetoEvent;
// b-tagging
Jets bjets;
- foreach (Jet jet, jets) {
- if (jet.bTagged(Cuts::pT > 5*GeV)) bjets +=jet;
+ for (const Jet& jet : jets) {
+ if (jet.bTagged(Cuts::pT > 5*GeV)) bjets +=jet;
}
if (bjets.empty()) vetoEvent;
- _h["pt"]->fill(photon.pT()/GeV, weight);
- _h["eta"]->fill(photon.abseta(), weight);
+ _h["pt"]->fill(photon.pT()/GeV);
+ _h["eta"]->fill(photon.abseta());
}
+
// Normalise histograms etc., after the run
void finalize() {
const double normto(crossSection() / femtobarn / sumOfWeights());
for (auto &hist : _h) { scale(hist.second, normto); }
}
+
private:
map<string, Histo1DPtr> _h;
+
};
// The hook for the plugin system
DECLARE_RIVET_PLUGIN(ATLAS_2017_I1604029);
}
diff --git a/analyses/pluginATLAS/ATLAS_2017_I1626105.cc b/analyses/pluginATLAS/ATLAS_2017_I1626105.cc
--- a/analyses/pluginATLAS/ATLAS_2017_I1626105.cc
+++ b/analyses/pluginATLAS/ATLAS_2017_I1626105.cc
@@ -1,130 +1,136 @@
+// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/VetoedFinalState.hh"
#include "Rivet/Projections/IdentifiedFinalState.hh"
#include "Rivet/Projections/PromptFinalState.hh"
#include "Rivet/Projections/DressedLeptons.hh"
#include "Rivet/Projections/FastJets.hh"
namespace Rivet {
+
/// @brief lepton differential ttbar analysis at 8 TeV
class ATLAS_2017_I1626105 : public Analysis {
public:
-
+
DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2017_I1626105);
-
+
void init() {
Cut eta_full = Cuts::abseta < 5.0 && Cuts::pT > 1.0*MeV;
// All final state particles
const FinalState fs;
// Get photons to dress leptons
IdentifiedFinalState photons(fs);
photons.acceptIdPair(PID::PHOTON);
// Projection to find the electrons
PromptFinalState prompt_el(Cuts::abspid == PID::ELECTRON, true);
DressedLeptons elecs(photons, prompt_el, 0.1, (Cuts::abseta < 2.5) && (Cuts::pT > 25*GeV));
DressedLeptons veto_elecs(photons, prompt_el, 0.1, eta_full, false);
declare(elecs, "elecs");
// Projection to find the muons
PromptFinalState prompt_mu(Cuts::abspid == PID::MUON, true);
DressedLeptons muons(photons, prompt_mu, 0.1, (Cuts::abseta < 2.5) && (Cuts::pT > 25*GeV));
DressedLeptons veto_muons(photons, prompt_mu, 0.1, eta_full, false);
declare(muons, "muons");
// Jet clustering.
VetoedFinalState vfs;
vfs.addVetoOnThisFinalState(veto_elecs);
vfs.addVetoOnThisFinalState(veto_muons);
FastJets jets(vfs, FastJets::ANTIKT, 0.4);
jets.useInvisibles();
declare(jets, "jets");
// Book histograms
bookHistos("lep_pt", 1);
bookHistos("lep_eta", 3);
bookHistos("dilep_pt", 5);
bookHistos("dilep_mass", 7);
bookHistos("dilep_rap", 9);
bookHistos("dilep_dphi", 11);
bookHistos("dilep_sumpt", 13);
bookHistos("dilep_sumE", 15);
}
+
void analyze(const Event& event) {
vector<DressedLepton> elecs = sortByPt(apply<DressedLeptons>(event, "elecs").dressedLeptons());
- vector<DressedLepton> muons = sortByPt(apply<DressedLeptons>(event, "muons").dressedLeptons());
+ vector<DressedLepton> muons = sortByPt(apply<DressedLeptons>(event, "muons").dressedLeptons());
Jets jets = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
-
+
// Check overlap of jets/leptons.
for (const Jet& jet : jets) {
for (const DressedLepton& el : elecs) {
if (deltaR(jet, el) < 0.4) delete el;
}
for (const DressedLepton& mu : muons) {
if (deltaR(jet, mu) < 0.4) delete mu;
}
}
if (elecs.empty() || muons.empty()) vetoEvent;
-
- if (elecs[0].charge() == muons[0].charge()) vetoEvent;
-
+
+ if (elecs[0].charge() == muons[0].charge()) vetoEvent;
+
FourMomentum el = elecs[0].momentum();
FourMomentum mu = muons[0].momentum();
FourMomentum ll = elecs[0].momentum() + muons[0].momentum();
-
+
// Fill histograms
- const double weight = event.weight();
- fillHistos("lep_pt", el.pT()/GeV, weight);
- fillHistos("lep_pt", mu.pT()/GeV, weight);
- fillHistos("lep_eta", el.abseta(), weight);
- fillHistos("lep_eta", mu.abseta(), weight);
- fillHistos("dilep_pt", ll.pT()/GeV, weight);
- fillHistos("dilep_mass", ll.mass()/GeV, weight);
- fillHistos("dilep_rap", ll.absrap(), weight);
- fillHistos("dilep_dphi", deltaPhi(el, mu), weight);
- fillHistos("dilep_sumpt", (el.pT() + mu.pT())/GeV, weight);
- fillHistos("dilep_sumE", (el.E() + mu.E())/GeV, weight);
+ fillHistos("lep_pt", el.pT()/GeV);
+ fillHistos("lep_pt", mu.pT()/GeV);
+ fillHistos("lep_eta", el.abseta());
+ fillHistos("lep_eta", mu.abseta());
+ fillHistos("dilep_pt", ll.pT()/GeV);
+ fillHistos("dilep_mass", ll.mass()/GeV);
+ fillHistos("dilep_rap", ll.absrap());
+ fillHistos("dilep_dphi", deltaPhi(el, mu));
+ fillHistos("dilep_sumpt", (el.pT() + mu.pT())/GeV);
+ fillHistos("dilep_sumE", (el.E() + mu.E())/GeV);
}
+
void finalize() {
// Normalize to cross-section
const double sf = crossSection()/femtobarn/sumOfWeights();
for (auto& hist : _h) {
const double norm = 1.0 / hist.second->integral();
- // add overflow to last bin
+ // add overflow to last bin
double overflow = hist.second->overflow().effNumEntries();
hist.second->fillBin(hist.second->numBins() - 1, overflow);
// histogram normalisation
if (hist.first.find("norm") != string::npos) scale(hist.second, norm);
else scale(hist.second, sf);
}
}
+
private:
/// @name Histogram helper functions
//@{
void bookHistos(const std::string name, unsigned int index) {
- _h[name] = bookHisto1D(index, 1, 1);
- _h["norm_" + name] = bookHisto1D(index + 1, 1, 1);
+ book(_h[name], index, 1, 1);
+ book(_h["norm_" + name], index + 1, 1, 1);
}
- void fillHistos(const std::string name, double value, double weight) {
- _h[name]->fill(value, weight);
- _h["norm_" + name]->fill(value, weight);
+ void fillHistos(const std::string name, double value) {
+ _h[name]->fill(value);
+ _h["norm_" + name]->fill(value);
}
map<string, Histo1DPtr> _h;
//@}
};
+
// Declare the class as a hook for the plugin system
DECLARE_RIVET_PLUGIN(ATLAS_2017_I1626105);
+
}
diff --git a/analyses/pluginATLAS/ATLAS_2017_I1644367.cc b/analyses/pluginATLAS/ATLAS_2017_I1644367.cc
--- a/analyses/pluginATLAS/ATLAS_2017_I1644367.cc
+++ b/analyses/pluginATLAS/ATLAS_2017_I1644367.cc
@@ -1,185 +1,187 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/IdentifiedFinalState.hh"
#include "Rivet/Projections/FastJets.hh"
namespace Rivet {
+
/// @brief Isolated triphotons at 8 TeV
class ATLAS_2017_I1644367 : public Analysis {
public:
// Constructor
DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2017_I1644367);
-
+
// Book histograms and initialise projections before the run
void init() {
-
+
const FinalState fs;
- declare(fs, "FS");
+ declare(fs, "FS");
FastJets fj(fs, FastJets::KT, 0.5);
- fj.useJetArea(new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec()));
- declare(fj, "KtJetsD05");
+ fj.useJetArea(new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec()));
+ declare(fj, "KtJetsD05");
- IdentifiedFinalState photonfs(Cuts::abspid == PID::PHOTON && Cuts::abseta < 2.37 && Cuts::pT > 15*GeV);
+ IdentifiedFinalState photonfs(Cuts::abspid == PID::PHOTON && Cuts::abseta < 2.37 && Cuts::pT > 15*GeV);
declare(photonfs, "Photon");
-
- // Histograms
- _h["etg1"] = bookHisto1D( 1, 1, 1);
- _h["etg2"] = bookHisto1D( 2, 1, 1);
- _h["etg3"] = bookHisto1D( 3, 1, 1);
- _h["dphig1g2"] = bookHisto1D( 4, 1, 1);
- _h["dphig1g3"] = bookHisto1D( 5, 1, 1);
- _h["dphig2g3"] = bookHisto1D( 6, 1, 1);
- _h["detag1g2"] = bookHisto1D( 7, 1, 1);
- _h["detag1g3"] = bookHisto1D( 8, 1, 1);
- _h["detag2g3"] = bookHisto1D( 9, 1, 1);
- _h["mg1g2"] = bookHisto1D(10, 1, 1);
- _h["mg1g3"] = bookHisto1D(11, 1, 1);
- _h["mg2g3"] = bookHisto1D(12, 1, 1);
- _h["mg1g2g3"] = bookHisto1D(13, 1, 1);
+
+ // Histograms
+ book(_h["etg1"] , 1, 1, 1);
+ book(_h["etg2"] , 2, 1, 1);
+ book(_h["etg3"] , 3, 1, 1);
+ book(_h["dphig1g2"], 4, 1, 1);
+ book(_h["dphig1g3"], 5, 1, 1);
+ book(_h["dphig2g3"], 6, 1, 1);
+ book(_h["detag1g2"], 7, 1, 1);
+ book(_h["detag1g3"], 8, 1, 1);
+ book(_h["detag2g3"], 9, 1, 1);
+ book(_h["mg1g2"] , 10, 1, 1);
+ book(_h["mg1g3"] , 11, 1, 1);
+ book(_h["mg2g3"] , 12, 1, 1);
+ book(_h["mg1g2g3"] , 13, 1, 1);
}
-
+
+
// Perform the per-event analysis
void analyze(const Event& event) {
- const double weight = event.weight();
- // Require at least 2 photons in final state
- const Particles photons = apply<IdentifiedFinalState>(event, "Photon").particlesByPt(Cuts::abseta < 1.37 || Cuts::abseta > 1.5);
- if (photons.size() < 3) vetoEvent;
+ // Require at least 2 photons in final state
+ const Particles photons = apply<IdentifiedFinalState>(event, "Photon").particlesByPt(Cuts::abseta < 1.37 || Cuts::abseta > 1.5);
+ if (photons.size() < 3) vetoEvent;
- // Get jets, and corresponding jet areas
- vector<vector<double> > ptDensities(ETA_BINS.size()-1);
+ // Get jets, and corresponding jet areas
+ vector<vector<double> > ptDensities(ETA_BINS.size()-1);
FastJets fastjets = apply<FastJets>(event, "KtJetsD05");
- const auto clust_seq_area = fastjets.clusterSeqArea();
- for (const Jet& jet : fastjets.jets()) {
- const double area = clust_seq_area->area(jet);
- if (area < 1e-3) continue;
- const int ieta = binIndex(jet.abseta(), ETA_BINS);
- if (ieta != -1) ptDensities[ieta].push_back(jet.pT()/area);
- }
+ const auto clust_seq_area = fastjets.clusterSeqArea();
+ for (const Jet& jet : fastjets.jets()) {
+ const double area = clust_seq_area->area(jet);
+ if (area < 1e-3) continue;
+ const int ieta = binIndex(jet.abseta(), ETA_BINS);
+ if (ieta != -1) ptDensities[ieta].push_back(jet.pT()/area);
+ }
- // Compute median jet properties over the jets in the event
- vector<double> ptDensity;
- for (size_t b = 0; b < ETA_BINS.size()-1; ++b) {
- double median = 0.0;
- if (ptDensities[b].size() > 0) {
- std::sort(ptDensities[b].begin(), ptDensities[b].end());
- int nDens = ptDensities[b].size();
- median = (nDens % 2 == 0) ? (ptDensities[b][nDens/2]+ptDensities[b][(nDens-2)/2])/2 : ptDensities[b][(nDens-1)/2];
+ // Compute median jet properties over the jets in the event
+ vector<double> ptDensity;
+ for (size_t b = 0; b < ETA_BINS.size()-1; ++b) {
+ double median = 0.0;
+ if (ptDensities[b].size() > 0) {
+ std::sort(ptDensities[b].begin(), ptDensities[b].end());
+ int nDens = ptDensities[b].size();
+ median = (nDens % 2 == 0) ? (ptDensities[b][nDens/2]+ptDensities[b][(nDens-2)/2])/2 : ptDensities[b][(nDens-1)/2];
}
ptDensity.push_back(median);
}
// Loop over photons and fill vector of isolated ones
Particles isolated_photons;
for (const Particle& photon : photons) {
- if (photon.fromDecay()) continue;
-
- // Remove photons in ECAL crack region
- const double eta_P = photon.eta();
- const double phi_P = photon.phi();
-
- // Compute isolation via particles within an R=0.4 cone of the photon
- const Particles fs = apply<FinalState>(event, "FS").particles();
- FourMomentum mom_in_EtCone;
- for (const Particle& p : fs) {
- // Reject if not in cone
- if (deltaR(photon.momentum(), p.momentum()) > 0.4) continue;
- // Reject if in the 5x7 cell central core
+ if (!photon.isPrompt()) continue;
+
+ // Remove photons in ECAL crack region
+ const double eta_P = photon.eta();
+ const double phi_P = photon.phi();
+
+ // Compute isolation via particles within an R=0.4 cone of the photon
+ const Particles fs = apply<FinalState>(event, "FS").particles();
+ FourMomentum mom_in_EtCone;
+ for (const Particle& p : fs) {
+ // Reject if not in cone
+ if (deltaR(photon.momentum(), p.momentum()) > 0.4) continue;
+ // Reject if in the 5x7 cell central core
if (fabs(eta_P - p.eta()) < 0.025 * 5 * 0.5 && fabs(phi_P - p.phi()) < PI/128. * 7 * 0.5) continue;
// Sum momentum
mom_in_EtCone += p.momentum();
}
// Now figure out the correction (area*density)
const double EtCone_area = M_PI*sqr(0.4) - (7*.025)*(5*M_PI/128.); // cone area - central core rectangle
const double correction = ptDensity[binIndex(fabs(eta_P), ETA_BINS)] * EtCone_area;
// Discard the photon if there is more than 11 GeV of cone activity
// NOTE: Shouldn't need to subtract photon itself (it's in the central core)
if (mom_in_EtCone.Et() - correction > 10*GeV) continue;
// Add isolated photon to list
isolated_photons.push_back(photon);
- }///loop over photons
-
- // Require at least two isolated photons
- if (isolated_photons.size() < 3) vetoEvent;
-
- // Select leading pT pair
- sortByPt(isolated_photons);
- const FourMomentum y1 = isolated_photons[0];
- const FourMomentum y2 = isolated_photons[1];
- const FourMomentum y3 = isolated_photons[2];
-
- // Leading photon should have pT > 40 GeV, subleading > 30 GeV
- if (y1.pT() < 27*GeV) vetoEvent;
- if (y2.pT() < 22*GeV) vetoEvent;
- if (y3.pT() < 15*GeV) vetoEvent;
-
+ }///loop over photons
+
+ // Require at least two isolated photons
+ if (isolated_photons.size() < 3) vetoEvent;
+
+ // Select leading pT pair
+ sortByPt(isolated_photons);
+ const FourMomentum y1 = isolated_photons[0];
+ const FourMomentum y2 = isolated_photons[1];
+ const FourMomentum y3 = isolated_photons[2];
+
+ // Leading photon should have pT > 40 GeV, subleading > 30 GeV
+ if (y1.pT() < 27*GeV) vetoEvent;
+ if (y2.pT() < 22*GeV) vetoEvent;
+ if (y3.pT() < 15*GeV) vetoEvent;
+
// Require the two photons to be separated (dR>0.4)
- if (deltaR(y1,y2) < 0.45) vetoEvent;
- if (deltaR(y1,y3) < 0.45) vetoEvent;
- if (deltaR(y2,y3) < 0.45) vetoEvent;
-
-
- const FourMomentum yyy = y1 + y2 + y3;
- const FourMomentum y1y2 = y1 + y2;
- const FourMomentum y1y3 = y1 + y3;
- const FourMomentum y2y3 = y2 + y3;
+ if (deltaR(y1,y2) < 0.45) vetoEvent;
+ if (deltaR(y1,y3) < 0.45) vetoEvent;
+ if (deltaR(y2,y3) < 0.45) vetoEvent;
- const double Myyy = yyy.mass() / GeV;
- const double dPhiy1y2 = mapAngle0ToPi(deltaPhi(y1, y2));
- const double dPhiy1y3 = mapAngle0ToPi(deltaPhi(y1, y3));
- const double dPhiy2y3 = mapAngle0ToPi(deltaPhi(y2, y3));
+ const FourMomentum yyy = y1 + y2 + y3;
+ const FourMomentum y1y2 = y1 + y2;
+ const FourMomentum y1y3 = y1 + y3;
+ const FourMomentum y2y3 = y2 + y3;
- const double dEtay1y2 = fabs(y1.eta() - y2.eta());
- const double dEtay1y3 = fabs(y1.eta() - y3.eta());
- const double dEtay2y3 = fabs(y2.eta() - y3.eta());
+ const double Myyy = yyy.mass() / GeV;
+
+ const double dPhiy1y2 = mapAngle0ToPi(deltaPhi(y1, y2));
+ const double dPhiy1y3 = mapAngle0ToPi(deltaPhi(y1, y3));
+ const double dPhiy2y3 = mapAngle0ToPi(deltaPhi(y2, y3));
+
+ const double dEtay1y2 = fabs(y1.eta() - y2.eta());
+ const double dEtay1y3 = fabs(y1.eta() - y3.eta());
+ const double dEtay2y3 = fabs(y2.eta() - y3.eta());
if(Myyy < 50.) vetoEvent;
// Fill histograms
- _h["etg1"]->fill(y1.pT() / GeV, weight);
- _h["etg2"]->fill(y2.pT() / GeV, weight);
- _h["etg3"]->fill(y3.pT() / GeV, weight);
+ _h["etg1"]->fill(y1.pT() / GeV);
+ _h["etg2"]->fill(y2.pT() / GeV);
+ _h["etg3"]->fill(y3.pT() / GeV);
- _h["dphig1g2"]->fill(dPhiy1y2);
- _h["dphig1g3"]->fill(dPhiy1y3);
- _h["dphig2g3"]->fill(dPhiy2y3);
+ _h["dphig1g2"]->fill(dPhiy1y2);
+ _h["dphig1g3"]->fill(dPhiy1y3);
+ _h["dphig2g3"]->fill(dPhiy2y3);
- _h["detag1g2"]->fill(dEtay1y2);
- _h["detag1g3"]->fill(dEtay1y3);
- _h["detag2g3"]->fill(dEtay2y3);
+ _h["detag1g2"]->fill(dEtay1y2);
+ _h["detag1g3"]->fill(dEtay1y3);
+ _h["detag2g3"]->fill(dEtay2y3);
- _h["mg1g2"]->fill(y1y2.mass() / GeV, weight);
- _h["mg1g3"]->fill(y1y3.mass() / GeV, weight);
- _h["mg2g3"]->fill(y2y3.mass() / GeV, weight);
+ _h["mg1g2"]->fill(y1y2.mass() / GeV);
+ _h["mg1g3"]->fill(y1y3.mass() / GeV);
+ _h["mg2g3"]->fill(y2y3.mass() / GeV);
- _h["mg1g2g3"]->fill(Myyy, weight);
- }
-
-
- // Normalise histograms etc., after the run
- void finalize() {
- const double sf = crossSection() / (femtobarn * sumOfWeights());
+ _h["mg1g2g3"]->fill(Myyy);
+ }
+
+
+ // Normalise histograms etc., after the run
+ void finalize() {
+ const double sf = crossSection() / (femtobarn * sumOfWeights());
for (auto &hist : _h) { scale(hist.second, sf); }
- }
-
-
- private:
-
- map<string, Histo1DPtr> _h;
- const vector<double> ETA_BINS = { 0.0, 1.5, 3.0 };
-
- };
-
- // The hook for the plugin system
- DECLARE_RIVET_PLUGIN(ATLAS_2017_I1644367);
+ }
+
+
+ private:
+
+ map<string, Histo1DPtr> _h;
+ const vector<double> ETA_BINS = { 0.0, 1.5, 3.0 };
+
+ };
+
+
+ DECLARE_RIVET_PLUGIN(ATLAS_2017_I1644367);
+
}
diff --git a/analyses/pluginATLAS/ATLAS_2017_I1645627.cc b/analyses/pluginATLAS/ATLAS_2017_I1645627.cc
--- a/analyses/pluginATLAS/ATLAS_2017_I1645627.cc
+++ b/analyses/pluginATLAS/ATLAS_2017_I1645627.cc
@@ -1,154 +1,153 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/VisibleFinalState.hh"
#include "Rivet/Projections/VetoedFinalState.hh"
#include "Rivet/Projections/PromptFinalState.hh"
#include "Rivet/Projections/FastJets.hh"
namespace Rivet {
/// @brief Isolated photon + jets at 13 TeV
class ATLAS_2017_I1645627 : public Analysis {
public:
// Constructor
DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2017_I1645627);
// Book histograms and initialise projections before the run
void init() {
const FinalState fs;
// calorimeter particles
VisibleFinalState visFS(fs);
VetoedFinalState calo_fs(visFS);
calo_fs.addVetoPairId(PID::MUON);
declare(calo_fs, "calo");
// Voronoi eta-phi tessellation with KT jets, for ambient energy density calculation
- FastJets fj(fs, FastJets::KT, 0.5, JetAlg::NO_MUONS, JetAlg::NO_INVISIBLES); // E-scheme used by default;
+ FastJets fj(fs, FastJets::KT, 0.5, JetAlg::Muons::NONE, JetAlg::Invisibles::NONE); // E-scheme used by default;
fj.useJetArea(new fastjet::AreaDefinition(fastjet::voronoi_area, fastjet::VoronoiAreaSpec(1.0)));
declare(fj, "KtJetsD05");
// photon
PromptFinalState photonfs(Cuts::abspid == PID::PHOTON && Cuts::abseta < 2.37 && Cuts::pT > 125*GeV);
declare(photonfs, "photons");
// Jets
- FastJets jetpro(fs, FastJets::ANTIKT, 0.4, JetAlg::NO_MUONS, JetAlg::NO_INVISIBLES);
+ FastJets jetpro(fs, FastJets::ANTIKT, 0.4, JetAlg::Muons::NONE, JetAlg::Invisibles::NONE);
declare(jetpro, "Jets");
-
- _h_photon_pt = bookHisto1D(1, 1, 1);
- _h_jet_pt = bookHisto1D(2, 1, 1);
- _h_phjet_dphi = bookHisto1D(3, 1, 1);
- _h_phjet_mass = bookHisto1D(4, 1, 1);
- _h_phjet_costheta = bookHisto1D(5, 1, 1);
+
+ // Histograms
+ book(_h_photon_pt , 1, 1, 1);
+ book(_h_jet_pt , 2, 1, 1);
+ book(_h_phjet_dphi , 3, 1, 1);
+ book(_h_phjet_mass , 4, 1, 1);
+ book(_h_phjet_costheta, 5, 1, 1);
}
- size_t getEtaBin(double eta) const {
- const double abseta = fabs(eta);
- return binIndex(abseta, _eta_bins_areaoffset);
+ size_t getEtaBin(double eta) const {
+ return binIndex(fabs(eta), _eta_bins_areaoffset);
}
-
+
// Perform the per-event analysis
void analyze(const Event& event) {
- const double weight = event.weight();
-
// Get the photon
const Particles& photons = apply<PromptFinalState>(event, "photons").particlesByPt(Cuts::abseta < 1.37 || Cuts::abseta > 1.56);
if (photons.empty()) vetoEvent;
const FourMomentum photon = photons[0].momentum();
// Get the jet
Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 100*GeV && Cuts::absrap < 2.37);
ifilter_discard(jets, deltaRLess(photon, 0.8));
if (jets.empty()) vetoEvent;
FourMomentum leadingJet = jets[0].momentum();
-
+
// Compute the jet pT densities
vector< vector<double> > ptDensities(_eta_bins_areaoffset.size()-1);
FastJets fastjets = apply<FastJets>(event, "KtJetsD05");
const auto clust_seq_area = fastjets.clusterSeqArea();
for (const Jet& jet : fastjets.jets()) {
- const double area = clust_seq_area->area(jet); // Implicit call to pseudojet().
- //const double area2 = (clust_seq_area->area_4vector(jet)).perp(); // Area definition used in egammaTruthParticles.
+ const double area = clust_seq_area->area(jet); // Implicit call to pseudojet().
+ //const double area2 = (clust_seq_area->area_4vector(jet)).perp(); // Area definition used in egammaTruthParticles.
if (area > 1e-3 && jet.abseta() < _eta_bins_areaoffset.back()) {
ptDensities.at(getEtaBin(jet.abseta())) += jet.pT()/area;
}
}
// Compute the median event energy density
vector<double> ptDensity;
for (size_t b = 0; b < _eta_bins_areaoffset.size()-1; ++b) {
ptDensity += ptDensities[b].empty() ? 0 : median(ptDensities[b]);
}
// Compute photon isolation with a standard ET cone
FourMomentum mom_in_EtCone;
const Particles calo_fs = apply<VetoedFinalState>(event, "calo").particles();
const double iso_dr = 0.4;
for (const Particle& p : calo_fs) {
// Check if it's in the cone of .4
if (sqrt(2.0*(cosh(p.eta()-photon.eta()) - cos(p.phi()-photon.phi()))) >= iso_dr) continue;
// Increment sum
mom_in_EtCone += p.momentum();
}
-
+
// Remove the photon energy from the isolation
mom_in_EtCone -= photon;
// Figure out the correction (area*density)
const double etcone_area = PI*iso_dr*iso_dr;
const double correction = ptDensity[getEtaBin(photon.abseta())] * etcone_area;
// Require photon to be isolated
if ((mom_in_EtCone.Et()-correction) > (0.0042*photon.pT() + 10*GeV)) vetoEvent;
-
+
// Fill histos
const double photon_pt = photon.pT()/GeV;
const double jet_pt = leadingJet.pT()/GeV;
const double phjet_dphi = deltaPhi(photon, leadingJet);
const double photon_eta = photon.eta();
const double jet_y = leadingJet.rapidity();
- _h_photon_pt->fill(photon_pt, weight);
- _h_jet_pt->fill(jet_pt, weight);
- _h_phjet_dphi->fill(phjet_dphi, weight);
+ _h_photon_pt->fill(photon_pt);
+ _h_jet_pt->fill(jet_pt);
+ _h_phjet_dphi->fill(phjet_dphi);
double dy = fabs(jet_y-photon_eta);
double phjet_costheta = tanh(dy/2.);
double phjet_mass= (photon+leadingJet).mass()/GeV;
if (phjet_mass <= 450.) vetoEvent;
if (fabs(photon_eta + jet_y) >= 2.37) vetoEvent;
if (phjet_costheta >= 0.83) vetoEvent;
- _h_phjet_costheta->fill(phjet_costheta, weight);
- _h_phjet_mass->fill(phjet_mass, weight);
+ _h_phjet_costheta->fill(phjet_costheta);
+ _h_phjet_mass->fill(phjet_mass);
}
/// Normalise histograms etc., after the run
void finalize() {
const double sf = crossSection() / sumOfWeights();
scale(_h_photon_pt, sf);
scale(_h_jet_pt, sf);
scale(_h_phjet_dphi, sf);
scale(_h_phjet_mass, sf);
scale(_h_phjet_costheta, sf);
}
private:
Histo1DPtr _h_photon_pt, _h_jet_pt;
Histo1DPtr _h_phjet_dphi, _h_phjet_mass, _h_phjet_costheta;
const vector<double> _eta_bins_areaoffset = {0.0, 1.5, 3.0, 4.0, 5.0};
};
// The hook for the plugin system
DECLARE_RIVET_PLUGIN(ATLAS_2017_I1645627);
+
}
diff --git a/analyses/pluginATLAS/ATLAS_2018_I1646686.cc b/analyses/pluginATLAS/ATLAS_2018_I1646686.cc
--- a/analyses/pluginATLAS/ATLAS_2018_I1646686.cc
+++ b/analyses/pluginATLAS/ATLAS_2018_I1646686.cc
@@ -1,299 +1,298 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/VetoedFinalState.hh"
#include "Rivet/Projections/IdentifiedFinalState.hh"
#include "Rivet/Projections/PromptFinalState.hh"
#include "Rivet/Projections/DressedLeptons.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/PartonicTops.hh"
#include "Rivet/Math/LorentzTrans.hh"
namespace Rivet {
/// @brief All-hadronic ttbar at 13 TeV
class ATLAS_2018_I1646686 : public Analysis {
public:
/// Constructor
DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2018_I1646686);
/// Book histograms and initialise projections before the run
void init() {
//histogram booking
book(_h["inclusive"],2,1,1);
bookHistograms("t_y", 1, true);
bookHistograms("t_pt", 2, true);
bookHistograms("t1_pt", 3);
bookHistograms("t1_y", 4);
bookHistograms("t2_pt", 5);
bookHistograms("t2_y", 6);
bookHistograms("tt_m", 7);
bookHistograms("tt_pt", 8);
bookHistograms("tt_y", 9);
bookHistograms("tt_chi", 10);
bookHistograms("tt_yboost", 11);
bookHistograms("tt_pout", 12);
bookHistograms("tt_dPhi", 13);
bookHistograms("tt_Ht", 14);
bookHistograms("tt_cosThStar", 15);
// Projections
Cut dressed_lep = (Cuts::abseta < 2.5) && (Cuts::pT >= 25*GeV);
Cut eta_full = (Cuts::abseta < 5.0);
// All final state particles
FinalState fs(eta_full);
// Get photons to dress leptons
IdentifiedFinalState photons(fs);
photons.acceptIdPair(PID::PHOTON);
// Projection to find the electrons
IdentifiedFinalState el_id(fs);
el_id.acceptIdPair(PID::ELECTRON);
PromptFinalState electrons(el_id);
electrons.acceptTauDecays(true);
DressedLeptons dressedelectrons(photons, electrons, 0.1, dressed_lep);
declare(dressedelectrons, "elecs");
DressedLeptons ewdressedelectrons(photons, electrons, 0.1, eta_full);
// Projection to find the muons
IdentifiedFinalState mu_id(fs);
mu_id.acceptIdPair(PID::MUON);
PromptFinalState muons(mu_id);
muons.acceptTauDecays(true);
DressedLeptons dressedmuons(photons, muons, 0.1, dressed_lep);
declare(dressedmuons, "muons");
DressedLeptons ewdressedmuons(photons, muons, 0.1, eta_full);
// Projection to find neutrinos for vetoing in jets
IdentifiedFinalState nu_id(fs);
nu_id.acceptNeutrinos();
PromptFinalState neutrinos(nu_id);
neutrinos.acceptTauDecays(true);
// Jet clustering.
VetoedFinalState lvfs(fs);
lvfs.addVetoOnThisFinalState(mu_id);
lvfs.addVetoOnThisFinalState(nu_id);
VetoedFinalState vfs;
vfs.addVetoOnThisFinalState(ewdressedelectrons);
vfs.addVetoOnThisFinalState(ewdressedmuons);
vfs.addVetoOnThisFinalState(neutrinos);
FastJets jets(vfs, FastJets::ANTIKT, 0.4);
jets.useInvisibles();
declare(jets, "jets");
FastJets ljets(lvfs, FastJets::ANTIKT, 1.0);
ljets.useInvisibles();
declare(ljets, "ljets" );
PartonicTops partonTops;
declare(partonTops, "partonicTops");
}
void analyze(const Event& event) {
// Parton-level top quarks
const Particles partonicTops = apply<PartonicTops>( event, "partonicTops").particlesByPt();
FourMomentum top, tbar;
bool foundT = false, foundTBar = false;
for (const Particle& ptop : partonicTops) {
const int pid = ptop.pid();
if (pid == PID::TQUARK) {
top = ptop.momentum();
foundT = true;
} else if (pid == -PID::TQUARK) {
tbar = ptop.momentum();
foundTBar = true;
}
}
FourMomentum t1_parton, t2_parton, ttbar_parton;
if ( foundT && foundTBar ) {
t1_parton = top.pT2() > tbar.pT2() ? top : tbar;
t2_parton = top.pT2() > tbar.pT2() ? tbar : top;
ttbar_parton = t1_parton + t2_parton;
if ( t1_parton.pT() > 500*GeV && t2_parton.pT() > 350*GeV) {
const double chi_parton = calcChi(t1_parton, t2_parton);
const double cosThetaStar_parton = abs(calcCosThetaStar(t1_parton, t2_parton));
const double pout_parton = abs(calcPout(t1_parton, t2_parton));
const double dPhi_parton = deltaPhi(t1_parton, t2_parton);
const int randomChoice = rand() % 2;
const FourMomentum& randomTopParton = (randomChoice == 0) ? t1_parton : t2_parton;
fillParton("t_pt", randomTopParton.pT()/GeV);
fillParton("t_y", randomTopParton.absrap());
fillParton("t1_pt", t1_parton.pT()/GeV);
fillParton("t1_y", t1_parton.absrap());
fillParton("t2_pt", t2_parton.pT()/GeV);
fillParton("t2_y", t2_parton.absrap());
fillParton("tt_m", ttbar_parton.mass()/GeV);
fillParton("tt_pt", ttbar_parton.pT()/GeV);
fillParton("tt_Ht", (t1_parton.pT() + t2_parton.pT())/GeV);
fillParton("tt_y", ttbar_parton.absrap());
fillParton("tt_yboost", 0.5 * abs(t1_parton.rapidity() + t2_parton.rapidity()));
fillParton("tt_chi", chi_parton);
fillParton("tt_cosThStar", cosThetaStar_parton);
fillParton("tt_pout", pout_parton/GeV);
fillParton("tt_dPhi", dPhi_parton);
}
}
// Get and veto on dressed leptons
const vector<DressedLepton> dressedElectrons = apply<DressedLeptons>(event, "elecs").dressedLeptons();
const vector<DressedLepton> dressedMuons = apply<DressedLeptons>(event, "muons").dressedLeptons();
if (!dressedElectrons.empty()) vetoEvent;
if (!dressedMuons.empty()) vetoEvent;
// Get jets
const Jets& all_jets = apply<FastJets>( event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
const FastJets& ljets_fj = apply<FastJets>( event, "ljets");
const Jets all_ljets = ljets_fj.jetsByPt();
// Trim the large-R jets
Jets trimmedJets;
fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::kt_algorithm, 0.2), fastjet::SelectorPtFractionMin(0.05));
for (const Jet& jet : all_ljets)
trimmedJets += ljets_fj.trimJet(jet, trimmer);
trimmedJets = sortByPt(trimmedJets);
// Check large-R jets
Jets ljets;
vector<bool> b_tagged;
for (const Jet& jet : trimmedJets) {
if (jet.pT() < 250 * GeV) continue;
if (jet.pT() > 3000 * GeV) continue;
if (jet.mass() > jet.pT()) continue;
if (jet.abseta() > 2.0 ) continue;
ljets += jet;
b_tagged += jet.bTagged();
}
if (all_jets.size() < 2) vetoEvent;
if (ljets.size() < 2) vetoEvent;
// Identify top and anti top, compute some event variables
const FourMomentum ttbar = ljets[0].momentum() + ljets[1].momentum();
const FourMomentum t1 = ljets[0].momentum();
const FourMomentum t2 = ljets[1].momentum();
const double chi = calcChi(t1, t2);
const double cosThetaStar = abs(calcCosThetaStar(t1, t2));
const double pout = abs(calcPout(t1, t2));
const double dPhi = deltaPhi(t1, t2);
if ( t2.pT() < 350*GeV) vetoEvent;
if ( t1.pT() < 500*GeV) vetoEvent;
// b-tagging for particle done on large-R jets
if (!(b_tagged[0] && b_tagged[1])) vetoEvent;
// Continues with signal region cuts
if ( abs(t1.mass() - 172.5 * GeV) > 50*GeV ) vetoEvent;
if ( abs(t2.mass() - 172.5 * GeV) > 50*GeV ) vetoEvent;
_h["inclusive"]->fill(0);
fillHistograms("t1_pt", t1.pT()/GeV);
fillHistograms("t1_y", t1.absrap());
fillHistograms("t2_pt", t2.pT()/GeV);
fillHistograms("t2_y", t2.absrap());
fillHistograms("tt_m", ttbar.mass()/TeV);
fillHistograms("tt_pt", ttbar.pT()/GeV);
fillHistograms("tt_Ht", (t1.pT() + t2.pT())/GeV);
fillHistograms("tt_y", ttbar.absrap());
fillHistograms("tt_yboost", 0.5 * abs(t1.rapidity() + t2.rapidity()));
fillHistograms("tt_chi", chi);
fillHistograms("tt_cosThStar", cosThetaStar);
fillHistograms("tt_pout", pout/GeV);
fillHistograms("tt_dPhi", dPhi);
}
void finalize() {
// Normalize histograms
const double sf = crossSection() / sumOfWeights();
for (auto &hist : _h) {
scale(hist.second, sf);
if ((hist.first.find("_norm") != string::npos) && hist.second->integral(false)>0) hist.second->normalize(1.0, false);
}
}
double calcChi(const FourMomentum& t1, const FourMomentum& t2) {
double ystar = 0.5 * (t1.rapidity()-t2.rapidity());
double chi = exp( 2 * abs(ystar));
return chi;
}
double calcCosThetaStar(const FourMomentum& t1, const FourMomentum& t2) {
FourMomentum ttbar = t1 + t2;
LorentzTransform centreOfMassTrans;
ttbar.setX(0);
ttbar.setY(0);
centreOfMassTrans.setBetaVec( -ttbar.betaVec() );
FourMomentum t1_star = centreOfMassTrans.transform(t1);
double cosThetaStar = t1_star.pz()/t1_star.p3().mod();
return cosThetaStar;
}
double calcPout(const FourMomentum& t1, const FourMomentum& t2) {
Vector3 t1V = t1.p3();
Vector3 t2V = t2.p3();
Vector3 zUnit = Vector3(0., 0., 1.);
Vector3 vPerp = zUnit.cross(t1V);
double pout = vPerp.dot(t2V)/vPerp.mod();
return pout;
}
private:
map<string, Histo1DPtr> _h;
//some functions for booking, filling and scaling the histograms
void fillHistograms(std::string name, double value) {
_h[name]->fill(value);
_h[name + "_norm"]->fill(value);
}
void fillParton(std::string name, double value) {
_h[name + "_parton"]->fill(value);
_h[name + "_parton_norm"]->fill(value);
}
void bookHistograms(const std::string name, unsigned int index, bool onlyParton = false) {
if (!onlyParton) {
book(_h[name], index, 1, 1 );
book(_h[name + "_norm"], index + 13, 1, 1 );
}
book(_h[name + "_parton"], index + 82, 1, 1 );
book(_h[name + "_parton_norm"], index + 97, 1, 1 );
}
};
DECLARE_RIVET_PLUGIN(ATLAS_2018_I1646686);
-
}
diff --git a/analyses/pluginCMS/CMS_2014_I1303894.cc b/analyses/pluginCMS/CMS_2014_I1303894.cc
--- a/analyses/pluginCMS/CMS_2014_I1303894.cc
+++ b/analyses/pluginCMS/CMS_2014_I1303894.cc
@@ -1,236 +1,235 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/VetoedFinalState.hh"
#include "Rivet/Projections/PromptFinalState.hh"
#include "Rivet/Projections/InvMassFinalState.hh"
#include "Rivet/Projections/MissingMomentum.hh"
#include "Rivet/Projections/WFinder.hh"
#include "Rivet/Projections/DressedLeptons.hh"
namespace Rivet {
/// @brief Differential cross-section of W bosons + jets in pp collisions at sqrt(s)=7 TeV
/// @author Darin Baumgartel (darinb@cern.ch)
///
/// Based on Rivet analysis originally created by Anil Singh (anil@cern.ch), Lovedeep Saini (lovedeep@cern.ch)
class CMS_2014_I1303894 : public Analysis {
public:
/// Constructor
CMS_2014_I1303894()
: Analysis("CMS_2014_I1303894")
{ }
// Book histograms and initialise projections before the run
void init() {
// Prompt leptons only, no test on nu flavour.
// Projections
const FinalState fs;
declare(fs, "FS");
MissingMomentum missing(fs);
declare(missing, "MET");
PromptFinalState pfs(fs);
IdentifiedFinalState bareMuons(pfs);
bareMuons.acceptIdPair(PID::MUON);
DressedLeptons muonClusters(fs, bareMuons, -1); //, Cuts::open(), false, false);
declare(muonClusters, "muonClusters");
IdentifiedFinalState neutrinos(pfs);
neutrinos.acceptIdPair(PID::NU_MU);
declare(neutrinos, "neutrinos");
VetoedFinalState jetFS(fs);
jetFS.addVetoOnThisFinalState(muonClusters);
jetFS.addVetoOnThisFinalState(neutrinos);
jetFS.vetoNeutrinos();
- FastJets JetProjection(jetFS, FastJets::ANTIKT, 0.5);
- JetProjection.useInvisibles(false);
- declare(JetProjection, "Jets");
+ FastJets jetprojection(jetFS, FastJets::ANTIKT, 0.5);
+ declare(jetprojection, "Jets");
// Histograms
book(_histDPhiMuJet1 ,1,1,1);
book(_histDPhiMuJet2 ,2,1,1);
book(_histDPhiMuJet3 ,3,1,1);
book(_histDPhiMuJet4 ,4,1,1);
book(_histEtaJet1 ,5,1,1);
book(_histEtaJet2 ,6,1,1);
book(_histEtaJet3 ,7,1,1);
book(_histEtaJet4 ,8,1,1);
book(_histHT1JetInc ,9,1,1);
book(_histHT2JetInc ,10,1,1);
book(_histHT3JetInc ,11,1,1);
book(_histHT4JetInc ,12,1,1);
book(_histJet30MultExc ,13,1,1);
book(_histJet30MultInc ,14,1,1);
book(_histPtJet1 ,15,1,1);
book(_histPtJet2 ,16,1,1);
book(_histPtJet3 ,17,1,1);
book(_histPtJet4 ,18,1,1);
// Counters
book(_n_1jet, "n_1jet");
book(_n_2jet, "n_2jet");
book(_n_3jet, "n_3jet");
book(_n_4jet, "n_4jet");
book(_n_inclusivebinsummation, "n_inclusivebinsummation");
}
void analyze(const Event& event) {
// Get the dressed muon
const DressedLeptons& muonClusters = apply<DressedLeptons>(event, "muonClusters");
int nmu = muonClusters.dressedLeptons().size();
if (nmu < 1) vetoEvent;
DressedLepton dressedmuon = muonClusters.dressedLeptons()[0];
if (dressedmuon.momentum().abseta() > 2.1) vetoEvent;
if (dressedmuon.momentum().pT() < 25.0*GeV) vetoEvent;
// Get the muon neutrino
//const Particles& neutrinos = apply<FinalState>(event, "neutrinos").particlesByPt();
// Check that the muon and neutrino are not decay products of tau
if (dressedmuon.constituentLepton().hasAncestor( PID::TAU)) vetoEvent;
if (dressedmuon.constituentLepton().hasAncestor(-PID::TAU)) vetoEvent;
// Get the missing momentum
const MissingMomentum& met = apply<MissingMomentum>(event, "MET");
const double ptmet = met.visibleMomentum().pT();
const double phimet = (-met.visibleMomentum()).phi();
// Calculate MET and MT(mu,MET), and remove events with MT < 50 GeV
const double ptmuon = dressedmuon.pT();
const double phimuon = dressedmuon.phi();
const double mt_mumet = sqrt(2*ptmuon*ptmet*(1.0 - cos(phimet-phimuon)));
// Remove events in MT < 50 region
if (mt_mumet < 50*GeV) vetoEvent;
// Loop over jets and fill pt/eta/phi quantities in vectors
const Jets& jets_filtered = apply<FastJets>(event, "Jets").jetsByPt(0.0*GeV);
vector<float> finaljet_pT_list, finaljet_eta_list, finaljet_phi_list;
double htjets = 0.0;
for (size_t ii = 0; ii < jets_filtered.size(); ++ii) {
// Jet pT/eta/phi
double jet_pt = jets_filtered[ii].pT();
double jet_eta = jets_filtered[ii].eta();
double jet_phi = jets_filtered[ii].phi();
// Kinemetic cuts for jet acceptance
if (fabs(jet_eta) > 2.4) continue;
if (jet_pt < 30.0*GeV) continue;
if (deltaR(dressedmuon, jets_filtered[ii]) < 0.5) continue;
// Add jet to jet list and increases the HT variable
finaljet_pT_list.push_back(jet_pt);
finaljet_eta_list.push_back(jet_eta);
finaljet_phi_list.push_back(jet_phi);
htjets += fabs(jet_pt);
}
// Filling of histograms:
// Fill as many jets as there are into the exclusive jet multiplicity
if (!finaljet_pT_list.empty())
_histJet30MultExc->fill(finaljet_pT_list.size());
for (size_t ij = 0; ij < finaljet_pT_list.size(); ++ij) {
_histJet30MultInc->fill(ij+1);
_n_inclusivebinsummation->fill();
}
if (finaljet_pT_list.size() >= 1) {
_histPtJet1->fill(finaljet_pT_list[0]);
_histEtaJet1->fill(fabs(finaljet_eta_list[0]));
_histDPhiMuJet1->fill(deltaPhi(finaljet_phi_list[0], phimuon));
_histHT1JetInc->fill(htjets);
_n_1jet->fill();
}
if (finaljet_pT_list.size() >= 2) {
_histPtJet2->fill(finaljet_pT_list[1]);
_histEtaJet2->fill(fabs(finaljet_eta_list[1]));
_histDPhiMuJet2->fill(deltaPhi(finaljet_phi_list[1], phimuon));
_histHT2JetInc->fill(htjets);
_n_2jet->fill();
}
if (finaljet_pT_list.size() >= 3) {
_histPtJet3->fill(finaljet_pT_list[2]);
_histEtaJet3->fill(fabs(finaljet_eta_list[2]));
_histDPhiMuJet3->fill(deltaPhi(finaljet_phi_list[2], phimuon));
_histHT3JetInc->fill(htjets);
_n_3jet->fill();
}
if (finaljet_pT_list.size() >=4 ) {
_histPtJet4->fill(finaljet_pT_list[3]);
_histEtaJet4->fill(fabs(finaljet_eta_list[3]));
_histDPhiMuJet4->fill(deltaPhi(finaljet_phi_list[3], phimuon));
_histHT4JetInc-> fill(htjets);
_n_4jet->fill();
}
}
// Finalize the histograms.
void finalize() {
const double inclusive_cross_section = crossSection();
const double norm_1jet_histo = inclusive_cross_section*dbl(*_n_1jet)/sumOfWeights();
const double norm_2jet_histo = inclusive_cross_section*dbl(*_n_2jet)/sumOfWeights();
const double norm_3jet_histo = inclusive_cross_section*dbl(*_n_3jet)/sumOfWeights();
const double norm_4jet_histo = inclusive_cross_section*dbl(*_n_4jet)/sumOfWeights();
const double norm_incmultiplicity = inclusive_cross_section*dbl(*_n_inclusivebinsummation)/sumOfWeights();
normalize(_histJet30MultExc, norm_1jet_histo);
normalize(_histJet30MultInc, norm_incmultiplicity);
normalize(_histPtJet1, norm_1jet_histo);
normalize(_histHT1JetInc, norm_1jet_histo);
normalize(_histEtaJet1, norm_1jet_histo);
normalize(_histDPhiMuJet1, norm_1jet_histo);
normalize(_histPtJet2, norm_2jet_histo);
normalize(_histHT2JetInc, norm_2jet_histo);
normalize(_histEtaJet2, norm_2jet_histo);
normalize(_histDPhiMuJet2, norm_2jet_histo);
normalize(_histPtJet3, norm_3jet_histo);
normalize(_histHT3JetInc, norm_3jet_histo);
normalize(_histEtaJet3, norm_3jet_histo);
normalize(_histDPhiMuJet3, norm_3jet_histo);
normalize(_histPtJet4, norm_4jet_histo);
normalize(_histHT4JetInc, norm_4jet_histo);
normalize(_histEtaJet4, norm_4jet_histo);
normalize(_histDPhiMuJet4, norm_4jet_histo);
}
private:
Histo1DPtr _histJet30MultExc, _histJet30MultInc;
Histo1DPtr _histPtJet1, _histPtJet2, _histPtJet3, _histPtJet4;
Histo1DPtr _histEtaJet1, _histEtaJet2, _histEtaJet3, _histEtaJet4;
Histo1DPtr _histDPhiMuJet1, _histDPhiMuJet2, _histDPhiMuJet3, _histDPhiMuJet4;
Histo1DPtr _histHT1JetInc, _histHT2JetInc, _histHT3JetInc, _histHT4JetInc;
CounterPtr _n_1jet, _n_2jet, _n_3jet, _n_4jet, _n_inclusivebinsummation;
};
DECLARE_RIVET_PLUGIN(CMS_2014_I1303894);
}
diff --git a/analyses/pluginCMS/CMS_2015_I1346843.cc b/analyses/pluginCMS/CMS_2015_I1346843.cc
--- a/analyses/pluginCMS/CMS_2015_I1346843.cc
+++ b/analyses/pluginCMS/CMS_2015_I1346843.cc
@@ -1,112 +1,113 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/ChargedLeptons.hh"
#include "Rivet/Projections/NeutralFinalState.hh"
#include "Rivet/Projections/IdentifiedFinalState.hh"
namespace Rivet {
/// Differential cross-section of FSR photons in Z decays
class CMS_2015_I1346843 : public Analysis {
public:
/// Constructor
DEFAULT_RIVET_ANALYSIS_CTOR(CMS_2015_I1346843);
/// Book histograms and initialise projections before the run
void init() {
Cut c_photons = Cuts::pT >= 5.0*GeV && (Cuts::etaIn(-2.5, 1.4) || Cuts::etaIn(1.6, 2.5));
IdentifiedFinalState photons(c_photons);
photons.acceptId(PID::PHOTON);
declare(photons, "PHOTFS");
Cut c_muons = Cuts::pT > 9*GeV && Cuts::abseta < 2.4;
IdentifiedFinalState muons(c_muons);
muons.acceptIdPair(PID::MUON);
declare(muons, "MUFS");
book(_hist_pho_et ,1, 1, 1); // photon transverse energy
book(_hist_pho_et_wide ,1, 2, 1); // photon transverse energy (0.5 < dr < 3.0)
book(_hist_pho_et_close ,1, 3, 1); // photon transverse energy (0.05 < dr < 0.5)
book(_hist_pho_et_lqt ,1, 4, 1); // photon transverse energy (q_T < 10)
book(_hist_pho_et_hqt ,1, 5, 1); // photon transverse energy (q_T > 50)
book(_hist_pho_dr ,2, 1, 1); // delta_R
book(_hist_pho_dr_lqt ,2, 2, 1); // delta_R (q_T < 10)
book(_hist_pho_dr_hqt ,2, 3, 1); // delta_R (q_T > 50)
}
// Perform the per-event analysis
void analyze(const Event& event) {
const Particles muons = apply<IdentifiedFinalState>(event, "MUFS").particlesByPt();
if (muons.size() < 2) vetoEvent;
if (muons[0].pT()/GeV < 31) vetoEvent;
if (muons[0].charge()*muons[1].charge() > 0) vetoEvent;
const double mZ = (muons[0].momentum() + muons[1].momentum()).mass();
if (!inRange(mZ, 30*GeV, 87*GeV)) vetoEvent;
const Particles photons = apply<IdentifiedFinalState>(event, "PHOTFS").particlesByPt();
// We want the photon with the highest pT that does not come from a decay
- for(const Particle& p : photons) {
- if (p.fromDecay() || !p.isStable()) continue;
+ for (const Particle& p : photons) {
+ if (!p.isDirect()) continue;
+ if (!p.isStable()) continue;
const double dR = std::min(deltaR(p, muons[0]), deltaR(p, muons[1]) );
if (!inRange(dR, 0.05, 3.0)) continue;
// Calculate the three-body (mu,mu,gamma) transverse momentum
const double qT = (muons[0].mom() + muons[1].mom() + p.mom()).pT();
// Fill the analysis histograms
_hist_pho_et->fill(p.pT()/GeV, 1.0);
_hist_pho_dr->fill(dR, 1.0);
(dR <= 0.5 ? _hist_pho_et_close : _hist_pho_et_wide)->fill(p.pT()/GeV, 1.0);
if (qT / GeV < 10.) {
_hist_pho_et_lqt->fill(p.pT()/GeV, 1.0);
_hist_pho_dr_lqt->fill(dR, 1.0);
}
if (qT / GeV > 50.) {
_hist_pho_et_hqt->fill(p.pT()/GeV, 1.0);
_hist_pho_dr_hqt->fill(dR, 1.0);
}
break; // Exit the loop since we found the highest pT lepton already
}
}
/// Normalise histograms etc., after the run
void finalize() {
scale(_hist_pho_et, crossSection() / sumOfWeights());
scale(_hist_pho_et_wide, crossSection() / sumOfWeights());
scale(_hist_pho_et_close, crossSection() / sumOfWeights());
scale(_hist_pho_et_lqt, crossSection() / sumOfWeights());
scale(_hist_pho_et_hqt, crossSection() / sumOfWeights());
scale(_hist_pho_dr, crossSection() / sumOfWeights());
scale(_hist_pho_dr_lqt, crossSection() / sumOfWeights());
scale(_hist_pho_dr_hqt, crossSection() / sumOfWeights());
}
private:
Histo1DPtr _hist_pho_et;
Histo1DPtr _hist_pho_et_wide, _hist_pho_et_close;
Histo1DPtr _hist_pho_et_lqt, _hist_pho_et_hqt;
Histo1DPtr _hist_pho_dr;
Histo1DPtr _hist_pho_dr_lqt, _hist_pho_dr_hqt;
};
DECLARE_RIVET_PLUGIN(CMS_2015_I1346843);
}
diff --git a/analyses/pluginCMS/CMS_2015_I1397174.cc b/analyses/pluginCMS/CMS_2015_I1397174.cc
--- a/analyses/pluginCMS/CMS_2015_I1397174.cc
+++ b/analyses/pluginCMS/CMS_2015_I1397174.cc
@@ -1,386 +1,386 @@
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/PartonicTops.hh"
#include "Rivet/Projections/VetoedFinalState.hh"
#include "Rivet/Projections/FastJets.hh"
namespace Rivet {
/// Fully leptonic partonic ttbar analysis
class CMS_2015_I1397174 : public Analysis {
public:
/// Minimal constructor
CMS_2015_I1397174()
: Analysis("CMS_2015_I1397174") { }
/// @name Analysis methods
//@{
/// Set up projections and book histograms
void init() {
// Parton level top quarks
declare(PartonicTops(PartonicTops::DecayMode::E_MU, false), "PartonTops");
// Find jets not related to the top/W decays
VetoedFinalState vfs;
vfs.addDecayProductsVeto(PID::WPLUSBOSON);
vfs.addDecayProductsVeto(PID::WMINUSBOSON);
FastJets fj(vfs, FastJets::ANTIKT, 0.5, JetAlg::Muons::ALL, JetAlg::Invisibles::ALL);
declare(fj, "Jets");
// Book histograms
book(_hVis_nJet30_abs , 1, 1, 1);
book(_hVis_nJet30 , 2, 1, 1);
book(_hVis_nJet60_abs , 3, 1, 1);
book(_hVis_nJet60 , 4, 1, 1);
book(_hVis_nJet100_abs , 5, 1, 1);
book(_hVis_nJet100 , 6, 1, 1);
book(_hVis_addJet1Pt_abs , 7, 1, 1);
book(_hVis_addJet1Pt , 8, 1, 1);
book(_hVis_addJet1Eta_abs , 9, 1, 1);
book(_hVis_addJet1Eta ,10, 1, 1);
book(_hVis_addJet2Pt_abs ,11, 1, 1);
book(_hVis_addJet2Pt ,12, 1, 1);
book(_hVis_addJet2Eta_abs ,13, 1, 1);
book(_hVis_addJet2Eta ,14, 1, 1);
book(_hVis_addJJMass_abs ,15, 1, 1);
book(_hVis_addJJMass ,16, 1, 1);
book(_hVis_addJJDR_abs ,17, 1, 1);
book(_hVis_addJJDR ,18, 1, 1);
book(_hVis_addJJHT_abs ,19, 1, 1);
book(_hVis_addJJHT ,20, 1, 1);
book(_hFull_addJet1Pt_abs ,21, 1, 1);
book(_hFull_addJet1Pt ,22, 1, 1);
book(_hFull_addJet1Eta_abs ,23, 1, 1);
book(_hFull_addJet1Eta ,24, 1, 1);
book(_hFull_addJet2Pt_abs ,25, 1, 1);
book(_hFull_addJet2Pt ,26, 1, 1);
book(_hFull_addJet2Eta_abs ,27, 1, 1);
book(_hFull_addJet2Eta ,28, 1, 1);
book(_hFull_addJJMass_abs ,29, 1, 1);
book(_hFull_addJJMass ,30, 1, 1);
book(_hFull_addJJDR_abs ,31, 1, 1);
book(_hFull_addJJDR ,32, 1, 1);
book(_hFull_addJJHT_abs ,33, 1, 1);
book(_hFull_addJJHT ,34, 1, 1);
book(_hVis_addBJet1Pt_abs ,35, 1, 1);
book(_hVis_addBJet1Pt ,36, 1, 1);
book(_hVis_addBJet1Eta_abs ,37, 1, 1);
book(_hVis_addBJet1Eta ,38, 1, 1);
book(_hVis_addBJet2Pt_abs ,39, 1, 1);
book(_hVis_addBJet2Pt ,40, 1, 1);
book(_hVis_addBJet2Eta_abs ,41, 1, 1);
book(_hVis_addBJet2Eta ,42, 1, 1);
book(_hVis_addBBMass_abs ,43, 1, 1);
book(_hVis_addBBMass ,44, 1, 1);
book(_hVis_addBBDR_abs ,45, 1, 1);
book(_hVis_addBBDR ,46, 1, 1);
book(_hFull_addBJet1Pt_abs ,47, 1, 1);
book(_hFull_addBJet1Pt ,48, 1, 1);
book(_hFull_addBJet1Eta_abs ,49, 1, 1);
book(_hFull_addBJet1Eta ,50, 1, 1);
book(_hFull_addBJet2Pt_abs ,51, 1, 1);
book(_hFull_addBJet2Pt ,52, 1, 1);
book(_hFull_addBJet2Eta_abs ,53, 1, 1);
book(_hFull_addBJet2Eta ,54, 1, 1);
book(_hFull_addBBMass_abs ,55, 1, 1);
book(_hFull_addBBMass ,56, 1, 1);
book(_hFull_addBBDR_abs ,57, 1, 1);
book(_hFull_addBBDR ,58, 1, 1);
book(_h_gap_addJet1Pt ,59, 1, 1);
book(_h_gap_addJet1Pt_eta0 ,60, 1, 1);
book(_h_gap_addJet1Pt_eta1 ,61, 1, 1);
book(_h_gap_addJet1Pt_eta2 ,62, 1, 1);
book(_h_gap_addJet2Pt ,63, 1, 1);
book(_h_gap_addJet2Pt_eta0 ,64, 1, 1);
book(_h_gap_addJet2Pt_eta1 ,65, 1, 1);
book(_h_gap_addJet2Pt_eta2 ,66, 1, 1);
book(_h_gap_addJetHT ,67, 1, 1);
book(_h_gap_addJetHT_eta0 ,68, 1, 1);
book(_h_gap_addJetHT_eta1 ,69, 1, 1);
book(_h_gap_addJetHT_eta2 ,70, 1, 1);
}
void analyze(const Event& event) {
// The objects used in the PAPER 12-041 are defined as follows (see p.16 for details):
//
// * Leptons : from the W boson decays after FSR
// * Jets : anti-kT R=0.5 to all stable particles
// exclude W->enu, munu, taunu
// * B jet : B-Ghost matched
// * B from top : B hadron from top->b decay
//
// Visible phase space definition:
//
// * Leptons : pT > 20, |eta| < 2.4
// * B jets from top : pT > 30, |eta| < 2.4
// Additional jets : pT > 20, |eta| < 2.4
// *
// Full phase space definition:
//
// * Correction to dilepton BR from W boson BR
// * No cut on top decay products
// * Additional jets : pT > 20, |eta| < 2.4
// Do the analysis only for the ttbar full leptonic channel, removing tau decays
const Particles partontops = apply<ParticleFinder>(event, "PartonTops").particlesByPt();
if (partontops.size() != 2) vetoEvent;
const Particle& t1 = partontops[0];
const Particle& t2 = partontops[1];
// Apply acceptance cuts on top-decay leptons (existence should be guaranteed)
- const auto isPromptChLepton = [](const Particle& p){return isChargedLepton(p) && !fromDecay(p);};
+ const auto isPromptChLepton = [](const Particle& p){return p.isPrompt() && isChargedLepton(p);};
const Particle lep1 = t1.allDescendants(lastParticleWith(isPromptChLepton)).front();
const Particle lep2 = t2.allDescendants(lastParticleWith(isPromptChLepton)).front();
if (lep1.pT() < 1e-9*GeV || lep2.pT() < 1e-9*GeV) vetoEvent; // sanity check?
const Jets jets = apply<JetAlg>(event, "Jets").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 2.4);
int nJet30 = 0, nJet60 = 0, nJet100 = 0;
Jets topBJets, addJets, addBJets, addJets_eta0, addJets_eta1, addJets_eta2;
for (const Jet& jet : jets) {
if (jet.pT() > 30*GeV) nJet30 += 1;
if (jet.pT() > 60*GeV) nJet60 += 1;
if (jet.pT() > 100*GeV) nJet100 += 1;
const bool isBtagged = jet.bTagged();
const bool isBFromTop = any(jet.bTags(), hasParticleAncestorWith(Cuts::abspid == PID::TQUARK));
if (isBFromTop) {
if (jet.pT() > 30*GeV) topBJets.push_back(jet);
} else {
addJets.push_back(jet);
if (isBtagged) addBJets.push_back(jet);
if (jet.abseta() < 0.8 ) addJets_eta0.push_back(jet);
else if (jet.abseta() < 1.5 ) addJets_eta1.push_back(jet);
else if (jet.abseta() < 2.4 ) addJets_eta2.push_back(jet);
}
}
const bool isVisiblePS = topBJets.size() >= 2
&& lep1.pT() > 20*GeV && lep1.abseta() < 2.4 && lep2.pT() > 20*GeV && lep2.abseta() < 2.4;
MSG_DEBUG(isVisiblePS << ": #b(top) = " << topBJets.size()
<< "; l1 = " << lep1.pT() << ", " << lep1.abseta()
<< "; l2 = " << lep2.pT() << ", " << lep2.abseta());
const double weight = 1.0;
if (isVisiblePS) {
fillWithOF(_hVis_nJet30_abs, nJet30, weight);
fillWithOF(_hVis_nJet30, nJet30, weight);
fillWithOF(_hVis_nJet60_abs, nJet60, weight);
fillWithOF(_hVis_nJet60, nJet60, weight);
fillWithOF(_hVis_nJet100_abs, nJet100, weight);
fillWithOF(_hVis_nJet100, nJet100, weight);
fillGapFractions(addJets, _h_gap_addJet1Pt, _h_gap_addJet2Pt, _h_gap_addJetHT, weight);
fillGapFractions(addJets_eta0, _h_gap_addJet1Pt_eta0, _h_gap_addJet2Pt_eta0, _h_gap_addJetHT_eta0, weight);
fillGapFractions(addJets_eta1, _h_gap_addJet1Pt_eta1, _h_gap_addJet2Pt_eta1, _h_gap_addJetHT_eta1, weight);
fillGapFractions(addJets_eta2, _h_gap_addJet1Pt_eta2, _h_gap_addJet2Pt_eta2, _h_gap_addJetHT_eta2, weight);
}
// Plots with two additional jets
if (addJets.size() >= 1) {
const double ht = sum(addJets, pT, 0.0);
_hFull_addJJHT_abs->fill(ht/GeV, weight);
_hFull_addJJHT ->fill(ht/GeV, weight);
if (isVisiblePS) {
_hVis_addJJHT_abs->fill(ht/GeV, weight);
_hVis_addJJHT ->fill(ht/GeV, weight);
}
const Jet& j1 = addJets[0];
_hFull_addJet1Pt_abs ->fill(j1.pT()/GeV, weight);
_hFull_addJet1Pt ->fill(j1.pT()/GeV, weight);
_hFull_addJet1Eta_abs->fill(j1.abseta(), weight);
_hFull_addJet1Eta ->fill(j1.abseta(), weight);
if (isVisiblePS) {
_hVis_addJet1Pt_abs ->fill(j1.pT()/GeV, weight);
_hVis_addJet1Pt ->fill(j1.pT()/GeV, weight);
_hVis_addJet1Eta_abs->fill(j1.abseta(), weight);
_hVis_addJet1Eta ->fill(j1.abseta(), weight);
}
if (addJets.size() >= 2) {
const Jet& j2 = addJets[1];
_hFull_addJet2Pt_abs ->fill(j2.pT()/GeV, weight);
_hFull_addJet2Pt ->fill(j2.pT()/GeV, weight);
_hFull_addJet2Eta_abs->fill(j2.abseta(), weight);
_hFull_addJet2Eta ->fill(j2.abseta(), weight);
if (isVisiblePS) {
_hVis_addJet2Pt_abs ->fill(j2.pT()/GeV, weight);
_hVis_addJet2Pt ->fill(j2.pT()/GeV, weight);
_hVis_addJet2Eta_abs->fill(j2.abseta(), weight);
_hVis_addJet2Eta ->fill(j2.abseta(), weight);
}
const double jjmass = (j1.mom() + j2.mom()).mass();
const double jjdR = deltaR(j1, j2);
_hFull_addJJMass_abs->fill(jjmass/GeV, weight);
_hFull_addJJMass ->fill(jjmass/GeV, weight);
_hFull_addJJDR_abs ->fill(jjdR, weight);
_hFull_addJJDR ->fill(jjdR, weight);
if (isVisiblePS) {
_hVis_addJJMass_abs->fill(jjmass/GeV, weight);
_hVis_addJJMass ->fill(jjmass/GeV, weight);
_hVis_addJJDR_abs ->fill(jjdR, weight);
_hVis_addJJDR ->fill(jjdR, weight);
}
}
}
// Same set of plots if there are additional b-jets
if (addBJets.size() >= 1) {
const Jet& b1 = addBJets[0];
_hFull_addBJet1Pt_abs ->fill(b1.pT()/GeV, weight);
_hFull_addBJet1Pt ->fill(b1.pT()/GeV, weight);
_hFull_addBJet1Eta_abs->fill(b1.abseta(), weight);
_hFull_addBJet1Eta ->fill(b1.abseta(), weight);
if (isVisiblePS) {
_hVis_addBJet1Pt_abs ->fill(b1.pT()/GeV, weight);
_hVis_addBJet1Pt ->fill(b1.pT()/GeV, weight);
_hVis_addBJet1Eta_abs->fill(b1.abseta(), weight);
_hVis_addBJet1Eta ->fill(b1.abseta(), weight);
}
if (addBJets.size() >= 2) {
const Jet& b2 = addBJets[1];
_hFull_addBJet2Pt_abs ->fill(b2.pT()/GeV, weight);
_hFull_addBJet2Pt ->fill(b2.pT()/GeV, weight);
_hFull_addBJet2Eta_abs->fill(b2.abseta(), weight);
_hFull_addBJet2Eta ->fill(b2.abseta(), weight);
if (isVisiblePS) {
_hVis_addBJet2Pt_abs ->fill(b2.pT()/GeV, weight);
_hVis_addBJet2Pt ->fill(b2.pT()/GeV, weight);
_hVis_addBJet2Eta_abs->fill(b2.abseta(), weight);
_hVis_addBJet2Eta ->fill(b2.abseta(), weight);
}
const double bbmass = (b1.mom() + b2.mom()).mass();
const double bbdR = deltaR(b1, b2);
_hFull_addBBMass_abs->fill(bbmass/GeV, weight);
_hFull_addBBMass ->fill(bbmass/GeV, weight);
_hFull_addBBDR_abs ->fill(bbdR, weight);
_hFull_addBBDR ->fill(bbdR, weight);
if (isVisiblePS) {
_hVis_addBBMass_abs->fill(bbmass/GeV, weight);
_hVis_addBBMass ->fill(bbmass/GeV, weight);
_hVis_addBBDR_abs ->fill(bbdR, weight);
_hVis_addBBDR ->fill(bbdR, weight);
}
}
}
}
void finalize() {
const double ttbarXS = !std::isnan(crossSectionPerEvent()) ? crossSection() : 252.89*picobarn;
if (std::isnan(crossSectionPerEvent()))
MSG_INFO("No valid cross-section given, using NNLO (arXiv:1303.6254; sqrt(s)=8 TeV, m_t=172.5 GeV): " << ttbarXS/picobarn << " pb");
normalize({_hVis_nJet30,_hVis_nJet60, _hVis_nJet100,
_hVis_addJet1Pt, _hVis_addJet1Eta, _hVis_addJet2Pt, _hVis_addJet2Eta,
_hVis_addJJMass, _hVis_addJJDR, _hVis_addJJHT,
_hFull_addJet1Pt, _hFull_addJet1Eta, _hFull_addJet2Pt, _hFull_addJet2Eta,
_hFull_addJJMass, _hFull_addJJDR, _hFull_addJJHT,
_hVis_addBJet1Pt, _hVis_addBJet1Eta, _hVis_addBJet2Pt, _hVis_addBJet2Eta,
_hVis_addBBMass, _hVis_addBBDR,
_hFull_addBJet1Pt, _hFull_addBJet1Eta, _hFull_addBJet2Pt, _hFull_addBJet2Eta,
_hFull_addBBMass, _hFull_addBBDR});
const double xsPerWeight = ttbarXS/picobarn / sumOfWeights();
scale({_hVis_nJet30_abs, _hVis_nJet60_abs, _hVis_nJet100_abs,
_hVis_addJet1Pt_abs, _hVis_addJet1Eta_abs, _hVis_addJet2Pt_abs, _hVis_addJet2Eta_abs,
_hVis_addJJMass_abs, _hVis_addJJDR_abs, _hVis_addJJHT_abs,
_hVis_addBJet1Pt_abs, _hVis_addBJet1Eta_abs, _hVis_addBJet2Pt_abs, _hVis_addBJet2Eta_abs,
_hVis_addBBMass_abs, _hVis_addBBDR_abs}, xsPerWeight);
const double sfull = xsPerWeight / 0.0454; //< correct for dilepton branching fraction
scale({_hFull_addJet1Pt_abs, _hFull_addJet1Eta_abs, _hFull_addJet2Pt_abs, _hFull_addJet2Eta_abs,
_hFull_addJJMass_abs, _hFull_addJJDR_abs, _hFull_addJJHT_abs,
_hFull_addBJet1Pt_abs, _hFull_addBJet1Eta_abs, _hFull_addBJet2Pt_abs, _hFull_addBJet2Eta_abs,
_hFull_addBBMass_abs, _hFull_addBBDR_abs}, sfull);
}
//@}
void fillWithOF(Histo1DPtr h, double x, double w) {
h->fill(std::min(x, h->xMax()-1e-9), w);
}
void fillGapFractions(const Jets& addJets, Profile1DPtr h_gap_addJet1Pt, Profile1DPtr h_gap_addJet2Pt, Profile1DPtr h_gap_addJetHT, double weight) {
const double j1pt = (addJets.size() > 0) ? addJets[0].pT() : 0;
for (size_t i = 0; i < h_gap_addJet1Pt->numBins(); ++i) {
const double binCenter = h_gap_addJet1Pt->bin(i).xMid();
h_gap_addJet1Pt->fillBin(i, int(j1pt/GeV < binCenter), weight);
}
const double j2pt = (addJets.size() > 1) ? addJets[1].pT() : 0;
for (size_t i = 0; i < h_gap_addJet2Pt->numBins(); ++i) {
const double binCenter = h_gap_addJet2Pt->bin(i).xMid();
h_gap_addJet2Pt->fillBin(i, int(j2pt/GeV < binCenter), weight);
}
const double ht = sum(addJets, pT, 0.);
for (size_t i = 0; i < h_gap_addJetHT->numBins(); ++i) {
const double binCenter = h_gap_addJetHT->bin(i).xMid();
h_gap_addJetHT->fillBin(i, int(ht/GeV < binCenter) , weight);
}
}
// @name Histogram data members
//@{
Histo1DPtr _hVis_nJet30_abs, _hVis_nJet60_abs, _hVis_nJet100_abs;
Histo1DPtr _hVis_addJet1Pt_abs, _hVis_addJet1Eta_abs, _hVis_addJet2Pt_abs, _hVis_addJet2Eta_abs;
Histo1DPtr _hVis_addJJMass_abs, _hVis_addJJDR_abs, _hVis_addJJHT_abs;
Histo1DPtr _hFull_addJet1Pt_abs, _hFull_addJet1Eta_abs, _hFull_addJet2Pt_abs, _hFull_addJet2Eta_abs;
Histo1DPtr _hFull_addJJMass_abs, _hFull_addJJDR_abs, _hFull_addJJHT_abs;
Histo1DPtr _hVis_addBJet1Pt_abs, _hVis_addBJet1Eta_abs, _hVis_addBJet2Pt_abs, _hVis_addBJet2Eta_abs;
Histo1DPtr _hVis_addBBMass_abs, _hVis_addBBDR_abs;
Histo1DPtr _hFull_addBJet1Pt_abs, _hFull_addBJet1Eta_abs, _hFull_addBJet2Pt_abs, _hFull_addBJet2Eta_abs;
Histo1DPtr _hFull_addBBMass_abs, _hFull_addBBDR_abs;
Histo1DPtr _hVis_nJet30, _hVis_nJet60, _hVis_nJet100;
Histo1DPtr _hVis_addJet1Pt, _hVis_addJet1Eta, _hVis_addJet2Pt, _hVis_addJet2Eta;
Histo1DPtr _hVis_addJJMass, _hVis_addJJDR, _hVis_addJJHT;
Histo1DPtr _hFull_addJet1Pt, _hFull_addJet1Eta, _hFull_addJet2Pt, _hFull_addJet2Eta;
Histo1DPtr _hFull_addJJMass, _hFull_addJJDR, _hFull_addJJHT;
Histo1DPtr _hVis_addBJet1Pt, _hVis_addBJet1Eta, _hVis_addBJet2Pt, _hVis_addBJet2Eta;
Histo1DPtr _hVis_addBBMass, _hVis_addBBDR;
Histo1DPtr _hFull_addBJet1Pt, _hFull_addBJet1Eta, _hFull_addBJet2Pt, _hFull_addBJet2Eta;
Histo1DPtr _hFull_addBBMass, _hFull_addBBDR;
Profile1DPtr _h_gap_addJet1Pt, _h_gap_addJet1Pt_eta0, _h_gap_addJet1Pt_eta1, _h_gap_addJet1Pt_eta2;
Profile1DPtr _h_gap_addJet2Pt, _h_gap_addJet2Pt_eta0, _h_gap_addJet2Pt_eta1, _h_gap_addJet2Pt_eta2;
Profile1DPtr _h_gap_addJetHT, _h_gap_addJetHT_eta0, _h_gap_addJetHT_eta1, _h_gap_addJetHT_eta2;
//@}
};
// The hook for the plugin system
DECLARE_RIVET_PLUGIN(CMS_2015_I1397174);
}
diff --git a/analyses/pluginMC/MC_DILEPTON.cc b/analyses/pluginMC/MC_DILEPTON.cc
--- a/analyses/pluginMC/MC_DILEPTON.cc
+++ b/analyses/pluginMC/MC_DILEPTON.cc
@@ -1,134 +1,132 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/PromptFinalState.hh"
namespace Rivet {
/// @brief Study of prompt dilepton kinematics sensitive to boson polarizations
class MC_DILEPTON : public Analysis {
public:
/// Constructor
DEFAULT_RIVET_ANALYSIS_CTOR(MC_DILEPTON);
/// @name Analysis methods
//@{
/// Book histograms and initialise projections before the run
void init() {
// Initialise and register projections
declare(PromptFinalState((Cuts::abspid == PID::ELECTRON || Cuts::abspid == PID::MUON)
&& Cuts::abseta < 5 && Cuts::pT > 10*GeV), "Leptons");
// Book histograms
- _h_pt_l1 = bookHisto1D("lep1_pt", logspace(40, 10, 400));
- _h_costheta_l1 = bookHisto1D("lep1_costheta", linspace(25, -1, 1));
- _h_ppara_l1 = bookHisto1D("lep1_ppara", linspace(40, -50, 350));
- _h_pperp_l1 = bookHisto1D("lep1_pperp", linspace(25, 0, 100));
+ book(_h_pt_l1, "lep1_pt", logspace(40, 10, 400));
+ book(_h_costheta_l1, "lep1_costheta", linspace(25, -1, 1));
+ book(_h_ppara_l1, "lep1_ppara", linspace(40, -50, 350));
+ book(_h_pperp_l1, "lep1_pperp", linspace(25, 0, 100));
//
- _h_pt_l2 = bookHisto1D("lep2_pt", logspace(40, 10, 400));
- _h_costheta_l2 = bookHisto1D("lep2_costheta", linspace(25, -1, 1));
- _h_ppara_l2 = bookHisto1D("lep2_ppara", linspace(40, -50, 350));
- _h_pperp_l2 = bookHisto1D("lep2_pperp", linspace(25, 0, 100));
+ book(_h_pt_l2, "lep2_pt", logspace(40, 10, 400));
+ book(_h_costheta_l2, "lep2_costheta", linspace(25, -1, 1));
+ book(_h_ppara_l2, "lep2_ppara", linspace(40, -50, 350));
+ book(_h_pperp_l2, "lep2_pperp", linspace(25, 0, 100));
//
- _h_costheta_com_l1 = bookHisto1D("com_costheta_l1", linspace(25, -1, 1));
- _h_costheta_com_l2 = bookHisto1D("com_costheta_l2", linspace(25, -1, 1));
- _h_ppara_com_l1 = bookHisto1D("com_ppara_l1", linspace(25, -50, 50));
- _h_ppara_com_l2 = bookHisto1D("com_ppara_l2", linspace(25, -50, 50));
+ book(_h_costheta_com_l1, "com_costheta_l1", linspace(25, -1, 1));
+ book(_h_costheta_com_l2, "com_costheta_l2", linspace(25, -1, 1));
+ book(_h_ppara_com_l1, "com_ppara_l1", linspace(25, -50, 50));
+ book(_h_ppara_com_l2, "com_ppara_l2", linspace(25, -50, 50));
//
- _h_costheta_com = bookHisto1D("com_costheta", linspace(25, -1, 1));
- _h_ppara_com = bookHisto1D("com_ppara", linspace(25, -50, 50));
- _h_pperp_com = bookHisto1D("com_pperp", linspace(25, 0, 100));
+ book(_h_costheta_com, "com_costheta", linspace(25, -1, 1));
+ book(_h_ppara_com, "com_ppara", linspace(25, -50, 50));
+ book(_h_pperp_com, "com_pperp", linspace(25, 0, 100));
}
/// Perform the per-event analysis
void analyze(const Event& event) {
- const double weight = event.weight();
-
const Particles& leptons = apply<FinalState>(event, "Leptons").particlesByPt();
if (leptons.size() != 2) vetoEvent;
const Particle& l1 = leptons[0];
const Particle& l2 = leptons[1];
- _h_pt_l1->fill(l1.pT()/GeV, weight);
- _h_pt_l2->fill(l2.pT()/GeV, weight);
+ _h_pt_l1->fill(l1.pT()/GeV);
+ _h_pt_l2->fill(l2.pT()/GeV);
const FourMomentum pcom = l1.mom() + l2.mom();
const Vector3 betacom = pcom.betaVec();
const Vector3 unitboostvec = betacom.unit();
const LorentzTransform comboost = LorentzTransform::mkFrameTransformFromBeta(betacom);
const double l1_costheta = cos(l1.p3().angle(unitboostvec));
const double l1_ppara = l2.p3().dot(unitboostvec);
const double l1_pperp = l2.p3().cross(unitboostvec).mod();
- _h_costheta_l1->fill(l1_costheta, weight);
- _h_ppara_l1->fill(l1_ppara, weight);
- _h_pperp_l1->fill(l1_pperp, weight);
+ _h_costheta_l1->fill(l1_costheta);
+ _h_ppara_l1->fill(l1_ppara);
+ _h_pperp_l1->fill(l1_pperp);
const double l2_costheta = cos(l2.p3().angle(unitboostvec));
const double l2_ppara = l2.p3().dot(unitboostvec);
const double l2_pperp = l2.p3().cross(unitboostvec).mod();
- _h_costheta_l2->fill(l2_costheta, weight);
- _h_ppara_l2->fill(l2_ppara, weight);
- _h_pperp_l2->fill(l2_pperp, weight);
+ _h_costheta_l2->fill(l2_costheta);
+ _h_ppara_l2->fill(l2_ppara);
+ _h_pperp_l2->fill(l2_pperp);
const FourMomentum p1com = comboost.transform(l1.mom());
const FourMomentum p2com = comboost.transform(l2.mom());
const double com_costheta1 = cos(p1com.p3().angle(unitboostvec));
const double com_costheta2 = cos(p2com.p3().angle(unitboostvec));
MSG_DEBUG("CoM cos(th)s: " << com_costheta1 << ", " << com_costheta2);
// assert(com_costheta1 == com_costheta2 && "CoM cos(th)s differ");
const double com_ppara1 = p1com.p3().dot(unitboostvec);
const double com_ppara2 = p2com.p3().dot(unitboostvec);
MSG_DEBUG("CoM p_paras: " << com_ppara1 << ", " << com_ppara2);
// assert(com_ppara1 == com_ppara2 && "CoM p_paras differ");
const double com_pperp1 = p1com.p3().cross(unitboostvec).mod();
const double com_pperp2 = p2com.p3().cross(unitboostvec).mod();
MSG_DEBUG("CoM p_pperps: " << com_pperp1 << ", " << com_pperp2);
// assert(com_pperp1 == com_pperp2 && "CoM p_perps differ");
- _h_costheta_com_l1->fill(com_costheta1, weight);
- _h_costheta_com_l2->fill(com_costheta2, weight);
- _h_costheta_com->fill(com_costheta1, 0.5*weight);
- _h_costheta_com->fill(com_costheta2, 0.5*weight);
- _h_ppara_com_l1->fill(com_ppara1, weight);
- _h_ppara_com_l2->fill(com_ppara2, weight);
- _h_ppara_com->fill(com_ppara1, 0.5*weight);
- _h_ppara_com->fill(com_ppara2, 0.5*weight);
- _h_pperp_com->fill(com_pperp1, weight);
+ _h_costheta_com_l1->fill(com_costheta1);
+ _h_costheta_com_l2->fill(com_costheta2);
+ _h_costheta_com->fill(com_costheta1, 0.5);
+ _h_costheta_com->fill(com_costheta2, 0.5);
+ _h_ppara_com_l1->fill(com_ppara1);
+ _h_ppara_com_l2->fill(com_ppara2);
+ _h_ppara_com->fill(com_ppara1, 0.5);
+ _h_ppara_com->fill(com_ppara2, 0.5);
+ _h_pperp_com->fill(com_pperp1);
}
/// Normalise histograms etc., after the run
void finalize() {
normalize({_h_pt_l1, _h_pt_l2});
normalize({_h_ppara_com, _h_pperp_com, _h_costheta_com});
normalize({_h_ppara_com_l1, _h_ppara_com_l2, _h_costheta_com_l1, _h_costheta_com_l2});
normalize({_h_ppara_l1, _h_pperp_l1, _h_costheta_l1});
normalize({_h_ppara_l2, _h_pperp_l2, _h_costheta_l2});
}
//@}
/// @name Histograms
//@{
Histo1DPtr _h_pt_l1, _h_pt_l2;
Histo1DPtr _h_ppara_com, _h_pperp_com, _h_costheta_com;
Histo1DPtr _h_ppara_com_l1, _h_ppara_com_l2, _h_costheta_com_l1, _h_costheta_com_l2;
Histo1DPtr _h_ppara_l1, _h_pperp_l1, _h_costheta_l1;
Histo1DPtr _h_ppara_l2, _h_pperp_l2, _h_costheta_l2;
//@}
};
// The hook for the plugin system
DECLARE_RIVET_PLUGIN(MC_DILEPTON);
}
diff --git a/analyses/pluginPetra/JADE_OPAL_2000_S4300807.cc b/analyses/pluginPetra/JADE_OPAL_2000_S4300807.cc
--- a/analyses/pluginPetra/JADE_OPAL_2000_S4300807.cc
+++ b/analyses/pluginPetra/JADE_OPAL_2000_S4300807.cc
@@ -1,186 +1,177 @@
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/FinalState.hh"
namespace Rivet {
/// @brief Jet rates in \f$ e^+ e^- \f$ at OPAL and JADE
/// @author Frank Siegert
class JADE_OPAL_2000_S4300807 : public Analysis {
public:
- /// @name Constructors etc.
- //@{
-
/// Constructor
- JADE_OPAL_2000_S4300807()
- : Analysis("JADE_OPAL_2000_S4300807")
- { }
-
- //@}
+ DEFAULT_RIVET_ANALYSIS_CTOR(JADE_OPAL_2000_S4300807)
/// @name Analysis methods
//@{
void init() {
// Projections
const FinalState fs;
declare(fs, "FS");
- FastJets jadeJets = FastJets(fs, FastJets::JADE, 0.7);
- FastJets durhamJets = FastJets(fs, FastJets::DURHAM, 0.7);
- jadeJets.useInvisibles(true);
- durhamJets.useInvisibles(true);
+ FastJets jadeJets = FastJets(fs, FastJets::JADE, 0.7, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY);
+ FastJets durhamJets = FastJets(fs, FastJets::DURHAM, 0.7, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY);
declare(jadeJets, "JadeJets");
declare(durhamJets, "DurhamJets");
// Histos
int offset = 0;
switch (int(sqrtS()/GeV + 0.5)) {
case 35: offset = 7; break;
case 44: offset = 8; break;
case 91: offset = 9; break;
case 133: offset = 10; break;
case 161: offset = 11; break;
case 172: offset = 12; break;
case 183: offset = 13; break;
case 189: offset = 14; break;
default: break;
}
for (size_t i = 0; i < 5; ++i) {
book(_h_R_Jade[i] ,offset, 1, i+1);
book(_h_R_Durham[i] ,offset+9, 1, i+1);
if (i < 4) book(_h_y_Durham[i] ,offset+17, 1, i+1);
}
}
void analyze(const Event& e) {
MSG_DEBUG("Num particles = " << apply<FinalState>(e, "FS").particles().size());
const FastJets& jadejet = apply<FastJets>(e, "JadeJets");
if (jadejet.clusterSeq()) {
const double y_23 = jadejet.clusterSeq()->exclusive_ymerge_max(2);
const double y_34 = jadejet.clusterSeq()->exclusive_ymerge_max(3);
const double y_45 = jadejet.clusterSeq()->exclusive_ymerge_max(4);
const double y_56 = jadejet.clusterSeq()->exclusive_ymerge_max(5);
for (size_t i = 0; i < _h_R_Jade[0]->numBins(); ++i) {
double ycut = _h_R_Jade[0]->bin(i).xMid();
double width = _h_R_Jade[0]->bin(i).xWidth();
if (y_23 < ycut) {
_h_R_Jade[0]->fillBin(i, width);
}
}
for (size_t i = 0; i < _h_R_Jade[1]->numBins(); ++i) {
double ycut = _h_R_Jade[1]->bin(i).xMid();
double width = _h_R_Jade[1]->bin(i).xWidth();
if (y_34 < ycut && y_23 > ycut) {
_h_R_Jade[1]->fillBin(i, width);
}
}
for (size_t i = 0; i < _h_R_Jade[2]->numBins(); ++i) {
double ycut = _h_R_Jade[2]->bin(i).xMid();
double width = _h_R_Jade[2]->bin(i).xWidth();
if (y_45 < ycut && y_34 > ycut) {
_h_R_Jade[2]->fillBin(i, width);
}
}
for (size_t i = 0; i < _h_R_Jade[3]->numBins(); ++i) {
double ycut = _h_R_Jade[3]->bin(i).xMid();
double width = _h_R_Jade[3]->bin(i).xWidth();
if (y_56 < ycut && y_45 > ycut) {
_h_R_Jade[3]->fillBin(i, width);
}
}
for (size_t i = 0; i < _h_R_Jade[4]->numBins(); ++i) {
double ycut = _h_R_Jade[4]->bin(i).xMid();
double width = _h_R_Jade[4]->bin(i).xWidth();
if (y_56 > ycut) {
_h_R_Jade[4]->fillBin(i, width);
}
}
}
const FastJets& durjet = apply<FastJets>(e, "DurhamJets");
if (durjet.clusterSeq()) {
const double y_23 = durjet.clusterSeq()->exclusive_ymerge_max(2);
const double y_34 = durjet.clusterSeq()->exclusive_ymerge_max(3);
const double y_45 = durjet.clusterSeq()->exclusive_ymerge_max(4);
const double y_56 = durjet.clusterSeq()->exclusive_ymerge_max(5);
_h_y_Durham[0]->fill(y_23);
_h_y_Durham[1]->fill(y_34);
_h_y_Durham[2]->fill(y_45);
_h_y_Durham[3]->fill(y_56);
for (size_t i = 0; i < _h_R_Durham[0]->numBins(); ++i) {
double ycut = _h_R_Durham[0]->bin(i).xMid();
double width = _h_R_Durham[0]->bin(i).xWidth();
if (y_23 < ycut) {
_h_R_Durham[0]->fillBin(i, width);
}
}
for (size_t i = 0; i < _h_R_Durham[1]->numBins(); ++i) {
double ycut = _h_R_Durham[1]->bin(i).xMid();
double width = _h_R_Durham[1]->bin(i).xWidth();
if (y_34 < ycut && y_23 > ycut) {
_h_R_Durham[1]->fillBin(i, width);
}
}
for (size_t i = 0; i < _h_R_Durham[2]->numBins(); ++i) {
double ycut = _h_R_Durham[2]->bin(i).xMid();
double width = _h_R_Durham[2]->bin(i).xWidth();
if (y_45 < ycut && y_34 > ycut) {
_h_R_Durham[2]->fillBin(i, width);
}
}
for (size_t i = 0; i < _h_R_Durham[3]->numBins(); ++i) {
double ycut = _h_R_Durham[3]->bin(i).xMid();
double width = _h_R_Durham[3]->bin(i).xWidth();
if (y_56 < ycut && y_45 > ycut) {
_h_R_Durham[3]->fillBin(i, width);
}
}
for (size_t i = 0; i < _h_R_Durham[4]->numBins(); ++i) {
double ycut = _h_R_Durham[4]->bin(i).xMid();
double width = _h_R_Durham[4]->bin(i).xWidth();
if (y_56 > ycut) {
_h_R_Durham[4]->fillBin(i, width);
}
}
}
}
/// Finalize
void finalize() {
for (size_t n = 0; n < 4; ++n) normalize(_h_y_Durham[n]);
for (size_t n = 0; n < 5; ++n) scale(_h_R_Jade[n], 100/sumOfWeights());
for (size_t n = 0; n < 5; ++n) scale(_h_R_Durham[n], 100/sumOfWeights());
}
//@}
private:
/// @name Histograms
//@{
Histo1DPtr _h_R_Jade[5];
Histo1DPtr _h_R_Durham[5];
Histo1DPtr _h_y_Durham[4];
//@}
};
// The hook for the plugin system
DECLARE_RIVET_PLUGIN(JADE_OPAL_2000_S4300807);
}
diff --git a/include/Rivet/Particle.hh b/include/Rivet/Particle.hh
--- a/include/Rivet/Particle.hh
+++ b/include/Rivet/Particle.hh
@@ -1,700 +1,700 @@
// -*- C++ -*-
#ifndef RIVET_Particle_HH
#define RIVET_Particle_HH
#include "Rivet/Particle.fhh"
#include "Rivet/ParticleBase.hh"
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/Tools/Cuts.hh"
#include "Rivet/Tools/Utils.hh"
#include "Rivet/Math/LorentzTrans.hh"
// NOTE: Rivet/Tools/ParticleUtils.hh included at the end
#include "fastjet/PseudoJet.hh"
namespace Rivet {
/// Particle representation, either from a HepMC::GenEvent or reconstructed.
class Particle : public ParticleBase {
public:
/// @name Constructors
//@{
/// Default constructor.
/// @note A particle without info is useless. This only exists to keep STL containers happy.
Particle()
: ParticleBase(),
_original(nullptr), _id(PID::ANY)
{ }
/// Constructor without GenParticle.
Particle(PdgId pid, const FourMomentum& mom, const FourVector& pos=FourVector())
: ParticleBase(),
_original(nullptr), _id(pid), _momentum(mom), _origin(pos)
{ }
/// Constructor from a HepMC GenParticle pointer.
Particle(const GenParticle* gp)
: ParticleBase(),
_original(gp), _id(gp->pdg_id()),
_momentum(gp->momentum())
{
const GenVertex* vprod = gp->production_vertex();
if (vprod != nullptr) {
setOrigin(vprod->position().t(), vprod->position().x(), vprod->position().y(), vprod->position().z());
}
}
/// Constructor from a HepMC GenParticle.
Particle(const GenParticle& gp)
: ParticleBase(),
_original(&gp), _id(gp.pdg_id()),
_momentum(gp.momentum())
{
const GenVertex* vprod = gp.production_vertex();
if (vprod != nullptr) {
setOrigin(vprod->position().t(), vprod->position().x(), vprod->position().y(), vprod->position().z());
}
}
//@}
/// @name Kinematic properties
//@{
/// The momentum.
const FourMomentum& momentum() const {
return _momentum;
}
/// Set the momentum.
Particle& setMomentum(const FourMomentum& momentum) {
_momentum = momentum;
return *this;
}
/// Set the momentum via components.
Particle& setMomentum(double E, double px, double py, double pz) {
_momentum = FourMomentum(E, px, py, pz);
return *this;
}
/// Apply an active Lorentz transform to this particle
Particle& transformBy(const LorentzTransform& lt);
//@
/// @name Positional properties
//@{
/// The origin position.
const FourVector& origin() const {
return _origin;
}
/// Set the origin position.
Particle& setOrigin(const FourVector& position) {
_origin = position;
return *this;
}
/// Set the origin position via components.
Particle& setOrigin(double t, double x, double y, double z) {
_origin = FourMomentum(t, x, y, z);
return *this;
}
//@}
/// @name Other representations and implicit casts to momentum-like objects
//@{
/// Converter to FastJet3 PseudoJet
virtual fastjet::PseudoJet pseudojet() const {
return fastjet::PseudoJet(mom().px(), mom().py(), mom().pz(), mom().E());
}
/// Cast operator to FastJet3 PseudoJet
operator PseudoJet () const { return pseudojet(); }
/// Get a const pointer to the original GenParticle
const GenParticle* genParticle() const {
return _original;
}
/// Cast operator for conversion to GenParticle*
operator const GenParticle* () const { return genParticle(); }
//@}
/// @name Particle ID code accessors
//@{
/// This Particle's PDG ID code.
PdgId pid() const { return _id; }
/// Absolute value of the PDG ID code.
PdgId abspid() const { return std::abs(_id); }
//@}
/// @name Charge
//@{
/// The charge of this Particle.
double charge() const { return PID::charge(pid()); }
/// The absolute charge of this Particle.
double abscharge() const { return PID::abscharge(pid()); }
/// Three times the charge of this Particle (i.e. integer multiple of smallest quark charge).
int charge3() const { return PID::charge3(pid()); }
/// Three times the absolute charge of this Particle (i.e. integer multiple of smallest quark charge).
int abscharge3() const { return PID::abscharge3(pid()); }
/// Is this Particle charged?
bool isCharged() const { return charge3() != 0; }
//@}
/// @name Particle species
//@{
/// Is this a hadron?
bool isHadron() const { return PID::isHadron(pid()); }
/// Is this a meson?
bool isMeson() const { return PID::isMeson(pid()); }
/// Is this a baryon?
bool isBaryon() const { return PID::isBaryon(pid()); }
/// Is this a lepton?
bool isLepton() const { return PID::isLepton(pid()); }
/// Is this a charged lepton?
bool isChargedLepton() const { return PID::isChargedLepton(pid()); }
/// Is this a neutrino?
bool isNeutrino() const { return PID::isNeutrino(pid()); }
/// Does this (hadron) contain a b quark?
bool hasBottom() const { return PID::hasBottom(pid()); }
/// Does this (hadron) contain a c quark?
bool hasCharm() const { return PID::hasCharm(pid()); }
// /// Does this (hadron) contain an s quark?
// bool hasStrange() const { return PID::hasStrange(pid()); }
/// Is this particle potentially visible in a detector?
bool isVisible() const;
//@}
/// @name Constituents (for composite particles)
//@{
/// Set direct constituents of this particle
virtual void setConstituents(const Particles& cs, bool setmom=false);
/// Add a single direct constituent to this particle
virtual void addConstituent(const Particle& c, bool addmom=false);
/// Add direct constituents to this particle
virtual void addConstituents(const Particles& cs, bool addmom=false);
/// Determine if this Particle is a composite of other Rivet Particles
bool isComposite() const { return !constituents().empty(); }
/// @brief Direct constituents of this particle, returned by reference
///
/// The returned vector will be empty if this particle is non-composite,
/// and its entries may themselves be composites.
const Particles& constituents() const { return _constituents; }
/// @brief Direct constituents of this particle, sorted by a functor
/// @note Returns a copy, thanks to the sorting
const Particles constituents(const ParticleSorter& sorter) const {
return sortBy(constituents(), sorter);
}
/// @brief Direct constituents of this particle, filtered by a Cut
/// @note Returns a copy, thanks to the filtering
const Particles constituents(const Cut& c) const {
return filter_select(constituents(), c);
}
/// @brief Direct constituents of this particle, sorted by a functor
/// @note Returns a copy, thanks to the filtering and sorting
const Particles constituents(const Cut& c, const ParticleSorter& sorter) const {
return sortBy(constituents(c), sorter);
}
/// @brief Direct constituents of this particle, filtered by a selection functor
/// @note Returns a copy, thanks to the filtering
const Particles constituents(const ParticleSelector& selector) const {
return filter_select(constituents(), selector);
}
/// @brief Direct constituents of this particle, filtered and sorted by functors
/// @note Returns a copy, thanks to the filtering and sorting
const Particles constituents(const ParticleSelector& selector, const ParticleSorter& sorter) const {
return sortBy(constituents(selector), sorter);
}
/// @brief Fundamental constituents of this particle
/// @note Returns {{*this}} if this particle is non-composite.
Particles rawConstituents() const;
/// @brief Fundamental constituents of this particle, sorted by a functor
/// @note Returns a copy, thanks to the sorting
const Particles rawConstituents(const ParticleSorter& sorter) const {
return sortBy(rawConstituents(), sorter);
}
/// @brief Fundamental constituents of this particle, filtered by a Cut
/// @note Returns a copy, thanks to the filtering
const Particles rawConstituents(const Cut& c) const {
return filter_select(rawConstituents(), c);
}
/// @brief Fundamental constituents of this particle, sorted by a functor
/// @note Returns a copy, thanks to the filtering and sorting
const Particles rawConstituents(const Cut& c, const ParticleSorter& sorter) const {
return sortBy(rawConstituents(c), sorter);
}
/// @brief Fundamental constituents of this particle, filtered by a selection functor
/// @note Returns a copy, thanks to the filtering
const Particles rawConstituents(const ParticleSelector& selector) const {
return filter_select(rawConstituents(), selector);
}
/// @brief Fundamental constituents of this particle, filtered and sorted by functors
/// @note Returns a copy, thanks to the filtering and sorting
const Particles rawConstituents(const ParticleSelector& selector, const ParticleSorter& sorter) const {
return sortBy(rawConstituents(selector), sorter);
}
//@}
/// @name Ancestry (for fundamental particles with a HepMC link)
//@{
/// Get a list of the direct parents of the current particle (with optional selection Cut)
///
/// @note This is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
Particles parents(const Cut& c=Cuts::OPEN) const;
/// Get a list of the direct parents of the current particle (with selector function)
///
/// @note This is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
Particles parents(const ParticleSelector& f) const {
return filter_select(parents(), f);
}
/// Check whether any particle in the particle's parent list has the requested property
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool hasParentWith(const ParticleSelector& f) const {
return !parents(f).empty();
}
/// Check whether any particle in the particle's parent list has the requested property
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool hasParentWith(const Cut& c) const;
/// Check whether any particle in the particle's parent list does not have the requested property
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool hasParentWithout(const ParticleSelector& f) const {
return hasParentWith([&](const Particle& p){ return !f(p); });
}
/// Check whether any particle in the particle's parent list does not have the requested property
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool hasParentWithout(const Cut& c) const;
/// Check whether a given PID is found in the particle's parent list
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
///
/// @deprecated Prefer e.g. hasParentWith(Cut::pid == 123)
bool hasParent(PdgId pid) const;
/// Get a list of the ancestors of the current particle (with optional selection Cut)
///
/// @note By default only physical ancestors, with status=2, are returned.
///
/// @note This is valid in MC, but may not be answerable experimentally --
/// use this function with care when replicating experimental analyses!
Particles ancestors(const Cut& c=Cuts::OPEN, bool only_physical=true) const;
/// Get a list of the direct parents of the current particle (with selector function)
///
/// @note By default only physical ancestors, with status=2, are returned.
///
/// @note This is valid in MC, but may not be answerable experimentally --
/// use this function with care when replicating experimental analyses!
Particles ancestors(const ParticleSelector& f, bool only_physical=true) const {
return filter_select(ancestors(Cuts::OPEN, only_physical), f);
}
/// Check whether any particle in the particle's ancestor list has the requested property
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool hasAncestorWith(const ParticleSelector& f, bool only_physical=true) const {
return !ancestors(f, only_physical).empty();
}
/// Check whether any particle in the particle's ancestor list has the requested property
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool hasAncestorWith(const Cut& c, bool only_physical=true) const;
/// Check whether any particle in the particle's ancestor list does not have the requested property
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool hasAncestorWithout(const ParticleSelector& f, bool only_physical=true) const {
return hasAncestorWith([&](const Particle& p){ return !f(p); }, only_physical);
}
/// Check whether any particle in the particle's ancestor list does not have the requested property
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool hasAncestorWithout(const Cut& c, bool only_physical=true) const;
/// Check whether a given PID is found in the particle's ancestor list
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
///
/// @deprecated Prefer hasAncestorWith(Cuts::pid == pid) etc.
bool hasAncestor(PdgId pid, bool only_physical=true) const;
/// @brief Determine whether the particle is from a b-hadron decay
///
/// @note This question is valid in MC, but may not be perfectly answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool fromBottom() const;
/// @brief Determine whether the particle is from a c-hadron decay
///
/// @note If a hadron contains b and c quarks it is considered a bottom
/// hadron and NOT a charm hadron.
///
/// @note This question is valid in MC, but may not be perfectly answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool fromCharm() const;
// /// @brief Determine whether the particle is from a s-hadron decay
// ///
// /// @note If a hadron contains b or c quarks as well as strange it is
// /// considered a b or c hadron, but NOT a strange hadron.
// ///
// /// @note This question is valid in MC, but may not be perfectly answerable
// /// experimentally -- use this function with care when replicating
// /// experimental analyses!
// bool fromStrange() const;
/// @brief Determine whether the particle is from a hadron decay
///
/// @note This question is valid in MC, but may not be perfectly answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool fromHadron() const;
/// @brief Determine whether the particle is from a tau decay
///
/// @note This question is valid in MC, but may not be perfectly answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool fromTau(bool prompt_taus_only=false) const;
/// @brief Determine whether the particle is from a prompt tau decay
///
/// @note This question is valid in MC, but may not be perfectly answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool fromPromptTau() const { return fromTau(true); }
/// @brief Determine whether the particle is from a tau which decayed hadronically
///
/// @note This question is valid in MC, but may not be perfectly answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool fromHadronicTau(bool prompt_taus_only=false) const;
/// @brief Determine whether the particle is from a hadron or tau decay
///
/// Specifically, walk up the ancestor chain until a status 2 hadron or
/// tau is found, if at all.
///
/// @note This question is valid in MC, but may not be perfectly answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
///
- DEPRECATED("Too vague: use fromHadron or fromHadronicTau")
+ DEPRECATED("Too vague: use fromHadron() || fromPromptTau(), or isDirect()")
bool fromDecay() const { return fromHadron() || fromPromptTau(); }
/// @brief Shorthand definition of 'promptness' based on set definition flags
///
/// A "direct" particle is one directly connected to the hard process. It is a
/// preferred alias for "prompt", since it has no confusing implications about
/// distinguishability by timing information.
///
/// The boolean arguments allow a decay lepton to be considered direct if
/// its parent was a "real" direct lepton.
///
/// @note This one doesn't make any judgements about final-stateness
bool isDirect(bool allow_from_direct_tau=false, bool allow_from_direct_mu=false) const;
/// Alias for isDirect
bool isPrompt(bool allow_from_prompt_tau=false, bool allow_from_prompt_mu=false) const {
return isDirect(allow_from_prompt_tau, allow_from_prompt_mu);
}
//@}
/// @name Decay info
//@{
/// Whether this particle is stable according to the generator
bool isStable() const;
/// @todo isDecayed? How to restrict to physical particles?
/// Get a list of the direct descendants from the current particle (with optional selection Cut)
Particles children(const Cut& c=Cuts::OPEN) const;
/// Get a list of the direct descendants from the current particle (with selector function)
Particles children(const ParticleSelector& f) const {
return filter_select(children(), f);
}
/// Check whether any direct child of this particle has the requested property
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool hasChildWith(const ParticleSelector& f) const {
return !children(f).empty();
}
/// Check whether any direct child of this particle has the requested property
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool hasChildWith(const Cut& c) const;
/// Check whether any direct child of this particle does not have the requested property
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool hasChildWithout(const ParticleSelector& f) const {
return hasChildWith([&](const Particle& p){ return !f(p); });
}
/// Check whether any direct child of this particle does not have the requested property
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool hasChildWithout(const Cut& c) const;
/// Get a list of all the descendants from the current particle (with optional selection Cut)
Particles allDescendants(const Cut& c=Cuts::OPEN, bool remove_duplicates=true) const;
/// Get a list of all the descendants from the current particle (with selector function)
Particles allDescendants(const ParticleSelector& f, bool remove_duplicates=true) const {
return filter_select(allDescendants(Cuts::OPEN, remove_duplicates), f);
}
/// Check whether any descendant of this particle has the requested property
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool hasDescendantWith(const ParticleSelector& f, bool remove_duplicates=true) const {
return !allDescendants(f, remove_duplicates).empty();
}
/// Check whether any descendant of this particle has the requested property
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool hasDescendantWith(const Cut& c, bool remove_duplicates=true) const;
/// Check whether any descendant of this particle does not have the requested property
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool hasDescendantWithout(const ParticleSelector& f, bool remove_duplicates=true) const {
return hasDescendantWith([&](const Particle& p){ return !f(p); }, remove_duplicates);
}
/// Check whether any descendant of this particle does not have the requested property
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool hasDescendantWithout(const Cut& c, bool remove_duplicates=true) const;
/// Get a list of all the stable descendants from the current particle (with optional selection Cut)
///
/// @todo Use recursion through replica-avoiding MCUtils functions to avoid bookkeeping duplicates
/// @todo Insist that the current particle is post-hadronization, otherwise throw an exception?
Particles stableDescendants(const Cut& c=Cuts::OPEN) const;
/// Get a list of all the stable descendants from the current particle (with selector function)
Particles stableDescendants(const ParticleSelector& f) const {
return filter_select(stableDescendants(), f);
}
/// Check whether any stable descendant of this particle has the requested property
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool hasStableDescendantWith(const ParticleSelector& f) const {
return !stableDescendants(f).empty();
}
/// Check whether any stable descendant of this particle has the requested property
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool hasStableDescendantWith(const Cut& c) const;
/// Check whether any stable descendant of this particle does not have the requested property
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool hasStableDescendantWithout(const ParticleSelector& f) const {
return hasStableDescendantWith([&](const Particle& p){ return !f(p); });
}
/// Check whether any stable descendant of this particle does not have the requested property
///
/// @note This question is valid in MC, but may not be answerable
/// experimentally -- use this function with care when replicating
/// experimental analyses!
bool hasStableDescendantWithout(const Cut& c) const;
/// Flight length (divide by mm or cm to get the appropriate units)
double flightLength() const;
//@}
/// @name Duplicate testing
//@{
/// @brief Determine whether a particle is the first in a decay chain to meet the function requirement
inline bool isFirstWith(const ParticleSelector& f) const {
if (!f(*this)) return false; //< This doesn't even meet f, let alone being the last to do so
if (any(parents(), f)) return false; //< If a direct parent has this property, this isn't the first
return true;
}
/// @brief Determine whether a particle is the first in a decay chain not to meet the function requirement
inline bool isFirstWithout(const ParticleSelector& f) const {
return isFirstWith([&](const Particle& p){ return !f(p); });
}
/// @brief Determine whether a particle is the last in a decay chain to meet the function requirement
inline bool isLastWith(const ParticleSelector& f) const {
if (!f(*this)) return false; //< This doesn't even meet f, let alone being the last to do so
if (any(children(), f)) return false; //< If a child has this property, this isn't the last
return true;
}
/// @brief Determine whether a particle is the last in a decay chain not to meet the function requirement
inline bool isLastWithout(const ParticleSelector& f) const {
return isLastWith([&](const Particle& p){ return !f(p); });
}
//@}
protected:
/// A pointer to the original GenParticle from which this Particle is projected (may be null)
const GenParticle* _original;
/// Constituent particles if this is a composite (may be empty)
Particles _constituents;
/// The PDG ID code for this Particle.
PdgId _id;
/// The momentum of this particle.
FourMomentum _momentum;
/// The creation position of this particle.
FourVector _origin;
};
/// @name String representation and streaming support
//@{
/// Allow a Particle to be passed to an ostream.
std::ostream& operator << (std::ostream& os, const Particle& p);
/// Allow ParticlePair to be passed to an ostream.
std::ostream& operator << (std::ostream& os, const ParticlePair& pp);
//@}
}
#include "Rivet/Tools/ParticleUtils.hh"
#endif
diff --git a/include/Rivet/Tools/ParticleUtils.hh b/include/Rivet/Tools/ParticleUtils.hh
--- a/include/Rivet/Tools/ParticleUtils.hh
+++ b/include/Rivet/Tools/ParticleUtils.hh
@@ -1,768 +1,768 @@
#ifndef RIVET_PARTICLEUTILS_HH
#define RIVET_PARTICLEUTILS_HH
#include "Rivet/Particle.hh"
#include "Rivet/Tools/ParticleBaseUtils.hh"
#include "Rivet/Tools/ParticleIdUtils.hh"
// Macros to map Rivet::Particle functions to PID:: functions of the same name
#define PARTICLE_TO_PID_BOOLFN(fname) inline bool fname (const Particle& p) { return PID:: fname (p.pid()); }
#define PARTICLE_TO_PID_INTFN(fname) inline int fname (const Particle& p) { return PID:: fname (p.pid()); }
#define PARTICLE_TO_PID_DBLFN(fname) inline double fname (const Particle& p) { return PID:: fname (p.pid()); }
namespace Rivet {
/// @name Particle classifier functions
//@{
/// Unbound function access to PID code
inline int pid(const Particle& p) { return p.pid(); }
/// Unbound function access to abs PID code
inline int abspid(const Particle& p) { return p.abspid(); }
/// Is this particle species charged?
PARTICLE_TO_PID_BOOLFN(isCharged)
/// Is this particle species neutral?
PARTICLE_TO_PID_BOOLFN(isNeutral)
/// Is this a neutrino?
PARTICLE_TO_PID_BOOLFN(isNeutrino)
/// Determine if the PID is that of a charged lepton
PARTICLE_TO_PID_BOOLFN(isChargedLepton)
PARTICLE_TO_PID_BOOLFN(isChLepton)
/// Determine if the PID is that of a lepton (charged or neutral)
PARTICLE_TO_PID_BOOLFN(isLepton)
/// Determine if the PID is that of a photon
PARTICLE_TO_PID_BOOLFN(isPhoton)
/// Determine if the PID is that of an electron or positron
PARTICLE_TO_PID_BOOLFN(isElectron)
/// Determine if the PID is that of an muon or antimuon
PARTICLE_TO_PID_BOOLFN(isMuon)
/// Determine if the PID is that of an tau or antitau
PARTICLE_TO_PID_BOOLFN(isTau)
/// Determine if the PID is that of a hadron
PARTICLE_TO_PID_BOOLFN(isHadron)
/// Determine if the PID is that of a meson
PARTICLE_TO_PID_BOOLFN(isMeson)
/// Determine if the PID is that of a baryon
PARTICLE_TO_PID_BOOLFN(isBaryon)
/// Determine if the PID is that of a quark
PARTICLE_TO_PID_BOOLFN(isQuark)
/// Determine if the PID is that of a parton (quark or gluon)
PARTICLE_TO_PID_BOOLFN(isParton)
/// Determine if the PID is that of a W+
PARTICLE_TO_PID_BOOLFN(isWplus)
/// Determine if the PID is that of a W-
PARTICLE_TO_PID_BOOLFN(isWminus)
/// Determine if the PID is that of a W+-
PARTICLE_TO_PID_BOOLFN(isW)
/// Determine if the PID is that of a Z0
PARTICLE_TO_PID_BOOLFN(isZ)
/// Determine if the PID is that of an SM/lightest SUSY Higgs
PARTICLE_TO_PID_BOOLFN(isHiggs)
/// Determine if the PID is that of an s/sbar
PARTICLE_TO_PID_BOOLFN(isStrange)
/// Determine if the PID is that of a c/cbar
PARTICLE_TO_PID_BOOLFN(isCharm)
/// Determine if the PID is that of a b/bbar
PARTICLE_TO_PID_BOOLFN(isBottom)
/// Determine if the PID is that of a t/tbar
PARTICLE_TO_PID_BOOLFN(isTop)
/// Determine if the particle is a heavy flavour hadron or parton
PARTICLE_TO_PID_BOOLFN(isHeavyFlavour)
/// Determine if the PID is that of a heavy parton (c,b,t)
PARTICLE_TO_PID_BOOLFN(isHeavyParton)
/// Determine if the PID is that of a light parton (u,d,s)
PARTICLE_TO_PID_BOOLFN(isLightParton)
/// Determine if the PID is that of a heavy flavour (b or c) meson
PARTICLE_TO_PID_BOOLFN(isHeavyMeson)
/// Determine if the PID is that of a heavy flavour (b or c) baryon
PARTICLE_TO_PID_BOOLFN(isHeavyBaryon)
/// Determine if the PID is that of a heavy flavour (b or c) hadron
PARTICLE_TO_PID_BOOLFN(isHeavyHadron)
/// Determine if the PID is that of a light flavour (not b or c) meson
PARTICLE_TO_PID_BOOLFN(isLightMeson)
/// Determine if the PID is that of a light flavour (not b or c) baryon
PARTICLE_TO_PID_BOOLFN(isLightBaryon)
/// Determine if the PID is that of a light flavour (not b or c) hadron
PARTICLE_TO_PID_BOOLFN(isLightHadron)
/// Determine if the PID is that of a b-meson.
PARTICLE_TO_PID_BOOLFN(isBottomMeson)
/// Determine if the PID is that of a b-baryon.
PARTICLE_TO_PID_BOOLFN(isBottomBaryon)
/// Determine if the PID is that of a b-hadron.
PARTICLE_TO_PID_BOOLFN(isBottomHadron)
/// @brief Determine if the PID is that of a c-meson.
///
/// Specifically, the _heaviest_ quark is a c: a B_c is a b-meson and NOT a c-meson.
/// Charmonia (closed charm) are counted as c-mesons here.
PARTICLE_TO_PID_BOOLFN(isCharmMeson)
/// @brief Determine if the PID is that of a c-baryon.
///
/// Specifically, the _heaviest_ quark is a c: a baryon containing a b & c
/// is a b-baryon and NOT a c-baryon. To test for the simpler case, just use
/// a combination of hasCharm() and isBaryon().
PARTICLE_TO_PID_BOOLFN(isCharmBaryon)
/// Determine if the PID is that of a c-hadron.
PARTICLE_TO_PID_BOOLFN(isCharmHadron)
// /// Determine if the PID is that of a strange meson
// PARTICLE_TO_PID_BOOLFN(isStrangeMeson)
// /// Determine if the PID is that of a strange baryon
// PARTICLE_TO_PID_BOOLFN(isStrangeBaryon)
// /// Determine if the PID is that of a strange hadron
// PARTICLE_TO_PID_BOOLFN(isStrangeHadron)
/// Is this a pomeron, odderon, or generic reggeon?
PARTICLE_TO_PID_BOOLFN(isReggeon)
/// Determine if the PID is that of a diquark (used in hadronization models)
PARTICLE_TO_PID_BOOLFN(isDiquark)
/// Determine if the PID is that of a pentaquark (hypothetical hadron)
PARTICLE_TO_PID_BOOLFN(isPentaquark)
/// Is this a fundamental SUSY particle?
PARTICLE_TO_PID_BOOLFN(isSUSY)
/// Is this an R-hadron?
PARTICLE_TO_PID_BOOLFN(isRhadron)
/// Is this a technicolor particle?
PARTICLE_TO_PID_BOOLFN(isTechnicolor)
/// Is this an excited (composite) quark or lepton?
PARTICLE_TO_PID_BOOLFN(isExcited)
/// Is this a Kaluza-Klein excitation?
PARTICLE_TO_PID_BOOLFN(isKK)
/// Is this a graviton?
PARTICLE_TO_PID_BOOLFN(isGraviton)
/// Is this a BSM particle (including graviton)?
PARTICLE_TO_PID_BOOLFN(isBSM)
/// Determine if the PID is in the generator-specific range
PARTICLE_TO_PID_BOOLFN(isGenSpecific)
/// Determine if the PID is that of an EW scale resonance
PARTICLE_TO_PID_BOOLFN(isResonance)
/// Check the PID for usability in transport codes like Geant4
PARTICLE_TO_PID_BOOLFN(isTransportable)
/// Does this particle contain an up quark?
PARTICLE_TO_PID_BOOLFN(hasUp)
/// Does this particle contain a down quark?
PARTICLE_TO_PID_BOOLFN(hasDown)
/// Does this particle contain a strange quark?
PARTICLE_TO_PID_BOOLFN(hasStrange)
/// Does this particle contain a charm quark?
PARTICLE_TO_PID_BOOLFN(hasCharm)
/// Does this particle contain a bottom quark?
PARTICLE_TO_PID_BOOLFN(hasBottom)
/// Does this particle contain a top quark?
PARTICLE_TO_PID_BOOLFN(hasTop)
/// jSpin returns 2J+1, where J is the total spin
PARTICLE_TO_PID_INTFN(jSpin)
/// sSpin returns 2S+1, where S is the spin
PARTICLE_TO_PID_INTFN(sSpin)
/// lSpin returns 2L+1, where L is the orbital angular momentum
PARTICLE_TO_PID_INTFN(lSpin)
/// Return the charge
PARTICLE_TO_PID_DBLFN(charge)
/// Return 3 times the charge (3 x quark charge is an int)
PARTICLE_TO_PID_INTFN(charge3)
/// Return the absolute charge
PARTICLE_TO_PID_DBLFN(abscharge)
/// Return 3 times the abs charge (3 x quark charge is an int)
PARTICLE_TO_PID_INTFN(abscharge3)
/// Get the atomic number (number of protons) in a nucleus/ion
PARTICLE_TO_PID_INTFN(nuclZ)
/// Get the atomic weight (number of nucleons) in a nucleus/ion
PARTICLE_TO_PID_INTFN(nuclA)
/// If this is a nucleus (ion), get nLambda
PARTICLE_TO_PID_INTFN(nuclNlambda)
//@}
/// @name Particle pair classifiers
/// @todo Make versions that work on ParticlePair?
//@{
inline bool isSameSign(const Particle& a, const Particle& b) { return PID::isSameSign(a.pid(), b.pid()); }
inline bool isOppSign(const Particle& a, const Particle& b) { return PID::isOppSign(a.pid(), b.pid()); }
inline bool isSameFlav(const Particle& a, const Particle& b) { return PID::isSameFlav(a.pid(), b.pid()); }
inline bool isOppFlav(const Particle& a, const Particle& b) { return PID::isOppFlav(a.pid(), b.pid()); }
inline bool isOSSF(const Particle& a, const Particle& b) { return PID::isOSSF(a.pid(), b.pid()); }
inline bool isSSSF(const Particle& a, const Particle& b) { return PID::isSSSF(a.pid(), b.pid()); }
inline bool isOSOF(const Particle& a, const Particle& b) { return PID::isOSOF(a.pid(), b.pid()); }
inline bool isSSOF(const Particle& a, const Particle& b) { return PID::isSSOF(a.pid(), b.pid()); }
//@}
/// @name Particle charge/sign comparison functions
//@{
/// @brief Return true if Particles @a a and @a b have the opposite charge sign
/// @note Two neutrals returns false
inline bool oppSign(const Particle& a, const Particle& b) {
return sign(a.charge3()) == -sign(b.charge3()) && sign(a.charge3()) != ZERO;
}
/// Return true if Particles @a a and @a b have the same charge sign
/// @note Two neutrals returns true
inline bool sameSign(const Particle& a, const Particle& b) {
return sign(a.charge3()) == sign(b.charge3());
}
/// Return true if Particles @a a and @a b have the exactly opposite charge
/// @note Two neutrals returns false
inline bool oppCharge(const Particle& a, const Particle& b) {
return a.charge3() == -b.charge3() && a.charge3() != 0;
}
/// Return true if Particles @a a and @a b have the same charge (including neutral)
/// @note Two neutrals returns true
inline bool sameCharge(const Particle& a, const Particle& b) {
return a.charge3() == b.charge3();
}
/// Return true if Particles @a a and @a b have a different (not necessarily opposite) charge
inline bool diffCharge(const Particle& a, const Particle& b) {
return a.charge3() != b.charge3();
}
//@}
//////////////////////////////////////
/// @name Non-PID particle properties, via unbound functions
//@{
/// @brief Determine whether a particle is the first in a decay chain to meet the function requirement
inline bool isFirstWith(const Particle& p, const ParticleSelector& f) {
return p.isFirstWith(f);
}
/// @brief Determine whether a particle is the first in a decay chain not to meet the function requirement
inline bool isFirstWithout(const Particle& p, const ParticleSelector& f) {
return p.isFirstWithout(f);
}
/// @brief Determine whether a particle is the last in a decay chain to meet the function requirement
inline bool isLastWith(const Particle& p, const ParticleSelector& f) {
return p.isLastWith(f);
}
/// @brief Determine whether a particle is the last in a decay chain not to meet the function requirement
inline bool isLastWithout(const Particle& p, const ParticleSelector& f) {
return p.isLastWithout(f);
}
/// @brief Determine whether a particle has an ancestor which meets the function requirement
inline bool hasAncestorWith(const Particle& p, const ParticleSelector& f) {
return p.hasAncestorWith(f);
}
/// @brief Determine whether a particle has an ancestor which doesn't meet the function requirement
inline bool hasAncestorWithout(const Particle& p, const ParticleSelector& f) {
return p.hasAncestorWithout(f);
}
/// @brief Determine whether a particle has a parent which meets the function requirement
inline bool hasParentWith(const Particle& p, const ParticleSelector& f) {
return p.hasParentWith(f);
}
/// @brief Determine whether a particle has a parent which doesn't meet the function requirement
inline bool hasParentWithout(const Particle& p, const ParticleSelector& f) {
return p.hasParentWithout(f);
}
/// @brief Determine whether a particle has a child which meets the function requirement
inline bool hasChildWith(const Particle& p, const ParticleSelector& f) {
return p.hasChildWith(f);
}
/// @brief Determine whether a particle has a child which doesn't meet the function requirement
inline bool hasChildWithout(const Particle& p, const ParticleSelector& f) {
return p.hasChildWithout(f);
}
/// @brief Determine whether a particle has a descendant which meets the function requirement
inline bool hasDescendantWith(const Particle& p, const ParticleSelector& f) {
return p.hasDescendantWith(f);
// return !p.allDescendants(f).empty();
}
/// @brief Determine whether a particle has a descendant which doesn't meet the function requirement
inline bool hasDescendantWithout(const Particle& p, const ParticleSelector& f) {
return p.hasDescendantWithout(f);
}
/// @brief Determine whether a particle has a stable descendant which meets the function requirement
inline bool hasStableDescendantWith(const Particle& p, const ParticleSelector& f) {
return p.hasStableDescendantWith(f);
}
/// @brief Determine whether a particle has a stable descendant which doesn't meet the function requirement
inline bool hasStableDescendantWithout(const Particle& p, const ParticleSelector& f) {
return p.hasStableDescendantWithout(f);
}
/// Is this particle potentially visible in a detector?
inline bool isVisible(const Particle& p) { return p.isVisible(); }
/// @brief Decide if a given particle is direct, via Particle::isDirect()
///
/// A "direct" particle is one directly connected to the hard process. It is a
/// preferred alias for "prompt", since it has no confusing implications about
/// distinguishability by timing information.
///
/// The boolean arguments allow a decay lepton to be considered direct if
/// its parent was a "real" direct lepton.
inline bool isDirect(const Particle& p, bool allow_from_direct_tau=false, bool allow_from_direct_mu=false) {
return p.isDirect(allow_from_direct_tau, allow_from_direct_mu);
}
/// @brief Decide if a given particle is prompt, via Particle::isPrompt()
///
/// The boolean arguments allow a decay lepton to be considered prompt if
/// its parent was a "real" prompt lepton.
inline bool isPrompt(const Particle& p, bool allow_from_prompt_tau=false, bool allow_from_prompt_mu=false) {
return p.isPrompt(allow_from_prompt_tau, allow_from_prompt_mu);
}
/// Decide if a given particle is stable, via Particle::isStable()
inline bool isStable(const Particle& p) { return p.isStable(); }
/// Decide if a given particle decays hadronically
inline bool hasHadronicDecay(const Particle& p) {
if (p.isStable()) return false;
if (p.hasChildWith(isHadron)) return true;
return false;
}
/// Decide if a given particle decays leptonically (decays, and no hadrons)
inline bool hasLeptonicDecay(const Particle& p) {
if (p.isStable()) return false;
if (p.hasChildWith(isHadron)) return false;
return true;
}
/// Check whether a given PID is found in the particle's ancestor list
/// @deprecated Prefer hasAncestorWith
inline bool hasAncestor(const Particle& p, PdgId pid) { return p.hasAncestor(pid); }
/// Determine whether the particle is from a b-hadron decay
inline bool fromBottom(const Particle& p) { return p.fromBottom(); }
/// @brief Determine whether the particle is from a c-hadron decay
inline bool fromCharm(const Particle& p) { return p.fromCharm(); }
/// @brief Determine whether the particle is from a hadron decay
inline bool fromHadron(const Particle& p) { return p.fromHadron(); }
/// @brief Determine whether the particle is from a tau decay
inline bool fromTau(const Particle& p, bool prompt_taus_only=false) {
return p.fromTau(prompt_taus_only);
}
/// @brief Determine whether the particle is from a prompt tau decay
inline bool fromPromptTau(const Particle& p) { return p.fromPromptTau(); }
- /// @brief Determine whether the particle is from a hadron or tau decay
- DEPRECATED("Too vague: use fromHadron or fromHadronicTau")
- inline bool fromDecay(const Particle& p) { return p.fromDecay(); }
+ // /// @brief Determine whether the particle is from a hadron or tau decay
+ // DEPRECATED("Too vague: use fromHadron or fromHadronicTau")
+ // inline bool fromDecay(const Particle& p) { return p.fromDecay(); }
//@}
/// @name Particle classifier -> bool functors
///
/// To be passed to any() or all() e.g. any(p.children(), HasPID(PID::MUON))
//@{
/// Base type for Particle -> bool functors
struct BoolParticleFunctor {
virtual bool operator()(const Particle& p) const = 0;
virtual ~BoolParticleFunctor() {}
};
/// Functor for and-combination of selector logic
struct BoolParticleAND : public BoolParticleFunctor {
BoolParticleAND(const std::vector<ParticleSelector>& sels) : selectors(sels) {}
BoolParticleAND(const ParticleSelector& a, const ParticleSelector& b) : selectors({a,b}) {}
BoolParticleAND(const ParticleSelector& a, const ParticleSelector& b, const ParticleSelector& c) : selectors({a,b,c}) {}
bool operator()(const Particle& p) const {
for (const ParticleSelector& sel : selectors) if (!sel(p)) return false;
return true;
}
std::vector<ParticleSelector> selectors;
};
/// Operator syntactic sugar for AND construction
inline BoolParticleAND operator && (const ParticleSelector& a, const ParticleSelector& b) {
return BoolParticleAND(a, b);
}
/// Functor for or-combination of selector logic
struct BoolParticleOR : public BoolParticleFunctor {
BoolParticleOR(const std::vector<ParticleSelector>& sels) : selectors(sels) {}
BoolParticleOR(const ParticleSelector& a, const ParticleSelector& b) : selectors({a,b}) {}
BoolParticleOR(const ParticleSelector& a, const ParticleSelector& b, const ParticleSelector& c) : selectors({a,b,c}) {}
bool operator()(const Particle& p) const {
for (const ParticleSelector& sel : selectors) if (sel(p)) return true;
return false;
}
std::vector<ParticleSelector> selectors;
};
/// Operator syntactic sugar for OR construction
inline BoolParticleOR operator || (const ParticleSelector& a, const ParticleSelector& b) {
return BoolParticleOR(a, b);
}
/// Functor for inverting selector logic
struct BoolParticleNOT : public BoolParticleFunctor {
BoolParticleNOT(const ParticleSelector& sel) : selector(sel) {}
bool operator()(const Particle& p) const { return !selector(p); }
ParticleSelector selector;
};
/// Operator syntactic sugar for NOT construction
inline BoolParticleNOT operator ! (const ParticleSelector& a) {
return BoolParticleNOT(a);
}
/// PID matching functor
struct HasPID : public BoolParticleFunctor {
HasPID(PdgId pid) : targetpids{pid} { }
HasPID(vector<PdgId> pids) : targetpids{pids} { }
HasPID(initializer_list<PdgId> pids) : targetpids{pids} { }
bool operator()(const Particle& p) const { return contains(targetpids, p.pid()); }
vector<PdgId> targetpids;
};
using hasPID = HasPID;
/// |PID| matching functor
struct HasAbsPID : public BoolParticleFunctor {
HasAbsPID(PdgId pid) : targetapids{abs(pid)} { }
HasAbsPID(vector<PdgId> pids) { for (PdgId pid : pids) targetapids.push_back(abs(pid)); }
HasAbsPID(initializer_list<PdgId> pids) { for (PdgId pid : pids) targetapids.push_back(abs(pid)); }
bool operator()(const Particle& p) const { return contains(targetapids, p.abspid()); }
vector<PdgId> targetapids;
};
using hasAbsPID = HasAbsPID;
/// Determine whether a particle is the first in a decay chain to meet the cut/function
struct FirstParticleWith : public BoolParticleFunctor {
FirstParticleWith(const ParticleSelector& f) : fn(f) { }
FirstParticleWith(const Cut& c);
bool operator()(const Particle& p) const { return isFirstWith(p, fn); }
ParticleSelector fn;
};
using firstParticleWith = FirstParticleWith;
/// Determine whether a particle is the first in a decay chain not to meet the cut/function
struct FirstParticleWithout : public BoolParticleFunctor {
FirstParticleWithout(const ParticleSelector& f) : fn(f) { }
FirstParticleWithout(const Cut& c);
bool operator()(const Particle& p) const { return isFirstWithout(p, fn); }
ParticleSelector fn;
};
using firstParticleWithout = FirstParticleWithout;
/// Determine whether a particle is the last in a decay chain to meet the cut/function
struct LastParticleWith : public BoolParticleFunctor {
template <typename FN>
LastParticleWith(const FN& f) : fn(f) { }
LastParticleWith(const Cut& c);
bool operator()(const Particle& p) const { return isLastWith(p, fn); }
std::function<bool(const Particle&)> fn;
};
using lastParticleWith = LastParticleWith;
/// Determine whether a particle is the last in a decay chain not to meet the cut/function
struct LastParticleWithout : public BoolParticleFunctor {
LastParticleWithout(const ParticleSelector& f) : fn(f) { }
LastParticleWithout(const Cut& c);
bool operator()(const Particle& p) const { return isLastWithout(p, fn); }
ParticleSelector fn;
};
using lastParticleWithout = LastParticleWithout;
/// Determine whether a particle has an ancestor which meets the cut/function
struct HasParticleAncestorWith : public BoolParticleFunctor {
HasParticleAncestorWith(const ParticleSelector& f) : fn(f) { }
HasParticleAncestorWith(const Cut& c);
bool operator()(const Particle& p) const { return hasAncestorWith(p, fn); }
ParticleSelector fn;
};
using hasParticleAncestorWith = HasParticleAncestorWith;
/// Determine whether a particle has an ancestor which doesn't meet the cut/function
struct HasParticleAncestorWithout : public BoolParticleFunctor {
HasParticleAncestorWithout(const ParticleSelector& f) : fn(f) { }
HasParticleAncestorWithout(const Cut& c);
bool operator()(const Particle& p) const { return hasAncestorWithout(p, fn); }
ParticleSelector fn;
};
using hasParticleAncestorWithout = HasParticleAncestorWithout;
/// Determine whether a particle has an parent which meets the cut/function
struct HasParticleParentWith : public BoolParticleFunctor {
HasParticleParentWith(const ParticleSelector& f) : fn(f) { }
HasParticleParentWith(const Cut& c);
bool operator()(const Particle& p) const { return hasParentWith(p, fn); }
ParticleSelector fn;
};
using hasParticleParentWith = HasParticleParentWith;
/// Determine whether a particle has an parent which doesn't meet the cut/function
struct HasParticleParentWithout : public BoolParticleFunctor {
HasParticleParentWithout(const ParticleSelector& f) : fn(f) { }
HasParticleParentWithout(const Cut& c);
bool operator()(const Particle& p) const { return hasParentWithout(p, fn); }
ParticleSelector fn;
};
using hasParticleParentWithout = HasParticleParentWithout;
/// Determine whether a particle has a child which meets the cut/function
struct HasParticleChildWith : public BoolParticleFunctor {
HasParticleChildWith(const ParticleSelector& f) : fn(f) { }
HasParticleChildWith(const Cut& c);
bool operator()(const Particle& p) const { return hasChildWith(p, fn); }
ParticleSelector fn;
};
using hasParticleChildWith = HasParticleChildWith;
/// Determine whether a particle has a child which doesn't meet the cut/function
struct HasParticleChildWithout : public BoolParticleFunctor {
HasParticleChildWithout(const ParticleSelector& f) : fn(f) { }
HasParticleChildWithout(const Cut& c);
bool operator()(const Particle& p) const { return hasChildWithout(p, fn); }
ParticleSelector fn;
};
using hasParticleChildWithout = HasParticleChildWithout;
/// Determine whether a particle has a descendant which meets the cut/function
struct HasParticleDescendantWith : public BoolParticleFunctor {
HasParticleDescendantWith(const ParticleSelector& f) : fn(f) { }
HasParticleDescendantWith(const Cut& c);
bool operator()(const Particle& p) const { return hasDescendantWith(p, fn); }
ParticleSelector fn;
};
using hasParticleDescendantWith = HasParticleDescendantWith;
/// Determine whether a particle has a descendant which doesn't meet the cut/function
struct HasParticleDescendantWithout : public BoolParticleFunctor {
HasParticleDescendantWithout(const ParticleSelector& f) : fn(f) { }
HasParticleDescendantWithout(const Cut& c);
bool operator()(const Particle& p) const { return hasDescendantWithout(p, fn); }
ParticleSelector fn;
};
using hasParticleDescendantWithout = HasParticleDescendantWithout;
//@}
/// @name Unbound functions for filtering particles
//@{
/// Filter a particle collection in-place to the subset that passes the supplied Cut
Particles& ifilter_select(Particles& particles, const Cut& c);
/// Alias for ifilter_select
/// @deprecated Use ifilter_select
inline Particles& ifilterBy(Particles& particles, const Cut& c) { return ifilter_select(particles, c); }
/// Filter a particle collection in-place to the subset that passes the supplied Cut
inline Particles filter_select(const Particles& particles, const Cut& c) {
Particles rtn = particles;
return ifilter_select(rtn, c);
}
/// Alias for ifilter_select
/// @deprecated Use filter_select
inline Particles filterBy(const Particles& particles, const Cut& c) { return filter_select(particles, c); }
/// Filter a particle collection in-place to the subset that passes the supplied Cut
inline Particles filter_select(const Particles& particles, const Cut& c, Particles& out) {
out = filter_select(particles, c);
return out;
}
/// Alias for ifilter_select
/// @deprecated Use filter_select
inline Particles filterBy(const Particles& particles, const Cut& c, Particles& out) { return filter_select(particles, c, out); }
/// Filter a particle collection in-place to the subset that fails the supplied Cut
Particles& ifilter_discard(Particles& particles, const Cut& c);
/// Filter a particle collection in-place to the subset that fails the supplied Cut
inline Particles filter_discard(const Particles& particles, const Cut& c) {
Particles rtn = particles;
return ifilter_discard(rtn, c);
}
/// Filter a particle collection in-place to the subset that fails the supplied Cut
inline Particles filter_discard(const Particles& particles, const Cut& c, Particles& out) {
out = filter_discard(particles, c);
return out;
}
// inline void ifilterIsolateDeltaR(Particles& particles, const FourMomenta& vecs) {
// ifilter_discard(particles,
// }
// inline Particles filterIsolateDeltaR(const Particles& particles, const FourMomenta& vecs) {
// }
//@}
/// @name Particle pair functions
//@{
/// Get the PDG ID codes of a ParticlePair
/// @todo Make ParticlePair a custom class instead?
inline PdgIdPair pids(const ParticlePair& pp) {
return make_pair(pp.first.pid(), pp.second.pid());
}
//@}
/// @name Operations on collections of Particle
/// @note This can't be done on generic collections of ParticleBase -- thanks, C++ :-/
//@{
namespace Kin {
inline double sumPt(const Particles& ps) {
return sum(ps, pT, 0.0);
}
inline FourMomentum sumP4(const Particles& ps) {
return sum(ps, p4, FourMomentum());
}
inline Vector3 sumP3(const Particles& ps) {
return sum(ps, p3, Vector3());
}
/// @todo Min dPhi, min dR?
/// @todo Isolation routines?
}
//@}
// Import Kin namespace into Rivet
using namespace Kin;
}
#endif
diff --git a/src/AnalysisTools/MC_ParticleAnalysis.cc b/src/AnalysisTools/MC_ParticleAnalysis.cc
--- a/src/AnalysisTools/MC_ParticleAnalysis.cc
+++ b/src/AnalysisTools/MC_ParticleAnalysis.cc
@@ -1,168 +1,168 @@
// -*- C++ -*-
#include "Rivet/Analyses/MC_ParticleAnalysis.hh"
namespace Rivet {
-
+
MC_ParticleAnalysis::MC_ParticleAnalysis(const string& name,
size_t nparticles,
const string& particle_name)
: Analysis(name),
_nparts(nparticles), _pname(particle_name),
_h_pt(nparticles),
_h_eta(nparticles), _h_eta_plus(nparticles), _h_eta_minus(nparticles),
_h_rap(nparticles), _h_rap_plus(nparticles), _h_rap_minus(nparticles),
tmpeta(nparticles), tmprap(nparticles)
{
}
// Book histograms
void MC_ParticleAnalysis::init() {
for (size_t i = 0; i < _nparts; ++i) {
book(tmpeta[i], _pname + "_eta_pmratio_" + to_str(i+1));
book(tmprap[i], _pname + "_y_pmratio_" + to_str(i+1));
const string ptname = _pname + "_pt_" + to_str(i+1);
const double ptmax = 1.0/(double(i)+2.0) * (sqrtS()>0.?sqrtS():14000.)/GeV/2.0;
const int nbins_pt = 100/(i+1);
book(_h_pt[i] ,ptname, logspace(nbins_pt, 1.0, ptmax));
const string etaname = _pname + "_eta_" + to_str(i+1);
book(_h_eta[i] ,etaname, i > 1 ? 25 : 50, -5.0, 5.0);
book(_h_eta_plus[i], "_" + etaname + "_plus", i > 1 ? 15 : 25, 0, 5);
book(_h_eta_minus[i], "_" + etaname + "_minus", i > 1 ? 15 : 25, 0, 5);
const string rapname = _pname + "_y_" + to_str(i+1);
book(_h_rap[i] ,rapname, i > 1 ? 25 : 50, -5.0, 5.0);
book(_h_rap_plus[i], "_" + rapname + "_plus", i > 1 ? 15 : 25, 0, 5);
book(_h_rap_minus[i], "_" + rapname + "_minus", i > 1 ? 15 : 25, 0, 5);
for (size_t j = i+1; j < min(size_t(3), _nparts); ++j) {
const pair<size_t, size_t> ij = std::make_pair(i, j);
string detaname = _pname + "s_deta_" + to_str(i+1) + to_str(j+1);
Histo1DPtr tmpeta;
book(tmpeta, detaname, 25, -5.0, 5.0);
_h_deta.insert(make_pair(ij, tmpeta));
string dphiname = _pname + "s_dphi_" + to_str(i+1) + to_str(j+1);
Histo1DPtr tmpphi;
book(tmpphi, dphiname, 25, 0.0, M_PI);
_h_dphi.insert(make_pair(ij, tmpphi));
string dRname = _pname + "s_dR_" + to_str(i+1) + to_str(j+1);
Histo1DPtr tmpR;
book(tmpR, dRname, 25, 0.0, 5.0);
_h_dR.insert(make_pair(ij, tmpR));
}
}
book(_h_multi_exclusive ,_pname + "_multi_exclusive", _nparts+3, -0.5, _nparts+3-0.5);
book(_h_multi_inclusive ,_pname + "_multi_inclusive", _nparts+3, -0.5, _nparts+3-0.5);
book(_h_multi_ratio, _pname + "_multi_ratio");
book(_h_multi_exclusive_prompt ,_pname + "_multi_exclusive_prompt", _nparts+3, -0.5, _nparts+3-0.5);
book(_h_multi_inclusive_prompt ,_pname + "_multi_inclusive_prompt", _nparts+3, -0.5, _nparts+3-0.5);
book(_h_multi_ratio_prompt, _pname + "_multi_ratio_prompt");
}
// Do the analysis
void MC_ParticleAnalysis::_analyze(const Event& event, const Particles& particles) {
Particles promptparticles;
for (const Particle& p : particles)
- if (!p.fromDecay()) promptparticles += p;
+ if (p.isPrompt()) promptparticles += p;
for (size_t i = 0; i < _nparts; ++i) {
if (particles.size() < i+1) continue;
_h_pt[i]->fill(particles[i].pt()/GeV);
// Eta
const double eta_i = particles[i].eta();
_h_eta[i]->fill(eta_i);
(eta_i > 0.0 ? _h_eta_plus : _h_eta_minus)[i]->fill(fabs(eta_i));
// Rapidity
const double rap_i = particles[i].rapidity();
_h_rap[i]->fill(rap_i);
(rap_i > 0.0 ? _h_rap_plus : _h_rap_minus)[i]->fill(fabs(rap_i));
// Inter-particle properties
for (size_t j = i+1; j < min(size_t(3),_nparts); ++j) {
if (particles.size() < j+1) continue;
std::pair<size_t, size_t> ij = std::make_pair(i, j);
double deta = particles[i].eta() - particles[j].eta();
double dphi = deltaPhi(particles[i].momentum(), particles[j].momentum());
double dR = deltaR(particles[i].momentum(), particles[j].momentum());
_h_deta[ij]->fill(deta);
_h_dphi[ij]->fill(dphi);
_h_dR[ij]->fill(dR);
}
}
// Multiplicities
_h_multi_exclusive->fill(particles.size());
_h_multi_exclusive_prompt->fill(promptparticles.size());
for (size_t i = 0; i < _nparts+2; ++i) {
if (particles.size() >= i) _h_multi_inclusive->fill(i);
if (promptparticles.size() >= i) _h_multi_inclusive_prompt->fill(i);
}
}
// Finalize
void MC_ParticleAnalysis::finalize() {
const double scaling = crossSection()/sumOfWeights();
for (size_t i = 0; i < _nparts; ++i) {
scale(_h_pt[i], scaling);
scale(_h_eta[i], scaling);
scale(_h_rap[i], scaling);
// Create eta/rapidity ratio plots
divide(_h_eta_plus[i], _h_eta_minus[i], tmpeta[i]);
divide(_h_rap_plus[i], _h_rap_minus[i], tmprap[i]);
}
// Scale the d{eta,phi,R} histograms
typedef map<pair<size_t, size_t>, Histo1DPtr> HistMap;
for (HistMap::value_type& it : _h_deta) scale(it.second, scaling);
for (HistMap::value_type& it : _h_dphi) scale(it.second, scaling);
for (HistMap::value_type& it : _h_dR) scale(it.second, scaling);
// Fill inclusive multi ratios
for (size_t i = 0; i < _h_multi_inclusive->numBins()-1; ++i) {
_h_multi_ratio->addPoint(i+1, 0, 0.5, 0);
if (_h_multi_inclusive->bin(i).sumW() > 0.0) {
const double ratio = _h_multi_inclusive->bin(i+1).sumW() / _h_multi_inclusive->bin(i).sumW();
const double relerr_i = _h_multi_inclusive->bin(i).relErr();
const double relerr_j = _h_multi_inclusive->bin(i+1).relErr();
const double err = ratio * (relerr_i + relerr_j);
_h_multi_ratio->point(i).setY(ratio, err);
}
}
for (size_t i = 0; i < _h_multi_inclusive_prompt->numBins()-1; ++i) {
_h_multi_ratio_prompt->addPoint(i+1, 0, 0.5, 0);
if (_h_multi_inclusive_prompt->bin(i).sumW() > 0.0) {
const double ratio = _h_multi_inclusive_prompt->bin(i+1).sumW() / _h_multi_inclusive_prompt->bin(i).sumW();
const double relerr_i = _h_multi_inclusive_prompt->bin(i).relErr();
const double relerr_j = _h_multi_inclusive_prompt->bin(i+1).relErr();
const double err = ratio * (relerr_i + relerr_j);
_h_multi_ratio_prompt->point(i).setY(ratio, err);
}
}
scale(_h_multi_exclusive, scaling);
scale(_h_multi_exclusive_prompt, scaling);
scale(_h_multi_inclusive, scaling);
scale(_h_multi_inclusive_prompt, scaling);
}
}
diff --git a/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc
--- a/src/Core/AnalysisHandler.cc
+++ b/src/Core/AnalysisHandler.cc
@@ -1,519 +1,520 @@
// -*- C++ -*-
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/AnalysisHandler.hh"
#include "Rivet/Analysis.hh"
#include "Rivet/Tools/ParticleName.hh"
#include "Rivet/Tools/BeamConstraint.hh"
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Projections/Beam.hh"
#include "YODA/ReaderYODA.h"
#include "YODA/WriterYODA.h"
#include <regex>
#include <iostream>
using std::cout;
using std::cerr;
namespace {
inline std::vector<std::string> split(const std::string& input, const std::string& regex) {
// passing -1 as the submatch index parameter performs splitting
std::regex re(regex);
std::sregex_token_iterator
first{input.begin(), input.end(), re, -1},
last;
return {first, last};
}
}
namespace Rivet {
AnalysisHandler::AnalysisHandler(const string& runname)
: _runname(runname),
_eventCounter(vector<string>(), Counter()), _xs(vector<string>(), Scatter1D()),
_initialised(false), _ignoreBeams(false),
_defaultWeightIdx(0)
{}
AnalysisHandler::~AnalysisHandler()
{}
Log& AnalysisHandler::getLog() const {
return Log::getLog("Rivet.Analysis.Handler");
}
/// http://stackoverflow.com/questions/4654636/how-to-determine-if-a-string-is-a-number-with-c
bool is_number(const std::string& s)
{
std::string::const_iterator it = s.begin();
while (it != s.end() && std::isdigit(*it)) ++it;
return !s.empty() && it == s.end();
}
/// Check if any of the weightnames is not a number
bool AnalysisHandler::haveNamedWeights() {
bool dec=false;
for (unsigned int i=0;i<_weightNames.size();++i) {
string s = _weightNames[i];
if (!is_number(s)) {
dec=true;
break;
}
}
return dec;
}
void AnalysisHandler::init(const GenEvent& ge) {
if (_initialised)
throw UserError("AnalysisHandler::init has already been called: cannot re-initialize!");
// TODO TODO
// should the rivet analysis objects know about weight names?
setRunBeams(Rivet::beams(ge));
MSG_DEBUG("Initialising the analysis handler");
_eventNumber = ge.event_number();
setWeightNames(ge);
if (haveNamedWeights())
MSG_INFO("Using named weights");
else
MSG_INFO("NOT using named weights. Using first weight as nominal weight");
_eventCounter = CounterPtr(weightNames(), Counter("_EVTCOUNT"));
// set the cross section based on what is reported by this event.
// if no cross section
if (ge.cross_section()) {
MSG_TRACE("getting cross section.");
double xs = ge.cross_section()->cross_section();
double xserr = ge.cross_section()->cross_section_error();
setCrossSection(xs, xserr);
}
// Check that analyses are beam-compatible, and remove those that aren't
const size_t num_anas_requested = analysisNames().size();
vector<string> anamestodelete;
for (const AnaHandle a : _analyses) {
if (!_ignoreBeams && !a->isCompatible(beams())) {
//MSG_DEBUG(a->name() << " requires beams " << a->requiredBeams() << " @ " << a->requiredEnergies() << " GeV");
anamestodelete.push_back(a->name());
}
}
for (const string& aname : anamestodelete) {
MSG_WARNING("Analysis '" << aname << "' is incompatible with the provided beams: removing");
removeAnalysis(aname);
}
if (num_anas_requested > 0 && analysisNames().empty()) {
cerr << "All analyses were incompatible with the first event's beams\n"
<< "Exiting, since this probably wasn't intentional!" << '\n';
exit(1);
}
// Warn if any analysis' status is not unblemished
for (const AnaHandle a : analyses()) {
if (toUpper(a->status()) == "PRELIMINARY") {
MSG_WARNING("Analysis '" << a->name() << "' is preliminary: be careful, it may change and/or be renamed!");
} else if (toUpper(a->status()) == "OBSOLETE") {
MSG_WARNING("Analysis '" << a->name() << "' is obsolete: please update!");
} else if (toUpper(a->status()).find("UNVALIDATED") != string::npos) {
MSG_WARNING("Analysis '" << a->name() << "' is unvalidated: be careful, it may be broken!");
}
}
// Initialize the remaining analyses
_stage = Stage::INIT;
for (AnaHandle a : _analyses) {
MSG_DEBUG("Initialising analysis: " << a->name());
try {
// Allow projection registration in the init phase onwards
a->_allowProjReg = true;
a->init();
//MSG_DEBUG("Checking consistency of analysis: " << a->name());
//a->checkConsistency();
} catch (const Error& err) {
cerr << "Error in " << a->name() << "::init method: " << err.what() << '\n';
exit(1);
}
MSG_DEBUG("Done initialising analysis: " << a->name());
}
_stage = Stage::OTHER;
_initialised = true;
MSG_DEBUG("Analysis handler initialised");
}
void AnalysisHandler::setWeightNames(const GenEvent& ge) {
/// reroute the print output to a std::stringstream and process
/// The iteration is done over a map in hepmc2 so this is safe
std::ostringstream stream;
ge.weights().print(stream); // Super lame, I know
string str = stream.str();
std::regex re("(([^()]+))"); // Regex for stuff enclosed by parentheses ()
size_t idx = 0;
for(std::sregex_iterator i = std::sregex_iterator(str.begin(), str.end(), re);
i != std::sregex_iterator(); ++i ) {
std::smatch m = *i;
vector<string> temp = ::split(m.str(), "[,]");
if (temp.size() ==2) {
MSG_DEBUG("Name of weight #" << _weightNames.size() << ": " << temp[0]);
// store the default weight based on weight names
if (temp[0] == "Weight" || temp[0] == "0" || temp[0] == "Default") {
MSG_DEBUG(_weightNames.size() << " is being used as the nominal.");
_weightNames.push_back("");
_defaultWeightIdx = idx;
} else
_weightNames.push_back(temp[0]);
idx++;
}
}
}
void AnalysisHandler::analyze(const GenEvent& ge) {
// Call init with event as template if not already initialised
if (!_initialised) init(ge);
assert(_initialised);
// Ensure that beam details match those from the first event (if we're checking beams)
if ( !_ignoreBeams ) {
const PdgIdPair beams = Rivet::beamIds(ge);
const double sqrts = Rivet::sqrtS(ge);
if (!compatible(beams, _beams) || !fuzzyEquals(sqrts, sqrtS())) {
cerr << "Event beams mismatch: "
<< PID::toBeamsString(beams) << " @ " << sqrts/GeV << " GeV" << " vs. first beams "
<< this->beams() << " @ " << this->sqrtS()/GeV << " GeV" << '\n';
exit(1);
}
}
// Create the Rivet event wrapper
/// @todo Filter/normalize the event here
Event event(ge);
// set the cross section based on what is reported by this event.
// if no cross section
MSG_TRACE("getting cross section.");
if (ge.cross_section()) {
MSG_TRACE("getting cross section from GenEvent.");
double xs = ge.cross_section()->cross_section();
double xserr = ge.cross_section()->cross_section_error();
setCrossSection(xs, xserr);
}
// won't happen for first event because _eventNumber is set in
// init()
if (_eventNumber != ge.event_number()) {
/// @todo
/// can we get away with not passing a matrix?
MSG_TRACE("AnalysisHandler::analyze(): Pushing _eventCounter to persistent.");
_eventCounter.get()->pushToPersistent(_subEventWeights);
// if this is indeed a new event, push the temporary
// histograms and reset
for (const AnaHandle& a : _analyses) {
for (auto ao : a->analysisObjects()) {
MSG_TRACE("AnalysisHandler::analyze(): Pushing " << a->name() << "'s " << ao->name() << " to persistent.");
ao.get()->pushToPersistent(_subEventWeights);
}
MSG_TRACE("AnalysisHandler::analyze(): finished pushing " << a->name() << "'s objects to persistent.");
}
_eventNumber = ge.event_number();
MSG_DEBUG("nominal event # " << _eventCounter.get()->_persistent[0]->numEntries());
MSG_DEBUG("nominal sum of weights: " << _eventCounter.get()->_persistent[0]->sumW());
MSG_DEBUG("Event has " << _subEventWeights.size() << " sub events.");
_subEventWeights.clear();
}
MSG_TRACE("starting new sub event");
_eventCounter.get()->newSubEvent();
for (const AnaHandle& a : _analyses) {
for (auto ao : a->analysisObjects()) {
ao.get()->newSubEvent();
}
}
_subEventWeights.push_back(event.weights());
MSG_DEBUG("Analyzing subevent #" << _subEventWeights.size() - 1 << ".");
_eventCounter->fill();
// Run the analyses
for (AnaHandle a : _analyses) {
MSG_TRACE("About to run analysis " << a->name());
try {
a->analyze(event);
} catch (const Error& err) {
cerr << "Error in " << a->name() << "::analyze method: " << err.what() << '\n';
exit(1);
}
MSG_TRACE("Finished running analysis " << a->name());
}
}
void AnalysisHandler::analyze(const GenEvent* ge) {
if (ge == NULL) {
MSG_ERROR("AnalysisHandler received null pointer to GenEvent");
//throw Error("AnalysisHandler received null pointer to GenEvent");
}
analyze(*ge);
}
void AnalysisHandler::finalize() {
if (!_initialised) return;
MSG_INFO("Finalising analyses");
MSG_TRACE("AnalysisHandler::finalize(): Pushing analysis objects to persistent.");
_eventCounter.get()->pushToPersistent(_subEventWeights);
for (const AnaHandle& a : _analyses) {
for (auto ao : a->analysisObjects())
ao.get()->pushToPersistent(_subEventWeights);
}
for (const AnaHandle& a : _analyses) {
for (size_t iW = 0; iW < numWeights(); iW++) {
_eventCounter.get()->setActiveWeightIdx(iW);
_xs.get()->setActiveWeightIdx(iW);
for (auto ao : a->analysisObjects())
ao.get()->setActiveWeightIdx(iW);
MSG_TRACE("running " << a->name() << "::finalize() for weight " << iW << ".");
try {
a->finalize();
} catch (const Error& err) {
cerr << "Error in " << a->name() << "::finalize method: " << err.what() << '\n';
exit(1);
}
}
}
// Print out number of events processed
MSG_INFO("Processed " << numEvents() << " event" << (numEvents() == 1 ? "" : "s"));
// Print out MCnet boilerplate
cout << '\n';
cout << "The MCnet usage guidelines apply to Rivet: see http://www.montecarlonet.org/GUIDELINES" << '\n';
cout << "Please acknowledge plots made with Rivet analyses, and cite arXiv:1003.0694 (http://arxiv.org/abs/1003.0694)" << '\n';
}
AnalysisHandler& AnalysisHandler::addAnalysis(const string& analysisname) {
// Check for a duplicate analysis
/// @todo Might we want to be able to run an analysis twice, with different params?
/// Requires avoiding histo tree clashes, i.e. storing the histos on the analysis objects.
for (const AnaHandle& a : _analyses) {
if (a->name() == analysisname) {
MSG_WARNING("Analysis '" << analysisname << "' already registered: skipping duplicate");
return *this;
}
}
AnaHandle analysis( AnalysisLoader::getAnalysis(analysisname) );
if (analysis.get() != 0) { // < Check for null analysis.
MSG_DEBUG("Adding analysis '" << analysisname << "'");
analysis->_analysishandler = this;
_analyses.insert(analysis);
} else {
MSG_WARNING("Analysis '" << analysisname << "' not found.");
}
// MSG_WARNING(_analyses.size());
// for (const AnaHandle& a : _analyses) MSG_WARNING(a->name());
return *this;
}
AnalysisHandler& AnalysisHandler::removeAnalysis(const string& analysisname) {
std::shared_ptr<Analysis> toremove;
for (const AnaHandle a : _analyses) {
if (a->name() == analysisname) {
toremove = a;
break;
}
}
if (toremove.get() != 0) {
MSG_DEBUG("Removing analysis '" << analysisname << "'");
_analyses.erase(toremove);
}
return *this;
}
void AnalysisHandler::addData(const std::vector<YODA::AnalysisObjectPtr>& aos) {
for (const YODA::AnalysisObjectPtr ao : aos) {
const string path = ao->path();
if (path.size() > 1) { // path > "/"
try {
const string ananame = ::split(path, "/")[0];
AnaHandle a = analysis(ananame);
//MultiweightAOPtr mao = ????; /// @todo generate right Multiweight object from ao
//a->addAnalysisObject(mao); /// @todo Need to statistically merge...
} catch (const Error& e) {
MSG_WARNING(e.what());
}
}
}
}
void AnalysisHandler::readData(const string& filename) {
vector<YODA::AnalysisObjectPtr> aos;
try {
/// @todo Use new YODA SFINAE to fill the smart ptr vector directly
vector<YODA::AnalysisObject*> aos_raw;
YODA::ReaderYODA::read(filename, aos_raw);
for (YODA::AnalysisObject* aor : aos_raw) aos.push_back(YODA::AnalysisObjectPtr(aor));
} catch (const YODA::ReadError & e) {
throw UserError("Unexpected error in reading file: " + filename);
}
if (!aos.empty()) addData(aos);
}
vector<MultiweightAOPtr> AnalysisHandler::getRivetAOs() const {
vector<MultiweightAOPtr> rtn;
for (AnaHandle a : _analyses) {
for (const auto & ao : a->analysisObjects()) {
rtn.push_back(ao);
}
}
rtn.push_back(_eventCounter);
rtn.push_back(_xs);
return rtn;
}
vector<YODA::AnalysisObjectPtr> AnalysisHandler::getYodaAOs() const {
vector<YODA::AnalysisObjectPtr> rtn;
for (auto rao : getRivetAOs()) {
// need to set the index
// before we can search the PATH
rao.get()->setActiveWeightIdx(_defaultWeightIdx);
if (rao->path().find("/TMP/") != string::npos)
continue;
for (size_t iW = 0; iW < numWeights(); iW++) {
rao.get()->setActiveWeightIdx(iW);
rtn.push_back(rao.get()->activeYODAPtr());
}
}
sort(rtn.begin(), rtn.end(),
[](YODA::AnalysisObjectPtr a, YODA::AnalysisObjectPtr b) {
return a->path() < b->path();
}
);
return rtn;
}
vector<YODA::AnalysisObjectPtr> AnalysisHandler::getData() const {
return getYodaAOs();
}
void AnalysisHandler::writeData(const string& filename) const {
const vector<YODA::AnalysisObjectPtr> aos = getData();
try {
YODA::WriterYODA::write(filename, aos.begin(), aos.end());
} catch ( YODA::WriteError ) {
throw UserError("Unexpected error in writing file: " + filename);
}
}
string AnalysisHandler::runName() const { return _runname; }
size_t AnalysisHandler::numEvents() const { return _eventCounter->numEntries(); }
std::vector<std::string> AnalysisHandler::analysisNames() const {
std::vector<std::string> rtn;
for (AnaHandle a : _analyses) {
rtn.push_back(a->name());
}
return rtn;
}
const AnaHandle AnalysisHandler::analysis(const std::string& analysisname) const {
for (const AnaHandle a : analyses())
if (a->name() == analysisname) return a;
throw Error("No analysis named '" + analysisname + "' registered in AnalysisHandler");
}
AnalysisHandler& AnalysisHandler::addAnalyses(const std::vector<std::string>& analysisnames) {
for (const string& aname : analysisnames) {
//MSG_DEBUG("Adding analysis '" << aname << "'");
addAnalysis(aname);
}
return *this;
}
AnalysisHandler& AnalysisHandler::removeAnalyses(const std::vector<std::string>& analysisnames) {
for (const string& aname : analysisnames) {
removeAnalysis(aname);
}
return *this;
}
AnalysisHandler& AnalysisHandler::setCrossSection(double xs, double xserr) {
_xs = Scatter1DPtr(weightNames(), Scatter1D("_XSEC"));
_eventCounter.get()->setActiveWeightIdx(_defaultWeightIdx);
double nomwgt = sumOfWeights();
- // the cross section of each weight variation is the nominal cross section
+ // The cross section of each weight variation is the nominal cross section
// times the sumOfWeights(variation) / sumOfWeights(nominal).
- // this way the cross section will work correctly
- for (int iW = 0; iW < numWeights(); iW++) {
- _eventCounter.get()->setActiveWeightIdx(iW);
- double s = sumOfWeights() / nomwgt;
- _xs.get()->setActiveWeightIdx(iW);
- _xs->addPoint(xs*s, xserr*s);
+ // This way the cross section will work correctly
+ for (size_t iW = 0; iW < numWeights(); iW++) {
+ _eventCounter.get()->setActiveWeightIdx(iW);
+ double s = sumOfWeights() / nomwgt;
+ _xs.get()->setActiveWeightIdx(iW);
+ _xs->addPoint(xs*s, xserr*s);
}
_eventCounter.get()->unsetActiveWeight();
_xs.get()->unsetActiveWeight();
return *this;
}
+
AnalysisHandler& AnalysisHandler::addAnalysis(Analysis* analysis) {
analysis->_analysishandler = this;
_analyses.insert(AnaHandle(analysis));
return *this;
}
PdgIdPair AnalysisHandler::beamIds() const {
return Rivet::beamIds(beams());
}
double AnalysisHandler::sqrtS() const {
return Rivet::sqrtS(beams());
}
void AnalysisHandler::setIgnoreBeams(bool ignore) {
_ignoreBeams=ignore;
}
}
diff --git a/src/Projections/DressedLeptons.cc b/src/Projections/DressedLeptons.cc
--- a/src/Projections/DressedLeptons.cc
+++ b/src/Projections/DressedLeptons.cc
@@ -1,137 +1,137 @@
// -*- C++ -*-
#include "Rivet/Projections/DressedLeptons.hh"
namespace Rivet {
// On DressedLepton helper class
//{
DressedLepton::DressedLepton(const Particle& dlepton)
: Particle(dlepton)
- {
+ {
setConstituents({{dlepton}}); //< bare lepton is first constituent
}
DressedLepton::DressedLepton(const Particle& lepton, const Particles& photons, bool momsum)
: Particle(lepton.pid(), lepton.momentum())
{
setConstituents({{lepton}}); //< bare lepton is first constituent
addConstituents(photons, momsum);
}
void DressedLepton::addPhoton(const Particle& p, bool momsum) {
if (p.pid() != PID::PHOTON) throw Error("Clustering a non-photon on to a DressedLepton:"+to_string(p.pid()));
addConstituent(p, momsum);
}
const Particle& DressedLepton::bareLepton() const {
const Particle& l = constituents().front();
if (!l.isChargedLepton()) throw Error("First constituent of a DressedLepton is not a bare lepton: oops");
return l;
}
//}
// Separate-FS version
DressedLeptons::DressedLeptons(const FinalState& photons, const FinalState& bareleptons,
double dRmax, const Cut& cut, bool useDecayPhotons)
: FinalState(cut),
_dRmax(dRmax), _fromDecay(useDecayPhotons)
{
setName("DressedLeptons");
IdentifiedFinalState photonfs(photons, PID::PHOTON);
declare(photonfs, "Photons");
IdentifiedFinalState leptonfs(bareleptons);
leptonfs.acceptIdPairs({PID::ELECTRON, PID::MUON, PID::TAU});
declare(leptonfs, "Leptons");
}
// Single-FS version
DressedLeptons::DressedLeptons(const FinalState& barefs,
double dRmax, const Cut& cut, bool useDecayPhotons)
: DressedLeptons(barefs, barefs, dRmax, cut, useDecayPhotons)
{ }
CmpState DressedLeptons::compare(const Projection& p) const {
// Compare the two as final states (for pT and eta cuts)
const DressedLeptons& other = dynamic_cast<const DressedLeptons&>(p);
CmpState fscmp = FinalState::compare(other);
if (fscmp != CmpState::EQ) return fscmp;
const PCmp phcmp = mkNamedPCmp(p, "Photons");
if (phcmp != CmpState::EQ) return phcmp;
const PCmp sigcmp = mkNamedPCmp(p, "Leptons");
if (sigcmp != CmpState::EQ) return sigcmp;
return (cmp(_dRmax, other._dRmax) ||
cmp(_fromDecay, other._fromDecay));
}
void DressedLeptons::project(const Event& e) {
_theParticles.clear();
// Get bare leptons
const FinalState& signal = applyProjection<FinalState>(e, "Leptons");
Particles bareleptons = signal.particles();
if (bareleptons.empty()) return;
// Initialise DL collection with bare leptons
vector<Particle> allClusteredLeptons;
allClusteredLeptons.reserve(bareleptons.size());
for (const Particle& bl : bareleptons) {
Particle dl(bl.pid(), bl.momentum());
dl.setConstituents({bl});
allClusteredLeptons += dl;
}
// If the radius is 0 or negative, don't even attempt to cluster
if (_dRmax > 0) {
// Match each photon to its closest charged lepton within the dR cone
const FinalState& photons = applyProjection<FinalState>(e, "Photons");
for (const Particle& photon : photons.particles()) {
// Ignore photon if it's from a hadron/tau decay and we're avoiding those
- if (!_fromDecay && photon.fromDecay()) continue;
+ if (!_fromDecay && !photon.isPrompt()) continue;
const FourMomentum& p_P = photon.momentum();
double dRmin = _dRmax;
int idx = -1;
for (size_t i = 0; i < bareleptons.size(); ++i) {
const Particle& bl = bareleptons[i];
// Only cluster photons around *charged* signal particles
if (bl.charge3() == 0) continue;
// Find the closest lepton
double dR = deltaR(bl, p_P);
if (dR < dRmin) {
dRmin = dR;
idx = i;
}
}
if (idx > -1) allClusteredLeptons[idx].addConstituent(photon, true);
}
}
// Fill the canonical particles collection with the composite DL Particles
for (const Particle& lepton : allClusteredLeptons) {
const bool acc = accept(lepton);
- MSG_TRACE("Clustered lepton " << lepton
- << " with constituents = " << lepton.constituents()
+ MSG_TRACE("Clustered lepton " << lepton
+ << " with constituents = " << lepton.constituents()
<< ", cut-pass = " << std::boolalpha << acc);
if (acc) _theParticles.push_back(lepton);
}
- MSG_DEBUG("#dressed leptons = " << allClusteredLeptons.size()
+ MSG_DEBUG("#dressed leptons = " << allClusteredLeptons.size()
<< " -> " << _theParticles.size() << " after cuts");
}
}
diff --git a/src/Projections/FastJets.cc b/src/Projections/FastJets.cc
--- a/src/Projections/FastJets.cc
+++ b/src/Projections/FastJets.cc
@@ -1,215 +1,215 @@
// -*- C++ -*-
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/HeavyHadrons.hh"
#include "Rivet/Projections/TauFinder.hh"
namespace Rivet {
void FastJets::_initBase() {
setName("FastJets");
declare(HeavyHadrons(), "HFHadrons");
declare(TauFinder(TauFinder::DecayMode::HADRONIC), "Taus");
}
void FastJets::_initJdef(Algo alg, double rparameter, double seed_threshold) {
MSG_DEBUG("JetAlg = " << static_cast<int>(alg));
MSG_DEBUG("R parameter = " << rparameter);
MSG_DEBUG("Seed threshold = " << seed_threshold);
if (alg == KT) {
_jdef = fastjet::JetDefinition(fastjet::kt_algorithm, rparameter, fastjet::E_scheme);
} else if (alg == CAM) {
_jdef = fastjet::JetDefinition(fastjet::cambridge_algorithm, rparameter, fastjet::E_scheme);
} else if (alg == ANTIKT) {
_jdef = fastjet::JetDefinition(fastjet::antikt_algorithm, rparameter, fastjet::E_scheme);
} else if (alg == DURHAM) {
_jdef = fastjet::JetDefinition(fastjet::ee_kt_algorithm, fastjet::E_scheme);
} else if (alg == GENKTEE) {
_jdef = fastjet::JetDefinition(fastjet::ee_genkt_algorithm, rparameter, -1);
} else {
// Plugins:
if (alg == SISCONE) {
const double OVERLAP_THRESHOLD = 0.75;
_plugin.reset(new fastjet::SISConePlugin(rparameter, OVERLAP_THRESHOLD));
// } else if (alg == PXCONE) {
// string msg = "PxCone currently not supported, since FastJet doesn't install it by default. ";
// msg += "Please notify the Rivet authors if this behaviour should be changed.";
// throw Error(msg);
// _plugin.reset(new fastjet::PxConePlugin(rparameter));
} else if (alg == ATLASCONE) {
const double OVERLAP_THRESHOLD = 0.5;
_plugin.reset(new fastjet::ATLASConePlugin(rparameter, seed_threshold, OVERLAP_THRESHOLD));
} else if (alg == CMSCONE) {
_plugin.reset(new fastjet::CMSIterativeConePlugin(rparameter, seed_threshold));
} else if (alg == CDFJETCLU) {
const double OVERLAP_THRESHOLD = 0.75;
_plugin.reset(new fastjet::CDFJetCluPlugin(rparameter, OVERLAP_THRESHOLD, seed_threshold));
} else if (alg == CDFMIDPOINT) {
const double OVERLAP_THRESHOLD = 0.5;
_plugin.reset(new fastjet::CDFMidPointPlugin(rparameter, OVERLAP_THRESHOLD, seed_threshold));
} else if (alg == D0ILCONE) {
const double min_jet_Et = 6.0;
_plugin.reset(new fastjet::D0RunIIConePlugin(rparameter, min_jet_Et));
} else if (alg == JADE) {
_plugin.reset(new fastjet::JadePlugin());
} else if (alg == TRACKJET) {
_plugin.reset(new fastjet::TrackJetPlugin(rparameter));
}
_jdef = fastjet::JetDefinition(_plugin.get());
}
}
CmpState FastJets::compare(const Projection& p) const {
const FastJets& other = dynamic_cast<const FastJets&>(p);
return \
cmp(_useMuons, other._useMuons) ||
cmp(_useInvisibles, other._useInvisibles) ||
mkNamedPCmp(other, "FS") ||
cmp(_jdef.jet_algorithm(), other._jdef.jet_algorithm()) ||
cmp(_jdef.recombination_scheme(), other._jdef.recombination_scheme()) ||
cmp(_jdef.plugin(), other._jdef.plugin()) ||
cmp(_jdef.R(), other._jdef.R()) ||
cmp(_adef, other._adef);
}
// STATIC
PseudoJets FastJets::mkClusterInputs(const Particles& fsparticles, const Particles& tagparticles) {
PseudoJets pjs;
/// @todo Use FastJet3's UserInfo system to store Particle pointers directly?
// Store 4 vector data about each particle into FastJet's PseudoJets
for (size_t i = 0; i < fsparticles.size(); ++i) {
fastjet::PseudoJet pj = fsparticles[i];
pj.set_user_index(i+1);
pjs.push_back(pj);
}
// And the same for ghost tagging particles (with negative user indices)
for (size_t i = 0; i < tagparticles.size(); ++i) {
fastjet::PseudoJet pj = tagparticles[i];
pj *= 1e-20; ///< Ghostify the momentum
pj.set_user_index(-i-1);
pjs.push_back(pj);
}
return pjs;
}
// STATIC
Jet FastJets::mkJet(const PseudoJet& pj, const Particles& fsparticles, const Particles& tagparticles) {
const PseudoJets pjconstituents = pj.constituents();
Particles constituents, tags;
constituents.reserve(pjconstituents.size());
for (const fastjet::PseudoJet& pjc : pjconstituents) {
// Pure ghosts don't have corresponding particles
if (pjc.has_area() && pjc.is_pure_ghost()) continue;
// Default user index = 0 doesn't give valid particle lookup
if (pjc.user_index() == 0) continue;
// Split by index sign into constituent & tag lookup
if (pjc.user_index() > 0) {
// Find constituents if index > 0
const size_t i = pjc.user_index() - 1;
if (i >= fsparticles.size()) throw RangeError("FS particle lookup failed in jet construction");
constituents.push_back(fsparticles.at(i));
} else if (!tagparticles.empty()) {
// Find tags if index < 0
const size_t i = abs(pjc.user_index()) - 1;
if (i >= tagparticles.size()) throw RangeError("Tag particle lookup failed in jet construction");
tags.push_back(tagparticles.at(i));
}
}
return Jet(pj, constituents, tags);
}
// STATIC
Jets FastJets::mkJets(const PseudoJets& pjs, const Particles& fsparticles, const Particles& tagparticles) {
Jets rtn; rtn.reserve(pjs.size());
for (const PseudoJet pj : pjs) {
rtn.push_back(FastJets::mkJet(pj, fsparticles, tagparticles));
}
return rtn;
}
void FastJets::project(const Event& e) {
// Assemble final state particles
const string fskey = (_useInvisibles == JetAlg::Invisibles::NONE) ? "VFS" : "FS";
Particles fsparticles = applyProjection<FinalState>(e, fskey).particles();
// Remove prompt invisibles if needed (already done by VFS if using NO_INVISIBLES)
if (_useInvisibles == JetAlg::Invisibles::DECAY) {
- ifilter_discard(fsparticles, [](const Particle& p) { return !(p.isVisible() || p.fromDecay()); });
+ ifilter_discard(fsparticles, [](const Particle& p) { return !p.isVisible() && p.isPrompt(); });
}
// Remove prompt/all muons if needed
if (_useMuons == JetAlg::Muons::DECAY) {
- ifilter_discard(fsparticles, [](const Particle& p) { return isMuon(p) && !p.fromDecay(); });
+ ifilter_discard(fsparticles, [](const Particle& p) { return isMuon(p) && p.isPrompt(); });
} else if (_useMuons == JetAlg::Muons::NONE) {
ifilter_discard(fsparticles, isMuon);
}
// Tagging particles
const Particles chadrons = applyProjection<HeavyHadrons>(e, "HFHadrons").cHadrons();
const Particles bhadrons = applyProjection<HeavyHadrons>(e, "HFHadrons").bHadrons();
const Particles taus = applyProjection<FinalState>(e, "Taus").particles();
calc(fsparticles, chadrons+bhadrons+taus);
}
void FastJets::calc(const Particles& fsparticles, const Particles& tagparticles) {
MSG_DEBUG("Finding jets from " << fsparticles.size() << " input particles + " << tagparticles.size() << " tagging particles");
_fsparticles = fsparticles;
_tagparticles = tagparticles;
// Make pseudojets, with mapping info to Rivet FS and tag particles
PseudoJets pjs = mkClusterInputs(_fsparticles, _tagparticles);
// Run either basic or area-calculating cluster sequence as reqd.
if (_adef) {
_cseq.reset(new fastjet::ClusterSequenceArea(pjs, _jdef, *_adef));
} else {
_cseq.reset(new fastjet::ClusterSequence(pjs, _jdef));
}
MSG_DEBUG("ClusterSequence constructed; Njets_tot = "
<< _cseq->inclusive_jets().size() << ", Njets(pT > 10 GeV) = "
<< _cseq->inclusive_jets(10*GeV).size());
}
void FastJets::reset() {
_yscales.clear();
_fsparticles.clear();
_tagparticles.clear();
/// @todo _cseq = fastjet::ClusterSequence();
}
Jets FastJets::_jets() const {
/// @todo Cache?
return mkJets(pseudojets(), _fsparticles, _tagparticles);
}
/// @todo "Automate" trimming as part of project() with pre-registered Filters
Jet FastJets::trimJet(const Jet& input, const fastjet::Filter& trimmer) const {
if (input.pseudojet().associated_cluster_sequence() != clusterSeq().get())
throw Error("To trim a Rivet::Jet, its associated PseudoJet must have come from this FastJets' ClusterSequence");
PseudoJet pj = trimmer(input);
return mkJet(pj, _fsparticles, _tagparticles);
}
PseudoJets FastJets::pseudoJets(double ptmin) const {
return clusterSeq() ? clusterSeq()->inclusive_jets(ptmin) : PseudoJets();
}
}
diff --git a/src/Projections/FinalPartons.cc b/src/Projections/FinalPartons.cc
--- a/src/Projections/FinalPartons.cc
+++ b/src/Projections/FinalPartons.cc
@@ -1,42 +1,42 @@
// -*- C++ -*-
#include "Rivet/Projections/FinalPartons.hh"
#include "Rivet/Tools/RivetHepMC.hh"
namespace Rivet {
bool FinalPartons::accept(const Particle& p) const {
// Reject if *not* a parton
if (!isParton(p))
return false;
// Accept partons if they end on a standard hadronization vertex
if (p.genParticle()->end_vertex() != nullptr && p.genParticle()->end_vertex()->id() == 5)
return true;
// Reject if p has a parton child.
for (const Particle& c : p.children())
if (isParton(c))
return false;
// Reject if from a hadron or tau decay.
- if (p.fromDecay())
+ if (p.fromHadron() || p.fromTau())
return false;
return _cuts->accept(p);
}
void FinalPartons::project(const Event& e) {
_theParticles.clear();
for (const GenParticle* gp : Rivet::particles(e.genEvent())) {
if (!gp) continue;
const Particle p(gp);
if (accept(p)) _theParticles.push_back(p);
}
}
}

File Metadata

Mime Type
text/x-diff
Expires
Wed, May 14, 11:50 AM (4 h, 42 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5111534
Default Alt Text
(186 KB)

Event Timeline