Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F7879447
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
88 KB
Subscribers
None
View Options
diff --git a/PDF/HwRemDecayer.cc b/PDF/HwRemDecayer.cc
--- a/PDF/HwRemDecayer.cc
+++ b/PDF/HwRemDecayer.cc
@@ -1,2021 +1,2018 @@
// -*- C++ -*-
//
// HwRemDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the HwRemDecayer class.
//
#include "HwRemDecayer.h"
#include "ThePEG/Utilities/DescribeClass.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->dataPtr()->coloured()) {
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( !parton->coloured()) {
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;
bool qed = child==22;
if(child==21||child==22) {
ids=content.flav;
for(unsigned int ix=0;ix<ids.size();++ix) ids[ix] *= content.sign;
}
else if(abs(child)<=6) {
ids.push_back(ParticleID::g);
}
else {
ids.push_back(ParticleID::gamma);
qed=true;
}
// 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 = !qed ?
_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 || ids[iflav]==ParticleID::gamma) {
// 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 &&
pit->first!=ParticleID::gamma) 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 {
+void HwRemDecayer::setRemMasses(PPair diquarks) 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;
+ tPPtr diquark = theprocessed[0] ? diquarks.first : diquarks.second;
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];
+ if(rem->children()[ix]==diquark) continue;
+ progenitors.push_back(rem->children()[ix]);
+ deltap -= rem->children()[ix]->momentum();
}
// 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;
//ofstream myfile2("softPt.txt", ios::app );
//myfile2 << pt2/GeV2 <<" "<<sqrt(pt2)/GeV<< endl;
//myfile2.close();
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) - (pt*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 << "pz1 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 -= g1;
r2 -= g2;
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());
// case 2:
oldrems.first->colourLine(anti.first)
->addColoured(gluons.second,anti.second);
cl2->addColoured(softRems_.first, anti.second);
cl2->addColoured(gluons.second, !anti.second);
oldrems.first->colourLine(anti.first)
->addColoured(gluons.second,anti.second);
oldrems.second->colourLine(anti.second)
->addColoured(gluons.first,anti.first);
cl1->addColoured(softRems_.second, anti.first);
cl1->addColoured(gluons.first, !anti.first);
cl2->addColoured(softRems_.first, anti.second);
cl2->addColoured(gluons.second, !anti.second);
//reset counter
tries = 1;
}
if(dbg)
cerr << "generated " << i << "th soft scatters\n";
}
// Solve the reshuffling equation to rescale the remnant momenta
double bisectReshuffling(const vector<PPtr>& particles,
Energy w,
double target = -16., double maxLevel = 80.) {
double level = 0;
double left = 0;
double right = 1;
double check = -1.;
double xi = -1;
while ( level < maxLevel ) {
xi = (left+right)*pow(0.5,level+1.);
check = 0.;
for (vector<PPtr>::const_iterator p = particles.begin(); p != particles.end(); ++p){
check += sqrt(sqr(xi)*((*p)->momentum().vect().mag2())+sqr((*p)->mass()))/w;
}
if ( check==1. || log10(abs(1.-check)) <= target )
break;
left *= 2.;
right *= 2.;
if ( check >= 1. ) {
right -= 1.;
++level;
}
if ( check < 1. ) {
left += 1.;
++level;
}
}
return xi;
}
LorentzRotation HwRemDecayer::rotate(const LorentzMomentum &p) const {
LorentzRotation R;
static const double ptcut = 1e-20;
Energy2 pt2 = sqr(p.x())+sqr(p.y());
Energy2 pp2 = sqr(p.z())+pt2;
double phi, theta;
if(pt2 <= pp2*ptcut) {
if(p.z() > ZERO) theta = 0.;
else theta = Constants::pi;
phi = 0.;
} else {
Energy pp = sqrt(pp2);
Energy pt = sqrt(pt2);
double ct = p.z()/pp;
double cf = p.x()/pt;
phi = -acos(cf);
theta = acos(ct);
}
// Rotate first around the z axis to put p in the x-z plane
// Then rotate around the Y axis to put p on the z axis
R.rotateZ(phi).rotateY(theta);
return R;
}
struct vectorSort{
bool operator() (Lorentz5Momentum i,Lorentz5Momentum j) {return(i.rapidity() < j.rapidity());}
} ySort;
void HwRemDecayer::doSoftInteractions_multiPeriph(unsigned int N) {
if(N == 0) return;
int Nmpi = N;
for(int j=0;j<Nmpi;j++){
///////////////////////
// TODO: parametrization of the ladder multiplicity (need to tune to 900GeV, 7Tev and 13Tev)
// Parameterize the ladder multiplicity to: ladderMult_ = A_0 * (s/1TeV^2)^alpha
// with the two tunable parameters A_0 =ladderNorm_ and alpha = ladderPower_
// Get the collision energy
//Energy energy(generator()->maximumCMEnergy());
//double reference = sqr(energy/TeV);
// double ladderMult_;
// Parametrization of the ladder multiplicity
// ladderMult_ = ladderNorm_ * pow( ( reference ) , ladderPower_ );
double avgN = 2.*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
N=UseRandom::rndPoisson(avgN);
valOfN_=N;
if(N <= 1){
// j--; //TODO: Do we want to make all Nmpi soft MPIs?
// Compare to MaxTryMPI for hard mpis.
continue;
}
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;
// The remnants
PPtr rem1 = softRems_.first;
PPtr rem2 = softRems_.second;
// Vector for the ladder particles
vector<Lorentz5Momentum> ladderMomenta;
// Remnant momenta
Lorentz5Momentum r1(softRems_.first->momentum()), r2(softRems_.second->momentum());
Lorentz5Momentum cm =r1+r2;
// Initialize partons in the ladder
// The toy masses are needed for the correct calculation of the available energy
Lorentz5Momentum sumMomenta;
for(unsigned int i = 0; i < N; i++) {
// choose constituents
Energy newMass = ZERO;
Energy toyMass;
if(i<2){
// u and d have the same mass so its enough to use u
toyMass = getParticleData(ParticleID::u)->constituentMass();
}
else{
toyMass = getParticleData(ParticleID::g)->constituentMass();
}
Lorentz5Momentum cp(ZERO,ZERO,ZERO,newMass,newMass);
// dummy container for the momentum that is used for momentum conservation
Lorentz5Momentum dummy(ZERO,ZERO,ZERO,toyMass,toyMass);
ladderMomenta.push_back(cp);
sumMomenta+=dummy;
}
// Get the beam energy
tcPPair beam(generator()->currentEventHandler()->currentCollision()->incoming());
//Lorentz5Momentum P1(beam.first->momentum()), P2(beam.second->momentum());
// Calculate available energy for the partons
double x1;//,x2;
double param = (1./(valOfN_+1.))*initTotRap_;
do{
// Need 1-x instead of x to get the proper final momenta
// TODO: physical to use different x's (see comment below)
x1 = UseRandom::rndGauss( gaussWidth_ , exp(-param) );
// x2 = UseRandom::rndGauss( gaussWidth_ , exp(-param) );
}while(x1 < 0 || x1>=1.0); // x2 < 0 || x2>=1.0);
// Remnants 1 and 2 need to be rescaled later by this amount
Lorentz5Momentum ig1 = x1*r1;
Lorentz5Momentum ig2 = x1*r2; //TODO: x2*r2
// requires boost of Ladder in x1/x2-dependent
// frame.
// If the remaining remnant energy is not sufficient for the restmass of the remnants
// then continue/try again
if ( cm.m() - (ig1+ig2).m() < r1.m()+r2.m() ){
continue;
}
// The available energy that is used to generate the ladder
// sumMomenta is the the sum of rest masses of the ladder partons
// the available energy goes all into the kinematics
Energy availableEnergy = (ig1+ig2).m() - sumMomenta.m();
// If not enough energy then continue
// The available energy has to be larger then the rest mass of the remnants
if ( availableEnergy < ZERO ) {
// j--; //TODO: Do we want to make all Nmpi soft MPIs?
continue;
}
unsigned int its(0);
// Generate the momenta of the partons in the ladder
if ( !(doPhaseSpaceGenerationGluons(ladderMomenta,availableEnergy,its)) ){
// j--; //TODO: Do we want to make all Nmpi soft MPIs?
continue;
}
// Add gluon mass and rescale
Lorentz5Momentum totalMomPartons;
Lorentz5Momentum totalMassLessPartons;
// Sort the ladder partons according to their rapidity and then choose which ones will be the quarks
sort(ladderMomenta.begin(),ladderMomenta.end(),ySort);
int countPartons=0;
long quarkID=0;
// Choose between up and down quarks
int choice = UseRandom::rnd2(1,1);
switch (choice) {
case 0: quarkID = ParticleID::u; break;
case 1: quarkID = ParticleID::d; break;
}
for (auto &p:ladderMomenta){
totalMomPartons+=p;
// Set the mass of the gluons and the two quarks in the ladder
if(countPartons==0 || countPartons==int(ladderMomenta.size()-1)){
p.setMass( getParticleData(quarkID)->constituentMass() );
}else{
p.setMass( getParticleData(ParticleID::g)->constituentMass() );
}
p.rescaleEnergy();
countPartons++;
}
// Continue if energy conservation is violated
if ( abs(availableEnergy - totalMomPartons.m()) > 1e-8*GeV){
// j--; //TODO: Do we want to make all Nmpi soft MPIs?
continue;
}
// Boost momenta into CM frame
const Boost boostv(-totalMomPartons.boostVector());
Lorentz5Momentum totalMomentumAfterBoost;
for ( unsigned int i=0; i<ladderMomenta.size(); i++){
ladderMomenta[i].boost(boostv);
totalMomentumAfterBoost +=ladderMomenta[i];
}
const Boost boostvR(-cm.boostVector());
r1.boost(boostvR);
r2.boost(boostvR);
// Remaining energy after generation of soft ladder
Energy remainingEnergy = cm.m() - totalMomentumAfterBoost.m();
// Continue if not enough energy
if(remainingEnergy<ZERO) {
// j--; //TODO: Do we want to make all Nmpi soft MPIs?
continue;
}
vector<PPtr> remnants;
rem1->set5Momentum(r1);
rem2->set5Momentum(r2);
remnants.push_back(rem1);
remnants.push_back(rem2);
vector<PPtr> reshuffledRemnants;
Lorentz5Momentum totalMomentumAll;
// Bisect reshuffling for rescaling of remnants
double xi_remnants = bisectReshuffling(remnants,remainingEnergy);
// Rescale remnants
for ( auto &rems: remnants ) {
Lorentz5Momentum reshuffledMomentum;
reshuffledMomentum = xi_remnants*rems->momentum();
reshuffledMomentum.setMass(getParticleData(softRems_.first->id())->constituentMass());
reshuffledMomentum.rescaleEnergy();
reshuffledMomentum.boost(-boostvR);
rems->set5Momentum(reshuffledMomentum);
totalMomentumAll+=reshuffledMomentum;
}
// Then the other particles
for ( auto &p:ladderMomenta ) {
p.boost(-boostvR);
totalMomentumAll+=p;
}
// sanity check
if ( abs(cm.m() - totalMomentumAll.m()) > 1e-8*GeV) {
continue;
}
// sort again
sort(ladderMomenta.begin(),ladderMomenta.end(),ySort);
// Do the colour connections
// Original rems are the ones which are connected to other parts of the event
PPair oldRems_ = softRems_;
pair<bool, bool> anti = make_pair(oldRems_.first->hasAntiColour(),
oldRems_.second->hasAntiColour());
// Replace first remnant
softRems_.first = addParticle(softRems_.first, softRems_.first->id(),
remnants[0]->momentum());
// Connect the old remnant to the new remnant
oldRems_.first->colourLine(anti.first)->addColoured(softRems_.first, anti.first);
// Replace second remnant
softRems_.second = addParticle(softRems_.second, softRems_.second->id(),
remnants[1]->momentum());
// This connects the old remnants to the new remnants
oldRems_.second->colourLine(anti.second)->addColoured(softRems_.second, anti.second);
// Add all partons to the first remnant for the event record
vector<PPtr> partons;
vector<PPtr> quarks;
int count=0;
// Choose the colour connections and position of quark antiquark
// Choose between R1-q-g..g-qbar-R2 or R1-qbar-g...g-q-R2
// (place of quark antiquarks in the ladder)
int quarkPosition = UseRandom::rnd2(1,1);
for (auto &p:ladderMomenta){
if(p.mass()==getParticleData(ParticleID::u)->constituentMass()){
if(count==0){
if(quarkPosition==0){
quarks.push_back(addParticle(softRems_.first, quarkID, p));
count++;
}else{
quarks.push_back(addParticle(softRems_.first, -quarkID, p));
count++;
}
}else{
if(quarkPosition==0){
quarks.push_back(addParticle(softRems_.first, -quarkID, p));
}else{
quarks.push_back(addParticle(softRems_.first, quarkID, p));
}
}
}else{
partons.push_back(addParticle(softRems_.first, ParticleID::g, p));
}
softRems_.first = addParticle(softRems_.first, softRems_.first->id(),
softRems_.first->momentum());
oldRems_.first->colourLine(anti.first)->addColoured(softRems_.first, anti.first);
}
// Need to differenciate between the two quark positions, this defines the
// colour connections to the new remnants and old remnants
if(quarkPosition==0){
// ladder self contained
if(partons.size()==0 && quarks.size()>0){
ColinePtr clq = new_ptr(ColourLine());
clq->addColoured(quarks[0]);
clq->addAntiColoured(quarks[1]);
}
ColinePtr clfirst = new_ptr(ColourLine());
ColinePtr cllast = new_ptr(ColourLine());
if(partons.size()>0){
clfirst->addColoured(quarks[0]);
clfirst->addAntiColoured(partons[0]);
cllast->addAntiColoured(quarks[1]);
cllast->addColoured(partons[partons.size()-1]);
//now the remaining gluons
for (unsigned int i=0; i<partons.size()-1; i++){
ColinePtr cl = new_ptr(ColourLine());
cl->addColoured(partons[i]);
cl->addAntiColoured(partons[i+1]);
}
}
} else {
if(partons.size()==0 && quarks.size()>0){
ColinePtr clq = new_ptr(ColourLine());
clq->addAntiColoured(quarks[0]);
clq->addColoured(quarks[1]);
}
ColinePtr clfirst = new_ptr(ColourLine());
ColinePtr cllast = new_ptr(ColourLine());
if(partons.size()>0){
clfirst->addAntiColoured(quarks[0]);
clfirst->addColoured(partons[0]);
cllast->addColoured(quarks[1]);
cllast->addAntiColoured(partons[partons.size()-1]);
//now the remaining gluons
for (unsigned int i=0; i<partons.size()-1; i++){
ColinePtr cl = new_ptr(ColourLine());
cl->addAntiColoured(partons[i]);
cl->addColoured(partons[i+1]);
}
}
}// end colour connection loop
}// end Nmpi loop
}//end function
// Do the phase space generation here is 1 to 1 the same from UA5 model
bool HwRemDecayer::doPhaseSpaceGenerationGluons(vector<Lorentz5Momentum> &softGluons, Energy CME, unsigned int &its)
const{
// Define the parameters
unsigned int _maxtries = 300;
double alog = log(CME*CME/GeV2);
unsigned int ncl = softGluons.size();
// calculate the slope parameters for the different clusters
// outside loop to save time
vector<Lorentz5Momentum> mom(ncl);
// Sets the slopes depending on the constituent quarks of the cluster
for(unsigned int ix=0;ix<ncl;++ix)
{
mom[ix]=softGluons[ix];
}
// generate the momenta
double eps = 1e-10/double(ncl);
vector<double> xi(ncl);
vector<Energy> tempEnergy(ncl);
Energy sum1(ZERO);
double yy(0.);
// We want to make sure that the first Pt is from the
// desired pt-distribution. If we select the first pt in the
// trial loop we introduce a bias.
Energy firstPt=softPt();
while(its < _maxtries) {
++its;
Energy sumx = ZERO;
Energy sumy = ZERO;
unsigned int iterations(0);
unsigned int _maxtriesNew = 100;
while(iterations < _maxtriesNew) {
iterations++;
Energy sumxIt = ZERO;
Energy sumyIt = ZERO;
bool success=false;
Energy pTmax=ZERO;
for(unsigned int i = 0; i<ncl; ++i) {
// Different options for soft pt sampling
//1) pT1>pT2...pTN
//2) pT1>pT2>..>pTN
//3) flat
//4) y dependent
//5) Frist then flat
int triesPt=0;
Energy pt;
//Energy ptTest;
switch(PtDistribution_) {
case 0: //default softPt()
pt=softPt();
break;
case 1: //pTordered
if(i==0){
pt=softPt();
pTmax=pt;
}else{
do{
pt=softPt();
}while(pt>pTmax);
}
break;
case 2: //strongly pT ordered
if ( i==0 ) {
pt=softPt();
pTmax=pt;
} else {
do {
if ( triesPt==20 ) {
pt=pTmax;
break;
}
pt=softPt();
triesPt++;
} while ( pt>pTmax );
pTmax=pt;
}
break;
case 3: //flat
pt = UseRandom::rnd(0.0,(double)(ptmin_/GeV))*GeV;
break;
case 4: //flat below first pT
if ( i==0 ) {
pt = firstPt;
} else {
pt = firstPt * UseRandom::rnd();
}
break;
case 5: //flat but rising below first pT
if ( i==0 ) {
pt=firstPt;
} else {
pt = firstPt * pow(UseRandom::rnd(),1/2);
}
}
Energy2 ptp = pt*pt;
if(ptp <= ZERO) pt = - sqrt(-ptp);
else pt = sqrt(ptp);
// randomize azimuth
Energy px,py;
//randomize the azimuth, but the last one should cancel all others
if(i<ncl-1){
randAzm(pt,px,py);
// set transverse momentum
mom[i].setX(px);
mom[i].setY(py);
sumxIt += px;
sumyIt += py;
}else{
//calculate azimuth angle s.t
// double factor;
Energy pTdummy;
pTdummy = sqrt(sumxIt*sumxIt+sumyIt*sumyIt);
if( pTdummy < ptmin_ ){
px=-sumxIt;
py=-sumyIt;
mom[i].setX(px);
mom[i].setY(py);
sumxIt+=px;
sumyIt+=py;
sumx = sumxIt;
sumy = sumyIt;
success=true;
}
}
}
if(success){
break;
}
}
sumx /= ncl;
sumy /= ncl;
// find the sum of the transverse mass
Energy sumtm=ZERO;
for(unsigned int ix = 0; ix<ncl; ++ix) {
mom[ix].setX(mom[ix].x()-sumx);
mom[ix].setY(mom[ix].y()-sumy);
Energy2 pt2 = mom[ix].perp2();
// Use the z component of the clusters momentum for temporary storage
mom[ix].setZ(sqrt(pt2+mom[ix].mass2()));
sumtm += mom[ix].z();
}
// if transverse mass greater the CMS try again
if(sumtm > CME) continue;
// randomize the mom vector to get the first and the compensating parton
// at all possible positions:
long (*p_irnd)(long) = UseRandom::irnd;
random_shuffle(mom.begin(),mom.end(),p_irnd);
for(unsigned int i = 0; i<ncl; i++) xi[i] = randUng(0.6,1.0);
// sort into ascending order
sort(xi.begin(), xi.end());
double ximin = xi[0];
double ximax = xi.back()-ximin;
xi[0] = 0.;
for(unsigned int i = ncl-2; i>=1; i--) xi[i+1] = (xi[i]-ximin)/ximax;
xi[1] = 1.;
yy= log(CME*CME/(mom[0].z()*mom[1].z()));
bool suceeded=false;
Energy sum2,sum3,sum4;
for(unsigned int j = 0; j<10; j++) {
sum1 = sum2 = sum3 = sum4 = ZERO;
for(unsigned int i = 0; i<ncl; i++) {
Energy tm = mom[i].z();
double ex = exp(yy*xi[i]);
sum1 += tm*ex;
sum2 += tm/ex;
sum3 += (tm*ex)*xi[i];
sum4 += (tm/ex)*xi[i];
}
double fy = alog-log(sum1*sum2/GeV2);
double dd = (sum3*sum2 - sum1*sum4)/(sum1*sum2);
double dyy = fy/dd;
if(abs(dyy/yy) < eps) {
yy += dyy;
suceeded=true;
break;
}
yy += dyy;
}
if(suceeded){
break;
}
if(its > 100) eps *= 10.;
}
if(its==_maxtries){
return false;
}
// throw Exception() << "Can't generate soft underlying event in "
// << "UA5Handler::generateCylindricalPS"
// << Exception::eventerror;
double zz = log(CME/sum1);
for(unsigned int i = 0; i<ncl; i++) {
Energy tm = mom[i].z();
double E1 = exp(zz + yy*xi[i]);
mom[i].setZ(0.5*tm*(1./E1-E1));
mom[i].setE( 0.5*tm*(1./E1+E1));
softGluons[i]=mom[i];
}
return true;
}
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();
+ setRemMasses(diquarks);
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_ << ladderPower_<< ladderNorm_ << ladderMult_ << ladderbFactor_ << pomeronStructure_
<< ounit(mg_,GeV) << ounit(ptmin_,GeV) << ounit(beta_,sqr(InvGeV))
<< allowTop_ << allowLeptons_
<< multiPeriph_ << valOfN_ << initTotRap_ << PtDistribution_;
}
void HwRemDecayer::persistentInput(PersistentIStream & is, int) {
is >> iunit(_kinCutoff, GeV) >> _range >> _zbin >> _ybin
>> _nbinmax >> _alphaS >> _alphaEM >> DISRemnantOpt_
>> maxtrySoft_ >> colourDisrupt_ >> ladderPower_ >> ladderNorm_ >> ladderMult_ >> ladderbFactor_ >> pomeronStructure_
>> iunit(mg_,GeV) >> iunit(ptmin_,GeV) >> iunit(beta_,sqr(InvGeV))
>> allowTop_ >> allowLeptons_
>> multiPeriph_ >> valOfN_ >> initTotRap_ >> PtDistribution_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<HwRemDecayer,RemnantDecayer>
describeHerwigHwRemDecayer("Herwig::HwRemDecayer", "HwShower.so");
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> interaceladderPower
("ladderPower",
"The power factor in the ladder parameterization.",
&HwRemDecayer::ladderPower_,
1.0, -5.0, 10.0,
false, false, Interface::limited);
static Parameter<HwRemDecayer,double> interfaceladderNorm
("ladderNorm",
"The normalization factor in the ladder parameterization",
&HwRemDecayer::ladderNorm_,
1.0, 0.0, 10.0,
false, false, Interface::limited);
static Parameter<HwRemDecayer,double> interfaceladderMult
("ladderMult",
"The ladder multiplicity factor ",
&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> interfaceAllowLeptons
("AllowLeptons",
"Allow charged leptons in the hadron",
&HwRemDecayer::allowLeptons_, false, false, false);
static SwitchOption interfaceAllowLeptonsNo
(interfaceAllowLeptons,
"No",
"Don't allow them",
false);
static SwitchOption interfaceAllowLeptonsYes
(interfaceAllowLeptons,
"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);
static Switch<HwRemDecayer,unsigned int> interfacePtDistribution
("PtDistribution",
"Options for different pT generation methods",
&HwRemDecayer::PtDistribution_, 0, false, false);
static SwitchOption interfacePtDistributionDefault
(interfacePtDistribution,
"Default",
"Default generation of pT",
0);
static SwitchOption interfacePtDistributionOrdered
(interfacePtDistribution,
"Ordered",
"Ordered generation of pT,where the first pT is the hardest",
1);
static SwitchOption interfacePtDistributionStronglyOrdered
(interfacePtDistribution,
"StronglyOrdered",
"Strongly ordered generation of pT",
2);
static SwitchOption interfacePtDistributionFlat
(interfacePtDistribution,
"Flat",
"Sample from a flat pT distribution",
3);
static SwitchOption interfacePtDistributionFlatOrdered
(interfacePtDistribution,
"FlatOrdered",
"First pT normal, then flat",
4);
static SwitchOption interfacePtDistributionFlatOrdered2
(interfacePtDistribution,
"FlatOrdered2",
"First pT normal, then flat but steep",
5);
}
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 allowed as a parton in the remnant 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 if(ChargedLeptonMatcher::Check(*parton)) {
if(!allowLeptons_)
throw Exception() << "Charged leptons not allowed as a parton in the remnant handling, please "
<< "use a PDF which does not contain top for the remnant"
<< " handling or allow leptons in the remnant using\n"
<< " set " << fullName() << ":AllowLeptons 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;
}
diff --git a/PDF/HwRemDecayer.h b/PDF/HwRemDecayer.h
--- a/PDF/HwRemDecayer.h
+++ b/PDF/HwRemDecayer.h
@@ -1,753 +1,753 @@
// -*- C++ -*-
//
// HwRemDecayer.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_HwRemDecayer_H
#define HERWIG_HwRemDecayer_H
//
// This is the declaration of the HwRemDecayer class.
//
#include "ThePEG/PDT/RemnantDecayer.h"
#include "ThePEG/Handlers/EventHandler.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/EventRecord/SubProcess.h"
#include "ThePEG/PDF/BeamParticleData.h"
#include "Herwig/Shower/ShowerAlpha.h"
#include "Herwig/PDT/StandardMatchers.h"
#include "ThePEG/PDT/StandardMatchers.h"
#include "HwRemDecayer.fh"
namespace Herwig {
using namespace ThePEG;
/**
* The HwRemDecayer class is responsible for the decay of the remnants. Additional
* secondary scatters have to be evolved backwards to a gluon, the
* first/hard interaction has to be evolved back to a valence quark.
* This is all generated inside this class,
* which main methods are then called by the ShowerHandler.
*
* A simple forced splitting algorithm is used.
* This takes the Remnant object produced from the PDF and backward
* evolution (hadron - parton) and produce partons with the remaining
* flavours and with the correct colour connections.
*
* The algorithim operates by starting with the parton which enters the hard process.
* If this is from the sea there is a forced branching to produce the antiparticle
* from a gluon branching. If the parton entering the hard process was a gluon, or
* a gluon was produced from the first step of the algorithm, there is then a further
* branching back to a valence parton. After these partons have been produced a quark or
* diquark is produced to give the remaining valence content of the incoming hadron.
*
* The forced branching are generated using a scale between QSpac and EmissionRange times
* the minimum scale. The energy fractions are then distributed using
* \f[\frac{\alpha_S}{2\pi}\frac{P(z)}{z}f(x/z,\tilde{q})\f]
* with the massless splitting functions.
*
* \author Manuel B\"ahr
*
* @see \ref HwRemDecayerInterfaces "The interfaces"
* defined for HwRemDecayer.
*/
class HwRemDecayer: public RemnantDecayer {
public:
/** Typedef to store information about colour partners */
typedef vector<pair<tPPtr, tPPtr> > PartnerMap;
public:
/**
* The default constructor.
*/
HwRemDecayer() : allowTop_(false), allowLeptons_(false),
multiPeriph_(true), quarkPair_(false),
ptmin_(-1.*GeV), beta_(ZERO),
maxtrySoft_(10),
colourDisrupt_(1.0),
ladderbFactor_(0.0),
ladderPower_(-0.08),
ladderNorm_(1.0),
ladderMult_(1.0),
gaussWidth_(0.1),
valOfN_(0),
initTotRap_(0),
_kinCutoff(0.75*GeV),
_forcedSplitScale(2.5*GeV),
_range(1.1), _zbin(0.05),_ybin(0.),
_nbinmax(100), DISRemnantOpt_(0),
PtDistribution_(0),
pomeronStructure_(0), mg_(ZERO) {}
/** @name Virtual functions required by the Decayer class. */
//@{
/**
* Check if this decayer can perfom the decay specified by the
* given decay mode.
* @return true if this decayer can handle the given mode, otherwise false.
*/
virtual bool accept(const DecayMode &) const {
return true;
}
/**
* Return true if this decayer can handle the extraction of the \a
* extracted parton from the given \a particle.
*/
virtual bool canHandle(tcPDPtr particle, tcPDPtr parton) const;
/**
* Return true if this decayed can extract more than one parton from
* a particle.
*/
virtual bool multiCapable() const {
return true;
}
/**
* Perform a decay for a given DecayMode and a given Particle instance.
* @param dm the DecayMode describing the decay.
* @param p the Particle instance to be decayed.
* @param step the step we are working on.
* @return a ParticleVector containing the decay products.
*/
virtual ParticleVector decay(const DecayMode & dm, const Particle & p, Step & step) const;
//@}
public:
/**
* struct that is used to catch exceptions which are thrown
* due to energy conservation issues of additional soft scatters
*/
struct ExtraSoftScatterVeto {};
/** @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();
/**
* Do several checks and initialization, for remnantdecay inside ShowerHandler.
*/
void initialize(pair<tRemPPtr, tRemPPtr> rems, tPPair beam, Step & step,
Energy forcedSplitScale);
/**
* Initialize the soft scattering machinery.
* @param ptmin = the pt cutoff used in the UE model
* @param beta = slope of the soft pt-spectrum
*/
void initSoftInteractions(Energy ptmin, InvEnergy2 beta);
/**
* Perform the acual forced splitting.
* @param partons is a pair of ThePEG::Particle pointers which store the final
* partons on which the shower ends.
* @param pdfs are pointers to the pdf objects for both beams
* @param first is a flage wether or not this is the first or a secondary interation
*/
void doSplit(pair<tPPtr, tPPtr> partons, pair<tcPDFPtr, tcPDFPtr> pdfs, bool first);
/**
* Perform the final creation of the diquarks. Set the remnant masses and do
* all colour connections.
* @param colourDisrupt = variable to control how many "hard" scatters
* are colour isolated
* @param softInt = parameter for the number of soft scatters
*/
void finalize(double colourDisrupt=0.0, unsigned int softInt=0);
/**
* Find the children
*/
void findChildren(tPPtr,vector<PPtr> &) const;
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const {return new_ptr(*this);}
/** 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 {return new_ptr(*this);}
//@}
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() {
Interfaced::doinit();
_ybin=0.25/_zbin;
mg_ = getParticleData(ParticleID::g)->constituentMass();
}
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
HwRemDecayer & operator=(const HwRemDecayer &) = delete;
public:
/**
* Simple struct to store info about baryon quark and di-quark
* constituents.
*/
struct HadronContent {
/**
* manually extract the valence flavour \a id.
*/
inline void extract(int id) {
for(unsigned int i=0; i<flav.size(); i++) {
if(id == sign*flav[i]){
if(hadron->id() == ParticleID::gamma ||
(hadron->id() == ParticleID::pomeron && pomeronStructure==1) ||
hadron->id() == ParticleID::reggeon) {
flav[0] = id;
flav[1] = -id;
extracted = 0;
flav.resize(2);
}
else if (hadron->id() == ParticleID::pomeron && pomeronStructure==0) {
extracted = 0;
}
else {
extracted = i;
}
break;
}
}
}
/**
* Return a proper particle ID assuming that \a id has been removed
* from the hadron.
*/
long RemID() const;
/**
* Method to determine whether \a parton is a quark from the sea.
* @return TRUE if \a parton is neither a valence quark nor a gluon.
*/
bool isSeaQuark(tcPPtr parton) const {
return ((parton->id() != ParticleID::g) && ( !isValenceQuark(parton) ) );
}
/**
* Method to determine whether \a parton is a valence quark.
*/
bool isValenceQuark(tcPPtr parton) const {
return isValenceQuark(parton->id());
}
/**
* Method to determine whether \a parton is a quark from the sea.
* @return TRUE if \a parton is neither a valence quark nor a gluon.
*/
bool isSeaQuarkData(tcPDPtr partonData) const {
return ((partonData->id() != ParticleID::g) && ( !isValenceQuarkData(partonData) ) );
}
/**
* Method to determine whether \a parton is a valence quark.
*/
bool isValenceQuarkData(tcPDPtr partonData) const {
int id(sign*partonData->id());
return find(flav.begin(),flav.end(),id) != flav.end();
}
/**
* Method to determine whether \a parton is a valence quark.
*/
bool isValenceQuark(int id) const {
return find(flav.begin(),flav.end(),sign*id) != flav.end();
}
/** The valence flavours of the corresponding baryon. */
vector<int> flav;
/** The array index of the extracted particle. */
int extracted;
/** -1 if the particle is an anti-particle. +1 otherwise. */
int sign;
/** The ParticleData objects of the hadron */
tcPDPtr hadron;
/** Pomeron treatment */
unsigned int pomeronStructure;
};
/**
* Return the hadron content objects for the incoming particles.
*/
const pair<HadronContent, HadronContent>& content() const {
return theContent;
}
/**
* Return a HadronContent struct from a PPtr to a hadron.
*/
HadronContent getHadronContent(tcPPtr hadron) const;
/**
* Set the hadron contents.
*/
void setHadronContent(tPPair beam) {
theContent.first = getHadronContent(beam.first);
theContent.second = getHadronContent(beam.second);
}
private:
/**
* Do the forced Splitting of the Remnant with respect to the
* extracted parton \a parton.
* @param parton = PPtr to the parton going into the subprocess.
* @param content = HadronContent struct to keep track of flavours.
* @param rem = Pointer to the ThePEG::RemnantParticle.
* @param used = Momentum vector to keep track of remaining momenta.
* @param partners = Vector of pairs filled with tPPtr to the particles
* which should be colour connected.
* @param pdf pointer to the PDF Object which is used for this particle
* @param first = Flag for the first interaction.
*/
void split(tPPtr parton, HadronContent & content, tRemPPtr rem,
Lorentz5Momentum & used, PartnerMap & partners, tcPDFPtr pdf, bool first);
/**
* Merge the colour lines of two particles
* @param p1 = Pointer to particle 1
* @param p2 = Pointer to particle 2
* @param anti = flag to indicate, if (anti)colour was extracted as first parton.
*/
void mergeColour(tPPtr p1, tPPtr p2, bool anti) const;
/**
* Set the colour connections.
* @param partners = Object that holds the information which particles to connect.
* @param anti = flag to indicate, if (anti)colour was extracted as first parton.
* @param disrupt parameter for disruption of the colour structure
*/
void fixColours(PartnerMap partners, bool anti, double disrupt) const;
/**
* Set the momenta of the Remnants properly and boost the decay particles.
*/
- void setRemMasses() const;
+ void setRemMasses(PPair diquarks) const;
/**
* This creates a parton from the remaining flavours of the hadron. The
* last parton used was a valance parton, so only 2 (or 1, if meson) flavours
* remain to be used.
*/
PPtr finalSplit(const tRemPPtr rem, long remID,
Lorentz5Momentum usedMomentum) const {
// Create the remnant and set its momentum, also reset all of the decay
// products from the hadron
PPtr remnant = new_ptr(Particle(getParticleData(remID)));
Lorentz5Momentum prem(rem->momentum()-usedMomentum);
prem.setMass(getParticleData(remID)->constituentMass());
prem.rescaleEnergy();
remnant->set5Momentum(prem);
// Add the remnant to the step, but don't do colour connections
thestep->addDecayProduct(rem,remnant,false);
return remnant;
}
/**
* This takes the particle and find a splitting for np -> p + child and
* creates the correct kinematics and connects for such a split. This
* Splitting has an upper bound on qtilde given by the energy argument
* @param rem The Remnant
* @param child The PDG code for the outgoing particle
* @param oldQ The maximum scale for the evolution
* @param oldx The fraction of the hadron's momentum carried by the last parton
* @param pf The momentum of the last parton at input and after branching at output
* @param p The total emitted momentum
* @param content The content of the hadron
*/
PPtr forceSplit(const tRemPPtr rem, long child, Energy &oldQ, double &oldx,
Lorentz5Momentum &pf, Lorentz5Momentum &p,
HadronContent & content) const;
/**
* Check if a particle is a parton from a hadron or not
* @param parton The parton to be tested
*/
bool isPartonic(tPPtr parton) const;
/** @name Soft interaction methods. */
//@{
/**
* Produce pt values according to dN/dp_T = N p_T exp(-beta_*p_T^2)
*/
Energy softPt() const;
/**
* Get the 2 pairs of 5Momenta for the scattering. Needs calling of
* initSoftInteractions.
*/
void softKinematics(Lorentz5Momentum &r1, Lorentz5Momentum &r2,
Lorentz5Momentum &g1, Lorentz5Momentum &g2) const;
/**
* Create N soft gluon interactions
*/
void doSoftInteractions(unsigned int N){
if(!multiPeriph_){
doSoftInteractions_old(N);} //outdated model for soft interactions
else{
doSoftInteractions_multiPeriph(N); // Multiperipheral model
}
}
/**
* Create N soft gluon interactions (old version)
*/
void doSoftInteractions_old(unsigned int N);
/**
* Create N soft gluon interactions with multiperhpheral kinematics
*/
void doSoftInteractions_multiPeriph(unsigned int N);
/**
* Phase space generation for the ladder partons
*/
bool doPhaseSpaceGenerationGluons(vector<Lorentz5Momentum> &softGluons, Energy energy, unsigned int &its)
const;
/**
* This returns the rotation matrix needed to rotate p into the z axis
*/
LorentzRotation rotate(const LorentzMomentum &p) const;
/**
* Methods to generate random distributions also all stolen form UA5Handler
**/
template <typename T>
inline T gaussDistribution(T mean, T stdev) const{
double x = rnd();
x = sqrt(-2.*log(x));
double y;
randAzm(x,x,y);
return mean + stdev*x;
}
/**
* This returns a random number with a flat distribution
* [-A,A] plus gaussian tail with stdev B
* TODO: Should move this to Utilities
* @param A The width of the flat part
* @param B The standard deviation of the gaussian tail
* @return the randomly generated value
*/
inline double randUng(double A, double B) const{
double prun;
if(A == 0.) prun = 0.;
else prun = 1./(1.+B*1.2533/A);
if(rnd() < prun) return 2.*(rnd()-0.5)*A;
else {
double temp = gaussDistribution(0.,B);
if(temp < 0) return temp - abs(A);
else return temp + abs(A);
}
}
template <typename T>
inline void randAzm(T pt, T &px, T &py) const{
double c,s,cs;
while(true) {
c = 2.*rnd()-1.;
s = 2.*rnd()-1.;
cs = c*c+s*s;
if(cs <= 1.&&cs!=0.) break;
}
T qt = pt/cs;
px = (c*c-s*s)*qt;
py = 2.*c*s*qt;
}
inline Energy randExt(Energy AM0,InvEnergy B) const{
double r = rnd();
// Starting value
Energy am = AM0-log(r)/B;
for(int i = 1; i<20; ++i) {
double a = exp(-B*(am-AM0))/(1.+B*AM0);
double f = (1.+B*am)*a-r;
InvEnergy df = -B*B*am*a;
Energy dam = -f/df;
am += dam;
if(am<AM0) am = AM0 + .001*MeV;
if(abs(dam) < .001*MeV) break;
}
return am;
}
/**
* Method to add a particle to the step
* @param parent = pointer to the parent particle
* @param id = Particle ID of the newly created particle
* @param p = Lorentz5Momentum of the new particle
*/
tPPtr addParticle(tcPPtr parent, long id, Lorentz5Momentum p) const;
//@}
/**
* A flag which indicates, whether the extracted valence quark was a
* anti particle.
*/
pair<bool, bool> theanti;
/**
* variable to sum up the x values of the extracted particles
*/
pair<double, double> theX;
/**Pair of HadronContent structs to know about the quark content of the beams*/
pair<HadronContent, HadronContent> theContent;
/**Pair of Lorentz5Momentum to keep track of the forced splitting product momenta*/
pair<Lorentz5Momentum, Lorentz5Momentum> theUsed;
/**
* Pair of PartnerMap's to store the particles, which will be colour
* connected in the end.
*/
pair<PartnerMap, PartnerMap> theMaps;
/**
* Variable to hold a pointer to the current step. The variable is used to
* determine, wether decay(const DecayMode & dm, const Particle & p, Step & step)
* has been called in this event or not.
*/
StepPtr thestep;
/**
* Pair of Remnant pointers. This is needed to boost
* in the Remnant-Remnant CMF after all have been decayed.
*/
pair<RemPPtr, RemPPtr> theRems;
/**
* The beam particle data for the current incoming hadron
*/
mutable tcPPtr theBeam;
/**
* the beam data
*/
mutable Ptr<BeamParticleData>::const_pointer theBeamData;
/**
* The PDF for the current initial-state shower
*/
mutable tcPDFPtr _pdf;
private:
/**
* Switch to control handling of top quarks in proton
*/
bool allowTop_;
/**
* Switch to control handling of charged leptons in proton
*/
bool allowLeptons_;
/**
* Switch to control using multiperipheral kinemaics
*/
bool multiPeriph_;
/**
* True if kinematics is to be calculated for quarks
*/
bool quarkPair_;
/** @name Soft interaction variables. */
//@{
/**
* Pair of soft Remnant pointers, i.e. Diquarks.
*/
tPPair softRems_;
/**
* ptcut of the UE model
*/
Energy ptmin_;
/**
* slope of the soft pt-spectrum: dN/dp_T = N p_T exp(-beta*p_T^2)
*/
InvEnergy2 beta_;
/**
* Maximum number of attempts for the regeneration of an additional
* soft scattering, before the number of scatters is reduced.
*/
unsigned int maxtrySoft_;
/**
* Variable to store the relative number of colour disrupted
* connections to additional soft subprocesses.
*/
double colourDisrupt_;
/**
* Variable to store the additive factor of the
multiperipheral ladder multiplicity.
*/
double ladderbFactor_;
/**
* Variable of the parameterization of the ladder multiplicity.
*/
double ladderPower_;
/**
* Variable of the parameterization of the ladder multiplicity.
*/
double ladderNorm_;
double ladderMult_;
/**
* Variable to store the gaussian width of the
* fluctuation of the longitudinal momentum
* fraction.
*/
double gaussWidth_;
/**
* Variable to store the current total multiplicity
of a ladder.
*/
double valOfN_;
/**
* Variable to store the initial total rapidity between
of the remnants.
*/
double initTotRap_;
//@}
/** @name Forced splitting variables. */
//@{
/**
* The kinematic cut-off
*/
Energy _kinCutoff;
/**
* The PDF freezing scale as set in ShowerHandler
*/
Energy _forcedSplitScale;
/**
* Range for emission
*/
double _range;
/**
* Size of the bins in z for the interpolation
*/
double _zbin;
/**
* Size of the bins in y for the interpolation
*/
double _ybin;
/**
* Maximum number of bins for the z interpolation
*/
int _nbinmax;
/**
* Pointer to the object calculating the QCD coupling
*/
ShowerAlphaPtr _alphaS;
/**
* Pointer to the object calculating the QED coupling
*/
ShowerAlphaPtr _alphaEM;
/**
* Option for the DIS remnant
*/
unsigned int DISRemnantOpt_;
/**
* Option for the pT generation
*/
unsigned int PtDistribution_;
/**
* Option for the treatment of the pomeron structure
*/
unsigned int pomeronStructure_;
//@}
/**
* The gluon constituent mass.
*/
Energy mg_;
};
}
#endif /* HERWIG_HwRemDecayer_H */
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Tue, Nov 19, 8:25 PM (1 d, 6 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3803549
Default Alt Text
(88 KB)
Attached To
rHERWIGHG herwighg
Event Timeline
Log In to Comment