#include "PHOTONS++/Main/Photons.H" #include "ATOOLS/Org/Message.H" #include "ATOOLS/Org/Run_Parameter.H" #include "ATOOLS/Org/My_Limits.H" #include "ATOOLS/Org/Scoped_Settings.H" #include "ATOOLS/Phys/Blob.H" #include "PHOTONS++/Main/Define_Dipole.H" #ifdef PHOTONS_DEBUG #include "ATOOLS/Math/Histogram_2D.H" #include "ATOOLS/Org/Shell_Tools.H" #endif using namespace PHOTONS; using namespace ATOOLS; using namespace std; std::ostream &PHOTONS::operator<<(std::ostream &str,const yfsmode::code &ym) { if (ym==yfsmode::off) return str<<"Off"; else if (ym==yfsmode::soft) return str<<"Soft"; else if (ym==yfsmode::full) return str<<"Full"; return str<<"unknown"; } std::istream &PHOTONS::operator>>(std::istream &str,yfsmode::code &ym) { std::string tag; str>>tag; ym=yfsmode::full; if (tag.find("None")!=std::string::npos) ym=yfsmode::off; else if (tag.find("Soft")!=std::string::npos) ym=yfsmode::soft; else if (tag.find("1")!=std::string::npos) ym=yfsmode::soft; else if (tag.find("Full")!=std::string::npos) ym=yfsmode::full; else if (tag.find("2")!=std::string::npos) ym=yfsmode::full; return str; } // define statics yfsmode::code PHOTONS::Photons::s_mode = yfsmode::full; bool PHOTONS::Photons::s_useme = true; double PHOTONS::Photons::s_ircutoff = 1E-3; double PHOTONS::Photons::s_uvcutoff = std::numeric_limits::max(); int PHOTONS::Photons::s_ircutoffframe = 0; double PHOTONS::Photons::s_accu = 1E-6; int PHOTONS::Photons::s_nmax = std::numeric_limits::max(); int PHOTONS::Photons::s_nmin = 0; double PHOTONS::Photons::s_drcut = 1000.; bool PHOTONS::Photons::s_strict = false; double PHOTONS::Photons::s_reducemaxenergy = 1.; double PHOTONS::Photons::s_increasemaxweight = 1.; bool PHOTONS::Photons::s_checkfirst = false; int PHOTONS::Photons::s_ffrecscheme = 0; int PHOTONS::Photons::s_firecscheme = 0; double PHOTONS::Photons::s_alpha = 0.; double PHOTONS::Photons::s_alpha_input = 0.; bool PHOTONS::Photons::s_userunningparameters = false; #ifdef PHOTONS_DEBUG std::string PHOTONS::Photons::s_histo_base_name("weights"); Histogram_2D PHOTONS::Photons::s_histo_dipole = Histogram_2D(101,1e-6,1e4,100,-0.5,10.5,11); Histogram_2D PHOTONS::Photons::s_histo_jacobianM = Histogram_2D(101,1e-6,1e4,100,-0.5,10.5,11); Histogram_2D PHOTONS::Photons::s_histo_jacobianL = Histogram_2D(101,1e-6,1e4,100,-0.5,10.5,11); Histogram_2D PHOTONS::Photons::s_histo_higher = Histogram_2D(101,1e-6,1e4,100,-0.5,10.5,11); Histogram_2D PHOTONS::Photons::s_histo_yfs = Histogram_2D(101,1e-6,1e4,100,-0.5,10.5,11); Histogram_2D PHOTONS::Photons::s_histo_total = Histogram_2D(101,1e-6,1e4,100,-0.5,10.5,11); Histogram_2D PHOTONS::Photons::s_histo_t_dipole = Histogram_2D(111,1e-6,1e4,100,1e-6,1e2,20); Histogram_2D PHOTONS::Photons::s_histo_t_jacobianM = Histogram_2D(111,1e-6,1e4,100,1e-6,1e2,20); Histogram_2D PHOTONS::Photons::s_histo_t_jacobianL = Histogram_2D(111,1e-6,1e4,100,1e-6,1e2,20); Histogram_2D PHOTONS::Photons::s_histo_t_higher = Histogram_2D(111,1e-6,1e4,100,1e-6,1e2,20); Histogram_2D PHOTONS::Photons::s_histo_t_yfs = Histogram_2D(111,1e-6,1e4,100,1e-6,1e2,20); Histogram_2D PHOTONS::Photons::s_histo_t_total = Histogram_2D(111,1e-6,1e4,100,1e-6,1e2,20); #endif // member functions of class Photons Photons::Photons() : m_name("Photons") { RegisterDefaults(); Scoped_Settings s{ Settings::GetMainSettings()["YFS"] }; rpa->gen.AddCitation(1, "Photons is published under \\cite{Schonherr:2008av}."); s_mode = s["MODE"].Get(); s_useme = (bool)s["USE_ME"].Get(); s_ircutoff = s["IR_CUTOFF"].Get(); s_uvcutoff = s["UV_CUTOFF"].Get(); s_alpha_input = s["1/ALPHAQED"].Get(); s_alpha_input = (s_alpha_input?1./s_alpha_input:MODEL::aqed->AqedThomson()); s_userunningparameters = (bool)s["USE_RUNNING_PARAMETERS"].Get(); std::string irframe = s["IR_CUTOFF_FRAME"].Get(); if (irframe == "Multipole_CMS") s_ircutoffframe = 0; else if (irframe == "Lab") s_ircutoffframe = 1; else if (irframe == "Decayer_Rest_Frame") s_ircutoffframe = 2; else { s_ircutoffframe = 0; msg_Info()<<"value '"<(); s_nmin = s["MINEM"].Get(); s_drcut = s["DRCUT"].Get(); s_strict = s["STRICTNESS"].Get(); s_reducemaxenergy = s["REDUCE_MAXIMUM_ENERGY"].Get(); s_increasemaxweight = s["INCREASE_MAXIMUM_WEIGHT"].Get(); s_checkfirst = (bool)s["CHECK_FIRST"].Get(); s_ffrecscheme = s["FF_RECOIL_SCHEME"].Get(); s_firecscheme = s["FI_RECOIL_SCHEME"].Get(); s_accu = sqrt(rpa->gen.Accu()); m_splitphotons = s["PHOTON_SPLITTER_MODE"].Get(); m_photonsplitter = Photon_Splitter(m_splitphotons); #ifdef PHOTONS_DEBUG s_histo_base_name = s["HISTO_BASE_NAME"].Get(); #endif m_success = true; m_photonsadded = false; msg_Debugging()<0) { msg_Debugging()<<" , MEs: "<<((int)s_mode>1?s_useme:0) <<" , nmax: "<0?s_ircutoff:0) <<" in frame "<::max()); s["1/ALPHAQED"].SetDefault(0.); s["USE_RUNNING_PARAMETERS"].SetDefault(0); s["IR_CUTOFF_FRAME"].SetDefault("Multipole_CMS"); s["MAXEM"].SetDefault(std::numeric_limits::max()); s["MINEM"].SetDefault(0); s["DRCUT"].SetDefault(std::numeric_limits::max()); s["STRICTNESS"].SetDefault(0); s["REDUCE_MAXIMUM_ENERGY"].SetDefault(1.); s["INCREASE_MAXIMUM_WEIGHT"].SetDefault(1.); s["CHECK_FIRST"].SetDefault(0); s["FF_RECOIL_SCHEME"].SetDefault(2); s["FI_RECOIL_SCHEME"].SetDefault(2); s["HISTO_BASE_NAME"].SetDefault("weights"); } bool Photons::AddRadiation(Blob * blob) { if (s_mode==yfsmode::off) return m_success=true; if (!CheckStateBeforeTreatment(blob)) { m_photonsadded=false; return m_success=false; } ResetAlphaQED(); Define_Dipole dress(blob); if (!dress.CheckMasses()) { msg_Error()<MomentumConserved()) { msg_Error()<CheckMomentumConservation()<MomentumConserved()) { msg_Error()<CheckMomentumConservation()<ShortProcessName()); if (!blob->MomentumConserved()) { msg_Error()<CheckMomentumConservation()<NOutP();++i) { if (blob->OutParticle(i)->FinalMass()==0. && blob->OutParticle(i)->Momentum().Mass()>1e-3) { msg_Debugging()<OutParticle(i)->Flav().IDName() <<" not onshell: " <OutParticle(i)->Momentum().Mass()<<" vs " <OutParticle(i)->FinalMass()<OutParticle(i)->FinalMass()!=0. && !IsEqual(blob->OutParticle(i)->Momentum().Mass(), blob->OutParticle(i)->FinalMass(),1e-3)) { msg_Debugging()<OutParticle(i)->Flav().IDName() <<" not onshell: " <OutParticle(i)->Momentum().Mass()<<" vs " <OutParticle(i)->FinalMass()<