diff --git a/src/ScaleFunction.cc b/src/ScaleFunction.cc index e2af3f4..5747f97 100644 --- a/src/ScaleFunction.cc +++ b/src/ScaleFunction.cc @@ -1,149 +1,149 @@ /** * \authors Jeppe Andersen, Tuomas Hapola, Marian Heil, Andreas Maier, Jennifer Smillie * \date 2019 * \copyright GPLv2 or later */ #include "HEJ/ScaleFunction.hh" #include #include "HEJ/Event.hh" #include "HEJ/exceptions.hh" namespace HEJ{ double H_T(Event const & ev){ double result = 0.; for(size_t i = 0; i < ev.outgoing().size(); ++i){ auto const decay_products = ev.decays().find(i); if(decay_products == end(ev.decays())){ result += ev.outgoing()[i].perp(); } else{ for(auto const & particle: decay_products->second){ result += particle.perp(); } } } return result; } double max_jet_pt(Event const & ev) { return sorted_by_pt(ev.jets()).front().pt(); } double jet_invariant_mass(Event const & ev) { fastjet::PseudoJet sum; for(const auto & jet: ev.jets()) sum+=jet; return sum.m(); } double m_j1j2(Event const & ev) { const auto jets = sorted_by_pt(ev.jets()); assert(jets.size() >= 2); return (jets[0] + jets[1]).m(); } ScaleFunction operator*(double factor, ScaleFunction base_scale) { - base_scale.name_.insert(0, std::to_string(factor) + '_times_'); + base_scale.name_.insert(0, std::to_string(factor) + "_times_"); auto new_fun = [factor,fun{std::move(base_scale.fun_)}](HEJ::Event const & ev) { return factor*fun(ev); }; base_scale.fun_ = std::move(new_fun); return base_scale; } ScaleFunction operator/(ScaleFunction base_scale, double denom) { - base_scale.name_.append('_over_' + std::to_string(denom)); + base_scale.name_.append("_over_" + std::to_string(denom)); auto new_fun = [denom,fun{std::move(base_scale.fun_)}](HEJ::Event const & ev) { return fun(ev)/denom; }; base_scale.fun_ = std::move(new_fun); return base_scale; } // TODO: significant logic duplication with operator() void ScaleGenerator::gen_descriptions() { if(scales_.empty()) { throw std::logic_error{"Need at least one scale"}; } descriptions_.emplace_back( std::make_shared(scales_.front().name(), 1., 1.) ); for(auto & base_scale: scales_){ const auto base_name = base_scale.name(); descriptions_.emplace_back( std::make_shared(base_name, 1., 1.) ); //multiplicative scale variation for(double mur_factor: scale_factors_){ for(double muf_factor: scale_factors_){ if(muf_factor == 1. && mur_factor == 1.) continue; if( mur_factor/muf_factor < 1/max_scale_ratio_ || mur_factor/muf_factor > max_scale_ratio_ ) continue; descriptions_.emplace_back( std::make_shared( base_name, mur_factor, muf_factor ) ); } } } } Event ScaleGenerator::operator()(Event ev) const { if(! ev.variations().empty()) { throw std::invalid_argument{"Input event already has scale variation"}; } assert(!scales_.empty()); assert(!descriptions_.empty()); size_t descr_idx = 0; const double mu_central = (scales_.front())(ev); ev.central().mur = mu_central; ev.central().muf = mu_central; assert(descr_idx < descriptions_.size()); assert(descriptions_[descr_idx]); ev.central().description = descriptions_[descr_idx++]; // check if we are doing scale variation at all if(scales_.size() == 1 && scale_factors_.empty()) return ev; for(auto & base_scale: scales_){ const double mu_base = base_scale(ev); assert(descr_idx < descriptions_.size()); assert(descriptions_[descr_idx]); ev.variations().emplace_back( EventParameters{ mu_base, mu_base, ev.central().weight, descriptions_[descr_idx++] } ); //multiplicative scale variation for(double mur_factor: scale_factors_){ const double mur = mu_base*mur_factor; for(double muf_factor: scale_factors_){ if(muf_factor == 1. && mur_factor == 1.) continue; const double muf = mu_base*muf_factor; if( mur/muf < 1/max_scale_ratio_ || mur/muf > max_scale_ratio_ ) continue; assert(descr_idx < descriptions_.size()); assert(descriptions_[descr_idx]); ev.variations().emplace_back( EventParameters{ mur, muf, ev.central().weight, descriptions_[descr_idx++] } ); } } }; return ev; } }