diff --git a/src/MatrixElement.cc b/src/MatrixElement.cc index 86c4539..a1d09d7 100644 --- a/src/MatrixElement.cc +++ b/src/MatrixElement.cc @@ -1,994 +1,994 @@ #include "RHEJ/MatrixElement.hh" #include <CLHEP/Random/Randomize.h> #include <CLHEP/Random/RanluxEngine.h> #include "RHEJ/Constants.hh" #include "RHEJ/currents.hh" #include "RHEJ/PDG_codes.hh" #include "RHEJ/uno.hh" #include "RHEJ/qqx.hh" #include "RHEJ/utility.hh" namespace RHEJ{ //cf. last line of eq. (22) in \ref Andersen:2011hs double MatrixElement::omega0( double alpha_s, double mur, fastjet::PseudoJet const & q_j, double lambda ) const { const double result = - alpha_s*N_C/M_PI*log(q_j.perp2()/(lambda*lambda)); if(! param_.log_correction) return result; // use alpha_s(sqrt(q_j*lambda)), evolved to mur return ( 1. + alpha_s/(4.*M_PI)*beta0*log(mur*mur/(q_j.perp()*lambda)) )*result; } double MatrixElement::virtual_corrections( double mur, std::array<Particle, 2> const & in, std::vector<Particle> const & out ) const{ fastjet::PseudoJet const & pa = in.front().p; #ifndef NDEBUG fastjet::PseudoJet const & pb = in.back().p; double const norm = (in.front().p + in.back().p).E(); #endif assert(std::is_sorted(out.begin(), out.end(), rapidity_less{})); assert(out.size() >= 2); assert(pa.pz() < pb.pz()); fastjet::PseudoJet q = pa - out[0].p; size_t first_idx = 0; size_t last_idx = out.size() - 1; // if there is a Higgs or unordered gluon outside the extremal partons // then it is not part of the FKL ladder and does not contribute // to the virtual corrections if(out.front().type == pid::Higgs || has_unob_gluon(in, out)){ q -= out[1].p; ++first_idx; } if(out.back().type == pid::Higgs || has_unof_gluon(in, out)){ --last_idx; } double exponent = 0; const double alpha_s = alpha_s_(mur); for(size_t j = first_idx; j < last_idx; ++j){ exponent += omega0(alpha_s, mur, q, CLAMBDA)*( out[j+1].rapidity() - out[j].rapidity() ); q -= out[j+1].p; } assert( nearby(q, -1*pb, norm) || out.back().type == pid::Higgs || has_unof_gluon(in, out) ); return exp(exponent); } } // namespace RHEJ namespace { //! Lipatov vertex for partons emitted into extremal jets double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2) { CLHEP::HepLorentzVector temptrans=-(qav+qbv); CLHEP::HepLorentzVector p5=qav-qbv; CLHEP::HepLorentzVector CL=temptrans + p1*(qav.m2()/p5.dot(p1) + 2.*p5.dot(p2)/p1.dot(p2)) - p2*(qbv.m2()/p5.dot(p2) + 2.*p5.dot(p1)/p1.dot(p2)); // cout << "#Fadin qa : "<<qav<<endl; // cout << "#Fadin qb : "<<qbv<<endl; // cout << "#Fadin p1 : "<<p1<<endl; // cout << "#Fadin p2 : "<<p2<<endl; // cout << "#Fadin p5 : "<<p5<<endl; // cout << "#Fadin Gauge Check : "<< CL.dot(p5)<<endl; // cout << "#Fadin C2L : "<< -CL.dot(CL)<<" "<<-CL.dot(CL)/(qav.m2()*qbv.m2())/(4./p5.perp2())<<endl; // TODO can this dead test go? // if (-CL.dot(CL)<0.) // if (fabs(CL.dot(p5))>fabs(CL.dot(CL))) // not sufficient! // return 0.; // else return -CL.dot(CL); } //! Lipatov vertex with soft subtraction for partons emitted into extremal jets double C2Lipatovots(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2) { double kperp=(qav-qbv).perp(); if (kperp>RHEJ::CLAMBDA) return C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2()); else { double Cls=(C2Lipatov(qav, qbv, p1, p2)/(qav.m2()*qbv.m2())); return Cls-4./(kperp*kperp); } } //! Lipatov vertex double C2Lipatov(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, CLHEP::HepLorentzVector pim, CLHEP::HepLorentzVector pip, CLHEP::HepLorentzVector pom, CLHEP::HepLorentzVector pop) // B { CLHEP::HepLorentzVector temptrans=-(qav+qbv); CLHEP::HepLorentzVector p5=qav-qbv; CLHEP::HepLorentzVector CL=temptrans + qav.m2()*(1./p5.dot(pip)*pip + 1./p5.dot(pop)*pop)/2. - qbv.m2()*(1./p5.dot(pim)*pim + 1./p5.dot(pom)*pom)/2. + ( pip*(p5.dot(pim)/pip.dot(pim) + p5.dot(pom)/pip.dot(pom)) + pop*(p5.dot(pim)/pop.dot(pim) + p5.dot(pom)/pop.dot(pom)) - pim*(p5.dot(pip)/pip.dot(pim) + p5.dot(pop)/pop.dot(pim)) - pom*(p5.dot(pip)/pip.dot(pom) + p5.dot(pop)/pop.dot(pom)) )/2.; return -CL.dot(CL); } //! Lipatov vertex with soft subtraction double C2Lipatovots(CLHEP::HepLorentzVector qav, CLHEP::HepLorentzVector qbv, CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector pb, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector p2) { double kperp=(qav-qbv).perp(); if (kperp>RHEJ::CLAMBDA) return C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2()); else { double Cls=(C2Lipatov(qav, qbv, pa, pb, p1, p2)/(qav.m2()*qbv.m2())); double temp=Cls-4./(kperp*kperp); return temp; } } /** Matrix element squared for tree-level current-current scattering * @param aptype Particle a PDG ID * @param bptype Particle b PDG ID * @param pn Particle n Momentum * @param pb Particle b Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @returns ME Squared for Tree-Level Current-Current Scattering */ double ME_current( int aptype, int bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa ){ if (aptype==21&&bptype==21) { return jM2gg(pn,pb,p1,pa); } else if (aptype==21&&bptype!=21) { if (bptype > 0) return jM2qg(pn,pb,p1,pa); else return jM2qbarg(pn,pb,p1,pa); } else if (bptype==21&&aptype!=21) { // ----- || ----- if (aptype > 0) return jM2qg(p1,pa,pn,pb); else return jM2qbarg(p1,pa,pn,pb); } else { // they are both quark if (bptype>0) { if (aptype>0) return jM2qQ(pn,pb,p1,pa); else return jM2qQbar(pn,pb,p1,pa); } else { if (aptype>0) return jM2qQbar(p1,pa,pn,pb); else return jM2qbarQbar(pn,pb,p1,pa); } } throw std::logic_error("unknown particle types"); } /** Matrix element squared for tree-level current-current scattering With W+Jets * @param aptype Particle a PDG ID * @param bptype Particle b PDG ID * @param pn Particle n Momentum * @param pb Particle b Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @returns ME Squared for Tree-Level Current-Current Scattering */ double ME_W_current( int aptype, int bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & plbar, CLHEP::HepLorentzVector const & pl ){ if (aptype==21&&bptype==21) { throw std::logic_error("gg incoming in W+jets, qqx not yet implemented"); } else if (aptype==21&&bptype!=21) { if (bptype > 0) return jMWqg(pn,pl,plbar,pb,p1,pa); else return jMWqbarg(pn,pl,plbar,pb,p1,pa); } else if (bptype==21&&aptype!=21) { // ----- || ----- if (aptype > 0) return jMWqg(p1,pl,plbar,pa,pn,pb); else return jMWqbarg(p1,pl,plbar,pa,pn,pb); } else { // they are both quark if (bptype>0) { if (aptype>0) return jMWqQ(pn,pl,plbar,pb,p1,pa); else return jMWqQbar(pn,pl,plbar,pb,p1,pa); } else { if (aptype>0) return jMWqQbar(p1,pl,plbar,pa,pn,pb); else return jMWqbarQbar(pn,pl,plbar,pb,p1,pa); } } throw std::logic_error("unknown particle types"); } /** \brief Matrix element squared for tree-level current-current scattering with Higgs * @param aptype Particle a PDG ID * @param bptype Particle b PDG ID * @param pn Particle n Momentum * @param pb Particle b Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @param qH t-channel momentum before Higgs * @param qHp1 t-channel momentum after Higgs * @returns ME Squared for Tree-Level Current-Current Scattering with Higgs */ double ME_Higgs_current( int aptype, int bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs double mt, bool include_bottom, double mb ){ if (aptype==21&&bptype==21) // gg initial state return MH2gg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else if (aptype==21&&bptype!=21) { if (bptype > 0) return MH2qg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.; else return MH2qbarg(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4./9.; } else if (bptype==21&&aptype!=21) { if (aptype > 0) return MH2qg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.; else return MH2qbarg(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4./9.; } else { // they are both quark if (bptype>0) { if (aptype>0) return MH2qQ(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.); else return MH2qQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.); } else { if (aptype>0) return MH2qQbar(p1,pa,pn,pb,-qH,-qHp1,mt,include_bottom,mb)*4.*4./(9.*9.); else return MH2qbarQbar(pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb)*4.*4./(9.*9.); } } throw std::logic_error("unknown particle types"); } /** \brief Current matrix element squared with Higgs and unordered forward emission * @param aptype Particle A PDG ID * @param bptype Particle B PDG ID * @param punof Unordered Particle Momentum * @param pn Particle n Momentum * @param pb Particle b Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @param qH t-channel momentum before Higgs * @param qHp1 t-channel momentum after Higgs * @returns ME Squared with Higgs and unordered forward emission */ double ME_Higgs_current_unof( int aptype, int bptype, CLHEP::HepLorentzVector const & punof, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs double mt, bool include_bottom, double mb ){ if (aptype==21&&bptype!=21) { if (bptype > 0) return jM2unogqHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unogqbarHg(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } else { // they are both quark if (bptype>0) { if (aptype>0) return jM2unogqHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unogqHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } else { if (aptype>0) return jM2unogqbarHQ(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unogqbarHQbar(punof,pn,pb,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } } throw std::logic_error("unknown particle types"); } /** \brief Current matrix element squared with Higgs and unordered backward emission * @param aptype Particle A PDG ID * @param bptype Particle B PDG ID * @param pn Particle n Momentum * @param pb Particle b Momentum * @param punob Unordered back Particle Momentum * @param p1 Particle 1 Momentum * @param pa Particle a Momentum * @param qH t-channel momentum before Higgs * @param qHp1 t-channel momentum after Higgs * @returns ME Squared with Higgs and unordered backward emission */ double ME_Higgs_current_unob( int aptype, int bptype, CLHEP::HepLorentzVector const & pn, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & punob, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & qH, // t-channel momentum before Higgs CLHEP::HepLorentzVector const & qHp1, // t-channel momentum after Higgs double mt, bool include_bottom, double mb ){ if (bptype==21&&aptype!=21) { if (aptype > 0) return jM2unobgHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unobgHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } else { // they are both quark if (aptype>0) { if (bptype>0) return jM2unobqHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unobqbarHQg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } else { if (bptype>0) return jM2unobqHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); else return jM2unobqbarHQbarg(pn,pb,punob,p1,pa,-qHp1,-qH,mt,include_bottom,mb); } } throw std::logic_error("unknown particle types"); } CLHEP::HepLorentzVector to_HepLorentzVector(RHEJ::Particle const & particle){ return {particle.p.px(), particle.p.py(), particle.p.pz(), particle.p.E()}; } void validate(RHEJ::MatrixElementConfig const & config) { #ifndef RHEJ_BUILD_WITH_QCDLOOP if(!config.Higgs_coupling.use_impact_factors) { throw std::invalid_argument{ "Invalid Higgs coupling settings.\n" "HEJ without QCDloop support can only use impact factors.\n" "Set use_impact_factors to true or recompile HEJ.\n" }; } #endif if(config.Higgs_coupling.use_impact_factors && config.Higgs_coupling.mt != std::numeric_limits<double>::infinity()) { throw std::invalid_argument{ "Conflicting settings: " "impact factors may only be used in the infinite top mass limit" }; } } } // namespace anonymous namespace RHEJ{ MatrixElement::MatrixElement( std::function<double (double)> alpha_s, MatrixElementConfig conf ): alpha_s_{std::move(alpha_s)}, param_{std::move(conf)} { validate(param_); } double MatrixElement::operator()( double mur, std::array<Particle, 2> const & incoming, std::vector<Particle> const & outgoing, std::unordered_map<int, std::vector<Particle>> const & decays, bool check_momenta ) const { return tree( mur, incoming, outgoing, decays, check_momenta )*virtual_corrections( mur, incoming, outgoing ); } double MatrixElement::tree_kin( std::array<Particle, 2> const & incoming, std::vector<Particle> const & outgoing, std::unordered_map<int, std::vector<Particle>> const & decays, bool check_momenta ) const { assert( std::is_sorted( incoming.begin(), incoming.end(), [](Particle o1, Particle o2){return o1.p.pz()<o2.p.pz();} ) ); assert(std::is_sorted(outgoing.begin(), outgoing.end(), rapidity_less{})); auto AWZH_boson = std::find_if( begin(outgoing), end(outgoing), [](Particle const & p){return is_AWZH_boson(p);} ); if(AWZH_boson == end(outgoing)){ return tree_kin_jets(incoming, outgoing, check_momenta); } switch(AWZH_boson->type){ case pid::Higgs: { return tree_kin_Higgs(incoming, outgoing, check_momenta); } // TODO case pid::Wp: { return tree_kin_W(incoming, outgoing, decays, true, check_momenta); } case pid::Wm: { return tree_kin_W(incoming, outgoing, decays, false, check_momenta); } case pid::photon: case pid::Z: default: throw std::logic_error("Emission of boson of unsupported type."); } } namespace{ constexpr int extremal_jet_idx = 1; constexpr int no_extremal_jet_idx = 0; bool treat_as_extremal(Particle const & parton){ return parton.p.user_index() == extremal_jet_idx; } template<class InputIterator> double FKL_ladder_weight( InputIterator begin_gluon, InputIterator end_gluon, CLHEP::HepLorentzVector const & q0, CLHEP::HepLorentzVector const & pa, CLHEP::HepLorentzVector const & pb, CLHEP::HepLorentzVector const & p1, CLHEP::HepLorentzVector const & pn ){ double wt = 1; auto qi = q0; for(auto gluon_it = begin_gluon; gluon_it != end_gluon; ++gluon_it){ assert(gluon_it->type == pid::gluon); const auto g = to_HepLorentzVector(*gluon_it); const auto qip1 = qi - g; if(treat_as_extremal(*gluon_it)){ wt *= C2Lipatovots(qip1, qi, pa, pb)*C_A; } else{ wt *= C2Lipatovots(qip1, qi, pa, pb, p1, pn)*C_A; } qi = qip1; } return wt; } } // namespace anonymous std::vector<Particle> MatrixElement::tag_extremal_jet_partons( std::array<Particle, 2> const & incoming, std::vector<Particle> out_partons, bool check_momenta ) const{ if(!check_momenta){ for(auto & parton: out_partons){ parton.p.set_user_index(no_extremal_jet_idx); } return out_partons; } fastjet::ClusterSequence cs(to_PseudoJet(out_partons), param_.jet_param.def); const auto jets = sorted_by_rapidity(cs.inclusive_jets(param_.jet_param.min_pt)); assert(jets.size() >= 2); auto most_backward = begin(jets); auto most_forward = end(jets) - 1; // skip jets caused by unordered emission if(has_unob_gluon(incoming, out_partons)){ assert(jets.size() >= 3); ++most_backward; } else if(has_unof_gluon(incoming, out_partons)){ assert(jets.size() >= 3); --most_forward; } const auto extremal_jet_indices = cs.particle_jet_indices( {*most_backward, *most_forward} ); assert(extremal_jet_indices.size() == out_partons.size()); for(size_t i = 0; i < out_partons.size(); ++i){ assert(RHEJ::is_parton(out_partons[i])); const int idx = (extremal_jet_indices[i]>=0)? extremal_jet_idx: no_extremal_jet_idx; out_partons[i].p.set_user_index(idx); } return out_partons; } double MatrixElement::tree_kin_jets( std::array<Particle, 2> const & incoming, std::vector<Particle> partons, bool check_momenta ) const { partons = tag_extremal_jet_partons(incoming, partons, check_momenta); if(has_unob_gluon(incoming, partons) || has_unof_gluon(incoming, partons)){ throw std::logic_error("unordered emission not implemented for pure jets"); } const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); const auto p1 = to_HepLorentzVector(partons.front()); const auto pn = to_HepLorentzVector(partons.back()); return ME_current( incoming[0].type, incoming[1].type, pn, pb, p1, pa )/(4*(N_C*N_C - 1))*FKL_ladder_weight( begin(partons) + 1, end(partons) - 1, pa - p1, pa, pb, p1, pn ); } double MatrixElement::tree_kin_W( std::array<Particle, 2> const & incoming, std::vector<Particle> const & outgoing, std::unordered_map<int, std::vector<Particle>> const & decays, bool WPlus, bool check_momenta ) const { if(has_unob_gluon(incoming, outgoing)){ throw std::logic_error("unordered emission not yet implemented for W+jets"); //return tree_kin_W_unob(incoming, outgoing, check_momenta); } else if(has_unof_gluon(incoming, outgoing)){ throw std::logic_error("unordered emission not yet implemented for W+jets"); // return tree_kin_W_unof(incoming, outgoing, check_momenta); } else if(has_Ex_qqx(incoming, outgoing)){ throw std::logic_error("Extremal qqx not yet implemented for W+jets"); // return tree_kin_W_Exqqx(incoming, outgoing, check_momenta); } else if(has_mid_qqx(outgoing)){ throw std::logic_error("Central qqx not yet implemented for W+jets"); // return tree_kin_W_qqxCentral(incoming, outgoing, check_momenta); } else{ return tree_kin_W_FKL(incoming, outgoing, decays, WPlus, check_momenta); } } double MatrixElement::tree_kin_W_FKL( std::array<Particle, 2> const & incoming, std::vector<Particle> const & outgoing, std::unordered_map<int, std::vector<Particle>> const & decays, bool WPlus, bool check_momenta ) const { const auto the_W = std::find_if( begin(outgoing), end(outgoing), [](Particle const & s){ return abs(s.type) == pid::Wp; } ); HLV plbar, pl; for (auto& x: decays) { if (x.second.at(0).type < 0){ plbar = to_HepLorentzVector(x.second.at(0)); pl = to_HepLorentzVector(x.second.at(1)); } else{ pl = to_HepLorentzVector(x.second.at(0)); plbar = to_HepLorentzVector(x.second.at(1)); } } const auto pW = to_HepLorentzVector(*the_W); std::vector<Particle> partons(begin(outgoing), the_W); partons.insert(end(partons), the_W + 1, end(outgoing)); partons = tag_extremal_jet_partons(incoming, partons, check_momenta); const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); auto p1 = to_HepLorentzVector(partons[0]); auto pn = to_HepLorentzVector(partons[partons.size() - 1]); auto first_after_W = begin(partons) + (the_W-begin(outgoing)); if(first_after_W == begin(partons)) ++first_after_W; else if(first_after_W == end(partons)) --first_after_W; // t-channel momentum before W auto qW = pa; for(auto parton_it = begin(partons); parton_it != first_after_W; ++parton_it){ qW -= to_HepLorentzVector(*parton_it); } auto q0 = pa - p1; auto begin_ladder = begin(partons) + 1; auto end_ladder = end(partons) - 1; double current_factor; if (WPlus){ current_factor = ME_W_current( incoming[0].type, incoming[1].type, pn, pb, p1, pa, pl, plbar ); } else{ current_factor = ME_W_current( incoming[0].type, incoming[1].type, pn, pb, p1, pa, plbar, pl ); } const double ladder_factor = FKL_ladder_weight( begin_ladder, first_after_W, q0, pa, pb, p1, pn )*FKL_ladder_weight( begin_ladder, end_ladder, pa - p1, pa, pb, p1, pn ); return current_factor*9./8.*ladder_factor; } double MatrixElement::tree_kin_Higgs( std::array<Particle, 2> const & incoming, std::vector<Particle> const & outgoing, bool check_momenta ) const { if(has_uno_gluon(incoming, outgoing)){ return tree_kin_Higgs_between(incoming, outgoing, check_momenta); } if(outgoing.front().type == pid::Higgs){ return tree_kin_Higgs_first(incoming, outgoing, check_momenta); } if(outgoing.back().type == pid::Higgs){ return tree_kin_Higgs_last(incoming, outgoing, check_momenta); } return tree_kin_Higgs_between(incoming, outgoing, check_momenta); } namespace { // Colour acceleration multipliers, for gluons see eq. (7) in arXiv:0910.5113 // TODO: code duplication with currents.cc double K_g(double p1minus, double paminus) { - return 1./2.*(p1minus/paminus + paminus/p1minus)*(C_A - 1/C_A) + 1/C_A; + return 1./2.*(p1minus/paminus + paminus/p1minus)*(C_A - 1./C_A) + 1./C_A; } double K_g( CLHEP::HepLorentzVector const & pout, CLHEP::HepLorentzVector const & pin ) { if(pin.z() > 0) return K_g(pout.plus(), pin.plus()); return K_g(pout.minus(), pin.minus()); } double K( ParticleID type, CLHEP::HepLorentzVector const & pout, CLHEP::HepLorentzVector const & pin ) { if(type == ParticleID::gluon) return K_g(pout, pin); return C_F; } // Colour factor in strict MRK limit double K_MRK(ParticleID type) { return (type == ParticleID::gluon)?C_A:C_F; } } double MatrixElement::MH2_forwardH( CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, ParticleID type2, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector pH, double t1, double t2 ) const{ ignore(p2out, p2in); const double shat = p1in.invariantMass2(p2in); // gluon case #ifdef RHEJ_BUILD_WITH_QCDLOOP if(!param_.Higgs_coupling.use_impact_factors){ return K(type2, p2out, p2in)*C_A*1./(16*M_PI*M_PI)*t1/t2*MH2gq_outsideH( p1out, p1in, p2out, p2in, pH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb )/(4*(N_C*N_C - 1)); } #endif return K_MRK(type2)/C_A*9./2.*shat*shat*( C2gHgp(p1in,p1out,pH) + C2gHgm(p1in,p1out,pH) )/(t1*t2); } double MatrixElement::tree_kin_Higgs_first( std::array<Particle, 2> const & incoming, std::vector<Particle> const & outgoing, bool check_momenta ) const { assert(outgoing.front().type == pid::Higgs); if(outgoing[1].type != pid::gluon) { assert(incoming.front().type == outgoing[1].type); return tree_kin_Higgs_between(incoming, outgoing, check_momenta); } const auto pH = to_HepLorentzVector(outgoing.front()); const auto partons = tag_extremal_jet_partons( incoming, std::vector<Particle>(begin(outgoing) + 1, end(outgoing)), check_momenta ); const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); const auto p1 = to_HepLorentzVector(partons.front()); const auto pn = to_HepLorentzVector(partons.back()); const auto q0 = pa - p1 - pH; const double t1 = q0.m2(); const double t2 = (pn - pb).m2(); return MH2_forwardH( p1, pa, incoming[1].type, pn, pb, pH, t1, t2 )*FKL_ladder_weight( begin(partons) + 1, end(partons) - 1, q0, pa, pb, p1, pn ); } double MatrixElement::tree_kin_Higgs_last( std::array<Particle, 2> const & incoming, std::vector<Particle> const & outgoing, bool check_momenta ) const { assert(outgoing.back().type == pid::Higgs); if(outgoing[outgoing.size()-2].type != pid::gluon) { assert(incoming.back().type == outgoing[outgoing.size()-2].type); return tree_kin_Higgs_between(incoming, outgoing, check_momenta); } const auto pH = to_HepLorentzVector(outgoing.back()); const auto partons = tag_extremal_jet_partons( incoming, std::vector<Particle>(begin(outgoing), end(outgoing) - 1), check_momenta ); const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); auto p1 = to_HepLorentzVector(partons.front()); const auto pn = to_HepLorentzVector(partons.back()); auto q0 = pa - p1; const double t1 = q0.m2(); const double t2 = (pn + pH - pb).m2(); return MH2_forwardH( pn, pb, incoming[0].type, p1, pa, pH, t2, t1 )*FKL_ladder_weight( begin(partons) + 1, end(partons) - 1, q0, pa, pb, p1, pn ); } double MatrixElement::tree_kin_Higgs_between( std::array<Particle, 2> const & incoming, std::vector<Particle> const & outgoing, bool check_momenta ) const { const auto the_Higgs = std::find_if( begin(outgoing), end(outgoing), [](Particle const & s){ return s.type == pid::Higgs; } ); assert(the_Higgs != end(outgoing)); const auto pH = to_HepLorentzVector(*the_Higgs); std::vector<Particle> partons(begin(outgoing), the_Higgs); partons.insert(end(partons), the_Higgs + 1, end(outgoing)); partons = tag_extremal_jet_partons(incoming, partons, check_momenta); const auto pa = to_HepLorentzVector(incoming[0]); const auto pb = to_HepLorentzVector(incoming[1]); auto p1 = to_HepLorentzVector( partons[has_unob_gluon(incoming, outgoing)?1:0] ); auto pn = to_HepLorentzVector( partons[partons.size() - (has_unof_gluon(incoming, outgoing)?2:1)] ); auto first_after_Higgs = begin(partons) + (the_Higgs-begin(outgoing)); assert( (first_after_Higgs == end(partons) && ( has_unob_gluon(incoming, outgoing) || partons.back().type != pid::gluon )) || first_after_Higgs->rapidity() >= the_Higgs->rapidity() ); assert( (first_after_Higgs == begin(partons) && ( has_unof_gluon(incoming, outgoing) || partons.front().type != pid::gluon )) || (first_after_Higgs-1)->rapidity() <= the_Higgs->rapidity() ); // always treat the Higgs as if it were in between the extremal FKL partons if(first_after_Higgs == begin(partons)) ++first_after_Higgs; else if(first_after_Higgs == end(partons)) --first_after_Higgs; // t-channel momentum before Higgs auto qH = pa; for(auto parton_it = begin(partons); parton_it != first_after_Higgs; ++parton_it){ qH -= to_HepLorentzVector(*parton_it); } auto q0 = pa - p1; auto begin_ladder = begin(partons) + 1; auto end_ladder = end(partons) - 1; double current_factor; if(has_unob_gluon(incoming, outgoing)){ current_factor = C_A*C_A/2.*ME_Higgs_current_unob( // 1/2 = "K_uno" incoming[0].type, incoming[1].type, pn, pb, to_HepLorentzVector(partons.front()), p1, pa, qH, qH - pH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb ); const auto p_unob = to_HepLorentzVector(partons.front()); q0 -= p_unob; p1 += p_unob; ++begin_ladder; } else if(has_unof_gluon(incoming, outgoing)){ current_factor = C_A*C_A/2.*ME_Higgs_current_unof( // 1/2 = "K_uno" incoming[0].type, incoming[1].type, to_HepLorentzVector(partons.back()), pn, pb, p1, pa, qH, qH - pH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb ); pn += to_HepLorentzVector(partons.back()); --end_ladder; } else{ current_factor = ME_Higgs_current( incoming[0].type, incoming[1].type, pn, pb, p1, pa, qH, qH - pH, param_.Higgs_coupling.mt, param_.Higgs_coupling.include_bottom, param_.Higgs_coupling.mb ); } const double ladder_factor = FKL_ladder_weight( begin_ladder, first_after_Higgs, q0, pa, pb, p1, pn )*FKL_ladder_weight( first_after_Higgs, end_ladder, qH - pH, pa, pb, p1, pn ); return current_factor*C_A*C_A/(N_C*N_C-1.)*ladder_factor; } double MatrixElement::tree_param_partons( double alpha_s, double mur, std::vector<Particle> const & partons ) const{ const double gs2 = 4.*M_PI*alpha_s; double wt = std::pow(gs2, partons.size()); if(param_.log_correction){ // use alpha_s(q_perp), evolved to mur assert(partons.size() >= 2); for(size_t i = 1; i < partons.size()-1; ++i){ wt *= 1 + alpha_s/(2*M_PI)*beta0*log(mur/partons[i].p.perp()); } } return wt; } double MatrixElement::tree_param( double mur, std::array<Particle, 2> const & incoming, std::vector<Particle> const & outgoing ) const{ const double alpha_s = alpha_s_(mur); auto AWZH_boson = std::find_if( begin(outgoing), end(outgoing), [](auto const & p){return is_AWZH_boson(p);} ); double AWZH_coupling = 1.; if(AWZH_boson != end(outgoing)){ switch(AWZH_boson->type){ case pid::Higgs: { AWZH_coupling = alpha_s*alpha_s; break; } // TODO case pid::Wp:{ AWZH_coupling = alpha_w*alpha_w/2; break; } case pid::Wm:{ AWZH_coupling = alpha_w*alpha_w/2; break; } case pid::photon: case pid::Z: default: throw std::logic_error("Emission of boson of unsupported type"); } } if(has_unob_gluon(incoming, outgoing)){ return AWZH_coupling*4*M_PI*alpha_s*tree_param_partons( alpha_s, mur, filter_partons({begin(outgoing) + 1, end(outgoing)}) ); } if(has_unof_gluon(incoming, outgoing)){ return AWZH_coupling*4*M_PI*alpha_s*tree_param_partons( alpha_s, mur, filter_partons({begin(outgoing), end(outgoing) - 1}) ); } return AWZH_coupling*tree_param_partons(alpha_s, mur, filter_partons(outgoing)); } double MatrixElement::tree( double mur, std::array<Particle, 2> const & incoming, std::vector<Particle> const & outgoing, std::unordered_map<int, std::vector<Particle>> const & decays, bool check_momenta ) const { return tree_param(mur, incoming, outgoing)*tree_kin( incoming, outgoing, decays, check_momenta ); } } // namespace RHEJ diff --git a/src/currents.cc b/src/currents.cc index c6377b3..9abcce8 100644 --- a/src/currents.cc +++ b/src/currents.cc @@ -1,3694 +1,3694 @@ ////////////////////////////////////////////////// ////////////////////////////////////////////////// // This source code is Copyright (2012) of // // Jeppe R. Andersen and Jennifer M. Smillie // // and is distributed under the // // Gnu Public License version 2 // // http://www.gnu.org/licenses/gpl-2.0.html // // You are allowed to distribute and alter the // // source under the conditions of the GPLv2 // // as long as this copyright notice // // is unaltered and distributed with the source // // Any use should comply with the // // MCNET GUIDELINES // // for Event Generator Authors and Users // // as distributed with this source code // ////////////////////////////////////////////////// ////////////////////////////////////////////////// #include "RHEJ/currents.hh" //#include "ZJets/Flags.h" #include "RHEJ/Constants.hh" #include "RHEJ/utility.hh" #include "RHEJ/PDG_codes.hh" const COM looprwfactor = (COM(0.,1.)*M_PI*M_PI)/pow((2.*M_PI),4); //const double HVE = 246.21845810181637; #ifdef RHEJ_BUILD_WITH_QCDLOOP #include "qcdloop/qcdloop.h" #endif #include <iostream> namespace { // Loop integrals #ifdef RHEJ_BUILD_WITH_QCDLOOP COM B0DD(CLHEP::HepLorentzVector q, double mq) { static std::vector<std::complex<double>> result(3); static auto ql_B0 = [](){ ql::Bubble<std::complex<double>,double,double> ql_B0; ql_B0.setCacheSize(100); return ql_B0; }(); static std::vector<double> masses(2); static std::vector<double> momenta(1); for(auto & m: masses) m = mq*mq; momenta.front() = q.m2(); ql_B0.integral(result, 1, masses, momenta); return result[0]; } COM C0DD(CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mq) { static std::vector<std::complex<double>> result(3); static auto ql_C0 = [](){ ql::Triangle<std::complex<double>,double,double> ql_C0; ql_C0.setCacheSize(100); return ql_C0; }(); static std::vector<double> masses(3); static std::vector<double> momenta(3); for(auto & m: masses) m = mq*mq; momenta[0] = q1.m2(); momenta[1] = q2.m2(); momenta[2] = (q1+q2).m2(); ql_C0.integral(result, 1, masses, momenta); return result[0]; } COM D0DD(CLHEP::HepLorentzVector q1,CLHEP::HepLorentzVector q2, CLHEP::HepLorentzVector q3, double mq) { static std::vector<std::complex<double>> result(3); static auto ql_D0 = [](){ ql::Box<std::complex<double>,double,double> ql_D0; ql_D0.setCacheSize(100); return ql_D0; }(); static std::vector<double> masses(4); static std::vector<double> momenta(6); for(auto & m: masses) m = mq*mq; momenta[0] = q1.m2(); momenta[1] = q2.m2(); momenta[2] = q3.m2(); momenta[3] = (q1+q2+q3).m2(); momenta[4] = (q1+q2).m2(); momenta[5] = (q2+q3).m2(); ql_D0.integral(result, 1, masses, momenta); return result[0]; } COM A1(CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt) // As given in Eq. (B.2) of VDD { double q12,q22,Q2; CLHEP::HepLorentzVector Q; double Delta3,mt2; COM ans(COM(0.,0.)); q12=q1.m2(); q22=q2.m2(); Q=-q1-q2; // Define all momenta ingoing as in appendix of VDD Q2=Q.m2(); // std::cout<<"Higgs mass? : "<<sqrt(Q2)<<std::endl; Delta3=q12*q12+q22*q22+Q2*Q2-2*q12*q22-2*q12*Q2-2*q22*Q2; if (mt < 0.) std::cerr<<"Problem in A1! mt = "<<mt<<std::endl; mt2=mt*mt; ans=looprwfactor*COM(0,-1)*C0DD(q1,q2,mt)*(4.*mt2/Delta3*(Q2-q12-q22)-1.-4.*q12*q22/Delta3-12.*q12*q22*Q2/Delta3/Delta3*(q12+q22-Q2)); ans=ans-looprwfactor*COM(0,-1)*(B0DD(q2,mt)-B0DD(Q,mt))*(2.*q22/Delta3+12.*q12*q22/Delta3/Delta3*(q22-q12+Q2)); ans=ans-looprwfactor*COM(0,-1)*(B0DD(q1,mt)-B0DD(Q,mt))*(2.*q12/Delta3+12.*q12*q22/Delta3/Delta3*(q12-q22+Q2)); ans=ans-2./Delta3/16/M_PI/M_PI*(q12+q22-Q2); //cout << "q12, q22= "<<q12<<" "<<q22<<" "<<endl; return ans; } COM A2(CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt) // As given in Eq. (B.2) of VDD, but with high energy limit // of invariants taken. { double q12,q22,Q2; CLHEP::HepLorentzVector Q; double Delta3,mt2; COM ans(COM(0.,0.)); if (mt < 0.) std::cerr<<"Problem in A2! mt = "<<mt<<std::endl; mt2=mt*mt; q12=q1.m2(); q22=q2.m2(); Q=-q1-q2; // Define all momenta ingoing as in appendix of VDD Q2=Q.m2(); // std::cout<<"Higgs mass Square? : "<<Q2<<std::endl; Delta3=q12*q12+q22*q22+Q2*Q2-2*q12*q22-2*q12*Q2-2*q22*Q2; ans=looprwfactor*COM(0,-1)*C0DD(q1,q2,mt)*(2.*mt2+1./2.*(q12+q22-Q2)+2.*q12*q22*Q2/Delta3); ans=ans+looprwfactor*COM(0,-1)*(B0DD(q2,mt)-B0DD(Q,mt))*q22*(q22-q12-Q2)/Delta3; ans=ans+looprwfactor*COM(0,-1)*(B0DD(q1,mt)-B0DD(Q,mt))*q12*(q12-q22-Q2)/Delta3+1./16/M_PI/M_PI; return ans; } #else // no QCDloop COM A1(CLHEP::HepLorentzVector, CLHEP::HepLorentzVector, double) { throw std::logic_error{"A1 called without QCDloop support"}; } COM A2(CLHEP::HepLorentzVector, CLHEP::HepLorentzVector, double) { throw std::logic_error{"A2 called without QCDloop support"}; } #endif void to_current(const CLHEP::HepLorentzVector & q, current & ret){ ret[0]=q.e(); ret[1]=q.x(); ret[2]=q.y(); ret[3]=q.z(); } constexpr double C_A = 3.; constexpr double C_F = 4./3.; using ParticleID = RHEJ::pid::ParticleID; // Colour acceleration multiplier for gluons see eq. (7) in arXiv:0910.5113 // TODO: this is not a current and should be moved somewhere else double K_g(double p1minus, double paminus) { - return 1./2.*(p1minus/paminus + paminus/p1minus)*(C_A - 1/C_A) + 1/C_A; + return 1./2.*(p1minus/paminus + paminus/p1minus)*(C_A - 1./C_A) + 1./C_A; } double K_g( CLHEP::HepLorentzVector const & pout, CLHEP::HepLorentzVector const & pin ) { if(pin.z() > 0) return K_g(pout.plus(), pin.plus()); return K_g(pout.minus(), pin.minus()); } } // namespace anonymous CCurrent CCurrent::operator+(const CCurrent& other) { COM result_c0=c0 + other.c0; COM result_c1=c1 + other.c1; COM result_c2=c2 + other.c2; COM result_c3=c3 + other.c3; return CCurrent(result_c0,result_c1,result_c2,result_c3); } CCurrent CCurrent::operator-(const CCurrent& other) { COM result_c0=c0 - other.c0; COM result_c1=c1 - other.c1; COM result_c2=c2 - other.c2; COM result_c3=c3 - other.c3; return CCurrent(result_c0,result_c1,result_c2,result_c3); } CCurrent CCurrent::operator*(const double x) { COM result_c0=x*CCurrent::c0; COM result_c1=x*CCurrent::c1; COM result_c2=x*CCurrent::c2; COM result_c3=x*CCurrent::c3; return CCurrent(result_c0,result_c1,result_c2,result_c3); } CCurrent CCurrent::operator/(const double x) { COM result_c0=CCurrent::c0/x; COM result_c1=CCurrent::c1/x; COM result_c2=CCurrent::c2/x; COM result_c3=CCurrent::c3/x; return CCurrent(result_c0,result_c1,result_c2,result_c3); } CCurrent CCurrent::operator*(const COM x) { COM result_c0=x*CCurrent::c0; COM result_c1=x*CCurrent::c1; COM result_c2=x*CCurrent::c2; COM result_c3=x*CCurrent::c3; return CCurrent(result_c0,result_c1,result_c2,result_c3); } CCurrent CCurrent::operator/(const COM x) { COM result_c0=(CCurrent::c0)/x; COM result_c1=(CCurrent::c1)/x; COM result_c2=(CCurrent::c2)/x; COM result_c3=(CCurrent::c3)/x; return CCurrent(result_c0,result_c1,result_c2,result_c3); } std::ostream& operator <<(std::ostream& os, const CCurrent& cur) { os << "("<<cur.c0<< " ; "<<cur.c1<<" , "<<cur.c2<<" , "<<cur.c3<<")"; return os; } CCurrent operator * ( double x, CCurrent& m) { return m*x; } CCurrent operator * ( COM x, CCurrent& m) { return m*x; } CCurrent operator / ( double x, CCurrent& m) { return m/x; } CCurrent operator / ( COM x, CCurrent& m) { return m/x; } COM CCurrent::dot(CLHEP::HepLorentzVector p1) { // Current goes (E,px,py,pz) // std::cout<<"current = ("<<c0<<","<<c1<<","<<c2<<","<<c3<<")\n"; // Vector goes (px,py,pz,E) // std::cout<<"vector = ("<<p1[0]<<","<<p1[1]<<","<<p1[2]<<","<<p1[3]<<")\n"; return p1[3]*c0-p1[0]*c1-p1[1]*c2-p1[2]*c3; } COM CCurrent::dot(CCurrent p1) { return p1.c0*c0-p1.c1*c1-p1.c2*c2-p1.c3*c3; } void j (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin,current &cur) { cur[0]=0.; cur[1]=0.; cur[2]=0.; cur[3]=0.; double sqpop=sqrt(pout.plus()); double sqpom=sqrt(pout.minus()); COM poperp=pout.x()+COM(0,1)*pout.y(); if (helout!=helin) { std::cerr<< "void j : Non-matching helicities at line " << __LINE__ << std::endl; } else if (helout==false) { // negative helicity if (pin.plus()>pin.minus()) { // if forward double sqpip=sqrt(pin.plus()); cur[0]=sqpop*sqpip; cur[1]=sqpom*sqpip*poperp/abs(poperp); cur[2]=-COM(0,1)*cur[1]; cur[3]=cur[0]; } else { // if backward double sqpim=sqrt(pin.minus()); cur[0]=-sqpom*sqpim*poperp/abs(poperp); cur[1]=-sqpim*sqpop; cur[2]=COM(0,1)*cur[1]; cur[3]=-cur[0]; } } else { // positive helicity if (pin.plus()>pin.minus()) { // if forward double sqpip=sqrt(pin.plus()); cur[0]=sqpop*sqpip; cur[1]=sqpom*sqpip*conj(poperp)/abs(poperp); cur[2]=COM(0,1)*cur[1]; cur[3]=cur[0]; } else { // if backward double sqpim=sqrt(pin.minus()); cur[0]=-sqpom*sqpim*conj(poperp)/abs(poperp); cur[1]=-sqpim*sqpop; cur[2]=-COM(0,1)*cur[1]; cur[3]=-cur[0]; } } } CCurrent j (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin) { COM cur[4]; cur[0]=0.; cur[1]=0.; cur[2]=0.; cur[3]=0.; double sqpop=sqrt(pout.plus()); double sqpom=sqrt(pout.minus()); COM poperp=pout.x()+COM(0,1)*pout.y(); if (helout!=helin) { std::cerr<< "void j : Non-matching helicities\n"; } else if (helout==false) { // negative helicity if (pin.plus()>pin.minus()) { // if forward double sqpip=sqrt(pin.plus()); cur[0]=sqpop*sqpip; cur[1]=sqpom*sqpip*poperp/abs(poperp); cur[2]=-COM(0,1)*cur[1]; cur[3]=cur[0]; } else { // if backward double sqpim=sqrt(pin.minus()); cur[0]=-sqpom*sqpim*poperp/abs(poperp); cur[1]=-sqpim*sqpop; cur[2]=COM(0,1)*cur[1]; cur[3]=-cur[0]; } } else { // positive helicity if (pin.plus()>pin.minus()) { // if forward double sqpip=sqrt(pin.plus()); cur[0]=sqpop*sqpip; cur[1]=sqpom*sqpip*conj(poperp)/abs(poperp); cur[2]=COM(0,1)*cur[1]; cur[3]=cur[0]; } else { // if backward double sqpim=sqrt(pin.minus()); cur[0]=-sqpom*sqpim*conj(poperp)/abs(poperp); cur[1]=-sqpim*sqpop; cur[2]=-COM(0,1)*cur[1]; cur[3]=-cur[0]; } } CCurrent temp(cur[0],cur[1],cur[2],cur[3]); return temp; } CCurrent jio (CLHEP::HepLorentzVector pin, bool helin, CLHEP::HepLorentzVector pout, bool helout) { COM cur[4]; cur[0]=0.; cur[1]=0.; cur[2]=0.; cur[3]=0.; double sqpop=sqrt(pout.plus()); double sqpom=sqrt(pout.minus()); COM poperp=pout.x()+COM(0,1)*pout.y(); if (helout!=helin) { std::cerr<< "void j : Non-matching helicities\n"; } else if (helout==false) { // negative helicity if (pin.plus()>pin.minus()) { // if forward double sqpip=sqrt(pin.plus()); cur[0]=sqpop*sqpip; cur[1]=sqpom*sqpip*conj(poperp)/abs(poperp); cur[2]=COM(0,1)*cur[1]; cur[3]=cur[0]; } else { // if backward double sqpim=sqrt(pin.minus()); cur[0]=-sqpom*sqpim*conj(poperp)/abs(poperp); cur[1]=-sqpim*sqpop; cur[2]=-COM(0,1)*cur[1]; cur[3]=-cur[0]; } } else { // positive helicity if (pin.plus()>pin.minus()) { // if forward double sqpip=sqrt(pin.plus()); cur[0]=sqpop*sqpip; cur[1]=sqpom*sqpip*poperp/abs(poperp); cur[2]=-COM(0,1)*cur[1]; cur[3]=cur[0]; } else { // if backward double sqpim=sqrt(pin.minus()); cur[0]=-sqpom*sqpim*poperp/abs(poperp); cur[1]=-sqpim*sqpop; cur[2]=COM(0,1)*cur[1]; cur[3]=-cur[0]; } } CCurrent temp(cur[0],cur[1],cur[2],cur[3]); return temp; } // Current for <incoming state | mu | outgoing state> void jio(HLV pin, bool helin, HLV pout, bool helout, current &cur) { cur[0] = 0.0; cur[1] = 0.0; cur[2] = 0.0; cur[3] = 0.0; if(helin!=helout){ std::cout<<__LINE__<<" "<<__FILE__<<std::endl; } double sqpop = sqrt(pout.plus()); double sqpom = sqrt(pout.minus()); COM poperp = pout.x() + COM(0, 1) * pout.y(); if (helout == false) { if (pin.plus() > pin.minus()) { // if forward double sqpip=sqrt(pin.plus()); cur[0] = sqpop * sqpip; cur[1] = sqpom * sqpip * conj(poperp) / abs(poperp); cur[2] = COM(0,1) * cur[1]; cur[3] = cur[0]; } else { double sqpim = sqrt(pin.minus()); cur[0] = -sqpom * sqpim * conj(poperp) / abs(poperp); cur[1] = -sqpim * sqpop; cur[2] = -COM(0,1) * cur[1]; cur[3] = -cur[0]; } } else { if (pin.plus() > pin.minus()) { // if forward double sqpip = sqrt(pin.plus()); cur[0] = sqpop * sqpip; cur[1] = sqpom * sqpip*poperp/abs(poperp); cur[2] = -COM(0,1)*cur[1]; cur[3] = cur[0]; } else { double sqpim = sqrt(pin.minus()); cur[0] = -sqpom * sqpim * poperp/abs(poperp); cur[1] = -sqpim * sqpop; cur[2] = COM(0,1)*cur[1]; cur[3] = -cur[0]; } } } // Current for <outgoing state | mu | outgoing state> void joo(HLV pi, bool heli, HLV pj, bool helj, current &cur) { // Zero our current cur[0] = 0.0; cur[1] = 0.0; cur[2] = 0.0; cur[3] = 0.0; if(helj){ std::cout<<__LINE__<<" "<<__FILE__<<std::endl; } // If positive helicity swap momenta if (heli == true) { HLV dummy; dummy = pi; pi = pj; pj = dummy; } double sqpjp = sqrt(pj.plus()); double sqpjm = sqrt(pj.minus()); double sqpip = sqrt(pi.plus()); double sqpim = sqrt(pi.minus()); COM piperp = pi.x() + COM(0,1) * pi.y(); COM pjperp = pj.x() + COM(0,1) * pj.y(); COM phasei = piperp / abs(piperp); COM phasej = pjperp / abs(pjperp); cur[0] = sqpim * sqpjm * phasei * conj(phasej) + sqpip * sqpjp; cur[1] = sqpim * sqpjp * phasei + sqpip * sqpjm * conj(phasej); cur[2] = -COM(0, 1) * (sqpim * sqpjp * phasei - sqpip * sqpjm * conj(phasej)); cur[3] = -sqpim * sqpjm * phasei * conj(phasej) + sqpip * sqpjp; } CCurrent joo (CLHEP::HepLorentzVector pi, bool heli, CLHEP::HepLorentzVector pj, bool helj) { COM cur[4]; if (heli!=helj) { std::cerr<< "void j : Non-matching helicities\n"; } else if (heli==true) { // negative helicity CLHEP::HepLorentzVector dummy; dummy=pi; pi=pj; pj=dummy; } double sqpjp=sqrt(pj.plus()); double sqpjm=sqrt(pj.minus()); double sqpip=sqrt(pi.plus()); double sqpim=sqrt(pi.minus()); COM piperp=pi.x()+COM(0,1)*pi.y(); COM pjperp=pj.x()+COM(0,1)*pj.y(); COM phasei=piperp/abs(piperp); COM phasej=pjperp/abs(pjperp); cur[0]=sqpim*sqpjm*phasei*conj(phasej)+sqpip*sqpjp; cur[1]=sqpim*sqpjp*phasei+sqpip*sqpjm*conj(phasej); cur[2]=-COM(0,1)*(sqpim*sqpjp*phasei-sqpip*sqpjm*conj(phasej)); cur[3]=-sqpim*sqpjm*phasei*conj(phasej)+sqpip*sqpjp; CCurrent temp(cur[0],cur[1],cur[2],cur[3]); return temp; } // Current Functions // Current for <outgoing state | mu | incoming state> void joi(HLV pout, bool helout, HLV pin, bool helin, current &cur) { cur[0] = 0.0; cur[1] = 0.0; cur[2] = 0.0; cur[3] = 0.0; if(helin){ std::cout<<__LINE__<<" "<<__FILE__<<std::endl; } double sqpop = sqrt(pout.plus()); double sqpom = sqrt(pout.minus()); COM poperp = pout.x() + COM(0, 1) * pout.y(); if (helout == false) { if (pin.plus() > pin.minus()) { // if forward double sqpip=sqrt(pin.plus()); cur[0] = sqpop * sqpip; cur[1] = sqpom * sqpip * poperp/abs(poperp); cur[2] = -COM(0,1)*cur[1]; cur[3] = cur[0]; } else { double sqpim = sqrt(pin.minus()); cur[0] = -sqpom*sqpim*poperp/abs(poperp); cur[1] = -sqpim*sqpop; cur[2] = COM(0,1)*cur[1]; cur[3] = -cur[0]; } } else { if (pin.plus() > pin.minus()) { // if forward double sqpip = sqrt(pin.plus()); cur[0] = sqpop * sqpip; cur[1] = sqpom * sqpip*conj(poperp)/abs(poperp); cur[2] = COM(0,1)*cur[1]; cur[3] = cur[0]; } else { double sqpim = sqrt(pin.minus()); cur[0] = -sqpom * sqpim * conj(poperp)/abs(poperp); cur[1] = -sqpim * sqpop; cur[2] = -COM(0,1)*cur[1]; cur[3] = -cur[0]; } } } namespace { /// @TODO unused function // double jM2 (CLHEP::HepLorentzVector p1out, bool hel1out, CLHEP::HepLorentzVector p1in, bool hel1in, CLHEP::HepLorentzVector p2out, bool hel2out, CLHEP::HepLorentzVector p2in, bool hel2in) // { // CLHEP::HepLorentzVector q1=p1in-p1out; // CLHEP::HepLorentzVector q2=-(p2in-p2out); // current C1,C2; // j (p1out,hel1out,p1in,hel1in, C1); // j (p2out,hel2out,p2in,hel2in, C2); // std::cout << "# From Currents, C1 : ("<<C1[0]<<","<<C1[1]<<","<<C1[2]<<","<<C1[3]<<"\n"; // std::cout << "# From Currents, C2 : ("<<C2[0]<<","<<C2[1]<<","<<C2[2]<<","<<C2[3]<<"\n"; // COM M=cdot(C1,C2); // return (M*conj(M)).real()/(q1.m2()*q2.m2()); // } void jW (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pe, bool hele, CLHEP::HepLorentzVector pnu, bool helnu, CLHEP::HepLorentzVector pin, bool helin, current cur) { // NOTA BENE: Conventions for W+ --> e+ nu, so that nu is lepton(6), e is anti-lepton(5) // Need to swap e and nu for events with W- --> e- nubar! if (helin==helout && hele==helnu) { CLHEP::HepLorentzVector qa=pout+pe+pnu; CLHEP::HepLorentzVector qb=pin-pe-pnu; double ta(qa.m2()),tb(qb.m2()); current t65,vout,vin,temp2,temp3,temp5; joo(pnu,helnu,pe,hele,t65); vout[0]=pout.e(); vout[1]=pout.x(); vout[2]=pout.y(); vout[3]=pout.z(); vin[0]=pin.e(); vin[1]=pin.x(); vin[2]=pin.y(); vin[3]=pin.z(); COM brac615=cdot(t65,vout); COM brac645=cdot(t65,vin); // prod1565 and prod6465 are zero for Ws (not Zs)!! // noalias(temp)=prod(trans(CurrentOutOut(pout,helout,pnu,helout)),metric); joo(pout,helout,pnu,helout,temp2); // noalias(temp2)=prod(temp,ctemp); COM prod1665=cdot(temp2,t65); // noalias(temp)=prod(trans(Current(pe,helin,pin,helin)),metric); // noalias(temp2)=prod(temp,ctemp); j(pe,helin,pin,helin,temp3); COM prod5465=cdot(temp3,t65); // noalias(temp)=prod(trans(Current(pnu,helin,pin,helin)),metric); // noalias(temp2)=prod(temp,ctemp); joo(pout,helout,pe,helout,temp2); j(pnu,helnu,pin,helin,temp3); j(pout,helout,pin,helin,temp5); current term1,term2,term3,sum; cmult(2.*brac615/ta+2.*brac645/tb,temp5,term1); cmult(prod1665/ta,temp3,term2); cmult(-prod5465/tb,temp2,term3); // cur=((2.*brac615*Current(pout,helout,pin,helin)+prod1565*Current(pe,helin,pin,helin)+prod1665*Current(pnu,helin,pin,helin))/ta + (2.*brac645*Current(pout,helout,pin,helin)-prod5465*CurrentOutOut(pout,helout,pe,helout)-prod6465*CurrentOutOut(pout,helout,pnu,helout))/tb); // cur=((2.*brac615*temp5+prod1565*temp3+prod1665*temp4)/ta + (2.*brac645*temp5-prod5465*temp1-prod6465*temp2)/tb); cadd(term1,term2,term3,sum); // std::cout<<"sum: ("<<sum[0]<<","<<sum[1]<<","<<sum[2]<<","<<sum[3]<<")\n"; cur[0]=sum[0]; cur[1]=sum[1]; cur[2]=sum[2]; cur[3]=sum[3]; } } void jWbar (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pe, bool hele, CLHEP::HepLorentzVector pnu, bool helnu, CLHEP::HepLorentzVector pin, bool helin, current cur) { // NOTA BENE: Conventions for W+ --> e+ nu, so that nu is lepton(6), e is anti-lepton(5) // Need to swap e and nu for events with W- --> e- nubar! if (helin==helout && hele==helnu) { CLHEP::HepLorentzVector qa=pout+pe+pnu; CLHEP::HepLorentzVector qb=pin-pe-pnu; double ta(qa.m2()),tb(qb.m2()); current t65,vout,vin,temp2,temp3,temp5; joo(pnu,helnu,pe,hele,t65); vout[0]=pout.e(); vout[1]=pout.x(); vout[2]=pout.y(); vout[3]=pout.z(); vin[0]=pin.e(); vin[1]=pin.x(); vin[2]=pin.y(); vin[3]=pin.z(); COM brac615=cdot(t65,vout); COM brac645=cdot(t65,vin); // prod1565 and prod6465 are zero for Ws (not Zs)!! joo(pe,helout,pout,helout,temp2); // temp2 is <5|alpha|1> COM prod5165=cdot(temp2,t65); jio(pin,helin,pnu,helin,temp3); // temp3 is <4|alpha|6> COM prod4665=cdot(temp3,t65); joo(pnu,helout,pout,helout,temp2); // temp2 is now <6|mu|1> jio(pin,helin,pe,helin,temp3); // temp3 is now <4|mu|5> jio(pin,helin,pout,helout,temp5); // temp5 is <4|mu|1> current term1,term2,term3,sum; cmult(-2.*brac615/ta-2.*brac645/tb,temp5,term1); cmult(-prod5165/ta,temp3,term2); cmult(prod4665/tb,temp2,term3); // cur=((2.*brac615*Current(pout,helout,pin,helin)+prod1565*Current(pe,helin,pin,helin)+prod1665*Current(pnu,helin,pin,helin))/ta + (2.*brac645*Current(pout,helout,pin,helin)-prod5465*CurrentOutOut(pout,helout,pe,helout)-prod6465*CurrentOutOut(pout,helout,pnu,helout))/tb); // cur=((2.*brac615*temp5+prod1565*temp3+prod1665*temp4)/ta + (2.*brac645*temp5-prod5465*temp1-prod6465*temp2)/tb); cadd(term1,term2,term3,sum); // std::cout<<"term1: ("<<temp5[0]<<" "<<temp5[1]<<" "<<temp5[2]<<" "<<temp5[3]<<")"<<std::endl; // std::cout<<"sum: ("<<sum[0]<<","<<sum[1]<<","<<sum[2]<<","<<sum[3]<<")\n"; cur[0]=sum[0]; cur[1]=sum[1]; cur[2]=sum[2]; cur[3]=sum[3]; } } } // namespace anonymous double jMWqQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in) // Calculates the square of the current contractions for qQ->qenuQ scattering // p1: quark (with W emittance) // p2: Quark { current mj1m,mj2p,mj2m; CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu; CLHEP::HepLorentzVector q2=-(p2in-p2out); jW(p1out,false,pe,false,pnu,false,p1in,false,mj1m); j(p2out,true,p2in,true,mj2p); j(p2out,false,p2in,false,mj2m); // std::cout<<"jMW1: ("<<mj1m[0]<<","<<mj1m[1]<<","<<mj1m[2]<<","<<mj1m[3]<<")\n"; // std::cout<<"jMW2: ("<<mj2p[0]<<","<<mj2p[1]<<","<<mj2p[2]<<","<<mj2p[3]<<")\n"; // std::cout<<"jMW3: ("<<mj2m[0]<<","<<mj2m[1]<<","<<mj2m[2]<<","<<mj2m[3]<<")\n"; // mj1m.mj2p COM Mmp=cdot(mj1m,mj2p); // mj1m.mj2m COM Mmm=cdot(mj1m,mj2m); // sum of spinor strings ||^2 double a2Mmp=abs2(Mmp); double a2Mmm=abs2(Mmm); // // Leave division by colour and Helicity avg until Tree files // Leave multi. of couplings to later // Multiply by Cf^2 return (4./3.)*(4./3.)*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2()); } double jMWqQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in) // Calculates the square of the current contractions for qQ->qenuQ scattering // p1: quark (with W emittance) // p2: Quark { current mj1m,mj2p,mj2m; CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu; CLHEP::HepLorentzVector q2=-(p2in-p2out); jW(p1out,false,pe,false,pnu,false,p1in,false,mj1m); jio(p2in,true,p2out,true,mj2p); jio(p2in,false,p2out,false,mj2m); // std::cout<<"jMW1: ("<<mj1m[0]<<","<<mj1m[1]<<","<<mj1m[2]<<","<<mj1m[3]<<")\n"; // std::cout<<"jMW2: ("<<mj2p[0]<<","<<mj2p[1]<<","<<mj2p[2]<<","<<mj2p[3]<<")\n"; // std::cout<<"jMW3: ("<<mj2m[0]<<","<<mj2m[1]<<","<<mj2m[2]<<","<<mj2m[3]<<")\n"; // mj1m.mj2p COM Mmp=cdot(mj1m,mj2p); // mj1m.mj2m COM Mmm=cdot(mj1m,mj2m); // sum of spinor strings ||^2 double a2Mmp=abs2(Mmp); double a2Mmm=abs2(Mmm); // // Leave division by colour and Helicity avg until Tree files // Leave multi. of couplings to later // Multiply by Cf^2 return (4./3.)*(4./3.)*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2()); } double jMWqbarQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in) // Calculates the square of the current contractions for qQ->qenuQ scattering // p1: quark (with W emittance) // p2: Quark { current mj1m,mj2p,mj2m; CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu; CLHEP::HepLorentzVector q2=-(p2in-p2out); jWbar(p1out,false,pe,false,pnu,false,p1in,false,mj1m); j(p2out,true,p2in,true,mj2p); j(p2out,false,p2in,false,mj2m); // std::cout<<"jMW1: ("<<mj1m[0]<<","<<mj1m[1]<<","<<mj1m[2]<<","<<mj1m[3]<<")\n"; // std::cout<<"jMW2: ("<<mj2p[0]<<","<<mj2p[1]<<","<<mj2p[2]<<","<<mj2p[3]<<")\n"; // std::cout<<"jMW3: ("<<mj2m[0]<<","<<mj2m[1]<<","<<mj2m[2]<<","<<mj2m[3]<<")\n"; // mj1m.mj2p COM Mmp=cdot(mj1m,mj2p); // mj1m.mj2m COM Mmm=cdot(mj1m,mj2m); // sum of spinor strings ||^2 double a2Mmp=abs2(Mmp); double a2Mmm=abs2(Mmm); // // Leave division by colour and Helicity avg until Tree files // Leave multi. of couplings to later // Multiply by Cf^2 return (4./3.)*(4./3.)*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2()); } double jMWqbarQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in) // Calculates the square of the current contractions for qQ->qenuQ scattering // p1: quark (with W emittance) // p2: Quark { current mj1m,mj2p,mj2m; CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu; CLHEP::HepLorentzVector q2=-(p2in-p2out); jWbar(p1out,false,pe,false,pnu,false,p1in,false,mj1m); jio(p2in,true,p2out,true,mj2p); jio(p2in,false,p2out,false,mj2m); // std::cout<<"jMW1: ("<<mj1m[0]<<","<<mj1m[1]<<","<<mj1m[2]<<","<<mj1m[3]<<")\n"; // std::cout<<"jMW2: ("<<mj2p[0]<<","<<mj2p[1]<<","<<mj2p[2]<<","<<mj2p[3]<<")\n"; // std::cout<<"jMW3: ("<<mj2m[0]<<","<<mj2m[1]<<","<<mj2m[2]<<","<<mj2m[3]<<")\n"; // mj1m.mj2p COM Mmp=cdot(mj1m,mj2p); // mj1m.mj2m COM Mmm=cdot(mj1m,mj2m); // sum of spinor strings ||^2 double a2Mmp=abs2(Mmp); double a2Mmm=abs2(Mmm); // // Leave division by colour and Helicity avg until Tree files // Leave multi. of couplings to later // Multiply by Cf^2 return (4./3.)*(4./3.)*(a2Mmp+a2Mmm)/(q1.m2()*q2.m2()); } double jMWqg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in) // Calculates the square of the current contractions for qg->qenug scattering // p1: quark // p2: gluon { CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu; CLHEP::HepLorentzVector q2=-(p2in-p2out); current mj1m,mj2p,mj2m; jW(p1out,false,pe,false,pnu,false,p1in,false,mj1m); j(p2out,true,p2in,true,mj2p); j(p2out,false,p2in,false,mj2m); // mj1m.mj2p COM Mmp=cdot(mj1m,mj2p); // mj1m.mj2m COM Mmm=cdot(mj1m,mj2m); const double K = K_g(p2out, p2in); // sum of spinor strings ||^2 double a2Mmp=abs2(Mmp); double a2Mmm=abs2(Mmm); double sst = K/C_A*(a2Mmp+a2Mmm); // double sstsave=sst; // // Leave division by colour and Helicity avg until Tree files // Leave multi. of couplings to later // Multiply by Cf*Ca=4 return 4.*sst/(q1.m2()*q2.m2()); } double jMWqbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector pe, CLHEP::HepLorentzVector pnu,CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in) // Calculates the square of the current contractions for qg->qenug scattering // p1: quark // p2: gluon { CLHEP::HepLorentzVector q1=p1in-p1out-pe-pnu; CLHEP::HepLorentzVector q2=-(p2in-p2out); current mj1m,mj2p,mj2m; jWbar(p1out,false,pe,false,pnu,false,p1in,false,mj1m); j(p2out,true,p2in,true,mj2p); j(p2out,false,p2in,false,mj2m); // mj1m.mj2p COM Mmp=cdot(mj1m,mj2p); // mj1m.mj2m COM Mmm=cdot(mj1m,mj2m); const double K = K_g(p2out, p2in); // sum of spinor strings ||^2 double a2Mmp=abs2(Mmp); double a2Mmm=abs2(Mmm); double sst = K/C_A*(a2Mmp+a2Mmm); // double sstsave=sst; // // Leave division by colour and Helicity avg until Tree files // Leave multi. of couplings to later // Multiply by Cf*Ca=4 return 4.*sst/(q1.m2()*q2.m2()); } double jM2qQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in) { // std::cerr<<"Current: "<<p1out<<" "<<p1in<<" "<<p2out<<" "<<p2in<<std::endl; CLHEP::HepLorentzVector q1=p1in-p1out; CLHEP::HepLorentzVector q2=-(p2in-p2out); // std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl; current mj1m,mj1p,mj2m,mj2p; j(p1out,true,p1in,true,mj1p); j(p1out,false,p1in,false,mj1m); j(p2out,true,p2in,true,mj2p); j(p2out,false,p2in,false,mj2m); COM Mmp=cdot(mj1m,mj2p); COM Mmm=cdot(mj1m,mj2m); COM Mpp=cdot(mj1p,mj2p); COM Mpm=cdot(mj1p,mj2m); double sst=abs2(Mmm)+abs2(Mmp)+abs2(Mpp)+abs2(Mpm); // Multiply by Cf^2 return RHEJ::C_F*RHEJ::C_F*(sst)/(q1.m2()*q2.m2()); } double jM2qQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in) { CLHEP::HepLorentzVector q1=p1in-p1out; CLHEP::HepLorentzVector q2=-(p2in-p2out); current mj1m,mj1p,mj2m,mj2p; j(p1out,true,p1in,true,mj1p); j(p1out,false,p1in,false,mj1m); jio(p2in,true,p2out,true,mj2p); jio(p2in,false,p2out,false,mj2m); COM Mmp=cdot(mj1m,mj2p); COM Mmm=cdot(mj1m,mj2m); COM Mpp=cdot(mj1p,mj2p); COM Mpm=cdot(mj1p,mj2m); double sumsq=abs2(Mmm)+abs2(Mmp)+abs2(Mpp)+abs2(Mpm); // Multiply by Cf^2 return (4./3.)*(4./3.)*(sumsq)/(q1.m2()*q2.m2()); } double jM2qbarQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in) { CLHEP::HepLorentzVector q1=p1in-p1out; CLHEP::HepLorentzVector q2=-(p2in-p2out); current mj1m,mj1p,mj2m,mj2p; jio(p1in,true,p1out,true,mj1p); jio(p1in,false,p1out,false,mj1m); jio(p2in,true,p2out,true,mj2p); jio(p2in,false,p2out,false,mj2m); COM Mmp=cdot(mj1m,mj2p); COM Mmm=cdot(mj1m,mj2m); COM Mpp=cdot(mj1p,mj2p); COM Mpm=cdot(mj1p,mj2m); double sumsq=abs2(Mmm)+abs2(Mmp)+abs2(Mpp)+abs2(Mpm); // Multiply by Cf^2 return (4./3.)*(4./3.)*(sumsq)/(q1.m2()*q2.m2()); } double jM2qg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in) // Calculates the square of the current contractions for qg scattering // p1: quark // p2: gluon { CLHEP::HepLorentzVector q1=p1in-p1out; CLHEP::HepLorentzVector q2=-(p2in-p2out); current mj1m,mj1p,mj2m,mj2p; j(p1out,true,p1in,true,mj1p); j(p1out,false,p1in,false,mj1m); j(p2out,true,p2in,true,mj2p); j(p2out,false,p2in,false,mj2m); COM Mmp=cdot(mj1m,mj2p); COM Mmm=cdot(mj1m,mj2m); COM Mpp=cdot(mj1p,mj2p); COM Mpm=cdot(mj1p,mj2m); const double K = K_g(p2out, p2in); // sum of spinor strings ||^2 double a2Mmp=abs2(Mmp); double a2Mmm=abs2(Mmm); double a2Mpp=abs2(Mpp); double a2Mpm=abs2(Mpm); double sst = K/C_A*(a2Mpp+a2Mpm+a2Mmp+a2Mmm); // double sstsave=sst; // std::cout <<"ratio: "<<sst/sstsave<<std::endl; // Cf*Ca=4 return 4.*sst/(q1.m2()*q2.m2()); } double jM2qbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in) // Calculates the square of the current contractions for qg scattering // p1: quark // p2: gluon { CLHEP::HepLorentzVector q1=p1in-p1out; CLHEP::HepLorentzVector q2=-(p2in-p2out); current mj1m,mj1p,mj2m,mj2p; jio(p1in,true,p1out,true,mj1p); jio(p1in,false,p1out,false,mj1m); j(p2out,true,p2in,true,mj2p); j(p2out,false,p2in,false,mj2m); COM Mmp=cdot(mj1m,mj2p); COM Mmm=cdot(mj1m,mj2m); COM Mpp=cdot(mj1p,mj2p); COM Mpm=cdot(mj1p,mj2m); const double K = K_g(p2out, p2in); // sum of spinor strings ||^2 double a2Mmp=abs2(Mmp); double a2Mmm=abs2(Mmm); double a2Mpp=abs2(Mpp); double a2Mpm=abs2(Mpm); double sst = K/C_A*(a2Mpp+a2Mpm+a2Mmp+a2Mmm); // double sstsave=sst; // std::cout <<"ratio: "<<sst/sstsave<<std::endl; // Cf*Ca=4 return 4.*sst/(q1.m2()*q2.m2()); } double jM2gg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in) // Calculates the square of the current contractions for gg scattering // p1: gluon // p2: gluon { CLHEP::HepLorentzVector q1=p1in-p1out; CLHEP::HepLorentzVector q2=-(p2in-p2out); current mj1m,mj1p,mj2m,mj2p; j(p1out,true,p1in,true,mj1p); j(p1out,false,p1in,false,mj1m); j(p2out,true,p2in,true,mj2p); j(p2out,false,p2in,false,mj2m); COM Mmp=cdot(mj1m,mj2p); COM Mmm=cdot(mj1m,mj2m); COM Mpp=cdot(mj1p,mj2p); COM Mpm=cdot(mj1p,mj2m); const double K_g1 = K_g(p1out, p1in); const double K_g2 = K_g(p2out, p2in); // sum of spinor strings ||^2 double a2Mmp=abs2(Mmp); double a2Mmm=abs2(Mmm); double a2Mpp=abs2(Mpp); double a2Mpm=abs2(Mpm); double sst = K_g1/C_A*K_g2/C_A*(a2Mpp+a2Mpm+a2Mmp+a2Mmm); // double sstsave=sst; // std::cout <<"ratio: "<<sst/sstsave<<std::endl; // Ca*Ca=9 return 9.*sst/(q1.m2()*q2.m2()); } namespace { /// @TODO what was this intended to do? // double MH2helper(current C1, current C2, current q1, current q2) // { // COM M; // COM temp1,temp2; // // First the C1.q2 * C2.q1 - part // temp1=cdot(C1,q2); // temp2=cdot(C2,q1); // M=temp1*temp2; // // Then the C1.C2 * q1.q2 // temp1=cdot(C1,C2); // temp2=cdot(q1,q2); // M-=temp1*temp2; // return (M*conj(M)).real(); // } /** * @brief Higgs vertex contracted with current @param C1 and @param C2 */ COM cHdot(const current & C1, const current & C2, const current & q1, const current & q2, double mt, bool incBot, double mb) { if (mt == infinity) { return (cdot(C1,C2)*cdot(q1,q2)-cdot(C1,q2)*cdot(C2,q1))/(6*M_PI*v); } else { CLHEP::HepLorentzVector vq1,vq2; vq1.set(q1[1].real(),q1[2].real(),q1[3].real(),q1[0].real()); vq2.set(q2[1].real(),q2[2].real(),q2[3].real(),q2[0].real()); // first minus sign obtained because of q1-difference to VDD // std::cout<<"A1 : " << A1(-vq1,vq2)<<std::endl; // std::cout<<"A2 : " << A2(-vq1,vq2)<<std::endl; if(!(incBot)) // Factor is because 4 mt^2 g^2/v A1 -> 16 pi mt^2/v alphas, // and we divide by a factor 4 at the amp sqaured level later // which I absorb here (i.e. I divide by 2) /// @TODO move factor 1/2 from S to |ME|^2 => consistent with general notation return 8.*M_PI*mt*mt/v*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mt)-cdot(C1,C2)*A2(-vq1,vq2,mt)); else return 8.*M_PI*mt*mt/v*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mt)-cdot(C1,C2)*A2(-vq1,vq2,mt)) + 8.*M_PI*mb*mb/v*(-cdot(C1,q2)*cdot(C2,q1)*A1(-vq1,vq2,mb)-cdot(C1,C2)*A2(-vq1,vq2,mb)); } } } // namespace anonymous double MH2qQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) { // CLHEP::HepLorentzVector q1=p1in-p1out; // CLHEP::HepLorentzVector q2=-(p2in-p2out); current j1p,j1m,j2p,j2m, q1v, q2v; j (p1out,true,p1in,true,j1p); j (p1out,false,p1in,false,j1m); j (p2out,true,p2in,true,j2p); j (p2out,false,p2in,false,j2m); to_current(q1, q1v); to_current(q2, q2v); COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb); COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb); COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb); COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb); double sst=abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm); // return (4./3.)*(4./3.)*sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); return sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); } double MH2qQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) { // CLHEP::HepLorentzVector q1=p1in-p1out; // CLHEP::HepLorentzVector q2=-(p2in-p2out); current j1p,j1m,j2p,j2m,q1v,q2v; j (p1out,true,p1in,true,j1p); j (p1out,false,p1in,false,j1m); jio (p2in,true,p2out,true,j2p); jio (p2in,false,p2out,false,j2m); to_current(q1, q1v); to_current(q2, q2v); COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb); COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb); COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb); COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb); double sst=abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm); // return (4./3.)*(4./3.)*sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); return sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); } double MH2qbarQ (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) { // CLHEP::HepLorentzVector q1=p1in-p1out; // CLHEP::HepLorentzVector q2=-(p2in-p2out); current j1p,j1m,j2p,j2m,q1v,q2v; jio (p1in,true,p1out,true,j1p); jio (p1in,false,p1out,false,j1m); j (p2out,true,p2in,true,j2p); j (p2out,false,p2in,false,j2m); to_current(q1, q1v); to_current(q2, q2v); COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb); COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb); COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb); COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb); double sst=abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm); // return (4./3.)*(4./3.)*sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); return sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); } double MH2qbarQbar (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) { // CLHEP::HepLorentzVector q1=p1in-p1out; // CLHEP::HepLorentzVector q2=-(p2in-p2out); current j1p,j1m,j2p,j2m,q1v,q2v; jio (p1in,true,p1out,true,j1p); jio (p1in,false,p1out,false,j1m); jio (p2in,true,p2out,true,j2p); jio (p2in,false,p2out,false,j2m); to_current(q1, q1v); to_current(q2, q2v); COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb); COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb); COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb); COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb); double sst=abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm); // return (4./3.)*(4./3.)*sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); return sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); } double MH2qg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) // q~p1 g~p2 (i.e. ALWAYS p1 for quark, p2 for gluon) // should be called with q1 meant to be contracted with p2 in first part of vertex // (i.e. if g is backward, q1 is forward) { current j1p,j1m,j2p,j2m,q1v,q2v; j (p1out,true,p1in,true,j1p); j (p1out,false,p1in,false,j1m); j (p2out,true,p2in,true,j2p); j (p2out,false,p2in,false,j2m); to_current(q1, q1v); to_current(q2, q2v); // First, calculate the non-flipping amplitudes: COM Mpp=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb); COM Mpm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb); COM Mmp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb); COM Mmm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb); //cout << "Bits in MH2qg: " << Mpp << " " << Mpm << " " << Mmp << " " << Mmm << endl; const double K = K_g(p2out, p2in); double sst=K/C_A*(abs2(Mmp)+abs2(Mmm)+abs2(Mpp)+abs2(Mpm)); // Cf*Ca=4 // return 4.*sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); return sst/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); } double MH2qbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) // qbar~p1 g~p2 (i.e. ALWAYS p1 for anti-quark, p2 for gluon) // should be called with q1 meant to be contracted with p2 in first part of vertex // (i.e. if g is backward, q1 is forward) { current j1p,j1m,j2p,j2m,q1v,q2v; jio (p1in,true,p1out,true,j1p); jio (p1in,false,p1out,false,j1m); j (p2out,true,p2in,true,j2p); j (p2out,false,p2in,false,j2m); to_current(q1, q1v); to_current(q2, q2v); // First, calculate the non-flipping amplitudes: COM amp,amm,apm,app; app=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb); apm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb); amp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb); amm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb); double MH2sum = abs2(app)+abs2(amm)+abs2(apm)+abs2(amp); const double K = K_g(p2out, p2in); MH2sum*=K/C_A; // Cf*Ca=4 // return 4.*MH2sum/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); return MH2sum/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); } double MH2gg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) // g~p1 g~p2 // should be called with q1 meant to be contracted with p2 in first part of vertex // (i.e. if g is backward, q1 is forward) { current j1p,j1m,j2p,j2m,q1v,q2v; j (p1out,true,p1in,true,j1p); j (p1out,false,p1in,false,j1m); j (p2out,true,p2in,true,j2p); j (p2out,false,p2in,false,j2m); to_current(q1, q1v); to_current(q2, q2v); // First, calculate the non-flipping amplitudes: COM amp,amm,apm,app; app=cHdot(j1p,j2p,q1v,q2v,mt, incBot, mb); apm=cHdot(j1p,j2m,q1v,q2v,mt, incBot, mb); amp=cHdot(j1m,j2p,q1v,q2v,mt, incBot, mb); amm=cHdot(j1m,j2m,q1v,q2v,mt, incBot, mb); double MH2sum = abs2(app)+abs2(amm)+abs2(apm)+abs2(amp); const double K_g1 = K_g(p1out, p1in); const double K_g2 = K_g(p2out, p2in); MH2sum*=K_g1/C_A*K_g2/C_A; // Ca*Ca=9 // return 9.*MH2sum/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); return MH2sum/((p1in-p1out).m2()*(p2in-p2out).m2()*q1.m2()*q2.m2()); } // // Z's stuff // void jZ(HLV pin, HLV pout, HLV pem, HLV pep, bool HelPartons, bool HelLeptons, current cur) { // // Init current to zero // cur[0] = 0.0; // cur[1] = 0.0; // cur[2] = 0.0; // cur[3] = 0.0; // // Temporary variables // COM temp; // current Term_1, Term_2, Term_3, Term_4, J_temp, TempCur1, TempCur2; // // Momentum of virtual gluons aroun weak boson emission site // HLV qa = pout + pep + pem; // HLV qb = pin - pep - pem; // double ta = qa.m2(); // double tb = qb.m2(); // // Out-Out currents: // current Em_Ep, Out_Em, Out_Ep; // // Other currents: // current Out_In, Em_In, Ep_In; // joi(pout, HelPartons, pin, HelPartons, Out_In); // joi(pem, HelLeptons, pin, HelPartons, Em_In); // joi(pep, HelLeptons, pin, HelPartons, Ep_In); // joo(pem, HelLeptons, pep, HelLeptons, Em_Ep); // joo(pout, HelPartons, pem, HelLeptons, Out_Em); // joo(pout, HelPartons, pep, HelLeptons, Out_Ep); // if (HelLeptons == HelPartons) { // temp = 2.0 * cdot(pout, Em_Ep); // cmult(temp / ta, Out_In, Term_1); // temp = cdot(Out_Em, Em_Ep); // cmult(temp / ta , Em_In, Term_2); // temp = 2.0 * cdot(pin, Em_Ep); // cmult(temp / tb, Out_In, Term_3); // temp = -cdot(Ep_In, Em_Ep); // cmult(temp / tb, Out_Ep, Term_4); // cadd(Term_1, Term_2, Term_3, Term_4, J_temp); // cur[0] = J_temp[0]; // cur[1] = J_temp[1]; // cur[2] = J_temp[2]; // cur[3] = J_temp[3]; // } // else { // if (HelPartons == true) { // temp = 2.0 * cdot(pout, Em_Ep); // cmult(temp / ta, Out_In, Term_1); // joo(pout, true, pep, true, TempCur1); // joi(pep, true, pin, true, TempCur2); // temp = cdot(TempCur1, Em_Ep); // cmult(temp / ta , TempCur2, Term_2); // temp = 2.0 * cdot(pin, Em_Ep); // cmult(temp / tb, Out_In, Term_3); // joo(pout, true, pem, true, TempCur1); // joi(pem, true, pin, true, TempCur2); // temp = -cdot(TempCur2, Em_Ep); // cmult(temp / tb, TempCur1, Term_4); // cadd(Term_1, Term_2, Term_3, Term_4, J_temp); // cur[0] = J_temp[0]; // cur[1] = J_temp[1]; // cur[2] = J_temp[2]; // cur[3] = J_temp[3]; // } // else { // temp = 2.0 * cdot(pout, Em_Ep); // cmult(temp / ta, Out_In, Term_1); // joo(pout, false, pep, false, TempCur1); // joi(pep, false, pin, false, TempCur2); // temp = cdot(TempCur1, Em_Ep); // cmult(temp / ta, TempCur2, Term_2); // temp = 2.0 * cdot(pin, Em_Ep); // cmult(temp / tb, Out_In, Term_3); // joo(pout, false, pem, false, TempCur1); // joi(pem, false, pin, false, TempCur2); // temp = -cdot(TempCur2, Em_Ep); // cmult(temp / tb, TempCur1, Term_4); // cadd(Term_1, Term_2, Term_3, Term_4, J_temp); // cur[0] = J_temp[0]; // cur[1] = J_temp[1]; // cur[2] = J_temp[2]; // cur[3] = J_temp[3]; // } // } // } // void jZbar(HLV pin, HLV pout, HLV pem, HLV pep, bool HelPartons, bool HelLeptons, current cur) { // // Init current to zero // cur[0] = 0.0; // cur[1] = 0.0; // cur[2] = 0.0; // cur[3] = 0.0; // // Temporary variables // COM temp; // current Term_1, Term_2, Term_3, Term_4, J_temp, TempCur1, TempCur2; // // Transfered 4-momenta // HLV qa = pout + pep + pem; // HLV qb = pin - pep - pem; // // The square of the transfered 4-momenta // double ta = qa.m2(); // double tb = qb.m2(); // // Out-Out currents: // current Em_Ep, Em_Out, Ep_Out; // // In-Out currents: // current In_Out, In_Em, In_Ep; // // Safe to use the currents since helicity structure is ok // if (HelPartons == HelLeptons) { // jio(pin, HelPartons, pout, HelPartons, In_Out); // joo(pem, HelLeptons, pep, HelLeptons, Em_Ep); // jio(pin, HelPartons, pem, HelLeptons, In_Em); // jio(pin, HelPartons, pep, HelLeptons, In_Ep); // joo(pem, HelLeptons, pout, HelPartons, Em_Out); // joo(pep, HelLeptons, pout, HelPartons, Ep_Out); // } // else { // jio(pin, HelPartons, pout, HelPartons, In_Out); // joo(pem, HelLeptons, pep, HelLeptons, Em_Ep); // In_Em[0] = 0.0; // In_Em[1] = 0.0; // In_Em[2] = 0.0; // In_Em[3] = 0.0; // In_Ep[0] = 0.0; // In_Ep[1] = 0.0; // In_Ep[2] = 0.0; // In_Ep[3] = 0.0; // Em_Out[0] = 0.0; // Em_Out[1] = 0.0; // Em_Out[2] = 0.0; // Em_Out[3] = 0.0; // Ep_Out[0] = 0.0; // Ep_Out[1] = 0.0; // Ep_Out[2] = 0.0; // Ep_Out[3] = 0.0; // } // if (HelLeptons == HelPartons) { // temp = 2.0 * cdot(pout, Em_Ep); // cmult(temp / ta, In_Out, Term_1); // temp = cdot(Ep_Out, Em_Ep); // cmult(temp / ta, In_Ep, Term_2); // temp = 2.0 * cdot(pin, Em_Ep); // cmult(temp / tb, In_Out, Term_3); // temp = - cdot(In_Em, Em_Ep); // cmult(temp / tb, Em_Out, Term_4); // cadd(Term_1, Term_2, Term_3, Term_4, J_temp); // cur[0] = J_temp[0]; // cur[1] = J_temp[1]; // cur[2] = J_temp[2]; // cur[3] = J_temp[3]; // } // else { // if (HelPartons == true) { // temp = 2.0 * cdot(pout, Em_Ep); // cmult(temp / ta, In_Out, Term_1); // joo(pem, true, pout, true, TempCur1); // jio(pin, true, pem, true, TempCur2); // temp = cdot(TempCur1, Em_Ep); // cmult(temp / ta , TempCur2, Term_2); // temp = 2.0 * cdot(pin, Em_Ep); // cmult(temp / tb, In_Out, Term_3); // joo(pep, true, pout, true, TempCur1); // jio(pin, true, pep, true, TempCur2); // temp = - cdot(TempCur2, Em_Ep); // cmult(temp / tb, TempCur1, Term_4); // cadd(Term_1, Term_2, Term_3, Term_4, J_temp); // cur[0] = J_temp[0]; // cur[1] = J_temp[1]; // cur[2] = J_temp[2]; // cur[3] = J_temp[3]; // } // else { // temp = 2.0 * cdot(pout, Em_Ep); // cmult(temp / ta, In_Out, Term_1); // joo(pem, false, pout, false, TempCur1); // jio(pin, false, pem, false, TempCur2); // temp = cdot(TempCur1, Em_Ep); // cmult(temp / ta , TempCur2, Term_2); // temp = 2.0 * cdot(pin, Em_Ep); // cmult(temp / tb, In_Out, Term_3); // joo(pep, false, pout, false, TempCur1); // jio(pin, false, pep, false, TempCur2); // temp = - cdot(TempCur2, Em_Ep); // cmult(temp / tb, TempCur1, Term_4); // cadd(Term_1, Term_2, Term_3, Term_4, J_temp); // cur[0] = J_temp[0]; // cur[1] = J_temp[1]; // cur[2] = J_temp[2]; // cur[3] = J_temp[3]; // } // } // } // // Progagators // COM PZ(double s) { // double MZ, GammaZ; // MZ = 9.118800e+01; // Mass of the mediating gauge boson // GammaZ = 2.441404e+00; // Z peak width // // Return Z Prop value // return 1.0 / (s - MZ * MZ + COM(0.0, 1.0) * GammaZ * MZ); // } // COM PG(double s) { // return 1.0 / s; // } // // Non-gluonic with pa emitting // std::vector <double> jMZqQ (HLV pa, HLV pb, HLV p1, HLV p2, HLV pep, HLV pem, std::vector <double> VProducts, std::vector < std::vector <double> > Virtuals, int aptype, int bptype, bool UseVirtuals, bool BottomLineEmit) { // std::vector <double> ScaledWeights; // double Sum; // // Propagator factors // COM PZs = PZ((pep + pem).m2()); // COM PGs = PG((pep + pem).m2()); // // Emitting current initialisation // current j1pptop, j1pmtop; // Emission from top line // current j1ppbot, j1pmbot; // Emission from bottom line // // Non-emitting current initialisation // current j2ptop, j2mtop; // Emission from top line // current j2pbot, j2mbot; // Emission from bottom line // // Currents for top emission // // Upper current calculations // // if a is a quark // if (aptype > 0) { // jZ(pa, p1, pem, pep, true, true, j1pptop); // jZ(pa, p1, pem, pep, true, false, j1pmtop); // } // // if a is an antiquark // else { // jZbar(pa, p1, pem, pep, true, true, j1pptop); // jZbar(pa, p1, pem, pep, true, false, j1pmtop); // } // // Lower current calculations // // if b is a quark // if (bptype > 0) { // joi(p2, true, pb, true, j2ptop); // joi(p2, false, pb, false, j2mtop); // } // // if b is an antiquark // else { // jio(pb, true, p2, true, j2ptop); // jio(pb, false, p2, false, j2mtop); // } // // Currents for bottom emission // // Lower current calculations // if (bptype > 0) { // jZ(pb, p2, pem, pep, true, true, j1ppbot); // jZ(pb, p2, pem, pep, true, false, j1pmbot); // } // else { // jZbar(pb, p2, pem, pep, true, true, j1ppbot); // jZbar(pb, p2, pem, pep, true, false, j1pmbot); // } // // Upper current calculations // if (aptype > 0) { // joi(p1, true, pa, true, j2pbot); // joi(p1, false, pa, false, j2mbot); // } // else { // jio(pa, true, p1, true, j2pbot); // jio(pa, false, p1, false, j2mbot); // } // COM Coeff[2][8]; // if (!Interference) { // double ZCharge_a_P = Zq(aptype, true); // double ZCharge_a_M = Zq(aptype, false); // double ZCharge_b_P = Zq(bptype, true); // double ZCharge_b_M = Zq(bptype, false); // if (BottomLineEmit) { // // Emission from top-line quark (pa/p1 line) // Coeff[0][0] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pptop, j2ptop); // Coeff[0][1] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pptop, j2mtop); // Coeff[0][2] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pmtop, j2ptop); // Coeff[0][3] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pmtop, j2mtop); // Coeff[0][4] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2ptop)); // Coeff[0][5] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2mtop)); // Coeff[0][6] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2ptop)); // Coeff[0][7] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2mtop)); // } // else { // // Emission from bottom-line quark (pb/p2 line) // Coeff[1][0] = (ZCharge_b_P * Zep * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1ppbot, j2pbot); // Coeff[1][7] = (ZCharge_b_P * Zep * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1ppbot, j2mbot); // Coeff[1][2] = (ZCharge_b_P * Zem * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1pmbot, j2pbot); // Coeff[1][5] = (ZCharge_b_P * Zem * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1pmbot, j2mbot); // Coeff[1][4] = (ZCharge_b_M * Zem * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1ppbot, j2pbot)); // Coeff[1][3] = (ZCharge_b_M * Zem * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1ppbot, j2mbot)); // Coeff[1][6] = (ZCharge_b_M * Zep * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1pmbot, j2pbot)); // Coeff[1][1] = (ZCharge_b_M * Zep * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1pmbot, j2mbot)); // } // } // // Else calculate all the possiblities // else { // double ZCharge_a_P = Zq(aptype, true); // double ZCharge_a_M = Zq(aptype, false); // double ZCharge_b_P = Zq(bptype, true); // double ZCharge_b_M = Zq(bptype, false); // // Emission from top-line quark (pa/p1 line) // Coeff[0][0] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pptop, j2ptop); // Coeff[0][1] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pptop, j2mtop); // Coeff[0][2] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pmtop, j2ptop); // Coeff[0][3] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pmtop, j2mtop); // Coeff[0][4] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2ptop)); // Coeff[0][5] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2mtop)); // Coeff[0][6] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2ptop)); // Coeff[0][7] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2mtop)); // // Emission from bottom-line quark (pb/p2 line) // Coeff[1][0] = (ZCharge_b_P * Zep * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1ppbot, j2pbot); // Coeff[1][7] = (ZCharge_b_P * Zep * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1ppbot, j2mbot); // Coeff[1][2] = (ZCharge_b_P * Zem * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1pmbot, j2pbot); // Coeff[1][5] = (ZCharge_b_P * Zem * PZs * RWeak + Gq(bptype) * PGs) * cdot(j1pmbot, j2mbot); // Coeff[1][4] = (ZCharge_b_M * Zem * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1ppbot, j2pbot)); // Coeff[1][3] = (ZCharge_b_M * Zem * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1ppbot, j2mbot)); // Coeff[1][6] = (ZCharge_b_M * Zep * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1pmbot, j2pbot)); // Coeff[1][1] = (ZCharge_b_M * Zep * PZs * RWeak + Gq(bptype) * PGs) * conj(cdot(j1pmbot, j2mbot)); // } // // Find the numbers of scales // int ScaleCount; // #if calcscaleunc // ScaleCount = 20; // #else // ScaleCount = 1; // #endif // // For each scale... // for (int j = 0; j < ScaleCount; j++) { // Sum = 0.0; // // If we want to compare back to the W's code only emit from one quark and only couple to left handed particles // // virtuals arent here since they are calculated and included in weight() call. // if (!Interference) { // if (BottomLineEmit) for (int i = 0; i < 8; i++) Sum += abs2(Coeff[1][i]) * VProducts.at(1); // else for (int i = 0; i < 8; i++) Sum += abs2(Coeff[0][i]) * VProducts.at(0); // } // // Else work out the full interference // else { // // For the full calculation... // if (UseVirtuals) { // for (int i = 0; i < 8; i++) { // Sum += abs2(Coeff[0][i]) * VProducts.at(0) * Virtuals.at(j).at(0) // + abs2(Coeff[1][i]) * VProducts.at(1) * Virtuals.at(j).at(1) // + 2.0 * real(Coeff[0][i] * conj(Coeff[1][i])) * VProducts.at(2) * Virtuals.at(j).at(2); // } // } // // For the tree level calculation... // else { // for (int i = 0; i < 8; i++) { // Sum += abs2(Coeff[0][i]) * VProducts.at(0) // + abs2(Coeff[1][i]) * VProducts.at(1) // + 2.0 * real(Coeff[0][i] * conj(Coeff[1][i])) * VProducts.at(2); // } // } // } // // Add this to the vector to be returned with the other factors of C_A and the helicity sum/average factors. // ScaledWeights.push_back(Sum / 18.0); // } // // Return all the scale values // return ScaledWeights; // } // // Semi-gluonic with pa emitting // std::vector <double> jMZqg (HLV pa, HLV pb, HLV p1, HLV p2, HLV pep, HLV pem, std::vector <double> VProducts, std::vector < std::vector <double> > Virtuals, int aptype, int bptype, bool UseVirtuals, bool BottomLineEmit) { // COM Coeff[8]; // double Sum; // std::vector <double> ScaledWeights; // COM PZs = PZ((pep + pem).m2()); // COM PGs = PG((pep + pem).m2()); // // Emitting current initialisation - Emission from top line // current j1pptop, j1pmtop; // // Non-emitting current initialisation - Emission from top line // current j2ptop, j2mtop; // // Currents for top emission // // Upper current calculations // if (aptype > 0) { // jZ (pa, p1, pem, pep, true, true, j1pptop); // jZ (pa, p1, pem, pep, true, false, j1pmtop); // } // else { // jZbar(pa, p1, pem, pep, true, true, j1pptop); // jZbar(pa, p1, pem, pep, true, false, j1pmtop); // } // // Lower current calculations // joi(p2, true, pb, true, j2ptop); // joi(p2, false, pb, false, j2mtop); // // Calculate all the possiblities // double ZCharge_a_P = Zq(aptype, true); // double ZCharge_a_M = Zq(aptype, false); // // Emission from top-line quark (pa/p1 line) // Coeff[0] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pptop, j2ptop); // Coeff[1] = (ZCharge_a_P * Zep * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pptop, j2mtop); // Coeff[2] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pmtop, j2ptop); // Coeff[3] = (ZCharge_a_P * Zem * PZs * RWeak + Gq(aptype) * PGs) * cdot(j1pmtop, j2mtop); // Coeff[4] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2ptop)); // Coeff[5] = (ZCharge_a_M * Zem * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pptop, j2mtop)); // Coeff[6] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2ptop)); // Coeff[7] = (ZCharge_a_M * Zep * PZs * RWeak + Gq(aptype) * PGs) * conj(cdot(j1pmtop, j2mtop)); // // Calculate gluon colour accelerated factor // double CAMFactor, z; // // If b is a forward moving gluon define z (C.F. multiple jets papers) // if (pb.pz() > 0) z = p2.plus() / pb.plus(); // else z = p2.minus() / pb.minus(); // CAMFactor = (1.0 - 1.0 / 9.0) / 2.0 * (z + 1.0 / z) + 1.0 / 9.0; // // Find the numbers of scales // int ScaleCount; // #if calcscaleunc // ScaleCount = 20; // #else // ScaleCount = 1; // #endif // // For each scale... // for (int j = 0; j < ScaleCount; j++) { // Sum = 0.0; // // If we dont want the interference // if (!Interference) for (int i = 0; i < 8; i++) Sum += abs2(Coeff[i]) * VProducts.at(0); // // Else work out the full interference // else { // if (UseVirtuals) { // for (int i = 0; i < 8; i++) Sum += abs2(Coeff[i]) * VProducts.at(0) * Virtuals.at(j).at(0); // } // else { // for (int i = 0; i < 8; i++) Sum += abs2(Coeff[i]) * VProducts.at(0); // } // } // // Add this to the vector to be returned with the other factors of C_A, the colour accelerated factor and the helicity sum/average factors.: (4/3)*3/32 // ScaledWeights.push_back(CAMFactor * Sum / 8.0); // } // return ScaledWeights; // } // // Electroweak Charge Functions // double Zq (int PID, bool Helcitiy) { // double temp; // // Positive Spin // if (Helcitiy == true) { // if (PID == 1 || PID == 3 || PID == 5) temp = (+ 1.0 * stw2 / 3.0) / ctw; // if (PID == 2 || PID == 4) temp = (- 2.0 * stw2 / 3.0) / ctw; // if (PID == -1 || PID == -3 || PID == -5) temp = (- 1.0 * stw2 / 3.0) / ctw; // if (PID == -2 || PID == -4) temp = (+ 2.0 * stw2 / 3.0) / ctw; // // If electron or positron // if (PID == 7 || PID == -7) temp = Zep; // } // // Negative Spin // else { // if (PID == 1 || PID == 3 || PID == 5) temp = (-0.5 + 1.0 * stw2 / 3.0) / ctw; // if (PID == 2 || PID == 4) temp = ( 0.5 - 2.0 * stw2 / 3.0) / ctw; // if (PID == -1 || PID == -3 || PID == -5) temp = ( 0.5 - 1.0 * stw2 / 3.0) / ctw; // if (PID == -2 || PID == -4) temp = (-0.5 + 2.0 * stw2 / 3.0) / ctw; // // If electron or positron // if (PID == 7 || PID == -7) temp = Zem; // } // return temp; // } // double Gq (int PID) { // if (!VirtualPhoton) return 0.0; // if (PID == -1) return 1.0 * ee / 3.0; // if (PID == -2) return -2.0 * ee / 3.0; // if (PID == -3) return 1.0 * ee / 3.0; // if (PID == -4) return -2.0 * ee / 3.0; // if (PID == -5) return 1.0 * ee / 3.0; // if (PID == 1) return -1.0 * ee / 3.0; // if (PID == 2) return 2.0 * ee / 3.0; // if (PID == 3) return -1.0 * ee / 3.0; // if (PID == 4) return 2.0 * ee / 3.0; // if (PID == 5) return -1.0 * ee / 3.0; // std::cout << "ERROR! No Electroweak Charge Found at line " << __LINE__ << "..." << std::endl; // return 0.0; // } namespace { //@{ /// @brief Higgs vertex contracted with one current CCurrent jH (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) { CCurrent j2 = j(pout,helout,pin,helin); CCurrent jq2(q2.e(),q2.px(),q2.py(),q2.pz()); if(mt == infinity) return ((q1.dot(q2))*j2 - j2.dot(q1)*jq2)/(3*M_PI*v); else { if(incBot) return (-16.*M_PI*mb*mb/v*j2.dot(q1)*jq2*A1(-q1,q2,mb)-16.*M_PI*mb*mb/v*j2*A2(-q1,q2,mb)) + (-16.*M_PI*mt*mt/v*j2.dot(q1)*jq2*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j2*A2(-q1,q2,mt)); else return (-16.*M_PI*mt*mt/v*j2.dot(q1)*jq2*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j2*A2(-q1,q2,mt)); } } CCurrent jioH (CLHEP::HepLorentzVector pin, bool helin, CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) { CCurrent j2 = jio(pin,helin,pout,helout); CCurrent jq2(q2.e(),q2.px(),q2.py(),q2.pz()); if(mt == infinity) return ((q1.dot(q2))*j2 - j2.dot(q1)*jq2)/(3*M_PI*v); else { if(incBot) return (-16.*M_PI*mb*mb/v*j2.dot(q1)*jq2*A1(-q1,q2,mb)-16.*M_PI*mb*mb/v*j2*A2(-q1,q2,mb)) + (-16.*M_PI*mt*mt/v*j2.dot(q1)*jq2*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j2*A2(-q1,q2,mt)); else return (-16.*M_PI*mt*mt/v*j2.dot(q1)*jq2*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j2*A2(-q1,q2,mt)); } } CCurrent jHtop (CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector pin, bool helin, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) { CCurrent j1 = j(pout,helout,pin,helin); CCurrent jq1(q1.e(),q1.px(),q1.py(),q1.pz()); if(mt == infinity) return ((q1.dot(q2))*j1 - j1.dot(q2)*jq1)/(3*M_PI*v); else { if(incBot) return (-16.*M_PI*mb*mb/v*j1.dot(q2)*jq1*A1(-q1,q2,mb)-16.*M_PI*mb*mb/v*j1*A2(-q1,q2,mb)) + (-16.*M_PI*mt*mt/v*j1.dot(q2)*jq1*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j1*A2(-q1,q2,mt)); else return (-16.*M_PI*mt*mt/v*j1.dot(q2)*jq1*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j1*A2(-q1,q2,mt)); } } CCurrent jioHtop (CLHEP::HepLorentzVector pin, bool helin, CLHEP::HepLorentzVector pout, bool helout, CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mt, bool incBot, double mb) { CCurrent j1 = jio(pin,helin,pout,helout); CCurrent jq1(q1.e(),q1.px(),q1.py(),q1.pz()); if(mt == infinity) return ((q1.dot(q2))*j1 - j1.dot(q2)*jq1)/(3*M_PI*v); else { if(incBot) return (-16.*M_PI*mb*mb/v*j1.dot(q2)*jq1*A1(-q1,q2,mb)-16.*M_PI*mb*mb/v*j1*A2(-q1,q2,mb)) + (-16.*M_PI*mt*mt/v*j1.dot(q2)*jq1*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j1*A2(-q1,q2,mt)); else return (-16.*M_PI*mt*mt/v*j1.dot(q2)*jq1*A1(-q1,q2,mt)-16.*M_PI*mt*mt/v*j1*A2(-q1,q2,mt)); } } //@} } // namespace anonymous double jM2unogqHQ (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { // This construction is taking rapidity order: pg > p1out >> p2out // std::cerr<<"This Uno Current: "<<p1out<<" "<<p1in<<" "<<p2out<<" "<<p2in<<" "<<pg<<std::endl; CLHEP::HepLorentzVector q1=p1in-p1out; // Top End CLHEP::HepLorentzVector q2=-(p2in-p2out); // Bottom End CLHEP::HepLorentzVector qg=p1in-p1out-pg; // Extra bit post-gluon // std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl; CCurrent mj1m,mj1p,mj2m,mj2p,mjH2m,mjH2p; mj1p=j(p1out,true,p1in,true); mj1m=j(p1out,false,p1in,false); mjH2p=jH(p2out,true,p2in,true,qH1,qH2, mt, incBot, mb); mjH2m=jH(p2out,false,p2in,false,qH1,qH2, mt, incBot, mb); // Dot products of these which occur again and again COM MHmp=mj1m.dot(mjH2p); // And now for the Higgs ones COM MHmm=mj1m.dot(mjH2m); COM MHpp=mj1p.dot(mjH2p); COM MHpm=mj1p.dot(mjH2m); // std::cout<< p1out.rapidity() << " " << p2out.rapidity()<< " " << qH1 << " " << qH2 << "\n" <<MHmm << " " << MHmp << " " << MHpm << " " << MHpp << std::endl; // Currents with pg CCurrent jgam,jgap,j2gm,j2gp; j2gp=joo(p1out,true,pg,true); j2gm=joo(p1out,false,pg,false); jgap=j(pg,true,p1in,true); jgam=j(pg,false,p1in,false); CCurrent qsum(q1+qg); CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in); CCurrent p1o(p1out); CCurrent p1i(p1in); Lmm=(qsum*(MHmm) + (-2.*mjH2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHmm/2.))/q1.m2(); Lmp=(qsum*(MHmp) + (-2.*mjH2p.dot(pg))*mj1m+2.*mj1m.dot(pg)*mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHmp/2.))/q1.m2(); Lpm=(qsum*(MHpm) + (-2.*mjH2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHpm/2.))/q1.m2(); Lpp=(qsum*(MHpp) + (-2.*mjH2p.dot(pg))*mj1p+2.*mj1p.dot(pg)*mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHpp/2.))/q1.m2(); U1mm=(jgam.dot(mjH2m)*j2gm+2.*p1o*MHmm)/(p1out+pg).m2(); U1mp=(jgam.dot(mjH2p)*j2gm+2.*p1o*MHmp)/(p1out+pg).m2(); U1pm=(jgap.dot(mjH2m)*j2gp+2.*p1o*MHpm)/(p1out+pg).m2(); U1pp=(jgap.dot(mjH2p)*j2gp+2.*p1o*MHpp)/(p1out+pg).m2(); U2mm=((-1.)*j2gm.dot(mjH2m)*jgam+2.*p1i*MHmm)/(p1in-pg).m2(); U2mp=((-1.)*j2gm.dot(mjH2p)*jgam+2.*p1i*MHmp)/(p1in-pg).m2(); U2pm=((-1.)*j2gp.dot(mjH2m)*jgap+2.*p1i*MHpm)/(p1in-pg).m2(); U2pp=((-1.)*j2gp.dot(mjH2p)*jgap+2.*p1i*MHpp)/(p1in-pg).m2(); const double cf=RHEJ::C_F; double amm,amp,apm,app; amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm); amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp); apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm); app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp); double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2()); // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH1.m2()*qg.m2(); ampsq/=th; ampsq/=16.; ampsq*=RHEJ::C_F*RHEJ::C_F/RHEJ::C_A/RHEJ::C_A; // Factor of (Cf/Ca) for each quark to match MH2qQ. //Higgs coupling is included in Hjets.C return ampsq; } double jM2unogqbarHQ (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { // This construction is taking rapidity order: pg > p1out >> p2out // std::cerr<<"This Uno Current: "<<p1out<<" "<<p1in<<" "<<p2out<<" "<<p2in<<" "<<pg<<std::endl; CLHEP::HepLorentzVector q1=p1in-p1out; // Top End CLHEP::HepLorentzVector q2=-(p2in-p2out); // Bottom End CLHEP::HepLorentzVector qg=p1in-p1out-pg; // Extra bit post-gluon // std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl; CCurrent mj1m,mj1p,mj2m,mj2p,mjH2m,mjH2p; mj1p=jio(p1in,true,p1out,true); mj1m=jio(p1in,false,p1out,false); mjH2p=jH(p2out,true,p2in,true,qH1,qH2, mt, incBot, mb); mjH2m=jH(p2out,false,p2in,false,qH1,qH2, mt, incBot, mb); // Dot products of these which occur again and again COM MHmp=mj1m.dot(mjH2p); // And now for the Higgs ones COM MHmm=mj1m.dot(mjH2m); COM MHpp=mj1p.dot(mjH2p); COM MHpm=mj1p.dot(mjH2m); // Currents with pg CCurrent jgam,jgap,j2gm,j2gp; j2gp=joo(pg,true,p1out,true); j2gm=joo(pg,false,p1out,false); jgap=jio(p1in,true,pg,true); jgam=jio(p1in,false,pg,false); CCurrent qsum(q1+qg); CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in); CCurrent p1o(p1out); CCurrent p1i(p1in); Lmm=(qsum*(MHmm) + (-2.*mjH2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHmm/2.))/q1.m2(); Lmp=(qsum*(MHmp) + (-2.*mjH2p.dot(pg))*mj1m+2.*mj1m.dot(pg)*mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHmp/2.))/q1.m2(); Lpm=(qsum*(MHpm) + (-2.*mjH2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHpm/2.))/q1.m2(); Lpp=(qsum*(MHpp) + (-2.*mjH2p.dot(pg))*mj1p+2.*mj1p.dot(pg)*mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHpp/2.))/q1.m2(); U1mm=(jgam.dot(mjH2m)*j2gm+2.*p1o*MHmm)/(p1out+pg).m2(); U1mp=(jgam.dot(mjH2p)*j2gm+2.*p1o*MHmp)/(p1out+pg).m2(); U1pm=(jgap.dot(mjH2m)*j2gp+2.*p1o*MHpm)/(p1out+pg).m2(); U1pp=(jgap.dot(mjH2p)*j2gp+2.*p1o*MHpp)/(p1out+pg).m2(); U2mm=((-1.)*j2gm.dot(mjH2m)*jgam+2.*p1i*MHmm)/(p1in-pg).m2(); U2mp=((-1.)*j2gm.dot(mjH2p)*jgam+2.*p1i*MHmp)/(p1in-pg).m2(); U2pm=((-1.)*j2gp.dot(mjH2m)*jgap+2.*p1i*MHpm)/(p1in-pg).m2(); U2pp=((-1.)*j2gp.dot(mjH2p)*jgap+2.*p1i*MHpp)/(p1in-pg).m2(); const double cf=RHEJ::C_F; double amm,amp,apm,app; amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm); amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp); apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm); app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp); double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2()); // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH1.m2()*qg.m2(); ampsq/=th; ampsq/=16.; ampsq*=4.*4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ. //Higgs coupling is included in Hjets.C return ampsq; } double jM2unogqHQbar (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { // This construction is taking rapidity order: pg > p1out >> p2out // std::cerr<<"This Uno Current: "<<p1out<<" "<<p1in<<" "<<p2out<<" "<<p2in<<" "<<pg<<std::endl; CLHEP::HepLorentzVector q1=p1in-p1out; // Top End CLHEP::HepLorentzVector q2=-(p2in-p2out); // Bottom End CLHEP::HepLorentzVector qg=p1in-p1out-pg; // Extra bit post-gluon // std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl; CCurrent mj1m,mj1p,mj2m,mj2p,mjH2m,mjH2p; mj1p=j(p1out,true,p1in,true); mj1m=j(p1out,false,p1in,false); mjH2p=jioH(p2in,true,p2out,true,qH1,qH2, mt, incBot, mb); mjH2m=jioH(p2in,false,p2out,false,qH1,qH2, mt, incBot, mb); // Dot products of these which occur again and again COM MHmp=mj1m.dot(mjH2p); // And now for the Higgs ones COM MHmm=mj1m.dot(mjH2m); COM MHpp=mj1p.dot(mjH2p); COM MHpm=mj1p.dot(mjH2m); // Currents with pg CCurrent jgam,jgap,j2gm,j2gp; j2gp=joo(p1out,true,pg,true); j2gm=joo(p1out,false,pg,false); jgap=j(pg,true,p1in,true); jgam=j(pg,false,p1in,false); CCurrent qsum(q1+qg); CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in); CCurrent p1o(p1out); CCurrent p1i(p1in); Lmm=(qsum*(MHmm) + (-2.*mjH2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHmm/2.))/q1.m2(); Lmp=(qsum*(MHmp) + (-2.*mjH2p.dot(pg))*mj1m+2.*mj1m.dot(pg)*mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHmp/2.))/q1.m2(); Lpm=(qsum*(MHpm) + (-2.*mjH2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHpm/2.))/q1.m2(); Lpp=(qsum*(MHpp) + (-2.*mjH2p.dot(pg))*mj1p+2.*mj1p.dot(pg)*mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHpp/2.))/q1.m2(); U1mm=(jgam.dot(mjH2m)*j2gm+2.*p1o*MHmm)/(p1out+pg).m2(); U1mp=(jgam.dot(mjH2p)*j2gm+2.*p1o*MHmp)/(p1out+pg).m2(); U1pm=(jgap.dot(mjH2m)*j2gp+2.*p1o*MHpm)/(p1out+pg).m2(); U1pp=(jgap.dot(mjH2p)*j2gp+2.*p1o*MHpp)/(p1out+pg).m2(); U2mm=((-1.)*j2gm.dot(mjH2m)*jgam+2.*p1i*MHmm)/(p1in-pg).m2(); U2mp=((-1.)*j2gm.dot(mjH2p)*jgam+2.*p1i*MHmp)/(p1in-pg).m2(); U2pm=((-1.)*j2gp.dot(mjH2m)*jgap+2.*p1i*MHpm)/(p1in-pg).m2(); U2pp=((-1.)*j2gp.dot(mjH2p)*jgap+2.*p1i*MHpp)/(p1in-pg).m2(); const double cf=RHEJ::C_F; double amm,amp,apm,app; amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm); amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp); apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm); app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp); double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2()); // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH1.m2()*qg.m2(); ampsq/=th; ampsq/=16.; ampsq*=4.*4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ. //Higgs coupling is included in Hjets.C return ampsq; } double jM2unogqbarHQbar (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { // This construction is taking rapidity order: pg > p1out >> p2out // std::cerr<<"This Uno Current: "<<p1out<<" "<<p1in<<" "<<p2out<<" "<<p2in<<" "<<pg<<std::endl; CLHEP::HepLorentzVector q1=p1in-p1out; // Top End CLHEP::HepLorentzVector q2=-(p2in-p2out); // Bottom End CLHEP::HepLorentzVector qg=p1in-p1out-pg; // Extra bit post-gluon // std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl; CCurrent mj1m,mj1p,mj2m,mj2p,mjH2m,mjH2p; mj1p=jio(p1in,true,p1out,true); mj1m=jio(p1in,false,p1out,false); mjH2p=jioH(p2in,true,p2out,true,qH1,qH2, mt, incBot, mb); mjH2m=jioH(p2in,false,p2out,false,qH1,qH2, mt, incBot, mb); // Dot products of these which occur again and again COM MHmp=mj1m.dot(mjH2p); // And now for the Higgs ones COM MHmm=mj1m.dot(mjH2m); COM MHpp=mj1p.dot(mjH2p); COM MHpm=mj1p.dot(mjH2m); // Currents with pg CCurrent jgam,jgap,j2gm,j2gp; j2gp=joo(pg,true,p1out,true); j2gm=joo(pg,false,p1out,false); jgap=jio(p1in,true,pg,true); jgam=jio(p1in,false,pg,false); CCurrent qsum(q1+qg); CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in); CCurrent p1o(p1out); CCurrent p1i(p1in); Lmm=(qsum*(MHmm) + (-2.*mjH2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHmm/2.))/q1.m2(); Lmp=(qsum*(MHmp) + (-2.*mjH2p.dot(pg))*mj1m+2.*mj1m.dot(pg)*mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHmp/2.))/q1.m2(); Lpm=(qsum*(MHpm) + (-2.*mjH2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHpm/2.))/q1.m2(); Lpp=(qsum*(MHpp) + (-2.*mjH2p.dot(pg))*mj1p+2.*mj1p.dot(pg)*mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHpp/2.))/q1.m2(); U1mm=(jgam.dot(mjH2m)*j2gm+2.*p1o*MHmm)/(p1out+pg).m2(); U1mp=(jgam.dot(mjH2p)*j2gm+2.*p1o*MHmp)/(p1out+pg).m2(); U1pm=(jgap.dot(mjH2m)*j2gp+2.*p1o*MHpm)/(p1out+pg).m2(); U1pp=(jgap.dot(mjH2p)*j2gp+2.*p1o*MHpp)/(p1out+pg).m2(); U2mm=((-1.)*j2gm.dot(mjH2m)*jgam+2.*p1i*MHmm)/(p1in-pg).m2(); U2mp=((-1.)*j2gm.dot(mjH2p)*jgam+2.*p1i*MHmp)/(p1in-pg).m2(); U2pm=((-1.)*j2gp.dot(mjH2m)*jgap+2.*p1i*MHpm)/(p1in-pg).m2(); U2pp=((-1.)*j2gp.dot(mjH2p)*jgap+2.*p1i*MHpp)/(p1in-pg).m2(); const double cf=RHEJ::C_F; double amm,amp,apm,app; amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm); amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp); apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm); app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp); double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2()); // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH1.m2()*qg.m2(); ampsq/=th; ampsq/=16.; //Higgs coupling is included in Hjets.C ampsq*=4.*4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ. return ampsq; } double jM2unogqHg (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { // This construction is taking rapidity order: pg > p1out >> p2out // std::cerr<<"This Uno Current: "<<p1out<<" "<<p1in<<" "<<p2out<<" "<<p2in<<" "<<pg<<std::endl; CLHEP::HepLorentzVector q1=p1in-p1out; // Top End CLHEP::HepLorentzVector q2=-(p2in-p2out); // Bottom End CLHEP::HepLorentzVector qg=p1in-p1out-pg; // Extra bit post-gluon CCurrent mj1m,mj1p,mj2m,mj2p,mjH2m,mjH2p; mj1p=j(p1out,true,p1in,true); mj1m=j(p1out,false,p1in,false); mjH2p=jH(p2out,true,p2in,true,qH1,qH2, mt, incBot, mb); mjH2m=jH(p2out,false,p2in,false,qH1,qH2, mt, incBot, mb); // Dot products of these which occur again and again COM MHmp=mj1m.dot(mjH2p); // And now for the Higgs ones COM MHmm=mj1m.dot(mjH2m); COM MHpp=mj1p.dot(mjH2p); COM MHpm=mj1p.dot(mjH2m); // Currents with pg CCurrent jgam,jgap,j2gm,j2gp; j2gp=joo(p1out,true,pg,true); j2gm=joo(p1out,false,pg,false); jgap=j(pg,true,p1in,true); jgam=j(pg,false,p1in,false); CCurrent qsum(q1+qg); CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in); CCurrent p1o(p1out); CCurrent p1i(p1in); Lmm=(qsum*(MHmm) + (-2.*mjH2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHmm/2.))/q1.m2(); Lmp=(qsum*(MHmp) + (-2.*mjH2p.dot(pg))*mj1m+2.*mj1m.dot(pg)*mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHmp/2.))/q1.m2(); Lpm=(qsum*(MHpm) + (-2.*mjH2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHpm/2.))/q1.m2(); Lpp=(qsum*(MHpp) + (-2.*mjH2p.dot(pg))*mj1p+2.*mj1p.dot(pg)*mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHpp/2.))/q1.m2(); U1mm=(jgam.dot(mjH2m)*j2gm+2.*p1o*MHmm)/(p1out+pg).m2(); U1mp=(jgam.dot(mjH2p)*j2gm+2.*p1o*MHmp)/(p1out+pg).m2(); U1pm=(jgap.dot(mjH2m)*j2gp+2.*p1o*MHpm)/(p1out+pg).m2(); U1pp=(jgap.dot(mjH2p)*j2gp+2.*p1o*MHpp)/(p1out+pg).m2(); U2mm=((-1.)*j2gm.dot(mjH2m)*jgam+2.*p1i*MHmm)/(p1in-pg).m2(); U2mp=((-1.)*j2gm.dot(mjH2p)*jgam+2.*p1i*MHmp)/(p1in-pg).m2(); U2pm=((-1.)*j2gp.dot(mjH2m)*jgap+2.*p1i*MHpm)/(p1in-pg).m2(); U2pp=((-1.)*j2gp.dot(mjH2p)*jgap+2.*p1i*MHpp)/(p1in-pg).m2(); const double cf=RHEJ::C_F; double amm,amp,apm,app; amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm); amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp); apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm); app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp); double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2()); // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH1.m2()*qg.m2(); ampsq/=th; ampsq/=16.; ampsq*=4./9.*4./9.; // Factor of (Cf/Ca) for each quark to match MH2qQ. // here we need 2 to match with the normalization // gq is 9./4. times the qQ //Higgs coupling is included in Hjets.C const double K = K_g(p2out, p2in); return ampsq*K/C_A*9./4.; //ca/cf = 9/4 } double jM2unogqbarHg (CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { // This construction is taking rapidity order: pg > p1out >> p2out // std::cerr<<"This Uno Current: "<<p1out<<" "<<p1in<<" "<<p2out<<" "<<p2in<<" "<<pg<<std::endl; CLHEP::HepLorentzVector q1=p1in-p1out; // Top End CLHEP::HepLorentzVector q2=-(p2in-p2out); // Bottom End CLHEP::HepLorentzVector qg=p1in-p1out-pg; // Extra bit post-gluon CCurrent mj1m,mj1p,mj2m,mj2p,mjH2m,mjH2p; mj1p=jio(p1in,true,p1out,true); mj1m=jio(p1in,false,p1out,false); mjH2p=jH(p2out,true,p2in,true,qH1,qH2, mt, incBot, mb); mjH2m=jH(p2out,false,p2in,false,qH1,qH2, mt, incBot, mb); // Dot products of these which occur again and again COM MHmp=mj1m.dot(mjH2p); // And now for the Higgs ones COM MHmm=mj1m.dot(mjH2m); COM MHpp=mj1p.dot(mjH2p); COM MHpm=mj1p.dot(mjH2m); // Currents with pg CCurrent jgam,jgap,j2gm,j2gp; j2gp=joo(pg,true,p1out,true); j2gm=joo(pg,false,p1out,false); jgap=jio(p1in,true,pg,true); jgam=jio(p1in,false,pg,false); CCurrent qsum(q1+qg); CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p2o(p2out),p2i(p2in); CCurrent p1o(p1out); CCurrent p1i(p1in); Lmm=(qsum*(MHmm) + (-2.*mjH2m.dot(pg))*mj1m+2.*mj1m.dot(pg)*mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHmm/2.))/q1.m2(); Lmp=(qsum*(MHmp) + (-2.*mjH2p.dot(pg))*mj1m+2.*mj1m.dot(pg)*mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHmp/2.))/q1.m2(); Lpm=(qsum*(MHpm) + (-2.*mjH2m.dot(pg))*mj1p+2.*mj1p.dot(pg)*mjH2m+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHpm/2.))/q1.m2(); Lpp=(qsum*(MHpp) + (-2.*mjH2p.dot(pg))*mj1p+2.*mj1p.dot(pg)*mjH2p+(p2o/pg.dot(p2out) + p2i/pg.dot(p2in))*(qg.m2()*MHpp/2.))/q1.m2(); U1mm=(jgam.dot(mjH2m)*j2gm+2.*p1o*MHmm)/(p1out+pg).m2(); U1mp=(jgam.dot(mjH2p)*j2gm+2.*p1o*MHmp)/(p1out+pg).m2(); U1pm=(jgap.dot(mjH2m)*j2gp+2.*p1o*MHpm)/(p1out+pg).m2(); U1pp=(jgap.dot(mjH2p)*j2gp+2.*p1o*MHpp)/(p1out+pg).m2(); U2mm=((-1.)*j2gm.dot(mjH2m)*jgam+2.*p1i*MHmm)/(p1in-pg).m2(); U2mp=((-1.)*j2gm.dot(mjH2p)*jgam+2.*p1i*MHmp)/(p1in-pg).m2(); U2pm=((-1.)*j2gp.dot(mjH2m)*jgap+2.*p1i*MHpm)/(p1in-pg).m2(); U2pp=((-1.)*j2gp.dot(mjH2p)*jgap+2.*p1i*MHpp)/(p1in-pg).m2(); const double cf=RHEJ::C_F; double amm,amp,apm,app; amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm); amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp); apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm); app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp); double ampsq=-(amm+amp+apm+app)/(q2.m2()*qH2.m2()); // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH1.m2()*qg.m2(); ampsq/=th; ampsq/=16.; ampsq*=4./9.*4./9.; // Factor of (Cf/Ca) for each quark to match MH2qQ. // here we need 2 to match with the normalization // gq is 9./4. times the qQ //Higgs coupling is included in Hjets.C const double K = K_g(p2out, p2in); return ampsq*K/C_F; } double jM2unobqHQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { // std::cout << "####################\n"; // std::cout << "# p1in : "<<p1in<< " "<<p1in.plus()<<" "<<p1in.minus()<<std::endl; // std::cout << "# p2in : "<<p2in<< " "<<p2in.plus()<<" "<<p2in.minus()<<std::endl; // std::cout << "# p1out : "<<p1out<< " "<<p1out.rapidity()<<std::endl; // std::cout << "# (qH1-qH2) : "<<(qH1-qH2)<< " "<<(qH1-qH2).rapidity()<<std::endl; // std::cout << "# pg : "<<pg<< " "<<pg.rapidity()<<std::endl; // std::cout << "# p2out : "<<p2out<< " "<<p2out.rapidity()<<std::endl; // std::cout << "####################\n"; CLHEP::HepLorentzVector q1=p1in-p1out; // Top End CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End // std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl; CCurrent mjH1m,mjH1p,mj2m,mj2p; mjH1p=jHtop(p1out,true,p1in,true,qH1,qH2, mt, incBot, mb); mjH1m=jHtop(p1out,false,p1in,false,qH1,qH2, mt, incBot, mb); mj2p=j(p2out,true,p2in,true); mj2m=j(p2out,false,p2in,false); // Dot products of these which occur again and again COM MHmp=mjH1m.dot(mj2p); // And now for the Higgs ones COM MHmm=mjH1m.dot(mj2m); COM MHpp=mjH1p.dot(mj2p); COM MHpm=mjH1p.dot(mj2m); // Currents with pg CCurrent jgbm,jgbp,j2gm,j2gp; j2gp=joo(p2out,true,pg,true); j2gm=joo(p2out,false,pg,false); jgbp=j(pg,true,p2in,true); jgbm=j(pg,false,p2in,false); CCurrent qsum(q2+q3); CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in); CCurrent p2o(p2out); CCurrent p2i(p2in); CCurrent pplus((p1in+p1out)/2.); CCurrent pminus((p2in+p2out)/2.); // COM test=pminus.dot(p1in); Lmm=((-1.)*qsum*(MHmm) + (-2.*mjH1m.dot(pg))*mj2m+2.*mj2m.dot(pg)*mjH1m + (p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHmm/2.))/q3.m2(); Lmp=((-1.)*qsum*(MHmp) + (-2.*mjH1m.dot(pg))*mj2p+2.*mj2p.dot(pg)*mjH1m +(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHmp/2.))/q3.m2(); Lpm=((-1.)*qsum*(MHpm) + (-2.*mjH1p.dot(pg))*mj2m+2.*mj2m.dot(pg)*mjH1p +(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHpm/2.))/q3.m2(); Lpp=((-1.)*qsum*(MHpp) + (-2.*mjH1p.dot(pg))*mj2p+2.*mj2p.dot(pg)*mjH1p +(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHpp/2.))/q3.m2(); U1mm=(jgbm.dot(mjH1m)*j2gm+2.*p2o*MHmm)/(p2out+pg).m2(); U1mp=(jgbp.dot(mjH1m)*j2gp+2.*p2o*MHmp)/(p2out+pg).m2(); U1pm=(jgbm.dot(mjH1p)*j2gm+2.*p2o*MHpm)/(p2out+pg).m2(); U1pp=(jgbp.dot(mjH1p)*j2gp+2.*p2o*MHpp)/(p2out+pg).m2(); U2mm=((-1.)*j2gm.dot(mjH1m)*jgbm+2.*p2i*MHmm)/(p2in-pg).m2(); U2mp=((-1.)*j2gp.dot(mjH1m)*jgbp+2.*p2i*MHmp)/(p2in-pg).m2(); U2pm=((-1.)*j2gm.dot(mjH1p)*jgbm+2.*p2i*MHpm)/(p2in-pg).m2(); U2pp=((-1.)*j2gp.dot(mjH1p)*jgbp+2.*p2i*MHpp)/(p2in-pg).m2(); const double cf=RHEJ::C_F; double amm,amp,apm,app; amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm); amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp); apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm); app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp); // 1/3. = 1/C_A ? double ampsq=-(amm+amp+apm+app)/(q1.m2()*qH1.m2()); // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs const double th=qH2.m2()*q2.m2(); ampsq/=th; ampsq/=16.; ampsq*=RHEJ::C_F*RHEJ::C_F/(RHEJ::C_A*RHEJ::C_A); // Factor of (Cf/Ca) for each quark to match MH2qQ. return ampsq; } double jM2unobqbarHQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { CLHEP::HepLorentzVector q1=p1in-p1out; // Top End CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End // std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl; CCurrent mjH1m,mjH1p,mj2m,mj2p; mjH1p=jioHtop(p1in,true,p1out,true,qH1,qH2, mt, incBot, mb); mjH1m=jioHtop(p1in,false,p1out,false,qH1,qH2, mt, incBot, mb); mj2p=j(p2out,true,p2in,true); mj2m=j(p2out,false,p2in,false); // Dot products of these which occur again and again COM MHmp=mjH1m.dot(mj2p); // And now for the Higgs ones COM MHmm=mjH1m.dot(mj2m); COM MHpp=mjH1p.dot(mj2p); COM MHpm=mjH1p.dot(mj2m); // Currents with pg CCurrent jgbm,jgbp,j2gm,j2gp; j2gp=joo(p2out,true,pg,true); j2gm=joo(p2out,false,pg,false); jgbp=j(pg,true,p2in,true); jgbm=j(pg,false,p2in,false); CCurrent qsum(q2+q3); CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in); CCurrent p2o(p2out); CCurrent p2i(p2in); CCurrent pplus((p1in+p1out)/2.); CCurrent pminus((p2in+p2out)/2.); // COM test=pminus.dot(p1in); Lmm=((-1.)*qsum*(MHmm) + (-2.*mjH1m.dot(pg))*mj2m+2.*mj2m.dot(pg)*mjH1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHmm/2.))/q3.m2(); Lmp=((-1.)*qsum*(MHmp) + (-2.*mjH1m.dot(pg))*mj2p+2.*mj2p.dot(pg)*mjH1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHmp/2.))/q3.m2(); Lpm=((-1.)*qsum*(MHpm) + (-2.*mjH1p.dot(pg))*mj2m+2.*mj2m.dot(pg)*mjH1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHpm/2.))/q3.m2(); Lpp=((-1.)*qsum*(MHpp) + (-2.*mjH1p.dot(pg))*mj2p+2.*mj2p.dot(pg)*mjH1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHpp/2.))/q3.m2(); U1mm=(jgbm.dot(mjH1m)*j2gm+2.*p2o*MHmm)/(p2out+pg).m2(); U1mp=(jgbp.dot(mjH1m)*j2gp+2.*p2o*MHmp)/(p2out+pg).m2(); U1pm=(jgbm.dot(mjH1p)*j2gm+2.*p2o*MHpm)/(p2out+pg).m2(); U1pp=(jgbp.dot(mjH1p)*j2gp+2.*p2o*MHpp)/(p2out+pg).m2(); U2mm=((-1.)*j2gm.dot(mjH1m)*jgbm+2.*p2i*MHmm)/(p2in-pg).m2(); U2mp=((-1.)*j2gp.dot(mjH1m)*jgbp+2.*p2i*MHmp)/(p2in-pg).m2(); U2pm=((-1.)*j2gm.dot(mjH1p)*jgbm+2.*p2i*MHpm)/(p2in-pg).m2(); U2pp=((-1.)*j2gp.dot(mjH1p)*jgbp+2.*p2i*MHpp)/(p2in-pg).m2(); const double cf=RHEJ::C_F; double amm,amp,apm,app; amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm); amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp); apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm); app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp); double ampsq=-(amm+amp+apm+app)/(q1.m2()*qH1.m2()); // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH2.m2()*q2.m2(); ampsq/=th; ampsq/=16.; ampsq*=4.*4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ. //Higgs coupling is included in Hjets.C return ampsq; } double jM2unobqHQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { CLHEP::HepLorentzVector q1=p1in-p1out; // Top End CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End // std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl; CCurrent mjH1m,mjH1p,mj2m,mj2p; mjH1p=jHtop(p1out,true,p1in,true,qH1,qH2,mt, incBot, mb); mjH1m=jHtop(p1out,false,p1in,false,qH1,qH2,mt, incBot, mb); mj2p=jio(p2in,true,p2out,true); mj2m=jio(p2in,false,p2out,false); // Dot products of these which occur again and again COM MHmp=mjH1m.dot(mj2p); // And now for the Higgs ones COM MHmm=mjH1m.dot(mj2m); COM MHpp=mjH1p.dot(mj2p); COM MHpm=mjH1p.dot(mj2m); // Currents with pg CCurrent jgbm,jgbp,j2gm,j2gp; j2gp=joo(pg,true,p2out,true); j2gm=joo(pg,false,p2out,false); jgbp=jio(p2in,true,pg,true); jgbm=jio(p2in,false,pg,false); CCurrent qsum(q2+q3); CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in); CCurrent p2o(p2out); CCurrent p2i(p2in); CCurrent pplus((p1in+p1out)/2.); CCurrent pminus((p2in+p2out)/2.); // COM test=pminus.dot(p1in); Lmm=((-1.)*qsum*(MHmm) + (-2.*mjH1m.dot(pg))*mj2m+2.*mj2m.dot(pg)*mjH1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHmm/2.))/q3.m2(); Lmp=((-1.)*qsum*(MHmp) + (-2.*mjH1m.dot(pg))*mj2p+2.*mj2p.dot(pg)*mjH1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHmp/2.))/q3.m2(); Lpm=((-1.)*qsum*(MHpm) + (-2.*mjH1p.dot(pg))*mj2m+2.*mj2m.dot(pg)*mjH1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHpm/2.))/q3.m2(); Lpp=((-1.)*qsum*(MHpp) + (-2.*mjH1p.dot(pg))*mj2p+2.*mj2p.dot(pg)*mjH1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHpp/2.))/q3.m2(); U1mm=(jgbm.dot(mjH1m)*j2gm+2.*p2o*MHmm)/(p2out+pg).m2(); U1mp=(jgbp.dot(mjH1m)*j2gp+2.*p2o*MHmp)/(p2out+pg).m2(); U1pm=(jgbm.dot(mjH1p)*j2gm+2.*p2o*MHpm)/(p2out+pg).m2(); U1pp=(jgbp.dot(mjH1p)*j2gp+2.*p2o*MHpp)/(p2out+pg).m2(); U2mm=((-1.)*j2gm.dot(mjH1m)*jgbm+2.*p2i*MHmm)/(p2in-pg).m2(); U2mp=((-1.)*j2gp.dot(mjH1m)*jgbp+2.*p2i*MHmp)/(p2in-pg).m2(); U2pm=((-1.)*j2gm.dot(mjH1p)*jgbm+2.*p2i*MHpm)/(p2in-pg).m2(); U2pp=((-1.)*j2gp.dot(mjH1p)*jgbp+2.*p2i*MHpp)/(p2in-pg).m2(); const double cf=RHEJ::C_F; double amm,amp,apm,app; amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm); amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp); apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm); app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp); double ampsq=-(amm+amp+apm+app)/(q1.m2()*qH1.m2()); // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH2.m2()*q2.m2(); ampsq/=th; ampsq/=16.; ampsq*=4.*4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ. //Higgs coupling is included in Hjets.C return ampsq; } double jM2unobqbarHQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { CLHEP::HepLorentzVector q1=p1in-p1out; // Top End CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End // std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl; CCurrent mjH1m,mjH1p,mj2m,mj2p; mjH1p=jioHtop(p1in,true,p1out,true,qH1,qH2,mt, incBot, mb); mjH1m=jioHtop(p1in,false,p1out,false,qH1,qH2,mt, incBot, mb); mj2p=jio(p2in,true,p2out,true); mj2m=jio(p2in,false,p2out,false); // Dot products of these which occur again and again COM MHmp=mjH1m.dot(mj2p); // And now for the Higgs ones COM MHmm=mjH1m.dot(mj2m); COM MHpp=mjH1p.dot(mj2p); COM MHpm=mjH1p.dot(mj2m); // Currents with pg CCurrent jgbm,jgbp,j2gm,j2gp; j2gp=joo(pg,true,p2out,true); j2gm=joo(pg,false,p2out,false); jgbp=jio(p2in,true,pg,true); jgbm=jio(p2in,false,pg,false); CCurrent qsum(q2+q3); CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in); CCurrent p2o(p2out); CCurrent p2i(p2in); CCurrent pplus((p1in+p1out)/2.); CCurrent pminus((p2in+p2out)/2.); // COM test=pminus.dot(p1in); Lmm=((-1.)*qsum*(MHmm) + (-2.*mjH1m.dot(pg))*mj2m+2.*mj2m.dot(pg)*mjH1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHmm/2.))/q3.m2(); Lmp=((-1.)*qsum*(MHmp) + (-2.*mjH1m.dot(pg))*mj2p+2.*mj2p.dot(pg)*mjH1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHmp/2.))/q3.m2(); Lpm=((-1.)*qsum*(MHpm) + (-2.*mjH1p.dot(pg))*mj2m+2.*mj2m.dot(pg)*mjH1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHpm/2.))/q3.m2(); Lpp=((-1.)*qsum*(MHpp) + (-2.*mjH1p.dot(pg))*mj2p+2.*mj2p.dot(pg)*mjH1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHpp/2.))/q3.m2(); U1mm=(jgbm.dot(mjH1m)*j2gm+2.*p2o*MHmm)/(p2out+pg).m2(); U1mp=(jgbp.dot(mjH1m)*j2gp+2.*p2o*MHmp)/(p2out+pg).m2(); U1pm=(jgbm.dot(mjH1p)*j2gm+2.*p2o*MHpm)/(p2out+pg).m2(); U1pp=(jgbp.dot(mjH1p)*j2gp+2.*p2o*MHpp)/(p2out+pg).m2(); U2mm=((-1.)*j2gm.dot(mjH1m)*jgbm+2.*p2i*MHmm)/(p2in-pg).m2(); U2mp=((-1.)*j2gp.dot(mjH1m)*jgbp+2.*p2i*MHmp)/(p2in-pg).m2(); U2pm=((-1.)*j2gm.dot(mjH1p)*jgbm+2.*p2i*MHpm)/(p2in-pg).m2(); U2pp=((-1.)*j2gp.dot(mjH1p)*jgbp+2.*p2i*MHpp)/(p2in-pg).m2(); const double cf=RHEJ::C_F; double amm,amp,apm,app; amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm); amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp); apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm); app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp); double ampsq=-(amm+amp+apm+app)/(q1.m2()*qH1.m2()); // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH2.m2()*q2.m2(); ampsq/=th; ampsq/=16.; ampsq*=4.*4./(9.*9.); // Factor of (Cf/Ca) for each quark to match MH2qQ. //Higgs coupling is included in Hjets.C return ampsq; } double jM2unobgHQg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { // std::cout << "####################\n"; // std::cout << "# p1in : "<<p1in<< " "<<p1in.plus()<<" "<<p1in.minus()<<std::endl; // std::cout << "# p2in : "<<p2in<< " "<<p2in.plus()<<" "<<p2in.minus()<<std::endl; // std::cout << "# p1out : "<<p1out<< " "<<p1out.rapidity()<<std::endl; // std::cout << "# (qH1-qH2) : "<<(qH1-qH2)<< " "<<(qH1-qH2).rapidity()<<std::endl; // std::cout << "# pg : "<<pg<< " "<<pg.rapidity()<<std::endl; // std::cout << "# p2out : "<<p2out<< " "<<p2out.rapidity()<<std::endl; // std::cout << "####################\n"; CLHEP::HepLorentzVector q1=p1in-p1out; // Top End CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End // std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl; CCurrent mjH1m,mjH1p,mj2m,mj2p; mjH1p=jHtop(p1out,true,p1in,true,qH1,qH2,mt, incBot, mb); mjH1m=jHtop(p1out,false,p1in,false,qH1,qH2,mt, incBot, mb); mj2p=j(p2out,true,p2in,true); mj2m=j(p2out,false,p2in,false); // Dot products of these which occur again and again COM MHmp=mjH1m.dot(mj2p); // And now for the Higgs ones COM MHmm=mjH1m.dot(mj2m); COM MHpp=mjH1p.dot(mj2p); COM MHpm=mjH1p.dot(mj2m); // Currents with pg CCurrent jgbm,jgbp,j2gm,j2gp; j2gp=joo(p2out,true,pg,true); j2gm=joo(p2out,false,pg,false); jgbp=j(pg,true,p2in,true); jgbm=j(pg,false,p2in,false); CCurrent qsum(q2+q3); CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in); CCurrent p2o(p2out); CCurrent p2i(p2in); CCurrent pplus((p1in+p1out)/2.); CCurrent pminus((p2in+p2out)/2.); // COM test=pminus.dot(p1in); Lmm=((-1.)*qsum*(MHmm) + (-2.*mjH1m.dot(pg))*mj2m+2.*mj2m.dot(pg)*mjH1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHmm/2.))/q3.m2(); Lmp=((-1.)*qsum*(MHmp) + (-2.*mjH1m.dot(pg))*mj2p+2.*mj2p.dot(pg)*mjH1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHmp/2.))/q3.m2(); Lpm=((-1.)*qsum*(MHpm) + (-2.*mjH1p.dot(pg))*mj2m+2.*mj2m.dot(pg)*mjH1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHpm/2.))/q3.m2(); Lpp=((-1.)*qsum*(MHpp) + (-2.*mjH1p.dot(pg))*mj2p+2.*mj2p.dot(pg)*mjH1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHpp/2.))/q3.m2(); U1mm=(jgbm.dot(mjH1m)*j2gm+2.*p2o*MHmm)/(p2out+pg).m2(); U1mp=(jgbp.dot(mjH1m)*j2gp+2.*p2o*MHmp)/(p2out+pg).m2(); U1pm=(jgbm.dot(mjH1p)*j2gm+2.*p2o*MHpm)/(p2out+pg).m2(); U1pp=(jgbp.dot(mjH1p)*j2gp+2.*p2o*MHpp)/(p2out+pg).m2(); U2mm=((-1.)*j2gm.dot(mjH1m)*jgbm+2.*p2i*MHmm)/(p2in-pg).m2(); U2mp=((-1.)*j2gp.dot(mjH1m)*jgbp+2.*p2i*MHmp)/(p2in-pg).m2(); U2pm=((-1.)*j2gm.dot(mjH1p)*jgbm+2.*p2i*MHpm)/(p2in-pg).m2(); U2pp=((-1.)*j2gp.dot(mjH1p)*jgbp+2.*p2i*MHpp)/(p2in-pg).m2(); const double cf=RHEJ::C_F; double amm,amp,apm,app; amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm); amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp); apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm); app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp); double ampsq=-(amm+amp+apm+app)/(q1.m2()*qH1.m2()); // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH2.m2()*q2.m2(); ampsq/=th; ampsq/=16.; ampsq*=4./9.*4./9.; // Factor of (Cf/Ca) for each quark to match MH2qQ. // need twice to match the normalization //Higgs coupling is included in Hjets.C const double K = K_g(p1out, p1in); return ampsq*K/C_F; } double jM2unobgHQbarg (CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector pg, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector qH1, CLHEP::HepLorentzVector qH2, double mt, bool incBot, double mb) { CLHEP::HepLorentzVector q1=p1in-p1out; // Top End CLHEP::HepLorentzVector q2=-(p2in-p2out-pg); // Extra bit pre-gluon CLHEP::HepLorentzVector q3=-(p2in-p2out); // Bottom End // std::cerr<<"Current: "<<q1.m2()<<" "<<q2.m2()<<std::endl; CCurrent mjH1m,mjH1p,mj2m,mj2p; mjH1p=jHtop(p1out,true,p1in,true,qH1,qH2,mt, incBot, mb); mjH1m=jHtop(p1out,false,p1in,false,qH1,qH2,mt, incBot, mb); mj2p=jio(p2in,true,p2out,true); mj2m=jio(p2in,false,p2out,false); // Dot products of these which occur again and again COM MHmp=mjH1m.dot(mj2p); // And now for the Higgs ones COM MHmm=mjH1m.dot(mj2m); COM MHpp=mjH1p.dot(mj2p); COM MHpm=mjH1p.dot(mj2m); // Currents with pg CCurrent jgbm,jgbp,j2gm,j2gp; j2gp=joo(pg,true,p2out,true); j2gm=joo(pg,false,p2out,false); jgbp=jio(p2in,true,pg,true); jgbm=jio(p2in,false,pg,false); CCurrent qsum(q2+q3); CCurrent Lmp,Lmm,Lpp,Lpm,U1mp,U1mm,U1pp,U1pm,U2mp,U2mm,U2pp,U2pm,p1o(p1out),p1i(p1in); CCurrent p2o(p2out); CCurrent p2i(p2in); CCurrent pplus((p1in+p1out)/2.); CCurrent pminus((p2in+p2out)/2.); // COM test=pminus.dot(p1in); Lmm=((-1.)*qsum*(MHmm) + (-2.*mjH1m.dot(pg))*mj2m+2.*mj2m.dot(pg)*mjH1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHmm/2.))/q3.m2(); Lmp=((-1.)*qsum*(MHmp) + (-2.*mjH1m.dot(pg))*mj2p+2.*mj2p.dot(pg)*mjH1m+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHmp/2.))/q3.m2(); Lpm=((-1.)*qsum*(MHpm) + (-2.*mjH1p.dot(pg))*mj2m+2.*mj2m.dot(pg)*mjH1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHpm/2.))/q3.m2(); Lpp=((-1.)*qsum*(MHpp) + (-2.*mjH1p.dot(pg))*mj2p+2.*mj2p.dot(pg)*mjH1p+(p1o/pg.dot(p1out) + p1i/pg.dot(p1in))*(q2.m2()*MHpp/2.))/q3.m2(); U1mm=(jgbm.dot(mjH1m)*j2gm+2.*p2o*MHmm)/(p2out+pg).m2(); U1mp=(jgbp.dot(mjH1m)*j2gp+2.*p2o*MHmp)/(p2out+pg).m2(); U1pm=(jgbm.dot(mjH1p)*j2gm+2.*p2o*MHpm)/(p2out+pg).m2(); U1pp=(jgbp.dot(mjH1p)*j2gp+2.*p2o*MHpp)/(p2out+pg).m2(); U2mm=((-1.)*j2gm.dot(mjH1m)*jgbm+2.*p2i*MHmm)/(p2in-pg).m2(); U2mp=((-1.)*j2gp.dot(mjH1m)*jgbp+2.*p2i*MHmp)/(p2in-pg).m2(); U2pm=((-1.)*j2gm.dot(mjH1p)*jgbm+2.*p2i*MHpm)/(p2in-pg).m2(); U2pp=((-1.)*j2gp.dot(mjH1p)*jgbp+2.*p2i*MHpp)/(p2in-pg).m2(); const double cf=RHEJ::C_F; double amm,amp,apm,app; amm=cf*(2.*vre(Lmm-U1mm,Lmm+U2mm))+2.*cf*cf/3.*vabs2(U1mm+U2mm); amp=cf*(2.*vre(Lmp-U1mp,Lmp+U2mp))+2.*cf*cf/3.*vabs2(U1mp+U2mp); apm=cf*(2.*vre(Lpm-U1pm,Lpm+U2pm))+2.*cf*cf/3.*vabs2(U1pm+U2pm); app=cf*(2.*vre(Lpp-U1pp,Lpp+U2pp))+2.*cf*cf/3.*vabs2(U1pp+U2pp); double ampsq=-(amm+amp+apm+app)/(q1.m2()*qH1.m2()); // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) > 1.0000001) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // if ((vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm))/(2*vre(Lmm-U1mm,Lmm+U2mm)) < 0.9999999) // std::cout << " Big Problem!! " << vabs2(Lmm-U1mm+U2mm)+vabs2(Lmm)-vabs2(U1mm)-vabs2(U2mm) << " " << 2*vre(Lmm-U1mm,Lmm+U2mm) << std::endl; // Now add the t-channels for the Higgs double th=qH2.m2()*q2.m2(); ampsq/=th; ampsq/=16.; ampsq*=4./9.*4./9.; // Factor of (Cf/Ca) for each quark to match MH2qQ. //Higgs coupling is included in Hjets.C const double K = K_g(p1out, p1in); return ampsq*K/C_F; //ca/cf = 9/4 } // Begin finite mass stuff #ifdef RHEJ_BUILD_WITH_QCDLOOP namespace { // All the stuff needed for the box functions in qg->qgH now... //COM E1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM E1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //CLHEP::HepLorentzVector q2=k3+k4; CLHEP::HepLorentzVector q2=-(k1+k2+kh); double Delta, Sigma, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); //s34 = 2.*k3.dot(k4); s34 = q2.m2(); Delta = s12*s34 - S1*S2; Sigma = 4.*s12*s34 - pow(S1+S2,2); return looprwfactor*(-s12*D0DD(k2, k1, q2, mq)*(1 - 8.*mq*mq/s12 + S2/(2.*s12) + S2*(s12 - 8.*mq*mq)*(s34 + S1)/(2.*s12*Delta) + 2.*(s34 + S1)*(s34 + S1)/Delta + S2*pow((s34 + S1),3)/Delta/Delta) - ((s12 + S2)*C0DD(k2, k1 + q2, mq) - s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) - S1*C0DD(k1, q2, mq))*(S2*(s12 - 4.*mq*mq)/(2.*s12*Delta) + 2.*(s34 + S1)/Delta + S2*pow((s34 + S1),2)/Delta/Delta) + (C0DD(k1, q2, mq) - C0DD(k1 + k2, q2, mq))*(1. - 4.*mq*mq/s12) - C0DD(k1 + k2, q2, mq)*2.*s34/ S1 - (B0DD(k1 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2.*s34*(s34 + S1)/(S1*Delta) + (B0DD(q2, mq) - B0DD(k1 + k2 + q2, mq) + s12*C0DD(k1 + k2, q2, mq))*(2.*s34*(s34 + S1)*(S1 - S2)/(Delta*Sigma) + 2.*s34*(s34 + S1)/(S1*Delta)) + (B0DD(k1 + k2, mq) - B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*2.*(s34 + S1)*(2.*s12*s34 - S2*(S1 + S2))/(Delta*Sigma)); } //COM F1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM F1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //CLHEP::HepLorentzVector q2=k3+k4; CLHEP::HepLorentzVector q2 = -(k1+k2+kh); double Delta, Sigma, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); //s34 = 2.*k3.dot(k4); s34 = q2.m2(); Delta = s12*s34 - S1*S2; Sigma = 4.*s12*s34 - pow(S1+S2,2); return looprwfactor*(-S2*D0DD(k1, k2, q2, mq)*(0.5 - (s12 - 8.*mq*mq)*(s34 + S2)/(2.*Delta) - s12*pow((s34 + S2),3)/Delta/Delta) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) - s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) - S2*C0DD(k2, q2, mq))*(S2*(s12 - 4.*mq*mq)/(2.*s12*Delta) + S2*pow((s34 + S2),2)/Delta/Delta) - (C0DD(k1 + k2, q2, mq) - C0DD(k1, k2 + q2, mq))*(1. - 4.*mq*mq/s12) - C0DD(k1, k2 + q2, mq) + (B0DD(k2 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2.*pow((s34 + S2),2)/((s12 + S1)*Delta) - (B0DD( q2, mq) - B0DD(k1 + k2 + q2, mq) + s12*C0DD(k1 + k2, q2, mq))*2.*s34*(s34 + S2)*(S2 - S1)/(Delta*Sigma) + (B0DD( k1 + k2, mq) - B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*2.*(s34 + S2)*(2.*s12*s34 - S2*(S1 + S2))/(Delta*Sigma)); } //COM G1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM G1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //CLHEP::HepLorentzVector q2=k3+k4; CLHEP::HepLorentzVector q2 = -(k1+k2+kh); double Delta, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); //s34 = 2.*k3.dot(k4); s34 = q2.m2(); Delta = s12*s34 - S1*S2; return looprwfactor*(S2*D0DD(k1, q2, k2, mq)*(Delta/s12/s12 - 4.*mq*mq/s12) - S2*((s12 + S1)*C0DD(k1, k2 + q2, mq) - S1*C0DD(k1, q2, mq))*(1./ s12/s12 - (s12 - 4.*mq*mq)/(2.*s12*Delta)) - S2*((s12 + S2)*C0DD(k1 + q2, k2, mq) - S2*C0DD(k2, q2, mq))*(1./ s12/s12 + (s12 - 4.*mq*mq)/(2.*s12*Delta)) - C0DD(k1, q2, mq) - (C0DD(k1, k2 + q2, mq) - C0DD(k1, q2, mq))*4.*mq*mq/ s12 + (B0DD(k1 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2./ s12 + (B0DD(k1 + q2, mq) - B0DD(q2, mq))*2.*s34/(s12*S1) + (B0DD(k2 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2.*(s34 + S2)/(s12*(s12 + S1))); } //COM E4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM E4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //CLHEP::HepLorentzVector q2=k3+k4; CLHEP::HepLorentzVector q2 = -(k1+k2+kh); double Delta, Sigma, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); //s34 = 2.*k3.dot(k4); s34 = q2.m2(); Delta = s12*s34 - S1*S2; Sigma = 4.*s12*s34 - pow(S1+S2,2); return looprwfactor* (-s12*D0DD(k2, k1, q2, mq)*(0.5 - (S1 - 8.*mq*mq)*(s34 + S1)/(2.*Delta) - s12*pow((s34 + S1),3)/Delta/Delta) + ((s12 + S2)*C0DD(k2, k1 + q2, mq) - s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) - S1*C0DD(k1, q2, mq))*((S1 - 4.*mq*mq)/(2.*Delta) + s12*pow((s34 + S1),2)/Delta/Delta) - C0DD(k1 + k2, q2, mq) + (B0DD(k1 + q2, mq) - B0DD(k1 + k2 + q2, mq))*(2.*s34/Delta + 2.*s12*(s34 + S1)/((s12 + S2)*Delta)) - (B0DD( q2, mq) - B0DD(k1 + k2 + q2, mq) + s12*C0DD(k1 + k2, q2, mq))*((2.*s34*(2.*s12*s34 - S2*(S1 + S2) + s12*(S1 - S2)))/(Delta*Sigma)) + (B0DD(k1 + k2, mq) - B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*((2.*s12*(2.*s12*s34 - S1*(S1 + S2) + s34*(S2 - S1)))/(Delta*Sigma))); } //COM F4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM F4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //CLHEP::HepLorentzVector q2=k3+k4; CLHEP::HepLorentzVector q2 = -(k1+k2+kh); double Delta, Sigma, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); //s34 = 2.*k3.dot(k4); s34 = q2.m2(); Delta = s12*s34 - S1*S2; Sigma = 4.*s12*s34 - pow(S1+S2,2); return looprwfactor* (-s12*D0DD(k1, k2, q2, mq)*(0.5 + (S1 - 8.*mq*mq)*(s34 + S2)/(2.*Delta) + s12*pow((s34 + S2),3)/Delta/Delta) - ((s12 + S1)*C0DD(k1, k2 + q2, mq) - s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) - S2*C0DD(k2, q2, mq))*((S1 - 4.*mq*mq)/(2.*Delta) + s12*pow((s34 + S2),2)/Delta/Delta) - C0DD(k1 + k2, q2, mq) - (B0DD(k2 + q2, mq) - B0DD(k1 + k2 + q2, mq))*2.*(s34 + S2)/Delta + (B0DD(q2, mq) - B0DD(k1 + k2 + q2, mq) + s12*C0DD(k1 + k2, q2, mq))*2.*s34*(2.*s12*s34 - S1*(S1 + S2) + s12*(S2 - S1))/(Delta*Sigma) - (B0DD(k1 + k2, mq) - B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*(2.*s12*(2.*s12*s34 - S2*(S1 + S2) + s34*(S1 - S2))/(Delta*Sigma))); } //COM G4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM G4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //CLHEP::HepLorentzVector q2=k3+k4; CLHEP::HepLorentzVector q2 = -(k1+k2+kh); double Delta, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); //s34 = 2.*k3.dot(k4); s34 = q2.m2(); Delta = s12*s34 - S1*S2; return looprwfactor* (-D0DD(k1, q2, k2, mq)*(Delta/s12 + (s12 + S1)/2. - 4.*mq*mq) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) - S1*C0DD(k1, q2, mq))*(1./ s12 - (S1 - 4.*mq*mq)/(2.*Delta)) + ((s12 + S2)*C0DD( k1 + q2, k2, mq) - S2*C0DD(k2, q2, mq))*(1./ s12 + (S1 - 4.*mq*mq)/(2.*Delta)) + (B0DD( k1 + k2 + q2, mq) - B0DD(k1 + q2, mq))*2./(s12 + S2)); } //COM E10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM E10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //CLHEP::HepLorentzVector q2=k3+k4; CLHEP::HepLorentzVector q2 = -(k1+k2+kh); double Delta, Sigma, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); //s34 = 2.*k3.dot(k4); s34 = q2.m2(); Delta = s12*s34 - S1*S2; Sigma = 4.*s12*s34 - pow(S1+S2,2); return looprwfactor*(-s12*D0DD(k2, k1, q2, mq)*((s34 + S1)/Delta + 12.*mq*mq*S1*(s34 + S1)/Delta/Delta - 4.*s12*S1*pow((s34 + S1),3)/Delta/Delta/Delta) - ((s12 + S2)*C0DD(k2, k1 + q2, mq) - s12*C0DD(k1, k2, mq) + (S1 - S2)*C0DD(k1 + k2, q2, mq) - S1*C0DD(k1, q2, mq))*(1./Delta + 4.*mq*mq*S1/Delta/Delta - 4.*s12*S1*pow((s34 + S1),2)/Delta/Delta/Delta) + C0DD(k1 + k2, q2, mq)*(4.*s12*s34*(S1 - S2)/(Delta*Sigma) - 4.*(s12 - 2.*mq*mq)*(2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma)) + (B0DD(k1 + q2, mq) - B0DD(k1 + k2 + q2, mq))*(4.*(s34 + S1)/((s12 + S2)*Delta) + 8.*S1*(s34 + S1)/Delta/Delta) + (B0DD(q2, mq) - B0DD(k1 + k2 + q2, mq) + s12*C0DD(k1 + k2, q2, mq))*(12.*s34*(2.*s12 + S1 + S2)*(2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma*Sigma) - 4.*s34*(4.*s12 + 3.*S1 + S2)/(Delta*Sigma) + 8.*s12*s34*(s34*(s12 + S2) - S1*(s34 + S1))/(Delta*Delta*Sigma)) + (B0DD(k1 + k2, mq) - B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*(12.*s12*(2.*s34 + S1 + S2)*(2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma*Sigma) + 8.*s12*S1*(s34*(s12 + S2) - S1*(s34 + S1))/(Delta*Delta*Sigma))) + (COM(0.,1.)/(4.*M_PI*M_PI))*((2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma)); } //COM F10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM F10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //CLHEP::HepLorentzVector q2=k3+k4; CLHEP::HepLorentzVector q2 = -(k1+k2+kh); double Delta, Sigma, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); //s34 = 2.*k3.dot(k4); s34 = q2.m2(); Delta = s12*s34 - S1*S2; Sigma = 4.*s12*s34 - pow(S1+S2,2); return looprwfactor* (s12*D0DD(k1, k2, q2, mq)*((s34 + S2)/Delta - 4.*mq*mq/Delta + 12.*mq*mq*s34*(s12 + S1)/Delta/Delta - 4.*s12*pow((s34 + S2),2)/Delta/Delta - 4.*s12*S1*pow((s34 + S2),3)/Delta/Delta/Delta) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) - s12*C0DD(k1, k2, mq) - (S1 - S2)*C0DD(k1 + k2, q2, mq) - S2*C0DD(k2, q2, mq))*(1./Delta + 4.*mq*mq*S1/Delta/Delta - 4.*s12*(s34 + S2)/Delta/Delta - 4.*s12*S1*pow((s34 + S2),2)/Delta/Delta/Delta) - C0DD(k1 + k2, q2, mq)*(4.*s12*s34/(S2*Delta) + 4.*s12*s34*(S2 - S1)/(Delta*Sigma) + 4.*(s12 - 2.*mq*mq)*(2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma)) - (B0DD( k2 + q2, mq) - B0DD(k1 + k2 + q2, mq))*(4.*s34/(S2*Delta) + 8.*s34*(s12 + S1)/Delta/Delta) - (B0DD(q2, mq) - B0DD(k1 + k2 + q2, mq) + s12*C0DD(k1 + k2, q2, mq))*(-12*s34*(2*s12 + S1 + S2)*(2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma*Sigma) - 4.*s12*s34*s34/(S2*Delta*Delta) + 4.*s34*S1/(Delta*Sigma) - 4.*s34*(s12*s34*(2.*s12 + S2) - S1*S1*(2.*s12 + S1))/(Delta*Delta*Sigma)) - (B0DD(k1 + k2, mq) - B0DD(k1 + k2 + q2, mq) - (s34 + S1 + S2)*C0DD(k1 + k2, q2, mq))*(-12.*s12*(2.*s34 + S1 + S2)*(2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma*Sigma) + 8.*s12*(2.*s34 + S1)/(Delta*Sigma) - 8.*s12*s34*(2.*s12*s34 - S1*(S1 + S2) + s12*(S2 - S1))/(Delta*Delta*Sigma))) + (COM(0.,1.)/(4.*M_PI*M_PI))*((2.*s12*s34 - S1*(S1 + S2))/(Delta*Sigma)); } //COM G10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM G10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //CLHEP::HepLorentzVector q2=k3+k4; CLHEP::HepLorentzVector q2 = -(k1+k2+kh); double Delta, S1, S2, s12, s34; S1 = 2.*k1.dot(q2); S2 = 2.*k2.dot(q2); s12 = 2.*k1.dot(k2); //s34 = 2.*k3.dot(k4); s34 = q2.m2(); Delta = s12*s34 - S1*S2; return looprwfactor* (-D0DD(k1, q2, k2, mq)*(1. + 4.*S1*mq*mq/Delta) + ((s12 + S1)*C0DD(k1, k2 + q2, mq) - S1*C0DD(k1, q2, mq))*(1./Delta + 4.*S1*mq*mq/Delta/Delta) - ((s12 + S2)*C0DD(k1 + q2, k2, mq) - S2*C0DD(k2, q2, mq))*(1./Delta + 4.*S1*mq*mq/Delta/Delta) + (B0DD(k1 + k2 + q2, mq) - B0DD(k1 + q2, mq))*4.*(s34 + S1)/(Delta*(s12 + S2)) + (B0DD(q2, mq) - B0DD(k2 + q2, mq))*4.*s34/(Delta*S2)); } //COM H1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM H1(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //return E1(k1,k2,k3,k4,mq)+F1(k1,k2,k3,k4,mq)+G1(k1,k2,k3,k4,mq); return E1(k1,k2,kh,mq)+F1(k1,k2,kh,mq)+G1(k1,k2,kh,mq); } //COM H4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM H4(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //return E4(k1,k2,k3,k4,mq)+F4(k1,k2,k3,k4,mq)+G4(k1,k2,k3,k4,mq); return E4(k1,k2,kh,mq)+F4(k1,k2,kh,mq)+G4(k1,k2,kh,mq); } //COM H10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM H10(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //return E10(k1,k2,k3,k4,mq)+F10(k1,k2,k3,k4,mq)+G10(k1,k2,k3,k4,mq); return E10(k1,k2,kh,mq)+F10(k1,k2,kh,mq)+G10(k1,k2,kh,mq); } //COM H2(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM H2(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //return -1.*H1(k2,k1,k3,k4,mq); return -1.*H1(k2,k1,kh,mq); } //COM H5(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM H5(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //return -1.*H4(k2,k1,k3,k4,mq); return -1.*H4(k2,k1,kh,mq); } //COM H12(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector k3, CLHEP::HepLorentzVector k4, double mq) COM H12(CLHEP::HepLorentzVector k1, CLHEP::HepLorentzVector k2, CLHEP::HepLorentzVector kh, double mq) { //return -1.*H10(k2,k1,k3,k4,mq); return -1.*H10(k2,k1,kh,mq); } // FL and FT functions COM FL(CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mq) { CLHEP::HepLorentzVector Q = q1 + q2; double detQ2 = q1.m2()*q2.m2() - q1.dot(q2)*q1.dot(q2); return -1./(2.*detQ2)*((2.- 3.*q1.m2()*q2.dot(Q)/detQ2)*(B0DD(q1, mq) - B0DD(Q, mq)) + (2. - 3.*q2.m2()*q1.dot(Q)/detQ2)*(B0DD(q2, mq) - B0DD(Q, mq)) - (4.*mq*mq + q1.m2() + q2.m2() + Q.m2() - 3.*q1.m2()*q2.m2()*Q.m2()/detQ2)*C0DD( q1, q2, mq) - 2.); } COM FT(CLHEP::HepLorentzVector q1, CLHEP::HepLorentzVector q2, double mq) { CLHEP::HepLorentzVector Q = q1 + q2; double detQ2 = q1.m2()*q2.m2() - q1.dot(q2)*q1.dot(q2); return -1./(2.*detQ2)*(Q.m2()*(B0DD(q1, mq) + B0DD(q2, mq) - 2.*B0DD(Q, mq) - 2.*q1.dot(q2)*C0DD(q1, q2, mq)) + (q1.m2() - q2.m2()) *(B0DD(q1, mq) - B0DD(q2, mq))) - q1.dot(q2)*FL(q1, q2, mq); } CLHEP::HepLorentzVector ParityFlip(CLHEP::HepLorentzVector p) { CLHEP::HepLorentzVector flippedVector; flippedVector.setE(p.e()); flippedVector.setX(-p.x()); flippedVector.setY(-p.y()); flippedVector.setZ(-p.z()); return flippedVector; } /// @brief HC amp for qg->qgH with finite top (i.e. j^{++}_H) void g_gH_HC(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pH, double mq, current &retAns) { current cura1,pacur,p1cur,pHcur,conjeps1,conjepsH1,epsa,epsHa,epsHapart1, epsHapart2,conjepsH1part1,conjepsH1part2; COM ang1a,sqa1; const double F = 4.*mq*mq/v; // Easier to have the whole thing as current object so I can use cdot functionality. // Means I need to write pa,p1 as current objects to_current(pa, pacur); to_current(p1,p1cur); to_current(pH,pHcur); bool gluonforward = true; if(pa.z() < 0) gluonforward = false; //HEJ gauge jio(pa,false,p1,false,cura1); if(gluonforward){ // sqrt(2pa_-/p1_-)*p1_perp/abs(p1_perp) ang1a = sqrt(pa.plus()*p1.minus())*(p1.x()+COM(0.,1.)*p1.y())/p1.perp(); // sqrt(2pa_-/p1_-)*p1_perp*/abs(p1_perp) sqa1 = sqrt(pa.plus()*p1.minus())*(p1.x()-COM(0.,1.)*p1.y())/p1.perp(); } else { ang1a = sqrt(pa.minus()*p1.plus()); sqa1 = sqrt(pa.minus()*p1.plus()); } const double prop = (pa-p1-pH).m2(); cmult(-1./sqrt(2)/ang1a,cura1,conjeps1); cmult(1./sqrt(2)/sqa1,cura1,epsa); const COM Fta = FT(-pa,pa-pH,mq)/(pa-pH).m2(); const COM Ft1 = FT(-p1-pH,p1,mq)/(p1+pH).m2(); const COM h4 = H4(p1,-pa,pH,mq); const COM h5 = H5(p1,-pa,pH,mq); const COM h10 = H10(p1,-pa,pH,mq); const COM h12 = H12(p1,-pa,pH,mq); cmult(Fta*pa.dot(pH), epsa, epsHapart1); cmult(-1.*Fta*cdot(pHcur,epsa), pacur, epsHapart2); cmult(Ft1*cdot(pHcur,conjeps1), p1cur, conjepsH1part1); cmult(-Ft1*p1.dot(pH), conjeps1, conjepsH1part2); cadd(epsHapart1, epsHapart2, epsHa); cadd(conjepsH1part1, conjepsH1part2, conjepsH1); const COM aH1 = cdot(pHcur, cura1); current T1,T2,T3,T4,T5,T6,T7,T8,T9,T10; if(gluonforward){ cmult(sqrt(2.)*sqrt(p1.plus()/pa.plus())*prop/sqa1, conjepsH1, T1); cmult(-sqrt(2.)*sqrt(pa.plus()/p1.plus())*prop/ang1a, epsHa, T2); } else{ cmult(-sqrt(2.)*sqrt(p1.minus()/pa.minus()) *((p1.x()-COM(0.,1.)*p1.y())/p1.perp())*prop/sqa1, conjepsH1, T1); cmult(sqrt(2.)*sqrt(pa.minus()/p1.minus()) *((p1.x()-COM(0.,1.)*p1.y())/p1.perp())*prop/ang1a, epsHa, T2); } cmult(sqrt(2.)/ang1a*aH1, epsHa, T3); cmult(sqrt(2.)/sqa1*aH1, conjepsH1, T4); cmult(-sqrt(2.)*Fta*pa.dot(p1)*aH1/sqa1, conjeps1, T5); cmult(-sqrt(2.)*Ft1*pa.dot(p1)*aH1/ang1a, epsa, T6); cmult(-aH1/sqrt(2.)/sqa1*h4*8.*COM(0.,1.)*M_PI*M_PI, conjeps1, T7); cmult(aH1/sqrt(2.)/ang1a*h5*8.*COM(0.,1.)*M_PI*M_PI, epsa, T8); cmult(aH1*aH1/2./ang1a/sqa1*h10*8.*COM(0.,1.)*M_PI*M_PI, pacur, T9); cmult(-aH1*aH1/2./ang1a/sqa1*h12*8.*COM(0.,1.)*M_PI*M_PI, p1cur, T10); current ans; for(int i=0;i<4;i++) { ans[i] = T1[i]+T2[i]+T3[i]+T4[i]+T5[i]+T6[i]+T7[i]+T8[i]+T9[i]+T10[i]; } retAns[0] = F/prop*ans[0]; retAns[1] = F/prop*ans[1]; retAns[2] = F/prop*ans[2]; retAns[3] = F/prop*ans[3]; } /// @brief HNC amp for qg->qgH with finite top (i.e. j^{+-}_H) void g_gH_HNC(CLHEP::HepLorentzVector pa, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pH, double mq, current &retAns) { const double F = 4.*mq*mq/v; COM ang1a,sqa1; current conjepsH1,epsHa,p1cur,pacur,pHcur,conjeps1,epsa,paplusp1cur, p1minuspacur,cur1a,cura1,epsHapart1,epsHapart2,conjepsH1part1, conjepsH1part2; // Find here if pa, meaning the gluon, is forward or backward bool gluonforward = true; if(pa.z() < 0) gluonforward = false; jio(pa,true,p1,true,cura1); j(p1,true,pa,true,cur1a); to_current(pa,pacur); to_current(p1,p1cur); to_current(pH,pHcur); to_current(pa+p1,paplusp1cur); to_current(p1-pa,p1minuspacur); const COM aH1 = cdot(pHcur,cura1); const COM oneHa = std::conj(aH1); // = cdot(pHcur,cur1a) if(gluonforward){ // sqrt(2pa_-/p1_-)*p1_perp/abs(p1_perp) ang1a = sqrt(pa.plus()*p1.minus())*(p1.x()+COM(0.,1.)*p1.y())/p1.perp(); // sqrt(2pa_-/p1_-)*p1_perp*/abs(p1_perp) sqa1 = sqrt(pa.plus()*p1.minus())*(p1.x()-COM(0.,1.)*p1.y())/p1.perp(); } else { ang1a = sqrt(pa.minus()*p1.plus()); sqa1 = sqrt(pa.minus()*p1.plus()); } const double prop = (pa-p1-pH).m2(); cmult(1./sqrt(2)/sqa1, cur1a, epsa); cmult(-1./sqrt(2)/sqa1, cura1, conjeps1); const COM phase = cdot(conjeps1, epsa); const COM Fta = FT(-pa,pa-pH,mq)/(pa-pH).m2(); const COM Ft1 = FT(-p1-pH,p1,mq)/(p1+pH).m2(); const COM Falpha = FT(p1-pa,pa-p1-pH,mq); const COM Fbeta = FL(p1-pa,pa-p1-pH,mq); const COM h1 = H1(p1,-pa, pH, mq); const COM h2 = H2(p1,-pa, pH, mq); const COM h4 = H4(p1,-pa, pH, mq); const COM h5 = H5(p1,-pa, pH, mq); const COM h10 = H10(p1,-pa, pH, mq); const COM h12 = H12(p1,-pa, pH, mq); cmult(Fta*pa.dot(pH), epsa, epsHapart1); cmult(-1.*Fta*cdot(pHcur,epsa), pacur, epsHapart2); cmult(Ft1*cdot(pHcur,conjeps1), p1cur, conjepsH1part1); cmult(-Ft1*p1.dot(pH), conjeps1, conjepsH1part2); cadd(epsHapart1, epsHapart2, epsHa); cadd(conjepsH1part1, conjepsH1part2, conjepsH1); current T1,T2,T3,T4,T5a,T5b,T6,T7,T8a,T8b,T9,T10,T11a, T11b,T12a,T12b,T13; if(gluonforward){ cmult(sqrt(2.)*sqrt(p1.plus()/pa.plus())*prop/sqa1, conjepsH1, T1); cmult(-sqrt(2.)*sqrt(pa.plus()/p1.plus())*prop/sqa1, epsHa, T2); } else{ cmult(-sqrt(2.)*sqrt(p1.minus()/pa.minus())*((p1.x()-COM(0.,1.)*p1.y())/p1.perp()) *prop/sqa1, conjepsH1, T1); cmult(sqrt(2.)*sqrt(pa.minus()/p1.minus())*((p1.x()+COM(0.,1.)*p1.y())/p1.perp()) *prop/sqa1, epsHa, T2); } const COM boxdiagFact = 8.*COM(0.,1.)*M_PI*M_PI; cmult(aH1*sqrt(2.)/sqa1, epsHa, T3); cmult(oneHa*sqrt(2.)/sqa1, conjepsH1, T4); cmult(-2.*phase*Fta*pa.dot(pH), p1cur, T5a); cmult(2.*phase*Ft1*p1.dot(pH), pacur, T5b); cmult(-sqrt(2.)*Fta*p1.dot(pa)*oneHa/sqa1, conjeps1, T6); cmult(-sqrt(2.)*Ft1*pa.dot(p1)*aH1/sqa1, epsa, T7); cmult(-boxdiagFact*phase*h2, pacur, T8a); cmult(boxdiagFact*phase*h1, p1cur, T8b); cmult(boxdiagFact*aH1/sqrt(2.)/sqa1*h5, epsa, T9); cmult(-boxdiagFact*oneHa/sqrt(2.)/sqa1*h4, conjeps1, T10); cmult(boxdiagFact*aH1*oneHa/2./sqa1/sqa1*h10, pacur, T11a); cmult(-boxdiagFact*aH1*oneHa/2./sqa1/sqa1*h12, p1cur, T11b); cmult(-phase/(pa-p1).m2()*Falpha*(p1-pa).dot(pa-p1-pH), paplusp1cur, T12a); cmult(phase/(pa-p1).m2()*Falpha*(pa+p1).dot(pa-p1-pH), p1minuspacur, T12b); cmult(-phase*Fbeta*(pa-p1-pH).m2(), paplusp1cur, T13); current ans; for(int i=0;i<4;i++) { ans[i] = T1[i]+T2[i]+T3[i]+T4[i]+T5a[i]+T5b[i]+T6[i]+T7[i]+T8a[i]+T8b[i]+T9[i]+T10[i]+T11a[i]+T11b[i]+T12a[i]+T12b[i]+T13[i]; } retAns[0] = F/prop*ans[0]; retAns[1] = F/prop*ans[1]; retAns[2] = F/prop*ans[2]; retAns[3] = F/prop*ans[3]; } } // namespace anonymous // JDC - new amplitude with Higgs emitted close to gluon with full mt effects. Keep usual HEJ-style function call double MH2gq_outsideH(CLHEP::HepLorentzVector p1out, CLHEP::HepLorentzVector p1in, CLHEP::HepLorentzVector p2out, CLHEP::HepLorentzVector p2in, CLHEP::HepLorentzVector pH, double mq, bool includeBottom, double mq2) { current cur2bplus,cur2bminus, cur2bplusFlip, cur2bminusFlip; current retAns,retAnsb; j(p2out,true,p2in,true,cur2bplus); j(p2out,false,p2in,false,cur2bminus); j(ParityFlip(p2out),true,ParityFlip(p2in),true,cur2bplusFlip); j(ParityFlip(p2out),false,ParityFlip(p2in),false,cur2bminusFlip); COM app1,app2,apm1,apm2; COM app3, app4, apm3, apm4; if(!includeBottom) { g_gH_HC(p1in,p1out,pH,mq,retAns); app1=cdot(retAns,cur2bplus); app2=cdot(retAns,cur2bminus); g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns); app3=cdot(retAns,cur2bplusFlip); app4=cdot(retAns,cur2bminusFlip); // And non-conserving bits g_gH_HNC(p1in,p1out,pH,mq,retAns); apm1=cdot(retAns,cur2bplus); apm2=cdot(retAns,cur2bminus); g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns); apm3=cdot(retAns,cur2bplusFlip); apm4=cdot(retAns,cur2bminusFlip); } else { g_gH_HC(p1in,p1out,pH,mq,retAns); g_gH_HC(p1in,p1out,pH,mq2,retAnsb); app1=cdot(retAns,cur2bplus) + cdot(retAnsb,cur2bplus); app2=cdot(retAns,cur2bminus) + cdot(retAnsb,cur2bminus); g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns); g_gH_HC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq2,retAnsb); app3=cdot(retAns,cur2bplusFlip) + cdot(retAnsb,cur2bplusFlip); app4=cdot(retAns,cur2bminusFlip) + cdot(retAnsb,cur2bminusFlip); // And non-conserving bits g_gH_HNC(p1in,p1out,pH,mq,retAns); g_gH_HNC(p1in,p1out,pH,mq2,retAnsb); apm1=cdot(retAns,cur2bplus) + cdot(retAnsb,cur2bplus); apm2=cdot(retAns,cur2bminus) + cdot(retAnsb,cur2bminus); g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq,retAns); g_gH_HNC(ParityFlip(p1in),ParityFlip(p1out),ParityFlip(pH),mq2,retAnsb); apm3=cdot(retAns,cur2bplusFlip) + cdot(retAnsb,cur2bplusFlip); apm4=cdot(retAns,cur2bminusFlip) + cdot(retAnsb,cur2bminusFlip); } return abs2(app1) + abs2(app2) + abs2(app3) + abs2(app4) + abs2(apm1) + abs2(apm2) + abs2(apm3) + abs2(apm4); } #endif // RHEJ_BUILD_WITH_QCDLOOP double C2gHgm(CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pH) { static double A=1./(3.*M_PI*v); // Implements Eq. (4.22) in hep-ph/0301013 with modifications to incoming plus momenta double s12,p1p,p2p; COM p1perp,p3perp,phperp; // Determine first whether this is the case p1p\sim php>>p3p og the opposite s12=p1.invariantMass2(-p2); if (p2.pz()>0.) { // case considered in hep-ph/0301013 p1p=p1.plus(); p2p=p2.plus(); } else { // opposite case p1p=p1.minus(); p2p=p2.minus(); } p1perp=p1.px()+COM(0,1)*p1.py(); phperp=pH.px()+COM(0,1)*pH.py(); p3perp=-(p1perp+phperp); COM temp=COM(0,1)*A/(2.*s12)*(p2p/p1p*conj(p1perp)*p3perp+p1p/p2p*p1perp*conj(p3perp)); temp=temp*conj(temp); return temp.real(); } double C2gHgp(CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pH) { static double A=1./(3.*M_PI*v); // Implements Eq. (4.23) in hep-ph/0301013 double s12,php,p1p,phm; COM p1perp,p3perp,phperp; // Determine first whether this is the case p1p\sim php>>p3p og the opposite s12=p1.invariantMass2(-p2); if (p2.pz()>0.) { // case considered in hep-ph/0301013 php=pH.plus(); phm=pH.minus(); p1p=p1.plus(); } else { // opposite case php=pH.minus(); phm=pH.plus(); p1p=p1.minus(); } p1perp=p1.px()+COM(0,1)*p1.py(); phperp=pH.px()+COM(0,1)*pH.py(); p3perp=-(p1perp+phperp); COM temp=-COM(0,1)*A/(2.*s12)*(conj(p1perp*p3perp)*pow(php/p1p,2)/(1.+php/p1p)+s12*(pow(conj(phperp),2)/(pow(abs(phperp),2)+p1p*phm)-pow(conj(p3perp)+(1.+php/p1p)*conj(p1perp),2)/((1.+php/p1p)*(pH.m2()+2.*p1.dot(pH))))); temp=temp*conj(temp); return temp.real(); } double C2qHqm(CLHEP::HepLorentzVector p2, CLHEP::HepLorentzVector p1, CLHEP::HepLorentzVector pH) { static double A=1./(3.*M_PI*v); // Implements Eq. (4.22) in hep-ph/0301013 double s12,p2p,p1p; COM p1perp,p3perp,phperp; // Determine first whether this is the case p1p\sim php>>p3p og the opposite s12=p1.invariantMass2(-p2); if (p2.pz()>0.) { // case considered in hep-ph/0301013 p2p=p2.plus(); p1p=p1.plus(); } else { // opposite case p2p=p2.minus(); p1p=p1.minus(); } p1perp=p1.px()+COM(0,1)*p1.py(); phperp=pH.px()+COM(0,1)*pH.py(); p3perp=-(p1perp+phperp); COM temp=A/(2.*s12)*(sqrt(p2p/p1p)*p3perp*conj(p1perp)+sqrt(p1p/p2p)*p1perp*conj(p3perp)); temp=temp*conj(temp); return temp.real(); } diff --git a/t/ME_data/ME_h_mt.dat b/t/ME_data/ME_h_mt.dat index 349d686..4fe4ab4 100644 --- a/t/ME_data/ME_h_mt.dat +++ b/t/ME_data/ME_h_mt.dat @@ -1,807 +1,807 @@ 3.999227056e-08 4.193300057e-05 3.89321526e-05 0.0168082963 4.057962883e-06 4.063245819e-07 6.411897603e-07 7.822040516e-07 6.064514954e-07 -2.338845036e-05 +2.323729555e-05 4.738208766e-05 1.624045121e-06 -0.0005243567971 +0.0005240707418 1.536443326e-05 9.5908372e-06 0.0339113482 8.152325116e-07 608.615459 5.962078524e-06 -0.01352350209 +0.01351013722 1.931630897 0.0005720381031 0.0008001259943 3.300880891e-05 -0.0007712160505 +0.0007251310095 0.0007037298905 0.03792377651 -1.079417185e-06 +1.031859624e-06 2.011798913 8.842231663e-08 0.0463578742 0.1064067277 9.357584012e-07 -0.007615558468 +0.007493578139 0.02272391876 3.606711294e-05 6.725111513e-07 0.005897687516 7.581255699e-05 0.005354681787 0.0006088947447 0.0083356406 0.003352568463 -3.95296533e-05 +3.873123108e-05 0.0001515363281 0.005531947882 -0.0388505972 +0.0388405017 2.312735565e-07 5.479259483e-05 8.480458133 -0.0001863702054 -7.349046553e-06 +0.0001862666068 +7.083344246e-06 5.091269317e-06 -0.0001452556797 -0.2136081233 +0.0001436528651 +0.2111893651 0.0238989834 4.446354691 -0.0002720132705 +0.0002706534845 1.267033388 0.005050109129 -2.848397879e-06 -0.07116846744 +2.768647989e-06 +0.06976587113 0.00280707581 -7.10418958 +7.013718037 7.216599144e-08 4.779195357e-06 -1.924410023e-05 +1.89454273e-05 1.522656761e-05 1.116772322e-06 0.0009539322219 9.279967349e-07 2.28564472 0.001331277201 0.001585734932 26972.25526 -1.2516302e-06 +1.207626582e-06 0.0006894623725 1.065245235e-07 0.01968534835 1.951627129e-06 4.401131811e-05 0.02561966521 49705.99132 5.383396351e-06 1.07104069e-06 2.808652999e-06 0.0001066898908 0.0002239236027 7.851273644e-08 0.002704063747 -0.0005842786082 +0.0005822864495 0.07550986067 -0.001335627914 -6.766626128e-05 +0.001330057569 +6.413753514e-05 0.1149401044 0.183403668 0.000113873958 -1.844065382e-05 -0.02689433062 +1.831402184e-05 +0.02660643026 0.001076022562 5.066109523e-06 -0.1097193795 +0.10939624 9.253312815e-05 -92.12237679 +92.02148669 0.01945125628 8.958261209e-05 -3.849412514e-06 -0.0003305165961 +3.608146549e-06 +0.0003265014563 0.0001299537598 1.170328977e-06 0.0001328201344 3034.923139 1.378838481 -0.08677679017 -36.20174579 +0.08663556522 +36.19909049 1.785612005e-06 -1.032516793e-05 +1.018401168e-05 0.001615788007 1.23208212 0.0003642605012 2.708222463e-06 4.719890307e-06 0.001710230342 0.3182780512 -0.0009491003377 -0.09165888221 +0.0009417630586 +0.09122795432 0.2629019547 0.02748771871 -3.179028152e-06 -0.002811209731 +3.029545829e-06 +0.002797757699 1.89343446e-06 1.77338727e-06 2.291261955e-05 1.31623411e-05 0.02865008699 -0.0001547359261 +0.0001521683687 5.992647136 427898745.7 0.0002908564713 180.1583563 1.149638817 -4.915348327e-06 +4.874677273e-06 0.02938178265 12.41064482 0.0005867842208 1.235864829e-06 -6.021062978e-05 +5.86484668e-05 0.7228122449 -0.0003703954772 +0.0003638994119 1.100546805e-07 0.0007569561799 3.456980613e-06 0.001902875586 1.766220045 0.0006636538355 435.7550934 1.468265667e-05 -0.003895447939 +0.003892402317 0.005375424637 9.895862027e-06 0.002634610055 0.04246482577 1.811370936e-05 0.0008533026641 -1.973172924 -2.216040365e-06 +1.969024996 +2.188193623e-06 0.0001940609469 6.432993714e-05 436.2325667 8.585463564e-05 0.003171908449 3.1453383e-07 0.10708771 1.134414955e-07 5.434578249e-06 6.681002905e-06 -0.002957334168 -0.005022320071 +0.002902943775 +0.004943649156 2.421882039e-05 3.783451257e-08 5.144794528e-06 -0.0001067748459 +0.000104411672 13660.98467 117.1131418 387.5973894 -0.2286644594 +0.2272332194 2.335478847e-06 0.0001022288321 0.004686932429 0.4278874185 -0.0001089502807 +0.0001086686824 1.985873409e-06 2.25291051e-07 2.945623146e-05 0.3581923626 -8.204917967e-06 +8.047235434e-06 0.004893941209 0.0005325322641 1.060152245e-07 1.118231806e-05 2.203633348e-06 -4.953456305e-06 +4.844095806e-06 0.007400958058 5.769269726e-05 1.899358848e-06 0.002291707993 0.007504512648 0.0946266303 -0.07677680285 -0.001406412108 +0.0767367567 +0.001382893492 0.000131740617 3.126413693e-07 6.676941465e-07 2.247981499e-07 1.394670123e-07 0.0003117083998 0.003124792631 3.196614452e-07 0.0005145950558 1.196095567 1.715539058e-05 2.527890627e-05 -4.198681685e-06 +4.028317881e-06 4.483788358e-07 3.597134253e-05 6.514469651e-07 0.002675227421 0.1700338633 0.0151159864 -2.549964327e-05 -0.00103462801 -0.0002092909388 +2.475242006e-05 +0.001032856771 +0.000205696546 3.068873861e-06 1176.630703 9.378378651e-06 7.10431401e-07 4.541139018e-05 3.954117424e-06 0.003271203855 0.01288628186 0.01758893926 0.0001040458517 0.08587993013 7.493297002e-07 0.0303353991 2.525780749e-07 1.075360445e-06 0.0005909996117 3.980399911e-06 67.71892274 -8.935469327e-06 +8.920686332e-06 3.997018327e-05 -3.48516788e-06 +3.396851165e-06 0.1177580565 1.362864165e-07 0.00719011227 3.379854611e-06 1.335802906e-08 0.0008096413047 8.017608649e-08 5.475646202e-07 9.89836826e-07 0.0001097524935 0.005006547178 0.01328382633 3.708077203e-10 1.7738466e-09 6.195338978e-10 1.375292134e-08 -1.198886719e-05 +1.151374303e-05 1.087672399e-10 7.612820481e-06 1.631362221e-06 7.195255533e-08 8.307139358e-09 2.261143452e-08 -3.176205908e-07 +3.075268256e-07 7.680504754e-08 1.175875499e-07 0.0006018013011 1.467361257e-05 1.350450721e-07 -4.010394331e-07 -7.773787384e-07 +3.850336065e-07 +7.767624725e-07 0.0170653268 0.0001792376906 3.072166197e-06 22266.55115 0.06729772649 -1.156967792 +1.14678661 7.103072447e-07 0.0008701526857 0.3457109502 0.001639388444 6.829709377e-07 103.0220399 1.724918219e-06 0.006442680294 0.0002924073612 0.0061407479 3.014208435e-07 -1.067129219e-05 +1.066948362e-05 5.819752501e-07 0.002933183567 1.390851143e-08 1.359009867e-06 0.02819012119 1.30006471e-09 78.00480652 3.27001628e-06 2.84175006e-07 -4.457488422e-07 +4.405951989e-07 0.0003605232041 1.871558115e-06 2.854501301e-07 0.02600693272 0.0004600409654 -0.001476738524 +0.001462458297 2.488693144e-08 2.395438283e-11 4.248528123e-07 1.389769895e-05 8.049135616e-07 2.12824677e-07 8.405933433e-07 2.36621812e-06 0.06259159217 1.831343487e-07 2.1663617e-08 6.323481085e-10 1.339773875e-07 4.883049792e-06 1.705369167e-05 0.0001013417062 0.01536206656 1.67093181 5.039245627e-10 1.141488157e-05 3.188719492e-07 3.24135522e-10 0.0002434524475 1.383868594e-06 0.498790572 9.03612161e-05 0.02764416817 0.01341805272 7.529041731e-07 0.0002257161473 0.0003717671499 1.858022341e-08 5.715357136e-07 2.477427879e-09 0.0001517676839 2.326424939e-11 1.729565889e-05 7.354412281e-09 -9.505701399e-06 +9.213659243e-06 0.0005289104709 0.02479967619 6.212490017e-06 1.714682348e-08 7.916835703e-09 0.00934088881 4.083237353e-07 0.004696390685 6.761095529e-10 0.003004586166 6.030026635e-05 8.803060787e-08 6.957661086e-08 9.612925872e-08 1.628346787e-05 5.143053469e-07 8.784807581e-08 1.024493405e-08 3.137825282e-07 42.00569276 0.004640201019 1.146513701e-05 -5.876691844e-08 +5.674508354e-08 3.381428808e-06 2.603443408e-07 -0.002350385674 +0.0023417509 1.276808422e-05 0.06226820245 1.778585408e-08 2.378758757e-07 9.7714181e-05 0.387678215 4.398556239e-07 3.0060972e-08 3.874778496e-05 -1.0694862e-05 +1.051980557e-05 1.402325829e-06 4.695875643e-06 2.383017224e-07 2.417306793e-06 0.0002708928034 0.01384277314 1.710909166e-07 1.217157554 177.0639344 0.03861638077 7.433640647e-05 0.0002263704383 0.002948183736 0.008773163623 5.353914571e-08 360.0073872 6.478760385e-07 2.750359008e-07 8.589493477e-08 4.263303371e-07 2.005861581e-05 5.503740271e-08 0.00180575553 16.03887797 7.922076143 5.036869078e-12 1.277093114e-06 0.003373959327 1.743675895e-05 1.778007998e-06 0.0002212911594 9.259789313e-05 8.539690283e-05 2.423945599e-06 0.0001233083611 -9.174270377e-09 +8.923119994e-09 2.570708518e-07 0.0008554774811 -0.3056275035 +0.3048484598 0.0006441353055 1.203219038e-07 4.75874322e-07 258.1003097 2.768523645e-09 1.270156521e-11 1.46097374e-09 34.87523283 8.506779165e-05 4.671813569e-06 1.907838682e-05 3.716521365e-07 1.888947852e-08 -0.02682389012 +0.02680923327 15.26083816 2.873162527e-10 1.360352114e-06 9.527755241e-07 0.0003084354836 0.0009838761157 1.302357534e-05 0.004623681617 0.01694286548 9.642283586e-07 8.499293031e-06 8.770386824e-06 0.001276821989 6.668543971e-07 5.652013261e-11 0.238779506 1.624068992e-07 0.0002438443448 0.002313588467 7.53704251e-09 7.248545824e-12 0.02160343338 6.465648118e-09 7.430463219e-08 0.000611761207 2.262200983e-07 8.810901655e-07 1.688360234e-07 -3.452136087e-06 +3.392772149e-06 5.480913194e-08 0.0004601409972 3.364979337e-10 8.967694023e-09 0.1345541103 2.492529354e-06 9.220762636e-06 7.386130159e-06 3.598519539 8.045779179e-06 4.730716431e-07 2.192459188e-07 2.194928052e-06 3.130712514e-07 1.43402212e-05 0.0002608840076 23.2295584 -1.820768009e-08 +1.778524358e-08 4.161369951e-06 0.001762466634 4.255417524e-09 3320585.541 0.0009563172933 0.03635196251 0.5161165209 0.004796521016 2.31018301e-06 0.06003700358 0.0004342970655 2.009638611e-08 0.0001808310494 0.0003251077547 -6.445298688e-07 +6.20919327e-07 2.882854303e-08 0.0422146142 8.524423251e-07 19.01825412 0.000339795268 1.575704657e-07 0.0004837774983 4.132343892e-07 5.931203648e-09 0.01135109493 5.877226405e-05 1.480445352e-08 0.3379228367 1.786471685e-06 1.021098483e-07 1.079160999e-08 -1.26575267e-08 +1.232033094e-08 0.002978083967 4.588139504e-06 6.6721226e-11 3.261090243e-06 0.06059954782 4.847476637e-09 1.104977758e-06 -0.0003977585716 -0.0006607921183 +0.0003884699782 +0.0006601268859 0.000225642662 4.294099642e-06 1.185934783e-07 0.001042912926 8.03811263e-06 1.755333168e-10 6.77273477e-08 598.9038787 1.572919365e-07 0.04266693618 2.739698989e-08 0.0005135423693 1.264426626e-10 2.608778512e-07 1.857091679e-06 1.075074651e-10 1.188966384e-07 0.03379476153 0.0006473640268 3.624390296e-09 1.237883305e-05 1.016117525e-05 9.687850834e-06 0.0002738964397 9.038889191e-09 7.786947927e-07 3.198680693e-05 1.608278345e-11 0.003353487144 6.107229663e-08 0.0003768704931 0.0002005945404 3.670401818e-11 9.030155872e-09 2.042723216e-08 2.386955262e-08 4.835084799e-09 3.685817471e-11 1.066443786e-07 2.652667148e-07 3.649486136e-11 0.02477462094 3.517378902e-05 3.080027841e-09 1.018709237e-09 2.870228082e-10 8.744977958e-06 5.273214292e-09 3.168385461e-11 0.001556612382 -7.188695072e-07 +6.736179068e-07 2.492191426e-09 7.615392241e-10 0.01004217426 1.101368692e-06 3.918013244e-06 3.713827343e-09 6.857767402e-09 1.290966146e-09 6.224822208e-12 3.376319521e-06 -4.097980204e-09 +4.080482631e-09 7.633357106e-05 6.164135969e-05 1.630673977e-06 7.08423096e-09 4.080955534e-10 232770.6483 0.008068531153 8.666549381e-10 2.246457597e-10 3.222301462e-08 1.395629818e-06 2.196466632e-08 1.592802071e-08 4.326999161e-05 0.0005042839931 0.00197308444 3.032370137e-05 0.0005857579152 7.171987925e-06 6.864328965e-06 4.00109422e-06 0.005161817275 2.375668536e-09 9.778482285e-06 1.181506681e-08 1.503749258e-08 3.203838424e-11 0.3676119133 1.281881299e-09 0.0004415836065 -1.205682724e-07 +1.203972243e-07 8.77428226e-07 7.210177682e-11 1.45886715e-06 0.02121994342 0.01168544154 3.982811877e-09 8.600471911e-09 1.434258247e-08 8.098070419e-10 0.0007724068717 2.490831815e-12 2.21358682e-06 -0.3039319263 -1.06879986e-05 +0.302087331 +1.012854812e-05 2.722304477e-09 3.122833335e-06 1.198571403e-06 4.489385847e-09 2.67381748e-08 3.471377557e-09 5.756644016e-06 6.839839041e-05 2.209846587e-07 8.291924237e-08 5.242886271e-07 6.844201737e-11 1.726048304e-06 9.225948144e-07 0.000185717301 0.002602315394 6.138200606e-11 0.0001818846094 7.733889705e-14 1.358383639e-08 1.07000594e-06 4.579302999e-07 7.849169727e-07 2.448333585 7.535064872e-11 3.779457663e-06 1.495153092e-07 3.954406596e-09 8.122956595e-09 0.0001228391098 2.269361403e-11 0.001844497418 0.01115645921 -2.597690032e-06 +2.540174677e-06 4.660852875e-11 4.199474743e-10 4.628821128e-05 0.0004980806998 7.598116953e-09 4.143456828e-13 0.002644658829 2.564668699e-08 6.50580785e-08 3.820040563e-11 0.2398554094 9.400695396e-08 0.002358652784 1.171877918e-11 1.637975849e-09 8.942959658e-08 3.408065038e-10 6.371021504e-07 8.862071381e-07 8.07800753e-13 1.275141633e-07 5.065134208e-08 3.723768129e-05 1.586955855e-06 4.630822765e-08 3.869149299e-05 0.0002373420357 8.091381421e-08 2.031825403e-10 1.493296083e-05 2.814058365e-06 4.293942142e-10 1.227720708e-09 9.668701181e-10 5.082776865e-06 0.001170366525 -3.153572708e-07 +3.13097726e-07 2.929772297e-06 3.908673483e-06 -1.837314879e-08 +1.815611522e-08 9.22544372e-06 0.004601199051 1.245696104e-11 1.623356072e-06 0.0004075079271 0.06223654452 2.229639676e-12 1.548827536e-10 0.2098995653 0.02739866019 1.922659779e-12 0.001586041269 1.845383825 0.8609993903 0.09603688052 5.904130818e-11 9.958673683e-12 0.05389921105 0.0001869127049 6.618563858e-08 3.845467031e-07 9.065990802e-07 0.0003895671306 2.782863325e-11 0.0002689153495 5.924260979e-06 0.0005747447785 4.117458599e-09 5.117012888e-08 3.396356612e-09 8.109913215e-12 4.329382871e-08 4.881990436e-07 1.276280038e-11 9.120675139e-12 7.192434752e-07 3.913876447e-11 1.9517939e-08 0.11670582 1.114874302e-06 -1.051244942e-06 +1.024310838e-06 2.596368516e-07 0.0001277437573 0.0004507490086 5.934514686e-10 9.045316154e-06 1.1680623e-08 4.405302098e-07 1.103299007e-15 7.105364025e-08 6.778113877e-05 2.587602347e-06 2.630707294e-07 2.256922483e-10 1.585065727e-05 -0.0001137727537 +0.0001112190986 1.439932327e-09 0.003315113132 4.417287152e-05 3.208255435e-09 7.222413113e-10 0.000300699488 4.624336745e-05 8.037199117e-10 6.718759329e-12 3.83547078e-05 7.94104255e-05 4.496517789e-05 1.351102043e-06 9.177810983e-10 8.693619098e-07 1.929511565e-06 4.773961763e-09 1.896551164e-08 1.552346596e-08 9.245502051e-06 0.0001037407773 6.368845512e-07 4.041904863e-06 1.506851808e-08 9.717844454e-12 2.297581007e-05 0.0234227485 0.00306449787 5.437659639e-08 0.001404494601 1.582456182e-12 10.74537802 0.0006876957837 3.080472252e-07 4.979458387e-09 2.191738006e-08 diff --git a/t/ME_data/ME_h_mtmb.dat b/t/ME_data/ME_h_mtmb.dat index 8370f81..0a6185c 100644 --- a/t/ME_data/ME_h_mtmb.dat +++ b/t/ME_data/ME_h_mtmb.dat @@ -1,807 +1,807 @@ 4.043634525e-08 3.947735809e-05 3.978333671e-05 0.01707962353 4.133534858e-06 4.12758231e-07 6.527379321e-07 7.807924441e-07 6.213247187e-07 -2.292064051e-05 +2.277250905e-05 4.929604661e-05 1.648537835e-06 -0.0005321070091 +0.0005318167258 1.551221705e-05 9.732109994e-06 0.03508523471 8.385697667e-07 608.0814044 6.0749939e-06 -0.01383234108 +0.01381867099 1.96358809 0.0005777778828 0.000827073758 3.396653466e-05 -0.0007925662403 +0.0007452053903 0.000712177849 0.03945983226 -1.109333656e-06 +1.060458019e-06 2.002990871 9.107975513e-08 0.04804489153 0.1078671772 9.503786206e-07 -0.007776569624 +0.007652010339 0.02341176253 3.684250185e-05 6.861318278e-07 0.005400516499 7.647108683e-05 0.005354950581 0.0006200616798 0.008448060377 0.003448403835 -3.93020183e-05 +3.850819386e-05 0.0001546133101 0.005591549046 -0.03958941813 +0.03957913065 2.350530066e-07 5.66847762e-05 8.602097332 -0.0001921261891 -7.502702823e-06 +0.000192019391 +7.231445125e-06 5.119802912e-06 -0.0001496819684 -0.213120324 +0.0001480303122 +0.2107070894 0.02428150541 4.539435975 -0.0002738016111 +0.0002724328853 1.314349378 0.005127584311 -2.816422817e-06 -0.07344901515 +2.73756817e-06 +0.07200147354 0.002877526755 -7.299557828 +7.206598278 7.320533172e-08 4.686272323e-06 -1.946859515e-05 +1.916643801e-05 1.561186199e-05 1.127468541e-06 0.0009956223126 9.322344484e-07 2.355779109 0.001354263217 0.001607437479 27300.17114 -1.282038387e-06 +1.236965707e-06 0.0006972261753 1.071457504e-07 0.02005686587 1.961498186e-06 4.460773766e-05 0.02601668696 50724.57307 5.480882962e-06 1.076584976e-06 2.869362531e-06 0.0001072518549 0.0002293795185 8.022114184e-08 0.002698012227 -0.0006015504033 +0.0005994993545 0.07725070316 -0.001377279885 -6.964246376e-05 +0.001371535827 +6.601068069e-05 0.1179283377 0.187466976 0.0001155356229 -1.895976068e-05 -0.02724871881 +1.882956399e-05 +0.02695702478 0.001087874101 5.154078743e-06 -0.1059710091 +0.1056589091 9.484855944e-05 -94.78499143 +94.6811853 0.02033391802 9.092704303e-05 -3.912572591e-06 -0.0003340065858 +3.667348002e-06 +0.0003299490494 0.0001322670062 1.193514887e-06 0.0001344144727 3095.568362 1.406347883 -0.08943984916 -35.93630323 +0.08929429021 +35.9336674 1.82576512e-06 -1.057243426e-05 +1.042789761e-05 0.001666287395 1.269161514 0.0003682229129 2.790042514e-06 4.841321328e-06 0.001729869633 0.3293946648 -0.0009589102693 -0.08469469459 +0.0009514971519 +0.08429650834 0.2664463909 0.02792004132 -3.241512757e-06 -0.002892655486 +3.089092321e-06 +0.002878813725 1.908915912e-06 1.7932868e-06 2.342012655e-05 1.3539045e-05 0.02918304927 -0.0001576463616 +0.0001550305109 6.099926061 438766540.3 0.0003006556126 182.0743425 1.166556788 -5.062168162e-06 +5.020282277e-06 0.02991092642 12.79082303 0.0005990054968 1.245221356e-06 -6.010815261e-05 +5.854864839e-05 0.729519575 -0.0003654608813 +0.0003590513598 1.113634813e-07 0.0007837766949 3.49282179e-06 0.001938345149 1.741954359 0.000679566605 447.8606334 1.502685811e-05 -0.004001081006 +0.003997952796 0.005572294955 1.00407234e-05 0.002697149585 0.04268640856 1.838639564e-05 0.0008683709605 -2.022642541 -2.281686632e-06 +2.018390621 +2.25301498e-06 0.000198987372 6.590604109e-05 451.0223058 8.870879716e-05 0.003243511862 3.187920803e-07 0.1107987679 1.145759171e-07 5.532565231e-06 6.903234805e-06 -0.002727714074 -0.005134520243 +0.002677546784 +0.005054091796 2.467146722e-05 3.863236208e-08 5.225855864e-06 -0.0001075855321 +0.0001052044158 13702.18793 118.8790931 398.6816605 -0.2300421985 +0.2286023351 2.355356756e-06 0.0001050726086 0.004854583584 0.4416609149 -0.0001117424705 +0.0001114536554 2.015137689e-06 2.270813172e-07 3.065037179e-05 0.3695610208 -8.012249806e-06 +7.85826998e-06 0.005034520677 0.0005496323042 1.077013565e-07 1.134883954e-05 2.217217266e-06 -5.058740304e-06 +4.947055385e-06 0.007630107217 5.870346284e-05 1.937770986e-06 0.002308560802 0.00756714829 0.09713373161 -0.07902736629 -0.001373738971 +0.07898614627 +0.001350766728 0.0001360267898 3.146930172e-07 6.721239549e-07 2.308081382e-07 1.402445764e-07 0.000317796034 0.003250721897 3.220983676e-07 0.0005361385069 1.1674911 1.701754998e-05 2.555659626e-05 -4.295773543e-06 +4.121470184e-06 4.575225674e-07 3.723503521e-05 6.625793946e-07 0.002713331741 0.1757705927 0.01528287892 -2.600058607e-05 -0.00104396898 -0.000215183149 +2.52386836e-05 +0.001042181749 +0.0002114875625 3.169381656e-06 1197.183436 9.459688271e-06 7.151860571e-07 4.709453403e-05 4.024117529e-06 0.003399032872 0.01326938543 0.0181234889 0.0001070400499 0.08816817298 7.557106621e-07 0.03150507468 2.546159018e-07 1.083218854e-06 0.0006029839411 4.109536026e-06 69.62539127 -9.190394104e-06 +9.175189357e-06 3.976334138e-05 -3.557434726e-06 +3.467286716e-06 0.1217200425 1.39072939e-07 0.007475099982 3.40898863e-06 1.355970335e-08 0.0008211141398 8.085185969e-08 5.540708628e-07 9.793920273e-07 0.0001116678867 0.005087064408 0.01371958748 3.747341913e-10 1.784590611e-09 6.285401932e-10 1.403711889e-08 -1.234662607e-05 +1.185732376e-05 1.099705834e-10 7.694717507e-06 1.645454731e-06 7.281682073e-08 8.437713425e-09 2.293298229e-08 -3.274272811e-07 +3.170218661e-07 7.748072516e-08 1.184163566e-07 0.0006048814428 1.496340286e-05 1.358047852e-07 -4.122637606e-07 -7.951549138e-07 +3.958099615e-07 +7.945245559e-07 0.01733476065 0.0001847336643 3.09120607e-06 22646.44954 0.06825326182 -1.192340601 +1.181848143 7.194716543e-07 0.0008831460117 0.3510995791 0.001703720233 6.923341702e-07 104.5021475 1.817413974e-06 0.006274458892 0.000298133205 0.006212019474 2.971106209e-07 -1.075176214e-05 +1.074993994e-05 6.001267547e-07 0.002926676713 1.394485679e-08 1.378850246e-06 0.02873430354 1.324231029e-09 78.68865757 3.323528263e-06 2.901477846e-07 -4.493480222e-07 +4.441527661e-07 0.0003762267637 1.888509349e-06 2.984779601e-07 0.02649146415 0.0004735978501 -0.001523297471 +0.001508567014 2.511800946e-08 2.406174629e-11 4.343797764e-07 1.413682848e-05 8.187236684e-07 2.137756399e-07 8.487602643e-07 2.382894999e-06 0.06385257472 1.843038529e-07 2.200340736e-08 6.500296312e-10 1.369543749e-07 5.034654856e-06 1.729239287e-05 0.0001029880478 0.01584013323 1.717937137 5.051793703e-10 1.169719308e-05 3.274892339e-07 3.257575264e-10 0.0002494213231 1.395452311e-06 0.503263349 9.183628731e-05 0.02652180806 0.01351083505 7.697372696e-07 0.0002345916253 0.0003785854477 1.911331681e-08 5.741981368e-07 2.496303102e-09 0.0001538309849 2.342705363e-11 1.733704579e-05 7.440708121e-09 -9.634074462e-06 +9.338088321e-06 0.0005366276141 0.02581853732 6.351131637e-06 1.725302693e-08 7.963798782e-09 0.009484858796 4.13192111e-07 0.004760868685 6.868174526e-10 0.003035550995 6.122424548e-05 8.912863951e-08 7.049197773e-08 9.784536423e-08 1.644442743e-05 5.17689779e-07 8.945719428e-08 1.043648257e-08 3.220171418e-07 42.8640364 0.004721640337 1.17318728e-05 -6.030229416e-08 +5.822763572e-08 3.442022366e-06 2.613636766e-07 -0.002417874611 +0.002408991899 1.311028403e-05 0.06454239294 1.788371853e-08 2.431129133e-07 9.887586539e-05 0.393826783 4.435492073e-07 3.017594495e-08 3.944030937e-05 -1.049807075e-05 +1.032623546e-05 1.480534305e-06 4.80208438e-06 2.278546025e-07 2.431285271e-06 0.0002795248679 0.01414217835 1.721740596e-07 1.225137873 184.4942235 0.03969271858 7.5889168e-05 0.0002307099487 0.002991500607 0.008823843341 5.377454625e-08 363.5972586 6.577363411e-07 2.768216121e-07 8.808374953e-08 4.320144869e-07 2.087813482e-05 5.563234756e-08 0.001828101352 16.31099206 8.047143214 5.054980175e-12 1.31165661e-06 0.003453525346 1.775158542e-05 1.735750599e-06 0.0002258034827 9.302350911e-05 8.722704625e-05 2.463340497e-06 0.0001266405389 -9.536618271e-09 +9.275548428e-09 2.616936561e-07 0.0008671500774 -0.3105556664 +0.3097640608 0.0006552874753 1.250994614e-07 4.836707652e-07 267.0859855 2.784590138e-09 1.277933942e-11 1.625022487e-09 35.42273613 8.653074513e-05 4.787553737e-06 1.988624796e-05 3.776413822e-07 1.896902599e-08 -0.02770700524 +0.02769186584 15.7257863 2.903221186e-10 1.383656135e-06 9.796122701e-07 0.000313466982 0.001008582166 1.326553759e-05 0.004678660451 0.01717544418 9.763673978e-07 8.769786833e-06 8.887711201e-06 0.001304503338 6.699929944e-07 5.668340919e-11 0.2432408562 1.65579253e-07 0.000246842309 0.002380529019 7.590644517e-09 7.382339234e-12 0.02191318318 6.595093824e-09 7.59955463e-08 0.0006241836255 2.306006909e-07 8.949234518e-07 1.745865901e-07 -3.531291359e-06 +3.470566244e-06 5.515085903e-08 0.0004672383281 3.41289611e-10 9.072028681e-09 0.1360437848 2.516296619e-06 9.347242718e-06 7.560996071e-06 3.650323724 8.084973509e-06 4.821049366e-07 2.262423437e-07 2.237305224e-06 3.165244113e-07 1.440989373e-05 0.0002661269589 24.06039199 -1.860370743e-08 +1.817208269e-08 4.189979743e-06 0.001821016076 4.36609224e-09 3372180.49 0.0009765702499 0.03671795044 0.5266849257 0.004945279364 2.349850855e-06 0.06096774779 0.000440582517 2.030747703e-08 0.0001828061992 0.0003290502891 -6.504778186e-07 +6.266493903e-07 2.945465959e-08 0.04314824277 8.695492983e-07 19.3255897 0.0003433917438 1.623604604e-07 0.0004985403516 4.206613628e-07 6.048976857e-09 0.01167821222 5.984173029e-05 1.495919742e-08 0.349858651 1.798434603e-06 1.039820674e-07 1.118772128e-08 -1.30086451e-08 +1.266209558e-08 0.003069395217 4.6681898e-06 6.697088362e-11 3.350850437e-06 0.06185962328 4.883448988e-09 1.110833085e-06 -0.0004052046295 -0.0006808854079 +0.0003957421533 +0.0006801999472 0.0002346905421 4.330939823e-06 1.219218909e-07 0.00107853716 8.281691391e-06 1.767444691e-10 6.900019008e-08 605.6626124 1.58440332e-07 0.04329139184 2.752138714e-08 0.0005185712214 1.286238525e-10 2.639855477e-07 1.869433683e-06 1.079551089e-10 1.229225358e-07 0.03468655392 0.0006584274735 3.689112319e-09 1.260350348e-05 1.040877611e-05 9.878538403e-06 0.00027695052 9.100037836e-09 8.033944341e-07 3.315013452e-05 1.633327055e-11 0.003414227642 6.245172231e-08 0.0003805246583 0.0002060494734 3.729836801e-11 9.13272934e-09 2.051216421e-08 2.436679394e-08 4.8855578e-09 3.735846979e-11 1.088675628e-07 2.683031236e-07 3.706588258e-11 0.02523141233 3.549809347e-05 3.149827076e-09 1.040176277e-09 2.910870047e-10 8.86122235e-06 5.308664537e-09 3.218598972e-11 0.001577139281 -7.379705718e-07 +6.915165923e-07 2.55834443e-09 7.691641781e-10 0.01021181751 1.105572103e-06 3.947538433e-06 3.734235002e-09 6.899637353e-09 1.304337436e-09 6.538762872e-12 3.449997548e-06 -4.212613391e-09 +4.194626357e-09 7.962318534e-05 6.295652823e-05 1.652753715e-06 7.186541052e-09 4.161402504e-10 240986.4156 0.00837085449 8.760965452e-10 2.273546864e-10 3.315418794e-08 1.416471539e-06 2.236464809e-08 1.633966009e-08 4.518163364e-05 0.0005095243003 0.001992318371 3.067811242e-05 0.000600451339 7.236523643e-06 6.932191963e-06 4.091606026e-06 0.00532740987 2.395091742e-09 9.944229055e-06 1.205282983e-08 1.590006487e-08 3.215700583e-11 0.3737229127 1.296635624e-09 0.0004534080133 -1.233610704e-07 +1.231860602e-07 8.898892246e-07 7.255073266e-11 1.472337758e-06 0.02150215784 0.01192727543 4.01006022e-09 8.731812041e-09 1.44194837e-08 8.229024075e-10 0.0007855372767 2.506773438e-12 2.241537936e-06 -0.3085183925 -1.101920302e-05 +0.3066459614 +1.044241605e-05 2.752322662e-09 3.158084387e-06 1.215204067e-06 4.56919432e-09 2.790752035e-08 3.511291959e-09 5.863328277e-06 6.94327717e-05 2.22257262e-07 8.350918192e-08 5.331043956e-07 6.969438883e-11 1.760059719e-06 9.217915003e-07 0.000186591209 0.002628593815 6.211758257e-11 0.0001834032575 7.768284292e-14 1.377546165e-08 1.091734431e-06 4.677996963e-07 7.943499964e-07 2.545817678 7.562150253e-11 3.799329956e-06 1.510534633e-07 3.983823634e-09 8.184966815e-09 0.000124713148 2.309165194e-11 0.001872087514 0.01153490791 -2.680798294e-06 +2.621442843e-06 4.693787525e-11 4.26518767e-10 4.684595348e-05 0.0005076579094 7.644830233e-09 4.193877203e-13 0.002676923954 2.654050699e-08 6.800195443e-08 3.888064789e-11 0.2446884347 9.628489292e-08 0.002537052711 1.189694289e-11 1.656137485e-09 9.040798991e-08 3.488787289e-10 6.557255554e-07 9.083676565e-07 8.150458329e-13 1.315321346e-07 5.098465051e-08 3.786448639e-05 1.600711068e-06 4.658140964e-08 3.906381479e-05 0.0002389891749 8.287668941e-08 2.051670765e-10 1.60134009e-05 2.847282071e-06 4.308493985e-10 1.238922307e-09 9.741854048e-10 5.146899686e-06 0.001209019717 -3.233996009e-07 +3.210824326e-07 2.974261102e-06 3.899610282e-06 -1.896885723e-08 +1.874478682e-08 9.332118489e-06 0.004694117038 1.266652132e-11 1.65446818e-06 0.0004100868952 0.06302457519 2.236229887e-12 1.584413907e-10 0.2179964362 0.0278380659 1.936978481e-12 0.001609711018 1.892891439 0.8712003892 0.09866862844 5.934210083e-11 1.002311696e-11 0.05494208606 0.000190347461 6.722151007e-08 4.030694324e-07 9.278979899e-07 0.0003978128465 2.803383506e-11 0.0002738201279 5.960007001e-06 0.0005814501443 4.157124101e-09 5.198894043e-08 3.465884394e-09 8.167365756e-12 4.450019241e-08 4.893968894e-07 1.288321375e-11 9.163290302e-12 7.451831479e-07 3.940383302e-11 1.996508632e-08 0.1236118645 1.182043812e-06 -1.081851533e-06 +1.054133253e-06 2.618934064e-07 0.0001307650885 0.0004571929478 6.07226713e-10 9.141057675e-06 1.18356058e-08 4.465375632e-07 1.107633947e-15 7.416788486e-08 6.870175273e-05 2.600618361e-06 2.647425108e-07 2.279494937e-10 1.611105573e-05 -0.0001167395509 +0.0001141193053 1.474187502e-09 0.003419123567 4.527389745e-05 3.271142207e-09 7.434473088e-10 0.0003102710133 4.710559791e-05 8.139484815e-10 6.749221433e-12 3.882836934e-05 8.016924727e-05 4.578651435e-05 1.374899461e-06 9.31258077e-10 9.107579154e-07 2.005597861e-06 4.853545363e-09 1.906194275e-08 1.569062034e-08 9.37420375e-06 0.0001051157618 6.464238034e-07 4.092696924e-06 1.538436694e-08 9.760057524e-12 2.626691189e-05 0.02379715614 0.00314703493 5.484543738e-08 0.00141690332 1.612112513e-12 11.09088062 0.0007046654234 3.111374555e-07 5.038033967e-09 2.21156295e-08