Page MenuHomeHEPForge
Files F8310550

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/Hadron/MEDiffraction.cc b/MatrixElement/Hadron/MEDiffraction.cc
--- a/MatrixElement/Hadron/MEDiffraction.cc
+++ b/MatrixElement/Hadron/MEDiffraction.cc
@@ -1,851 +1,853 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEDiffraction class.
//
#include "MEDiffraction.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Handlers/StandardXComb.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) {}
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, -3)));
+ 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, -3)));
+ 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, -3)));
+ 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, -7)));
- //d -- uu_1 left u -- ud_1 right
- //add(new_ptr((Tree2toNDiagram(5), prt1, q2, pom, q1, prt1, 1, dq2, 2, q2, 3, q1, 4, dq11, -8)));
+ 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, -9)));
+ 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<pair<Energy2,Energy2>,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( getParticleData(2212)->mass() );
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 mqq2(sqr(mqq())), mq2(sqr(mq()));
- double sintheta = sqrt(1 - sqr(costheta));
- Energy pqq;
+
Energy2 Mx2;
switch(diffDirection){
case 0:
Mx2=M12;
break;
case 1:
Mx2=M22;
break;
}
- //total momentum
- pqq = ((Mx2-mq2+mqq2)*pprime*costheta+ sqrt((Mx2+sqr(pprime))*(kallen(Mx2,mq2,mqq2)
- -4*mqq2*sqr(pprime)*sqr(sintheta))))/(2*Mx2+2*sqr(pprime)*sqr(sintheta));
+
/* Select between left/right single diffraction and double diffraction */
//check if we want only delta for the excited state
- //new boost vectors
- //const Boost p3boost = p3.findBoostToCM();
- //const Boost p4boost = p4.findBoostToCM();
- //Axis norm;
//pair of momenta for double decay for a two cluster case
pair<Lorentz5Momentum,Lorentz5Momentum> momPair, momPair1;
//fraction of momenta
double frac = UseRandom::rnd();
switch(dissociationDecay){
case 0:
if(!deltaOnly)
{
pair<Lorentz5Momentum,Lorentz5Momentum> decayMomenta;
pair<Lorentz5Momentum,Lorentz5Momentum> decayMomentaTwo;
const double phiprime = UseRandom::rnd() * Constants::twopi;
//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
{
meMomenta()[2+diffDirection].setVect(p3.vect());
meMomenta()[2+diffDirection].setT(p3.t());
meMomenta()[3-diffDirection].setVect(p4.vect());
meMomenta()[3-diffDirection].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<pair<Energy2,Energy2>,Energy2> MEDiffraction::diffractiveMassAndMomentumTransfer() const {
Energy2 theM12(ZERO),theM22(ZERO), thet(ZERO);
int count = 0;
//proton mass squared
const Energy2 m2 = sqr(getParticleData(2212)->mass());
//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;
thet = randomt(md2);
break;
case 1:
theM22 = md2;
theM12 = m2;
M2 = md2;
thet = randomt(md2);
break;
case 2:
theM12 = md2;
theM22 = md2;
M2 = md2;
thet = doublediffrandomt(theM12,theM22);
break;
}
}else{
switch (diffDirection){
case 0:
M2=randomM2();
thet = randomt(M2);
theM12=M2;
theM22=m2;
break;
case 1:
theM12=m2;
M2=randomM2();
thet = randomt(M2);
theM22=M2;
break;
case 2:
theM12=randomM2();
theM22=randomM2();
M2=(theM12>theM22) ? theM12: theM22;
thet = doublediffrandomt(theM12,theM22);
break;
}
}
count++;
const Energy cmEnergy = generator()->maximumCMEnergy();
const Energy2 s = sqr(cmEnergy);
InvEnergy2 slope;
if(diffDirection==2){
slope = 2*softPomeronSlope()*log(.1+(sqr(cmEnergy)/softPomeronSlope())/(theM12*theM22));
}else{
slope = protonPomeronSlope()
+ 2*softPomeronSlope()*log(sqr(cmEnergy)/M2);
}
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()*GeV2)*(expmax-expmin)/(slope*GeV2))
||((theM12/GeV2)*(theM22/GeV2)>=(sqr(cmEnergy)/GeV2)/(softPomeronSlope()*GeV2));
}while(condition);
return make_pair (make_pair(theM12,theM22),thet);
}
//Decay of the excited proton to quark-diquark
pair<Lorentz5Momentum,Lorentz5Momentum> 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( getParticleData(2212)->mass() );
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);
return log( exp(slope*ttmin) +
UseRandom::rnd()*(exp(slope*ttmax) - exp(slope*ttmin)) ) / slope;
}
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);
return log( exp(slope*ttmin) +
UseRandom::rnd()*(exp(slope*ttmax) - exp(slope*ttmin)) ) / slope;
}
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( getParticleData(2212)->mass() );
return (1/(2*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( getParticleData(2212)->mass() );
return (1/(2*s))*(sqrt(kallen(s, m2, m2)*kallen(s, M12, M22))-sqr(s)+2*s*m2+s*M12+s*M22);
}
double MEDiffraction::me2() const{
return theme2;
}
CrossSection MEDiffraction::dSigHatDR() const {
return me2()*jacobian()/sHat()*sqr(hbarc);
}
unsigned int MEDiffraction::orderInAlphaS() const {
return 0;
}
unsigned int MEDiffraction::orderInAlphaEW() const {
return 0;
}
Selector<MEBase::DiagramIndex>
-MEDiffraction::diagrams(const DiagramVector & ) const {
+MEDiffraction::diagrams(const DiagramVector & diags) const {
Selector<DiagramIndex> sel;
-
if(!deltaOnly){
if(diffDirection<2){
- sel.insert(2/3,0);
- sel.insert(1/3,1);
+
+ 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{
-
- sel.insert(4/9,0);
- sel.insert(2/9,1);
- sel.insert(2/9,2);
- sel.insert(1/9,3);
+ 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<const ColourLines *>
MEDiffraction::colourGeometries(tcDiagPtr ) const {
Selector<const ColourLines *> 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::doinitrun() {
HwMEBase::doinitrun();
isInRunPhase = true;
}
IBPtr MEDiffraction::clone() const {
return new_ptr(*this);
}
IBPtr MEDiffraction::fullclone() const {
return new_ptr(*this);
}
ClassDescription<MEDiffraction> MEDiffraction::initMEDiffraction;
// Definition of the static class description member.
void MEDiffraction::persistentOutput(PersistentOStream & os) const {
os << theme2 << deltaOnly << diffDirection << theprotonPomeronSlope
<< thesoftPomeronIntercept << thesoftPomeronSlope << dissociationDecay;
}
void MEDiffraction::persistentInput(PersistentIStream & is, int) {
is >> theme2 >> deltaOnly >> diffDirection >> theprotonPomeronSlope
>> thesoftPomeronIntercept >> thesoftPomeronSlope >> dissociationDecay;
}
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<MEDiffraction> documentation
("There is no documentation for the MEDiffraction class");
static Parameter<MEDiffraction,double> 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<MEDiffraction,double> interfaceprotonPomeronSlope
("ProtonPomeronSlope",
"The proton-pomeron slope parameter.",
&MEDiffraction::theprotonPomeronSlope, 10.1, 0.00001, 100.0,
false, false, Interface::limited);
static Parameter<MEDiffraction,double> interfacesoftPomeronIntercept
("SoftPomeronIntercept",
"The soft pomeron intercept.",
&MEDiffraction::thesoftPomeronIntercept, 1.08, 0.00001, 100.0,
false, false, Interface::limited);
static Parameter<MEDiffraction,double> interfacesoftPomeronSlope
("SoftPomeronSlope",
"The soft pomeron slope parameter.",
&MEDiffraction::thesoftPomeronSlope, 0.25, 0.00001, 100.0,
false, false, Interface::limited);
static Switch<MEDiffraction, bool> 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<MEDiffraction, unsigned int> 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<MEDiffraction, unsigned int> 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);
}
diff --git a/MatrixElement/Hadron/MEDiffraction.h b/MatrixElement/Hadron/MEDiffraction.h
--- a/MatrixElement/Hadron/MEDiffraction.h
+++ b/MatrixElement/Hadron/MEDiffraction.h
@@ -1,334 +1,336 @@
// -*- C++ -*-
#ifndef HERWIG_MEDiffraction_H
#define HERWIG_MEDiffraction_H
//
// This is the declaration of the MEDiffraction class.
//
#include "Herwig/MatrixElement/HwMEBase.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;
/**
* 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<DiagramIndex> 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<const ColourLines *>
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:
/**
* 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;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
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<pair<Energy2,Energy2>,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<Lorentz5Momentum,Lorentz5Momentum> 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(0.938+0.325+2*0.65)*GeV2;}
+ 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()-1*GeV);}//TODO:modify to get proper parameters
+ Energy2 M2max() const{
+ return sqr(generator()->maximumCMEnergy()-getParticleData(2212)->mass());
+ }//TODO:modify to get proper parameters
InvEnergy2 softPomeronSlope() const;
/* Kallen function */
template<int L, int E, int Q, int DL, int DE, int DQ>
Qty<2*L,2*E,2*Q,DL,DE,DQ> kallen(Qty<L,E,Q,DL,DE,DQ> a,
Qty<L,E,Q,DL,DE,DQ> b,
Qty<L,E,Q,DL,DE,DQ> c) const {
return a*a + b*b + c*c - 2.0*(a*b + b*c + c*a);
}
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<MEDiffraction> initMEDiffraction;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MEDiffraction & operator=(const MEDiffraction &);
bool isInRunPhase;
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of MEDiffraction. */
template <>
struct BaseClassTrait<Herwig::MEDiffraction,1> {
/** Typedef of the first base class of MEDiffraction. */
typedef Herwig::HwMEBase NthBase;
};
/** This template specialization informs ThePEG about the name of
* the MEDiffraction class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::MEDiffraction>
: public ClassTraitsBase<Herwig::MEDiffraction> {
/** Return a platform-independent class name */
static string className() { return "Herwig::MEDiffraction"; }
/**
* The name of a file containing the dynamic library where the class
* MEDiffraction is implemented. It may also include several, space-separated,
* libraries if the class MEDiffraction depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() {return "HwMEHadron.so";}
};
/** @endcond */
}
#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,152 +1,164 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEMinBias class.
//
#include "MEMinBias.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/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(10*GeV);
}
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 ) {
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::me2() const {
//tuned so it gives the correct normalization for xmin = 0.11
- return 4.5584*(sqr(generator()->maximumCMEnergy())/GeV2);
+ return csNorm_*(sqr(generator()->maximumCMEnergy())/GeV2);
}
CrossSection MEMinBias::dSigHatDR() const {
return me2()*jacobian()/sHat()*sqr(hbarc);
}
unsigned int MEMinBias::orderInAlphaS() const {
return 2;
}
unsigned int MEMinBias::orderInAlphaEW() const {
return 0;
}
Selector<MEBase::DiagramIndex>
MEMinBias::diagrams(const DiagramVector & diags) const {
Selector<DiagramIndex> sel;
for ( DiagramIndex i = 0; i < diags.size(); ++i )
sel.insert(1.0, i);
return sel;
}
Selector<const ColourLines *>
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<const ColourLines *> 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);
}
ClassDescription<MEMinBias> MEMinBias::initMEMinBias;
// Definition of the static class description member.
+void MEMinBias::persistentOutput(PersistentOStream & os) const {
+ os << csNorm_;
+}
+void MEMinBias::persistentInput(PersistentIStream & is, int) {
+ is >> csNorm_;
+}
void MEMinBias::Init() {
static ClassDocumentation<MEMinBias> documentation
("There is no documentation for the MEMinBias class");
+ static Parameter<MEMinBias,double> interfacecsNorm
+ ("csNorm",
+ "Normalization of the min-bias cross section.",
+ &MEMinBias::csNorm_,
+ 1.0, 0.0, 100.0,
+ false, false, Interface::limited);
}
diff --git a/MatrixElement/Hadron/MEMinBias.h b/MatrixElement/Hadron/MEMinBias.h
--- a/MatrixElement/Hadron/MEMinBias.h
+++ b/MatrixElement/Hadron/MEMinBias.h
@@ -1,199 +1,224 @@
// -*- C++ -*-
#ifndef HERWIG_MEMinBias_H
#define HERWIG_MEMinBias_H
//
// This is the declaration of the MEMinBias class.
//
#include "Herwig/MatrixElement/HwMEBase.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:
/** @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;
/**
* 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<DiagramIndex> 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<const ColourLines *>
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
+ */
+ double csNorm_;
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<MEMinBias> initMEMinBias;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MEMinBias & operator=(const MEMinBias &);
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of MEMinBias. */
template <>
struct BaseClassTrait<Herwig::MEMinBias,1> {
/** Typedef of the first base class of MEMinBias. */
typedef Herwig::HwMEBase NthBase;
};
/** This template specialization informs ThePEG about the name of
* the MEMinBias class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::MEMinBias>
: public ClassTraitsBase<Herwig::MEMinBias> {
/** Return a platform-independent class name */
static string className() { return "Herwig::MEMinBias"; }
/**
* The name of a file containing the dynamic library where the class
* MEMinBias is implemented. It may also include several, space-separated,
* libraries if the class MEMinBias depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() {return "HwMEHadron.so";}
};
/** @endcond */
}
#endif /* HERWIG_MEMinBias_H */
diff --git a/PDF/HwRemDecayer.cc b/PDF/HwRemDecayer.cc
--- a/PDF/HwRemDecayer.cc
+++ b/PDF/HwRemDecayer.cc
@@ -1,1512 +1,1514 @@
// -*- C++ -*-
//
// HwRemDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig is licenced under version 2 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 HwRemDecayer class.
//
#include "HwRemDecayer.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Utilities/UtilityBase.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
#include "ThePEG/Utilities/Throw.h"
#include "Herwig/Shower/ShowerHandler.h"
using namespace Herwig;
namespace{
const bool dbg = false;
void reShuffle(Lorentz5Momentum &p1, Lorentz5Momentum &p2, Energy m1, Energy m2){
Lorentz5Momentum ptotal(p1+p2);
ptotal.rescaleMass();
if( ptotal.m() < m1+m2 ) {
if(dbg)
cerr << "Not enough energy to perform reshuffling \n";
throw HwRemDecayer::ExtraSoftScatterVeto();
}
Boost boostv = -ptotal.boostVector();
ptotal.boost(boostv);
p1.boost(boostv);
// set the masses and energies,
p1.setMass(m1);
p1.setE(0.5/ptotal.m()*(ptotal.m2()+sqr(m1)-sqr(m2)));
p1.rescaleRho();
// boost back to the lab
p1.boost(-boostv);
p2.boost(boostv);
// set the masses and energies,
p2.setMass(m2);
p2.setE(0.5/ptotal.m()*(ptotal.m2()+sqr(m2)-sqr(m1)));
p2.rescaleRho();
// boost back to the lab
p2.boost(-boostv);
}
}
void HwRemDecayer::initialize(pair<tRemPPtr, tRemPPtr> rems, tPPair beam, Step & step,
Energy forcedSplitScale) {
// the step
thestep = &step;
// valence content of the hadrons
theContent.first = getHadronContent(beam.first);
theContent.second = getHadronContent(beam.second);
// momentum extracted from the hadrons
theUsed.first = Lorentz5Momentum();
theUsed.second = Lorentz5Momentum();
theMaps.first.clear();
theMaps.second.clear();
theX.first = 0.0;
theX.second = 0.0;
theRems = rems;
_forcedSplitScale = forcedSplitScale;
// check remnants attached to the right hadrons
if( (theRems.first && parent(theRems.first ) != beam.first ) ||
(theRems.second && parent(theRems.second) != beam.second) )
throw Exception() << "Remnant order wrong in "
<< "HwRemDecayer::initialize(...)"
<< Exception::runerror;
return;
}
void HwRemDecayer::split(tPPtr parton, HadronContent & content,
tRemPPtr rem, Lorentz5Momentum & used,
PartnerMap &partners, tcPDFPtr pdf, bool first) {
theBeam = parent(rem);
theBeamData = dynamic_ptr_cast<Ptr<BeamParticleData>::const_pointer>
(theBeam->dataPtr());
double currentx = parton->momentum().rho()/theBeam->momentum().rho();
double check = rem==theRems.first ? theX.first : theX.second;
check += currentx;
if(1.0-check < 1e-3) throw ShowerHandler::ExtraScatterVeto();
bool anti;
Lorentz5Momentum lastp(parton->momentum());
int lastID(parton->id());
Energy oldQ(_forcedSplitScale);
_pdf = pdf;
//do nothing if already valence quark
if(first && content.isValenceQuark(parton)) {
//set the extracted value, because otherwise no RemID could be generated.
content.extract(lastID);
// add the particle to the colour partners
partners.push_back(make_pair(parton, tPPtr()));
//set the sign
anti = parton->hasAntiColour() && parton->id()!=ParticleID::g;
if(rem==theRems.first) theanti.first = anti;
else theanti.second = anti;
// add the x and return
if(rem==theRems.first) theX.first += currentx;
else theX.second += currentx;
return;
}
//or gluon for secondaries
else if(!first && lastID == ParticleID::g) {
partners.push_back(make_pair(parton, tPPtr()));
// add the x and return
if(rem==theRems.first) theX.first += currentx;
else theX.second += currentx;
return;
}
// if a sea quark.antiquark forced splitting to a gluon
// Create the new parton with its momentum and parent/child relationship set
PPtr newSea;
if( !(lastID == ParticleID::g ||
lastID == ParticleID::gamma) ) {
newSea = forceSplit(rem, -lastID, oldQ, currentx, lastp, used,content);
ColinePtr cl = new_ptr(ColourLine());
if(newSea->id() > 0) cl-> addColoured(newSea);
else cl->addAntiColoured(newSea);
// if a secondard scatter finished so return
if(!first || content.isValenceQuark(ParticleID::g) ){
partners.push_back(make_pair(parton, newSea));
// add the x and return
if(rem==theRems.first) theX.first += currentx;
else theX.second += currentx;
if(first) content.extract(ParticleID::g);
return;
}
}
// otherwise evolve back to valence
// final valence splitting
PPtr newValence = forceSplit(rem,
lastID!=ParticleID::gamma ?
ParticleID::g : ParticleID::gamma,
oldQ, currentx , lastp, used, content);
// extract from the hadron to allow remnant to be determined
content.extract(newValence->id());
// case of a gluon going into the hard subprocess
if( lastID == ParticleID::g ) {
partners.push_back(make_pair(parton, tPPtr()));
anti = newValence->hasAntiColour();
if(rem==theRems.first) theanti.first = anti;
else theanti.second = anti;
parton->colourLine(!anti)->addColoured(newValence, anti);
return;
}
else if( lastID == ParticleID::gamma) {
partners.push_back(make_pair(parton, newValence));
anti = newValence->hasAntiColour();
ColinePtr newLine(new_ptr(ColourLine()));
newLine->addColoured(newValence, anti);
if(rem==theRems.first) theanti.first = anti;
else theanti.second = anti;
// add the x and return
if(rem==theRems.first) theX.first += currentx;
else theX.second += currentx;
return;
}
//The valence quark will always be connected to the sea quark with opposite sign
tcPPtr particle;
if(lastID*newValence->id() < 0){
particle = parton;
partners.push_back(make_pair(newSea, tPPtr()));
}
else {
particle = newSea;
partners.push_back(make_pair(parton, tPPtr()));
}
anti = newValence->hasAntiColour();
if(rem==theRems.first) theanti.first = anti;
else theanti.second = anti;
if(particle->colourLine())
particle->colourLine()->addAntiColoured(newValence);
if(particle->antiColourLine())
particle->antiColourLine()->addColoured(newValence);
// add the x and return
if(rem==theRems.first) theX.first += currentx;
else theX.second += currentx;
return;
}
void HwRemDecayer::doSplit(pair<tPPtr, tPPtr> partons,
pair<tcPDFPtr, tcPDFPtr> pdfs,
bool first) {
if(theRems.first) {
ParticleVector children=theRems.first->children();
for(unsigned int ix=0;ix<children.size();++ix) {
if(children[ix]->dataPtr()==theRems.first->dataPtr())
theRems.first = dynamic_ptr_cast<RemPPtr>(children[ix]);
}
}
if(theRems.second) {
ParticleVector children=theRems.second->children();
for(unsigned int ix=0;ix<children.size();++ix) {
if(children[ix]->dataPtr()==theRems.second->dataPtr())
theRems.second = dynamic_ptr_cast<RemPPtr>(children[ix]);
}
}
// forced splitting for first parton
if(isPartonic(partons.first )) {
try {
split(partons.first, theContent.first, theRems.first,
theUsed.first, theMaps.first, pdfs.first, first);
}
catch(ShowerHandler::ExtraScatterVeto) {
throw ShowerHandler::ExtraScatterVeto();
}
}
// forced splitting for second parton
if(isPartonic(partons.second)) {
try {
split(partons.second, theContent.second, theRems.second,
theUsed.second, theMaps.second, pdfs.second, first);
// additional check for the remnants
// if can't do the rescale veto the emission
if(!first&&partons.first->data().coloured()&&
partons.second->data().coloured()) {
Lorentz5Momentum pnew[2]=
{theRems.first->momentum() - theUsed.first - partons.first->momentum(),
theRems.second->momentum() - theUsed.second - partons.second->momentum()};
pnew[0].setMass(getParticleData(theContent.first.RemID())->constituentMass());
pnew[0].rescaleEnergy();
pnew[1].setMass(getParticleData(theContent.second.RemID())->constituentMass());
pnew[1].rescaleEnergy();
for(unsigned int iy=0; iy<theRems.first->children().size(); ++iy)
pnew[0] += theRems.first->children()[iy]->momentum();
for(unsigned int iy=0; iy<theRems.second->children().size(); ++iy)
pnew[1] += theRems.second->children()[iy]->momentum();
Lorentz5Momentum ptotal=
theRems.first ->momentum()-partons.first ->momentum()+
theRems.second->momentum()-partons.second->momentum();
// add x limits
if(ptotal.m() < (pnew[0].m() + pnew[1].m()) ) {
if(partons.second->id() != ParticleID::g){
if(partons.second==theMaps.second.back().first)
theUsed.second -= theMaps.second.back().second->momentum();
else
theUsed.second -= theMaps.second.back().first->momentum();
thestep->removeParticle(theMaps.second.back().first);
thestep->removeParticle(theMaps.second.back().second);
}
theMaps.second.pop_back();
theX.second -= partons.second->momentum().rho()/
parent(theRems.second)->momentum().rho();
throw ShowerHandler::ExtraScatterVeto();
}
}
}
catch(ShowerHandler::ExtraScatterVeto){
if(!partons.first||!partons.second||
!theRems.first||!theRems.second)
throw ShowerHandler::ExtraScatterVeto();
//case of the first forcedSplitting worked fine
theX.first -= partons.first->momentum().rho()/
parent(theRems.first)->momentum().rho();
//case of the first interaction
//throw veto immediately, because event get rejected anyway.
if(first) throw ShowerHandler::ExtraScatterVeto();
//secondary interactions have to end on a gluon, if parton
//was NOT a gluon, the forced splitting particles must be removed
if(partons.first->id() != ParticleID::g) {
if(partons.first==theMaps.first.back().first)
theUsed.first -= theMaps.first.back().second->momentum();
else
theUsed.first -= theMaps.first.back().first->momentum();
thestep->removeParticle(theMaps.first.back().first);
thestep->removeParticle(theMaps.first.back().second);
}
theMaps.first.pop_back();
throw ShowerHandler::ExtraScatterVeto();
}
}
// veto if not enough energy for extraction
if( !first &&(theRems.first ->momentum().e() -
partons.first ->momentum().e() < 1.0e-3*MeV ||
theRems.second->momentum().e() -
partons.second->momentum().e() < 1.0e-3*MeV )) {
if(partons.first->id() != ParticleID::g) {
if(partons.first==theMaps.first.back().first)
theUsed.first -= theMaps.first.back().second->momentum();
else
theUsed.first -= theMaps.first.back().first->momentum();
thestep->removeParticle(theMaps.first.back().first);
thestep->removeParticle(theMaps.first.back().second);
}
theMaps.first.pop_back();
if(partons.second->id() != ParticleID::g) {
if(partons.second==theMaps.second.back().first)
theUsed.second -= theMaps.second.back().second->momentum();
else
theUsed.second -= theMaps.second.back().first->momentum();
thestep->removeParticle(theMaps.second.back().first);
thestep->removeParticle(theMaps.second.back().second);
}
theMaps.second.pop_back();
throw ShowerHandler::ExtraScatterVeto();
}
}
void HwRemDecayer::mergeColour(tPPtr pold, tPPtr pnew, bool anti) const {
ColinePtr clnew, clold;
//save the corresponding colour lines
clold = pold->colourLine(anti);
clnew = pnew->colourLine(!anti);
assert(clold);
// There is already a colour line (not the final diquark)
if(clnew){
if( (clnew->coloured().size() + clnew->antiColoured().size()) > 1 ){
if( (clold->coloured().size() + clold->antiColoured().size()) > 1 ){
//join the colour lines
//I don't use the join method, because potentially only (anti)coloured
//particles belong to one colour line
if(clold!=clnew){//procs are not already connected
while ( !clnew->coloured().empty() ) {
tPPtr p = clnew->coloured()[0];
clnew->removeColoured(p);
clold->addColoured(p);
}
while ( !clnew->antiColoured().empty() ) {
tPPtr p = clnew->antiColoured()[0];
clnew->removeAntiColoured(p);
clold->addAntiColoured(p);
}
}
}else{
//if pold is the only member on it's
//colour line, remove it.
clold->removeColoured(pold, anti);
//and add it to clnew
clnew->addColoured(pold, anti);
}
} else{//pnnew is the only member on it's colour line.
clnew->removeColoured(pnew, !anti);
clold->addColoured(pnew, !anti);
}
} else {//there is no coline at all for pnew
clold->addColoured(pnew, !anti);
}
}
void HwRemDecayer::fixColours(PartnerMap partners, bool anti,
double colourDisrupt) const {
PartnerMap::iterator prev;
tPPtr pnew, pold;
assert(partners.size()>=2);
PartnerMap::iterator it=partners.begin();
while(it != partners.end()) {
//skip the first one to have a partner
if(it==partners.begin()){
it++;
continue;
}
prev = it - 1;
//determine the particles to work with
pold = prev->first;
if(prev->second) {
if(!pold->coloured())
pold = prev->second;
else if(pold->hasAntiColour() != anti)
pold = prev->second;
}
assert(pold);
pnew = it->first;
if(it->second) {
if(it->second->colourLine(!anti)) //look for the opposite colour
pnew = it->second;
}
assert(pnew);
// Implement the disruption of colour connections
if( it != partners.end()-1 ) {//last one is diquark-has to be connected
//has to be inside the if statement, so that the probability is
//correctly counted:
if( UseRandom::rnd() < colourDisrupt ){
if(!it->second){//check, whether we have a gluon
mergeColour(pnew, pnew, anti);
}else{
if(pnew==it->first)//be careful about the order
mergeColour(it->second, it->first, anti);
else
mergeColour(it->first, it->second, anti);
}
it = partners.erase(it);
continue;
}
}
// regular merging
mergeColour(pold, pnew, anti);
//end of loop
it++;
}
return;
}
PPtr HwRemDecayer::forceSplit(const tRemPPtr rem, long child, Energy &lastQ,
double &lastx, Lorentz5Momentum &pf,
Lorentz5Momentum &p,
HadronContent & content) const {
static const double eps=1e-6;
// beam momentum
Lorentz5Momentum beam = theBeam->momentum();
// the last scale is minimum of last value and upper limit
Energy minQ=_range*_kinCutoff*sqrt(lastx)/(1-lastx);
if(minQ>lastQ) lastQ=minQ;
// generate the new value of qtilde
// weighted towards the lower value: dP/dQ = 1/Q -> Q(R) =
// Q0 (Qmax/Q0)^R
Energy q;
unsigned int ntry=0,maxtry=100;
double xExtracted = rem==theRems.first ? theX.first : theX.second;
double zmin= lastx/(1.-xExtracted) ,zmax,yy;
if(1-lastx<eps) throw ShowerHandler::ExtraScatterVeto();
do {
q = minQ*pow(lastQ/minQ,UseRandom::rnd());
yy = 1.+0.5*sqr(_kinCutoff/q);
zmax = yy - sqrt(sqr(yy)-1.);
++ntry;
}
while(zmax<zmin&&ntry<maxtry);
if(ntry==maxtry) throw ShowerHandler::ExtraScatterVeto();
if(zmax-zmin<eps) throw ShowerHandler::ExtraScatterVeto();
// now generate z as in FORTRAN HERWIG
// use y = ln(z/(1-z)) as integration variable
double ymin=log(zmin/(1.-zmin));
double ymax=log(zmax/(1.-zmax));
double dely=ymax-ymin;
unsigned int nz=_nbinmax;
dely/=nz;
yy=ymin+0.5*dely;
vector<int> ids;
if(child==21||child==22) {
ids=content.flav;
for(unsigned int ix=0;ix<ids.size();++ix) ids[ix] *= content.sign;
}
else {
ids.push_back(ParticleID::g);
}
// probabilities of the different types of possible splitting
map<long,pair<double,vector<double> > > partonprob;
double ptotal(0.);
for(unsigned int iflav=0;iflav<ids.size();++iflav) {
// only do each parton once
if(partonprob.find(ids[iflav])!=partonprob.end()) continue;
// particle data object
tcPDPtr in = getParticleData(ids[iflav]);
double psum(0.);
vector<double> prob;
for(unsigned int iz=0;iz<nz;++iz) {
double ez=exp(yy);
double wr=1.+ez;
double zr=wr/ez;
double wz=1./wr;
double zz=wz*ez;
double coup = child!=22 ?
_alphaS ->value(sqr(max(wz*q,_kinCutoff))) :
_alphaEM->value(sqr(max(wz*q,_kinCutoff)));
double az=wz*zz*coup;
// g -> q qbar
if(ids[iflav]==ParticleID::g) {
// calculate splitting function
// SP as q is always less than forcedSplitScale, the pdf scale is fixed
// pdfval = _pdf->xfx(theBeamData,in,sqr(q),lastx*zr);
double pdfval=_pdf->xfx(theBeamData,in,sqr(_forcedSplitScale),lastx*zr);
if(pdfval>0.) psum += pdfval*az*0.5*(sqr(zz)+sqr(wz));
}
// q -> q g
else {
// calculate splitting function
// SP as q is always less than forcedSplitScale, the pdf scale is fixed
// pdfval = _pdf->xfx(theBeamData,in,sqr(q),lastx*zr);
double pdfval=_pdf->xfx(theBeamData,in,sqr(_forcedSplitScale),lastx*zr);
if(pdfval>0.) psum += pdfval*az*4./3.*(1.+sqr(wz))*zr;
}
if(psum>0.) prob.push_back(psum);
yy+=dely;
}
if(psum>0.) partonprob[ids[iflav]] = make_pair(psum,prob);
ptotal+=psum;
}
// select the flavour
if(ptotal==0.) throw ShowerHandler::ExtraScatterVeto();
ptotal *= UseRandom::rnd();
map<long,pair<double,vector<double> > >::const_iterator pit;
for(pit=partonprob.begin();pit!=partonprob.end();++pit) {
if(pit->second.first>=ptotal) break;
else ptotal -= pit->second.first;
}
if(pit==partonprob.end())
throw Exception() << "Can't select parton for forced backward evolution in "
<< "HwRemDecayer::forceSplit" << Exception::eventerror;
// select z
unsigned int iz=0;
for(;iz<pit->second.second.size();++iz) {
if(pit->second.second[iz]>ptotal) break;
}
if(iz==pit->second.second.size()) --iz;
double ey=exp(ymin+dely*(float(iz+1)-UseRandom::rnd()));
double z=ey/(1.+ey);
Energy2 pt2=sqr((1.-z)*q)- z*sqr(_kinCutoff);
// create the particle
if(pit->first!=ParticleID::g) child=pit->first;
PPtr parton = getParticleData(child)->produceParticle();
Energy2 emittedm2 = sqr(parton->dataPtr()->constituentMass());
// Now boost pcm and pf to z only frame
Lorentz5Momentum p_ref = Lorentz5Momentum(ZERO, beam.vect());
Lorentz5Momentum n_ref = Lorentz5Momentum(ZERO, -beam.vect());
// generate phi and compute pt of branching
double phi = Constants::twopi*UseRandom::rnd();
Energy pt=sqrt(pt2);
Lorentz5Momentum qt = LorentzMomentum(pt*cos(phi), pt*sin(phi), ZERO, ZERO);
Axis axis(p_ref.vect().unit());
if(axis.perp2()>0.) {
LorentzRotation rot;
double sinth(sqrt(sqr(axis.x())+sqr(axis.y())));
rot.setRotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.));
qt.transform(rot);
}
// compute alpha for previous particle
Energy2 p_dot_n = p_ref*n_ref;
double lastalpha = pf*n_ref/p_dot_n;
Lorentz5Momentum qtout=qt;
Energy2 qtout2=-qt*qt;
double alphaout=(1.-z)/z*lastalpha;
double betaout=0.5*(emittedm2+qtout2)/alphaout/p_dot_n;
Lorentz5Momentum k=alphaout*p_ref+betaout*n_ref+qtout;
k.rescaleMass();
parton->set5Momentum(k);
pf+=k;
lastQ=q;
lastx/=z;
p += parton->momentum();
thestep->addDecayProduct(rem,parton,false);
return parton;
}
void HwRemDecayer::setRemMasses() const {
// get the masses of the remnants
Energy mrem[2];
Lorentz5Momentum ptotal,pnew[2];
vector<tRemPPtr> theprocessed;
theprocessed.push_back(theRems.first);
theprocessed.push_back(theRems.second);
// one remnant in e.g. DIS
if(!theprocessed[0]||!theprocessed[1]) {
tRemPPtr rem = theprocessed[0] ? theprocessed[0] : theprocessed[1];
Lorentz5Momentum deltap(rem->momentum());
// find the diquark and momentum we still need in the energy
tPPtr diquark;
vector<PPtr> progenitors;
for(unsigned int ix=0;ix<rem->children().size();++ix) {
if(!DiquarkMatcher::Check(rem->children()[ix]->data())) {
progenitors.push_back(rem->children()[ix]);
deltap -= rem->children()[ix]->momentum();
}
else
diquark = rem->children()[ix];
}
// now find the total momentum of the hadronic final-state to
// reshuffle against
// find the hadron for this remnant
tPPtr hadron=rem;
do hadron=hadron->parents()[0];
while(!hadron->parents().empty());
// find incoming parton to hard process from this hadron
tPPtr hardin =
generator()->currentEvent()->primaryCollision()->incoming().first==hadron ?
generator()->currentEvent()->primarySubProcess()->incoming().first :
generator()->currentEvent()->primarySubProcess()->incoming().second;
tPPtr parent=hardin;
vector<PPtr> tempprog;
// find the outgoing particles emitted from the backward shower
do {
assert(!parent->parents().empty());
tPPtr newparent=parent->parents()[0];
if(newparent==hadron) break;
for(unsigned int ix=0;ix<newparent->children().size();++ix) {
if(newparent->children()[ix]!=parent)
findChildren(newparent->children()[ix],tempprog);
}
parent=newparent;
}
while(parent!=hadron);
// add to list of potential particles to reshuffle against in right order
for(unsigned int ix=tempprog.size();ix>0;--ix) progenitors.push_back(tempprog[ix-1]);
// final-state particles which are colour connected
tColinePair lines = make_pair(hardin->colourLine(),hardin->antiColourLine());
vector<PPtr> others;
for(ParticleVector::const_iterator
cit = generator()->currentEvent()->primarySubProcess()->outgoing().begin();
cit!= generator()->currentEvent()->primarySubProcess()->outgoing().end();++cit) {
// colour connected
if(lines.first&&lines.first==(**cit).colourLine()) {
findChildren(*cit,progenitors);
continue;
}
// anticolour connected
if(lines.second&&lines.second==(**cit).antiColourLine()) {
findChildren(*cit,progenitors);
continue;
}
// not connected
for(unsigned int ix=0;ix<(**cit).children().size();++ix)
others.push_back((**cit).children()[ix]);
}
// work out how much of the system needed for rescaling
unsigned int iloc=0;
Lorentz5Momentum psystem,ptotal;
do {
psystem+=progenitors[iloc]->momentum();
ptotal = psystem + deltap;
ptotal.rescaleMass();
psystem.rescaleMass();
++iloc;
if(ptotal.mass() > psystem.mass() + diquark->mass() &&
psystem.mass()>1*MeV && DISRemnantOpt_<2 && ptotal.e() > 0.*GeV ) break;
}
while(iloc<progenitors.size());
if(ptotal.mass() > psystem.mass() + diquark->mass()) --iloc;
if(iloc==progenitors.size()) {
// try touching the lepton in dis as a last restort
for(unsigned int ix=0;ix<others.size();++ix) {
progenitors.push_back(others[ix]);
psystem+=progenitors[iloc]->momentum();
ptotal = psystem + deltap;
ptotal.rescaleMass();
psystem.rescaleMass();
++iloc;
}
--iloc;
if(ptotal.mass() > psystem.mass() + diquark->mass()) {
if(DISRemnantOpt_==0||DISRemnantOpt_==2)
Throw<Exception>() << "Warning had to adjust the momentum of the"
<< " non-colour connected"
<< " final-state, e.g. the scattered lepton in DIS"
<< Exception::warning;
else
throw Exception() << "Can't set remnant momentum without adjusting "
<< "the momentum of the"
<< " non-colour connected"
<< " final-state, e.g. the scattered lepton in DIS"
<< " vetoing event"
<< Exception::eventerror;
}
else {
throw Exception() << "Can't put the remnant on-shell in HwRemDecayer::setRemMasses()"
<< Exception::eventerror;
}
}
psystem.rescaleMass();
LorentzRotation R = Utilities::getBoostToCM(make_pair(psystem, deltap));
Energy pz = SimplePhaseSpace::getMagnitude(sqr(ptotal.mass()),
psystem.mass(), diquark->mass());
LorentzRotation Rs(-(R*psystem).boostVector());
Rs.boost(0.0, 0.0, pz/sqrt(sqr(pz) + sqr(psystem.mass())));
Rs = Rs*R;
// put remnant on shell
deltap.transform(R);
deltap.setMass(diquark->mass());
deltap.setE(sqrt(sqr(diquark->mass())+sqr(pz)));
deltap.rescaleRho();
R.invert();
deltap.transform(R);
Rs = R*Rs;
// apply transformation to required particles to absorb recoil
for(unsigned int ix=0;ix<=iloc;++ix) {
progenitors[ix]->deepTransform(Rs);
}
diquark->set5Momentum(deltap);
}
// two remnants
else {
for(unsigned int ix=0;ix<2;++ix) {
if(!theprocessed[ix]) continue;
pnew[ix]=Lorentz5Momentum();
for(unsigned int iy=0;iy<theprocessed[ix]->children().size();++iy) {
pnew[ix]+=theprocessed[ix]->children()[iy]->momentum();
}
mrem[ix]=sqrt(pnew[ix].m2());
}
// now find the remnant remnant cmf frame
Lorentz5Momentum prem[2]={theprocessed[0]->momentum(),
theprocessed[1]->momentum()};
ptotal=prem[0]+prem[1];
ptotal.rescaleMass();
// boost momenta to this frame
if(ptotal.m()< (pnew[0].m()+pnew[1].m()))
throw Exception() << "Not enough energy in both remnants in "
<< "HwRemDecayer::setRemMasses() "
<< Exception::eventerror;
Boost boostv(-ptotal.boostVector());
ptotal.boost(boostv);
for(unsigned int ix=0;ix<2;++ix) {
prem[ix].boost(boostv);
// set the masses and energies,
prem[ix].setMass(mrem[ix]);
prem[ix].setE(0.5/ptotal.m()*(sqr(ptotal.m())+sqr(mrem[ix])-sqr(mrem[1-ix])));
prem[ix].rescaleRho();
// boost back to the lab
prem[ix].boost(-boostv);
// set the momenta of the remnants
theprocessed[ix]->set5Momentum(prem[ix]);
}
// boost the decay products
Lorentz5Momentum ptemp;
for(unsigned int ix=0;ix<2;++ix) {
Boost btorest(-pnew[ix].boostVector());
Boost bfmrest( prem[ix].boostVector());
for(unsigned int iy=0;iy<theprocessed[ix]->children().size();++iy) {
ptemp=theprocessed[ix]->children()[iy]->momentum();
ptemp.boost(btorest);
ptemp.boost(bfmrest);
theprocessed[ix]->children()[iy]->set5Momentum(ptemp);
}
}
}
}
void HwRemDecayer::initSoftInteractions(Energy ptmin, InvEnergy2 beta){
ptmin_ = ptmin;
beta_ = beta;
}
Energy HwRemDecayer::softPt() const {
Energy2 pt2(ZERO);
double xmin(0.0), xmax(1.0), x(0);
if(beta_ == ZERO){
return UseRandom::rnd(0.0,(double)(ptmin_/GeV))*GeV;
}
if(beta_ < ZERO){
xmin = 1.0;
xmax = exp( -beta_*sqr(ptmin_) );
}else{
xmin = exp( -beta_*sqr(ptmin_) );
xmax = 1.0;
}
x = UseRandom::rnd(xmin, xmax);
pt2 = 1.0/beta_ * log(1/x);
if( pt2 < ZERO || pt2 > sqr(ptmin_) )
throw Exception() << "HwRemDecayer::softPt generation of pt "
<< "outside allowed range [0," << ptmin_/GeV << "]."
<< Exception::runerror;
return sqrt(pt2);
}
void HwRemDecayer::softKinematics(Lorentz5Momentum &r1, Lorentz5Momentum &r2,
Lorentz5Momentum &g1, Lorentz5Momentum &g2) const {
g1 = Lorentz5Momentum();
g2 = Lorentz5Momentum();
//All necessary variables for the two soft gluons
Energy pt(softPt()), pz(ZERO);
Energy2 pz2(ZERO);
double phi(UseRandom::rnd(2.*Constants::pi));
double x_g1(0.0), x_g2(0.0);
//Get the external momenta
tcPPair beam(generator()->currentEventHandler()->currentCollision()->incoming());
Lorentz5Momentum P1(beam.first->momentum()), P2(beam.second->momentum());
if(dbg){
cerr << "new event --------------------\n"
<< *(beam.first) << *(softRems_.first)
<< "-------------------\n"
<< *(beam.second) << *(softRems_.second) << endl;
}
//parton mass
Energy mp;
if(quarkPair_){
mp = getParticleData(ParticleID::u)->constituentMass();
}else{
mp = mg_;
}
//Get x_g1 and x_g2
//first limits
double xmin = sqr(ptmin_)/4.0/(P1+P2).m2();
double x1max = (r1.e()+abs(r1.z()))/(P1.e() + abs(P1.z()));
double x2max = (r2.e()+abs(r2.z()))/(P2.e() + abs(P2.z()));
double x1;
if(!multiPeriph_){
//now generate according to 1/x
x_g1 = xmin * exp(UseRandom::rnd(log(x1max/xmin)));
x_g2 = xmin * exp(UseRandom::rnd(log(x2max/xmin)));
}else{
if(valOfN_==0) return;
double param = (1/(2*valOfN_+1))*initTotRap_;
do{
// need 1-x instead of x to get the proper final momenta
x1 = UseRandom::rndGauss(gaussWidth_, 1 - (exp(param)-1)/exp(param));
}while(x1 < 0 || x1>=1.0);
x_g1 = x1max*x1;
x_g2 = x2max*x1;
}
if(dbg)
cerr << x1max << " " << x_g1 << endl << x2max << " " << x_g2 << endl;
Lorentz5Momentum ig1, ig2, cmf;
ig1 = x_g1*P1;
ig2 = x_g2*P2;
ig1.setMass(mp);
ig2.setMass(mp);
ig1.rescaleEnergy();
ig2.rescaleEnergy();
cmf = ig1 + ig2;
//boost vector from cmf to lab
Boost boostv(cmf.boostVector());
//outgoing gluons in cmf
g1.setMass(mp);
g2.setMass(mp);
g1.setX(pt*cos(phi));
g2.setX(-pt*cos(phi));
g1.setY(pt*sin(phi));
g2.setY(-pt*sin(phi));
pz2 = cmf.m2()/4 - sqr(mp) - sqr(pt);
if(pz2/GeV2 < 0.0){
if(dbg)
cerr << "EXCEPTION not enough energy...." << endl;
throw ExtraSoftScatterVeto();
}
if(!multiPeriph_){
if(UseRandom::rndbool())
pz = sqrt(pz2);
else
pz = -sqrt(pz2);
}else{
pz = pz2 > ZERO ? sqrt(pz2) : ZERO;
}
if(dbg)
cerr << "pz has been calculated to: " << pz/GeV << endl;
g1.setZ(pz);
g2.setZ(-pz);
g1.rescaleEnergy();
g2.rescaleEnergy();
if(dbg){
cerr << "check inv mass in cmf frame: " << (g1+g2).m()/GeV
<< " vs. lab frame: " << (ig1+ig2).m()/GeV << endl;
}
g1.boost(boostv);
g2.boost(boostv);
//recalc the remnant momenta
Lorentz5Momentum r1old(r1), r2old(r2);
r1 -= ig1;
r2 -= ig2;
try{
reShuffle(r1, r2, r1old.m(), r2old.m());
}catch(ExtraSoftScatterVeto){
r1 = r1old;
r2 = r2old;
throw ExtraSoftScatterVeto();
}
if(dbg){
cerr << "remnant 1,2 momenta: " << r1/GeV << "--" << r2/GeV << endl;
cerr << "remnant 1,2 masses: " << r1.m()/GeV << " " << r2.m()/GeV << endl;
cerr << "check momenta in the lab..." << (-r1old-r2old+r1+r2+g1+g2)/GeV << endl;
}
}
void HwRemDecayer::doSoftInteractions_old(unsigned int N) {
if(N == 0) return;
if(!softRems_.first || !softRems_.second)
throw Exception() << "HwRemDecayer::doSoftInteractions: no "
<< "Remnants available."
<< Exception::runerror;
if( ptmin_ == -1.*GeV )
throw Exception() << "HwRemDecayer::doSoftInteractions: init "
<< "code has not been called! call initSoftInteractions."
<< Exception::runerror;
Lorentz5Momentum g1, g2;
Lorentz5Momentum r1(softRems_.first->momentum()), r2(softRems_.second->momentum());
unsigned int tries(1), i(0);
for(i=0; i<N; i++){
//check how often this scattering has been regenerated
if(tries > maxtrySoft_) break;
if(dbg){
cerr << "new try \n" << *softRems_.first << *softRems_.second << endl;
}
try{
softKinematics(r1, r2, g1, g2);
}catch(ExtraSoftScatterVeto){
tries++;
i--;
continue;
}
PPair oldrems = softRems_;
PPair gluons = make_pair(addParticle(softRems_.first, ParticleID::g, g1),
addParticle(softRems_.second, ParticleID::g, g2));
//now reset the remnants with the new ones
softRems_.first = addParticle(softRems_.first, softRems_.first->id(), r1);
softRems_.second = addParticle(softRems_.second, softRems_.second->id(), r2);
//do the colour connections
pair<bool, bool> anti = make_pair(oldrems.first->hasAntiColour(),
oldrems.second->hasAntiColour());
ColinePtr cl1 = new_ptr(ColourLine());
ColinePtr cl2 = new_ptr(ColourLine());
if( UseRandom::rnd() < colourDisrupt_ ){//this is the member variable, i.e. SOFT colour disruption
//connect the remnants independent of the gluons
oldrems.first->colourLine(anti.first)->addColoured(softRems_.first, anti.first);
oldrems.second->colourLine(anti.second)->addColoured(softRems_.second, anti.second);
//connect the gluons to each other
cl1->addColoured(gluons.first);
cl1->addAntiColoured(gluons.second);
cl2->addColoured(gluons.second);
cl2->addAntiColoured(gluons.first);
}else{
//connect the remnants to the gluons
oldrems.first->colourLine(anti.first)->addColoured(gluons.first, anti.first);
oldrems.second->colourLine(anti.second)->addColoured(gluons.second, anti.second);
//and the remaining colour line to the final remnant
cl1->addColoured(softRems_.first, anti.first);
cl1->addColoured(gluons.first, !anti.first);
cl2->addColoured(softRems_.second, anti.second);
cl2->addColoured(gluons.second, !anti.second);
}
//reset counter
tries = 1;
}
if(dbg)
cerr << "generated " << i << "th soft scatters\n";
}
void HwRemDecayer::doSoftInteractions_multiPeriph(unsigned int N) {
if(N == 0) return;
int Nmpi = N;
for(int j=0;j<Nmpi;j++){
///////////////////////
// Average number of partons in the ladder
int avgN = floor(ladderMult_*log((softRems_.first->momentum()
+softRems_.second->momentum()).m()/mg_) + ladderbFactor_);
initTotRap_ = abs(softRems_.first->momentum().rapidity())
+abs(softRems_.second->momentum().rapidity());
//generate the poisson distribution with mean avgN
double L = exp(-double(avgN));
int k = 0;
double p = 1;
do {
k++;
p *= UseRandom::rnd();
} while( p > L);
N=k-1;
valOfN_=N;
if(N == 0) return;
if(!softRems_.first || !softRems_.second)
throw Exception() << "HwRemDecayer::doSoftInteractions: no "
<< "Remnants available."
<< Exception::runerror;
if( ptmin_ == -1.*GeV )
throw Exception() << "HwRemDecayer::doSoftInteractions: init "
<< "code has not been called! call initSoftInteractions."
<< Exception::runerror;
Lorentz5Momentum g1, g2, q1, q2;
//intermediate gluons; erased in if there is
//another step
Lorentz5Momentum gint1, gint2;
//momenta of the gluon pair generated in
//each step
vector< pair<Lorentz5Momentum,Lorentz5Momentum> > gluonMomPairs;
//momenta of remnants
Lorentz5Momentum r1(softRems_.first->momentum()), r2(softRems_.second->momentum());
unsigned int tries(1), i(0);
//generate the momenta of particles in the ladder
for(i = 0; i < N; i++){
//check how often this scattering has been regenerated
//and break if exeeded maximal number
- if(tries > maxtrySoft_) break;
-
+ if(tries > maxtrySoft_){
+ if(quarkPair_)return;
+ break;
+ }
if(dbg){
cerr << "new try \n" << *softRems_.first << *softRems_.second << endl;
}
try{
if(i==0){
//generated partons in the ladder are quark-antiquark
quarkPair_ = true;
//first splitting: remnant -> remnant + quark
softKinematics(r1, r2, q1, q2);
}else if(i==1){
//generated pair is gluon-gluon
quarkPair_ = false;
//second splitting; quark -> quark + gluon
softKinematics(q1, q2, g1, g2);
//record generated gluon pair
gluonMomPairs.push_back(make_pair(g1,g2));
}else{
//consequent splittings gluon -> gluon+gluon
//but, the previous gluon gets deleted
quarkPair_ = false;
//save first the previous gluon momentum
g1=gluonMomPairs.back().first;
g2=gluonMomPairs.back().second;
//split gluon momentum
softKinematics(g1, g2, gint1, gint2);
//erase the last entry
gluonMomPairs.pop_back();
//add new gluons
gluonMomPairs.push_back(make_pair(g1,g2));
gluonMomPairs.push_back(make_pair(gint1,gint2));
}
}catch(ExtraSoftScatterVeto){
tries++;
i--;
continue;
}
//reset counter
tries = 1;
}
//if no gluons discard the ladder
- if(gluonMomPairs.size()==0) return;
+ //if(gluonMomPairs.size()==0) return;
if(dbg)
cerr << "generated " << i << "th soft scatters\n";
//gluons from previous step
PPair oldgluons;
//quark-antiquark pair
PPair quarks;
//colour direction of quark-antiquark (true if anticolour)
pair<bool, bool> anti_q;
//generate particles (quark-antiquark and gluons) in the
//ladder from momenta generated above
for(i = 0; i <= gluonMomPairs.size(); i++){
//remnants before splitting
PPair oldrems = softRems_;
//current gluon pair
PPair gluons;
//add quarks
if(i==0){
quarks = make_pair(addParticle(softRems_.first, ParticleID::u, q1),
addParticle(softRems_.second, ParticleID::ubar, q2));
anti_q = make_pair(quarks.first->hasAntiColour(),
quarks.second->hasAntiColour());
}
if(i>0){
gluons = make_pair(addParticle(softRems_.first, ParticleID::g,
gluonMomPairs[i-1].first),
addParticle(softRems_.second, ParticleID::g,
gluonMomPairs[i-1].second));
}
//now reset the remnants with the new ones
softRems_.first = addParticle(softRems_.first, softRems_.first->id(), r1);
softRems_.second = addParticle(softRems_.second, softRems_.second->id(), r2);
//colour direction of remnats
pair<bool, bool> anti = make_pair(oldrems.first->hasAntiColour(),
oldrems.second->hasAntiColour());
ColinePtr cl1 = new_ptr(ColourLine());
ColinePtr cl2 = new_ptr(ColourLine());
//first colour connect remnants and if no
//gluons the quark-antiquark pair
if(i==0){
oldrems.first->colourLine(anti.first)->addColoured(softRems_.first,
anti.first);
oldrems.second->colourLine(anti.second)->addColoured(softRems_.second,
anti.second);
if(gluonMomPairs.size()==0){
ColinePtr clq = new_ptr(ColourLine());
clq->addColoured(quarks.first, anti_q.first);
clq->addColoured(quarks.second, anti_q.second);
}
}//colour connect remnants from previous step and gluons to quarks or gluons
if(i!=0){
oldrems.first->colourLine(anti.first)->addColoured(softRems_.first,
anti.first);
oldrems.second->colourLine(anti.second)->addColoured(softRems_.second,
anti.second);
if(i==1){
cl1->addColoured(quarks.first, anti_q.first);
cl1->addColoured(gluons.first, !anti_q.first);
cl2->addColoured(quarks.second, anti_q.second);
cl2->addColoured(gluons.second, !anti_q.second);
}else{
cl1->addColoured(oldgluons.first, anti_q.first);
cl1->addColoured(gluons.first, !anti_q.first);
cl2->addColoured(oldgluons.second, anti_q.second);
cl2->addColoured(gluons.second, !anti_q.second);
}
} //last step; connect last gluons
if(i > 0 && i == gluonMomPairs.size()){
ColinePtr clg = new_ptr(ColourLine());
clg->addColoured(gluons.first, anti_q.first);
clg->addColoured(gluons.second, anti_q.second);
}
//save gluons for the next step
if(i < gluonMomPairs.size()) oldgluons = gluons;
}
////////
}
////////
}
void HwRemDecayer::finalize(double colourDisrupt, unsigned int softInt){
PPair diquarks;
//Do the final Rem->Diquark or Rem->quark "decay"
if(theRems.first) {
diquarks.first = finalSplit(theRems.first, theContent.first.RemID(),
theUsed.first);
theMaps.first.push_back(make_pair(diquarks.first, tPPtr()));
}
if(theRems.second) {
diquarks.second = finalSplit(theRems.second, theContent.second.RemID(),
theUsed.second);
theMaps.second.push_back(make_pair(diquarks.second, tPPtr()));
}
setRemMasses();
if(theRems.first) {
fixColours(theMaps.first, theanti.first, colourDisrupt);
if(theContent.first.hadron->id()==ParticleID::pomeron&&
pomeronStructure_==0) fixColours(theMaps.first, !theanti.first, colourDisrupt);
}
if(theRems.second) {
fixColours(theMaps.second, theanti.second, colourDisrupt);
if(theContent.second.hadron->id()==ParticleID::pomeron&&
pomeronStructure_==0) fixColours(theMaps.second, !theanti.second, colourDisrupt);
}
if( !theRems.first || !theRems.second ) return;
//stop here if we don't have two remnants
softRems_ = diquarks;
doSoftInteractions(softInt);
}
HwRemDecayer::HadronContent
HwRemDecayer::getHadronContent(tcPPtr hadron) const {
HadronContent hc;
hc.hadron = hadron->dataPtr();
long id(hadron->id());
// baryon
if(BaryonMatcher::Check(hadron->data())) {
hc.sign = id < 0? -1: 1;
hc.flav.push_back((id = abs(id)/10)%10);
hc.flav.push_back((id /= 10)%10);
hc.flav.push_back((id /= 10)%10);
hc.extracted = -1;
}
else if(hadron->data().id()==ParticleID::gamma ||
(hadron->data().id()==ParticleID::pomeron && pomeronStructure_==1)) {
hc.sign = 1;
for(int ix=1;ix<6;++ix) {
hc.flav.push_back( ix);
hc.flav.push_back(-ix);
}
}
else if(hadron->data().id()==ParticleID::pomeron ) {
hc.sign = 1;
hc.flav.push_back(ParticleID::g);
hc.flav.push_back(ParticleID::g);
}
else if(hadron->data().id()==ParticleID::reggeon ) {
hc.sign = 1;
for(int ix=1;ix<3;++ix) {
hc.flav.push_back( ix);
hc.flav.push_back(-ix);
}
}
hc.pomeronStructure = pomeronStructure_;
return hc;
}
long HwRemDecayer::HadronContent::RemID() const{
if(extracted == -1)
throw Exception() << "Try to build a Diquark id without "
<< "having extracted something in "
<< "HwRemDecayer::RemID(...)"
<< Exception::runerror;
//the hadron was a meson or photon
if(flav.size()==2) return sign*flav[(extracted+1)%2];
long remId;
int id1(sign*flav[(extracted+1)%3]),
id2(sign*flav[(extracted+2)%3]),
sign(0), spin(0);
if (abs(id1) > abs(id2)) swap(id1, id2);
sign = (id1 < 0) ? -1 : 1; // Needed for the spin 0/1 part
remId = id2*1000+id1*100;
// Now decide if we have spin 0 diquark or spin 1 diquark
if(id1 == id2) spin = 3; // spin 1
else spin = 1; // otherwise spin 0
remId += sign*spin;
return remId;
}
tPPtr HwRemDecayer::addParticle(tcPPtr parent, long id, Lorentz5Momentum p) const {
PPtr newp = new_ptr(Particle(getParticleData(id)));
newp->set5Momentum(p);
// Add the new remnant to the step, but don't do colour connections
thestep->addDecayProduct(parent,newp,false);
return newp;
}
void HwRemDecayer::findChildren(tPPtr part,vector<PPtr> & particles) const {
if(part->children().empty()) particles.push_back(part);
else {
for(unsigned int ix=0;ix<part->children().size();++ix)
findChildren(part->children()[ix],particles);
}
}
ParticleVector HwRemDecayer::decay(const DecayMode &,
const Particle &, Step &) const {
throw Exception() << "HwRemDecayer::decay(...) "
<< "must not be called explicitely."
<< Exception::runerror;
}
void HwRemDecayer::persistentOutput(PersistentOStream & os) const {
os << ounit(_kinCutoff, GeV) << _range << _zbin << _ybin
<< _nbinmax << _alphaS << _alphaEM << DISRemnantOpt_
<< maxtrySoft_ << colourDisrupt_ << ladderMult_ << ladderbFactor_ << pomeronStructure_
<< ounit(mg_,GeV) << ounit(ptmin_,GeV) << ounit(beta_,sqr(InvGeV))
<< allowTop_ << multiPeriph_ << valOfN_ << initTotRap_;
}
void HwRemDecayer::persistentInput(PersistentIStream & is, int) {
is >> iunit(_kinCutoff, GeV) >> _range >> _zbin >> _ybin
>> _nbinmax >> _alphaS >> _alphaEM >> DISRemnantOpt_
>> maxtrySoft_ >> colourDisrupt_ >> ladderMult_ >> ladderbFactor_ >> pomeronStructure_
>> iunit(mg_,GeV) >> iunit(ptmin_,GeV) >> iunit(beta_,sqr(InvGeV))
>> allowTop_ >> multiPeriph_ >> valOfN_ >> initTotRap_;
}
ClassDescription<HwRemDecayer> HwRemDecayer::initHwRemDecayer;
// Definition of the static class description member.
void HwRemDecayer::Init() {
static ClassDocumentation<HwRemDecayer> documentation
("The HwRemDecayer class decays the remnant for Herwig");
static Parameter<HwRemDecayer,double> interfaceZBinSize
("ZBinSize",
"The size of the vbins in z for the interpolation of the splitting function.",
&HwRemDecayer::_zbin, 0.05, 0.001, 0.1,
false, false, Interface::limited);
static Parameter<HwRemDecayer,int> interfaceMaxBin
("MaxBin",
"Maximum number of z bins",
&HwRemDecayer::_nbinmax, 100, 10, 1000,
false, false, Interface::limited);
static Reference<HwRemDecayer,ShowerAlpha> interfaceAlphaS
("AlphaS",
"Pointer to object to calculate the strong coupling",
&HwRemDecayer::_alphaS, false, false, true, false, false);
static Reference<HwRemDecayer,ShowerAlpha> interfaceAlphaEM
("AlphaEM",
"Pointer to object to calculate the electromagnetic coupling",
&HwRemDecayer::_alphaEM, false, false, true, false, false);
static Parameter<HwRemDecayer,Energy> interfaceKinCutoff
("KinCutoff",
"Parameter kinCutoff used to constrain qtilde",
&HwRemDecayer::_kinCutoff, GeV, 0.75*GeV, 0.5*GeV, 10.0*GeV,
false, false, Interface::limited);
static Parameter<HwRemDecayer,double> interfaceEmissionRange
("EmissionRange",
"Factor above the minimum possible value in which the forced splitting is allowed.",
&HwRemDecayer::_range, 1.1, 1.0, 10.0,
false, false, Interface::limited);
static Switch<HwRemDecayer,unsigned int> interfaceDISRemnantOption
("DISRemnantOption",
"Options for the treatment of the remnant in DIS",
&HwRemDecayer::DISRemnantOpt_, 0, false, false);
static SwitchOption interfaceDISRemnantOptionDefault
(interfaceDISRemnantOption,
"Default",
"Use the minimum number of particles needed to take the recoil"
" and allow the lepton to be used if needed",
0);
static SwitchOption interfaceDISRemnantOptionNoLepton
(interfaceDISRemnantOption,
"NoLepton",
"Use the minimum number of particles needed to take the recoil but"
" veto events where the lepton kinematics would need to be altered",
1);
static SwitchOption interfaceDISRemnantOptionAllParticles
(interfaceDISRemnantOption,
"AllParticles",
"Use all particles in the colour connected system to take the recoil"
" and use the lepton if needed.",
2);
static SwitchOption interfaceDISRemnantOptionAllParticlesNoLepton
(interfaceDISRemnantOption,
"AllParticlesNoLepton",
"Use all the particles in the colour connected system to take the"
" recoil but don't use the lepton.",
3);
static Parameter<HwRemDecayer,unsigned int> interfaceMaxTrySoft
("MaxTrySoft",
"The maximum number of regeneration attempts for an additional soft scattering",
&HwRemDecayer::maxtrySoft_, 10, 0, 100,
false, false, Interface::limited);
static Parameter<HwRemDecayer,double> interfacecolourDisrupt
("colourDisrupt",
"Fraction of connections to additional soft subprocesses, which are colour disrupted.",
&HwRemDecayer::colourDisrupt_,
1.0, 0.0, 1.0,
false, false, Interface::limited);
static Parameter<HwRemDecayer,double> interfaceladderMult
("ladderMult",
"The multiplicity factor in the multiperipheral ladder.",
&HwRemDecayer::ladderMult_,
1.0, 0.0, 10.0,
false, false, Interface::limited);
static Parameter<HwRemDecayer,double> interfaceladderbFactor
("ladderbFactor",
"The additive factor in the multiperipheral ladder multiplicity.",
&HwRemDecayer::ladderbFactor_,
1.0, 0.0, 10.0,
false, false, Interface::limited);
static Parameter<HwRemDecayer,double> interfacegaussWidth
("gaussWidth",
"The gaussian width of the fluctuation of longitudinal momentum fraction.",
&HwRemDecayer::gaussWidth_,
0.1, 0.0, 1.0,
false, false, Interface::limited);
static Switch<HwRemDecayer,unsigned int> interfacePomeronStructure
("PomeronStructure",
"Option for the treatment of the valance structure of the pomeron",
&HwRemDecayer::pomeronStructure_, 0, false, false);
static SwitchOption interfacePomeronStructureGluon
(interfacePomeronStructure,
"Gluon",
"Assume the pomeron is a two gluon state",
0);
static SwitchOption interfacePomeronStructureQQBar
(interfacePomeronStructure,
"QQBar",
"Assumne the pomeron is q qbar as for the photon,"
" this option is not recommended and is provide for compatiblity with POMWIG",
1);
static Switch<HwRemDecayer,bool> interfaceAllowTop
("AllowTop",
"Allow top quarks in the hadron",
&HwRemDecayer::allowTop_, false, false, false);
static SwitchOption interfaceAllowTopNo
(interfaceAllowTop,
"No",
"Don't allow them",
false);
static SwitchOption interfaceAllowTopYes
(interfaceAllowTop,
"Yes",
"Allow them",
true);
static Switch<HwRemDecayer,bool> interfaceMultiPeriph
("MultiPeriph",
"Use multiperipheral kinematics",
&HwRemDecayer::multiPeriph_, false, false, false);
static SwitchOption interfaceMultiPeriphNo
(interfaceMultiPeriph,
"No",
"Don't use multiperipheral",
false);
static SwitchOption interfaceMultiPeriphYes
(interfaceMultiPeriph,
"Yes",
"Use multiperipheral kinematics",
true);
}
bool HwRemDecayer::canHandle(tcPDPtr particle, tcPDPtr parton) const {
if(! (StandardQCDPartonMatcher::Check(*parton) || parton->id()==ParticleID::gamma) ) {
if(abs(parton->id())==ParticleID::t) {
if(!allowTop_)
throw Exception() << "Top is not allow as a parton in the remant handling, please "
<< "use a PDF which does not contain top for the remnant"
<< " handling (preferred) or allow top in the remnant using\n"
<< " set " << fullName() << ":AllowTop Yes\n"
<< Exception::runerror;
}
else
return false;
}
return HadronMatcher::Check(*particle) || particle->id()==ParticleID::gamma
|| particle->id()==ParticleID::pomeron || particle->id()==ParticleID::reggeon;
}
bool HwRemDecayer::isPartonic(tPPtr parton) const {
if(parton->parents().empty()) return false;
tPPtr parent = parton->parents()[0];
bool partonic = false;
for(unsigned int ix=0;ix<parent->children().size();++ix) {
if(dynamic_ptr_cast<tRemPPtr>(parent->children()[ix])) {
partonic = true;
break;
}
}
return partonic;
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 6:47 PM (5 h, 42 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023868
Default Alt Text
(102 KB)

Event Timeline