#include "SHRiMPS/Ladders/Ladder_Generator_Base.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 using namespace SHRIMPS; using namespace MODEL; using namespace ATOOLS; using namespace std; Ladder_Generator_Base::Ladder_Generator_Base() : m_partonic(Sigma_Partonic(xs_mode::perturbative)), m_Ymax(MBpars.GetEikonalParameters().Ymax), m_deltaY(MBpars.GetEikonalParameters().cutoffY), m_beamY(m_Ymax+m_deltaY), m_qt2min(MBpars.GetLadderParameters().Q02), m_qt2minFF(0.), m_kt2min(MBpars.GetLadderParameters().Q02), m_kt2minShower(MBpars.GetShowerLinkParameters().KT2min), m_density(MBpars.GetEikonalParameters().Delta, MBpars.GetEikonalParameters().lambda,m_Ymax, MBpars.GetEikonalParameters().absorp), p_ladder(0) { Running_AlphaS * as = static_cast(s_model->GetScalarFunction(string("alpha_S"))); p_alphaS = new Strong_Coupling(as,asform::smooth,MBpars.GetLadderParameters().Qas2); m_partonic.SetAlphaS(p_alphaS); m_me.SetPartonic(&m_partonic); } Ladder_Generator_Base::~Ladder_Generator_Base() { delete p_alphaS; } void Ladder_Generator_Base::Initialise(Remnant_Handler * remnants) { m_partonic.Initialise(remnants); } void Ladder_Generator_Base::InitLadder(const Vec4D & pos) { m_shat = 4.*m_E[0]*m_E[1]; m_weight = 1.; p_ladder = new Ladder(pos); p_emissions = p_ladder->GetEmissions(); p_props = p_ladder->GetProps(); } void Ladder_Generator_Base::ConstructSimpleLadder() { size_t dir = (dabs(p_emissions->begin()->first) > dabs(p_emissions->rbegin()->first) ? 0 : 1); double qt2max = sqr(m_E[dir]/ (cosh(dir==0 ? p_ladder->GetEmissions()->begin()->first : p_ladder->GetEmissions()->rbegin()->first))); do { m_qt2 = p_eikonal->FF(dir)->SelectQT2(qt2max); } while (ReggeWeight(m_qt2,m_ylimits[0],m_ylimits[1]) < ran->Get()); MakeTransverseUnitVector(); Vec4D k[2]; k[0] = sqrt(m_qt2)*(Vec4D(cosh(m_ylimits[0]),0.,0.,sinh(m_ylimits[0])) + m_eqt); k[1] = sqrt(m_qt2)*(Vec4D(cosh(m_ylimits[1]),0.,0.,sinh(m_ylimits[1])) - m_eqt); p_emissions->begin()->second.SetMomentum(k[0]); p_emissions->rbegin()->second.SetMomentum(k[1]); T_Prop & prop = *p_props->begin(); prop.SetQT2(m_qt2); prop.SetQ02(m_qt2min); prop.SetQ(sqrt(m_qt2)*m_eqt); } void Ladder_Generator_Base::ConstructISKinematics() { Vec4D Ksum = p_ladder->FSMomentum(); for (size_t i=0;i<2;i++) { p_ladder->InPart(i)->SetMomentum(i==0 ? Ksum.PPlus()/2. * Vec4D(1.,0.,0., 1.) : Ksum.PMinus()/2. * Vec4D(1.,0.,0.,-1.) ); p_ladder->InPart(i)->SetBeam(i); } p_emissions->begin()->second.SetBeam(p_ladder->InPart(0)->Beam()); p_emissions->rbegin()->second.SetBeam(p_ladder->InPart(1)->Beam()); p_ladder->UpdatePropagatorKinematics(); } double Ladder_Generator_Base::RescaleLadder(Ladder * ladder,const Vec4D & P_in) { Vec4D P_out = ladder->FSMomentum(); double factor = sqrt(P_in.Abs2()/P_out.Abs2()), weight = 1.; Poincare out = Poincare(P_out), in = Poincare(P_in); for (LadderMap::iterator lit=ladder->GetEmissions()->begin(); lit!=ladder->GetEmissions()->end();lit++) { Vec4D oldp = lit->second.Momentum(); out.Boost(oldp); in.BoostBack(oldp); Vec4D newp = factor * oldp; lit->second.SetMomentum(newp); if (dabs(newp.Y())GetProps()->begin(); pit!=ladder->GetProps()->end();pit++) { Vec4D oldq = pit->Q(), newq = factor*oldq; double oldqt2 = pit->QT2(), newqt2 = sqr(factor)*oldqt2; pit->SetQ(newq); pit->SetQT2(newqt2); weight *= oldqt2/newqt2; } return weight; } double Ladder_Generator_Base::TotalReggeWeight(Ladder * ladder) { LadderMap::iterator lit1=ladder->GetEmissions()->begin(), lit2=lit1; lit2++; TPropList::iterator pit=ladder->GetProps()->begin(); double weight = 1.; while (lit2!=ladder->GetEmissions()->end() && pit!=ladder->GetProps()->end()) { weight *= ReggeWeight(dabs(pit->Q2()),lit1->first,lit2->first); pit++; lit1++; lit2++; } return weight; } void Ladder_Generator_Base::MakeTransverseUnitVector() { double phi = 2.*M_PI*ran->Get(); m_eqt = Vec4D(0.,cos(phi),sin(phi),0.); } void Ladder_Generator_Base::ResetFSFlavours() { for (LadderMap::iterator lit=p_emissions->begin(); lit!=p_emissions->end();lit++) { lit->second.SetFlavour(Flavour(kf_gluon)); } for (TPropList::iterator pit=p_props->begin();pit!=p_props->end();pit++) { if (pit->Col()==colour_type::triplet) pit->SetCol(colour_type::octet); } } void Ladder_Generator_Base::QuarkReplace() { LadderMap::iterator lit1, lit2; TPropList::iterator pit; double pp, ppold; int dir; if (ran->Get() > .5) dir = 0; else dir = 1; if (dir == 0) { lit1=p_emissions->begin(); lit2=lit1; lit2++; pit=p_props->begin(); pp = p_ladder->InPart(0)->Momentum().PPlus(); } else { lit1=p_emissions->end(); lit1--; lit2=lit1; lit2--; pit=p_props->end(); pit--; pp = p_ladder->InPart(1)->Momentum().PMinus(); } ppold=pp; bool last = false, stop = false; do { if (dir == 0) pp = pit->Q().PPlus(); else pp = pit->Q().PMinus(); if (last) last=false; else { if (pit->Col()!=colour_type::singlet) { if (pp > ran->Get()*ppold) { pit->SetCol(colour_type::triplet); Flavour flav = Flavour(int(1+ran->Get()*3.)); if (dir==0){//ran->Get()>.5) { lit1->second.SetFlavour(flav); lit2->second.SetFlavour(flav.Bar()); lit1->second.SetFlow(2,0); lit2->second.SetFlow(1,0); } else { lit2->second.SetFlavour(flav); lit1->second.SetFlavour(flav.Bar()); lit2->second.SetFlow(2,0); lit1->second.SetFlow(1,0); } last = true; } } } ppold = pp; if(dir == 0) { lit1++; lit2++; pit++; stop = pit==p_props->end(); } else { stop = pit == p_props->begin(); if(!stop) {lit1--; lit2--; pit--;} } } while (!stop); } double Ladder_Generator_Base::AlphaSWeight(const double & kt2) { return AlphaS(kt2)/AlphaSMax(); } double Ladder_Generator_Base:: ReggeWeight(const double & qt2, const double & y1,const double y2) { return exp(-3.*AlphaS(qt2)/M_PI * dabs(y1-y2) * log(1.+qt2/m_qt2min)); } double Ladder_Generator_Base:: LDCWeight(const double & qt2, const double & qt2prev) { return qt2/Max(qt2,qt2prev); } void Ladder_Generator_Base::Test() { vector > * eikonals(MBpars.GetEikonals()); double b1, b2, y, asym12, asym34, d1, d2, d3, d4; for (size_t i=0;isize();i++) { for (size_t j=i;j<(*eikonals)[i].size();j++) { msg_Out()<<"=================================\n" <<"Testing eikonals["<