diff --git a/MatrixElement/Hadron/MEDiffraction.cc b/MatrixElement/Hadron/MEDiffraction.cc --- a/MatrixElement/Hadron/MEDiffraction.cc +++ b/MatrixElement/Hadron/MEDiffraction.cc @@ -1,958 +1,1026 @@ // -*- 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/Utilities/UnitIO.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), + 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+ - + 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 + + //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 + 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 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 + //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))); + 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, + 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, + 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, + 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, + 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; - } + } + 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()); + 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 + + //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: + 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; + 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] = p4; - meMomenta()[3] = decayMomenta. first; - meMomenta()[4] = decayMomenta.second; - break; + meMomenta()[4].setT(sqrt(mq2+sqr(meMomenta()[4].x())+sqr(meMomenta()[4].y())+sqr(meMomenta()[4].z()))); + //////////////////////////////////////////////////// + // Proton decays along z axis in rest frame + if (smearAngle == 0){ + do{} + while(!Kinematics::twoBodyDecay(p3,mqq(),mq(),-dir,decayMomenta.first,decayMomenta.second)); + } + // Otherwise it has randomly sampled azimuth and polar angles. + else{ + do{} + while(!twoBodyDecayMomenta(p3,mqq(),mq(),decayMomenta.first,decayMomenta.second)); + } + /////////// + + meMomenta()[2] = p4; + meMomenta()[3] = decayMomenta. first; + meMomenta()[4] = decayMomenta.second; + 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] = p3; - meMomenta()[3] = decayMomenta. first; - meMomenta()[4] = decayMomenta.second; - break; + meMomenta()[4].setT(sqrt(mq2+sqr(meMomenta()[4].x())+sqr(meMomenta()[4].y())+sqr(meMomenta()[4].z()))); + //////////////////////////////////////////////////// + // Proton decays along z axis in rest frame + if (smearAngle == 0){ + do{} + while(!Kinematics::twoBodyDecay(p4,mqq(),mq(),dir,decayMomenta.first,decayMomenta.second)); + } + // Otherwise it has randomly sampled azimuth and polar angles. + else{ + do{} + while(!twoBodyDecayMomenta(p4,mqq(),mq(),decayMomenta.first,decayMomenta.second)); + } + meMomenta()[2] = p3; + meMomenta()[3] = decayMomenta. first; + meMomenta()[4] = decayMomenta.second; + 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] = decayMomenta. first; - meMomenta()[3] = decayMomenta.second; - meMomenta()[4] = decayMomentaTwo.second; - meMomenta()[5] = decayMomentaTwo. first; - break; + // Proton decays along z axis in rest frame + if (smearAngle == 0){ + do{} + while(!Kinematics::twoBodyDecay(p3,mqq(),mq(),-dir,decayMomenta.first,decayMomenta.second)); + + do{} + while(!Kinematics::twoBodyDecay(p4,mqq(),mq(),dir,decayMomentaTwo.first,decayMomentaTwo.second)); + } + // Otherwise it has randomly sampled azimuth and polar angles. + else{ + do{} + while(!twoBodyDecayMomenta(p3,mqq(),mq(),decayMomenta.first,decayMomenta.second)); + + do{} + while(!twoBodyDecayMomenta(p4,mqq(),mq(),decayMomentaTwo.first,decayMomentaTwo.second)); + } + + meMomenta()[2] = decayMomenta. first; + meMomenta()[3] = decayMomenta.second; + meMomenta()[4] = decayMomentaTwo.second; + meMomenta()[5] = decayMomentaTwo. first; + break; } - + } else { meMomenta()[2] = p3; meMomenta()[3] = p4; } break; case 1: switch(diffDirection){ - case 0: + 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 + //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: + 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 + + //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: + 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)); - + 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) +//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 + 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); } +// Implementation of two body decay with random angle chosen +// Used for the decay of the excited proton into quark and diquark +bool MEDiffraction::twoBodyDecayMomenta(const Lorentz5Momentum & pIn, + const Energy m1, const Energy m2, + Lorentz5Momentum & p1, Lorentz5Momentum & p2){ + Energy min = pIn.mass(); + // Check that the decay is kinematically possible + if (min >= m1 + m2 && m1 >= ZERO && m2 >= ZERO){ + Energy2 M2 = sqr(min); + // pstar squared + const Energy2 psq = ((M2-sqr(m1+m2)*(M2-sqr(m1-m2))/(4*M2); + // Shouldn't have to check (since the ), but no harm + assert(psq/GeV2 > 0); + const Energy p(sqrt(psq)); + // Sample the angle + const double phi = UseRandom::rnd() * Constants::twopi; + const double costheta =1-2*UseRandom::rnd(); + const double sintheta = sqrt(1-sqr(costheta)); + + // Generate the momenta + Lorentz5Momentum k1=Lorentz5Momentum(p*sintheta*cos(phi), p*sintheta*sin(phi), p*costheta, sqrt(sqr(m1)+psq)); + Lorentz5Momentum k2=Lorentz5Momentum(-p*sintheta*cos(phi), -p*sintheta*sin(phi), -p*costheta,sqrt(sqr(m2)+psq)); + + //find boost to pp center of mass + const Boost betap3 = (pIn).findBoostToCM(); + + //k1 and k2 calculated at p3 center of mass, so boost back + k1.boost(-betap3); + k2.boost(-betap3); + p1 = k1; + p2 = k2; + return true; + } + return false; + } +/* //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); + + 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. + // 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. + // 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; - + // Cross section without reweighting: XS_norew // XS_have = avcsNorm2*XS_norew (for large N) // We want to determine the rew that allows to get the wanted XS. // In deltaN points we want (left) and we get (right): // XS_wanted*(countN+deltaN) = XS_have*countN + rew*deltaN*XS_norew // Solve for rew: - + rew_=avRew*(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. + //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; } CrossSection MEDiffraction::dSigHatDR() const { // Apply the correction weight switch if turned on if ( correctionWeightSwitch ) { return me2()*jacobian()/sHat()*sqr(hbarc)*correctionweight(); } else { // for single diffraction return me2()*jacobian()/sHat()*sqr(hbarc); } } 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) + if(diags[0]->id()==-1) sel.insert(2./3.,i); else - sel.insert(1./3.,i); + sel.insert(1./3.,i); } - + }else{ for(unsigned int i = 0; i < diags.size(); i++){ - if(diags[0]->id()==-1) + if(diags[0]->id()==-1) sel.insert(4./9.,i); else if(diags[0]->id()==-2) - sel.insert(2./9.,i); + sel.insert(2./9.,i); else if(diags[0]->id()==-3) sel.insert(2./9.,i); else - sel.insert(1./9.,i); + 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); + + sel.insert(1.0, &cl); } break; case 1: switch(diffDirection){ - case 0: + case 0: static ColourLines clleft("-6 2 3 8, -8 -3 7"); sel.insert(1.0, &clleft); break; - case 1: + case 1: static ColourLines clright("-9 4 3 7, -7 -3 8"); sel.insert(1.0, &clright); break; - case 2: + 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) << MPIHandler_ << correctionWeightSwitch - << thefixedM2min; + << thefixedM2min << smearAngle; } void MEDiffraction::persistentInput(PersistentIStream & is, int) { is >> theme2 >> deltaOnly >> diffDirection >> theprotonPomeronSlope >> thesoftPomeronIntercept >> thesoftPomeronSlope >> dissociationDecay >> iunit(theProtonMass,GeV)>> MPIHandler_ >> correctionWeightSwitch - >> thefixedM2min; + >> thefixedM2min >> smearAngle; } InvEnergy2 MEDiffraction::protonPomeronSlope() const{ return theprotonPomeronSlope/GeV2; } double MEDiffraction::softPomeronIntercept() const { return thesoftPomeronIntercept; } InvEnergy2 MEDiffraction::softPomeronSlope() const { return thesoftPomeronSlope/GeV2; } // return fixedM2min if in allowed range, else return kinematic limits Energy2 MEDiffraction::M2min() const{ if (fixedM2min() < sqr(getParticleData(2212)->mass()+mq()+mqq())) return sqr(getParticleData(2212)->mass()+mq()+mqq()); else return fixedM2min() < M2max() ? fixedM2min(): M2max(); } Energy2 MEDiffraction::M2max() const{ return sqr(generator()->maximumCMEnergy()-getParticleData(2212)->mass()); } Energy2 MEDiffraction::fixedM2min() const { return thefixedM2min*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); - + false, false, Interface::limited); + static Parameter interfacesoftPomeronIntercept ("SoftPomeronIntercept", "The soft pomeron intercept.", &MEDiffraction::thesoftPomeronIntercept, 1.08, 0.00001, 100.0, - false, false, Interface::limited); - + 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); - - + 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); - + (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); static Switch interfaceCorrectionWeightSwitch ("CorrectionWeight", "use correction weight for diffraction crosssection", &MEDiffraction::correctionWeightSwitch, 1, false, false); static SwitchOption interfaceCorrectionWeight0 (interfaceCorrectionWeightSwitch,"No","Correction weight is not used", 0); static SwitchOption interfaceCorrectionWeight1 (interfaceCorrectionWeightSwitch,"Yes","Correction weight is used", 1); static Parameter interfacesoftM2min ("MinDiffMass", "Manually set the minimum diffractive mass(GeV^2), bound by kinematics", &MEDiffraction::thefixedM2min, 0.0, 0.0, 100.0, false, false, Interface::nolimits); + + static Switch interfaceSmearAngle + ("SmearAngle", + "Smear the angle of the 1 to 2 decay in the rest frame of the decaying diffractive mass", + &MEDiffraction::smearAngle, 0, false, false); + static SwitchOption interfaceSmearAngle0 + (interfaceSmeaAngle,"No","Decay along z axis (in rest frame)", 0); + static SwitchOption interfaceSmearAngle1 + (interfaceSmearAngle,"Yes","Sample theta and phi (in rest frame)", 1); } diff --git a/MatrixElement/Hadron/MEDiffraction.h b/MatrixElement/Hadron/MEDiffraction.h --- a/MatrixElement/Hadron/MEDiffraction.h +++ b/MatrixElement/Hadron/MEDiffraction.h @@ -1,324 +1,330 @@ // -*- 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 + * 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; + bool deltaOnly; /* use correction weight for diffraction compared to inelastic weights */ bool correctionWeightSwitch; + /* allow for non-zero azimuthal and polar angles in the 1 to 2 decay */ + bool smearAngle; + /* manual set min mass */ double thefixedM2min; - + /* Direction of the excited proton */ - unsigned int diffDirection; - + unsigned int diffDirection; + /* Number of clusters the dissociated proton decays into */ - unsigned int dissociationDecay; - + 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; + //pair twoBodyDecayMomenta(Lorentz5Momentum pp) const; + bool twoBodyDecayMomenta(const Lorentz5Momentum & p, + const Energy m1, const Energy m2, + Lorentz5Momentum & p1, Lorentz5Momentum & p2); /** * Returns the proton-pomeron slope */ - InvEnergy2 protonPomeronSlope() const; - + InvEnergy2 protonPomeronSlope() const; + /** * Returns the soft pomeron intercept - */ - double softPomeronIntercept() const; - - //M12 and M22 are masses squared of + */ + 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; + 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; + Energy2 M2min() const; /** * Returns the maximal possible value of diffractive mass */ Energy2 M2max() const; //TODO:modify to get proper parameters InvEnergy2 softPomeronSlope() const; /** * Allows for a fixed m2 value */ Energy2 fixedM2min() 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; - - + bool isInRunPhase; + + /* The proton mass */ Energy theProtonMass; /** - * a MPIHandler to administer the creation of several (semihard) + * a MPIHandler to administer the creation of several (semihard) * partonic interactions. */ UEBasePtr MPIHandler_; - + }; } #endif /* HERWIG_MEDiffraction_H */