#include "SHRiMPS/Cross_Sections/Sigma_Partonic.H" #include "SHRiMPS/Beam_Remnants/Remnant_Handler.H" #include "SHRiMPS/Tools/MinBias_Parameters.H" #include "MODEL/Main/Model_Base.H" #include "ATOOLS/Math/Random.H" #include "ATOOLS/Org/Run_Parameter.H" #include "ATOOLS/Org/Message.H" using namespace SHRIMPS; using namespace MODEL; using namespace ATOOLS; using namespace std; Sigma_Partonic::Sigma_Partonic(const xs_mode::code & mode) : p_alphaS(NULL), m_mode(mode), m_fixflavour(true), m_Ymax(MBpars.GetEikonalParameters().originalY), m_S(sqr(rpa->gen.Ecms())), m_eta(0.08), m_smin(1.), m_tmin(1.), m_accu(0.01), m_sigma(0.), m_maxdsigma(0.), m_Nmaxtrials(100), m_kinX_fails(0), m_ana(true) { if (m_ana) { m_histos[string("Yhat_sigma")] = new Histogram(0, -12.5, 12.5, 50); m_histos[string("Yhat_asym_sigma")] = new Histogram(0, 0.0, 12.5, 50); m_histos[string("Y1_sigma")] = new Histogram(0, -12.5, 12.5, 50); m_histos[string("Y2_sigma")] = new Histogram(0, -12.5, 12.5, 50); } } Sigma_Partonic::~Sigma_Partonic() { if (m_ana) { std::string name = std::string("Ladder_Analysis/"); for (std::map::iterator hit=m_histos.begin(); hit!=m_histos.end();hit++) { hit->second->Finalize(); hit->second->Output(name+hit->first); delete hit->second; } } } void Sigma_Partonic::Initialise(Remnant_Handler * remnants) { for (size_t beam=0;beam<2;beam++) p_pdf[beam]=remnants->GetPDF(beam); m_smin = Max(m_smin, m_S*p_pdf[0]->XMin()*p_pdf[1]->XMin()); if (!Calculate()) { msg_Error()<m_maxdsigma*ran->Get()) { m_phi = 2.*M_PI*ran->Get(); SelectFlavours(m_fixflavour); if (m_ana) { msg_Out()<<" - "<Insert(m_yhat); m_histos[string("Yhat_asym_sigma")]->Insert(dabs(m_yhat),(m_yhat>0.?1.:-1.)); m_histos[string("Y1_sigma")]->Insert(m_y[0]); m_histos[string("Y2_sigma")]->Insert(m_y[1]); } return m_shat; } } return -1.; } void Sigma_Partonic::SelectFlavours(const bool & fixflavour) { m_flavs[0] = m_flavs[1] = Flavour(kf_gluon); if (fixflavour) return; for (size_t i=0;i<2;i++) { double disc = 0.999999 * m_xpdf[i]*ran->Get(); for (list::const_iterator flit=p_pdf[i]->GetFlavours().begin(); flit!=p_pdf[i]->GetFlavours().end();flit++) { if (p_pdf[i]->XPDF((*flit))<1.e-4) continue; disc -= p_pdf[i]->XPDF((*flit)) * ColourFactor((*flit)); if (disc<=0.) { m_flavs[i] = (*flit); break; } } } } const double Sigma_Partonic:: PDF(const size_t beam,const double & x,const double & Q2,const Flavour & flav) { if (x<1.e-6) return 0.; p_pdf[beam]->Calculate(x,Q2); return p_pdf[beam]->XPDF(flav); } const bool Sigma_Partonic::Calculate() { size_t iter = 0, Npoints = 10000; long int N = 0; double dsigma, sigma, res = 0., res2 = 0., accu = 1.e99; do { for (size_t i=0;im_maxdsigma) m_maxdsigma = dsigma; } sigma = res/double(N); accu = sqrt((res2/double(N) - sqr(sigma))/double(N))/sigma; if (accu " <m_accu); msg_Out()<Get(); switch (m_mode) { case xs_mode::perturbative: for (size_t i=0;i<2;i++) { m_y[i] = m_Ymax*(-1.+2.*ran->Get()); } if (m_y[0]==m_y[1]) { m_y[0]+=1.e-6; m_y[1]-=1.e-6; } else if (m_y[0]Get()); ymax = 1./2.*log(m_S/m_shat); m_yhat = ymax * (-1. + 2.*rand); return log(m_S/m_smin) * (2.*ymax); default: break; } return 0.; } const double Sigma_Partonic::dSigma() { double flux = 1/(2.*m_shat), Enorm = sqrt(m_shat/m_S), scale = m_that; for (size_t i=0;i<2;i++) { m_x[i] = Enorm * (i==0 ? exp(m_yhat) : exp(-m_yhat)); if (m_x[i]XMin() || m_x[i]>0.9999999) { m_kinX_fails++; return 0.; } m_xpdf[i] = 0.; p_pdf[i]->Calculate(m_x[i],scale); for (list::const_iterator flit=p_pdf[i]->GetFlavours().begin(); flit!=p_pdf[i]->GetFlavours().end();flit++) { if (p_pdf[i]->XPDF((*flit))>1.e-4) m_xpdf[i] += p_pdf[i]->XPDF((*flit)) * ColourFactor((*flit)); } } return flux * m_xpdf[0]*m_xpdf[1] * ME2(m_shat,m_that,scale); } const double Sigma_Partonic::dSigma(const double & shat,const double & yhat) { double flux = 1/(2.*m_shat), Enorm = sqrt(m_shat/m_S), scale = 0.; for (size_t i=0;i<2;i++) { m_x[i] = Enorm * (i==0 ? exp(m_yhat) : exp(-m_yhat)); if (m_x[i]<1.e-6 || m_x[i]>0.9999999) return 0.; m_xpdf[i] = 0.; p_pdf[i]->Calculate(m_x[i],scale); for (list::const_iterator flit=p_pdf[i]->GetFlavours().begin(); flit!=p_pdf[i]->GetFlavours().end();flit++) { if (p_pdf[i]->XPDF((*flit))<1.e-4) continue; m_xpdf[i] += p_pdf[i]->XPDF((*flit)) * ColourFactor((*flit)); } } return flux * m_xpdf[0]*m_xpdf[1] * ME2(m_shat,m_that,scale); } const double Sigma_Partonic::ME2(const double & shat,const double & that,const double & scale) { switch (m_mode) { case xs_mode::perturbative: return sqr(4.*M_PI*(*p_alphaS)(scale))*(sqr(shat)+sqr(shat+that))/sqr(-that+m_tmin); case xs_mode::integrated: return 4.*M_PI*sqr((*p_alphaS)(scale))*sqr(shat)/(m_tmin*(shat+m_tmin)); case xs_mode::Regge: return pow( (shat/m_smin), 1.+m_eta); default: break; } return 0.; }