#include "SHRiMPS/Beam_Remnants/Hadron_Dissociation.H" #include "SHRiMPS/Tools/MinBias_Parameters.H" #include "ATOOLS/Math/Random.H" #include "ATOOLS/Org/Message.H" #include using namespace SHRIMPS; using namespace ATOOLS; using namespace std; Hadron_Dissociation::Hadron_Dissociation(Continued_PDF *const pdf) : p_pdf(pdf), m_bunch(p_pdf->Bunch()), //m_ycut(MBpars("originalY")), //-MBpars("deltaY")), m_analyse(true) { if (m_analyse) { m_histomap[string("KT_remn_orig")] = new Histogram(0,0.0,5.0,50); m_histomap[string("KT_remn_resc")] = new Histogram(0,0.0,5.0,50); m_histomap[string("X_quark")] = new Histogram(0,0.0,1.0,1000); m_histomap[string("X_gluon")] = new Histogram(0,0.0,1.0,1000); m_histomap[string("X_diquark")] = new Histogram(0,0.0,1.0,1000); m_histomap2D[string("X_quark_2D")] = new Histogram_2D(0,0.,25.,25,0.0,1.0,100); m_histomap2D[string("X_gluon_2D")] = new Histogram_2D(0,0.,25.,25,0.0,1.0,100); m_histomap2D[string("X_diquark_2D")] = new Histogram_2D(0,0.,25.,25,0.0,1.0,100); } } Hadron_Dissociation::~Hadron_Dissociation() { if (m_analyse) { msg_Info()<<"Initial kt's: "<Average() <<" and "<Average() <<" after rescaling.\n"; if (!m_histomap.empty()) { Histogram * histo; string name; for (map::iterator hit=m_histomap.begin();hit!=m_histomap.end();hit++) { histo = hit->second; name = string("Ladder_Analysis/")+hit->first+string(".dat"); histo->Finalize(); histo->Output(name); delete histo; } m_histomap.clear(); } if (!m_histomap2D.empty()) { Histogram_2D * histo; string name; for (map::iterator hit=m_histomap2D.begin();hit!=m_histomap2D.end();hit++) { histo = hit->second; name = string("Ladder_Analysis/")+hit->first+string(".dat"); histo->Finalize(); histo->Output(name); delete histo; } m_histomap.clear(); } } } void Hadron_Dissociation::FixFlavourConstituents() { double random(ran->Get()); if (m_bunch==Flavour(kf_p_plus)) { if (random<1./3.) { m_quark = Flavour(kf_d); m_remnant = Flavour(kf_uu_1); } else if (random<1./2.) { m_quark = Flavour(kf_u); m_remnant = Flavour(kf_ud_1); } else { m_quark = Flavour(kf_u); m_remnant = Flavour(kf_ud_0); } } else if (m_bunch==Flavour(kf_p_plus).Bar()) { if (random<1./3.) { m_quark = Flavour(kf_d).Bar(); m_remnant = Flavour(kf_uu_1).Bar(); } else if (random<1./2.) { m_quark = Flavour(kf_u).Bar(); m_remnant = Flavour(kf_ud_1).Bar(); } else { m_quark = Flavour(kf_u).Bar(); m_remnant = Flavour(kf_ud_0).Bar(); } } else { msg_Error()<<"Error in "<SetNumber(0); m_particles.push_back(particle); for (int i=0;iSetNumber(0); m_particles.push_back(particle); } random_shuffle(m_particles.begin(),m_particles.end(),*ran); particle = new Particle(0,m_remnant,defmom,'B'); particle->SetNumber(0); m_particles.push_back(particle); } bool Hadron_Dissociation:: DefineDissociation(const int & Nladders,const double B, const double & xcut,const double & eta, Form_Factor * ff) { //msg_Out()<XMin()*double(npart+1)),xave(1./double(npart+1)); double startweight(pow(1.25,npart-1)*pow(xave,-(npart+1)*eta)); if (xminGet(); m_xs.push_back(x); } for (int i=0;iFlav(); if (flav.IsDiQuark()) { if (x/2XMin()) { weight = 0.; break; } weight *= wt = exp((x-1.)/Nladders); //pow(1.+xcut-x,-npart+1);// } else if (flav.IsQuark()) { if (xXMin()) { weight = 0.; break; } p_pdf->Calculate(x,0.); weight *= wt = p_pdf->XPDF(flav)/p_pdf->XPDFMax(flav)/x; } // the extra x in xpdf compensates for the incoming flux. if (i!=npart && !IsZero(eta)) weight *= pow(x,eta); } if (weight>maxwt) maxwt=weight; if (weight>1.) { msg_Tracking()<<"\n"; msg_Tracking()<<" Potential Error in "<1 " <<"(start = "<ran->Get()) { for (size_t i=0;iSetMomentum(m_xs[i] * m_bunchmom); if (m_analyse) { if (m_particles[i]->Flav().IsGluon()){ m_histomap["X_gluon"]->Insert(m_xs[i]); m_histomap2D["X_gluon_2D"]->Insert(npart-1,m_xs[i]); } else if (m_particles[i]->Flav().IsQuark()){ m_histomap["X_quark"]->Insert(m_xs[i]); m_histomap2D["X_quark_2D"]->Insert(npart-1,m_xs[i]); } else if (m_particles[i]->Flav().IsDiQuark()){ m_histomap["X_diquark"]->Insert(m_xs[i]); m_histomap2D["X_diquark_2D"]->Insert(npart-1,m_xs[i]); } } } DefineTransverseMomenta(ff); return true; } } msg_Tracking()<<"\n"; msg_Error()< qts, qls; m_qtvecs.clear(); for (size_t i=0;iSelectQT2(QT2max,QT2min); } while (qt2>QT2max || qt2Insert(qt); } phi = ran->Get()*2.*M_PI; m_qtvecs.push_back(qt*Vec4D(0.,cos(phi),sin(phi),0.)); sumvec += m_qtvecs.back(); KL += qls.back(); } for (size_t i=0;iInsert(qt); } } } void Hadron_Dissociation::FillBeamBlob() { p_blob->SetType(btp::Beam); p_blob->SetTypeSpec("Shrimps"); p_blob->SetStatus(blob_status::inactive); if (!m_elastic) { for (size_t i=0;iAddToOutParticles(m_particles[i]); } } //msg_Out()<SetBeam(beam); blob->AddToInParticles(m_particles[i]); } } void Hadron_Dissociation::PrintParticles() const { msg_Out()<GetFlow(pos)==c1) { m_particles[i]->SetFlow(pos,c2); return true; } } return false; }