diff --git a/MatrixElement/Hadron/MEDiffraction.cc b/MatrixElement/Hadron/MEDiffraction.cc --- a/MatrixElement/Hadron/MEDiffraction.cc +++ b/MatrixElement/Hadron/MEDiffraction.cc @@ -1,890 +1,953 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the MEDiffraction class. // #include "MEDiffraction.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" +#include "ThePEG/Interface/Reference.h" #include "ThePEG/Utilities/SimplePhaseSpace.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Handlers/StandardXComb.h" - +#include "ThePEG/Handlers/SamplerBase.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace Herwig; #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/MatrixElement/Tree2toNDiagram.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "Herwig/Utilities/Kinematics.h" MEDiffraction::MEDiffraction() : HwMEBase(), deltaOnly(false), isInRunPhase(false), theProtonMass(-MeV) // to be set in doinit {} void MEDiffraction::getDiagrams() const { //incoming particles cPDPair incomingHardons = generator()->eventHandler()->incoming(); tcPDPtr pom = getParticleData(990); //get incoming particles tcPDPtr prt11 = getParticleData(incomingHardons.first->id()); tcPDPtr prt12 = getParticleData(incomingHardons.second->id()); //get sign of id int sign1=0, sign2=0; sign1 = (incomingHardons.first->id() > 0) ? 1 : -1; sign2 = (incomingHardons.second->id() > 0) ? 1 : -1; tcPDPtr prt21 = getParticleData(sign1*2214);//Delta+ tcPDPtr prt22 = getParticleData(sign2*2214);//Delta+ //for the left side tcPDPtr q11 = getParticleData(sign1*2); //u tcPDPtr q21 = getParticleData(sign1*1); //d //for the right side tcPDPtr q12 = getParticleData(sign2*2); //u tcPDPtr q22 = getParticleData(sign2*1); //d //for the left side tcPDPtr dq11 = getParticleData(sign1*2101); //ud_0 tcPDPtr dq111 = getParticleData(sign1*2103); //ud_1 tcPDPtr dq21 = getParticleData(sign1*2203); //uu_1 //for the right side tcPDPtr dq12 = getParticleData(sign2*2101); //ud_0 tcPDPtr dq112 = getParticleData(sign2*2103); //ud_1 tcPDPtr dq22 = getParticleData(sign2*2203); //uu_1 tcPDPtr gl = getParticleData(21);//gluon //switch between dissociation decays to different //number of clusters or dissociation into delta only //(Maybe can be automated???) //(Should be generalized to ppbar, for example!!!) switch(dissociationDecay){ case 0: //one cluster or only delta in the final state if(deltaOnly) //only delta in the final state { switch (diffDirection){ case 0: add(new_ptr((Tree2toNDiagram(3), prt11, pom, prt12, 1, prt21, 2, prt12, -1))); break; case 1: add(new_ptr((Tree2toNDiagram(3), prt11, pom, prt12, 1, prt11, 2, prt22, -1))); break; case 2: add(new_ptr((Tree2toNDiagram(3), prt11, pom, prt12, 1, prt21, 2, prt22, -1))); break; } }else { //switch between direction of dissociated proton for single diffraction or //double diffraction switch (diffDirection){ case 0: //left //u -- ud_0 add(new_ptr((Tree2toNDiagram(4), prt11, q11, pom, prt12, 3, prt12, 1, dq11, 2, q11, -1))); //d -- uu_1 add(new_ptr((Tree2toNDiagram(4), prt11, q21, pom, prt12, 3, prt12, 1, dq21, 2, q21, -2))); break; case 1: //right //u -- ud_0 add(new_ptr((Tree2toNDiagram(4), prt11, pom, q12, prt12, 1, prt11, 3, dq12, 2, q12, -1))); //d -- uu_1 add(new_ptr((Tree2toNDiagram(4), prt11, pom, q22, prt12, 1, prt11, 3, dq22, 2, q22, -2))); break; case 2: //double //u -- ud_0 left u -- ud_0 right add(new_ptr((Tree2toNDiagram(5), prt11, q11, pom, q12, prt12, 1, dq11, 2, q11, 3, q12, 4, dq12, -1))); //u -- ud_0 left d -- uu_1 right add(new_ptr((Tree2toNDiagram(5), prt11, q11, pom, q22, prt12, 1, dq11, 2, q11, 3, q22, 4, dq22, -2))); //d -- uu_1 left u -- ud_0 right add(new_ptr((Tree2toNDiagram(5), prt11, q21, pom, q12, prt12, 1,dq21, 2, q21, 3, q12, 4, dq12, -3))); //d -- uu_1 left d -- uu_1 right add(new_ptr((Tree2toNDiagram(5), prt11, q21, pom, q22, prt12, 1, dq21, 2, q21, 3, q22, 4, dq22, -4))); break; } } break; case 1: //two clusters (cases with ud_1 not included) switch (diffDirection){ case 0: //left //u -- ud_0 add(new_ptr((Tree2toNDiagram(5), prt11, q11, gl, pom, prt12, 1, dq11, 2, q11, 3, gl, 4, prt12, -1))); //d -- uu_1 add(new_ptr((Tree2toNDiagram(5), prt11, q21, gl, pom, prt12, 1, dq21, 2, q21, 3, gl, 4, prt12, -2))); break; case 1: //right //u -- ud_0 add(new_ptr((Tree2toNDiagram(5), prt11, pom, gl, q12, prt12, 1, prt11, 2, gl, 3, q12, 4, dq12, -1))); //d -- ud_1 add(new_ptr((Tree2toNDiagram(5), prt11, pom, gl, q22, prt12, 1, prt11, 2, gl, 3, q22, 4, dq22, -2))); break; case 2: //double //u -- ud_0 left u -- ud_0 right add(new_ptr((Tree2toNDiagram(7), prt11, q11, gl, pom, gl, q12, prt12, 1, dq11, 2, q11, 3, gl, 4, gl, 5, q12, 6, dq12, -1))); //u -- ud_0 left d -- uu_1 right add(new_ptr((Tree2toNDiagram(7), prt11, q11, gl, pom, gl, q22, prt12, 1, dq11, 2, q11, 3, gl, 4, gl, 5, q22, 6, dq22, -2))); //d -- uu_1 left u -- ud_0 right add(new_ptr((Tree2toNDiagram(7), prt11, q21, gl, pom, gl, q12, prt12, 1, dq21, 2, q21, 3, gl, 4, gl, 5, q12, 6, dq12, -3))); //d -- uu_1 left d -- uu_1 right add(new_ptr((Tree2toNDiagram(7), prt11, q21, gl, pom, gl, q22, prt12, 1, dq21, 2, q21, 3, gl, 4, gl, 5, q22, 6, dq22, -4))); break; } break; } } Energy2 MEDiffraction::scale() const { return sqr(10*GeV); } int MEDiffraction::nDim() const { return 0; } void MEDiffraction::setKinematics() { HwMEBase::setKinematics(); // Always call the base class method first } bool MEDiffraction::generateKinematics(const double * ) { // generate the masses of the particles for (size_t i = 2; i < meMomenta().size(); ++i) meMomenta()[i] = Lorentz5Momentum(mePartonData()[i]->generateMass()); /* sample M12, M22 and t, characterizing the diffractive final state */ const pair,Energy2> point = diffractiveMassAndMomentumTransfer(); const Energy2 M12 (point.first.first); const Energy2 M22 (point.first.second); const Energy2 t(point.second); /* construct the hadronic momenta in the lab frame */ const double phi = UseRandom::rnd() * Constants::twopi; const Energy cmEnergy = generator()->maximumCMEnergy(); const Energy2 s = sqr(cmEnergy); //proton mass const Energy2 m2 = sqr( theProtonMass ); const Energy E3 = (s - M22 + M12) / (2.*cmEnergy); const Energy E4 = (s + M22 - M12) / (2.*cmEnergy); //Momentum of outgoing proton and dissociated proton const Energy pprime = sqrt(kallen(s, M12, M22)) / (2.*cmEnergy); //costheta of scattering angle double costheta = s*(s + 2*t - 2*m2 - M12 - M22) / sqrt( kallen(s, M12, M22)*kallen(s, m2, m2) ); assert(abs(costheta)<=1.); const Energy pzprime = pprime*costheta; const Energy pperp = pprime*sqrt(1 - sqr(costheta)); /* momenta in the lab frame */ const Lorentz5Momentum p3 = Lorentz5Momentum(pperp*cos(phi), pperp*sin(phi), pzprime, E3); const Lorentz5Momentum p4 = Lorentz5Momentum(-pperp*cos(phi), -pperp*sin(phi), -pzprime, E4); /* decay dissociated proton into quark-diquark */ //squares of constituent masses of quark and diquark const Energy2 mq2(sqr(mq())); Energy2 Mx2; switch(diffDirection){ case 0: Mx2=M12; break; case 1: Mx2=M22; break; } /* Select between left/right single diffraction and double diffraction */ //check if we want only delta for the excited state //pair of momenta for double decay for a two cluster case pair momPair, momPair1; //fraction of momenta double frac = UseRandom::rnd(); switch(dissociationDecay){ case 0: if(!deltaOnly) { pair decayMomenta; pair decayMomentaTwo; //absolute collinear dissociation of the hadron const double phiprime = phi; //aligned with outgoing dissociated proton const double costhetaprime = costheta; const double sinthetaprime=sqrt(1-sqr(costhetaprime)); //axis along which diquark from associated proton is aligned Axis dir = Axis(sinthetaprime*cos(phiprime), sinthetaprime*sin(phiprime), costhetaprime); switch (diffDirection){ case 0://Left single diffraction meMomenta()[4].setT(sqrt(mq2+sqr(meMomenta()[4].x())+sqr(meMomenta()[4].y())+sqr(meMomenta()[4].z()))); //////////////////////////////////////////////////// do{} while(!Kinematics::twoBodyDecay(p3,mqq(),mq(),-dir,decayMomenta.first,decayMomenta.second)); /////////// meMomenta()[2].setVect(p4.vect()); meMomenta()[2].setT(p4.t()); meMomenta()[3].setVect(decayMomenta.first.vect()); meMomenta()[3].setT(decayMomenta.first.t()); meMomenta()[4].setVect(decayMomenta.second.vect()); meMomenta()[4].setT(decayMomenta.second.t()); meMomenta()[2].rescaleEnergy(); meMomenta()[3].rescaleEnergy(); meMomenta()[4].rescaleEnergy(); break; case 1://Right single diffraction meMomenta()[4].setT(sqrt(mq2+sqr(meMomenta()[4].x())+sqr(meMomenta()[4].y())+sqr(meMomenta()[4].z()))); //////////////////////////////////////////////////// do{} while(!Kinematics::twoBodyDecay(p4,mqq(),mq(),dir,decayMomenta.first,decayMomenta.second)); meMomenta()[2].setVect(p3.vect()); meMomenta()[2].setT(p3.t()); meMomenta()[3].setVect(decayMomenta.first.vect()); meMomenta()[3].setT(decayMomenta.first.t()); meMomenta()[4].setVect(decayMomenta.second.vect()); meMomenta()[4].setT(decayMomenta.second.t()); meMomenta()[2].rescaleEnergy(); meMomenta()[3].rescaleEnergy(); meMomenta()[4].rescaleEnergy(); break; case 2://double diffraction do{} while(!Kinematics::twoBodyDecay(p3,mqq(),mq(),-dir,decayMomenta.first,decayMomenta.second)); do{} while(!Kinematics::twoBodyDecay(p4,mqq(),mq(),dir,decayMomentaTwo.first,decayMomentaTwo.second)); meMomenta()[2].setVect(decayMomenta.first.vect()); meMomenta()[2].setT(decayMomenta.first.t()); meMomenta()[3].setVect(decayMomenta.second.vect()); meMomenta()[3].setT(decayMomenta.second.t()); meMomenta()[4].setVect(decayMomentaTwo.second.vect()); meMomenta()[4].setT(decayMomentaTwo.second.t()); meMomenta()[5].setVect(decayMomentaTwo.first.vect()); meMomenta()[5].setT(decayMomentaTwo.first.t()); meMomenta()[2].rescaleEnergy(); meMomenta()[3].rescaleEnergy(); meMomenta()[4].rescaleEnergy(); meMomenta()[5].rescaleEnergy(); break; } }else { const auto tmp=diffDirection==1?1:0; meMomenta()[2+tmp].setVect(p3.vect()); meMomenta()[2+tmp].setT(p3.t()); meMomenta()[3-tmp].setVect(p4.vect()); meMomenta()[3-tmp].setT(p4.t()); meMomenta()[2].rescaleEnergy(); meMomenta()[3].rescaleEnergy(); } break; case 1: switch(diffDirection){ case 0: //quark and diquark masses meMomenta()[2].setMass(mqq()); meMomenta()[3].setMass(mq()); //gluon constituent mass meMomenta()[4].setMass(getParticleData(21)->constituentMass()); //outgoing proton meMomenta()[5].setVect(p4.vect()); meMomenta()[5].setT(p4.t()); //two body decay of the outgoing dissociation proton do{} while(!Kinematics::twoBodyDecay(p3,mqq()+mq(),getParticleData(21)->constituentMass(), p3.vect().unit(),momPair.first,momPair.second)); //put gluon back-to-back with quark-diquark //set momenta of quark and diquark frac = mqq()/(mqq()+mq()); meMomenta()[2].setVect(frac*momPair.first.vect()); meMomenta()[2].setT(sqrt(sqr(frac)*momPair.first.vect().mag2()+sqr(mqq()))); meMomenta()[3].setVect((1-frac)*momPair.first.vect()); meMomenta()[3].setT(sqrt(sqr(1-frac)*momPair.first.vect().mag2()+sqr(mq()))); //set momentum of gluon meMomenta()[4].setVect(momPair.second.vect()); meMomenta()[4].setT(momPair.second.t()); break; case 1: //quark and diquark masses meMomenta()[5].setMass(mqq()); meMomenta()[4].setMass(mq()); //gluon constituent mass meMomenta()[3].setMass(getParticleData(21)->constituentMass()); //outgoing proton meMomenta()[2].setVect(p3.vect()); meMomenta()[2].setT(p3.t()); //two body decay of the outgoing dissociation proton do{} while(!Kinematics::twoBodyDecay(p4,mqq()+mq(),getParticleData(21)->constituentMass(), p4.vect().unit(),momPair.first,momPair.second)); //put gluon back-to-back with quark-diquark //set momenta of quark and diquark frac = mqq()/(mqq()+mq()); meMomenta()[5].setVect(frac*momPair.first.vect()); meMomenta()[5].setT(sqrt(sqr(frac)*momPair.first.vect().mag2()+sqr(mqq()))); meMomenta()[4].setVect((1-frac)*momPair.first.vect()); meMomenta()[4].setT(sqrt(sqr(1-frac)*momPair.first.vect().mag2()+sqr(mq()))); //set momentum of gluon meMomenta()[3].setVect(momPair.second.vect()); meMomenta()[3].setT(momPair.second.t()); break; case 2: //first dissociated proton constituents meMomenta()[2].setMass(mqq()); meMomenta()[3].setMass(mq()); meMomenta()[4].setMass(getParticleData(21)->constituentMass()); //second dissociated proton constituents meMomenta()[5].setMass(getParticleData(21)->constituentMass()); meMomenta()[6].setMass(mq()); meMomenta()[7].setMass(mqq()); //two body decay of the outgoing dissociation proton do{} while(!Kinematics::twoBodyDecay(p3,mqq()+mq(),getParticleData(21)->constituentMass(), p3.vect().unit(),momPair.first,momPair.second)); do{} while(!Kinematics::twoBodyDecay(p4,mqq()+mq(),getParticleData(21)->constituentMass(), p4.vect().unit(),momPair1.first,momPair1.second)); //put gluon back-to-back with quark-diquark frac = mqq()/(mqq()+mq()); //first dissociated proton //set momenta of quark and diquark meMomenta()[2].setVect(frac*momPair.first.vect()); meMomenta()[2].setT(sqrt(sqr(frac)*momPair.first.vect().mag2()+sqr(mqq()))); meMomenta()[3].setVect((1-frac)*momPair.first.vect()); meMomenta()[3].setT(sqrt(sqr(1-frac)*momPair.first.vect().mag2()+sqr(mq()))); //set momentum of gluon meMomenta()[4].setVect(momPair.second.vect()); meMomenta()[4].setT(momPair.second.t()); //first dissociated proton //set momenta of quark and diquark meMomenta()[7].setVect(frac*momPair1.first.vect()); meMomenta()[7].setT(sqrt(sqr(frac)*momPair1.first.vect().mag2()+sqr(mqq()))); meMomenta()[6].setVect((1-frac)*momPair1.first.vect()); meMomenta()[6].setT(sqrt(sqr(1-frac)*momPair1.first.vect().mag2()+sqr(mq()))); //set momentum of gluon meMomenta()[5].setVect(momPair1.second.vect()); meMomenta()[5].setT(momPair1.second.t()); break; } meMomenta()[2].rescaleEnergy(); meMomenta()[3].rescaleEnergy(); meMomenta()[4].rescaleEnergy(); meMomenta()[5].rescaleEnergy(); if(diffDirection==2){ meMomenta()[6].rescaleEnergy(); meMomenta()[7].rescaleEnergy(); } break; } jacobian(sqr(cmEnergy)/GeV2); return true; } //Generate masses of dissociated protons and momentum transfer from probability f(M2,t) //(for single diffraction). Sample according to f(M2,t)=f(M2)f(t|M2). pair,Energy2> MEDiffraction::diffractiveMassAndMomentumTransfer() const { Energy2 theM12(ZERO),theM22(ZERO), thet(ZERO); //proton mass squared const Energy2 m2 = sqr(theProtonMass); //delta mass squared const Energy2 md2 = sqr(getParticleData(2214)->mass()); Energy2 M2; bool condition = true; do { //check if we want only delta if(deltaOnly) { switch(diffDirection){ case 0: theM12 = md2; theM22 = m2; M2 = md2; if(generator()->maximumCMEnergy()maximumCMEnergy()maximumCMEnergy()maximumCMEnergy()maximumCMEnergy()theM22) ? theM12: theM22; if(generator()->maximumCMEnergy()maximumCMEnergy(); const Energy2 s = sqr(cmEnergy); if(generator()->maximumCMEnergy()= sqr(cmEnergy)/softPomeronSlope()) { condition = true; continue; } const double expmax = exp(slope*tmaxfun(s,m2,M2)); const double expmin = exp(slope*tminfun(s,m2,M2)); // without (1-M2/s) constraint condition = (UseRandom::rnd() > protonPomeronSlope()*(expmax-expmin)/slope ); } while(condition); return make_pair (make_pair(theM12,theM22),thet); } //Decay of the excited proton to quark-diquark pair MEDiffraction::twoBodyDecayMomenta(Lorentz5Momentum pp) const{ //Decay of the excited proton const Energy2 Mx2(sqr(pp.mass())),mq2(sqr(mq())),mqq2(sqr(mqq())); const Energy2 psq = ((Mx2-sqr(mq()+mqq()))*(Mx2-sqr(mq()-mqq())))/(4*Mx2); assert(psq/GeV2>0); const Energy p(sqrt(psq)); const double phi = UseRandom::rnd() * Constants::twopi; const double costheta =1-2*UseRandom::rnd(); const double sintheta = sqrt(1-sqr(costheta)); Lorentz5Momentum k1=Lorentz5Momentum(p*sintheta*cos(phi), p*sintheta*sin(phi), p*costheta, sqrt(mq2+psq)); Lorentz5Momentum k2=Lorentz5Momentum(-p*sintheta*cos(phi), -p*sintheta*sin(phi), -p*costheta,sqrt(mqq2+psq)); //find boost to pp center of mass const Boost betap3 = (pp).findBoostToCM(); //k1 and k2 calculated at p3 center of mass, so boost back k1.boost(-betap3); k2.boost(-betap3); //first is quark, second diquark return make_pair(k1,k2); } Energy2 MEDiffraction::randomt(Energy2 M2) const { assert(protonPomeronSlope()*GeV2 > 0); //proton mass const Energy2 m2 = sqr( theProtonMass ); const Energy cmEnergy = generator()->maximumCMEnergy(); const Energy2 ttmin = tminfun(sqr(cmEnergy),m2,M2); const Energy2 ttmax = tmaxfun(sqr(cmEnergy),m2,M2); const InvEnergy2 slope = protonPomeronSlope() + 2*softPomeronSlope()*log(sqr(cmEnergy)/M2); double r = UseRandom::rnd(); Energy2 newVal; if(slope*ttmax>slope*ttmin) { newVal = ttmax + log( r + (1.-r)*exp(slope*(ttmin-ttmax)) ) / slope; } else { newVal = ttmin + log( 1. - r + r*exp(slope*(ttmax-ttmin))) / slope; } return newVal; } Energy2 MEDiffraction::doublediffrandomt(Energy2 M12, Energy2 M22) const { const Energy cmEnergy = generator()->maximumCMEnergy(); const double shift = 0.1; const InvEnergy2 slope = 2*softPomeronSlope()*log(shift+(sqr(cmEnergy)/softPomeronSlope())/(M12*M22)); const Energy2 ttmin = tminfun(sqr(cmEnergy),M12,M22); const Energy2 ttmax = tmaxfun(sqr(cmEnergy),M12,M22); double r = UseRandom::rnd(); Energy2 newVal; if(slope*ttmax>slope*ttmin) { newVal = ttmax + log( r + (1.-r)*exp(slope*(ttmin-ttmax)) ) / slope; } else { newVal = ttmin + log( 1. - r + r*exp(slope*(ttmax-ttmin))) / slope; } return newVal; } Energy2 MEDiffraction::randomM2() const { const double tmp = 1 - softPomeronIntercept(); const Energy cmEnergy = generator()->maximumCMEnergy(); return sqr(cmEnergy) * pow( pow(M2min()/sqr(cmEnergy),tmp) + UseRandom::rnd() * (pow(M2max()/sqr(cmEnergy),tmp) - pow(M2min()/sqr(cmEnergy),tmp)), 1.0/tmp ); } Energy2 MEDiffraction::tminfun(Energy2 s, Energy2 M12, Energy2 M22) const { const Energy2 m2 = sqr( theProtonMass ); return 0.5/s*(-sqrt(kallen(s, m2, m2)*kallen(s, M12, M22))-sqr(s)+2*s*m2+s*M12+s*M22); } Energy2 MEDiffraction::tmaxfun(Energy2 s, Energy2 M12, Energy2 M22) const { const Energy2 m2 = sqr( theProtonMass ); return 0.5/s*(sqrt(kallen(s, m2, m2)*kallen(s, M12, M22))-sqr(s)+2*s*m2+s*M12+s*M22); } + +double MEDiffraction::correctionweight() const { + + + // Here we calculate the weight to restore the diffractiveXSec + // given by the MPIHandler. + + // First get the eventhandler to get the current cross sections. + static Ptr::tptr eh = + dynamic_ptr_cast::tptr>(generator()->eventHandler()); + + // All diffractive processes make use of this ME. + // The static map can be used to collect all the sumOfWeights. + static map weightsmap; + weightsmap[lastXCombPtr()]=lastXComb().stats().sumWeights(); + + + // Define static variable to keep trac of reweighting + static double rew_=1.; + static int countUpdateWeight=50; + static double sumRew=0.; + static double countN=0; + + // if we produce events we count + if(eh->integratedXSec()>ZERO)sumRew+=rew_; + if(eh->integratedXSec()>ZERO)countN+=1.; + + + + if(countUpdateWeightsampler()->maxXSec()/eh->sampler()->attempts()*sum; + CrossSection XS_wanted=MPIHandler_->diffractiveXSec(); + double deltaN=50; + rew_=avcsNorm2*(XS_wanted*(countN+deltaN)-XS_have*countN)/(XS_have*deltaN); + countUpdateWeight+=deltaN; + } + //Make sure we dont produce negative weights. + // TODO: write finalize method that checks if reweighting was performed correctly. + rew_=max(rew_,0.000001); + rew_=min(rew_,10000.0); + + return rew_; + +} + double MEDiffraction::me2() const{ - return theme2; + return theme2; } CrossSection MEDiffraction::dSigHatDR() const { - return me2()*jacobian()/sHat()*sqr(hbarc); + return me2()*jacobian()/sHat()*sqr(hbarc)*correctionweight(); } unsigned int MEDiffraction::orderInAlphaS() const { return 0; } unsigned int MEDiffraction::orderInAlphaEW() const { return 0; } Selector MEDiffraction::diagrams(const DiagramVector & diags) const { Selector sel; if(!deltaOnly){ if(diffDirection<2){ for(unsigned int i = 0; i < diags.size(); i++){ if(diags[0]->id()==-1) sel.insert(2./3.,i); else sel.insert(1./3.,i); } }else{ for(unsigned int i = 0; i < diags.size(); i++){ if(diags[0]->id()==-1) sel.insert(4./9.,i); else if(diags[0]->id()==-2) sel.insert(2./9.,i); else if(diags[0]->id()==-3) sel.insert(2./9.,i); else sel.insert(1./9.,i); } } }else{ sel.insert(1.0,0); } return sel; } Selector MEDiffraction::colourGeometries(tcDiagPtr ) const { Selector sel; int sign1=0, sign2=0; sign1 = (generator()->eventHandler()->incoming().first->id() > 0) ? 1 : -1; sign2 = (generator()->eventHandler()->incoming().second->id() > 0) ? 1 : -1; switch(dissociationDecay){ case 0: if(!deltaOnly) { if(diffDirection!=2){ if (diffDirection == 0){ if(sign1>0){ static ColourLines dqq0=ColourLines("-6 2 7"); sel.insert(1.0,&dqq0); }else{ static ColourLines dqq0=ColourLines("6 -2 -7"); sel.insert(1.0,&dqq0); } } else{ if(sign2>0){ static ColourLines dqq1=ColourLines("-6 3 7"); sel.insert(1.0,&dqq1); }else{ static ColourLines dqq1=ColourLines("6 -3 -7"); sel.insert(1.0,&dqq1); } } }else{ if(sign1>0 && sign2>0){ static ColourLines ddqq0=ColourLines("-6 2 7, -9 4 8"); sel.insert(1.0,&ddqq0); }else if(sign1<0 && sign2>0){ static ColourLines ddqq0=ColourLines("6 -2 -7, -9 4 8"); sel.insert(1.0,&ddqq0); }else if(sign1>0&& sign2<0){ static ColourLines ddqq0=ColourLines("-6 2 7, 9 -4 -8"); sel.insert(1.0,&ddqq0); }else{ static ColourLines ddqq0=ColourLines("6 -2 -7, 9 -4 -8"); sel.insert(1.0,&ddqq0); } } }else { static ColourLines cl(""); sel.insert(1.0, &cl); } break; case 1: switch(diffDirection){ case 0: static ColourLines clleft("-6 2 3 8, -8 -3 7"); sel.insert(1.0, &clleft); break; case 1: static ColourLines clright("-9 4 3 7, -7 -3 8"); sel.insert(1.0, &clright); break; case 2: static ColourLines cldouble("-8 2 3 10, -10 -3 9, -13 6 5 11, -11 -5 12"); sel.insert(1.0, &cldouble); break; } break; } return sel; } void MEDiffraction::doinit() { HwMEBase::doinit(); theProtonMass = getParticleData(2212)->mass(); } void MEDiffraction::doinitrun() { HwMEBase::doinitrun(); isInRunPhase = true; } IBPtr MEDiffraction::clone() const { return new_ptr(*this); } IBPtr MEDiffraction::fullclone() const { return new_ptr(*this); } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigMEDiffraction("Herwig::MEDiffraction", "HwMEHadron.so"); void MEDiffraction::persistentOutput(PersistentOStream & os) const { os << theme2 << deltaOnly << diffDirection << theprotonPomeronSlope << thesoftPomeronIntercept << thesoftPomeronSlope << dissociationDecay - << ounit(theProtonMass,GeV); + << ounit(theProtonMass,GeV) << MPIHandler_; } void MEDiffraction::persistentInput(PersistentIStream & is, int) { is >> theme2 >> deltaOnly >> diffDirection >> theprotonPomeronSlope >> thesoftPomeronIntercept >> thesoftPomeronSlope >> dissociationDecay - >> iunit(theProtonMass,GeV); + >> iunit(theProtonMass,GeV)>> MPIHandler_; } InvEnergy2 MEDiffraction::protonPomeronSlope() const{ return theprotonPomeronSlope/GeV2; } double MEDiffraction::softPomeronIntercept() const { return thesoftPomeronIntercept; } InvEnergy2 MEDiffraction::softPomeronSlope() const { return thesoftPomeronSlope/GeV2; } void MEDiffraction::Init() { static ClassDocumentation documentation ("There is no documentation for the MEDiffraction class"); static Parameter interfaceme2 ("DiffractionAmplitude", "The square of the diffraction amplitude used to determine the " "cross section.", &MEDiffraction::theme2, 1.0, 0.00001, 100.0, false, false, Interface::limited); static Parameter interfaceprotonPomeronSlope ("ProtonPomeronSlope", "The proton-pomeron slope parameter.", &MEDiffraction::theprotonPomeronSlope, 10.1, 0.00001, 100.0, false, false, Interface::limited); static Parameter interfacesoftPomeronIntercept ("SoftPomeronIntercept", "The soft pomeron intercept.", &MEDiffraction::thesoftPomeronIntercept, 1.08, 0.00001, 100.0, false, false, Interface::limited); static Parameter interfacesoftPomeronSlope ("SoftPomeronSlope", "The soft pomeron slope parameter.", &MEDiffraction::thesoftPomeronSlope, 0.25, 0.00001, 100.0, false, false, Interface::limited); static Switch interfaceDeltaOnly ("DeltaOnly", "proton-proton to proton-delta only", &MEDiffraction::deltaOnly, 0, false, false); static SwitchOption interfaceDeltaOnly0 (interfaceDeltaOnly,"No","Final state with Delta only is OFF", 0); static SwitchOption interfaceDeltaOnly1 (interfaceDeltaOnly,"Yes","Final state with Delta only is ON", 1); //Select if the left, right or both protons are excited static Switch interfaceDiffDirection ("DiffDirection", "Direction of the excited proton", &MEDiffraction::diffDirection, 0, false, false); static SwitchOption left (interfaceDiffDirection,"Left","Proton moving in the positive z direction", 0); static SwitchOption right (interfaceDiffDirection,"Right","Proton moving in the negative z direction", 1); static SwitchOption both (interfaceDiffDirection,"Both","Both protons", 2); //Select if two or three body decay static Switch interfaceDissociationDecay ("DissociationDecay", "Number of clusters the dissociated proton decays", &MEDiffraction::dissociationDecay, 0, false, false); static SwitchOption one (interfaceDissociationDecay,"One","Dissociated proton decays into one cluster", 0); static SwitchOption two (interfaceDissociationDecay,"Two","Dissociated proton decays into two clusters", 1); + + + static Reference interfaceMPIHandler + ("MPIHandler", + "The object that administers all additional scatterings.", + &MEDiffraction::MPIHandler_, false, false, true, true); + + + + } diff --git a/MatrixElement/Hadron/MEDiffraction.h b/MatrixElement/Hadron/MEDiffraction.h --- a/MatrixElement/Hadron/MEDiffraction.h +++ b/MatrixElement/Hadron/MEDiffraction.h @@ -1,302 +1,315 @@ // -*- C++ -*- #ifndef HERWIG_MEDiffraction_H #define HERWIG_MEDiffraction_H // // This is the declaration of the MEDiffraction class. // #include "Herwig/MatrixElement/HwMEBase.h" +#include "Herwig/Shower/UEBase.h" namespace Herwig { using namespace ThePEG; /** * The MEDiffraction class provides a simple colour singlet exchange matrix element * to be used in the soft component of the multiple scattering model of the * underlying event * * @see \ref MEDiffractionInterfaces "The interfaces" * defined for MEDiffraction. */ class MEDiffraction: public HwMEBase { public: MEDiffraction(); /** @name Virtual functions required by the MEBase class. */ //@{ /** * Return the order in \f$\alpha_S\f$ in which this matrix * element is given. */ virtual unsigned int orderInAlphaS() const; /** * Return the order in \f$\alpha_{EW}\f$ in which this matrix * element is given. */ virtual unsigned int orderInAlphaEW() const; /** * The matrix element for the kinematical configuration * previously provided by the last call to setKinematics(), suitably * scaled by sHat() to give a dimension-less number. * @return the matrix element scaled with sHat() to give a * dimensionless number. */ virtual double me2() const; + + /** + * Correction weight to reweight the cross section to the diffractive + * cross section. + */ + double correctionweight() const; /** * Return the scale associated with the last set phase space point. */ virtual Energy2 scale() const; /** * Set the typed and momenta of the incoming and outgoing partons to * be used in subsequent calls to me() and colourGeometries() * according to the associated XComb object. If the function is * overridden in a sub class the new function must call the base * class one first. */ virtual void setKinematics(); /** * The number of internal degrees of freedom used in the matrix * element. */ virtual int nDim() const; /** * Generate internal degrees of freedom given nDim() uniform * random numbers in the interval \f$ ]0,1[ \f$. To help the phase space * generator, the dSigHatDR should be a smooth function of these * numbers, although this is not strictly necessary. * @param r a pointer to the first of nDim() consecutive random numbers. * @return true if the generation succeeded, otherwise false. */ virtual bool generateKinematics(const double * r); /** * Return the matrix element squared differential in the variables * given by the last call to generateKinematics(). */ virtual CrossSection dSigHatDR() const; /** * Add all possible diagrams with the add() function. */ virtual void getDiagrams() const; /** * Get diagram selector. With the information previously supplied with the * setKinematics method, a derived class may optionally * override this method to weight the given diagrams with their * (although certainly not physical) relative probabilities. * @param dv the diagrams to be weighted. * @return a Selector relating the given diagrams to their weights. */ virtual Selector diagrams(const DiagramVector & dv) const; /** * Return a Selector with possible colour geometries for the selected * diagram weighted by their relative probabilities. * @param diag the diagram chosen. * @return the possible colour geometries weighted by their * relative probabilities. */ virtual Selector colourGeometries(tcDiagPtr diag) const; //@} /** * Expect the incoming partons in the laboratory frame */ /* virtual bool wantCMS() const { return false; } */ public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); /** * Initialize this object. Called in the run phase just before a run begins. */ virtual void doinitrun(); //@} /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} private: /* The matrix element squared */ double theme2; /* Use only delta as excited state */ bool deltaOnly; /* Direction of the excited proton */ unsigned int diffDirection; /* Number of clusters the dissociated proton decays into */ unsigned int dissociationDecay; /* The mass of the consitutent quark */ Energy mq() const {return Energy(0.325*GeV);} /* The mass of the constituent diquark */ Energy mqq() const {return Energy(0.650*GeV);} /* The proton-pomeron slope */ double theprotonPomeronSlope; /* The soft pomeron intercept */ double thesoftPomeronIntercept; /* The soft pomeron slope */ double thesoftPomeronSlope; /** * Sample the diffractive mass squared M2 and the momentum transfer t */ pair,Energy2> diffractiveMassAndMomentumTransfer() const; /** * Random value for the diffractive mass squared M2 according to (M2/s0)^(-intercept) */ Energy2 randomM2() const; /** * Random value for t according to exp(diffSlope*t) */ Energy2 randomt(Energy2 M2) const; /** * Random value for t according to exp(diffSlope*t) for double diffraction */ Energy2 doublediffrandomt(Energy2 M12, Energy2 M22) const; /** * Returns the momenta of the two-body decay of momentum pp */ pair twoBodyDecayMomenta(Lorentz5Momentum pp) const; /** * Returns the proton-pomeron slope */ InvEnergy2 protonPomeronSlope() const; /** * Returns the soft pomeron intercept */ double softPomeronIntercept() const; //M12 and M22 are masses squared of //outgoing particles /** * Returns the minimal possible value of momentum transfer t given the center * of mass energy and diffractive masses */ Energy2 tminfun(Energy2 s, Energy2 M12, Energy2 M22) const; /** * Returns the maximal possible value of momentum transfer t given the center * of mass energy and diffractive masses */ Energy2 tmaxfun(Energy2 s , Energy2 M12, Energy2 M22) const; /** * Returns the minimal possible value of diffractive mass */ //lowest possible mass given the constituent masses of quark and diquark Energy2 M2min() const{return sqr(getParticleData(2212)->mass()+mq()+mqq());} /** * Returns the maximal possible value of diffractive mass */ Energy2 M2max() const{ return sqr(generator()->maximumCMEnergy()-getParticleData(2212)->mass()); }//TODO:modify to get proper parameters InvEnergy2 softPomeronSlope() const; /* Kallen function */ template auto kallen(A a, B b, C c) const -> decltype(a*a) { return a*a + b*b + c*c - 2.0*(a*b + b*c + c*a); } /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ MEDiffraction & operator=(const MEDiffraction &) = delete; bool isInRunPhase; /* The proton mass */ Energy theProtonMass; + + /** + * a MPIHandler to administer the creation of several (semihard) + * partonic interactions. + */ + UEBasePtr MPIHandler_; }; } #endif /* HERWIG_MEDiffraction_H */ diff --git a/MatrixElement/Hadron/MEMinBias.cc b/MatrixElement/Hadron/MEMinBias.cc --- a/MatrixElement/Hadron/MEMinBias.cc +++ b/MatrixElement/Hadron/MEMinBias.cc @@ -1,174 +1,244 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the MEMinBias class. // #include "MEMinBias.h" #include "ThePEG/Utilities/DescribeClass.h" +#include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Utilities/SimplePhaseSpace.h" //#include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Handlers/StandardXComb.h" #include "ThePEG/Interface/Parameter.h" +#include "ThePEG/Handlers/SamplerBase.h" + #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace Herwig; #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/MatrixElement/Tree2toNDiagram.h" void MEMinBias::getDiagrams() const { int maxflav(2); // Pomeron data tcPDPtr pom = getParticleData(990); for ( int i = 1; i <= maxflav; ++i ) { for( int j=1; j <= i; ++j){ tcPDPtr q1 = getParticleData(i); tcPDPtr q1b = q1->CC(); tcPDPtr q2 = getParticleData(j); tcPDPtr q2b = q2->CC(); // For each flavour we add: //qq -> qq add(new_ptr((Tree2toNDiagram(3), q1, pom, q2, 1, q1, 2, q2, -1))); //qqb -> qqb add(new_ptr((Tree2toNDiagram(3), q1, pom, q2b, 1, q1, 2, q2b, -2))); //qbqb -> qbqb add(new_ptr((Tree2toNDiagram(3), q1b, pom, q2b, 1, q1b, 2, q2b, -3))); } } } Energy2 MEMinBias::scale() const { return sqr(Scale_); } int MEMinBias::nDim() const { return 0; } void MEMinBias::setKinematics() { HwMEBase::setKinematics(); // Always call the base class method first. } bool MEMinBias::generateKinematics(const double *) { // generate the masses of the particles for ( int i = 2, N = meMomenta().size(); i < N; ++i ) { meMomenta()[i] = Lorentz5Momentum(mePartonData()[i]->generateMass()); } Energy q = ZERO; try { q = SimplePhaseSpace:: getMagnitude(sHat(), meMomenta()[2].mass(), meMomenta()[3].mass()); } catch ( ImpossibleKinematics & e ) { return false; } Energy pt = ZERO; meMomenta()[2].setVect(Momentum3( pt, pt, q)); meMomenta()[3].setVect(Momentum3(-pt, -pt, -q)); meMomenta()[2].rescaleEnergy(); meMomenta()[3].rescaleEnergy(); jacobian(1.0); return true; } + +double MEMinBias::correctionweight() const { + + + + // Here we calculate the weight to restore the inelastic-diffractiveXSec + // given by the MPIHandler. + + // First get the eventhandler to get the current cross sections. + static Ptr::tptr eh = + dynamic_ptr_cast::tptr>(generator()->eventHandler()); + + // All diffractive processes make use of this ME. + // The static map can be used to collect all the sumOfWeights. + static map weightsmap; + weightsmap[lastXCombPtr()]=lastXComb().stats().sumWeights(); + + + // Define static variable to keep trac of reweighting + static double rew_=1.; + static int countUpdateWeight=50; + static double sumRew=0.; + static double countN=0; + + // if we produce events we count + if(eh->integratedXSec()>ZERO)sumRew+=rew_; + if(eh->integratedXSec()>ZERO)countN+=1.; + + + + if(countUpdateWeightsampler()->maxXSec()/eh->sampler()->attempts()*sum; + CrossSection XS_wanted=MPIHandler_->inelasticXSec()-MPIHandler_->diffractiveXSec(); + double deltaN=50; + rew_=avcsNorm2*(XS_wanted*(countN+deltaN)-XS_have*countN)/(XS_have*deltaN); + countUpdateWeight+=deltaN; + } + //Make sure we dont produce negative weights. + // TODO: write finalize method that checks if reweighting was performed correctly. + rew_=max(rew_,0.000001); + rew_=min(rew_,10000.0); + + return rew_; + + + + + + + +} + + + double MEMinBias::me2() const { //tuned so it gives the correct normalization for xmin = 0.11 return csNorm_*(sqr(generator()->maximumCMEnergy())/GeV2); } CrossSection MEMinBias::dSigHatDR() const { - return me2()*jacobian()/sHat()*sqr(hbarc); + return me2()*jacobian()/sHat()*sqr(hbarc)*correctionweight(); } unsigned int MEMinBias::orderInAlphaS() const { return 2; } unsigned int MEMinBias::orderInAlphaEW() const { return 0; } Selector MEMinBias::diagrams(const DiagramVector & diags) const { Selector sel; for ( DiagramIndex i = 0; i < diags.size(); ++i ) sel.insert(1.0, i); return sel; } Selector MEMinBias::colourGeometries(tcDiagPtr diag) const { static ColourLines qq("1 4, 3 5"); static ColourLines qqb("1 4, -3 -5"); static ColourLines qbqb("-1 -4, -3 -5"); Selector sel; switch(diag->id()){ case -1: sel.insert(1.0, &qq); break; case -2: sel.insert(1.0, &qqb); break; case -3: sel.insert(1.0, &qbqb); break; } return sel; } IBPtr MEMinBias::clone() const { return new_ptr(*this); } IBPtr MEMinBias::fullclone() const { return new_ptr(*this); } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigMEMinBias("Herwig::MEMinBias", "HwMEHadron.so"); void MEMinBias::persistentOutput(PersistentOStream & os) const { - os << csNorm_ << ounit(Scale_,GeV);; + os << csNorm_ << ounit(Scale_,GeV) << MPIHandler_; } void MEMinBias::persistentInput(PersistentIStream & is, int) { - is >> csNorm_ >> iunit(Scale_,GeV); + is >> csNorm_ >> iunit(Scale_,GeV) >> MPIHandler_; } void MEMinBias::Init() { static ClassDocumentation documentation ("There is no documentation for the MEMinBias class"); static Parameter interfacecsNorm ("csNorm", "Normalization of the min-bias cross section.", &MEMinBias::csNorm_, 1.0, 0.0, 100.0, false, false, Interface::limited); static Parameter interfaceScale ("Scale", "Scale for the Min Bias matrix element.", &MEMinBias::Scale_,GeV, 2.0*GeV, 0.0*GeV, 100.0*GeV, false, false, Interface::limited); + static Reference interfaceMPIHandler + ("MPIHandler", + "The object that administers all additional scatterings.", + &MEMinBias::MPIHandler_, false, false, true, true); + + } diff --git a/MatrixElement/Hadron/MEMinBias.h b/MatrixElement/Hadron/MEMinBias.h --- a/MatrixElement/Hadron/MEMinBias.h +++ b/MatrixElement/Hadron/MEMinBias.h @@ -1,195 +1,212 @@ // -*- C++ -*- #ifndef HERWIG_MEMinBias_H #define HERWIG_MEMinBias_H // // This is the declaration of the MEMinBias class. // #include "Herwig/MatrixElement/HwMEBase.h" +#include "Herwig/Shower/UEBase.h" namespace Herwig { using namespace ThePEG; /** * The MEMinBias class provides a simple colour singlet exchange matrix element * to be used in the soft component of the multiple scattering model of the * underlying event * * @see \ref MEMinBiasInterfaces "The interfaces" * defined for MEMinBias. */ class MEMinBias: public HwMEBase { public: /** * The default constructor. */ MEMinBias() : csNorm_(1.), Scale_(2.*GeV) {} public: /** @name Virtual functions required by the MEBase class. */ //@{ /** * Return the order in \f$\alpha_S\f$ in which this matrix * element is given. */ virtual unsigned int orderInAlphaS() const; /** * Return the order in \f$\alpha_{EW}\f$ in which this matrix * element is given. */ virtual unsigned int orderInAlphaEW() const; /** * The matrix element for the kinematical configuration * previously provided by the last call to setKinematics(), suitably * scaled by sHat() to give a dimension-less number. * @return the matrix element scaled with sHat() to give a * dimensionless number. */ virtual double me2() const; /** + * Correction weight to reweight the cross section to the inelastic cross + * section subtracted by the diffractive cross section. + */ + double correctionweight() const; + + /** * Return the scale associated with the last set phase space point. */ virtual Energy2 scale() const; /** * Set the typed and momenta of the incoming and outgoing partons to * be used in subsequent calls to me() and colourGeometries() * according to the associated XComb object. If the function is * overridden in a sub class the new function must call the base * class one first. */ virtual void setKinematics(); /** * The number of internal degrees of freedom used in the matrix * element. */ virtual int nDim() const; /** * Generate internal degrees of freedom given nDim() uniform * random numbers in the interval \f$ ]0,1[ \f$. To help the phase space * generator, the dSigHatDR should be a smooth function of these * numbers, although this is not strictly necessary. * @param r a pointer to the first of nDim() consecutive random numbers. * @return true if the generation succeeded, otherwise false. */ virtual bool generateKinematics(const double * r); /** * Return the matrix element squared differential in the variables * given by the last call to generateKinematics(). */ virtual CrossSection dSigHatDR() const; /** * Add all possible diagrams with the add() function. */ virtual void getDiagrams() const; /** * Get diagram selector. With the information previously supplied with the * setKinematics method, a derived class may optionally * override this method to weight the given diagrams with their * (although certainly not physical) relative probabilities. * @param dv the diagrams to be weighted. * @return a Selector relating the given diagrams to their weights. */ virtual Selector diagrams(const DiagramVector & dv) const; /** * Return a Selector with possible colour geometries for the selected * diagram weighted by their relative probabilities. * @param diag the diagram chosen. * @return the possible colour geometries weighted by their * relative probabilities. */ virtual Selector colourGeometries(tcDiagPtr diag) const; //@} public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} // If needed, insert declarations of virtual function defined in the // InterfacedBase class here (using ThePEG-interfaced-decl in Emacs). private: /** - * Normalization of the min-bias cross section - */ + * Normalization of the min-bias cross section. + * Note that the cross section is reweighted in addition to produce the + * non-diffractive cross section given by the MPIHandler + * csNorm can be modified to improve the unweighting effiency. + */ double csNorm_; /** * Scale for the Min Bias matrix element */ Energy Scale_; /** + * a MPIHandler to administer the creation of several (semihard) + * partonic interactions. + * Needed to comunicate the non-diffractive cross section. + */ + UEBasePtr MPIHandler_; + + /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ MEMinBias & operator=(const MEMinBias &) = delete; }; } #endif /* HERWIG_MEMinBias_H */ diff --git a/Shower/UEBase.h b/Shower/UEBase.h --- a/Shower/UEBase.h +++ b/Shower/UEBase.h @@ -1,131 +1,141 @@ // -*- C++ -*- #ifndef HERWIG_UEBase_H #define HERWIG_UEBase_H // // This is the declaration of the UEBase class. // #include "ThePEG/Interface/Interfaced.h" #include "ThePEG/Handlers/StandardXComb.fh" #include "UEBase.fh" namespace Herwig { using namespace ThePEG; /** * Abstract base class used to minimize the dependence between * MPIHandler and all Shower classes. * * \author Manuel B\"ahr * * @see \ref UEBaseInterfaces "The interfaces" * defined for UEBase. */ class UEBase: public Interfaced { public: /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); /** @name Virtual functions used for the generation of additional interactions . */ //@{ /** * Some initialization code eventually. */ virtual void initialize() {} /** * Return true or false depending on the generator setup. */ virtual bool beamOK() const = 0; /** * Return true or false depending on whether soft interactions are enabled. */ virtual bool softInt() const {return false;} /** * Return the value of the pt cutoff. */ virtual Energy Ptmin() const = 0; /** * Return the slope of the soft pt spectrum. Only necessary when the * soft part is modelled. */ virtual InvEnergy2 beta() const {return ZERO;} /** * Some finalize code eventually. */ virtual void finalize() {} /** * Clean up method called after each event. */ virtual void clean() {} /** * Return the number of different hard processes. Use 0 as default to * not require implementation. */ virtual unsigned int additionalHardProcs() const {return 0;} /** * return the hard multiplicity of process i. Can't be constant in my * case because drawing from the probability distribution also * specifies the soft multiplicity that has to be stored.... */ virtual unsigned int multiplicity(unsigned int i=0) = 0; /** * Generate a additional interaction for ProcessHandler sel. Method * can't be const because it saves the state of the underlying XComb * object on it's way. */ virtual tStdXCombPtr generate(unsigned int sel=0) = 0; /** * Return the type of algorithm. */ virtual int Algorithm() const = 0; /** * Return the value of the hard Process pt cutoff for vetoing. */ virtual Energy PtForVeto() const = 0; /** * Return the fraction of colour disrupted subprocesses. Use default 0 * so that it is not required to implement. */ virtual double colourDisrupt() const {return 0.0;} /** * Return the soft multiplicity. Use 0 as default to not require * implementation. */ virtual unsigned int softMultiplicity() const {return 0;} //@} + /** + * Return the inelastic cross section ( sigmaND + sigmaDiff ) + */ + virtual CrossSection inelasticXSec() const =0; + + /** + * Return the diffractiv cross section assumed by the model. + */ + virtual CrossSection diffractiveXSec() const =0; + private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ UEBase & operator=(const UEBase &) = delete; }; } #endif /* HERWIG_UEBase_H */ diff --git a/UnderlyingEvent/MPIHandler.cc b/UnderlyingEvent/MPIHandler.cc --- a/UnderlyingEvent/MPIHandler.cc +++ b/UnderlyingEvent/MPIHandler.cc @@ -1,857 +1,874 @@ // -*- C++ -*- // // MPIHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the MPIHandler class. // #include "MPIHandler.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Handlers/StandardXComb.h" #include "ThePEG/Handlers/SubProcessHandler.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/MatrixElement/MEBase.h" #include "ThePEG/Handlers/CascadeHandler.h" #include "ThePEG/Cuts/Cuts.h" #include "ThePEG/Cuts/JetCuts.h" #include "ThePEG/Cuts/SimpleKTCut.h" #include "ThePEG/PDF/PartonExtractor.h" #include "gsl/gsl_sf_bessel.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace Herwig; MPIHandler * MPIHandler::currentHandler_ = 0; bool MPIHandler::beamOK() const { return (HadronMatcher::Check(*eventHandler()->incoming().first) && HadronMatcher::Check(*eventHandler()->incoming().second) ); } tStdXCombPtr MPIHandler::generate(unsigned int sel) { //generate a certain process if(sel+1 > processHandlers().size()) throw Exception() << "MPIHandler::generate called with argument out of range" << Exception::runerror; return processHandlers()[sel]->generate(); } IBPtr MPIHandler::clone() const { return new_ptr(*this); } IBPtr MPIHandler::fullclone() const { return new_ptr(*this); } void MPIHandler::finalize() { if( beamOK() ){ statistics(); } } void MPIHandler::initialize() { currentHandler_ = this; useMe(); theHandler = generator()->currentEventHandler(); //stop if the EventHandler is not present: assert(theHandler); //check if MPI is wanted if( !beamOK() ){ throw Exception() << "You have requested multiple parton-parton scattering,\n" << "but the model is not forseen for the beam setup you chose.\n" << "You should therefore disable that by setting XXXGenerator:EventHandler:" << "CascadeHandler:MPIHandler to NULL" << Exception::runerror; } numSubProcs_ = subProcesses().size(); if( numSubProcs_ != cuts().size() ) throw Exception() << "MPIHandler::each SubProcess needs a Cuts Object" << "ReferenceVectors are not equal in size" << Exception::runerror; if( additionalMultiplicities_.size()+1 != numSubProcs_ ) throw Exception() << "MPIHandler: each additional SubProcess needs " << "a multiplicity assigned. This can be done in with " << "insert MPIHandler:additionalMultiplicities 0 1" << Exception::runerror; //identicalToUE_ = 0 hard process is identical to ue, -1 no one if( identicalToUE_ > (int)numSubProcs_ || identicalToUE_ < -1 ) throw Exception() << "MPIHandler:identicalToUE has disallowed value" << Exception::runerror; // override the cuts for the additional scatters if energyExtrapolation_ is // set if (energyExtrapolation_ != 0 ) { overrideUECuts(); } tcPDPtr gluon=getParticleData(ParticleID::g); //determine ptmin Ptmin_ = cuts()[0]->minKT(gluon); if(identicalToUE_ == -1){ algorithm_ = 2; }else{ if(identicalToUE_ == 0){ //Need to work a bit, in case of LesHouches events for QCD2to2 Ptr::pointer eH = dynamic_ptr_cast::pointer>(eventHandler()); if( eH ) { PtOfQCDProc_ = eH->cuts()->minKT(gluon); // find the jet cut in the new style cuts for(unsigned int ix=0;ixcuts()->multiCuts().size();++ix) { Ptr::pointer jetCuts = dynamic_ptr_cast::pointer>(eH->cuts()->multiCuts()[ix]); if(jetCuts) { Energy ptMin=1e30*GeV; for(unsigned int iy=0;iyjetRegions().size();++iy) { ptMin = min(ptMin,jetCuts->jetRegions()[iy]->ptMin()); } if(ptMin<1e29*GeV&&ptMin>PtOfQCDProc_) PtOfQCDProc_ = ptMin; } } } else { if(PtOfQCDProc_ == -1.0*GeV) throw Exception() << "MPIHandler: You need to specify the pt cutoff " << "used to in the LesHouches file for QCD2to2 events" << Exception::runerror; } } else { PtOfQCDProc_ = cuts()[identicalToUE_]->minKT(gluon); } if(PtOfQCDProc_ > 2*Ptmin_) algorithm_ = 1; else algorithm_ = 0; if(PtOfQCDProc_ == ZERO)//pure MinBias mode algorithm_ = -1; } //Init all subprocesses for(unsigned int i=0; iinitialize(subProcesses()[i], cuts()[i], eventHandler()); processHandlers().back()->initrun(); } //now calculate the individual Probabilities XSVector UEXSecs; UEXSecs.push_back(processHandlers()[0]->integratedXSec()); //save the hard cross section hardXSec_ = UEXSecs.front(); //determine sigma_soft and beta if(softInt_){//check that soft ints are requested GSLBisection rootFinder; if(twoComp_){ //two component model /* GSLMultiRoot eqSolver; slopeAndTotalXSec eq(this); pair res = eqSolver.value(eq, 10*millibarn, 0.6*GeV2); softXSec_ = res.first; softMu2_ = res.second; */ slopeBisection fs(this); try{ softMu2_ = rootFinder.value(fs, 0.3*GeV2, 1.*GeV2); softXSec_ = fs.softXSec(); }catch(GSLBisection::IntervalError){ try{ // Very low energies (e.g. 200 GeV) need this. // Very high energies (e.g. 100 TeV) produce issues with // this choice. This is a temp. fix. softMu2_ = rootFinder.value(fs, 0.25*GeV2, 1.3*GeV2); softXSec_ = fs.softXSec(); }catch(GSLBisection::IntervalError){ throw Exception() << "\n**********************************************************\n" "* Inconsistent MPI parameter choice for this beam energy *\n" "**********************************************************\n" "MPIHandler parameter choice is unable to reproduce\n" "the total cross section. Please check arXiv:0806.2949\n" "for the allowed parameter space." << Exception::runerror; } } }else{ //single component model TotalXSecBisection fn(this); try{ softXSec_ = rootFinder.value(fn, 0*millibarn, 5000*millibarn); }catch(GSLBisection::IntervalError){ throw Exception() << "\n**********************************************************\n" "* Inconsistent MPI parameter choice for this beam energy *\n" "**********************************************************\n" "MPIHandler parameter choice is unable to reproduce\n" "the total cross section. Please check arXiv:0806.2949\n" "for the allowed parameter space." << Exception::runerror; } } //now get the differential cross section at ptmin ProHdlPtr qcd = new_ptr(ProcessHandler()); Energy eps = 0.1*GeV; Energy ptminPlus = Ptmin_ + eps; Ptr::pointer ktCut = new_ptr(SimpleKTCut(ptminPlus)); ktCut->init(); ktCut->initrun(); CutsPtr qcdCut = new_ptr(Cuts(2*ptminPlus)); qcdCut->add(dynamic_ptr_cast(ktCut)); qcdCut->init(); qcdCut->initrun(); qcd->initialize(subProcesses()[0], qcdCut, eventHandler()); qcd->initrun(); // ds/dp_T^2 = 1/2/p_T ds/dp_T DiffXSec hardPlus = (hardXSec_-qcd->integratedXSec())/(2*Ptmin_*eps); betaBisection fn2(softXSec_, hardPlus, Ptmin_); try{ beta_ = rootFinder.value(fn2, -10/GeV2, 2/GeV2); }catch(GSLBisection::IntervalError){ try{ // Very low energies (e.g. 200 GeV) need this. // Very high energies (e.g. 100 TeV) produce issues with // this choice. This is a temp. fix. beta_ = rootFinder.value(fn2, -5/GeV2, 8./GeV2); }catch(GSLBisection::IntervalError){ throw Exception() << "MPIHandler: slope of soft pt spectrum couldn't be " << "determined." << Exception::runerror; } } } Probs(UEXSecs); //MultDistribution("probs.test"); UEXSecs.clear(); } void MPIHandler::MultDistribution(string filename) const { ofstream file; double p(0.0), pold(0.0); file.open(filename.c_str()); //theMultiplicities Selector::const_iterator it = theMultiplicities.begin(); while(it != theMultiplicities.end()){ p = it->first; file << it->second.first << " " << it->second.second << " " << p-pold << '\n'; it++; pold = p; } file << "sum of all probabilities: " << theMultiplicities.sum() << endl; file.close(); } void MPIHandler::statistics() const { ostream & file = generator()->misc(); string line = "=======================================" "=======================================\n"; for(unsigned int i=0; istatistics(file, tot); file << "\n"; } if(softInt_){ file << line << "Eikonalized and soft cross sections:\n\n" << "Model parameters: " << "ptmin: " << Ptmin_/GeV << " GeV" << ", mu2: " << invRadius_/sqr(1.*GeV) << " GeV2\n" << " " << "DL mode: " << DLmode_ << ", CMenergy: " << generator()->maximumCMEnergy()/GeV << " GeV" << '\n' << "hard inclusive cross section (mb): " << hardXSec_/millibarn << '\n' << "soft inclusive cross section (mb): " << softXSec_/millibarn << '\n' << "total cross section (mb): " << totalXSecExp()/millibarn << '\n' << "inelastic cross section (mb): " << inelXSec_/millibarn << '\n' + << "diffractive cross section (mb): " + << diffratio_*totalXSecExp()/millibarn << '\n' << "soft inv radius (GeV2): " << softMu2_/GeV2 << '\n' << "slope of soft pt spectrum (1/GeV2): " << beta_*sqr(1.*GeV) << '\n' << "Average hard multiplicity: " << avgNhard_ << '\n' << "Average soft multiplicity: " << avgNsoft_ << '\n' << line << endl; }else{ file << line << "Eikonalized and soft cross sections:\n\n" << "Model parameters: " << "ptmin: " << Ptmin_/GeV << " GeV" << ", mu2: " << invRadius_/sqr(1.*GeV) << " GeV2\n" << " " << ", CMenergy: " << generator()->maximumCMEnergy()/GeV << " GeV" << '\n' << "hard inclusive cross section (mb): " << hardXSec_/millibarn << '\n' << "Average hard multiplicity: " << avgNhard_ << '\n' << line << endl; } } unsigned int MPIHandler::multiplicity(unsigned int sel){ if(sel==0){//draw from the pretabulated distribution MPair m = theMultiplicities.select(UseRandom::rnd()); softMult_ = m.second; return m.first; } else{ //fixed multiplicities for the additional hard scatters if(additionalMultiplicities_.size() < sel) throw Exception() << "MPIHandler::multiplicity: process index " << "is out of range" << Exception::runerror; return additionalMultiplicities_[sel-1]; } } void MPIHandler::Probs(XSVector UEXSecs) { GSLIntegrator integrator; unsigned int iH(1), iS(0); double P(0.0); double P0(0.0);//the probability for i hard and zero soft scatters Length bmax(500.0*sqrt(millibarn)); //only one UE process will be drawn from a probability distribution, //so check that. assert(UEXSecs.size() == 1); for ( XSVector::const_iterator it = UEXSecs.begin(); it != UEXSecs.end(); ++it ) { iH = 0; //get the inel xsec Eikonalization inelint(this, -1, *it, softXSec_, softMu2_); - inelXSec_ = integrator.value(inelint, ZERO, bmax); + + + auto nondiffXS=integrator.value(inelint, ZERO, bmax); + inelXSec_ = nondiffXS+diffratio_*totalXSecExp(); avgNhard_ = 0.0; avgNsoft_ = 0.0; bmax = 10.0*sqrt(millibarn); do{//loop over hard ints if(Algorithm()>-1 && iH==0){ iH++; continue; } iS = 0; do{//loop over soft ints if( ( Algorithm() == -1 && iS==0 && iH==0 ) ){ iS++; continue; } Eikonalization integrand(this, iH*100+iS, *it, softXSec_, softMu2_); if(Algorithm() > 0){ P = integrator.value(integrand, ZERO, bmax) / *it; }else{ - P = integrator.value(integrand, ZERO, bmax) / inelXSec_; + P = integrator.value(integrand, ZERO, bmax) / nondiffXS; } //store the probability if(Algorithm()>-1){ theMultiplicities.insert(P, make_pair(iH-1, iS)); avgNhard_ += P*(iH-1); }else{ theMultiplicities.insert(P, make_pair(iH, iS)); avgNhard_ += P*(iH); } avgNsoft_ += P*iS; if(iS==0) P0 = P; iS++; } while ( (iS < maxScatters_) && (iS < 5 || P > 1.e-15 ) && softInt_ ); iH++; } while ( (iH < maxScatters_) && (iH < 5 || P0 > 1.e-15) ); } } // calculate the integrand Length Eikonalization::operator() (Length b) const { unsigned int Nhard(0), Nsoft(0); //fac is just: d^2b=fac*db despite that large number Length fac(Constants::twopi*b); double chiTot(( theHandler->OverlapFunction(b)*hardXSec_ + theHandler->OverlapFunction(b, softMu2_)*softXSec_ ) / 2.0); //total cross section wanted if(theoption == -2) return 2 * fac * ( 1 - exp(-chiTot) ); //inelastic cross section if(theoption == -1) return fac * ( 1 - exp(- 2.0 * chiTot) ); if(theoption >= 0){ //encode multiplicities as: N_hard*100 + N_soft Nhard = theoption/100; Nsoft = theoption%100; if(theHandler->Algorithm() > 0){ //P_n*sigma_hard: n-1 extra scatters + 1 hard scatterer != hardXSec_ return fac * Nhard * theHandler->poisson(b, hardXSec_, Nhard) * theHandler->poisson(b, softXSec_, Nsoft, softMu2_); }else{ //P_n*sigma_inel: n extra scatters return fac * theHandler->poisson(b, hardXSec_, Nhard) * theHandler->poisson(b, softXSec_, Nsoft, softMu2_); } }else{ throw Exception() << "Parameter theoption in Struct Eikonalization in " << "MPIHandler.cc has not allowed value" << Exception::runerror; return 0.0*meter; } } InvEnergy2 slopeBisection::operator() (Energy2 softMu2) const { GSLBisection root; //single component model TotalXSecBisection fn(handler_, softMu2); try{ softXSec_ = root.value(fn, 0*millibarn, 5000*millibarn); }catch(GSLBisection::IntervalError){ throw Exception() << "MPIHandler 2-Component model didn't work out." << Exception::runerror; } return handler_->slopeDiff(softXSec_, softMu2); } LengthDiff slopeInt::operator() (Length b) const { //fac is just: d^2b=fac*db Length fac(Constants::twopi*b); double chiTot(( handler_->OverlapFunction(b)*hardXSec_ + handler_->OverlapFunction(b, softMu2_)*softXSec_ ) / 2.0); InvEnergy2 b2 = sqr(b/hbarc); //B*sigma_tot return fac * b2 * ( 1 - exp(-chiTot) ); } double MPIHandler::factorial (unsigned int n) const { static double f[] = {1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880.,3.6288e6, 3.99168e7,4.790016e8,6.2270208e9,8.71782912e10,1.307674368e12, 2.0922789888e13,3.55687428096e14,6.402373705728e15,1.21645100408832e17, 2.43290200817664e18,5.10909421717094e19,1.12400072777761e21, 2.5852016738885e22,6.20448401733239e23,1.5511210043331e25, 4.03291461126606e26,1.08888694504184e28,3.04888344611714e29, 8.8417619937397e30,2.65252859812191e32,8.22283865417792e33, 2.63130836933694e35,8.68331761881189e36,2.95232799039604e38, 1.03331479663861e40,3.71993326789901e41,1.37637530912263e43, 5.23022617466601e44,2.03978820811974e46,8.15915283247898e47, 3.34525266131638e49,1.40500611775288e51,6.04152630633738e52, 2.65827157478845e54,1.1962222086548e56,5.50262215981209e57, 2.58623241511168e59,1.24139155925361e61,6.08281864034268e62, 3.04140932017134e64,1.55111875328738e66,8.06581751709439e67, 4.27488328406003e69,2.30843697339241e71,1.26964033536583e73, 7.10998587804863e74,4.05269195048772e76,2.35056133128288e78, 1.3868311854569e80,8.32098711274139e81,5.07580213877225e83, 3.14699732603879e85,1.98260831540444e87,1.26886932185884e89, 8.24765059208247e90,5.44344939077443e92,3.64711109181887e94, 2.48003554243683e96,1.71122452428141e98,1.19785716699699e100, 8.50478588567862e101,6.12344583768861e103,4.47011546151268e105, 3.30788544151939e107,2.48091408113954e109,1.88549470166605e111, 1.45183092028286e113,1.13242811782063e115,8.94618213078298e116, 7.15694570462638e118,5.79712602074737e120,4.75364333701284e122, 3.94552396972066e124,3.31424013456535e126,2.81710411438055e128, 2.42270953836727e130,2.10775729837953e132,1.85482642257398e134, 1.65079551609085e136,1.48571596448176e138,1.3520015276784e140, 1.24384140546413e142,1.15677250708164e144,1.08736615665674e146, 1.03299784882391e148,9.9167793487095e149,9.61927596824821e151, 9.42689044888325e153,9.33262154439442e155,9.33262154439442e157}; if(n > maxScatters_) throw Exception() << "MPIHandler::factorial called with too large argument" << Exception::runerror; else return f[n]; } InvArea MPIHandler::OverlapFunction(Length b, Energy2 mu2) const { if(mu2 == ZERO) mu2 = invRadius_; InvLength mu = sqrt(mu2)/hbarc; return (sqr(mu)/96/Constants::pi)*pow(mu*b, 3)*(gsl_sf_bessel_Kn(3, mu*b)); } double MPIHandler::poisson(Length b, CrossSection sigma, unsigned int N, Energy2 mu2) const { if(sigma > 0*millibarn){ return pow(OverlapFunction(b, mu2)*sigma, (double)N)/factorial(N) *exp(-OverlapFunction(b, mu2)*sigma); }else{ return (N==0) ? 1.0 : 0.0; } } CrossSection MPIHandler::totalXSecDiff(CrossSection softXSec, Energy2 softMu2) const { GSLIntegrator integrator; Eikonalization integrand(this, -2, hardXSec_, softXSec, softMu2); Length bmax = 500.0*sqrt(millibarn); - CrossSection tot = integrator.value(integrand, ZERO, bmax); - return (tot-totalXSecExp()); + CrossSection minbiasXS = integrator.value(integrand, ZERO, bmax); + + auto totalXS=totalXSecExp(); + CrossSection tot=minbiasXS+diffratio_*totalXS; + + return (tot-totalXS); } InvEnergy2 MPIHandler::slopeDiff(CrossSection softXSec, Energy2 softMu2) const { GSLIntegrator integrator; Eikonalization integrand(this, -2, hardXSec_, softXSec, softMu2); Length bmax = 500.0*sqrt(millibarn); CrossSection tot = integrator.value(integrand, ZERO, bmax); slopeInt integrand2(this, hardXSec_, softXSec, softMu2); return integrator.value(integrand2, ZERO, bmax)/tot - slopeExp(); } CrossSection MPIHandler::totalXSecExp() const { if(totalXSecExp_ != 0*millibarn) return totalXSecExp_; double pom_old = 0.0808; CrossSection coef_old = 21.7*millibarn; double pom_new_hard = 0.452; CrossSection coef_new_hard = 0.0139*millibarn; double pom_new_soft = 0.0667; CrossSection coef_new_soft = 24.22*millibarn; Energy energy(generator()->maximumCMEnergy()); switch(DLmode_){ case 1://old DL extrapolation return coef_old * pow(energy/GeV, 2*pom_old); break; case 2://old DL extrapolation fixed to CDF return 81.8*millibarn * pow(energy/1800.0/GeV, 2*pom_old); break; case 3://new DL extrapolation - return 0.8*coef_new_hard * pow(energy/GeV, 2*pom_new_hard) + + return coef_new_hard * pow(energy/GeV, 2*pom_new_hard) + coef_new_soft * pow(energy/GeV, 2*pom_new_soft); break; default: throw Exception() << "MPIHandler::totalXSecExp non-existing mode selected" << Exception::runerror; } } InvEnergy2 MPIHandler::slopeExp() const{ //Currently return the slope as calculated in the pomeron fit by //Donnachie & Landshoff Energy energy(generator()->maximumCMEnergy()); //slope at Energy e_0 = 1800*GeV; InvEnergy2 b_0 = 17/GeV2; return b_0 + log(energy/e_0)/GeV2; } void MPIHandler::overrideUECuts() { if(energyExtrapolation_==1) Ptmin_ = EEparamA_ * log(generator()->maximumCMEnergy() / EEparamB_); else if(energyExtrapolation_==2) Ptmin_ = pT0_*pow(double(generator()->maximumCMEnergy()/refScale_),b_); else assert(false); // create a new SimpleKTCut object with the calculated ptmin value Ptr::pointer newUEktCut = new_ptr(SimpleKTCut(Ptmin_)); newUEktCut->init(); newUEktCut->initrun(); // create a new Cuts object with MHatMin = 2 * Ptmin_ CutsPtr newUEcuts = new_ptr(Cuts(2*Ptmin_)); newUEcuts->add(dynamic_ptr_cast(newUEktCut)); newUEcuts->init(); newUEcuts->initrun(); // replace the old Cuts object cuts()[0] = newUEcuts; } void MPIHandler::persistentOutput(PersistentOStream & os) const { os << theMultiplicities << theHandler << theSubProcesses << theCuts << theProcessHandlers << additionalMultiplicities_ << identicalToUE_ << ounit(PtOfQCDProc_, GeV) << ounit(Ptmin_, GeV) << ounit(hardXSec_, millibarn) << ounit(softXSec_, millibarn) << ounit(beta_, 1/GeV2) << algorithm_ << ounit(invRadius_, GeV2) << numSubProcs_ << colourDisrupt_ << softInt_ << twoComp_ << DLmode_ << ounit(totalXSecExp_, millibarn) << energyExtrapolation_ << ounit(EEparamA_, GeV) << ounit(EEparamB_, GeV) << ounit(refScale_,GeV) << ounit(pT0_,GeV) << b_ << avgNhard_ << avgNsoft_ << softMult_ << ounit(inelXSec_, millibarn) - << ounit(softMu2_, GeV2); + << ounit(softMu2_, GeV2) << diffratio_; } void MPIHandler::persistentInput(PersistentIStream & is, int) { is >> theMultiplicities >> theHandler >> theSubProcesses >> theCuts >> theProcessHandlers >> additionalMultiplicities_ >> identicalToUE_ >> iunit(PtOfQCDProc_, GeV) >> iunit(Ptmin_, GeV) >> iunit(hardXSec_, millibarn) >> iunit(softXSec_, millibarn) >> iunit(beta_, 1/GeV2) >> algorithm_ >> iunit(invRadius_, GeV2) >> numSubProcs_ >> colourDisrupt_ >> softInt_ >> twoComp_ >> DLmode_ >> iunit(totalXSecExp_, millibarn) >> energyExtrapolation_ >> iunit(EEparamA_, GeV) >> iunit(EEparamB_, GeV) >> iunit(refScale_,GeV) >> iunit(pT0_,GeV) >> b_ >> avgNhard_ >> avgNsoft_ >> softMult_ >> iunit(inelXSec_, millibarn) - >> iunit(softMu2_, GeV2); + >> iunit(softMu2_, GeV2) >> diffratio_; currentHandler_ = this; } void MPIHandler::clean() { // ThePEG's event handler's usual event cleanup doesn't reach these // XCombs. Need to do it by hand here. for ( size_t i = 0; i < theSubProcesses.size(); ++i ) { theSubProcesses[i]->pExtractor()->lastXCombPtr()->clean(); } } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigMPIHandler("Herwig::MPIHandler", "JetCuts.so SimpleKTCut.so HwMPI.so"); void MPIHandler::Init() { static ClassDocumentation documentation ("The MPIHandler class is the main administrator of the multiple interaction model", "The underlying event was simulated with an eikonal model for multiple partonic interactions." "Details can be found in Ref.~\\cite{Bahr:2008dy,Bahr:2009ek}.", "%\\cite{Bahr:2008dy}\n" "\\bibitem{Bahr:2008dy}\n" " M.~Bahr, S.~Gieseke and M.~H.~Seymour,\n" " ``Simulation of multiple partonic interactions in Herwig,''\n" " JHEP {\\bf 0807}, 076 (2008)\n" " [arXiv:0803.3633 [hep-ph]].\n" " %%CITATION = JHEPA,0807,076;%%\n" "\\bibitem{Bahr:2009ek}\n" " M.~Bahr, J.~M.~Butterworth, S.~Gieseke and M.~H.~Seymour,\n" " ``Soft interactions in Herwig,''\n" " arXiv:0905.4671 [hep-ph].\n" " %%CITATION = ARXIV:0905.4671;%%\n" ); static RefVector interfaceSubhandlers ("SubProcessHandlers", "The list of sub-process handlers used in this EventHandler. ", &MPIHandler::theSubProcesses, -1, false, false, true, false, false); static RefVector interfaceCuts ("Cuts", "List of cuts used for the corresponding list of subprocesses. These cuts " "should not be overidden in individual sub-process handlers.", &MPIHandler::theCuts, -1, false, false, true, false, false); static Parameter interfaceInvRadius ("InvRadius", "The inverse hadron radius squared used in the overlap function", &MPIHandler::invRadius_, GeV2, 2.0*GeV2, 0.2*GeV2, 4.0*GeV2, true, false, Interface::limited); static ParVector interfaceadditionalMultiplicities ("additionalMultiplicities", "specify the multiplicities of secondary hard processes (multiple parton scattering)", &MPIHandler::additionalMultiplicities_, -1, 0, 0, 3, false, false, true); static Parameter interfaceIdenticalToUE ("IdenticalToUE", "Specify which of the hard processes is identical to the UE one (QCD dijets)", &MPIHandler::identicalToUE_, -1, 0, 0, false, false, Interface::nolimits); static Parameter interfacePtOfQCDProc ("PtOfQCDProc", "Specify the value of the pt cutoff for the process that is identical to the UE one", &MPIHandler::PtOfQCDProc_, GeV, -1.0*GeV, ZERO, ZERO, false, false, Interface::nolimits); static Parameter interfacecolourDisrupt ("colourDisrupt", "Fraction of connections to additional subprocesses, which are colour disrupted.", &MPIHandler::colourDisrupt_, 0.0, 0.0, 1.0, false, false, Interface::limited); + + static Parameter interfaceDiffRatio + ("DiffractiveRatio", + "Fraction of diffractive cross section in total cross section.", + &MPIHandler::diffratio_, + 0.2, 0.0, 1.0, + false, false, Interface::limited); + static Switch interfacesoftInt ("softInt", "Switch to enable soft interactions", &MPIHandler::softInt_, true, false, false); static SwitchOption interfacesoftIntYes (interfacesoftInt, "Yes", "enable the two component model", true); static SwitchOption interfacesoftIntNo (interfacesoftInt, "No", "disable the model", false); static Switch interEnergyExtrapolation ("EnergyExtrapolation", "Switch to ignore the cuts object at MPIHandler:Cuts[0]. " "Instead, extrapolate the pt cut.", &MPIHandler::energyExtrapolation_, 2, false, false); static SwitchOption interEnergyExtrapolationLog (interEnergyExtrapolation, "Log", "Use logarithmic dependence, ptmin = A * log (sqrt(s) / B).", 1); static SwitchOption interEnergyExtrapolationPower (interEnergyExtrapolation, "Power", "Use power law, ptmin = pt_0 * (sqrt(s) / E_0)^b.", 2); static SwitchOption interEnergyExtrapolationNo (interEnergyExtrapolation, "No", "Use manually set value for the minimal pt, " "specified in MPIHandler:Cuts[0]:OneCuts[0]:MinKT.", 0); static Parameter interfaceEEparamA ("EEparamA", "Parameter A in the empirical parametrization " "ptmin = A * log (sqrt(s) / B)", &MPIHandler::EEparamA_, GeV, 0.6*GeV, ZERO, Constants::MaxEnergy, false, false, Interface::limited); static Parameter interfaceEEparamB ("EEparamB", "Parameter B in the empirical parametrization " "ptmin = A * log (sqrt(s) / B)", &MPIHandler::EEparamB_, GeV, 39.0*GeV, ZERO, Constants::MaxEnergy, false, false, Interface::limited); static Switch interfacetwoComp ("twoComp", "switch to enable the model with a different radius for soft interactions", &MPIHandler::twoComp_, true, false, false); static SwitchOption interfacetwoCompYes (interfacetwoComp, "Yes", "enable the two component model", true); static SwitchOption interfacetwoCompNo (interfacetwoComp, "No", "disable the model", false); static Parameter interfaceMeasuredTotalXSec ("MeasuredTotalXSec", "Value for the total cross section (assuming rho=0). If non-zero, this " "overwrites the Donnachie-Landshoff parametrizations.", &MPIHandler::totalXSecExp_, millibarn, 0.0*millibarn, 0.0*millibarn, 0*millibarn, false, false, Interface::lowerlim); static Switch interfaceDLmode ("DLmode", "Choice of Donnachie-Landshoff parametrization for the total cross section.", &MPIHandler::DLmode_, 2, false, false); static SwitchOption interfaceDLmodeStandard (interfaceDLmode, "Standard", "Standard parametrization with s**0.08", 1); static SwitchOption interfaceDLmodeCDF (interfaceDLmode, "CDF", "Standard parametrization but normalization fixed to CDF's measured value", 2); static SwitchOption interfaceDLmodeNew (interfaceDLmode, "New", "Parametrization taking hard and soft pomeron contributions into account", 3); static Parameter interfaceReferenceScale ("ReferenceScale", "The reference energy for power law energy extrapolation of pTmin", &MPIHandler::refScale_, GeV, 7000.0*GeV, 0.0*GeV, 20000.*GeV, false, false, Interface::limited); static Parameter interfacepTmin0 ("pTmin0", "The pTmin at the reference scale for power law extrapolation of pTmin.", &MPIHandler::pT0_, GeV, 3.11*GeV, 0.0*GeV, 10.0*GeV, false, false, Interface::limited); static Parameter interfacePower ("Power", "The power for power law extrapolation of the pTmin cut-off.", &MPIHandler::b_, 0.21, 0.0, 10.0, false, false, Interface::limited); } diff --git a/UnderlyingEvent/MPIHandler.h b/UnderlyingEvent/MPIHandler.h --- a/UnderlyingEvent/MPIHandler.h +++ b/UnderlyingEvent/MPIHandler.h @@ -1,878 +1,895 @@ // -*- C++ -*- // // MPIHandler.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef HERWIG_MPIHandler_H #define HERWIG_MPIHandler_H // // This is the declaration of the MPIHandler class. // #include "ThePEG/Interface/Interfaced.h" #include "ThePEG/Handlers/StandardEventHandler.h" #include "ThePEG/Repository/EventGenerator.h" #include "Herwig/PDT/StandardMatchers.h" #include "Herwig/Utilities/GSLBisection.h" //#include "Herwig/Utilities/GSLMultiRoot.h" #include "Herwig/Utilities/GSLIntegrator.h" #include "Herwig/Shower/UEBase.h" #include #include "ProcessHandler.h" #include "MPIHandler.fh" namespace Herwig { using namespace ThePEG; /** \ingroup UnderlyingEvent * \class MPIHandler * This class is responsible for generating additional * semi hard partonic interactions. * * \author Manuel B\"ahr * * @see \ref MPIHandlerInterfaces "The interfaces" * defined for MPIHandler. * @see ProcessHandler * @see ShowerHandler * @see HwRemDecayer */ class MPIHandler: public UEBase { /** * Maximum number of scatters */ static const unsigned int maxScatters_ = 99; /** * Class for the integration is a friend to access private members */ friend struct Eikonalization; friend struct TotalXSecBisection; friend struct slopeAndTotalXSec; friend struct slopeInt; friend struct slopeBisection; public: /** A vector of SubProcessHandlers. */ typedef vector SubHandlerList; /** A vector of Cuts. */ typedef vector CutsList; /** A vector of ProcessHandlers. */ typedef vector ProcessHandlerList; /** A vector of cross sections. */ typedef vector XSVector; /** A pair of multiplicities: hard, soft. */ typedef pair MPair; /** @name Standard constructors and destructors. */ //@{ /** * The default constructor. */ MPIHandler(): softMult_(0), identicalToUE_(-1), PtOfQCDProc_(-1.0*GeV), Ptmin_(-1.0*GeV), hardXSec_(0*millibarn), softXSec_(0*millibarn), totalXSecExp_(0*millibarn), softMu2_(ZERO), beta_(100.0/GeV2), algorithm_(2), numSubProcs_(0), colourDisrupt_(0.0), softInt_(true), twoComp_(true), DLmode_(2), avgNhard_(0.0), avgNsoft_(0.0), energyExtrapolation_(2), EEparamA_(0.6*GeV), EEparamB_(37.5*GeV), refScale_(7000.*GeV), pT0_(3.11*GeV), b_(0.21) {} /** * The destructor. */ virtual ~MPIHandler(){} //@} public: /** @name Methods for the MPI generation. */ //@{ /* * @return true if for this beam setup MPI can be generated */ virtual bool beamOK() const; /** * Return true or false depending on whether soft interactions are enabled. */ virtual bool softInt() const {return softInt_;} /** * Get the soft multiplicity from the pretabulated multiplicity * distribution. Generated in multiplicity in the first place. * @return the number of extra soft events in this collision */ virtual unsigned int softMultiplicity() const {return softMult_;} /** * Sample from the pretabulated multiplicity distribution. * @return the number of extra events in this collision */ virtual unsigned int multiplicity(unsigned int sel=0); /** * Select a StandardXComb according to it's weight * @return that StandardXComb Object * @param sel is the subprocess that should be returned, * if more than one is specified. */ virtual tStdXCombPtr generate(unsigned int sel=0); //@} /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); /** * Initialize this Multiple Interaction handler and all related objects needed to * generate additional events. */ virtual void initialize(); /** * Finalize this Multiple Interaction handler and all related objects needed to * generate additional events. */ virtual void finalize(); /** * Clean up the XCombs from our subprocesses after each event. * ThePEG cannot see them, so the usual cleaning misses these. */ virtual void clean(); /** * Write out accumulated statistics about integrated cross sections. */ void statistics() const; /** * The level of statistics. Controlls the amount of statistics * written out after each run to the EventGenerators * .out file. Simply the EventHandler method is called here. */ int statLevel() const {return eventHandler()->statLevel();} /** * Return the hard cross section above ptmin */ CrossSection hardXSec() const { return hardXSec_; } /** * Return the soft cross section below ptmin */ CrossSection softXSec() const { return softXSec_; } /** * Return the inelastic cross section */ CrossSection inelasticXSec() const { return inelXSec_; } + + /** + * Return the diffractive cross section assumed by the model. + * The diffractive cross section is seen as part of the + * inelastic cross section. + */ + CrossSection diffractiveXSec() const { + static auto totalXS=totalXSecExp(); + return diffratio_*totalXS; + } + + /** @name Simple access functions. */ //@{ /** * Return the ThePEG::EventHandler assigned to this handler. * This methods shadows ThePEG::StepHandler::eventHandler(), because * it is not virtual in ThePEG::StepHandler. This is ok, because this * method would give a null-pointer at some stages, whereas this method * gives access to the explicitely copied pointer (in initialize()) * to the ThePEG::EventHandler. */ tEHPtr eventHandler() const {return theHandler;} /** * Return the current handler */ static const MPIHandler * currentHandler() { return currentHandler_; } /** * Return theAlgorithm. */ virtual int Algorithm() const {return algorithm_;} /** * Return the ptmin parameter of the model */ virtual Energy Ptmin() const { if(Ptmin_ > ZERO) return Ptmin_; else throw Exception() << "MPIHandler::Ptmin called without initialize before" << Exception::runerror; } /** * Return the slope of the soft pt spectrum as calculated. */ virtual InvEnergy2 beta() const { if(beta_ != 100.0/GeV2) return beta_; else throw Exception() << "MPIHandler::beta called without initialization" << Exception::runerror; } /** * Return the pt Cutoff of the Interaction that is identical to the UE * one. */ virtual Energy PtForVeto() const {return PtOfQCDProc_;} /** * Return the number of additional "hard" processes ( = multiple * parton scattering) */ virtual unsigned int additionalHardProcs() const {return numSubProcs_-1;} /** * Return the fraction of colour disrupted connections to the * suprocesses. */ virtual double colourDisrupt() const {return colourDisrupt_;} protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} private: /** * Access the list of sub-process handlers. */ const SubHandlerList & subProcesses() const {return theSubProcesses;} /** * Access the list of sub-process handlers. */ SubHandlerList & subProcesses() {return theSubProcesses;} /** * Access the list of cuts. */ const CutsList & cuts() const {return theCuts;} /** * Access the list of cuts. */ CutsList & cuts() {return theCuts;} /** * Access the list of sub-process handlers. */ const ProcessHandlerList & processHandlers() const {return theProcessHandlers;} /** * Access the list of sub-process handlers. */ ProcessHandlerList & processHandlers() {return theProcessHandlers;} /** * Method to calculate the individual probabilities for N scatters in the event. * @param UEXSecs is(are) the inclusiv cross section(s) for the UE process(es). */ void Probs(XSVector UEXSecs); /** * Debug method to check the individual probabilities. * @param filename is the file the output gets written to */ void MultDistribution(string filename) const; /** * Return the value of the Overlap function A(b) for a given impact * parameter \a b. * @param b impact parameter * @param mu2 = inv hadron radius squared. 0 will use the value of * invRadius_ * @return inverse area. */ InvArea OverlapFunction(Length b, Energy2 mu2=ZERO) const; /** * Method to calculate the poisson probability for expectation value * \f$ = A(b)\sigma\f$, and multiplicity N. */ double poisson(Length b, CrossSection sigma, unsigned int N, Energy2 mu2=ZERO) const; /** * Return n! */ double factorial (unsigned int n) const; /** * Returns the total cross section for the current CMenergy. The * decision which parametrization will be used is steered by a * external parameter of this class. */ CrossSection totalXSecExp() const; /** * Difference of the calculated total cross section and the * experimental one from totalXSecExp. * @param softXSec = the soft cross section that is used * @param softMu2 = the soft radius, if 0 the hard radius will be used */ CrossSection totalXSecDiff(CrossSection softXSec, Energy2 softMu2=ZERO) const; /** * Difference of the calculated elastic slope and the * experimental one from slopeExp. * @param softXSec = the soft cross section that is used * @param softMu2 = the soft radius, if 0 the hard radius will be used */ InvEnergy2 slopeDiff(CrossSection softXSec, Energy2 softMu2=ZERO) const; /** * Returns the value of the elastic slope for the current CMenergy. * The decision which parametrization will be used is steered by a * external parameter of this class. */ InvEnergy2 slopeExp() const; /** * Calculate the minimal transverse momentum from the extrapolation */ void overrideUECuts(); private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ MPIHandler & operator=(const MPIHandler &) = delete; /** * A pointer to the EventHandler that calls us. Has to be saved, because the * method eventHandler() inherited from ThePEG::StepHandler returns a null-pointer * sometimes. Leif changed that in r1053 so that a valid pointer is present, when * calling doinitrun(). */ tEHPtr theHandler; /** * The list of SubProcessHandlers. */ SubHandlerList theSubProcesses; /** * The kinematical cuts used for this collision handler. */ CutsList theCuts; /** * List of ProcessHandler used to sample different processes independently */ ProcessHandlerList theProcessHandlers; /** * A ThePEG::Selector where the individual Probabilities P_N are stored * and the actual Multiplicities can be selected. */ Selector theMultiplicities; /** * Variable to store the soft multiplicity generated for a event. This * has to be stored as it is generated at the time of the hard * additional interactions but used later on. */ unsigned int softMult_; /** * Variable to store the multiplicity of the second hard process */ vector additionalMultiplicities_; /** * Variable to store the information, which process is identical to * the UE one (QCD dijets). * 0 means "real" hard one * n>0 means the nth additional hard scatter * -1 means no one! */ int identicalToUE_; /** * Variable to store the minimal pt of the process that is identical * to the UE one. This only has to be set, if it can't be determined * automatically (i.e. when reading QCD LesHouches files in). */ Energy PtOfQCDProc_; /** * Variable to store the parameter ptmin */ Energy Ptmin_; /** * Variable to store the hard cross section above ptmin */ CrossSection hardXSec_; /** * Variable to store the final soft cross section below ptmin */ CrossSection softXSec_; /** * Variable to store the inelastic cross section */ CrossSection inelXSec_; /** * Variable to store the total pp cross section (assuming rho=0!) as * measured at LHC. If this variable is set, this value is used in the * subsequent run instead of any of the Donnachie-Landshoff * parametrizations. */ CrossSection totalXSecExp_; /** * Variable to store the soft radius, that is calculated during * initialization for the two-component model. */ Energy2 softMu2_; /** * slope to the non-perturbative pt spectrum: \f$d\sigma/dp_T^2 = A \exp * (- beta p_T^2)\f$. Its value is determined durint initialization. */ InvEnergy2 beta_; /** * Switch to be set from outside to determine the algorithm used for * UE activity. */ int algorithm_; /** * Inverse hadron Radius squared \f$ (\mu^2) \f$. Used inside the overlap function. */ Energy2 invRadius_; /** * Member variable to store the actual number of separate SubProcesses */ unsigned int numSubProcs_; /** * Variable to store the relative number of colour disrupted * connections to additional subprocesses. This variable is used in * Herwig::HwRemDecayer but store here, to have access to all * parameters through one Object. */ double colourDisrupt_; /** * Flag to store whether soft interactions, i.e. pt < ptmin should be * simulated. */ bool softInt_; /** * Flag to steer wheather the soft part has a different radius, that * will be dynamically fixed. */ bool twoComp_; /** * Switch to determine which Donnachie & Landshoff parametrization * should be used. */ unsigned int DLmode_; /** * Variable to store the average hard multiplicity. */ double avgNhard_; /** * Variable to store the average soft multiplicity. */ double avgNsoft_; /** * The current handler */ static MPIHandler * currentHandler_; /** * Flag to store whether to calculate the minimal UE pt according to an * extrapolation formula or whether to use MPIHandler:Cuts[0]:OneCuts[0]:MinKT */ unsigned int energyExtrapolation_; /** * Parameters for the energy extrapolation formula */ Energy EEparamA_; Energy EEparamB_; Energy refScale_; Energy pT0_; double b_; + /** + * Parameters to set the fraction of diffractive cross section in the total cross section. + */ + double diffratio_=0.2; + protected: /** @cond EXCEPTIONCLASSES */ /** * Exception class used by the MultipleInteractionHandler, when something * during initialization went wrong. * \todo understand!!! */ class InitError: public Exception {}; /** @endcond */ }; } namespace Herwig { /** * A struct for the 2D root finding that is necessary to determine the * soft cross section and the soft radius that is needed to describe * the total cross section correctly. * NOT IN USE CURRENTLY */ struct slopeAndTotalXSec : public GSLHelper { public: /** * Constructor */ slopeAndTotalXSec(tcMPIHPtr handler): handler_(handler) {} /** second argument type */ typedef Energy2 ArgType2; /** second value type */ typedef InvEnergy2 ValType2; /** first element of the vector like function to find root for * @param softXSec soft cross-section * @param softMu2 \f$\mu^2\f$ */ CrossSection f1(ArgType softXSec, ArgType2 softMu2) const { return handler_->totalXSecDiff(softXSec, softMu2); } /** second element of the vector like function to find root for * @param softXSec soft cross-section * @param softMu2 \f$\mu^2\f$ */ InvEnergy2 f2(ArgType softXSec, ArgType2 softMu2) const { return handler_->slopeDiff(softXSec, softMu2); } /** provide the actual units of use */ virtual ValType vUnit() const {return 1.0*millibarn;} /** otherwise rounding errors may get significant */ virtual ArgType aUnit() const {return 1.0*millibarn;} /** provide the actual units of use */ ValType2 vUnit2() const {return 1.0/GeV2;} /** otherwise rounding errors may get significant */ ArgType2 aUnit2() const {return GeV2;} private: /** * Pointer to the handler */ tcMPIHPtr handler_; }; /** * A struct for the root finding that is necessary to determine the * slope of the soft pt spectrum to match the soft cross section */ struct betaBisection : public GSLHelper{ public: /** * Constructor. * @param soft = soft cross section, i.e. the integral of the soft * pt spectrum f(u=p_T^2) = dsig exp(-beta*u/u_min) * @param dsig = dsigma_hard/dp_T^2 at the p_T cutoff * @param ptmin = p_T cutoff */ betaBisection(CrossSection soft, DiffXSec dsig, Energy ptmin) : softXSec_(soft), dsig_(dsig), ptmin_(ptmin) {} /** * Operator that is used inside the GSLBisection class */ virtual Energy2 operator ()(InvEnergy2 beta) const { if( fabs(beta*GeV2) < 1.E-4 ) beta = (beta > ZERO) ? 1.E-4/GeV2 : -1.E-4/GeV2; return (exp(beta*sqr(ptmin_)) - 1.0)/beta - softXSec_/dsig_; } /** provide the actual units of use */ virtual ValType vUnit() const {return 1.0*GeV2;} /** provide the actual units of use */ virtual ArgType aUnit() const {return 1.0/GeV2;} private: /** soft cross section */ CrossSection softXSec_; /** dsigma/dp_T^2 at ptmin */ DiffXSec dsig_; /** pt cutoff */ Energy ptmin_; }; /** * A struct for the root finding that is necessary to determine the * soft cross section and soft mu2 that are needed to describe the * total cross section AND elastic slope correctly. */ struct slopeBisection : public GSLHelper { public: /** Constructor */ slopeBisection(tcMPIHPtr handler) : handler_(handler) {} /** * Return the difference of the calculated elastic slope to the * experimental one for a given value of the soft mu2. During that, * the soft cross section get fixed. */ InvEnergy2 operator ()(Energy2 arg) const; /** Return the soft cross section that has been calculated */ CrossSection softXSec() const {return softXSec_;} private: /** const pointer to the MPIHandler to give access to member functions.*/ tcMPIHPtr handler_; /** soft cross section that is determined on the fly.*/ mutable CrossSection softXSec_; }; /** * A struct for the root finding that is necessary to determine the * soft cross section that is needed to describe the total cross * section correctly. */ struct TotalXSecBisection : public GSLHelper { public: /** * Constructor * @param handler The handler * @param softMu2 \f$\mu^2\f$ */ TotalXSecBisection(tcMPIHPtr handler, Energy2 softMu2=ZERO): handler_(handler), softMu2_(softMu2) {} /** * operator to return the cross section * @param argument input cross section */ CrossSection operator ()(CrossSection argument) const { return handler_->totalXSecDiff(argument, softMu2_); } /** provide the actual units of use */ virtual ValType vUnit() const {return 1.0*millibarn;} /** otherwise rounding errors may get significant */ virtual ArgType aUnit() const {return 1.0*millibarn;} private: /** * The handler */ tcMPIHPtr handler_; /** * \f$\mu^2\f$ */ Energy2 softMu2_; }; /** * Typedef for derivative of the length */ typedef decltype(mm/GeV2) LengthDiff; /** * A struct for the integrand for the slope */ struct slopeInt : public GSLHelper{ public: /** Constructor * @param handler The handler * @param hard The hard cross section * @param soft The soft cross section * @param softMu2 \f$\mu^2\f$ */ slopeInt(tcMPIHPtr handler, CrossSection hard, CrossSection soft=0*millibarn, Energy2 softMu2=ZERO) : handler_(handler), hardXSec_(hard), softXSec_(soft), softMu2_(softMu2) {} /** * Operator to return the answer * @param arg The argument */ ValType operator ()(ArgType arg) const; private: /** * Pointer to the Handler that calls this integrand */ tcMPIHPtr handler_; /** * The hard cross section to be eikonalized */ CrossSection hardXSec_; /** * The soft cross section to be eikonalized. Default is zero */ CrossSection softXSec_; /** * The inv radius^2 of the soft interactions. */ Energy2 softMu2_; }; /** * A struct for the eikonalization of the inclusive cross section. */ struct Eikonalization : public GSLHelper{ /** * The constructor * @param handler is the pointer to the MPIHandler to get access to * MPIHandler::OverlapFunction and member variables of the MPIHandler. * @param option is a flag, whether the inelastic or the total * @param handler The handler * @param hard The hard cross section * @param soft The soft cross section * @param softMu2 \f$\mu^2\f$ * cross section should be returned (-2 or -1). For option = N > 0 the integrand * is N*(A(b)*sigma)^N/N! exp(-A(b)*sigma) this is the P_N*sigma where * P_N is the Probability of having exactly N interaction (including the hard one) * This is equation 14 from "Jimmy4: Multiparton Interactions in HERWIG for the LHC" */ Eikonalization(tcMPIHPtr handler, int option, CrossSection hard, CrossSection soft=0*millibarn, Energy2 softMu2=ZERO) : theHandler(handler), theoption(option), hardXSec_(hard), softXSec_(soft), softMu2_(softMu2) {} /** * Get the function value */ Length operator ()(Length argument) const; private: /** * Pointer to the Handler that calls this integrand */ tcMPIHPtr theHandler; /** * A flag to switch between the calculation of total and inelastic cross section * or calculations for the individual probabilities. See the constructor */ int theoption; /** * The hard cross section to be eikonalized */ CrossSection hardXSec_; /** * The soft cross section to be eikonalized. Default is zero */ CrossSection softXSec_; /** * The inv radius^2 of the soft interactions. */ Energy2 softMu2_; }; } #endif /* HERWIG_MPIHandler_H */ diff --git a/src/LHC-MB.in b/src/LHC-MB.in --- a/src/LHC-MB.in +++ b/src/LHC-MB.in @@ -1,102 +1,103 @@ # -*- ThePEG-repository -*- ################################################################################ # This file contains our best tune to UE data from ATLAS at 7 TeV. More recent # tunes and tunes for other centre-of-mass energies as well as more usage # instructions can be obtained from this Herwig wiki page: # http://projects.hepforge.org/herwig/trac/wiki/MB_UE_tunes # The model for soft interactions and diffractions is explained in # [S. Gieseke, P. Kirchgaesser, F. Loshaj, arXiv:1612.04701] ################################################################################ read snippets/PPCollider.in ################################################## # Technical parameters for this run ################################################## cd /Herwig/Generators ################################################## # LHC physics parameters (override defaults here) ################################################## set EventGenerator:EventHandler:LuminosityFunction:Energy 7000.0 # Intrinsic pT tune extrapolated to LHC energy set /Herwig/Shower/ShowerHandler:IntrinsicPtGaussian 2.2*GeV # Minimum Bias read snippets/MB.in # Read in parameters of the soft model recommended for MB/UE simulations read snippets/SoftTune.in #Diffraction model read snippets/Diffraction.in #Turn on Baryonic reconnection read snippets/BaryonicReconnection.in # Normalization of the Min bias cross section for correct diffractive cross section set /Herwig/MatrixElements/MEMinBias:csNorm 0.01 # Use LHC parametrization of the cross section -set /Herwig/UnderlyingEvent/MPIHandler:DLmode 3 +set /Herwig/UnderlyingEvent/MPIHandler:DLmode 2 #some preliminary parameters for the MPI model which need to be tuned #TODO set /Herwig/UnderlyingEvent/MPIHandler:pTmin0 3.02 -set /Herwig/UnderlyingEvent/MPIHandler:InvRadius 0.9 +set /Herwig/UnderlyingEvent/MPIHandler:InvRadius 1.4 set /Herwig/UnderlyingEvent/MPIHandler:Power 0.308 set /Herwig/Partons/RemnantDecayer:ladderMult 0.45 set /Herwig/Partons/RemnantDecayer:ladderbFactor 1.0 +set /Herwig/UnderlyingEvent/MPIHandler:DiffractiveRatio 0.17 #ordered pTs -set /Herwig/Partons/RemnantDecayer:PtDistribution 5 +set /Herwig/Partons/RemnantDecayer:PtDistribution 4 #use random remnant connections set /Herwig/Partons/RemnantDecayer:RandomConnection No # set the correct PDFs set /Herwig/Partons/HardLOPDF:PDFName CT14lo set /Herwig/Partons/ShowerLOPDF:PDFName CT14lo set /Herwig/Partons/MPIPDF:PDFName CT14lo set /Herwig/Partons/RemnantPDF:PDFName CT14lo #set /Herwig/Shower/ShowerHandler:CascadeHandler NULL # Switch off parton shower #set /Herwig/Hadronization/HadronizationHandler:HadronizationHandler NULL # do Soft interactions #set /Herwig/Hadronization/ColourReconnector:ColourReconnection No ################################################# # Analyses ################################################## #Comment these lines out in order to use rivet analyses -#cd /Herwig/Analysis -#create ThePEG::RivetAnalysis RivetAnalysis RivetAnalysis.so +cd /Herwig/Analysis +create ThePEG::RivetAnalysis RivetAnalysis RivetAnalysis.so -#cd /Herwig/Generators -#insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/RivetAnalysis +cd /Herwig/Generators +insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/RivetAnalysis #Some example analyses -#insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2012_I1084540 -#insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2010_S8918562 +insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2012_I1084540 +insert /Herwig/Analysis/RivetAnalysis:Analyses 0 ATLAS_2010_S8918562 #set /Herwig/Analysis/Plot:EventNumber 54 #cd /Herwig/Generators #insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/Plot #insert EventGenerator:AnalysisHandlers 0 /Herwig/Analysis/HepMCFile #set /Herwig/Analysis/HepMCFile:PrintEvent 1000000 #set /Herwig/Analysis/HepMCFile:Format GenEvent #set /Herwig/Analysis/HepMCFile:Units GeV_mm #set /Herwig/Analysis/HepMCFile:Filename events.fifo ################################################## # Save run for later usage with 'Herwig run' ################################################## cd /Herwig/Generators saverun LHC-MB EventGenerator diff --git a/src/snippets/Diffraction.in b/src/snippets/Diffraction.in --- a/src/snippets/Diffraction.in +++ b/src/snippets/Diffraction.in @@ -1,104 +1,114 @@ ################################################## # Create separate SubProcessHandler for Diffraction ################################################## cd /Herwig/MatrixElements create Herwig::MEDiffraction MEDiffractionLeft set MEDiffractionLeft:DiffDirection Left create Herwig::MEDiffraction MEDiffractionRight set MEDiffractionRight:DiffDirection Right create Herwig::MEDiffraction MEDiffractionDouble set MEDiffractionDouble:DiffDirection Both create Herwig::MEDiffraction MEDiffractionDeltaLeft set MEDiffractionDeltaLeft:DiffDirection Left create Herwig::MEDiffraction MEDiffractionDeltaRight set MEDiffractionDeltaRight:DiffDirection Right create Herwig::MEDiffraction MEDiffractionDeltaDouble set MEDiffractionDeltaDouble:DiffDirection Both # Make a parton extractor for diffraction cd /Herwig/Partons cp PPExtractor DiffPPExtractor set DiffPPExtractor:FirstPDF /Herwig/Partons/NoPDF set DiffPPExtractor:SecondPDF /Herwig/Partons/NoPDF cd /Herwig/MatrixElements/ # Create Diffraction SubProcessHandler cp SubProcess QCDDiffraction # Assign the PartonExtractor to the SubProcessHandler set QCDDiffraction:PartonExtractor /Herwig/Partons/DiffPPExtractor # Use only Delta as final excited state (Yes/No) set MEDiffractionLeft:DeltaOnly No set MEDiffractionRight:DeltaOnly No set MEDiffractionDouble:DeltaOnly No set MEDiffractionDeltaLeft:DeltaOnly Yes set MEDiffractionDeltaRight:DeltaOnly Yes set MEDiffractionDeltaDouble:DeltaOnly Yes # Set weight for Diffraction set MEDiffractionLeft:DiffractionAmplitude 12 set MEDiffractionRight:DiffractionAmplitude 12 set MEDiffractionDouble:DiffractionAmplitude 8 set MEDiffractionDeltaLeft:DiffractionAmplitude 4 set MEDiffractionDeltaRight:DiffractionAmplitude 4 set MEDiffractionDeltaDouble:DiffractionAmplitude 2 # Set soft diffraction paramters # Parameter values from arxiv/0709.0395 set MEDiffractionLeft:ProtonPomeronSlope 10.1 set MEDiffractionLeft:SoftPomeronIntercept 1.08 set MEDiffractionLeft:SoftPomeronSlope 0.25 set MEDiffractionRight:ProtonPomeronSlope 10.1 set MEDiffractionRight:SoftPomeronIntercept 1.08 set MEDiffractionRight:SoftPomeronSlope 0.25 set MEDiffractionDouble:ProtonPomeronSlope 10.1 set MEDiffractionDouble:SoftPomeronIntercept 1.08 set MEDiffractionDouble:SoftPomeronSlope 0.25 set MEDiffractionDeltaLeft:ProtonPomeronSlope 10.1 set MEDiffractionDeltaLeft:SoftPomeronIntercept 1.08 set MEDiffractionDeltaLeft:SoftPomeronSlope 0.25 set MEDiffractionDeltaRight:ProtonPomeronSlope 10.1 set MEDiffractionDeltaRight:SoftPomeronIntercept 1.08 set MEDiffractionDeltaRight:SoftPomeronSlope 0.25 set MEDiffractionDeltaDouble:ProtonPomeronSlope 10.1 set MEDiffractionDeltaDouble:SoftPomeronIntercept 1.08 set MEDiffractionDeltaDouble:SoftPomeronSlope 0.25 # Set number of clusters for dissociation set MEDiffractionLeft:DissociationDecay One set MEDiffractionRight:DissociationDecay One set MEDiffractionDouble:DissociationDecay One set MEDiffractionDeltaLeft:DissociationDecay One set MEDiffractionDeltaRight:DissociationDecay One set MEDiffractionDeltaDouble:DissociationDecay One # Insert matrix elements insert QCDDiffraction:MatrixElements[0] MEDiffractionLeft insert QCDDiffraction:MatrixElements[0] MEDiffractionRight insert QCDDiffraction:MatrixElements[0] MEDiffractionDouble insert QCDDiffraction:MatrixElements[0] MEDiffractionDeltaLeft insert QCDDiffraction:MatrixElements[0] MEDiffractionDeltaRight #insert QCDDiffraction:MatrixElements[0] MEDiffractionDeltaDouble # No cluster fission (set to E/2, if needed. Here E=7TeV.) #set /Herwig/Hadronization/ClusterFissioner:ClMaxLight 3500 + +# The DiffractionMEs need to know the MPIHandler for the cross section. +set MEDiffractionLeft:MPIHandler /Herwig/UnderlyingEvent/MPIHandler +set MEDiffractionRight:MPIHandler /Herwig/UnderlyingEvent/MPIHandler +set MEDiffractionDouble:MPIHandler /Herwig/UnderlyingEvent/MPIHandler +set MEDiffractionDeltaLeft:MPIHandler /Herwig/UnderlyingEvent/MPIHandler +set MEDiffractionDeltaRight:MPIHandler /Herwig/UnderlyingEvent/MPIHandler +set MEDiffractionDeltaDouble:MPIHandler /Herwig/UnderlyingEvent/MPIHandler + + set QCDDiffraction:CascadeHandler NULL set /Herwig/Generators/EventGenerator:EventHandler:CascadeHandler NULL insert /Herwig/Generators/EventGenerator:EventHandler:SubProcessHandlers[0] QCDDiffraction ########################################################################## ########################################################################## diff --git a/src/snippets/MB-DipoleShower.in b/src/snippets/MB-DipoleShower.in --- a/src/snippets/MB-DipoleShower.in +++ b/src/snippets/MB-DipoleShower.in @@ -1,47 +1,54 @@ ################################################## # MEMinBias Matrix Element ################################################## ### Note to users - Release 7.1: ### This currently uses parameters tuned for the ### default shower. # MPI model settings set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0 ## Report the correct cross section cd /Herwig/Generators -create Herwig::MPIXSecReweighter MPIXSecReweighter -insert EventGenerator:EventHandler:PostSubProcessHandlers 0 MPIXSecReweighter + +# TODO: remove the MPIXSecReweighter once validated. +# create Herwig::MPIXSecReweighter MPIXSecReweighter +# insert EventGenerator:EventHandler:PostSubProcessHandlers 0 MPIXSecReweighter set EventGenerator:EventHandler:CascadeHandler NULL clear EventGenerator:EventHandler:SubProcessHandlers[0] ################################################## # Create separate SubProcessHandler for MinBias ################################################## cd /Herwig/MatrixElements/ cp SubProcess QCDMinBias set QCDMinBias:CascadeHandler /Herwig/DipoleShower/DipoleShowerHandler set QCDMinBias:CascadeHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler set QCDMinBias:DecayHandler /Herwig/Decays/DecayHandler # Due to numerics the pomeron could be seen as timelike. set /Herwig/Shower/ShowerHandler:SplitHardProcess No set /Herwig/DipoleShower/DipoleShowerHandler:SplitHardProcess No +set MEMinBias:MPIHandler /Herwig/UnderlyingEvent/MPIHandler + insert QCDMinBias:MatrixElements[0] MEMinBias cd /Herwig/Generators # MinBias parameters used for the new kinematics of soft MPI -set /Herwig/Cuts/MinBiasCuts:X1Min 0.11 -set /Herwig/Cuts/MinBiasCuts:X2Min 0.11 +set /Herwig/Cuts/MinBiasCuts:X1Min 0.011 +set /Herwig/Cuts/MinBiasCuts:X2Min 0.011 + # Needed to get the correct fraction of diffractive events -set /Herwig/MatrixElements/MEMinBias:csNorm 4.5584 +set /Herwig/MatrixElements/MEMinBias:csNorm 0.01 +set /Herwig/MatrixElements/MEMinBias:Scale 2.0 + set EventGenerator:EventHandler:Cuts /Herwig/Cuts/MinBiasCuts cd /Herwig/MatrixElements/ insert /Herwig/Generators/EventGenerator:EventHandler:SubProcessHandlers[0] QCDMinBias diff --git a/src/snippets/MB.in b/src/snippets/MB.in --- a/src/snippets/MB.in +++ b/src/snippets/MB.in @@ -1,55 +1,63 @@ ################################################## # MEMinBias Matrix Element ################################################## # MPI model settings set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0 ## Report the correct cross section cd /Herwig/Generators -create Herwig::MPIXSecReweighter MPIXSecReweighter -insert EventGenerator:EventHandler:PostSubProcessHandlers 0 MPIXSecReweighter + +# TODO: remove the MPIXSecReweighter once validated. +#create Herwig::MPIXSecReweighter MPIXSecReweighter +#insert EventGenerator:EventHandler:PostSubProcessHandlers 0 MPIXSecReweighter + + set EventGenerator:EventHandler:CascadeHandler NULL clear EventGenerator:EventHandler:SubProcessHandlers[0] ################################################## # Create separate SubProcessHandler for MinBias ################################################## cd /Herwig/MatrixElements/ cp SubProcess QCDMinBias set QCDMinBias:CascadeHandler /Herwig/Shower/ShowerHandler set QCDMinBias:CascadeHandler:MPIHandler /Herwig/UnderlyingEvent/MPIHandler set QCDMinBias:DecayHandler /Herwig/Decays/DecayHandler # Due to numerics the pomeron could be seen as timelike. set /Herwig/Shower/ShowerHandler:SplitHardProcess No set /Herwig/DipoleShower/DipoleShowerHandler:SplitHardProcess No insert QCDMinBias:MatrixElements[0] MEMinBias +# MEMinBias is automatically reweighted to give the cross section +# comunicated by the MPIHandler. +set MEMinBias:MPIHandler /Herwig/UnderlyingEvent/MPIHandler + cd /Herwig/Generators # MinBias parameters used for the new kinematics of soft MPI set /Herwig/Cuts/MinBiasCuts:X1Min 0.011 set /Herwig/Cuts/MinBiasCuts:X2Min 0.011 #PDFs for MPI,Underlying Event, Min Bias # set the correct PDFs set /Herwig/Partons/HardLOPDF:PDFName CT14lo set /Herwig/Partons/ShowerLOPDF:PDFName CT14lo set /Herwig/Partons/MPIPDF:PDFName CT14lo set /Herwig/Partons/RemnantPDF:PDFName CT14lo # Needed to get the correct fraction of diffractive events set /Herwig/MatrixElements/MEMinBias:csNorm 0.01 set /Herwig/MatrixElements/MEMinBias:Scale 2.0 set EventGenerator:EventHandler:Cuts /Herwig/Cuts/MinBiasCuts cd /Herwig/MatrixElements/ insert /Herwig/Generators/EventGenerator:EventHandler:SubProcessHandlers[0] QCDMinBias