#include "PHASIC++/Channels/Exponential_Channels.H" #include "PHASIC++/Channels/Channel_Elements.H" #include "ATOOLS/Org/MyStrStream.H" #include "ATOOLS/Org/Message.H" #include <stdio.h> using namespace PHASIC; using namespace ATOOLS; using namespace std; Exponential_RelicDensity:: Exponential_RelicDensity(const double exponent,const double mass1,const double mass2, const std::string cinfo,ATOOLS::Integration_Info *info): ISR_Channel_Base(info), m_exponent(exponent) { m_mass[0]=mass1; m_mass[1]=mass2; m_name="Exponential_"+ATOOLS::ToString(exponent)+"_RelicDensity"; m_spkey.SetInfo(std::string("Exponential_")+ATOOLS::ToString(exponent)); m_spkey.Assign(cinfo+std::string("::s'"),5,0,info); m_sgridkey.Assign(m_spkey.Info(),1,0,info); m_zchannel = m_spkey.Name().find("z-channel")!=std::string::npos; m_rannum = 1; p_vegas = new Vegas(m_rannum,100,m_name); p_rans = new double[m_rannum]; } void Exponential_RelicDensity::GeneratePoint(const double *rns) { double *ran = p_vegas->GeneratePoint(rns); for(int i=0;i<m_rannum;i++) p_rans[i]=ran[i]; m_spkey[3] = CE.ExponentialMomenta(m_exponent,m_spkey[0],m_spkey[1],m_mass,p_rans[0]); } void Exponential_RelicDensity::GenerateWeight(const int & mode) { if (m_spkey.Weight()==ATOOLS::UNDEFINED_WEIGHT) { if (m_spkey[3]>=m_spkey[0] && m_spkey[3]<=m_spkey[1]) { m_spkey<<1./CE.ExponentialWeight(m_exponent,m_spkey[0],m_spkey[1],m_mass, m_spkey[3],m_sgridkey[0]); } } if (m_spkey[4]>0.0) { p_vegas->ConstChannel(0); m_spkey<<M_PI*2.0; } p_rans[0] = m_sgridkey[0]; double pw = p_vegas->GenerateWeight(p_rans); m_weight=pw*m_spkey.Weight(); } /////////////////////////////////////////////////////////////////////////////// Exponential_DM_Annihilation:: Exponential_DM_Annihilation(const double exponent,const double mass1,const double mass2, const std::string cinfo,ATOOLS::Integration_Info *info): ISR_Channel_Base(info), m_exponent(exponent) { m_mass[0]=mass1; m_mass[1]=mass2; m_name="Exponential_"+ATOOLS::ToString(exponent)+"_DM_Annihilation"; m_spkey.SetInfo(std::string("Exponential_")+ATOOLS::ToString(exponent)); m_spkey.Assign(cinfo+std::string("::s'"),5,0,info); m_sgridkey.Assign(m_spkey.Info(),1,0,info); m_xkey.Assign(cinfo+std::string("::xDM"),3,0,info); m_cosxikey.Assign(cinfo+std::string("::cosXi"),3,0,info); m_sgridkey.Assign(m_spkey.Info(),1,0,info); m_xgridkey.Assign(m_xkey.Info(),1,0,info); m_cosgridkey.Assign(m_cosxikey.Info(),1,0,info); m_zchannel = m_spkey.Name().find("z-channel")!=std::string::npos; m_rannum = 2; p_vegas = new Vegas(m_rannum,100,m_name); p_rans = new double[m_rannum]; } void Exponential_DM_Annihilation::GeneratePoint(const double *rns) { double *ran = p_vegas->GeneratePoint(rns); for(int i=0;i<m_rannum;i++) p_rans[i]=ran[i]; m_spkey[3] = m_spkey[2] = CE.ExponentialMomenta(m_exponent,m_spkey[0],m_spkey[1],m_mass,p_rans[0]); m_cosxikey[2] = CE.GenerateDMAngleUniform(p_rans[1],3); m_xkey[2] = CE.GenerateDMRapidityUniform(m_mass,m_spkey.Doubles(),m_xkey.Doubles(), m_cosxikey[2], p_rans[0], 3); } void Exponential_DM_Annihilation::GenerateWeight(const int & mode) { if (m_spkey.Weight()==ATOOLS::UNDEFINED_WEIGHT) { if (m_spkey[3]>=m_spkey[0] && m_spkey[3]<=m_spkey[1]) { m_spkey<<1./CE.ExponentialWeight(m_exponent,m_spkey[0],m_spkey[1],m_mass, m_spkey[3],m_sgridkey[0]); } } if (m_spkey[4]>0.0) { p_vegas->ConstChannel(0); m_spkey<<M_PI*2.0; } p_rans[0] = m_sgridkey[0]; double pw = p_vegas->GenerateWeight(p_rans); m_weight=pw*m_spkey.Weight()/m_spkey[2]; // msg_Out()<<"s="<<m_spkey[3]<<", weight="<<m_weight<<"\n"; //debugging }