diff --git a/FixedOrderGen/t/W_2j_classify.cc b/FixedOrderGen/t/W_2j_classify.cc index 9e7c2d8..525903b 100644 --- a/FixedOrderGen/t/W_2j_classify.cc +++ b/FixedOrderGen/t/W_2j_classify.cc @@ -1,114 +1,147 @@ #ifdef NDEBUG #undef NDEBUG #endif #include "PhaseSpacePoint.hh" #include "Process.hh" #include "JetParameters.hh" #include "HiggsProperties.hh" #include "RHEJ/Mixmax.hh" #include "RHEJ/PDF.hh" #include "RHEJ/utility.hh" using namespace HEJFOG; using namespace RHEJ; namespace { void print_event(PhaseSpacePoint const & psp){ std::cerr << "Process:\n" << psp.incoming()[0].type << " + "<< psp.incoming()[1].type << " -> "; for(auto const & out: psp.outgoing()){ std::cerr << out.type << " "; } std::cerr << "\n"; } void bail_out(PhaseSpacePoint const & psp, std::string msg){ print_event(psp); throw std::logic_error{msg}; } bool is_up_type(Particle const & part){ return RHEJ::is_anyquark(part) && !(abs(part.type)%2); } bool is_down_type(Particle const & part){ return RHEJ::is_anyquark(part) && abs(part.type)%2; } + + bool check_W2j(PhaseSpacePoint const & psp, ParticleID W_type){ + bool found_quark = false; + bool found_anti = false; + std::vector out_partons; + std::vector Wp; + for(auto const & p: psp.outgoing()){ + if(p.type == W_type) Wp.push_back(p); + else if(is_parton(p)) out_partons.push_back(p); + else bail_out(psp, "Found particle with is not " + +std::to_string(int(W_type))+" or parton"); + } + if(Wp.size() != 1 || out_partons.size() != 2){ + bail_out(psp, "Found wrong number of outgoing partons"); + } + for(size_t j=0; j<2; ++j){ + auto const & in = psp.incoming()[j]; + auto const & out = out_partons[j]; + if(is_quark(in) || is_antiquark(in)) { + found_quark = true; + if(in.type != out.type) { // switch in quark type -> Wp couples to it + if(found_anti){ // already found qq for coupling to W + bail_out(psp, "Found second up/down pair"); + } + found_anti = true; + if( is_up_type(in)) { // "up" in + if(W_type > 0){ + // -> only allowed u -> Wp + d + if(in.type < 0 || is_up_type(out) || out.type < 0) + bail_out(psp, "u -/> Wp + d"); + + } else { + // -> only allowed ux -> Wm + dx + if(in.type > 0 || is_up_type(out) || out.type > 0) + bail_out(psp, "ux -/> Wm + dx"); + } + } else { // "down" in + if(W_type > 0){ + // -> only allowed dx -> Wp + ux + if(in.type > 0 || is_down_type(out) || out.type > 0) + bail_out(psp, "dx -/> Wp + ux"); + } else { + // -> only allowed d -> Wm + u + if(in.type < 0 || is_down_type(out) || out.type < 0) + bail_out(psp, "d -/> Wm + u"); + } + } + } + } + } + if(!found_quark) { + bail_out(psp, "Found no initial quarks"); + } else if(!found_anti){ + bail_out(psp, "Found no up/down pair"); + } + return true; + } } int main(){ constexpr size_t n_psp_base = 1337; const JetParameters jet_para{ fastjet::JetDefinition(fastjet::JetAlgorithm::antikt_algorithm, 0.4), 20,5,30}; PDF pdf(11000, pid::proton, pid::proton); constexpr double E_cms = 13000.; constexpr double subl_change = 0.5; const HiggsProperties higgs_para{125., 0.001, {Decay{{pid::photon, pid::photon},1.}}}; RHEJ::Mixmax ran{}; - // W2j - const Process proc {{pid::proton,pid::proton}, 2, pid::Wp}; + // Wp2j + Process proc {{pid::proton,pid::proton}, 2, pid::Wp}; size_t n_psp = n_psp_base; size_t t_fail = 0; for( size_t i = 0; i out_partons; - std::vector Wp; - for(auto const & p: psp.outgoing()){ - if(p.type == pid::Wp) Wp.push_back(p); - else if(is_parton(p)) out_partons.push_back(p); - else bail_out(psp, "Found particle with is not Wp or parton"); - } - if(Wp.size() != 1 || out_partons.size() != 2){ - bail_out(psp, "Found wrong number of outgoing partons"); - } - for(size_t j=0; j<2; ++j){ - auto const & in = psp.incoming()[j]; - auto const & out = out_partons[j]; - if(is_quark(in) || is_antiquark(in)) { - found_quark = true; - if(in.type != out.type) { // switch in quark type -> Wp couples to it - if(found_anti){ // already found qq for coupling to W - bail_out(psp, "Found second up/down pair"); - } - found_anti = true; - if( is_up_type(in)) { // "up" in - // -> only allowed u -> Wp + d - if(in.type < 0 || is_up_type(out) || out.type < 0){ - // not in "q" or out "up" or out "qx" -> bail out - bail_out(psp, "u -/> Wp + d"); - } - } else { // "down" in - // -> only allowed dx -> Wp + ux - if(in.type > 0 || is_down_type(out) || out.type > 0){ - // not in "qx" or out "down" or out "q" -> bail out - bail_out(psp, "dx -/> Wp + ux"); - } - } - } - } - } - if(!found_quark) { - bail_out(psp, "Found no initial quarks"); - } else if(!found_anti){ - bail_out(psp, "Found no up/down pair"); - } + check_W2j(psp, *proc.boson); } else { // bad process -> try again ++n_psp; if(psp.status()==wrong_final_state) ++t_fail; } } - // @TODO test final state for 3j - - std::cout << "All processes passed.\nTook " << n_psp << " to generate " + std::cout << "Wp+2j: Took " << n_psp << " to generate " + << n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl; + std::cout << "Failed for wrong final state " << t_fail + << " (" << 100.*t_fail/(n_psp-n_psp_base) << " % of fails)" << std::endl; + // Wm2j + proc = Process{{pid::proton,pid::proton}, 2, pid::Wm}; + n_psp = n_psp_base; + t_fail = 0; + for( size_t i = 0; i try again + ++n_psp; + if(psp.status()==wrong_final_state) ++t_fail; + } + } + std::cout << "Wm+2j: Took " << n_psp << " to generate " << n_psp_base << " successfully PSP (" << 1.*n_psp/n_psp_base << " trials/PSP)" << std::endl; std::cout << "Failed for wrong final state " << t_fail << " (" << 100.*t_fail/(n_psp-n_psp_base) << " % of fails)" << std::endl; + // @TODO test final state for 3j + + std::cout << "All processes passed." << std::endl; return EXIT_SUCCESS; }