Page MenuHomeHEPForge
Files F7879915

No OneTemporary
Actions

View Options

diff --git a/analyses/pluginMC/MC_WJETS.cc b/analyses/pluginMC/MC_WJETS.cc
--- a/analyses/pluginMC/MC_WJETS.cc
+++ b/analyses/pluginMC/MC_WJETS.cc
@@ -1,130 +1,130 @@
// -*- C++ -*-
#include "Rivet/Analyses/MC_JetAnalysis.hh"
#include "Rivet/Projections/WFinder.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Analysis.hh"
namespace Rivet {
/// @brief MC validation analysis for W + jets events
class MC_WJETS : public MC_JetAnalysis {
public:
/// Default constructor
MC_WJETS(string name="MC_WJETS")
: MC_JetAnalysis(name, 4, "Jets")
{
_dR=0.2;
_lepton=PID::ELECTRON;
}
/// @name Analysis methods
//@{
/// Book histograms
void init() {
FinalState fs;
WFinder wfinder(fs, Cuts::abseta < 3.5 && Cuts::pT > 25*GeV, _lepton, 60.0*GeV, 100.0*GeV, 25.0*GeV, _dR);
declare(wfinder, "WFinder");
FastJets jetpro(wfinder.remainingFinalState(), FastJets::ANTIKT, 0.4);
declare(jetpro, "Jets");
book(_h_W_jet1_deta ,"W_jet1_deta", 50, -5.0, 5.0);
book(_h_W_jet1_dR ,"W_jet1_dR", 25, 0.5, 7.0);
MC_JetAnalysis::init();
}
/// Do the analysis
void analyze(const Event & e) {
const double weight = 1.0;
const WFinder& wfinder = apply<WFinder>(e, "WFinder");
if (wfinder.bosons().size() != 1) {
vetoEvent;
}
FourMomentum wmom(wfinder.bosons().front().momentum());
const Jets& jets = apply<FastJets>(e, "Jets").jetsByPt(_jetptcut);
if (jets.size() > 0) {
- _h_W_jet1_deta->fill(wmom.eta()-jets[0].eta(), weight);
- _h_W_jet1_dR->fill(deltaR(wmom, jets[0].momentum()), weight);
+ _h_W_jet1_deta->fill(wmom.eta()-jets[0].eta());
+ _h_W_jet1_dR->fill(deltaR(wmom, jets[0].momentum()));
}
MC_JetAnalysis::analyze(e);
}
/// Finalize
void finalize() {
scale(_h_W_jet1_deta, crossSection()/picobarn/sumOfWeights());
scale(_h_W_jet1_dR, crossSection()/picobarn/sumOfWeights());
MC_JetAnalysis::finalize();
}
//@}
protected:
/// @name Parameters for specialised e/mu and dressed/bare subclassing
//@{
double _dR;
PdgId _lepton;
//@}
private:
/// @name Histograms
//@{
Histo1DPtr _h_W_jet1_deta;
Histo1DPtr _h_W_jet1_dR;
//@}
};
struct MC_WJETS_EL : public MC_WJETS {
MC_WJETS_EL() : MC_WJETS("MC_WJETS_EL") {
_dR = 0.2;
_lepton = PID::ELECTRON;
}
};
struct MC_WJETS_EL_BARE : public MC_WJETS {
MC_WJETS_EL_BARE() : MC_WJETS("MC_WJETS_EL_BARE") {
_dR = 0;
_lepton = PID::ELECTRON;
}
};
struct MC_WJETS_MU : public MC_WJETS {
MC_WJETS_MU() : MC_WJETS("MC_WJETS_MU") {
_dR = 0.2;
_lepton = PID::MUON;
}
};
struct MC_WJETS_MU_BARE : public MC_WJETS {
MC_WJETS_MU_BARE() : MC_WJETS("MC_WJETS_MU_BARE") {
_dR = 0;
_lepton = PID::MUON;
}
};
// The hooks for the plugin system
DECLARE_RIVET_PLUGIN(MC_WJETS);
DECLARE_RIVET_PLUGIN(MC_WJETS_EL);
DECLARE_RIVET_PLUGIN(MC_WJETS_EL_BARE);
DECLARE_RIVET_PLUGIN(MC_WJETS_MU);
DECLARE_RIVET_PLUGIN(MC_WJETS_MU_BARE);
}
diff --git a/include/Rivet/Tools/Cmp.fhh b/include/Rivet/Tools/Cmp.fhh
--- a/include/Rivet/Tools/Cmp.fhh
+++ b/include/Rivet/Tools/Cmp.fhh
@@ -1,18 +1,27 @@
// -*- C++ -*-
#ifndef RIVET_Cmp_FHH
#define RIVET_Cmp_FHH
namespace Rivet {
// Forward-declare the Cmp template class
template <typename T>
class Cmp;
enum class CmpState {
UNDEF, LT, EQ, GT, NEQ
};
+ /*std::ostream& operator << (std::ostream& os, const CmpState& obj) {
+ if (obj == CmpState::UNDEF) os << "UNDEF";
+ if (obj == CmpState::LT) os << "LT";
+ if (obj == CmpState::GT) os << "GT";
+ if (obj == CmpState::NEQ) os << "NEQ";
+ return os;
+ }*/
+
+
}
#endif
diff --git a/src/Projections/WFinder.cc b/src/Projections/WFinder.cc
--- a/src/Projections/WFinder.cc
+++ b/src/Projections/WFinder.cc
@@ -1,163 +1,166 @@
// -*- C++ -*-
#include "Rivet/Projections/WFinder.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
#include "Rivet/Projections/PromptFinalState.hh"
#include "Rivet/Projections/InvMassFinalState.hh"
#include "Rivet/Projections/MergedFinalState.hh"
#include "Rivet/Projections/DressedLeptons.hh"
#include "Rivet/Projections/VetoedFinalState.hh"
#include "Rivet/Projections/Beam.hh"
namespace Rivet {
WFinder::WFinder(const FinalState& inputfs,
const Cut& leptoncuts,
PdgId pid,
double minmass, double maxmass,
double missingET,
double dRmax,
ChargedLeptons chLeptons,
ClusterPhotons clusterPhotons,
AddPhotons trackPhotons,
MassWindow masstype,
double masstarget) {
setName("WFinder");
_etMissMin = missingET;
_minmass = minmass;
_maxmass = maxmass;
_masstarget = masstarget;
_pid = abs(pid);
_trackPhotons = trackPhotons;
_useTransverseMass = (masstype == MassWindow::MT);
// Check that the arguments are legal
if (_pid != PID::ELECTRON && _pid != PID::MUON)
throw Error("Invalid charged lepton PID given to WFinder");
// Identify bare leptons for dressing
// Bit of a code nightmare -- FS projection copy constructors don't work?
/// @todo Fix FS copy constructors!!
if (chLeptons == ChargedLeptons::PROMPT) {
PromptFinalState inputfs_prompt(inputfs);
IdentifiedFinalState bareleptons(inputfs_prompt);
bareleptons.acceptIdPair(_pid);
declare(bareleptons, "BareLeptons");
} else {
IdentifiedFinalState bareleptons(inputfs);
bareleptons.acceptIdPair(_pid);
declare(bareleptons, "BareLeptons");
}
// Dress the bare leptons
const bool doClustering = (clusterPhotons != ClusterPhotons::NONE);
const bool useDecayPhotons = (clusterPhotons == ClusterPhotons::ALL);
DressedLeptons leptons(inputfs, get<FinalState>("BareLeptons"), (doClustering ? dRmax : -1.), leptoncuts, useDecayPhotons);
declare(leptons, "DressedLeptons");
// Add MissingMomentum proj to calc MET
MissingMomentum vismom(inputfs);
declare(vismom, "MissingET");
// Identify the non-Z part of the event
VetoedFinalState remainingFS;
remainingFS.addVetoOnThisFinalState(*this);
declare(remainingFS, "RFS");
}
/////////////////////////////////////////////////////
const VetoedFinalState& WFinder::remainingFinalState() const {
return getProjection<VetoedFinalState>("RFS");
}
const MissingMomentum& WFinder::missingMom() const {
return getProjection<MissingMomentum>("MissingET");
}
CmpState WFinder::compare(const Projection& p) const {
PCmp dlcmp = mkNamedPCmp(p, "DressedLeptons");
if (dlcmp != CmpState::EQ) return dlcmp;
+ PCmp metcmp = mkNamedPCmp(p, "MissingET");
+ if (metcmp != CmpState::EQ) return metcmp;
+
const WFinder& other = dynamic_cast<const WFinder&>(p);
return (cmp(_minmass, other._minmass) ||
cmp(_maxmass, other._maxmass) ||
cmp(_useTransverseMass, other._useTransverseMass) ||
cmp(_etMissMin, other._etMissMin) ||
cmp(_pid, other._pid) ||
cmp(_trackPhotons, other._trackPhotons));
}
void WFinder::project(const Event& e) {
clear();
_leptons.clear();
_neutrinos.clear();
// Check missing ET
const MissingMomentum& missmom = applyProjection<MissingMomentum>(e, "MissingET");
const double met = missmom.vectorEt().mod();
MSG_TRACE("MET = " << met/GeV << " GeV vs. required > " << _etMissMin/GeV << " GeV");
if (met < _etMissMin) {
MSG_DEBUG("Not enough missing ET: " << met/GeV << " GeV vs. required > " << _etMissMin/GeV << " GeV");
return;
}
// Get lepton
const DressedLeptons& leptons = applyProjection<DressedLeptons>(e, "DressedLeptons");
if ( leptons.dressedLeptons().empty() ) {
MSG_DEBUG("No dressed leptons");
return;
}
MSG_DEBUG("Found at least one dressed lepton: " << leptons.dressedLeptons().front().momentum() );
// Get missing momentum 4-vector, assuming a massless invisible particle
const FourMomentum pmiss = missmom.missingMomentum(0*GeV);
MSG_DEBUG("Found missing 4-momentum: " << pmiss);
// Compute an invariant mass final state for the W decay leptons (using pseudo-neutrinos from ETmiss)
PdgId _nu_pid = _pid + 1;
assert(_nu_pid == PID::NU_E || _nu_pid == PID::NU_MU);
vector<pair<PdgId, PdgId> > l_nu_ids;
l_nu_ids += make_pair(_pid, -_nu_pid);
l_nu_ids += make_pair(-_pid, _nu_pid);
InvMassFinalState imfs(l_nu_ids, _minmass, _maxmass, _masstarget);
imfs.useTransverseMass(_useTransverseMass);
Particles tmp = leptons.particles();
tmp += { Particle( _nu_pid, pmiss), Particle(-_nu_pid, pmiss) }; // fake (anti)neutrinos from ETmiss vector
imfs.calc(tmp);
if (imfs.particlePairs().size() < 1) return;
// Assemble a pseudo-W particle
const ParticlePair Wconstituents = imfs.particlePairs().front();
const Particle& p1(Wconstituents.first), p2(Wconstituents.second);
const FourMomentum pW = p1.momentum() + p2.momentum();
const int wcharge3 = p1.charge3() + p2.charge3();
assert(abs(wcharge3) == 3);
const int wcharge = wcharge3/3;
const PdgId wpid = (wcharge == 1) ? PID::WPLUSBOSON : PID::WMINUSBOSON;
Particle w(wpid, pW);
MSG_DEBUG(w << " reconstructed from: " << p1 << " + " << p2);
// Add (dressed) lepton constituents to the W (skipping photons if requested)
/// @todo Do we need to add all used invisibles to _theParticles ?
const Particle l = p1.isChargedLepton() ? p1 : p2;
_leptons += (_trackPhotons == AddPhotons::YES) ? l : l.constituents().front();
w.addConstituent(_leptons.back());
const Particle nu = p1.isNeutrino() ? p1 : p2;
_neutrinos += nu;
w.addConstituent(nu);
// Register the completed W
_theParticles.push_back(w);
}
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 9:22 PM (1 d, 1 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3806262
Default Alt Text
(10 KB)

Event Timeline