diff --git a/analyses/pluginMC/EXAMPLE_SMEAR.cc b/analyses/pluginMC/EXAMPLE_SMEAR.cc --- a/analyses/pluginMC/EXAMPLE_SMEAR.cc +++ b/analyses/pluginMC/EXAMPLE_SMEAR.cc @@ -1,238 +1,238 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/TauFinder.hh" #include "Rivet/Projections/SmearedJets.hh" #include "Rivet/Projections/SmearedParticles.hh" #include "Rivet/Projections/SmearedMET.hh" namespace Rivet { class EXAMPLE_SMEAR : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(EXAMPLE_SMEAR); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { MissingMomentum mm(Cuts::abseta < 5); declare(mm, "MET0"); SmearedMET smm1(mm, MET_SMEAR_IDENTITY); declare(smm1, "MET1"); SmearedMET smm2(mm, [](const Vector3& met, double){ return P3_SMEAR_LEN_GAUSS(met, 0.1*met.mod()); }); declare(smm2, "MET2"); FastJets fj(FinalState(Cuts::abseta < 5), FastJets::ANTIKT, 0.4); declare(fj, "Jets0"); SmearedJets sj1(fj, JET_SMEAR_IDENTITY); declare(sj1, "Jets1"); SmearedJets sj2(fj, JET_SMEAR_ATLAS_RUN1, [](const Jet& j){ return j.bTagged() ? 0.7*(1 - exp(-j.pT()/(10*GeV))) : 0.01; } ); declare(sj2, "Jets2"); SmearedJets sj3(fj, [](const Jet& j){ return j; }, - [](const Jet& j){ return j.bTagged() ? 0.7*(1 - exp(-j.pT()/(10*GeV))) : 0.01; }, + JET_BTAG_EFFS(0.7, 0.1, 0.01), JET_CTAG_PERFECT, [](const Jet& j){ return 0.8; }); declare(sj3, "Jets3"); IdentifiedFinalState photons(Cuts::abseta < 5, PID::PHOTON); IdentifiedFinalState truthelectrons(Cuts::abseta < 5 && Cuts::pT > 10*GeV, {{PID::ELECTRON, PID::POSITRON}}); declare(truthelectrons, "Electrons0"); DressedLeptons dressedelectrons(photons, truthelectrons, 0.2); declare(dressedelectrons, "Electrons1"); SmearedParticles recoelectrons(truthelectrons, ELECTRON_EFF_ATLAS_RUN1, ELECTRON_SMEAR_ATLAS_RUN1); //< @note Can't use dressedelectrons yet... declare(recoelectrons, "Electrons2"); IdentifiedFinalState truthmuons(Cuts::abseta < 5 && Cuts::pT > 10*GeV, {{PID::MUON, PID::ANTIMUON}}); declare(truthmuons, "Muons0"); DressedLeptons dressedmuons(photons, truthmuons, 0.2); declare(dressedmuons, "Muons1"); SmearedParticles recomuons(truthmuons, MUON_EFF_ATLAS_RUN1, MUON_SMEAR_ATLAS_RUN1); //< @note Can't use dressedmuons yet... declare(recomuons, "Muons2"); TauFinder truthtaus(TauFinder::ANY, Cuts::abseta < 5 && Cuts::pT > 10*GeV); declare(truthtaus, "Taus0"); DressedLeptons dressedtaus(photons, truthtaus, 0.2); declare(dressedtaus, "Taus1"); SmearedParticles recotaus(truthtaus, TAU_EFF_ATLAS_RUN1, TAU_SMEAR_ATLAS_RUN1); //< @note Can't use dressedtaus yet... declare(recotaus, "Taus2"); _h_met_true = bookHisto1D("met_true", 30, 0.0, 120); _h_met_reco = bookHisto1D("met_reco", 30, 0.0, 120); _h_nj_true = bookHisto1D("jet_N_true", 10, -0.5, 9.5); _h_nj_reco = bookHisto1D("jet_N_reco", 10, -0.5, 9.5); _h_j1pt_true = bookHisto1D("jet_pt1_true", 30, 0.0, 120); _h_j1pt_reco = bookHisto1D("jet_pt1_reco", 30, 0.0, 120); _h_j1eta_true = bookHisto1D("jet_eta1_true", 20, -5.0, 5.0); _h_j1eta_reco = bookHisto1D("jet_eta1_reco", 20, -5.0, 5.0); _h_ne_true = bookHisto1D("elec_N_true", 5, -0.5, 4.5); _h_ne_reco = bookHisto1D("elec_N_reco", 5, -0.5, 4.5); _h_e1pt_true = bookHisto1D("elec_pt1_true", 30, 0, 120); _h_e1pt_reco = bookHisto1D("elec_pt1_reco", 30, 0, 120); _h_e1eta_true = bookHisto1D("elec_eta1_true", 20, -5.0, 5.0); _h_e1eta_reco = bookHisto1D("elec_eta1_reco", 20, -5.0, 5.0); _h_nm_true = bookHisto1D("muon_N_true", 5, -0.5, 4.5); _h_nm_reco = bookHisto1D("muon_N_reco", 5, -0.5, 4.5); _h_m1pt_true = bookHisto1D("muon_pt1_true", 30, 0, 120); _h_m1pt_reco = bookHisto1D("muon_pt1_reco", 30, 0, 120); _h_m1eta_true = bookHisto1D("muon_eta1_true", 20, -5.0, 5.0); _h_m1eta_reco = bookHisto1D("muon_eta1_reco", 20, -5.0, 5.0); _h_nt_true = bookHisto1D("tau_N_true", 5, -0.5, 4.5); _h_nt_reco = bookHisto1D("tau_N_reco", 5, -0.5, 4.5); _h_t1pt_true = bookHisto1D("tau_pt1_true", 30, 0, 120); _h_t1pt_reco = bookHisto1D("tau_pt1_reco", 30, 0, 120); _h_t1eta_true = bookHisto1D("tau_eta1_true", 20, -5.0, 5.0); _h_t1eta_reco = bookHisto1D("tau_eta1_reco", 20, -5.0, 5.0); } /// Perform the per-event analysis void analyze(const Event& event) { const double weight = event.weight(); const Vector3 met0 = apply(event, "MET0").vectorEt(); const Vector3 met1 = apply(event, "MET1").vectorEt(); const Vector3 met2 = apply(event, "MET2").vectorEt(); MSG_DEBUG("MET = " << met0.mod()/GeV << ", " << met1.mod()/GeV << ", " << met2.mod()/GeV << " GeV"); _h_met_true->fill(met0.mod()/GeV, weight); _h_met_reco->fill(met2.mod()/GeV, weight); const Jets jets0 = apply(event, "Jets0").jetsByPt(Cuts::pT > 10*GeV); const Jets jets1 = apply(event, "Jets1").jetsByPt(Cuts::pT > 10*GeV); const Jets jets2 = apply(event, "Jets2").jetsByPt(Cuts::pT > 10*GeV); const Jets jets3 = apply(event, "Jets3").jetsByPt(Cuts::pT > 10*GeV); MSG_DEBUG("Numbers of jets = " << jets0.size() << " true; " << jets1.size() << ", " << jets2.size() << ", " << jets3.size()); _h_nj_true->fill(jets0.size(), weight); _h_nj_reco->fill(jets2.size(), weight); if (!jets0.empty()) { _h_j1pt_true->fill(jets0.front().pT()/GeV, weight); _h_j1eta_true->fill(jets0.front().eta(), weight); } if (!jets2.empty()) { _h_j1pt_reco->fill(jets2.front().pT()/GeV, weight); _h_j1eta_reco->fill(jets2.front().eta(), weight); } const Particles& elecs1 = apply(event, "Electrons1").particlesByPt(); const Particles& elecs2 = apply(event, "Electrons2").particlesByPt(); MSG_DEBUG("Numbers of electrons = " << elecs1.size() << " true; " << elecs2.size() << " reco"); _h_ne_true->fill(elecs1.size(), weight); _h_ne_reco->fill(elecs2.size(), weight); if (!elecs1.empty()) { _h_e1pt_true->fill(elecs1.front().pT()/GeV, weight); _h_e1eta_true->fill(elecs1.front().eta(), weight); } if (!elecs2.empty()) { _h_e1pt_reco->fill(elecs2.front().pT()/GeV, weight); _h_e1eta_reco->fill(elecs2.front().eta(), weight); } const Particles& muons1 = apply(event, "Muons1").particlesByPt(); const Particles& muons2 = apply(event, "Muons2").particlesByPt(); MSG_DEBUG("Numbers of muons = " << muons1.size() << " true; " << muons2.size() << " reco"); _h_nm_true->fill(muons1.size(), weight); _h_nm_reco->fill(muons2.size(), weight); if (!muons1.empty()) { _h_m1pt_true->fill(muons1.front().pT()/GeV, weight); _h_m1eta_true->fill(muons1.front().eta(), weight); } if (!muons2.empty()) { _h_m1pt_reco->fill(muons2.front().pT()/GeV, weight); _h_m1eta_reco->fill(muons2.front().eta(), weight); } const Particles& taus1 = apply(event, "Taus1").particlesByPt(); const Particles& taus2 = apply(event, "Taus2").particlesByPt(); MSG_DEBUG("Numbers of taus = " << taus1.size() << " true; " << taus2.size() << " reco"); _h_nt_true->fill(taus1.size(), weight); _h_nt_reco->fill(taus2.size(), weight); if (!taus1.empty()) { _h_t1pt_true->fill(taus1.front().pT()/GeV, weight); _h_t1eta_true->fill(taus1.front().eta(), weight); } if (!taus2.empty()) { _h_t1pt_reco->fill(taus2.front().pT()/GeV, weight); _h_t1eta_reco->fill(taus2.front().eta(), weight); } } /// Normalise histograms etc., after the run void finalize() { normalize(_h_nj_true); normalize(_h_nj_reco); normalize(_h_j1pt_true, 1-_h_nj_true->bin(0).area()); normalize(_h_j1pt_reco, 1-_h_nj_reco->bin(0).area()); normalize(_h_j1eta_true, 1-_h_nj_true->bin(0).area()); normalize(_h_j1eta_reco, 1-_h_nj_reco->bin(0).area()); normalize(_h_ne_true); normalize(_h_ne_reco); normalize(_h_e1pt_true, 1-_h_ne_true->bin(0).area()); normalize(_h_e1pt_reco, 1-_h_ne_reco->bin(0).area()); normalize(_h_e1eta_true, 1-_h_ne_true->bin(0).area()); normalize(_h_e1eta_reco, 1-_h_ne_reco->bin(0).area()); normalize(_h_nm_true); normalize(_h_nm_reco); normalize(_h_m1pt_true, 1-_h_nm_true->bin(0).area()); normalize(_h_m1pt_reco, 1-_h_nm_reco->bin(0).area()); normalize(_h_m1eta_true, 1-_h_nm_true->bin(0).area()); normalize(_h_m1eta_reco, 1-_h_nm_reco->bin(0).area()); normalize(_h_nt_true); normalize(_h_nt_reco); normalize(_h_t1pt_true, 1-_h_nt_true->bin(0).area()); normalize(_h_t1pt_reco, 1-_h_nt_reco->bin(0).area()); normalize(_h_t1eta_true, 1-_h_nt_true->bin(0).area()); normalize(_h_t1eta_reco, 1-_h_nt_reco->bin(0).area()); } //@} private: /// @name Histograms //@{ Histo1DPtr _h_met_true, _h_met_reco; Histo1DPtr _h_nj_true, _h_nj_reco, _h_ne_true, _h_ne_reco, _h_nm_true, _h_nm_reco, _h_nt_true, _h_nt_reco; Histo1DPtr _h_j1pt_true, _h_j1pt_reco, _h_e1pt_true, _h_e1pt_reco, _h_m1pt_true, _h_m1pt_reco, _h_t1pt_true, _h_t1pt_reco; Histo1DPtr _h_j1eta_true, _h_j1eta_reco, _h_e1eta_true, _h_e1eta_reco, _h_m1eta_true, _h_m1eta_reco, _h_t1eta_true, _h_t1eta_reco; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(EXAMPLE_SMEAR); } diff --git a/include/Rivet/Tools/JetSmearingFunctions.hh b/include/Rivet/Tools/JetSmearingFunctions.hh --- a/include/Rivet/Tools/JetSmearingFunctions.hh +++ b/include/Rivet/Tools/JetSmearingFunctions.hh @@ -1,135 +1,135 @@ // -*- C++ -*- #ifndef RIVET_JetSmearingFunctions_HH #define RIVET_JetSmearingFunctions_HH #include "Rivet/Jet.hh" #include "Rivet/Tools/MomentumSmearingFunctions.hh" #include "Rivet/Tools/ParticleSmearingFunctions.hh" #include "Rivet/Tools/Random.hh" namespace Rivet { /// @name Jet filtering, efficiency and smearing utils //@{ /// @name Typedef for Jet smearing functions/functors typedef std::function JetSmearFn; /// @name Typedef for Jet efficiency functions/functors typedef std::function JetEffFn; /// Return a constant 0 given a Jet as argument inline double JET_EFF_ZERO(const Jet& p) { return 0; } /// Return a constant 1 given a Jet as argument inline double JET_EFF_ONE(const Jet& p) { return 1; } /// Take a Jet and return a constant efficiency struct JET_EFF_CONST { JET_EFF_CONST(double eff) : _eff(eff) {} double operator () (const Jet& ) const { return _eff; } double _eff; }; /// Return 1 if the given Jet contains a b, otherwise 0 inline double JET_BTAG_PERFECT(const Jet& j) { return j.bTagged() ? 1 : 0; } /// Return 1 if the given Jet contains a c, otherwise 0 inline double JET_CTAG_PERFECT(const Jet& j) { return j.cTagged() ? 1 : 0; } /// @brief b-tagging efficiency functor, for more readable b-tag effs and mistag rates /// Note several constructors, allowing for optional specification of charm, tau, and light jet mistag rates struct JET_BTAG_EFFS { - JET_BTAG_EFFS(double eff_b, double eff_light=0) : _eff_b(eff_b), _eff_c(-1), _eff_t(-1), _eff_l(eff_light) { } - JET_BTAG_EFFS(double eff_b, double eff_c, double eff_light=0) : _eff_b(eff_b), _eff_c(eff_c), _eff_t(-1), _eff_l(eff_light) { } - JET_BTAG_EFFS(double eff_b, double eff_c, double eff_tau, double eff_light=0) : _eff_b(eff_b), _eff_c(eff_c), _eff_t(eff_tau), _eff_l(eff_light) { } + JET_BTAG_EFFS(double eff_b, double eff_light=0) : _eff_b(eff_b), _eff_c(0), _eff_t(0), _eff_l(eff_light) { } + JET_BTAG_EFFS(double eff_b, double eff_c, double eff_light) : _eff_b(eff_b), _eff_c(eff_c), _eff_t(0), _eff_l(eff_light) { } + JET_BTAG_EFFS(double eff_b, double eff_c, double eff_tau, double eff_light) : _eff_b(eff_b), _eff_c(eff_c), _eff_t(eff_tau), _eff_l(eff_light) { } inline double operator () (const Jet& j) { if (j.bTagged()) return _eff_b; if (_eff_c >= 0 && j.cTagged()) return _eff_c; if (_eff_t >= 0 && j.tauTagged()) return _eff_t; return _eff_l; } double _eff_b, _eff_c, _eff_t, _eff_l; }; /// Take a jet and return an unmodified copy /// @todo Modify constituent particle vectors for consistency /// @todo Set a null PseudoJet if the Jet is smeared? inline Jet JET_SMEAR_IDENTITY(const Jet& j) { return j; } /// @brief Functor for simultaneous efficiency-filtering and smearing of Jets /// /// A central element of the SmearedJets system /// /// @todo Include tagging efficiency functions? struct JetEffSmearFn { JetEffSmearFn(const JetSmearFn& s, const JetEffFn& e) : sfn(s), efn(e) { } JetEffSmearFn(const JetEffFn& e, const JetSmearFn& s) : sfn(s), efn(e) { } JetEffSmearFn(const JetSmearFn& s) : sfn(s), efn(JET_EFF_ONE) { } JetEffSmearFn(const JetEffFn& e) : sfn(JET_SMEAR_IDENTITY), efn(e) { } JetEffSmearFn(double eff) : JetEffSmearFn(JET_EFF_CONST(eff)) { } /// Smear and calculate an efficiency for the given jet pair operator() (const Jet& j) const { return make_pair(sfn(j), efn(j)); } /// Compare to another, for use in the projection system int cmp(const JetEffSmearFn& other) const { // cout << "Eff hashes = " << get_address(efn) << "," << get_address(other.efn) << "; " // << "smear hashes = " << get_address(sfn) << "," << get_address(other.sfn) << endl; if (get_address(sfn) == 0 || get_address(other.sfn) == 0) return UNDEFINED; if (get_address(efn) == 0 || get_address(other.efn) == 0) return UNDEFINED; return Rivet::cmp(get_address(sfn), get_address(other.sfn)) || Rivet::cmp(get_address(efn), get_address(other.efn)); } /// Automatic conversion to a smearing function operator JetSmearFn () { return sfn; } /// Automatic conversion to an efficiency function /// @todo Ambiguity re. whether reco eff or a tagging efficiency... // operator JetEffFn () { return efn; } // Stored functions/functors JetSmearFn sfn; JetEffFn efn; }; /// Return true if Jet @a j is chosen to survive a random efficiency selection template inline bool efffilt(const Jet& j, FN& feff) { return rand01() < feff(j); } /// A functor to return true if Jet @a j survives a random efficiency selection struct JetEffFilter { template JetEffFilter(const FN& feff) : _feff(feff) {} JetEffFilter(double eff) : JetEffFilter( [&](const Jet& j){return eff;} ) {} bool operator () (const Jet& j) const { return efffilt(j, _feff); } private: const std::function _feff; }; using jetEffFilter = JetEffFilter; //@} } #endif