#include "PHOTONS++/PhaseSpace/Avarage_Photon_Number.H" #include "ATOOLS/Math/Random.H" #include "ATOOLS/Math/Vector.H" #include "ATOOLS/Org/Message.H" #include "ATOOLS/Phys/Particle.H" #include "PHOTONS++/Main/Photons.H" #include "PHOTONS++/Main/Dipole_Type.H" using namespace PHOTONS; using namespace ATOOLS; using namespace std; Avarage_Photon_Number::Avarage_Photon_Number (const Particle_Vector& dip, const double& wmax, const double& wmin) : m_omegaMax(wmax), m_omegaMin(wmin), m_dipole(dip), m_nbar(0.) { DEBUG_FUNC(dip.size()); CalculateAvaragePhotonNumber(); msg_Debugging()<<"all dipoles:\n"; for (unsigned int i=0; iFlav()<<"," <Flav()<<"]"< "<Flav().Charge(); double Zj = m_dipole[j]->Flav().Charge(); double titj = TiTj(i,j); double betai = CalculateBeta(m_dipole[i]->Momentum()); double betaj = CalculateBeta(m_dipole[j]->Momentum()); Vec3D pi = m_dipole[i]->Momentum(); Vec3D pj = m_dipole[j]->Momentum(); double alpha = 0.; if (!IsZero(pi.Abs()) && !IsZero(pj.Abs())) { double arg = pi*pj/(pi.Abs()*pj.Abs()); if (IsEqual(arg,-1.)) alpha = 0.; else if (IsEqual(arg,1.)) alpha = M_PI/2.; else alpha = (M_PI - acos(arg))/2.; } double ai = betai*sin(alpha); double aj = betaj*sin(alpha); double bi = betai*cos(alpha); double bj = betaj*cos(alpha); double intterm = 0.; if (IsZero(alpha,2.e-4)) { intterm = log(((1.+betai)*(1.+betaj))/((1.-betai)*(1.-betaj))) /(betai+betaj); msg_Debugging()<<"back-to-back pair discovered: "< "<0) || !(2.*sqrt(Cj-Dj+Ej)+(B*(2.*Cj-Dj)-A*(Dj-2.*Ej))/rootj>0) || !(2.*sqrt(Ci+Di+Ei)+(B*(2.*Ci+Di)-A*(Di+2.*Ei))/rooti>0) || !(2.*sqrt(Cj+Dj+Ej)+(B*(2.*Cj+Dj)-A*(Dj+2.*Ej))/rootj>0)) return 0.; else return bi*bj*(1./(bj*rooti) *log(abs((A+B)/(A-B)) *(2.*sqrt(Ci-Di+Ei)+(B*(2.*Ci-Di)-A*(Di-2.*Ei))/rooti) /(2.*sqrt(Ci+Di+Ei)+(B*(2.*Ci+Di)-A*(Di+2.*Ei))/rooti)) -1./(bi*rootj) *log(abs((A+B)/(A-B)) *(2.*sqrt(Cj-Dj+Ej)+(B*(2.*Cj-Dj)-A*(Dj-2.*Ej))/rootj) /(2.*sqrt(Cj+Dj+Ej)+(B*(2.*Cj+Dj)-A*(Dj+2.*Ej))/rootj))); } double Avarage_Photon_Number::TiTj(const size_t& i, const size_t& j) { if ((m_dipole[i]->ProductionBlob() == m_dipole[j]->ProductionBlob()) && (m_dipole[i]->ProductionBlob() != NULL)) return 1.; else if ((m_dipole[i]->DecayBlob() == m_dipole[j]->ProductionBlob()) && (m_dipole[i]->DecayBlob() != NULL)) return -1.; else if ((m_dipole[i]->ProductionBlob() == m_dipole[j]->DecayBlob()) && (m_dipole[i]->ProductionBlob() != NULL)) return -1.; else if ((m_dipole[i]->DecayBlob() == m_dipole[j]->DecayBlob()) && (m_dipole[i]->DecayBlob() != NULL)) return 1.; else return 0.; }