#include "PHOTONS++/Tools/Dress_Blob_Base.H" #include "ATOOLS/Math/MathTools.H" #include "ATOOLS/Math/Vector.H" #include "ATOOLS/Phys/Particle.H" #include "ATOOLS/Phys/Blob.H" #include "ATOOLS/Math/Random.H" #include "ATOOLS/Math/Poincare.H" #include "PHOTONS++/Main/Photons.H" #include "PHOTONS++/Tools/Generate_One_Photon.H" #include "PHOTONS++/Tools/Weight_Dipole.H" #include "PHOTONS++/Tools/Weight_Jacobian.H" #include "PHOTONS++/Tools/Weight_YFS.H" #include "PHOTONS++/Tools/Weight_Higher_Order_Corrections.H" #include "PHOTONS++/Tools/Generate_One_Photon.H" #include "PHOTONS++/PhaseSpace/Avarage_Photon_Number.H" #include "ATOOLS/Org/Run_Parameter.H" #include "ATOOLS/Org/Message.H" #include #ifdef PHOTONS_DEBUG #include "ATOOLS/Math/Histogram_2D.H" #endif using namespace PHOTONS; using namespace ATOOLS; using namespace std; Dress_Blob_Base::Dress_Blob_Base() : m_photonsadded(false), m_success(false), m_soft((PHOTONS::Photons::s_mode==yfsmode::soft)?true:false), p_newinitialstate(NULL), m_accu(PHOTONS::Photons::s_accu), m_genweight(1.), m_genmaxweight(1.), m_nbar(0.), m_n(0.), m_omegaMax(PHOTONS::Photons::s_uvcutoff), m_omegaMin(PHOTONS::Photons::s_ircutoff) {} Dress_Blob_Base::~Dress_Blob_Base() { if (m_dtype == Dipole_Type::ff) delete p_newinitialstate; } void Dress_Blob_Base::BuildNewParticleVectorVector() { if ((m_dtype == Dipole_Type::ff) && (p_newinitialstate != NULL)) delete p_newinitialstate; m_pvv_new.clear(); Particle_Vector cinit, ninit; if (m_dtype == Dipole_Type::ff) { p_newinitialstate = new Particle(*m_neutralinparticles[0]); p_newinitialstate->SetMomentum(m_P+m_PN+m_K); ninit.push_back(p_newinitialstate); } else if (m_dtype == Dipole_Type::fi) { p_newinitialstate = m_newdipole[0]; cinit.push_back(p_newinitialstate); } m_pvv_new.push_back(cinit); m_pvv_new.push_back(ninit); m_pvv_new.push_back(m_newdipole); m_pvv_new.push_back(m_newspectator); m_pvv_new.push_back(m_softphotons); #ifdef PHOTONS_DEBUG msg_Debugging()<<"new particles:\n"; for (unsigned int i=0; iMomentum(); return sum; } void Dress_Blob_Base::CalculateWeights() { BuildNewParticleVectorVector(); Weight_Dipole dip(m_olddipole,m_newdipole,m_softphotons,m_dtype); Weight_YFS yfs(m_newdipole,m_olddipole,m_omegaMin,m_nbar); Weight_Jacobian_Mapping jacobM(m_newdipole,m_newspectator,m_olddipole, m_oldspectator,m_K,m_M,m_u,m_dtype); Weight_Jacobian_Lorentz jacobL(m_newdipole,m_newspectator,m_olddipole, m_oldspectator,m_softphotons,m_dtype); double wdipole = dip.Get(); double mwdipole = dip.GetMax(); double wyfs = yfs.Get(); double mwyfs = yfs.GetMax(); double wjacobianM = jacobM.Get(); double mwjacobianM= jacobM.GetMax(); double wjacobianL = jacobL.Get(); double mwjacobianL= jacobL.GetMax(); double whigher = 1.; double mwhigher = 1.; if (!m_soft) { Weight_Higher_Order_Corrections c(m_pvv,m_pvv_new,m_dtype); whigher = c.Get(); mwhigher = c.GetMax(); } m_genweight = wdipole * wjacobianM * wjacobianL * whigher * wyfs; m_genmaxweight = mwdipole * mwjacobianM * mwjacobianL * mwhigher * mwyfs; #ifdef PHOTONS_DEBUG if (IsNan(wdipole) || IsNan(mwdipole) || IsNan(wjacobianM) || IsNan(mwjacobianM) || IsNan(wjacobianL) || IsNan(mwjacobianL) || IsNan(whigher) || IsNan(mwhigher) || IsNan(wyfs) || IsNan(mwyfs) || IsNan(m_genweight) || IsNan(m_genmaxweight) || IsNan((double)m_n)) { msg_Error()<Momentum(); Poincare dec(m_newdipole[0]->Momentum()); dec.Boost(sump); Vec4D test(m_newdipole[0]->Momentum()); dec.Boost(test); double t(sump[0]); Photons::s_histo_t_dipole.Insert(wdipole/mwdipole,t); Photons::s_histo_t_jacobianM.Insert(wjacobianM/mwjacobianM,t); Photons::s_histo_t_jacobianL.Insert(wjacobianL/mwjacobianL,t); Photons::s_histo_t_higher.Insert(whigher/mwhigher,t); Photons::s_histo_t_yfs.Insert(wyfs/mwyfs,t); Photons::s_histo_t_total.Insert(m_genweight/m_genmaxweight,t); msg_Debugging()<<"weights: "< m_genmaxweight*Photons::s_increasemaxweight) { msg_Tracking()<<"weight: "< maxweight: " < mC2; std::vector mN2; std::vector q; for (unsigned int i=0; iMomentum())); for (unsigned int i=0; iMomentum())); double c = Func(M2,mC2,mN2,q,1.); if (abs(c) < 1E-12) m_u = 1.; else { double i = 0.; while ((c*Func(M2,mC2,mN2,q,1.-i) > 0.) && (i<=1)) i = i + 1E-2; if (abs(Func(M2,mC2,mN2,q,1.-i)) < 1E-14) m_u = 1.-i; else { i = i - 1E-2; double j = 0.; while ((c*Func(M2,mC2,mN2,q,1.-i-j) > 0.) && (j<=1E-2)) j = j + 1E-4; if (abs(Func(M2,mC2,mN2,q,1.-i-j)) < 1E-14) m_u = 1.-i-j; else { j = j - 1E-4; double k = 0.; while ((c*Func(M2,mC2,mN2,q,1.-i-j-k) > 0.) && (k<=1E-4)) k = k + 1E-6; if (abs(Func(M2,mC2,mN2,q,1.-i-j-k)) < 1E-14) m_u = 1.-i-j-k; else { k = k - 1E-6; double l = 0.; while ((c*Func(M2,mC2,mN2,q,1.-i-j-k-l) > 0.) && (l<=1E-6)) l = l + 1E-8; if (abs(Func(M2,mC2,mN2,q,1.-i-j-k-l)) < 1E-14) m_u = 1.-i-j-k-l; else { l = l - 1E-8; double m = 0.; while ((c*Func(M2,mC2,mN2,q,1.-i-j-k-l-m) > 0.) && (m<=1E-8)) m = m + 1E-10; if (abs(Func(M2,mC2,mN2,q,1.-i-j-k-l-m)) < 1E-14) m_u = 1.-i-j-k-l-m; else { m = m - 1E-10; double n = 0.; while ((c*Func(M2,mC2,mN2,q,1.-i-j-k-l-m-n) > 0.) && (n<=1E-10)) n = n + 1E-12; if (abs(Func(M2,mC2,mN2,q,1.-i-j-k-l-m-n)) < 1E-14) m_u = 1.-i-j-k-l-m-n; else { n = n - 1E-12; double o = 0.; while ((c*Func(M2,mC2,mN2,q,1.-i-j-k-l-m-n-o) > 0.) && (o<=1E-12)) o = o + 1E-14; if (abs(Func(M2,mC2,mN2,q,1.-i-j-k-l-m-n-o)) < 1E-14) m_u = 1.-i-j-k-l-m-n-o; else { o = o - 1E-14; double p = 0.; while ((c*Func(M2,mC2,mN2,q,1.-i-j-k-l-m-n-o-p) > 0.) && (p < 1E-14)) p = p + 1E-16; if (p < 1E-14) m_u = 1.-i-j-k-l-m-n-o-p+0.5E-16; else m_u = -1.; } } } } } } } } msg_Debugging()<<"u: "< Dress_Blob_Base::GenerateNumberAndEnergies() { // std::vector R, k; // R.push_back(0); // for (unsigned int i=0; R.back()Get())); // k.push_back(m_omegaMax*pow(m_omegaMin/m_omegaMax,R.back()/m_nbar)); // if (k.size()>Photons::s_nmax) break; // } // k.pop_back(); // m_n = k.size(); m_n = -1; double expnbar = exp(-m_nbar); double prod = 1.; while (true) { m_n++; prod = prod*ran->Get(); if (prod <= expnbar) break; } m_n = Min(Max(m_n,Photons::s_nmin),Photons::s_nmax); std::vector k; for (unsigned int i=0; iGet())); } if (m_n && msg_LevelIsDebugging()) { msg_Debugging()<<"k_i = ["< k; k = GenerateNumberAndEnergies(); for (unsigned int i=0; i k; k = GenerateNumberAndEnergies(); for (unsigned int i=0; iFlav().IDName()+" "; for (size_t i(0);iFlav().IDName()+" "; str+=" -> "; for (size_t i(0);iFlav().IDName()+" "; for (size_t i(0);iFlav().IDName()+" "; return str; }