diff --git a/MatrixElement/MEMinBias.cc b/MatrixElement/MEMinBias.cc
--- a/MatrixElement/MEMinBias.cc
+++ b/MatrixElement/MEMinBias.cc
@@ -1,279 +1,295 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEMinBias class.
//
#include "MEMinBias.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
//#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Handlers/SamplerBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/MatrixElement/Tree2toNDiagram.h"
inline bool checkValence(int i,int side,Ptr<StandardEventHandler>::tptr eh){
// Inline function to check for valence quarks of the beam.
// i: pdgid of quark
// side: beam side
// eh: pointer to the eventhandler
int beam= ( side == 0 ) ? eh->incoming().first->id() : eh->incoming().second->id();
vector<int> val;
- if( beam == ParticleID::pplus || beam == ParticleID::n0 ) val = {1,2};
- if( beam == ParticleID::pbarminus || beam == ParticleID::nbar0 ) val = { -1 , -2 };
- if( val.size() == 0 )
- {cerr<<"\n\n MEMinBias: Valence Quarks not defined for pid "<<beam;assert(false);}
- for(auto v:val)if(v==i)return true;
+ switch (beam) {
+ case ParticleID::pplus :
+ case ParticleID::n0 :
+ val = { 1, 2 };
+ break;
+ case ParticleID::pbarminus:
+ case ParticleID::nbar0:
+ val = { -1, -2 };
+ break;
+ case ParticleID::piplus :
+ val = { -1, 2 };
+ break;
+ case ParticleID::piminus :
+ val = { 1, -2 };
+ break;
+ default :
+ cerr << "\n\n MEMinBias: Valence Quarks not defined for pid " << beam;
+ assert(false);
+ };
+ for(const int & v : val )
+ if(v==i) return true;
return false;
}
void MEMinBias::getDiagrams() const {
int maxflav(2);
// Pomeron data
tcPDPtr pom = getParticleData(990);
Ptr<StandardEventHandler>::tptr eh = dynamic_ptr_cast<Ptr<StandardEventHandler>::tptr>(generator()->eventHandler());
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
if(!onlyValQuarks_) add(new_ptr((Tree2toNDiagram(3), q1, pom, q2, 1, q1, 2, q2, -1)));
else if(checkValence(i,0,eh) && checkValence(j,1,eh) ) add(new_ptr((Tree2toNDiagram(3), q1, pom, q2, 1, q1, 2, q2, -1)));
//qqb -> qqb
if(!onlyValQuarks_) add(new_ptr((Tree2toNDiagram(3), q1, pom, q2b, 1, q1, 2, q2b, -2)));
else if(checkValence(i,0,eh) && checkValence(-j,1,eh) ) add(new_ptr((Tree2toNDiagram(3), q1, pom, q2b, 1, q1, 2, q2b, -2)));
//qbqb -> qbqb
if(!onlyValQuarks_) add(new_ptr((Tree2toNDiagram(3), q1b, pom, q2b, 1, q1b, 2, q2b, -3)));
else if(checkValence(-i,0,eh) && checkValence(-j,1,eh) ) add(new_ptr((Tree2toNDiagram(3), q1b, pom, q2b, 1, q1b, 2, q2b, -3)));
}
}
}
Energy2 MEMinBias::scale() const {
return sqr(Scale_);
}
int MEMinBias::nDim() const {
return 0;
}
void MEMinBias::setKinematics() {
HwMEBase::setKinematics(); // Always call the base class method first.
}
bool MEMinBias::generateKinematics(const double *) {
// generate the masses of the particles
for ( int i = 2, N = meMomenta().size(); i < N; ++i ) {
meMomenta()[i] = Lorentz5Momentum(mePartonData()[i]->constituentMass());
}
Energy q = ZERO;
try {
q = SimplePhaseSpace::
getMagnitude(sHat(), meMomenta()[2].mass(), meMomenta()[3].mass());
} catch ( ImpossibleKinematics & e ) {
return false;
}
Energy pt = ZERO;
meMomenta()[2].setVect(Momentum3( pt, pt, q));
meMomenta()[3].setVect(Momentum3(-pt, -pt, -q));
meMomenta()[2].rescaleEnergy();
meMomenta()[3].rescaleEnergy();
jacobian(1.0);
return true;
}
double MEMinBias::correctionweight() const {
// Here we calculate the weight to restore the inelastic-diffractiveXSec
// given by the MPIHandler.
// First get the eventhandler to get the current cross sections.
Ptr<StandardEventHandler>::tptr eh =
dynamic_ptr_cast<Ptr<StandardEventHandler>::tptr>(generator()->eventHandler());
// All diffractive processes make use of this ME.
// The static map can be used to collect all the sumOfWeights.
static map<XCombPtr,double> weightsmap;
weightsmap[lastXCombPtr()]=lastXComb().stats().sumWeights();
// Define static variable to keep trac of reweighting
static double rew_=1.;
static int countUpdateWeight=50;
static double sumRew=0.;
static double countN=0;
// if we produce events we count
if(eh->integratedXSec()>ZERO)sumRew+=rew_;
if(eh->integratedXSec()>ZERO)countN+=1.;
if(countUpdateWeight<countN){
// Summing all diffractive processes (various initial states)
double sum=0.;
for(auto xx:weightsmap){
sum+=xx.second;
}
double avRew=sumRew/countN;
CrossSection XS_have =eh->sampler()->maxXSec()/eh->sampler()->attempts()*sum;
CrossSection XS_wanted=MPIHandler_->nonDiffractiveXSec();
double deltaN=50;
// Cross section without reweighting: XS_norew
// XS_have = avcsNorm2*XS_norew (for large N)
// We want to determine the rew that allows to get the wanted XS.
// In deltaN points we want (left) and we get (right):
// XS_wanted*(countN+deltaN) = XS_have*countN + rew*deltaN*XS_norew
// Solve for rew:
rew_=avRew*(XS_wanted*(countN+deltaN)-XS_have*countN)/(XS_have*deltaN);
countUpdateWeight+=deltaN;
}
//Make sure we dont produce negative weights.
// TODO: write finalize method that checks if reweighting was performed correctly.
rew_=max(rew_,0.000001);
rew_=min(rew_,10000.0);
return rew_;
}
double MEMinBias::me2() const {
//tuned so it gives the correct normalization for xmin = 0.11
return csNorm_*(sqr(generator()->maximumCMEnergy())/GeV2);
}
CrossSection MEMinBias::dSigHatDR() const {
return me2()*jacobian()/sHat()*sqr(hbarc)*correctionweight();
}
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);
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<MEMinBias,HwMEBase>
describeHerwigMEMinBias("Herwig::MEMinBias", "HwMEHadron.so");
void MEMinBias::persistentOutput(PersistentOStream & os) const {
os << csNorm_ << ounit(Scale_,GeV) << MPIHandler_;
}
void MEMinBias::persistentInput(PersistentIStream & is, int) {
is >> csNorm_ >> iunit(Scale_,GeV) >> MPIHandler_;
}
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);
static Parameter<MEMinBias,Energy> interfaceScale
("Scale",
"Scale for the Min Bias matrix element.",
&MEMinBias::Scale_,GeV,
2.0*GeV, 0.0*GeV, 100.0*GeV,
false, false, Interface::limited);
static Reference<MEMinBias,UEBase> interfaceMPIHandler
("MPIHandler",
"The object that administers all additional scatterings.",
&MEMinBias::MPIHandler_, false, false, true, true);
static Switch<MEMinBias , bool> interfaceOnlyVal
("OnlyValence" ,
"Allow the dummy process to only extract valence quarks." ,
&MEMinBias::onlyValQuarks_ , false , false , false );
static SwitchOption interfaceOnlyValYes
( interfaceOnlyVal , "Yes" , "" , true );
static SwitchOption interfaceOnlyValNo
( interfaceOnlyVal , "No" , "" , false );
}
diff --git a/PDF/HwRemDecayer.cc b/PDF/HwRemDecayer.cc
--- a/PDF/HwRemDecayer.cc
+++ b/PDF/HwRemDecayer.cc
@@ -1,2020 +1,2028 @@
// -*- 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/Utilities/EnumParticles.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
+ // 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 )) {
+ 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(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 = theprocessed[0] ? diquarks.first : diquarks.second;
vector<PPtr> progenitors;
for(unsigned int ix=0;ix<rem->children().size();++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;
hadron=rem->parents()[0];
while(hadron->id()==ParticleID::Remnant) hadron=hadron->parents()[0];
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(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);
}
}
+ else {
+ hc.sign = id < 0? -1: 1;
+ hc.flav.push_back(-(id = abs(id)/10)%10);
+ hc.flav.push_back((id = abs(id)/10)%10);
+ if(hc.flav.back()%2!=0) {
+ for(int & i : hc.flav) i *=-1;
+ }
+ }
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/Tests/python/make_input_files.py.in b/Tests/python/make_input_files.py.in
--- a/Tests/python/make_input_files.py.in
+++ b/Tests/python/make_input_files.py.in
@@ -1,1614 +1,1625 @@
#! @PYTHON@
# -*- mode: python -*-
from __future__ import print_function
import logging,sys,os
from string import Template
from HerwigInputs import *
import sys
if sys.version_info[:3] < (2,4,0):
print ("rivet scripts require Python version >= 2.4.0... exiting")
sys.exit(1)
if __name__ == "__main__":
import logging
from optparse import OptionParser, OptionGroup
parser = OptionParser(usage="%prog name [...]")
simulation=""
numberOfAddedProcesses=0
def addProcess(thefactory,theProcess,Oas,Oew,scale,mergedlegs,NLOprocesses):
global numberOfAddedProcesses
global simulation
numberOfAddedProcesses+=1
res ="set "+thefactory+":OrderInAlphaS "+Oas+"\n"
res+="set "+thefactory+":OrderInAlphaEW "+Oew+"\n"
res+="do "+thefactory+":Process "+theProcess+" "
if ( mergedlegs != 0 ):
if simulation!="Merging":
print ("simulation is not Merging, trying to add merged legs.")
sys.exit(1)
res+="["
for j in range(mergedlegs):
res+=" j "
res+="]"
res+="\n"
if (NLOprocesses!=0):
if simulation!="Merging":
print ("simulation is not Merging, trying to add NLOProcesses.")
sys.exit(1)
res+="set MergingFactory:NLOProcesses %s \n" % NLOprocesses
if ( scale != "" ):
res+="set "+thefactory+":ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/"+scale+"\n"
return res
addedBRReweighter=False
def addBRReweighter():
global addedBRReweighter
if(addedBRReweighter):
logging.error("Can only add BRReweighter once.")
sys.exit(1)
res="create Herwig::BranchingRatioReweighter /Herwig/Generators/BRReweighter\n"
res+="insert /Herwig/Generators/EventGenerator:EventHandler:PostHadronizationHandlers 0 /Herwig/Generators/BRReweighter\n"
addedBRReweighter=True
return res
selecteddecaymode=False
def selectDecayMode(particle,decaymodes):
global selecteddecaymode
res="do /Herwig/Particles/"+particle+":SelectDecayModes"
for decay in decaymodes:
res+=" /Herwig/Particles/"+particle+"/"+decay
res+="\n"
selecteddecaymode=True
return res
ME_Upsilon = """\
create Herwig::MEee2VectorMeson /Herwig/MatrixElements/MEUpsilon HwMELepton.so
set /Herwig/MatrixElements/MEUpsilon:VectorMeson /Herwig/Particles/Upsilon(4S)
set /Herwig/MatrixElements/MEUpsilon:Coupling 96.72794
""" + insert_ME("MEUpsilon")
(opts, args) = parser.parse_args()
## Check args
if len(args) != 1:
logging.error("Must specify at least input file")
sys.exit(1)
name = args[0]
print (name)
# work out name and type of the collider
(collider,have_hadronic_collider) = identifyCollider(name)
# workout the type of simulation
(simulation,templateName,parameterName,parameters)=identifySimulation(name,collider,have_hadronic_collider)
if simulation=="Merging" :
thefactory="MergingFactory"
else :
thefactory="Factory"
# settings for four flavour scheme
fourFlavour="""
read Matchbox/FourFlavourScheme.in
{bjetgroup}
set /Herwig/Cuts/MatchboxJetMatcher:Group bjet
""".format(bjetgroup=particlegroup(thefactory,'bjet','b','bbar','c', 'cbar',
's','sbar','d','dbar','u','ubar','g'))
# work out the process and parameters
process=StringBuilder()
# DIS
if(collider=="DIS" and "Photo" not in name) :
if(simulation=="") :
if "NoME" in name :
process = StringBuilder("set /Herwig/Shower/ShowerHandler:HardEmission None")
parameterName=parameterName.replace("NoME-","")
parameterName=parameterName.replace("DIS-" ,"")
if "CC" in parameterName :
process += insert_ME("MEDISCC")
else :
process += insert_ME("MEDISNC")
elif(simulation=="Powheg") :
if "CC" in parameterName :
process = StringBuilder(insert_ME("PowhegMEDISCC"))
else :
process = StringBuilder(insert_ME("PowhegMEDISNC"))
elif(simulation=="Matchbox" ) :
if "CC" in name :
if "e-" in parameterName :
process = StringBuilder(addProcess(thefactory,"e- p -> nu_e j","0","2","",0,0))
else :
process = StringBuilder(addProcess(thefactory,"e+ p -> nu_ebar j","0","2","",0,0))
else :
if "e-" in parameterName :
process = StringBuilder(addProcess(thefactory,"e- p -> e- j","0","2","",0,0))
else :
process = StringBuilder(addProcess(thefactory,"e+ p -> e+ j","0","2","",0,0))
elif(simulation=="Merging" ) :
if "CC" in name :
if "e-" in parameterName :
process = StringBuilder(addProcess(thefactory,"e- p -> e- j","0","2","",2,2))
else :
process = StringBuilder(addProcess(thefactory,"e+ p -> e+ j","0","2","",2,2))
else :
if "e-" in parameterName :
process = StringBuilder(addProcess(thefactory,"e- p -> nu_e j","0","2","",2,2))
else :
process = StringBuilder(addProcess(thefactory,"e+ p -> nu_ebar j","0","2","",2,2))
Q2Min=1.
Q2Max=1000000.
if "VeryLow" in name :
Q2Max=20.
parameterName=parameterName.replace("-VeryLowQ2","")
elif "Low" in name :
Q2Min=20.
Q2Max=100.
parameterName=parameterName.replace("-LowQ2","")
elif "Med" in name :
Q2Min=100.
Q2Max=1000.
parameterName=parameterName.replace("-MedQ2","")
elif "High" in name :
Q2Min=1000.
parameterName=parameterName.replace("-HighQ2","")
if "CC" in name :
process+="set /Herwig/Cuts/ChargedCurrentCut:MaxQ2 %s\nset /Herwig/Cuts/ChargedCurrentCut:MinQ2 %s\n" %(Q2Max,Q2Min)
else :
process+="set /Herwig/Cuts/NeutralCurrentCut:MaxQ2 %s\nset /Herwig/Cuts/NeutralCurrentCut:MinQ2 %s\n" %(Q2Max,Q2Min)
# DIS photoproduction
elif(collider=="DIS" and "Photo" in name) :
assert(simulation=="")
ecms=float(parameterName.split("-")[1])
cuts=[3.,6.,10.,ecms]
if "Direct" in parameterName :
parameterName=parameterName.replace("Direct-","")
process = StringBuilder(insert_ME("MEGammaP2Jets",None,"Process","SubProcess"))
elif "Resolved" in parameterName :
process = StringBuilder(insert_ME("MEQCD2to2"))
parameterName=parameterName.replace("Resolved-","")
for i in range(1,len(cuts)) :
tstring = "-Jets-%s"%i
if tstring in parameterName :
process+=jet_kt_cut(cuts[i-1],cuts[i])
parameterName=parameterName.replace(tstring,"")
# EE
elif(collider=="EE") :
if(simulation=="") :
if "gg" in parameterName :
process = StringBuilder("create Herwig::MEee2Higgs2SM /Herwig/MatrixElements/MEee2Higgs2SM\n")
process+=insert_ME("MEee2Higgs2SM","Gluon","Allowed")
elif "LL" in parameterName :
process = StringBuilder(insert_ME("MEee2gZ2ll"))
process += "set /Herwig/MatrixElements/MEee2gZ2ll:Allowed Charged\n"
elif "WW" in parameterName :
process = StringBuilder(insert_ME("MEee2VV"))
process += "set /Herwig/MatrixElements/MEee2VV:Process WW\n"
else :
process = StringBuilder(insert_ME("MEee2gZ2qq"))
try :
ecms = float(parameterName)
if(ecms<=3.75) :
process+= "set /Herwig/MatrixElements/MEee2gZ2qq:MaximumFlavour 3\n"
elif(ecms<=10.6) :
process+= "set /Herwig/MatrixElements/MEee2gZ2qq:MaximumFlavour 4\n"
except :
pass
elif(simulation=="Powheg") :
if "LL" in parameterName :
process = StringBuilder(insert_ME("PowhegMEee2gZ2ll"))
process += "set /Herwig/MatrixElements/PowhegMEee2gZ2ll:Allowed Charged\n"
else :
process = StringBuilder(insert_ME("PowhegMEee2gZ2qq"))
try :
ecms = float(parameterName)
if(ecms<=3.75) :
process+= "set /Herwig/MatrixElements/PowhegMEee2gZ2qq:MaximumFlavour 3\n"
elif(ecms<=10.6) :
process+= "set /Herwig/MatrixElements/PowhegMEee2gZ2qq:MaximumFlavour 4\n"
except :
pass
elif(simulation=="Matchbox" ) :
try :
ecms = float(parameterName)
if(ecms<=3.75) :
process = StringBuilder(addProcess(thefactory,"e- e+ -> u ubar","0","2","",0,0))
process+=addProcess(thefactory,"e- e+ -> d dbar","0","2","",0,0)
process+=addProcess(thefactory,"e- e+ -> s sbar","0","2","",0,0)
elif(ecms<=10.6) :
process = StringBuilder(addProcess(thefactory,"e- e+ -> u ubar","0","2","",0,0))
process+=addProcess(thefactory,"e- e+ -> d dbar","0","2","",0,0)
process+=addProcess(thefactory,"e- e+ -> c cbar","0","2","",0,0)
process+=addProcess(thefactory,"e- e+ -> s sbar","0","2","",0,0)
else :
process = StringBuilder(addProcess(thefactory,"e- e+ -> j j","0","2","",0,0))
except:
process = StringBuilder(addProcess(thefactory,"e- e+ -> j j","0","2","",0,0))
elif(simulation=="Merging" ) :
try :
ecms = float(parameterName)
if(ecms<=10.1) :
process = StringBuilder(addProcess(thefactory,"e- e+ -> j j","0","2","",2,2))
process+="read Matchbox/FourFlavourScheme.in"
else :
process = StringBuilder(addProcess(thefactory,"e- e+ -> j j","0","2","",2,2))
except:
process = StringBuilder(addProcess(thefactory,"e- e+ -> j j","0","2","",2,2))
# EE-Gamma
elif(collider=="EE-Gamma") :
if(simulation=="") :
if("mumu" in parameterName) :
process = StringBuilder(insert_ME("MEgg2ff","Muon"))
process +="set /Herwig/Cuts/Cuts:MHatMin 3.\n"
parameterName=parameterName.replace("Direct-","")
elif( "tautau" in parameterName) :
process = StringBuilder(insert_ME("MEgg2ff","Tau"))
process +="set /Herwig/Cuts/Cuts:MHatMin 3.\n"
parameterName=parameterName.replace("Direct-","")
elif( "Jets" in parameterName) :
if("Direct" in parameterName ) :
process = StringBuilder(insert_ME("MEgg2ff","Quarks"))
parameterName=parameterName.replace("Direct-","")
elif("Single-Resolved" in parameterName ) :
process = StringBuilder(insert_ME("MEGammaP2Jets",None,"Process","SubProcess"))
process+= insert_ME("MEGammaP2Jets",None,"Process","SubProcess2")
parameterName=parameterName.replace("Single-Resolved-","")
else :
process = StringBuilder(insert_ME("MEQCD2to2"))
parameterName=parameterName.replace("Double-Resolved-","")
process+="insert /Herwig/Cuts/Cuts:OneCuts[0] /Herwig/Cuts/JetKtCut"
process+="set /Herwig/Cuts/JetKtCut:MinKT 3."
else :
print ("process not supported for Gamma Gamma processes at EE")
quit()
else :
print ("Only internal matrix elements currently supported for Gamma Gamma processes at EE")
quit()
elif(collider=="GammaGamma") :
if(simulation=="") :
if("mumu" in parameterName) :
process = StringBuilder(insert_ME("MEgg2ff"))
process +="set /Herwig/MatrixElements/MEgg2ff:Process Muon\n"
process +="set /Herwig/Cuts/Cuts:MHatMin 3.\n"
else :
print ("process not supported for Gamma Gamma processes at EE")
quit()
else :
print ("Only internal matrix elements currently supported for Gamma Gamma processes at EE")
quit()
# TVT
elif(collider=="TVT") :
process = StringBuilder("set /Herwig/Generators/EventGenerator:EventHandler:BeamB /Herwig/Particles/pbar-\n")
ecms=1960.
if "Run-II" in parameterName : ecms = 1960.0
elif "Run-I" in parameterName : ecms = 1800.0
elif "900" in parameterName : ecms = 900.0
elif "630" in parameterName : ecms = 630.0
elif "300" in parameterName : ecms = 300.0
process+=collider_lumi(ecms)
if(simulation=="") :
if "PromptPhoton" in parameterName :
process+=insert_ME("MEGammaJet")
process+="set /Herwig/Cuts/PhotonKtCut:MinKT 15.\n"
elif "DiPhoton-GammaGamma" in parameterName :
process+=insert_ME("MEGammaGamma")
process+="set /Herwig/Cuts/PhotonKtCut:MinKT 5.\n"
parameterName=parameterName.replace("-GammaGamma","")
elif "DiPhoton-GammaJet" in parameterName :
process+=insert_ME("MEGammaJet")
process+="set /Herwig/Cuts/PhotonKtCut:MinKT 5.\n"
parameterName=parameterName.replace("-GammaJet","")
elif "UE" in parameterName :
if "Dipole" in parameters["shower"]:
process+="read snippets/MB-DipoleShower.in\n"
else:
process+="read snippets/MB.in\n"
process+="read snippets/Diffraction.in\n"
process += "set /Herwig/Decays/DecayHandler:LifeTimeOption 0\n"
process += "set /Herwig/Decays/DecayHandler:MaxLifeTime 10*mm\n"
elif "Jets" in parameterName :
process+=insert_ME("MEQCD2to2")
process+="set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0\n"
if "DiJets" in name :
process +=jet_kt_cut( 30.)
cuts=[100.,300.,600.,900.,ecms]
for i in range(1,len(cuts)) :
tstring = "-DiJets-%s"%i
if tstring in parameterName :
process+=mhat_cut(cuts[i-1],cuts[i])
parameterName=parameterName.replace(tstring,"-DiJets")
else :
if "Run" in parameterName :
cuts=[5.,20.,40.,80.,160.,320.]
elif "300" in parameterName :
cuts=[5.,]
elif "630" in parameterName :
cuts=[5.,20.,40.]
elif "900" in parameterName :
cuts=[5.,]
cuts.append(ecms)
for i in range(1,len(cuts)) :
tstring = "-Jets-%s"%i
if tstring in parameterName :
process+=jet_kt_cut(cuts[i-1],cuts[i])
parameterName=parameterName.replace(tstring,"-Jets")
elif "Run-I-WZ" in parameterName :
process+=insert_ME("MEqq2W2ff","Electron")
process+=insert_ME("MEqq2gZ2ff","Electron")
elif "Run-II-W" in parameterName or "Run-I-W" in parameterName :
process+=insert_ME("MEqq2W2ff","Electron")
elif "Run-II-Z-e" in parameterName or "Run-I-Z" in parameterName :
process +=insert_ME("MEqq2gZ2ff","Electron")
elif "Run-II-Z-LowMass-mu" in parameterName :
process +=insert_ME("MEqq2gZ2ff","Muon")
process+=addLeptonPairCut("25","70")
elif "Run-II-Z-HighMass-mu" in parameterName :
process +=insert_ME("MEqq2gZ2ff","Muon")
process+=addLeptonPairCut("150","600")
elif "Run-II-Z-mu" in parameterName :
process +=insert_ME("MEqq2gZ2ff","Muon")
elif(simulation=="Powheg") :
if "Run-I-WZ" in parameterName :
process+=insert_ME("PowhegMEqq2W2ff","Electron")
process+=insert_ME("PowhegMEqq2gZ2ff","Electron")
elif "Run-II-W" in parameterName or "Run-I-W" in parameterName :
process+=insert_ME("PowhegMEqq2W2ff","Electron")
elif "Run-II-Z-e" in parameterName or "Run-I-Z" in parameterName :
process+=insert_ME("PowhegMEqq2gZ2ff","Electron")
elif "Run-II-Z-LowMass-mu" in parameterName :
process+=insert_ME("PowhegMEqq2gZ2ff","Muon")
process+=addLeptonPairCut("25","70")
elif "Run-II-Z-HighMass-mu" in parameterName :
process+=insert_ME("PowhegMEqq2gZ2ff","Muon")
process+=addLeptonPairCut("150","600")
elif "Run-II-Z-mu" in parameterName :
process+=insert_ME("PowhegMEqq2gZ2ff","Muon")
elif "DiPhoton-GammaGamma" in parameterName :
process+=insert_ME("MEGammaGammaPowheg","GammaGamma")
process+=insert_ME("MEGammaGamma","gg")
process+="set /Herwig/Cuts/PhotonKtCut:MinKT 5.\n"
process+=jet_kt_cut(5.)
parameterName=parameterName.replace("-GammaGamma","")
elif "DiPhoton-GammaJet" in parameterName :
process+=insert_ME("MEGammaGammaPowheg","VJet")
process+="set /Herwig/Cuts/PhotonKtCut:MinKT 5.\n"
process+=jet_kt_cut(5.)
parameterName=parameterName.replace("-GammaJet","")
elif(simulation=="Matchbox" or simulation=="Merging" ) :
if "Jets" in parameterName :
if(simulation=="Matchbox"):
process+=addProcess(thefactory,"p p -> j j","2","0","MaxJetPtScale",0,0)
elif(simulation=="Merging"):
process+=addProcess(thefactory,"p p -> j j","2","0","MaxJetPtScale",1,0)
if "DiJets" in parameterName :
process+=addFirstJet("30")+addSecondJet("25")
cuts=[100.,300.,600.,900.,ecms]
for i in range(1,len(cuts)) :
tstring = "-DiJets-%s"%i
if tstring in parameterName :
process+=addJetPairCut(cuts[i-1],cuts[i])
parameterName=parameterName.replace(tstring,"-DiJets")
else :
if "Run" in parameterName :
cuts=[5.,20.,40.,80.,160.,320.]
elif "300" in parameterName :
cuts=[5.,]
elif "630" in parameterName :
cuts=[5.,20.,40.]
elif "900" in parameterName :
cuts=[5.,]
cuts.append(ecms)
for i in range(1,len(cuts)) :
tstring = "-Jets-%s"%i
if tstring in parameterName :
process+=addFirstJet(cuts[i-1],cuts[i])
parameterName=parameterName.replace(tstring,"-Jets")
elif "Run-I-WZ" in parameterName :
if(simulation=="Matchbox"):
process+=addProcess(thefactory,"p pbar e+ e-","0","2","LeptonPairMassScale",0,0)
process+=addProcess(thefactory,"p pbar e+ nu","0","2","LeptonPairMassScale",0,0)
process+=addProcess(thefactory,"p pbar e- nu","0","2","LeptonPairMassScale",0,0)
elif(simulation=="Merging"):
process+=particlegroup(thefactory,'epm','e+','e-')
process+=particlegroup(thefactory,'epmnu','e+','e-','nu_e','nu_ebar')
process+=addProcess(thefactory,"p pbar epm epmnu","0","2","LeptonPairMassScale",2,2)
process+=addLeptonPairCut("60","120")
elif "Run-II-W" in parameterName or "Run-I-W" in parameterName :
if(simulation=="Matchbox"):
process+=addProcess(thefactory,"p pbar e+ nu","0","2","LeptonPairMassScale",0,0)
process+=addProcess(thefactory,"p pbar e- nu","0","2","LeptonPairMassScale",0,0)
elif(simulation=="Merging"):
process+=particlegroup(thefactory,'epm','e+','e-')
process+=addProcess(thefactory,"p pbar epm nu","0","2","LeptonPairMassScale",2,2)
process+=addLeptonPairCut("60","120")
elif "Run-II-Z-e" in parameterName or "Run-I-Z" in parameterName :
if(simulation=="Matchbox"):
process+=addProcess(thefactory,"p pbar e+ e-","0","2","LeptonPairMassScale",0,0)
elif(simulation=="Merging"):
process+=addProcess(thefactory,"p pbar e+ e-","0","2","LeptonPairMassScale",2,2)
process+=addLeptonPairCut("60","120")
elif "Run-II-Z-LowMass-mu" in parameterName :
if(simulation=="Matchbox"):
process+=addProcess(thefactory,"p pbar mu+ mu-","0","2","LeptonPairMassScale",0,0)
elif(simulation=="Merging"):
process+=addProcess(thefactory,"p pbar mu+ mu-","0","2","LeptonPairMassScale",2,2)
process+=addLeptonPairCut("25","70")
elif "Run-II-Z-HighMass-mu" in parameterName :
if(simulation=="Matchbox"):
process+=addProcess(thefactory,"p pbar mu+ mu-","0","2","LeptonPairMassScale",0,0)
elif(simulation=="Merging"):
process+=addProcess(thefactory,"p pbar mu+ mu-","0","2","LeptonPairMassScale",2,2)
process+=addLeptonPairCut("150","600")
elif "Run-II-Z-mu" in parameterName :
if(simulation=="Matchbox"):
process+=addProcess(thefactory,"p pbar mu+ mu-","0","2","LeptonPairMassScale",0,0)
elif(simulation=="Merging"):
process+=addProcess(thefactory,"p pbar mu+ mu-","0","2","LeptonPairMassScale",2,2)
process+=addLeptonPairCut("60","120")
# Star
elif(collider=="Star" ) :
process = StringBuilder("set /Herwig/Decays/DecayHandler:LifeTimeOption 0\n")
process+= "set /Herwig/Decays/DecayHandler:MaxLifeTime 10*mm\n"
process+= "set /Herwig/Generators/EventGenerator:EventHandler:BeamB /Herwig/Particles/p+\n"
process+= collider_lumi(200.0)
process+= "set /Herwig/Cuts/Cuts:X2Min 0.01\n"
if(simulation=="") :
if "UE" in parameterName :
if "Dipole" in parameters["shower"]:
process+="read snippets/MB-DipoleShower.in\n"
else:
process+="read snippets/MB.in\n"
process+="read snippets/Diffraction.in\n"
else :
process+=insert_ME("MEQCD2to2")
process+="set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0\n"
if "Jets-1" in parameterName : process+=jet_kt_cut(2.)
elif "Jets-2" in parameterName : process+=jet_kt_cut(5.)
elif "Jets-3" in parameterName : process+=jet_kt_cut(20.)
elif "Jets-4" in parameterName : process+=jet_kt_cut(25.)
else :
logging.error("Star not supported for %s " % simulation)
sys.exit(1)
# ISR and SppS
elif ( collider=="ISR" or collider =="SppS" or collider == "SPS" or collider == "Fermilab" ) :
process = StringBuilder("set /Herwig/Decays/DecayHandler:LifeTimeOption 0\n")
process+="set /Herwig/Decays/DecayHandler:MaxLifeTime 10*mm\n"
if(collider=="SppS") :
process = StringBuilder("set /Herwig/Generators/EventGenerator:EventHandler:BeamB /Herwig/Particles/pbar-\n")
if "17.4" in parameterName : process+=collider_lumi( 17.4)
elif "27.4" in parameterName : process+=collider_lumi( 27.4)
elif "30" in parameterName : process+=collider_lumi( 30.4)
elif "38.8" in parameterName : process+=collider_lumi( 38.8)
elif "44" in parameterName : process+=collider_lumi( 44.4)
elif "53" in parameterName : process+=collider_lumi( 53.0)
elif "62" in parameterName : process+=collider_lumi( 62.2)
elif "63" in parameterName : process+=collider_lumi( 63.0)
elif "200" in parameterName : process+=collider_lumi(200.0)
elif "500" in parameterName : process+=collider_lumi(500.0)
elif "546" in parameterName : process+=collider_lumi(546.0)
elif "900" in parameterName : process+=collider_lumi(900.0)
if "UE" in parameterName :
if(simulation=="") :
if "Dipole" in parameters["shower"]:
process+="read snippets/MB-DipoleShower.in\n"
else:
process+="read snippets/MB.in\n"
+ if "EHS" not in name :
process+="read snippets/Diffraction.in\n"
else :
logging.error(" SppS and ISR not supported for %s " % simulation)
sys.exit(1)
elif "Z-mu" in parameterName :
if simulation == "" :
process+=insert_ME("MEqq2gZ2ff","Muon")
process+=mhat_minm_maxm(2,2,20)
elif simulation == "Powheg" :
process+=insert_ME("PowhegMEqq2gZ2ff","Muon")
process+=mhat_minm_maxm(2,2,20)
elif(simulation=="Matchbox"):
process+=addProcess(thefactory,"p p mu+ mu-","0","2","LeptonPairMassScale",0,0)
process+=addLeptonPairCut("2","20")
elif(simulation=="Merging"):
process+=addProcess(thefactory,"p p mu+ mu-","0","2","LeptonPairMassScale",2,2)
process+=addLeptonPairCut("2","20")
else :
logging.error(" SppS and ISR not supported for %s " % simulation)
sys.exit(1)
else :
logging.error(" Process not supported for SppS and ISR %s " % parameterName )
sys.exit(1)
# LHC
elif(collider=="LHC") :
ecms=7000.0
if parameterName.startswith("7-") : ecms = 7000.0
elif parameterName.startswith("8-") : ecms = 8000.0
elif parameterName.startswith("13-") : ecms = 13000.0
elif parameterName.startswith("900") : ecms = 900.0
elif parameterName.startswith("2360") : ecms = 2360.0
elif parameterName.startswith("2760") : ecms = 2760.0
elif parameterName.startswith("5-") : ecms = 5020.0
else : ecms = 7000.0
process = StringBuilder(collider_lumi(ecms))
if(simulation=="") :
if "VBF" in parameterName :
process+=insert_ME("MEPP2HiggsVBF")
if "GammaGamma" in parameterName :
process+=selectDecayMode("h0",["h0->gamma,gamma;"])
addedBRReweighter = True
elif "WW" in parameterName :
process+=selectDecayMode("h0",["h0->W+,W-;"])
addedBRReweighter = True
elif "ZZ" in parameterName :
process+=selectDecayMode("h0",["h0->Z0,Z0;"])
addedBRReweighter = True
elif "8-" not in parameterName :
process+=selectDecayMode("h0",["h0->tau-,tau+;"])
addedBRReweighter = True
process+="set /Herwig/Particles/tau-:Stable Stable\n"
elif "ggHJet" in parameterName :
process+=selectDecayMode("h0",["h0->tau-,tau+;"])
addedBRReweighter = True
process+="set /Herwig/Particles/tau-:Stable Stable\n"
process+=insert_ME("MEHiggsJet")
process+=jet_kt_cut(20.)
elif "ggH" in parameterName :
process+=insert_ME("MEHiggs")
process+=insert_ME("MEHiggsJet","qqbar")
process+=jet_kt_cut(0.0)
if "GammaGamma" in parameterName :
process+=selectDecayMode("h0",["h0->gamma,gamma;"])
addedBRReweighter = True
elif "WW" in parameterName :
process+=selectDecayMode("h0",["h0->W+,W-;"])
addedBRReweighter = True
elif "ZZ" in parameterName :
process+=selectDecayMode("h0",["h0->Z0,Z0;"])
addedBRReweighter = True
elif "8-" not in parameterName :
process+=selectDecayMode("h0",["h0->tau-,tau+;"])
addedBRReweighter = True
process+="set /Herwig/Particles/tau-:Stable Stable\n"
elif "PromptPhoton" in parameterName :
process+=insert_ME("MEGammaJet")
if "PromptPhoton-1" in parameterName :
process+="set /Herwig/Cuts/PhotonKtCut:MinKT 5.\n"
process+="set /Herwig/Cuts/PhotonKtCut:MaxKT 25.\n"
parameterName=parameterName.replace("-1","")
elif "PromptPhoton-2" in parameterName :
process+="set /Herwig/Cuts/PhotonKtCut:MinKT 25.\n"
process+="set /Herwig/Cuts/PhotonKtCut:MaxKT 80.\n"
parameterName=parameterName.replace("-2","")
elif "PromptPhoton-3" in parameterName :
process+="set /Herwig/Cuts/PhotonKtCut:MinKT 80.\n"
process+="set /Herwig/Cuts/PhotonKtCut:MaxKT 150.\n"
parameterName=parameterName.replace("-3","")
elif "PromptPhoton-4" in parameterName :
process+="set /Herwig/Cuts/PhotonKtCut:MinKT 150.\n"
process+="set /Herwig/Cuts/PhotonKtCut:MaxKT 500.\n"
parameterName=parameterName.replace("-4","")
elif "PromptPhoton-5" in parameterName :
process+="set /Herwig/Cuts/PhotonKtCut:MinKT 500.\n"
parameterName=parameterName.replace("-5","")
elif "DiPhoton-GammaGamma" in parameterName :
process+=insert_ME("MEGammaGamma")
process+="set /Herwig/Cuts/PhotonKtCut:MinKT 5.\n"
parameterName=parameterName.replace("-GammaGamma","")
elif "DiPhoton-GammaJet" in parameterName :
process+=insert_ME("MEGammaJet")
process+="set /Herwig/Cuts/PhotonKtCut:MinKT 5.\n"
parameterName=parameterName.replace("-GammaJet","")
elif "8-WH" in parameterName :
process+=insert_ME("MEPP2WH")
process+=jet_kt_cut(0.0)
if "GammaGamma" in parameterName :
process+=selectDecayMode("h0",["h0->gamma,gamma;"])
addedBRReweighter = True
elif "WW" in parameterName :
process+=selectDecayMode("h0",["h0->W+,W-;"])
addedBRReweighter = True
elif "ZZ" in parameterName :
process+=selectDecayMode("h0",["h0->Z0,Z0;"])
addedBRReweighter = True
elif "8-ZH" in parameterName :
process+=insert_ME("MEPP2ZH")
process+=jet_kt_cut(0.0)
if "GammaGamma" in parameterName :
process+=selectDecayMode("h0",["h0->gamma,gamma;"])
addedBRReweighter = True
elif "WW" in parameterName :
process+=selectDecayMode("h0",["h0->W+,W-;"])
addedBRReweighter = True
elif "ZZ" in parameterName :
process+=selectDecayMode("h0",["h0->Z0,Z0;"])
addedBRReweighter = True
elif "WH" in parameterName :
process+=selectDecayMode("h0",["h0->b,bbar;"])
process+=selectDecayMode("W+",["W+->nu_e,e+;",
"W+->nu_mu,mu+;"])
addedBRReweighter = True
process+=insert_ME("MEPP2WH")
process+=jet_kt_cut(0.0)
elif "ZH" in parameterName :
process+=selectDecayMode("h0",["h0->b,bbar;"])
process+=selectDecayMode("Z0",["Z0->e-,e+;",
"Z0->mu-,mu+;"])
addedBRReweighter = True
process+=insert_ME("MEPP2ZH")
process+=jet_kt_cut(0.0)
elif "UE" in parameterName or "Cent" in parameterName :
if "Dipole" in parameters["shower"]:
process+="read snippets/MB-DipoleShower.in\n"
else:
process+="set /Herwig/Shower/ShowerHandler:IntrinsicPtGaussian 2.2*GeV\n"
process+="read snippets/MB.in\n"
process+="read snippets/Diffraction.in\n"
if "Long" in parameterName :
process += "set /Herwig/Decays/DecayHandler:MaxLifeTime 100*mm\n"
elif "8-DiJets" in parameterName or "7-DiJets" in parameterName or "13-DiJets" in parameterName :
process+=insert_ME("MEQCD2to2")
process+="set MEQCD2to2:MaximumFlavour 5\n"
process+="set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0\n"
if "13-DiJets" not in parameterName :
if "-A" in parameterName :
process+=jet_kt_cut(45.)
process+="set /Herwig/Cuts/JetKtCut:MinEta -3.\n"
process+="set /Herwig/Cuts/JetKtCut:MaxEta 3.\n"
elif "-B" in parameterName :
process+=jet_kt_cut(20.)
process+="set /Herwig/Cuts/JetKtCut:MinEta -2.7\n"
process+="set /Herwig/Cuts/JetKtCut:MaxEta 2.7\n"
elif "-C" in parameterName :
process+=jet_kt_cut(20.)
process+="set /Herwig/Cuts/JetKtCut:MinEta -4.8\n"
process+="set /Herwig/Cuts/JetKtCut:MaxEta 4.8\n"
else :
if "-A" in parameterName :
process+=jet_kt_cut(60.)
process+="set /Herwig/Cuts/JetKtCut:MinEta -3.\n"
process+="set /Herwig/Cuts/JetKtCut:MaxEta 3.\n"
elif "-B" in parameterName :
process+=jet_kt_cut(180.)
process+="set /Herwig/Cuts/JetKtCut:MinEta -3.\n"
process+="set /Herwig/Cuts/JetKtCut:MaxEta 3.\n"
if "DiJets-1" in parameterName : process+=mhat_cut(90.)
elif "DiJets-2" in parameterName : process+=mhat_cut(200.)
elif "DiJets-3" in parameterName : process+=mhat_cut(450.)
elif "DiJets-4" in parameterName : process+=mhat_cut(750.)
elif "DiJets-5" in parameterName : process+=mhat_cut(950.)
elif "DiJets-6" in parameterName : process+=mhat_cut(1550.)
elif "DiJets-7" in parameterName : process+=mhat_cut(2150.)
elif "DiJets-8" in parameterName : process+=mhat_cut(2750.)
elif "DiJets-9" in parameterName : process+=mhat_cut(3750.)
elif "DiJets-10" in parameterName : process+=mhat_cut(4750.)
elif "DiJets-11" in parameterName : process+=mhat_cut(5750.)
elif( "7-Jets" in parameterName
or "8-Jets" in parameterName
or "13-Jets" in parameterName
or "2760-Jets" in parameterName
) :
process+=insert_ME("MEQCD2to2")
process+="set MEQCD2to2:MaximumFlavour 5\n"
process+="set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0\n"
if "Jets-10" in parameterName : process+=jet_kt_cut(1800.)
elif "Jets-0" in parameterName : process+=jet_kt_cut(5.)
elif "Jets-1" in parameterName : process+=jet_kt_cut(10.)
elif "Jets-2" in parameterName : process+=jet_kt_cut(20.)
elif "Jets-3" in parameterName : process+=jet_kt_cut(40.)
elif "Jets-4" in parameterName : process+=jet_kt_cut(70.)
elif "Jets-5" in parameterName : process+=jet_kt_cut(150.)
elif "Jets-6" in parameterName : process+=jet_kt_cut(200.)
elif "Jets-7" in parameterName : process+=jet_kt_cut(300.)
elif "Jets-8" in parameterName : process+=jet_kt_cut(500.)
elif "Jets-9" in parameterName : process+=jet_kt_cut(800.)
elif( "-Charm" in parameterName or "-Bottom" in parameterName ) :
if("8-Bottom" in parameterName) :
addBRReweighter()
process+=selectDecayMode("Jpsi",["Jpsi->mu-,mu+;"])
if "Bottom" in parameterName :
process+="cp MEHeavyQuark MEBottom\n"
process+="set MEBottom:QuarkType Bottom\n"
process+=insert_ME("MEBottom")
else :
process+="cp MEHeavyQuark MECharm\n"
process+="set MECharm:QuarkType Charm\n"
process+=insert_ME("MECharm")
process+="set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0\n"
if "-0" in parameterName :
if "Bottom" in parameterName :
process+="set MEBottom:Process Pair\n"
process+=jet_kt_cut(0.)
else :
process+=jet_kt_cut(1.)
elif "-1" in parameterName : process+=jet_kt_cut(5.)
elif "-2" in parameterName : process+=jet_kt_cut(15.)
elif "-3" in parameterName : process+=jet_kt_cut(20.)
elif "-4" in parameterName : process+=jet_kt_cut(50.)
elif "-5" in parameterName : process+=jet_kt_cut(80.)
elif "-6" in parameterName : process+=jet_kt_cut(110.)
elif "-7" in parameterName : process+=jet_kt_cut(30.)+mhat_cut(90.)
elif "-8" in parameterName : process+=jet_kt_cut(30.)+mhat_cut(340.)
elif "-9" in parameterName : process+=jet_kt_cut(30.)+mhat_cut(500.)
elif "Top-L" in parameterName :
process+="set MEHeavyQuark:QuarkType Top\n"
process+=insert_ME("MEHeavyQuark")
process+=selectDecayMode("t",["t->nu_e,e+,b;",
"t->nu_mu,mu+,b;"])
process+=addBRReweighter()
elif "Top-SL" in parameterName :
process+="set MEHeavyQuark:QuarkType Top\n"
process+=insert_ME("MEHeavyQuark")
process+="set /Herwig/Particles/t:Synchronized Not_synchronized\n"
process+="set /Herwig/Particles/tbar:Synchronized Not_synchronized\n"
process+=selectDecayMode("t",["t->nu_e,e+,b;","t->nu_mu,mu+,b;"])
process+=selectDecayMode("tbar",["tbar->b,bbar,cbar;",
"tbar->bbar,cbar,d;",
"tbar->bbar,cbar,s;",
"tbar->bbar,s,ubar;",
"tbar->bbar,ubar,d;"])
process+=addBRReweighter()
elif "Top-All" in parameterName :
process+="set MEHeavyQuark:QuarkType Top\n"
process+=insert_ME("MEHeavyQuark")
elif "WZ" in parameterName :
process+=insert_ME("MEPP2VV","WZ")
process+=selectDecayMode("W+",["W+->nu_e,e+;",
"W+->nu_mu,mu+;"])
process+=selectDecayMode("W-",["W-->nu_ebar,e-;",
"W-->nu_mubar,mu-;"])
process+=selectDecayMode("Z0",["Z0->e-,e+;",
"Z0->mu-,mu+;"])
addedBRReweighter = True
elif "WW-emu" in parameterName :
process+=insert_ME("MEPP2VV","WW")
process+="set /Herwig/Particles/W+:Synchronized 0\n"
process+="set /Herwig/Particles/W-:Synchronized 0\n"
process+=selectDecayMode("W+",["W+->nu_e,e+;"])
process+=selectDecayMode("W-",["W-->nu_mubar,mu-;"])
addedBRReweighter = True
elif "WW-ll" in parameterName :
process+=insert_ME("MEPP2VV","WW")
process+=selectDecayMode("W+",["W+->nu_e,e+;","W+->nu_mu,mu+;","W+->nu_tau,tau+;"])
addedBRReweighter = True
elif "ZZ-ll" in parameterName :
process+=insert_ME("MEPP2VV","ZZ")
process+=selectDecayMode("Z0",["Z0->e-,e+;",
"Z0->mu-,mu+;",
"Z0->tau-,tau+;"])
addedBRReweighter = True
elif "ZZ-lv" in parameterName :
process+=insert_ME("MEPP2VV","ZZ")
process+=selectDecayMode("Z0",["Z0->e-,e+;",
"Z0->mu-,mu+;",
"Z0->tau-,tau+;",
"Z0->nu_e,nu_ebar;",
"Z0->nu_mu,nu_mubar;",
"Z0->nu_tau,nu_taubar;"])
addedBRReweighter = True
elif "W-e" in parameterName :
process+=insert_ME("MEqq2W2ff","Electron")
elif "W-mu" in parameterName :
process+=insert_ME("MEqq2W2ff","Muon")
elif "Z-e" in parameterName or "Z-mu" in parameterName :
if "Z-e" in parameterName:
process+=insert_ME("MEqq2gZ2ff","Electron")
else :
process+=insert_ME("MEqq2gZ2ff","Muon")
mcuts=[10,35,75,110,400,ecms]
for i in range(1,6) :
tstring = "-Mass%s"%i
if tstring in parameterName :
process+=mhat_minm_maxm(mcuts[i-1],mcuts[i-1],mcuts[i])
parameterName=parameterName.replace(tstring,"")
elif "Z-nu" in parameterName :
process+=insert_ME("MEqq2gZ2ff","Neutrinos")
elif "W-Jet" in parameterName :
process+=insert_ME("MEWJet","Electron","WDecay")
if "W-Jet-1-e" in parameterName :
process+="set /Herwig/Cuts/WBosonKtCut:MinKT 100.0*GeV\n"
parameterName=parameterName.replace("W-Jet-1-e","W-Jet-e")
elif "W-Jet-2-e" in parameterName :
process+="set /Herwig/Cuts/WBosonKtCut:MinKT 190.0*GeV\n"
parameterName=parameterName.replace("W-Jet-2-e","W-Jet-e")
elif "W-Jet-3-e" in parameterName :
process+="set /Herwig/Cuts/WBosonKtCut:MinKT 270.0*GeV\n"
parameterName=parameterName.replace("W-Jet-3-e","W-Jet-e")
elif "Z-Jet" in parameterName :
if "-e" in parameterName :
process+=insert_ME("MEZJet","Electron","ZDecay")
if "Z-Jet-0-e" in parameterName :
process+="set /Herwig/Cuts/ZBosonKtCut:MinKT 35.0*GeV\n"
parameterName=parameterName.replace("Z-Jet-0-e","Z-Jet-e")
elif "Z-Jet-1-e" in parameterName :
process+="set /Herwig/Cuts/ZBosonKtCut:MinKT 100.0*GeV\n"
parameterName=parameterName.replace("Z-Jet-1-e","Z-Jet-e")
elif "Z-Jet-2-e" in parameterName :
process+="set /Herwig/Cuts/ZBosonKtCut:MinKT 190.0*GeV\n"
parameterName=parameterName.replace("Z-Jet-2-e","Z-Jet-e")
elif "Z-Jet-3-e" in parameterName :
process+="set /Herwig/Cuts/ZBosonKtCut:MinKT 270.0*GeV\n"
parameterName=parameterName.replace("Z-Jet-3-e","Z-Jet-e")
else :
process+=insert_ME("MEZJet","Muon","ZDecay")
process+="set /Herwig/Cuts/ZBosonKtCut:MinKT 35.0*GeV\n"
parameterName=parameterName.replace("Z-Jet-0-mu","Z-Jet-mu")
elif "WGamma" in parameterName :
process+=insert_ME("MEPP2VGamma","1")
process+="set MEPP2VGamma:MassOption 1"
process+="set /Herwig/Cuts/PhotonKtCut:MinKT 10.\n"
if "-e" in parameterName :
process+=selectDecayMode("W+",["W+->nu_e,e+;"])
addedBRReweighter=True
else :
process+=selectDecayMode("W+",["W+->nu_mu,mu+;"])
addedBRReweighter=True
elif "ZGamma" in parameterName :
process+=insert_ME("MEPP2VGamma","2")
process+="set /Herwig/Cuts/PhotonKtCut:MinKT 10.\n"
if "-e" in parameterName :
process+=selectDecayMode("Z0",["Z0->e-,e+;"])
addedBRReweighter=True
elif "-mu" in parameterName :
process+=selectDecayMode("Z0",["Z0->mu-,mu+;"])
addedBRReweighter=True
elif "-nu" in parameterName :
process+=selectDecayMode("Z0",["Z0->nu_e,nu_ebar;","Z0->nu_mu,nu_mubar;","Z0->nu_tau,nu_taubar;"])
addedBRReweighter=True
else :
logging.error(" Process %s not supported for internal matrix elements" % name)
sys.exit(1)
elif(simulation=="Powheg") :
if "VBF" in parameterName :
process+=insert_ME("PowhegMEPP2HiggsVBF")
if "GammaGamma" in parameterName :
process+=selectDecayMode("h0",["h0->gamma,gamma;"])
addedBRReweighter = True
elif "WW" in parameterName :
process+=selectDecayMode("h0",["h0->W+,W-;"])
addedBRReweighter = True
elif "ZZ" in parameterName :
process+=selectDecayMode("h0",["h0->Z0,Z0;"])
addedBRReweighter = True
elif "8-" not in parameterName :
process+=selectDecayMode("h0",["h0->tau-,tau+;"])
addedBRReweighter = True
process+="set /Herwig/Particles/tau-:Stable Stable\n"
elif "ggHJet" in parameterName :
logging.error(" Process %s not supported for POWHEG matrix elements" % name)
sys.exit(1)
elif "ggH" in parameterName :
process+=insert_ME("PowhegMEHiggs")
if "GammaGamma" in parameterName :
process+=selectDecayMode("h0",["h0->gamma,gamma;"])
addedBRReweighter = True
elif "WW" in parameterName :
process+=selectDecayMode("h0",["h0->W+,W-;"])
addedBRReweighter = True
elif "ZZ" in parameterName :
process+=selectDecayMode("h0",["h0->Z0,Z0;"])
addedBRReweighter = True
elif "8-" not in parameterName :
process+=selectDecayMode("h0",["h0->tau-,tau+;"])
addedBRReweighter = True
process+="set /Herwig/Particles/tau-:Stable Stable\n"
elif "8-WH" in parameterName :
process+=insert_ME("PowhegMEPP2WH")
process+=jet_kt_cut(0.0)
if "GammaGamma" in parameterName :
process+=selectDecayMode("h0",["h0->gamma,gamma;"])
addedBRReweighter = True
elif "WW" in parameterName :
process+=selectDecayMode("h0",["h0->W+,W-;"])
addedBRReweighter = True
elif "ZZ" in parameterName :
process+=selectDecayMode("h0",["h0->Z0,Z0;"])
addedBRReweighter = True
elif "8-ZH" in parameterName :
process+=insert_ME("PowhegMEPP2ZH")
process+=jet_kt_cut(0.0)
if "GammaGamma" in parameterName :
process+=selectDecayMode("h0",["h0->gamma,gamma;"])
addedBRReweighter = True
elif "WW" in parameterName :
process+=selectDecayMode("h0",["h0->W+,W-;"])
addedBRReweighter = True
elif "ZZ" in parameterName :
process+=selectDecayMode("h0",["h0->Z0,Z0;"])
addedBRReweighter = True
elif "WH" in parameterName :
process+=selectDecayMode("h0",["h0->b,bbar;"])
process+=selectDecayMode("W+",["W+->nu_e,e+;",
"W+->nu_mu,mu+;"])
addedBRReweighter = True
process+=insert_ME("PowhegMEPP2WH")
process+=jet_kt_cut(0.0)
elif "ZH" in parameterName :
process+=selectDecayMode("h0",["h0->b,bbar;"])
process+=selectDecayMode("Z0",["Z0->e-,e+;",
"Z0->mu-,mu+;"])
addedBRReweighter = True
process+=insert_ME("PowhegMEPP2ZH")
process+=jet_kt_cut(0.0)
elif "UE" in parameterName :
logging.error(" Process %s not supported for powheg matrix elements" % name)
sys.exit(1)
elif "WZ" in parameterName :
process+="create Herwig::HwDecayHandler /Herwig/NewPhysics/DecayHandler\n"
process+="set /Herwig/NewPhysics/DecayHandler:NewStep No\n"
process+="set /Herwig/Shower/ShowerHandler:SplitHardProcess No\n";
process+="set /Herwig/Decays/ZDecayer:PhotonGenerator NULL\n";
process+="set /Herwig/Decays/WDecayer:PhotonGenerator NULL\n";
process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 0 /Herwig/Particles/tau-\n"
process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 1 /Herwig/Particles/tau+\n"
process+="insert /Herwig/Generators/EventGenerator:EventHandler:PreCascadeHandlers 0 /Herwig/NewPhysics/DecayHandler\n"
process+=insert_ME("PowhegMEPP2VV","WZ")
process+=selectDecayMode("W+",["W+->nu_e,e+;",
"W+->nu_mu,mu+;"])
process+=selectDecayMode("W-",["W-->nu_ebar,e-;",
"W-->nu_mubar,mu-;"])
process+=selectDecayMode("Z0",["Z0->e-,e+;",
"Z0->mu-,mu+;"])
addedBRReweighter = True
elif "WW-emu" in parameterName :
process+="create Herwig::HwDecayHandler /Herwig/NewPhysics/DecayHandler\n"
process+="set /Herwig/NewPhysics/DecayHandler:NewStep No\n"
process+="set /Herwig/Shower/ShowerHandler:SplitHardProcess No\n";
process+="set /Herwig/Decays/ZDecayer:PhotonGenerator NULL\n";
process+="set /Herwig/Decays/WDecayer:PhotonGenerator NULL\n";
process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 0 /Herwig/Particles/tau-\n"
process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 1 /Herwig/Particles/tau+\n"
process+="insert /Herwig/Generators/EventGenerator:EventHandler:PreCascadeHandlers 0 /Herwig/NewPhysics/DecayHandler\n"
process+=insert_ME("PowhegMEPP2VV","WW")
process+="set /Herwig/Particles/W+:Synchronized 0\n"
process+="set /Herwig/Particles/W-:Synchronized 0\n"
process+=selectDecayMode("W+",["W+->nu_e,e+;"])
process+=selectDecayMode("W-",["W-->nu_mubar,mu-;"])
addedBRReweighter = True
elif "WW-ll" in parameterName :
process+="create Herwig::HwDecayHandler /Herwig/NewPhysics/DecayHandler\n"
process+="set /Herwig/NewPhysics/DecayHandler:NewStep No\n"
process+="set /Herwig/Shower/ShowerHandler:SplitHardProcess No\n";
process+="set /Herwig/Decays/ZDecayer:PhotonGenerator NULL\n";
process+="set /Herwig/Decays/WDecayer:PhotonGenerator NULL\n";
process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 0 /Herwig/Particles/tau-\n"
process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 1 /Herwig/Particles/tau+\n"
process+="insert /Herwig/Generators/EventGenerator:EventHandler:PreCascadeHandlers 0 /Herwig/NewPhysics/DecayHandler\n"
process+=insert_ME("PowhegMEPP2VV","WW")
process+=selectDecayMode("W+",["W+->nu_e,e+;",
"W+->nu_mu,mu+;",
"W+->nu_tau,tau+;"])
addedBRReweighter = True
elif "ZZ-ll" in parameterName :
process+="create Herwig::HwDecayHandler /Herwig/NewPhysics/DecayHandler\n"
process+="set /Herwig/NewPhysics/DecayHandler:NewStep No\n"
process+="set /Herwig/Shower/ShowerHandler:SplitHardProcess No\n";
process+="set /Herwig/Decays/ZDecayer:PhotonGenerator NULL\n";
process+="set /Herwig/Decays/WDecayer:PhotonGenerator NULL\n";
process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 0 /Herwig/Particles/tau-\n"
process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 1 /Herwig/Particles/tau+\n"
process+="insert /Herwig/Generators/EventGenerator:EventHandler:PreCascadeHandlers 0 /Herwig/NewPhysics/DecayHandler\n"
process+=insert_ME("PowhegMEPP2VV","ZZ")
process+=selectDecayMode("Z0",["Z0->e-,e+;",
"Z0->mu-,mu+;",
"Z0->tau-,tau+;"])
addedBRReweighter = True
elif "ZZ-lv" in parameterName :
process+="create Herwig::HwDecayHandler /Herwig/NewPhysics/DecayHandler\n"
process+="set /Herwig/NewPhysics/DecayHandler:NewStep No\n"
process+="set /Herwig/Shower/ShowerHandler:SplitHardProcess No\n";
process+="set /Herwig/Decays/ZDecayer:PhotonGenerator NULL\n";
process+="set /Herwig/Decays/WDecayer:PhotonGenerator NULL\n";
process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 0 /Herwig/Particles/tau-\n"
process+="insert /Herwig/NewPhysics/DecayHandler:Excluded 1 /Herwig/Particles/tau+\n"
process+="insert /Herwig/Generators/EventGenerator:EventHandler:PreCascadeHandlers 0 /Herwig/NewPhysics/DecayHandler\n"
process+=insert_ME("PowhegMEPP2VV","ZZ")
process+=selectDecayMode("Z0",["Z0->e-,e+;",
"Z0->mu-,mu+;",
"Z0->tau-,tau+;",
"Z0->nu_e,nu_ebar;",
"Z0->nu_mu,nu_mubar;",
"Z0->nu_tau,nu_taubar;"])
addedBRReweighter = True
elif "W-e" in parameterName :
process+=insert_ME("PowhegMEqq2W2ff","Electron")
elif "W-mu" in parameterName :
process+=insert_ME("PowhegMEqq2W2ff","Muon")
elif "Z-e" in parameterName or "Z-mu" in parameterName :
if "Z-e" in parameterName:
process+=insert_ME("PowhegMEqq2gZ2ff","Electron")
else :
process+=insert_ME("PowhegMEqq2gZ2ff","Muon")
mcuts=[10,35,75,110,400,ecms]
for i in range(1,6) :
tstring = "-Mass%s"%i
if tstring in parameterName :
process+=mhat_minm_maxm(mcuts[i-1],mcuts[i-1],mcuts[i])
parameterName=parameterName.replace(tstring,"")
elif "Z-nu" in parameterName :
process+=insert_ME("PowhegMEqq2gZ2ff","Neutrinos")
elif "DiPhoton-GammaGamma" in parameterName :
process+=insert_ME("MEGammaGammaPowheg","GammaGamma")
process+=insert_ME("MEGammaGamma","gg")
process+="set /Herwig/Cuts/PhotonKtCut:MinKT 5.\n"
process+=jet_kt_cut(5.)
parameterName=parameterName.replace("-GammaGamma","")
elif "DiPhoton-GammaJet" in parameterName :
process+=insert_ME("MEGammaGammaPowheg","VJet")
process+="set /Herwig/Cuts/PhotonKtCut:MinKT 5.\n"
process+=jet_kt_cut(5.)
parameterName=parameterName.replace("-GammaJet","")
else :
logging.error(" Process %s not supported for internal POWHEG matrix elements" % name)
sys.exit(1)
elif( simulation=="Matchbox" or simulation=="Merging" ) :
if "VBF" in parameterName :
parameters["nlo"] = "read Matchbox/VBFNLO.in\n"
if(simulation=="Merging"):
process+="cd /Herwig/Merging/\n"
process+="insert "+thefactory+":DiagramGenerator:RestrictLines 0 /Herwig/Particles/Z0\n"
process+="insert "+thefactory+":DiagramGenerator:RestrictLines 0 /Herwig/Particles/W+\n"
process+="insert "+thefactory+":DiagramGenerator:RestrictLines 0 /Herwig/Particles/W-\n"
process+="insert "+thefactory+":DiagramGenerator:RestrictLines 0 /Herwig/Particles/gamma\n"
process+="do "+thefactory+":DiagramGenerator:TimeLikeRange 0 0\n"
if(simulation=="Matchbox"):
process+=addProcess(thefactory,"p p h0 j j","0","3","FixedScale",0,0)
elif(simulation=="Merging"):
process+=addProcess(thefactory,"p p h0 j j","0","3","FixedScale",1,1)
process+=setHardProcessWidthToZero(["h0"])
process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 125.7\n"
if "GammaGamma" in parameterName :
process+=selectDecayMode("h0",["h0->gamma,gamma;"])
process+=addBRReweighter()
elif "WW" in parameterName :
process+=selectDecayMode("h0",["h0->W+,W-;"])
process+=addBRReweighter()
elif "ZZ" in parameterName :
process+=selectDecayMode("h0",["h0->Z0,Z0;"])
process+=addBRReweighter()
elif "8-" not in parameterName :
process+=selectDecayMode("h0",["h0->tau-,tau+;"])
process+=addBRReweighter()
process+="set /Herwig/Particles/tau-:Stable Stable\n"
elif "ggHJet" in parameterName :
if(simulation=="Merging"):
logging.warning("ggHJet not explicitly tested for %s " % simulation)
sys.exit(0)
parameters["nlo"] = "read Matchbox/MadGraph-GoSam.in\nread Matchbox/HiggsEffective.in\n"
process+=selectDecayMode("h0",["h0->tau-,tau+;"])
process+=addBRReweighter()
process+="set /Herwig/Particles/tau-:Stable Stable\n"
process+=setHardProcessWidthToZero(["h0"])
process+=addProcess(thefactory,"p p h0 j","3","1","FixedScale",0,0)
process+=addFirstJet("20")
process+="set "+thefactory+":ScaleChoice /Herwig/MatrixElements/Matchbox/Scales/FixedScale\n"
process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 125.7\n"
elif "ggH" in parameterName :
parameters["nlo"] = "read Matchbox/MadGraph-GoSam.in\nread Matchbox/HiggsEffective.in\n"
if(simulation=="Merging"):
process+= "cd /Herwig/MatrixElements/Matchbox/Amplitudes\nset OpenLoops:HiggsEff Yes\nset MadGraph:Model heft\n"
process+="cd /Herwig/Merging/\n"
process+=setHardProcessWidthToZero(["h0"])
if(simulation=="Matchbox"):
process+=addProcess(thefactory,"p p h0","2","1","FixedScale",0,0)
elif(simulation=="Merging"):
process+=addProcess(thefactory,"p p h0","2","1","FixedScale",2,2)
process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 125.7\n"
if "GammaGamma" in parameterName :
process+=selectDecayMode("h0",["h0->gamma,gamma;"])
process+=addBRReweighter()
elif "WW" in parameterName :
process+=selectDecayMode("h0",["h0->W+,W-;"])
process+=addBRReweighter()
elif "ZZ" in parameterName :
process+=selectDecayMode("h0",["h0->Z0,Z0;"])
process+=addBRReweighter()
elif "8-" not in parameterName :
process+=selectDecayMode("h0",["h0->tau-,tau+;"])
process+=addBRReweighter()
process+="set /Herwig/Particles/tau-:Stable Stable\n"
elif "8-WH" in parameterName :
if(simulation=="Merging"):
logging.warning("8-WH not explicitly tested for %s " % simulation)
sys.exit(0)
process+=setHardProcessWidthToZero(["h0","W+","W-"])
process+=addProcess(thefactory,"p p W+ h0","0","2","FixedScale",0,0)
process+=addProcess(thefactory,"p p W- h0","0","2","FixedScale",0,0)
process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 125.7\n"
if "GammaGamma" in parameterName :
process+=selectDecayMode("h0",["h0->gamma,gamma;"])
process+=addBRReweighter()
elif "WW" in parameterName :
process+=selectDecayMode("h0",["h0->W+,W-;"])
process+=addBRReweighter()
elif "ZZ" in parameterName :
process+=selectDecayMode("h0",["h0->Z0,Z0;"])
process+=addBRReweighter()
elif "8-ZH" in parameterName :
if(simulation=="Merging"):
logging.warning("8-ZH not explicitly tested for %s " % simulation)
sys.exit(0)
process+=setHardProcessWidthToZero(["h0","Z0"])
process+=addProcess(thefactory,"p p Z0 h0","0","2","FixedScale",0,0)
process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 125.7\n"
if "GammaGamma" in parameterName :
process+=selectDecayMode("h0",["h0->gamma,gamma;"])
process+=addBRReweighter()
elif "WW" in parameterName :
process+=selectDecayMode("h0",["h0->W+,W-;"])
process+=addBRReweighter()
elif "ZZ" in parameterName :
process+=selectDecayMode("h0",["h0->Z0,Z0;"])
process+=addBRReweighter()
elif "WH" in parameterName :
if(simulation=="Merging"):
logging.warning("WH not explicitly tested for %s " % simulation)
sys.exit(0)
process+=selectDecayMode("h0",["h0->b,bbar;"])
process+=addBRReweighter()
process+=setHardProcessWidthToZero(["h0"])
process+=addProcess(thefactory,"p p e+ nu h0","0","3","LeptonPairMassScale",0,0)
process+=addProcess(thefactory,"p p e- nu h0","0","3","LeptonPairMassScale",0,0)
process+=addProcess(thefactory,"p p mu+ nu h0","0","3","LeptonPairMassScale",0,0)
process+=addProcess(thefactory,"p p mu- nu h0","0","3","LeptonPairMassScale",0,0)
process+=addLeptonPairCut("60","120")
elif "ZH" in parameterName :
if(simulation=="Merging"):
logging.warning("ZH not explicitly tested for %s " % simulation)
sys.exit(0)
process+=selectDecayMode("h0",["h0->b,bbar;"])
process+=addBRReweighter()
process+=setHardProcessWidthToZero(["h0"])
process+=addProcess(thefactory,"p p e+ e- h0","0","3","LeptonPairMassScale",0,0)
process+=addProcess(thefactory,"p p mu+ mu- h0","0","3","LeptonPairMassScale",0,0)
process+=addLeptonPairCut("60","120")
elif "UE" in parameterName :
logging.error(" Process %s not supported for Matchbox matrix elements" % name)
sys.exit(1)
elif "8-DiJets" in parameterName or "7-DiJets" in parameterName or "13-DiJets" in parameterName :
if(simulation=="Matchbox"):
process+=addProcess(thefactory,"p p j j","2","0","MaxJetPtScale",0,0)
elif(simulation=="Merging"):
process+=addProcess(thefactory,"p p j j","2","0","MaxJetPtScale",1,1)
process+="set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0\n"
if "13-DiJets" not in parameterName :
if "-A" in parameterName :
process+=addFirstJet("45")
process+=addSecondJet("25")
process+="set /Herwig/Cuts/FirstJet:YRange -3. 3.\n"
process+="set /Herwig/Cuts/SecondJet:YRange -3. 3.\n"
elif "-B" in parameterName :
process+=addFirstJet("20")
process+=addSecondJet("15")
process+="set /Herwig/Cuts/FirstJet:YRange -2.7 2.7\n"
process+="set /Herwig/Cuts/SecondJet:YRange -2.7 2.7\n"
elif "-C" in parameterName :
process+=addFirstJet("20")
process+=addSecondJet("15")
process+="set /Herwig/Cuts/FirstJet:YRange -4.8 4.8\n"
process+="set /Herwig/Cuts/SecondJet:YRange -4.8 4.8\n"
else :
logging.error("Exit 00001")
sys.exit(1)
else :
if "-A" in parameterName :
process+= addFirstJet("75.")
process+=addSecondJet("60.")
process+="set /Herwig/Cuts/JetKtCut:MinEta -3.\n"
process+="set /Herwig/Cuts/JetKtCut:MaxEta 3.\n"
elif "-B" in parameterName :
process+= addFirstJet("220.")
process+=addSecondJet("180.")
process+="set /Herwig/Cuts/JetKtCut:MinEta -3.\n"
process+="set /Herwig/Cuts/JetKtCut:MaxEta 3.\n"
else :
logging.error("Exit 00001")
sys.exit(1)
if "DiJets-1" in parameterName : process+=addJetPairCut("90")
elif "DiJets-2" in parameterName : process+=addJetPairCut("200")
elif "DiJets-3" in parameterName : process+=addJetPairCut("450")
elif "DiJets-4" in parameterName : process+=addJetPairCut("750")
elif "DiJets-5" in parameterName : process+=addJetPairCut("950")
elif "DiJets-6" in parameterName : process+=addJetPairCut("1550")
elif "DiJets-7" in parameterName : process+=addJetPairCut("2150")
elif "DiJets-8" in parameterName : process+=addJetPairCut("2750")
elif "DiJets-9" in parameterName : process+=mhat_cut(3750.)
elif "DiJets-10" in parameterName : process+=mhat_cut(4750.)
elif "DiJets-11" in parameterName : process+=mhat_cut(5750.)
else :
logging.error("Exit 00002")
sys.exit(1)
elif( "7-Jets" in parameterName
or "8-Jets" in parameterName
or "13-Jets" in parameterName
or "2760-Jets" in parameterName
) :
if(simulation=="Matchbox"):
process+=addProcess(thefactory,"p p j j","2","0","MaxJetPtScale",0,0)
elif(simulation=="Merging"):
process+=addProcess(thefactory,"p p j j","2","0","MaxJetPtScale",1,1)
process+="set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0\n"
if "Jets-10" in parameterName : process+=addFirstJet("1800")
elif "Jets-0" in parameterName : process+=addFirstJet("5")
elif "Jets-1" in parameterName : process+=addFirstJet("10")
elif "Jets-2" in parameterName : process+=addFirstJet("20")
elif "Jets-3" in parameterName : process+=addFirstJet("40")
elif "Jets-4" in parameterName : process+=addFirstJet("70")
elif "Jets-5" in parameterName : process+=addFirstJet("150")
elif "Jets-6" in parameterName : process+=addFirstJet("200")
elif "Jets-7" in parameterName : process+=addFirstJet("300")
elif "Jets-8" in parameterName : process+=addFirstJet("500")
elif "Jets-9" in parameterName : process+=addFirstJet("800")
else :
logging.error("Exit 00003")
sys.exit(1)
elif( "-Charm" in parameterName or "-Bottom" in parameterName) :
parameters["bscheme"]=fourFlavour
process+="set /Herwig/Particles/b:HardProcessMass 4.2*GeV\n"
process+="set /Herwig/Particles/bbar:HardProcessMass 4.2*GeV\n"
if("8-Bottom" in parameterName) :
addBRReweighter()
process+=selectDecayMode("Jpsi",["Jpsi->mu-,mu+;"])
if "Bottom" in parameterName :
if(simulation=="Matchbox"):
process+=addProcess(thefactory,"p p b bbar","2","0","MaxJetPtScale",0,0)
elif(simulation=="Merging"):
process+=addProcess(thefactory,"p p b bbar","2","0","MaxJetPtScale",1,0)
else:
if(simulation=="Matchbox"):
process+=addProcess(thefactory,"p p c cbar","2","0","MaxJetPtScale",0,0)
elif(simulation=="Merging"):
process+=addProcess(thefactory,"p p c cbar","2","0","MaxJetPtScale",1,0)
process+="set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0\n"
if "-0" in parameterName : process+=addFirstJet("0")
elif "-1" in parameterName : process+=addFirstJet("5")
elif "-2" in parameterName : process+=addFirstJet("15")
elif "-3" in parameterName : process+=addFirstJet("20")
elif "-4" in parameterName : process+=addFirstJet("50")
elif "-5" in parameterName : process+=addFirstJet("80")
elif "-6" in parameterName : process+=addFirstJet("110")
elif "-7" in parameterName :
process+=addFirstJet("30")
process+=addSecondJet("25")
process+=addJetPairCut("90")
elif "-8" in parameterName :
process+=addFirstJet("30")
process+=addSecondJet("25")
process+=addJetPairCut("340")
elif "-9" in parameterName :
process+=addFirstJet("30")
process+=addSecondJet("25")
process+=addJetPairCut("500")
else :
logging.error("Exit 00004")
sys.exit(1)
elif "Top-L" in parameterName :
process+=setHardProcessWidthToZero(["t","tbar"])
if(simulation=="Matchbox"):
process+=addProcess(thefactory,"p p t tbar","2","0","TopPairMTScale",0,0)
elif(simulation=="Merging"):
process+=addProcess(thefactory,"p p t tbar","2","0","TopPairMTScale",2,2)
process+=selectDecayMode("t",["t->nu_e,e+,b;",
"t->nu_mu,mu+,b;"])
process+=addBRReweighter()
elif "Top-SL" in parameterName :
process+=setHardProcessWidthToZero(["t","tbar"])
if(simulation=="Matchbox"):
process+=addProcess(thefactory,"p p t tbar","2","0","TopPairMTScale",0,0)
elif(simulation=="Merging"):
process+=addProcess(thefactory,"p p t tbar","2","0","TopPairMTScale",2,2)
process+="set /Herwig/Particles/t:Synchronized Not_synchronized\n"
process+="set /Herwig/Particles/tbar:Synchronized Not_synchronized\n"
process+=selectDecayMode("t",["t->nu_e,e+,b;",
"t->nu_mu,mu+,b;"])
process+=selectDecayMode("tbar",["tbar->b,bbar,cbar;",
"tbar->bbar,cbar,d;",
"tbar->bbar,cbar,s;",
"tbar->bbar,s,ubar;",
"tbar->bbar,ubar,d;"])
process+=addBRReweighter()
elif "Top-All" in parameterName :
process+=setHardProcessWidthToZero(["t","tbar"])
if(simulation=="Matchbox"):
process+=addProcess(thefactory,"p p t tbar","2","0","TopPairMTScale",0,0)
elif(simulation=="Merging"):
process+=addProcess(thefactory,"p p t tbar","2","0","TopPairMTScale",2,2)
elif "WZ" in parameterName :
if(simulation=="Merging"):
logging.warning("WZ not explicitly tested for %s " % simulation)
sys.exit(0)
process+=setHardProcessWidthToZero(["W+","W-","Z0"])
process+=addProcess(thefactory,"p p W+ Z0","0","2","FixedScale",0,0)
process+=addProcess(thefactory,"p p W- Z0","0","2","FixedScale",0,0)
process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 171.6*GeV\n\n"
process+=selectDecayMode("W+",["W+->nu_e,e+;",
"W+->nu_mu,mu+;"])
process+=selectDecayMode("W-",["W-->nu_ebar,e-;",
"W-->nu_mubar,mu-;"])
process+=selectDecayMode("Z0",["Z0->e-,e+;",
"Z0->mu-,mu+;"])
process+=addBRReweighter()
process+=addLeptonPairCut("60","120")
elif "WW-emu" in parameterName :
if(simulation=="Merging"):
logging.warning("WW-emu not explicitly tested for %s " % simulation)
sys.exit(0)
process+=setHardProcessWidthToZero(["W+","W-","Z0"])
process+=addProcess(thefactory,"p p W+ W-","0","2","FixedScale",0,0)
process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 160.8*GeV\n"
process+="set /Herwig/Particles/W+:Synchronized 0\n"
process+="set /Herwig/Particles/W-:Synchronized 0\n"
process+=selectDecayMode("W+",["W+->nu_e,e+;"])
process+=selectDecayMode("W-",["W-->nu_mubar,mu-;"])
process+=addBRReweighter()
parameters["bscheme"] = "read Matchbox/FourFlavourScheme.in\n"
process+=addLeptonPairCut("60","120")
elif "WW-ll" in parameterName :
if(simulation=="Merging"):
logging.warning("WW-ll not explicitly tested for %s " % simulation)
sys.exit(0)
process+=setHardProcessWidthToZero(["W+","W-","Z0"])
process+=addProcess(thefactory,"p p W+ W-","0","2","FixedScale",0,0)
process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 160.8*GeV\n"
process+=selectDecayMode("W+",["W+->nu_e,e+;",
"W+->nu_mu,mu+;",
"W+->nu_tau,tau+;"])
process+=addBRReweighter()
process+=addLeptonPairCut("60","120")
parameters["bscheme"] = "read Matchbox/FourFlavourScheme.in\n"
elif "ZZ-ll" in parameterName :
if(simulation=="Merging"):
logging.warning("ZZ-ll not explicitly tested for %s " % simulation)
sys.exit(0)
process+=setHardProcessWidthToZero(["W+","W-","Z0"])
process+=addProcess(thefactory,"p p Z0 Z0","0","2","FixedScale",0,0)
process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 182.2*GeV\n"
process+=selectDecayMode("Z0",["Z0->e-,e+;",
"Z0->mu-,mu+;",
"Z0->tau-,tau+;"])
process+=addBRReweighter()
process+=addLeptonPairCut("60","120")
elif "ZZ-lv" in parameterName :
if(simulation=="Merging"):
logging.warning("ZZ-lv not explicitly tested for %s " % simulation)
sys.exit(0)
process+=setHardProcessWidthToZero(["W+","W-","Z0"])
process+=addProcess(thefactory,"p p Z0 Z0","0","2","FixedScale",0,0)
process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 182.2*GeV\n"
process+=selectDecayMode("Z0",["Z0->e-,e+;",
"Z0->mu-,mu+;",
"Z0->tau-,tau+;",
"Z0->nu_e,nu_ebar;",
"Z0->nu_mu,nu_mubar;",
"Z0->nu_tau,nu_taubar;"])
process+=addBRReweighter()
process+=addLeptonPairCut("60","120")
elif "W-e" in parameterName :
if(simulation=="Matchbox"):
process+=addProcess(thefactory,"p p e+ nu","0","2","LeptonPairMassScale",0,0)
process+=addProcess(thefactory,"p p e- nu","0","2","LeptonPairMassScale",0,0)
elif(simulation=="Merging"):
process+=particlegroup(thefactory,'epm','e+','e-')
process+=addProcess(thefactory,"p p epm nu","0","2","LeptonPairMassScale",2,2)
process+=addLeptonPairCut("60","120")
elif "W-mu" in parameterName :
if(simulation=="Matchbox"):
process+=addProcess(thefactory,"p p mu+ nu","0","2","LeptonPairMassScale",0,0)
process+=addProcess(thefactory,"p p mu- nu","0","2","LeptonPairMassScale",0,0)
elif(simulation=="Merging"):
process+=particlegroup(thefactory,'mupm','mu+','mu-')
process+=addProcess(thefactory,"p p mupm nu","0","2","LeptonPairMassScale",2,2)
process+=addLeptonPairCut("60","120")
elif "Z-e" in parameterName or "Z-mu" in parameterName :
if "Z-e" in parameterName :
if(simulation=="Matchbox"):
process+=addProcess(thefactory,"p p e+ e-","0","2","LeptonPairMassScale",0,0)
elif(simulation=="Merging"):
process+=addProcess(thefactory,"p p e+ e-","0","2","LeptonPairMassScale",2,2)
elif "Z-mu" in parameterName :
if(simulation=="Matchbox"):
process+=addProcess(thefactory,"p p mu+ mu-","0","2","LeptonPairMassScale",0,0)
elif(simulation=="Merging"):
process+=addProcess(thefactory,"p p mu+ mu-","0","2","LeptonPairMassScale",2,2)
mcuts=[10,35,75,110,400,ecms]
for i in range(1,6) :
tstring = "-Mass%s"%i
if tstring in parameterName :
process+=addLeptonPairCut(mcuts[i-1],mcuts[i])
parameterName=parameterName.replace(tstring,"")
elif "Z-nu" in parameterName :
if(simulation=="Matchbox"):
process+=addProcess(thefactory,"p p nu nu","0","2","LeptonPairMassScale",0,0)
elif(simulation=="Merging"):
process+=addProcess(thefactory,"p p nu nu","0","2","LeptonPairMassScale",2,2)
elif "Z-jj" in parameterName :
if(simulation=="Merging"):
logging.warning("Z-jj not explicitly tested for %s " % simulation)
sys.exit(0)
process+=addProcess(thefactory,"p p e+ e- j j","2","2","LeptonPairMassScale",0,0)
process+=addFirstJet("40")
process+=addSecondJet("30")
process+=addLeptonPairCut("60","120")
elif "W-Jet" in parameterName :
if(simulation=="Merging"):
logging.warning("W-Jet not explicitly tested for %s " % simulation)
sys.exit(0)
process+=addProcess(thefactory,"p p e+ nu j","1","2","HTScale",0,0)
process+=addProcess(thefactory,"p p e- nu j","1","2","HTScale",0,0)
process+=addLeptonPairCut("60","120")
if "W-Jet-1-e" in parameterName :
process+=addFirstJet("100")
parameterName=parameterName.replace("W-Jet-1-e","W-Jet-e")
elif "W-Jet-2-e" in parameterName :
process+=addFirstJet("190")
parameterName=parameterName.replace("W-Jet-2-e","W-Jet-e")
elif "W-Jet-3-e" in parameterName :
process+=addFirstJet("270")
parameterName=parameterName.replace("W-Jet-3-e","W-Jet-e")
else :
logging.error("Exit 00005")
sys.exit(1)
elif "Z-Jet" in parameterName :
if(simulation=="Merging"):
logging.warning("Z-Jet not explicitly tested for %s " % simulation)
sys.exit(0)
if "-e" in parameterName :
process+=addProcess(thefactory,"p p e+ e- j","1","2","HTScale",0,0)
if "Z-Jet-0-e" in parameterName :
process+=addFirstJet("35")
parameterName=parameterName.replace("Z-Jet-0-e","Z-Jet-e")
elif "Z-Jet-1-e" in parameterName :
process+=addFirstJet("100")
parameterName=parameterName.replace("Z-Jet-1-e","Z-Jet-e")
elif "Z-Jet-2-e" in parameterName :
process+=addFirstJet("190")
parameterName=parameterName.replace("Z-Jet-2-e","Z-Jet-e")
elif "Z-Jet-3-e" in parameterName :
process+=addFirstJet("270")
parameterName=parameterName.replace("Z-Jet-3-e","Z-Jet-e")
else :
logging.error("Exit 00006")
sys.exit(1)
else :
process+=addProcess(thefactory,"p p mu+ mu- j","1","2","HTScale",0,0)
process+=addFirstJet("35")
parameterName=parameterName.replace("Z-Jet-0-mu","Z-Jet-mu")
process+=addLeptonPairCut("60","120")
elif "Z-bb" in parameterName :
if(simulation=="Merging"):
logging.warning("Z-bb not explicitly tested for %s " % simulation)
sys.exit(0)
parameters["bscheme"]=fourFlavour
process+="set /Herwig/Particles/b:HardProcessMass 4.2*GeV\nset /Herwig/Particles/bbar:HardProcessMass 4.2*GeV\n"
process+=addProcess(thefactory,"p p e+ e- b bbar","2","2","FixedScale",0,0)
process+=addLeptonPairCut("66","116")
process+=addFirstJet("18")
process+=addSecondJet("15")
process+=addLeptonPairCut("60","120")
elif "Z-b" in parameterName :
if(simulation=="Merging"):
logging.warning("Z-b not explicitly tested for %s " % simulation)
sys.exit(0)
process+=particlegroup(thefactory,'bjet','b','bbar')
process+=addProcess(thefactory,"p p e+ e- bjet","1","2","FixedScale",0,0)
process+="set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 91.2*GeV\n"
process+=addLeptonPairCut("60","120")
process+=addFirstJet("15")
elif "W-b" in parameterName :
if(simulation=="Merging"):
logging.warning("W-b not explicitly tested for %s " % simulation)
sys.exit(0)
parameters["bscheme"]=fourFlavour
process += "set /Herwig/Particles/b:HardProcessMass 4.2*GeV\nset /Herwig/Particles/bbar:HardProcessMass 4.2*GeV\n"
process+=addProcess(thefactory,"p p e- nu b bbar","2","2","FixedScale",0,0)
process+=addProcess(thefactory,"p p mu+ nu b bbar","2","2","FixedScale",0,0)
process += "set /Herwig/MatrixElements/Matchbox/Scales/FixedScale:FixedScale 80.4*GeV\n"
process+=addFirstJet("30")
process+=addLeptonPairCut("60","120")
else :
logging.error(" Process %s not supported for Matchbox matrix elements" % name)
sys.exit(1)
# LHC-GammaGamma
elif(collider=="LHC-GammaGamma" ) :
if "7" in parameterName : process = StringBuilder(collider_lumi( 7000.0))
elif "8" in parameterName : process = StringBuilder(collider_lumi( 8000.0))
elif "13" in parameterName : process = StringBuilder(collider_lumi(13000.0))
else : process = StringBuilder(collider_lumi( 7000.0))
if(simulation=="") :
process += "cp MEgg2ff MEgg2ee\n"
process += insert_ME("MEgg2ee","Electron")
process += "cp MEgg2ff MEgg2mm\n"
process += insert_ME("MEgg2mm","Muon")
else :
logging.error("LHC-GammaGamma not supported for %s " % simulation)
sys.exit(1)
+pion=False
+if "Pion-" in parameterName and have_hadronic_collider:
+ parameterName = parameterName.replace("Pion-","")
+ pion = True
+
if "EHS" in name :
pFile = os.path.join(collider,"{c}-{pn}.in".format(c="EHS", pn=parameterName))
elif "Photo" == name[0:5] :
pFile = os.path.join(collider,"{c}-{pn}.in".format(c="Photo", pn=parameterName))
else :
pFile = os.path.join(collider,"{c}-{pn}.in".format(c=collider, pn=parameterName))
with open(os.path.join("Rivet",pFile), 'r') as f:
parameters['parameterFile'] = f.read()
parameters['runname'] = 'Rivet-%s' % name
parameters['process'] = str(process)
if have_hadronic_collider :
if "EHS" in name :
parameters['collider'] = "PPCollider.in\nread snippets/FixedTarget-PP.in"
else :
parameters['collider'] = "PPCollider.in"
+if pion :
+ parameters['collider'] +="\nread snippets/PionPDF.in\nset /Herwig/Generators/EventGenerator:EventHandler:BeamA /Herwig/Particles/pi+\n"
+ parameters['collider'] +="set /Herwig/Shower/ShowerHandler:PDFA /Herwig/Partons/PionPDF\nset /Herwig/Shower/ShowerHandler:PDFARemnant /Herwig/Partons/PionPDF\n"
+ parameters['collider'] +="set /Herwig/Partons/PPExtractor:FirstPDF /Herwig/Partons/PionPDF\nset /Herwig/Partons/MPIExtractor:FirstPDF /Herwig/Partons/PionPDF\n"
+ parameters['collider'] +="set /Herwig/EventHandlers/FixedTargetLuminosity:BeamParticle /Herwig/Particles/pi+\n"
#check if selecteddecaymode and addedBRReweighter is consistent
if selecteddecaymode and not addedBRReweighter:
logging.error("Decaymode was selected but no BRReweighter was added.")
sys.exit(1)
if addedBRReweighter and not selecteddecaymode:
logging.error("BRReweighter was added but no Decaymode was selected.")
sys.exit(1)
# check that we only add one process if in merging mode:
if numberOfAddedProcesses > 1 and simulation =="Merging":
logging.error("In Merging only one process is allowed at the moment. See ticket #403.")
sys.exit(1)
# Check if a process was added for Merging or Matchbox:
if numberOfAddedProcesses == 0 and (simulation =="Merging" or simulation =="Matchbox"):
logging.error("No process was selected.")
sys.exit(1)
# get template and write the file
with open(os.path.join("Rivet/Templates",templateName), 'r') as f:
templateText = f.read()
template = Template( templateText )
with open(os.path.join("Rivet",name+".in"), 'w') as f:
f.write( template.substitute(parameters) )
diff --git a/src/defaults/mesons.in b/src/defaults/mesons.in
--- a/src/defaults/mesons.in
+++ b/src/defaults/mesons.in
@@ -1,893 +1,893 @@
# -*- ThePEG-repository -*-
#
# file containing the particle data for the mesons
#
#
# the 1^1S_0 multiplet
#
-create ThePEG::ParticleData pi+
+create ThePEG::BeamParticleData pi+
setup pi+ 211 pi+ 0.13957039 2.5283740149914E-17 0 7804.5 3 0 1 1
create ThePEG::ParticleData pi0
setup pi0 111 pi0 0.1349768 7.7994841897233E-9 0 2.53E-5 0 0 1 0
-create ThePEG::ParticleData pi-
+create ThePEG::BeamParticleData pi-
setup pi- -211 pi- 0.13957039 2.5283740149914E-17 0 7804.5 -3 0 1 1
makeanti pi- pi+
create ThePEG::ParticleData eta
setup eta 221 eta 0.547862 1.31E-6 1.31E-5 0 0 0 1 0
create ThePEG::ParticleData eta'
setup eta' 331 eta' 0.95778 0.000203 0.00203 0 0 0 1 0
create ThePEG::ParticleData eta_c
setup eta_c 441 eta_c 2.9839 0.032 0.16 0 0 0 1 0
create ThePEG::ParticleData eta_b
setup eta_b 551 eta_b 9.3987 0.01 0.1 0 0 0 1 0
create ThePEG::ParticleData K+
setup K+ 321 K+ 0.493677 5.3173524656427E-17 0 3711 3 0 1 1
create ThePEG::ParticleData K0
setup K0 311 K0 0.497611 1.0E+27 0 0 0 0 1 0
create ThePEG::ParticleData K-
setup K- -321 K- 0.493677 5.3173524656427E-17 0 3711 -3 0 1 1
makeanti K- K+
create ThePEG::ParticleData Kbar0
setup Kbar0 -311 Kbar0 0.497611 1.0E+27 0 0 0 0 1 0
makeanti Kbar0 K0
create ThePEG::ParticleData D0
setup D0 421 D0 1.86484 1.6042841463415E-12 0 0.123 0 0 1 0
create ThePEG::ParticleData D+
setup D+ 411 D+ 1.86966 6.3286385503528E-13 0 0.3118 3 0 1 0
create ThePEG::ParticleData Dbar0
setup Dbar0 -421 Dbar0 1.86484 1.6042841463415E-12 0 0.123 0 0 1 0
makeanti Dbar0 D0
create ThePEG::ParticleData D-
setup D- -411 D- 1.86966 6.3286385503528E-13 0 0.3118 -3 0 1 0
makeanti D- D+
create ThePEG::ParticleData D_s+
setup D_s+ 431 D_s+ 1.9682 1.315513E-12 0 0.15 3 0 1 0
create ThePEG::ParticleData D_s-
setup D_s- -431 D_s- 1.9682 1.315513E-12 0 0.15 -3 0 1 0
makeanti D_s- D_s+
create ThePEG::MixedParticleData B0
setup B0 511 B0 5.27965 4.2897163043478E-13 0 0.46 0 0 1 0
newdef B0:DeltaM 3.337134e-13*GeV
newdef B0:DeltaGamma 0.*GeV
newdef B0:PQMagnitude 1.
newdef B0:PQPhase 1.
newdef B0:ZMagnitude 0.
newdef B0:ZPhase 0.
create ThePEG::ParticleData B+
setup B+ 521 B+ 5.27934 3.9386616766467E-13 0 0.501 3 0 1 0
create ThePEG::MixedParticleData Bbar0
setup Bbar0 -511 Bbar0 5.27965 4.2897163043478E-13 0 0.46 0 0 1 0
makeanti Bbar0 B0
newdef Bbar0:Synchronized 0
newdef Bbar0:DeltaM 3.337134e-13*GeV
newdef Bbar0:DeltaGamma 0.*GeV
newdef Bbar0:PQMagnitude 1.
newdef Bbar0:PQPhase 1.
newdef Bbar0:ZMagnitude 0.
newdef Bbar0:ZPhase 0.
create ThePEG::ParticleData B-
setup B- -521 B- 5.27934 3.9386616766467E-13 0 0.501 -3 0 1 0
makeanti B- B+
newdef B-:Synchronized 0
create ThePEG::MixedParticleData B_s0
setup B_s0 531 B_s0 5.36688 4.5051815068493E-13 0 0.438 0 0 1 0
newdef B_s0:DeltaM 1.169642e-11*GeV
newdef B_s0:DeltaGamma 6.676243e-14*GeV
newdef B_s0:PQMagnitude 1.
newdef B_s0:PQPhase 1.
newdef B_s0:ZMagnitude 0.
newdef B_s0:ZPhase 0.
create ThePEG::MixedParticleData B_sbar0
setup B_sbar0 -531 B_sbar0 5.36688 4.5051815068493E-13 0 0.438 0 0 1 0
makeanti B_sbar0 B_s0
newdef B_sbar0:Synchronized 0
newdef B_sbar0:DeltaM 1.169642e-11*GeV
newdef B_sbar0:DeltaGamma 6.676243e-14*GeV
newdef B_sbar0:PQMagnitude 1.
newdef B_sbar0:PQPhase 1.
newdef B_sbar0:ZMagnitude 0.
newdef B_sbar0:ZPhase 0.
create ThePEG::ParticleData B_c+
setup B_c+ 541 B_c+ 6.286 1.4299054347826E-12 0 0.138 3 0 1 0
create ThePEG::ParticleData B_c-
setup B_c- -541 B_c- 6.286 1.4299054347826E-12 0 0.138 -3 0 1 0
makeanti B_c- B_c+
#
# the 1^3S_1 multiplet
#
create ThePEG::ParticleData rho+
setup rho+ 213 rho+ 0.77526 0.1491 0.45 0 3 0 3 0
newdef rho+:WidthLoCut 0.4
newdef rho+:WidthUpCut 0.5
create ThePEG::ParticleData rho0
setup rho0 113 rho0 0.77526 0.1491 0.45 0 0 0 3 0
newdef rho0:WidthLoCut 0.4
newdef rho0:WidthUpCut 0.5
create ThePEG::ParticleData rho-
setup rho- -213 rho- 0.77526 0.1491 0.45 0 -3 0 3 0
makeanti rho- rho+
newdef rho-:WidthLoCut 0.4
newdef rho-:WidthUpCut 0.5
create ThePEG::ParticleData omega
setup omega 223 omega 0.78266 0.00868 0.0868 0 0 0 3 0
create ThePEG::ParticleData phi
setup phi 333 phi 1.019461 0.004249 0.036245 0 0 0 3 0
newdef phi:WidthLoCut 0.03
newdef phi:WidthUpCut 0.04249
create ThePEG::ParticleData Jpsi
setup Jpsi 443 Jpsi 3.096916 9.34E-5 0.000934 0 0 0 3 0
create ThePEG::ParticleData Upsilon
setup Upsilon 553 Upsilon 9.4603 5.402E-5 0.00054 0 0 0 3 0
create ThePEG::ParticleData K*+
setup K*+ 323 K*+ 0.89167 0.0514 0.382 0 3 0 3 0
newdef K*+:WidthLoCut 0.25
newdef K*+:WidthUpCut 0.514
create ThePEG::ParticleData K*0
setup K*0 313 K*0 0.89555 0.0473 0.3615 0 0 0 3 0
newdef K*0:WidthLoCut 0.25
newdef K*0:WidthUpCut 0.473
create ThePEG::ParticleData K*-
setup K*- -323 K*- 0.89167 0.0514 0.382 0 -3 0 3 0
makeanti K*- K*+
newdef K*-:WidthLoCut 0.25
newdef K*-:WidthUpCut 0.514
create ThePEG::ParticleData K*bar0
setup K*bar0 -313 K*bar0 0.89555 0.0473 0.3615 0 0 0 3 0
makeanti K*bar0 K*0
newdef K*bar0:WidthLoCut 0.25
newdef K*bar0:WidthUpCut 0.473
create ThePEG::ParticleData D*0
setup D*0 423 D*0 2.00685 6.8E-5 0.00068 0 0 0 3 0
create ThePEG::ParticleData D*+
setup D*+ 413 D*+ 2.01026 8.34E-5 0.000834 0 3 0 3 0
create ThePEG::ParticleData D*bar0
setup D*bar0 -423 D*bar0 2.00685 6.8E-5 0.00068 0 0 0 3 0
makeanti D*bar0 D*0
create ThePEG::ParticleData D*-
setup D*- -413 D*- 2.01026 8.34E-5 0.000834 0 -3 0 3 0
makeanti D*- D*+
create ThePEG::ParticleData D_s*+
setup D_s*+ 433 D_s*+ 2.1122 4.4E-5 0.00044 0 3 0 3 0
create ThePEG::ParticleData D_s*-
setup D_s*- -433 D_s*- 2.1122 4.4E-5 0.00044 0 -3 0 3 0
makeanti D_s*- D_s*+
create ThePEG::ParticleData B*0
setup B*0 513 B*0 5.32471 2.4E-7 2.4E-6 0 0 0 3 0
create ThePEG::ParticleData B*+
setup B*+ 523 B*+ 5.32471 7.8E-7 7.8E-6 0 3 0 3 0
create ThePEG::ParticleData B*bar0
setup B*bar0 -513 B*bar0 5.32471 2.4E-7 2.4E-6 0 0 0 3 0
makeanti B*bar0 B*0
create ThePEG::ParticleData B*-
setup B*- -523 B*- 5.32471 7.8E-7 7.8E-6 0 -3 0 3 0
makeanti B*- B*+
create ThePEG::ParticleData B_s*0
setup B_s*0 533 B_s*0 5.4154 1.5E-7 1.5E-6 0 0 0 3 0
create ThePEG::ParticleData B_s*bar0
setup B_s*bar0 -533 B_s*bar0 5.4154 1.5E-7 1.5E-6 0 0 0 3 0
makeanti B_s*bar0 B_s*0
create ThePEG::ParticleData B_c*+
setup B_c*+ 543 B_c*+ 6.321 8.0E-8 8.0E-7 0 3 0 3 0
create ThePEG::ParticleData B_c*-
setup B_c*- -543 B_c*- 6.321 8.0E-8 8.0E-7 0 -3 0 3 0
makeanti B_c*- B_c*+
#
# the 1^1P_1 multiplet
#
create ThePEG::ParticleData b_1+
setup b_1+ 10213 b_1+ 1.2295 0.142 0.505 0 3 0 3 0
newdef b_1+:WidthLoCut 0.3
newdef b_1+:WidthUpCut 0.71
create ThePEG::ParticleData b_10
setup b_10 10113 b_10 1.2295 0.142 0.505 0 0 0 3 0
newdef b_10:WidthLoCut 0.3
newdef b_10:WidthUpCut 0.71
create ThePEG::ParticleData b_1-
setup b_1- -10213 b_1- 1.2295 0.142 0.505 0 -3 0 3 0
makeanti b_1- b_1+
newdef b_1-:WidthLoCut 0.3
newdef b_1-:WidthUpCut 0.71
create ThePEG::ParticleData h_1
setup h_1 10223 h_1 1.166 0.375 0.375 0 0 0 3 0
create ThePEG::ParticleData h'_1
setup h'_1 10333 h'_1 1.416 0.09 0.18 0 0 0 3 0
create ThePEG::ParticleData h_c
setup h_c 10443 h_c 3.52538 0.0007 0.007 0 0 0 3 0
create ThePEG::ParticleData h_b
setup h_b 10553 h_b 9.8993 8.94E-5 0.000894 0 0 0 3 0
create ThePEG::ParticleData K_1+
setup K_1+ 10323 K_1+ 1.27 0.09 0.18 0 3 0 3 0
newdef K_1+:VariableRatio 1
create ThePEG::ParticleData K_10
setup K_10 10313 K_10 1.272 0.09 0.18 0 0 0 3 0
newdef K_10:VariableRatio 1
create ThePEG::ParticleData K_1-
setup K_1- -10323 K_1- 1.272 0.09 0.18 0 -3 0 3 0
makeanti K_1- K_1+
newdef K_1-:VariableRatio 1
create ThePEG::ParticleData K_1bar0
setup K_1bar0 -10313 K_1bar0 1.272 0.09 0.18 0 0 0 3 0
makeanti K_1bar0 K_10
newdef K_1bar0:VariableRatio 1
create ThePEG::ParticleData D_10
setup D_10 10423 D_10 2.4221 0.0313 0.0939 0 0 0 3 0
create ThePEG::ParticleData D_1+
setup D_1+ 10413 D_1+ 2.4221 0.0313 0.0939 0 3 0 3 0
create ThePEG::ParticleData D_1bar0
setup D_1bar0 -10423 D_1bar0 2.4221 0.0313 0.0939 0 0 0 3 0
makeanti D_1bar0 D_10
create ThePEG::ParticleData D_1-
setup D_1- -10413 D_1- 2.4221 0.0313 0.0939 0 -3 0 3 0
makeanti D_1- D_1+
create ThePEG::ParticleData D_s1+
setup D_s1+ 10433 D_s1+ 2.53511 0.00092 0.0092 0 3 0 3 0
create ThePEG::ParticleData D_s1-
setup D_s1- -10433 D_s1- 2.53511 0.00092 0.0092 0 -3 0 3 0
makeanti D_s1- D_s1+
create ThePEG::ParticleData B_10
setup B_10 10513 B_10 5.7261 0.0275 0.0825 0 0 0 3 0
create ThePEG::ParticleData B_1+
setup B_1+ 10523 B_1+ 5.7259 0.031 0.093 0 3 0 3 0
create ThePEG::ParticleData B_1bar0
setup B_1bar0 -10513 B_1bar0 5.7261 0.0275 0.0825 0 0 0 3 0
makeanti B_1bar0 B_10
create ThePEG::ParticleData B_1-
setup B_1- -10523 B_1- 5.7259 0.031 0.093 0 -3 0 3 0
makeanti B_1- B_1+
create ThePEG::ParticleData B_s10
setup B_s10 10533 B_s10 5.8287 0.0005 0.005 0 0 0 3 0
create ThePEG::ParticleData B_s1bar0
setup B_s1bar0 -10533 B_s1bar0 5.8287 0.0005 0.005 0 0 0 3 0
makeanti B_s1bar0 B_s10
create ThePEG::ParticleData B_c1+
setup B_c1+ 10543 B_c1+ 6.743 7.3E-5 0.00073 0 3 0 3 0
create ThePEG::ParticleData B_c1-
setup B_c1- -10543 B_c1- 6.743 7.3E-5 0.00073 0 -3 0 3 0
makeanti B_c1- B_c1+
#
# the 1^3P_0 multiplet
#
create ThePEG::ParticleData a'_0+
setup a'_0+ 10211 a'_0+ 1.474 0.265 0.265 0 3 0 1 0
create ThePEG::ParticleData a'_00
setup a'_00 10111 a'_00 1.474 0.265 0.265 0 0 0 1 0
create ThePEG::ParticleData a'_0-
setup a'_0- -10211 a'_0- 1.474 0.265 0.265 0 -3 0 1 0
makeanti a'_0- a'_0+
create ThePEG::ParticleData f'_0
setup f'_0 10221 f'_0 1.395 0.275 0.275 0 0 0 1 0
create ThePEG::ParticleData f_0(1710)
setup f_0(1710) 10331 f_0(1710) 1.704 0.123 0.369 0 0 0 1 0
create ThePEG::ParticleData chi_c0
setup chi_c0 10441 chi_c0 3.41476 0.0104 0.104 0 0 0 1 0
create ThePEG::ParticleData chi_b0
setup chi_b0 10551 chi_b0 9.85944 0.000817 0.00817 0 0 0 1 0
create ThePEG::ParticleData K*_0+
setup K*_0+ 10321 K*_0+ 1.425 0.27 0.79 0 3 0 1 0
newdef K*_0+:WidthLoCut 0.77
newdef K*_0+:WidthUpCut 0.81
create ThePEG::ParticleData K*_00
setup K*_00 10311 K*_00 1.425 0.27 0.79 0 0 0 1 0
newdef K*_00:WidthLoCut 0.77
newdef K*_00:WidthUpCut 0.81
create ThePEG::ParticleData K*_0-
setup K*_0- -10321 K*_0- 1.425 0.27 0.79 0 -3 0 1 0
makeanti K*_0- K*_0+
newdef K*_0-:WidthLoCut 0.77
newdef K*_0-:WidthUpCut 0.81
create ThePEG::ParticleData K*_0bar0
setup K*_0bar0 -10311 K*_0bar0 1.425 0.27 0.79 0 0 0 1 0
makeanti K*_0bar0 K*_00
newdef K*_0bar0:WidthLoCut 0.77
newdef K*_0bar0:WidthUpCut 0.81
create ThePEG::ParticleData D_0*+
setup D_0*+ 10411 D_0*+ 2.343 0.229 0.458 0 3 0 1 0
newdef D_0*+:WidthLoCut 0.229
newdef D_0*+:WidthUpCut 0.687
create ThePEG::ParticleData D_0*0
setup D_0*0 10421 D_0*0 2.343 0.229 0.458 0 0 0 1 0
newdef D_0*0:WidthLoCut 0.229
newdef D_0*0:WidthUpCut 0.687
create ThePEG::ParticleData D_0*-
setup D_0*- -10411 D_0*- 2.343 0.229 0.458 0 -3 0 1 0
makeanti D_0*- D_0*+
newdef D_0*-:WidthLoCut 0.229
newdef D_0*-:WidthUpCut 0.687
create ThePEG::ParticleData D_0*bar0
setup D_0*bar0 -10421 D_0*bar0 2.343 0.229 0.458 0 0 0 1 0
makeanti D_0*bar0 D_0*0
newdef D_0*bar0:WidthLoCut 0.229
newdef D_0*bar0:WidthUpCut 0.687
create ThePEG::ParticleData D_s0+
setup D_s0+ 10431 D_s0+ 2.3178 2.32E-5 0.000232 0 3 0 1 0
create ThePEG::ParticleData D_s0-
setup D_s0- -10431 D_s0- 2.3178 2.32E-5 0.000232 0 -3 0 1 0
makeanti D_s0- D_s0+
create ThePEG::ParticleData B*_00
setup B*_00 10511 B*_00 5.726 0.14 0.28 0 0 0 1 0
create ThePEG::ParticleData B*_0+
setup B*_0+ 10521 B*_0+ 5.726 0.14 0.28 0 3 0 1 0
create ThePEG::ParticleData B*_0bar0
setup B*_0bar0 -10511 B*_0bar0 5.726 0.14 0.28 0 0 0 1 0
makeanti B*_0bar0 B*_00
create ThePEG::ParticleData B*_0-
setup B*_0- -10521 B*_0- 5.726 0.14 0.28 0 -3 0 1 0
makeanti B*_0- B*_0+
create ThePEG::ParticleData B*_s00
setup B*_s00 10531 B*_s00 5.818 0.07 0.0875 0 0 0 1 0
newdef B*_s00:WidthLoCut 0.035
newdef B*_s00:WidthUpCut 0.14
create ThePEG::ParticleData B*_s0bar0
setup B*_s0bar0 -10531 B*_s0bar0 5.818 0.07 0.0875 0 0 0 1 0
makeanti B*_s0bar0 B*_s00
newdef B*_s0bar0:WidthLoCut 0.035
newdef B*_s0bar0:WidthUpCut 0.14
create ThePEG::ParticleData B*_c0+
setup B*_c0+ 10541 B*_c0+ 6.727 5.5E-5 0.00055 0 3 0 1 0
create ThePEG::ParticleData B*_c0-
setup B*_c0- -10541 B*_c0- 6.727 5.5E-5 0.00055 0 -3 0 1 0
makeanti B*_c0- B*_c0+
#
# the 1^3P_1 multiplet
#
create ThePEG::ParticleData a_1+
setup a_1+ 20213 a_1+ 1.23 0.42 0.56 0 3 0 3 0
newdef a_1+:VariableRatio 1
create ThePEG::ParticleData a_10
setup a_10 20113 a_10 1.23 0.42 0.56 0 0 0 3 0
newdef a_10:VariableRatio 1
create ThePEG::ParticleData a_1-
setup a_1- -20213 a_1- 1.23 0.42 0.56 0 -3 0 3 0
makeanti a_1- a_1+
newdef a_1-:VariableRatio 1
create ThePEG::ParticleData f_1
setup f_1 20223 f_1 1.2819 0.0227 0.242 0 0 0 3 0
create ThePEG::ParticleData f'_1
setup f'_1 20333 f'_1 1.4263 0.0545 0.40875 0 0 0 3 0
newdef f'_1:WidthLoCut 0.2725
newdef f'_1:WidthUpCut 0.545
create ThePEG::ParticleData chi_c1
setup chi_c1 20443 chi_c1 3.51067 0.00084 0.0089 0 0 0 3 0
create ThePEG::ParticleData chi_b1
setup chi_b1 20553 chi_b1 9.89278 7.1E-5 0.00071 0 0 0 3 0
create ThePEG::ParticleData K'_1+
setup K'_1+ 20323 K'_1+ 1.403 0.174 0.348 0 3 0 3 0
newdef K'_1+:VariableRatio 1
create ThePEG::ParticleData K'_10
setup K'_10 20313 K'_10 1.403 0.174 0.348 0 0 0 3 0
newdef K'_10:VariableRatio 1
create ThePEG::ParticleData K'_1-
setup K'_1- -20323 K'_1- 1.403 0.174 0.348 0 -3 0 3 0
makeanti K'_1- K'_1+
newdef K'_1-:VariableRatio 1
create ThePEG::ParticleData K'_1bar0
setup K'_1bar0 -20313 K'_1bar0 1.403 0.174 0.348 0 0 0 3 0
makeanti K'_1bar0 K'_10
newdef K'_1bar0:VariableRatio 1
create ThePEG::ParticleData D'_10
setup D'_10 20423 D'_10 2.412 0.314 0.439 0 0 0 3 0
newdef D'_10:WidthLoCut 0.25
newdef D'_10:WidthUpCut 0.628
create ThePEG::ParticleData D'_1+
setup D'_1+ 20413 D'_1+ 2.412 0.314 0.439 0 3 0 3 0
newdef D'_1+:WidthLoCut 0.25
newdef D'_1+:WidthUpCut 0.628
create ThePEG::ParticleData D'_1bar0
setup D'_1bar0 -20423 D'_1bar0 2.412 0.314 0.439 0 0 0 3 0
makeanti D'_1bar0 D'_10
newdef D'_1bar0:WidthLoCut 0.25
newdef D'_1bar0:WidthUpCut 0.628
create ThePEG::ParticleData D'_1-
setup D'_1- -20413 D'_1- 2.412 0.314 0.439 0 -3 0 3 0
makeanti D'_1- D'_1+
newdef D'_1-:WidthLoCut 0.25
newdef D'_1-:WidthUpCut 0.628
create ThePEG::ParticleData D'_s1+
setup D'_s1+ 20433 D'_s1+ 2.4595 3.82E-5 0.000382 0 3 0 3 0
create ThePEG::ParticleData D'_s1-
setup D'_s1- -20433 D'_s1- 2.4595 3.82E-5 0.000382 0 -3 0 3 0
makeanti D'_s1- D'_s1+
create ThePEG::ParticleData B'_10
setup B'_10 20513 B'_10 5.762 0.13 0.26 0 0 0 3 0
create ThePEG::ParticleData B'_1+
setup B'_1+ 20523 B'_1+ 5.762 0.13 0.26 0 3 0 3 0
create ThePEG::ParticleData B'_1bar0
setup B'_1bar0 -20513 B'_1bar0 5.762 0.13 0.26 0 0 0 3 0
makeanti B'_1bar0 B'_10
create ThePEG::ParticleData B'_1-
setup B'_1- -20523 B'_1- 5.762 0.13 0.26 0 -3 0 3 0
makeanti B'_1- B'_1+
create ThePEG::ParticleData B'_s10
setup B'_s10 20533 B'_s10 5.856 0.06 0.075 0 0 0 3 0
newdef B'_s10:WidthLoCut 0.03
newdef B'_s10:WidthUpCut 0.12
create ThePEG::ParticleData B'_s1bar0
setup B'_s1bar0 -20533 B'_s1bar0 5.856 0.06 0.075 0 0 0 3 0
makeanti B'_s1bar0 B'_s10
newdef B'_s1bar0:WidthLoCut 0.03
newdef B'_s1bar0:WidthUpCut 0.12
create ThePEG::ParticleData B'_c1+
setup B'_c1+ 20543 B'_c1+ 6.765 9.1E-5 0.00091 0 3 0 3 0
create ThePEG::ParticleData B'_c1-
setup B'_c1- -20543 B'_c1- 6.765 9.1E-5 0.00091 0 -3 0 3 0
makeanti B'_c1- B'_c1+
#
# the 1^3P_2 multiplet
#
create ThePEG::ParticleData a_2+
setup a_2+ 215 a_2+ 1.3182 0.105 0.21 0 3 0 5 0
create ThePEG::ParticleData a_20
setup a_20 115 a_20 1.3182 0.105 0.21 0 0 0 5 0
create ThePEG::ParticleData a_2-
setup a_2- -215 a_2- 1.3182 0.105 0.21 0 -3 0 5 0
makeanti a_2- a_2+
create ThePEG::ParticleData f_2
setup f_2 225 f_2 1.2755 0.1867 0.18085 0 0 0 5 0
newdef f_2:WidthLoCut 0.175
newdef f_2:WidthUpCut 0.1867
create ThePEG::ParticleData f'_2
setup f'_2 335 f'_2 1.5174 0.086 0.344 0 0 0 5 0
create ThePEG::ParticleData chi_c2
setup chi_c2 445 chi_c2 3.55617 0.00197 0.0206 0 0 0 5 0
create ThePEG::ParticleData chi_b2
setup chi_b2 555 chi_b2 9.91221 0.00017 0.0017 0 0 0 5 0
create ThePEG::ParticleData K*_2+
setup K*_2+ 325 K*_2+ 1.4273 0.1 0.2 0 3 0 5 0
create ThePEG::ParticleData K*_20
setup K*_20 315 K*_20 1.4324 0.109 0.218 0 0 0 5 0
create ThePEG::ParticleData K*_2-
setup K*_2- -325 K*_2- 1.4273 0.1 0.2 0 -3 0 5 0
makeanti K*_2- K*_2+
create ThePEG::ParticleData K*_2bar0
setup K*_2bar0 -315 K*_2bar0 1.4324 0.109 0.218 0 0 0 5 0
makeanti K*_2bar0 K*_20
create ThePEG::ParticleData D*_20
setup D*_20 425 D*_20 2.4611 0.0473 0.2365 0 0 0 5 0
create ThePEG::ParticleData D*_2+
setup D*_2+ 415 D*_2+ 2.4611 0.0473 0.2365 0 3 0 5 0
create ThePEG::ParticleData D*_2bar0
setup D*_2bar0 -425 D*_2bar0 2.4611 0.0473 0.2365 0 0 0 5 0
makeanti D*_2bar0 D*_20
create ThePEG::ParticleData D*_2-
setup D*_2- -415 D*_2- 2.4611 0.0473 0.2365 0 -3 0 5 0
makeanti D*_2- D*_2+
create ThePEG::ParticleData D_s2+
setup D_s2+ 435 D_s2+ 2.5691 0.0169 0.0507 0 3 0 5 0
create ThePEG::ParticleData D_s2-
setup D_s2- -435 D_s2- 2.5691 0.0169 0.0507 0 -3 0 5 0
makeanti D_s2- D_s2+
create ThePEG::ParticleData B_20
setup B_20 515 B_20 5.7395 0.0242 0.121 0 0 0 5 0
create ThePEG::ParticleData B_2+
setup B_2+ 525 B_2+ 5.7372 0.0242 0.121 0 3 0 5 0
create ThePEG::ParticleData B_2bar0
setup B_2bar0 -515 B_2bar0 5.7395 0.0242 0.121 0 0 0 5 0
makeanti B_2bar0 B_20
create ThePEG::ParticleData B_2-
setup B_2- -525 B_2- 5.7372 0.0242 0.121 0 -3 0 5 0
makeanti B_2- B_2+
create ThePEG::ParticleData B_s20
setup B_s20 535 B_s20 5.83986 0.00149 0.0149 0 0 0 5 0
create ThePEG::ParticleData B_s2bar0
setup B_s2bar0 -535 B_s2bar0 5.83986 0.00149 0.0149 0 0 0 5 0
makeanti B_s2bar0 B_s20
create ThePEG::ParticleData B_c2+
setup B_c2+ 545 B_c2+ 6.783 8.3E-5 0.00083 0 3 0 5 0
create ThePEG::ParticleData B_c2-
setup B_c2- -545 B_c2- 6.783 8.3E-5 0.00083 0 -3 0 5 0
makeanti B_c2- B_c2+
#
# the 1^1D_2 multiplet
#
create ThePEG::ParticleData pi_2+
setup pi_2+ 10215 pi_2+ 1.6706 0.258 0.258 0 3 0 5 0
create ThePEG::ParticleData pi_20
setup pi_20 10115 pi_20 1.6706 0.258 0.258 0 0 0 5 0
create ThePEG::ParticleData pi_2-
setup pi_2- -10215 pi_2- 1.6706 0.258 0.258 0 -3 0 5 0
makeanti pi_2- pi_2+
create ThePEG::ParticleData eta_2
setup eta_2 10225 eta_2 1.617 0.181 0.362 0 0 0 5 0
create ThePEG::ParticleData eta'_2
setup eta'_2 10335 eta'_2 1.842 0.225 0.225 0 0 0 5 0
create ThePEG::ParticleData eta_b2
setup eta_b2 10555 eta_b2 10.158 3.2E-5 0.00032 0 0 0 5 0
create ThePEG::ParticleData K_2(1770)+
setup K_2(1770)+ 10325 K_2(1770)+ 1.773 0.186 0.558 0 3 0 5 0
create ThePEG::ParticleData K_2(1770)0
setup K_2(1770)0 10315 K_2(1770)0 1.773 0.186 0.558 0 0 0 5 0
create ThePEG::ParticleData K_2(1770)-
setup K_2(1770)- -10325 K_2(1770)- 1.773 0.186 0.558 0 -3 0 5 0
makeanti K_2(1770)- K_2(1770)+
create ThePEG::ParticleData K_2(1770)bar0
setup K_2(1770)bar0 -10315 K_2(1770)bar0 1.773 0.186 0.558 0 0 0 5 0
makeanti K_2(1770)bar0 K_2(1770)0
create ThePEG::ParticleData B_c2(L)+
setup B_c2(L)+ 10545 B_c2(L)+ 7.036 8.3E-5 0.00083 0 3 0 5 0
create ThePEG::ParticleData B_c2(L)-
setup B_c2(L)- -10545 B_c2(L)- 7.036 8.3E-5 0.00083 0 -3 0 5 0
makeanti B_c2(L)- B_c2(L)+
#
# the 1^3D_1 multiplet
#
create ThePEG::ParticleData rho''+
setup rho''+ 30213 rho''+ 1.72 0.25 0.5 0 3 0 3 0
create ThePEG::ParticleData rho''0
setup rho''0 30113 rho''0 1.72 0.25 0.5 0 0 0 3 0
create ThePEG::ParticleData rho''-
setup rho''- -30213 rho''- 1.72 0.25 0.5 0 -3 0 3 0
makeanti rho''- rho''+
create ThePEG::ParticleData omega''
setup omega'' 30223 omega'' 1.67 0.315 0.63 0 0 0 3 0
create ThePEG::ParticleData phi(2170)
setup phi(2170) 30333 phi(2170) 2.162 0.1 0.3 0 0 0 3 0
create ThePEG::ParticleData psi(3770)
setup psi(3770) 30443 psi(3770) 3.7737 0.0272 0.0816 0 0 0 3 0
newdef psi(3770):WidthLoCut 0.0272
newdef psi(3770):WidthUpCut 0.136
create ThePEG::ParticleData Upsilon_1(1D)
setup Upsilon_1(1D) 30553 Upsilon_1(1D) 10.1551 3.56E-5 0.000356 0 0 0 3 0
create ThePEG::ParticleData K''*+
setup K''*+ 30323 K''*+ 1.718 0.322 0.322 0 3 0 3 0
create ThePEG::ParticleData K''*0
setup K''*0 30313 K''*0 1.718 0.322 0.322 0 0 0 3 0
create ThePEG::ParticleData K''*-
setup K''*- -30323 K''*- 1.718 0.322 0.322 0 -3 0 3 0
makeanti K''*- K''*+
create ThePEG::ParticleData K''*bar0
setup K''*bar0 -30313 K''*bar0 1.718 0.322 0.322 0 0 0 3 0
makeanti K''*bar0 K''*0
create ThePEG::ParticleData D_1*(2760)0
setup D_1*(2760)0 30423 D_1*(2760)0 2.781 0.177 0.354 0 0 0 3 0
create ThePEG::ParticleData D_1*(2760)+
setup D_1*(2760)+ 30413 D_1*(2760)+ 2.781 0.177 0.354 0 3 0 3 0
create ThePEG::ParticleData Dbar_1*(2760)0
setup Dbar_1*(2760)0 -30423 Dbar_1*(2760)0 2.781 0.177 0.354 0 0 0 3 0
makeanti Dbar_1*(2760)0 D_1*(2760)0
create ThePEG::ParticleData D_1*(2760)-
setup D_1*(2760)- -30413 D_1*(2760)- 2.781 0.177 0.354 0 -3 0 3 0
makeanti D_1*(2760)- D_1*(2760)+
create ThePEG::ParticleData D_s1*(2860)+
setup D_s1*(2860)+ 30433 D_s1*(2860)+ 2.859 0.159 0.318 0 3 0 3 0
create ThePEG::ParticleData D_s1*(2860)-
setup D_s1*(2860)- -30433 D_s1*(2860)- 2.859 0.159 0.318 0 -3 0 3 0
makeanti D_s1*(2860)- D_s1*(2860)+
create ThePEG::ParticleData B_c(1D)*+
setup B_c(1D)*+ 30543 B_c(1D)*+ 7.028 9.35E-5 0.000935 0 3 0 3 0
create ThePEG::ParticleData B_c(1D)*-
setup B_c(1D)*- -30543 B_c(1D)*- 7.028 9.35E-5 0.000935 0 -3 0 3 0
makeanti B_c(1D)*- B_c(1D)*+
#
# the 1^3D_2 multiplet
#
create ThePEG::ParticleData psi_2(1D)
setup psi_2(1D) 20445 psi_2(1D) 3.8237 0.00074 0.0074 0 0 0 5 0
create ThePEG::ParticleData Upsilon_2(1D)
setup Upsilon_2(1D) 20555 Upsilon_2(1D) 10.1637 2.8E-5 0.00028 0 0 0 5 0
create ThePEG::ParticleData K_2(1820)+
setup K_2(1820)+ 20325 K_2(1820)+ 1.819 0.264 0.528 0 3 0 5 0
create ThePEG::ParticleData K_2(1820)0
setup K_2(1820)0 20315 K_2(1820)0 1.819 0.264 0.528 0 0 0 5 0
create ThePEG::ParticleData K_2(1820)-
setup K_2(1820)- -20325 K_2(1820)- 1.819 0.264 0.528 0 -3 0 5 0
makeanti K_2(1820)- K_2(1820)+
create ThePEG::ParticleData K_2(1820)bar0
setup K_2(1820)bar0 -20315 K_2(1820)bar0 1.819 0.264 0.528 0 0 0 5 0
makeanti K_2(1820)bar0 K_2(1820)0
create ThePEG::ParticleData D_2(2740)0
setup D_2(2740)0 20425 D_2(2740)0 2.747 0.088 0.264 0 0 0 5 0
create ThePEG::ParticleData D_2(2740)+
setup D_2(2740)+ 20415 D_2(2740)+ 2.747 0.088 0.264 0 3 0 5 0
create ThePEG::ParticleData Dbar_2(2740)0
setup Dbar_2(2740)0 -20425 Dbar_2(2740)0 2.747 0.088 0.264 0 0 0 5 0
makeanti Dbar_2(2740)0 D_2(2740)0
create ThePEG::ParticleData D_2(2740)-
setup D_2(2740)- -20415 D_2(2740)- 2.747 0.088 0.264 0 -3 0 5 0
makeanti D_2(2740)- D_2(2740)+
create ThePEG::ParticleData D_s2(2850)+
setup D_s2(2850)+ 20435 D_s2(2850)+ 2.851 0.052 0.156 0 3 0 5 0
create ThePEG::ParticleData D_s2(2850)-
setup D_s2(2850)- -20435 D_s2(2850)- 2.851 0.052 0.156 0 -3 0 5 0
makeanti D_s2(2850)- D_s2(2850)+
create ThePEG::ParticleData B_c2(H)+
setup B_c2(H)+ 20545 B_c2(H)+ 7.041 9.3E-5 0.00093 0 3 0 5 0
create ThePEG::ParticleData B_c2(H)-
setup B_c2(H)- -20545 B_c2(H)- 7.041 9.3E-5 0.00093 0 -3 0 5 0
makeanti B_c2(H)- B_c2(H)+
#
# the 1^3D_3 multiplet
#
create ThePEG::ParticleData rho_3+
setup rho_3+ 217 rho_3+ 1.6888 0.161 0.592 0 3 0 7 0
newdef rho_3+:WidthLoCut 0.54
newdef rho_3+:WidthUpCut 0.644
create ThePEG::ParticleData rho_30
setup rho_30 117 rho_30 1.6888 0.161 0.592 0 0 0 7 0
newdef rho_30:WidthLoCut 0.54
newdef rho_30:WidthUpCut 0.644
create ThePEG::ParticleData rho_3-
setup rho_3- -217 rho_3- 1.6888 0.161 0.592 0 -3 0 7 0
makeanti rho_3- rho_3+
newdef rho_3-:WidthLoCut 0.54
newdef rho_3-:WidthUpCut 0.644
create ThePEG::ParticleData omega_3
setup omega_3 227 omega_3 1.667 0.168 0.336 0 0 0 7 0
create ThePEG::ParticleData phi_3
setup phi_3 337 phi_3 1.854 0.087 0.435 0 0 0 7 0
create ThePEG::ParticleData psi_3(1D)
setup psi_3(1D) 447 psi_3(1D) 3.84271 0.0028 0.014 0 0 0 7 0
create ThePEG::ParticleData Upsilon_3(1D)
setup Upsilon_3(1D) 557 Upsilon_3(1D) 10.1651 2.55E-5 0.000255 0 0 0 7 0
create ThePEG::ParticleData K_3*+
setup K_3*+ 327 K_3*+ 1.779 0.161 0.477 0 3 0 7 0
create ThePEG::ParticleData K_3*0
setup K_3*0 317 K_3*0 1.779 0.161 0.477 0 0 0 7 0
create ThePEG::ParticleData K_3*-
setup K_3*- -327 K_3*- 1.779 0.161 0.477 0 -3 0 7 0
makeanti K_3*- K_3*+
create ThePEG::ParticleData K_3*bar0
setup K_3*bar0 -317 K_3*bar0 1.779 0.161 0.477 0 0 0 7 0
makeanti K_3*bar0 K_3*0
create ThePEG::ParticleData D_3*(2750)0
setup D_3*(2750)0 427 D_3*(2750)0 2.7631 0.066 0.192 0 0 0 7 0
create ThePEG::ParticleData D_3*(2750)+
setup D_3*(2750)+ 417 D_3*(2750)+ 2.7631 0.066 0.192 0 3 0 7 0
create ThePEG::ParticleData Dbar_3*(2750)0
setup Dbar_3*(2750)0 -427 Dbar_3*(2750)0 2.7631 0.066 0.192 0 0 0 7 0
makeanti Dbar_3*(2750)0 D_3*(2750)0
create ThePEG::ParticleData D_3*(2750)-
setup D_3*(2750)- -417 D_3*(2750)- 2.7631 0.066 0.192 0 -3 0 7 0
makeanti D_3*(2750)- D_3*(2750)+
create ThePEG::ParticleData D_s3*(2860)+
setup D_s3*(2860)+ 437 D_s3*(2860)+ 2.86 0.053 0.159 0 3 0 7 0
create ThePEG::ParticleData D_s3*(2860)-
setup D_s3*(2860)- -437 D_s3*(2860)- 2.86 0.053 0.159 0 -3 0 7 0
makeanti D_s3*(2860)- D_s3*(2860)+
create ThePEG::ParticleData B_c3(1D)*+
setup B_c3(1D)*+ 547 B_c3(1D)*+ 7.045 8.23E-5 0.000823 0 3 0 7 0
create ThePEG::ParticleData B_c3(1D)*-
setup B_c3(1D)*- -547 B_c3(1D)*- 7.045 8.23E-5 0.000823 0 -3 0 7 0
makeanti B_c3(1D)*- B_c3(1D)*+
#
# the 1^3F_4 multiplet
#
#
# the 2^1S_0 multiplet
#
create ThePEG::ParticleData pi'+
setup pi'+ 100211 pi'+ 1.3 0.4 0.4 0 3 0 1 0
create ThePEG::ParticleData pi'0
setup pi'0 100111 pi'0 1.3 0.4 0.4 0 0 0 1 0
create ThePEG::ParticleData pi'-
setup pi'- -100211 pi'- 1.3 0.4 0.4 0 -3 0 1 0
makeanti pi'- pi'+
create ThePEG::ParticleData eta(1295)
setup eta(1295) 100221 eta(1295) 1.294 0.055 0.395 0 0 0 1 0
newdef eta(1295):WidthLoCut 0.24
newdef eta(1295):WidthUpCut 0.55
create ThePEG::ParticleData eta(1475)
setup eta(1475) 100331 eta(1475) 1.475 0.09 0.174 0 0 0 1 0
create ThePEG::ParticleData eta_c(2S)
setup eta_c(2S) 100441 eta_c(2S) 3.6375 0.0113 0.113 0 0 0 1 0
create ThePEG::ParticleData eta_b(2S)
setup eta_b(2S) 100551 eta_b(2S) 9.999 0.0041 0.041 0 0 0 1 0
create ThePEG::ParticleData K'+
setup K'+ 100321 K'+ 1.46 0.26 0.26 0 3 0 1 0
create ThePEG::ParticleData K'0
setup K'0 100311 K'0 1.46 0.26 0.26 0 0 0 1 0
create ThePEG::ParticleData K'-
setup K'- -100321 K'- 1.46 0.26 0.26 0 -3 0 1 0
makeanti K'- K'+
create ThePEG::ParticleData K'bar0
setup K'bar0 -100311 K'bar0 1.46 0.26 0.26 0 0 0 1 0
makeanti K'bar0 K'0
create ThePEG::ParticleData D_0(2550)0
setup D_0(2550)0 100421 D_0(2550)0 2.549 0.165 0.33 0 0 0 1 0
create ThePEG::ParticleData D_0(2550)+
setup D_0(2550)+ 100411 D_0(2550)+ 2.549 0.165 0.33 0 3 0 1 0
create ThePEG::ParticleData Dbar_0(2550)0
setup Dbar_0(2550)0 -100421 Dbar_0(2550)0 2.549 0.165 0.33 0 0 0 1 0
makeanti Dbar_0(2550)0 D_0(2550)0
create ThePEG::ParticleData D_0(2550)-
setup D_0(2550)- -100411 D_0(2550)- 2.549 0.165 0.33 0 -3 0 1 0
makeanti D_0(2550)- D_0(2550)+
create ThePEG::ParticleData D_s0(2590)+
setup D_s0(2590)+ 100431 D_s0(2590)+ 2.591 0.089 0.267 0 3 0 1 0
create ThePEG::ParticleData D_s0(2590)-
setup D_s0(2590)- -100431 D_s0(2590)- 2.591 0.089 0.267 0 -3 0 1 0
makeanti D_s0(2590)- D_s0(2590)+
create ThePEG::ParticleData B_c(2S)+
setup B_c(2S)+ 100541 B_c(2S)+ 6.8712 6.5E-5 0.00065 0 3 0 1 0
create ThePEG::ParticleData B_c(2S)-
setup B_c(2S)- -100541 B_c(2S)- 6.8712 6.5E-5 0.00065 0 -3 0 1 0
makeanti B_c(2S)- B_c(2S)+
#
# the 2^3S_1 multiplet
#
create ThePEG::ParticleData rho'+
setup rho'+ 100213 rho'+ 1.465 0.4 0.4 0 3 0 3 0
create ThePEG::ParticleData rho'0
setup rho'0 100113 rho'0 1.465 0.4 0.4 0 0 0 3 0
create ThePEG::ParticleData rho'-
setup rho'- -100213 rho'- 1.465 0.4 0.4 0 -3 0 3 0
makeanti rho'- rho'+
create ThePEG::ParticleData omega'
setup omega' 100223 omega' 1.41 0.29 0.3225 0 0 0 3 0
newdef omega':WidthLoCut 0.215
newdef omega':WidthUpCut 0.43
create ThePEG::ParticleData phi'
setup phi' 100333 phi' 1.68 0.15 0.3 0 0 0 3 0
create ThePEG::ParticleData psi(2S)
setup psi(2S) 100443 psi(2S) 3.6861 0.000294 0.00337 0 0 0 3 0
create ThePEG::ParticleData Upsilon(2S)
setup Upsilon(2S) 100553 Upsilon(2S) 10.02326 3.198E-5 0.0003198 0 0 0 3 0
create ThePEG::ParticleData K'*+
setup K'*+ 100323 K'*+ 1.414 0.232 0.464 0 3 0 3 0
create ThePEG::ParticleData K'*0
setup K'*0 100313 K'*0 1.414 0.232 0.464 0 0 0 3 0
create ThePEG::ParticleData K'*-
setup K'*- -100323 K'*- 1.414 0.232 0.464 0 -3 0 3 0
makeanti K'*- K'*+
create ThePEG::ParticleData K'*bar0
setup K'*bar0 -100313 K'*bar0 1.414 0.232 0.464 0 0 0 3 0
makeanti K'*bar0 K'*0
create ThePEG::ParticleData D_1*(2600)0
setup D_1*(2600)0 100423 D_1*(2600)0 2.627 0.141 0.282 0 0 0 3 0
create ThePEG::ParticleData D_1*(2600)+
setup D_1*(2600)+ 100413 D_1*(2600)+ 2.627 0.141 0.282 0 3 0 3 0
create ThePEG::ParticleData Dbar_1*(2600)0
setup Dbar_1*(2600)0 -100423 Dbar_1*(2600)0 2.627 0.141 0.282 0 0 0 3 0
makeanti Dbar_1*(2600)0 D_1*(2600)0
create ThePEG::ParticleData D_1*(2600)-
setup D_1*(2600)- -100413 D_1*(2600)- 2.627 0.141 0.282 0 -3 0 3 0
makeanti D_1*(2600)- D_1*(2600)+
create ThePEG::ParticleData D_s1*(2700)+
setup D_s1*(2700)+ 100433 D_s1*(2700)+ 2.714 0.122 0.244 0 3 0 3 0
create ThePEG::ParticleData D_s1*(2700)-
setup D_s1*(2700)- -100433 D_s1*(2700)- 2.714 0.122 0.244 0 -3 0 3 0
makeanti D_s1*(2700)- D_s1*(2700)+
create ThePEG::ParticleData B_c(2S)*+
setup B_c(2S)*+ 100543 B_c(2S)*+ 6.8752 7.2E-5 0.00072 0 3 0 3 0
create ThePEG::ParticleData B_c(2S)*-
setup B_c(2S)*- -100543 B_c(2S)*- 6.8752 7.2E-5 0.00072 0 -3 0 3 0
makeanti B_c(2S)*- B_c(2S)*+
#
# the 2^1P_1 multiplet
#
create ThePEG::ParticleData h_b(2P)
setup h_b(2P) 110553 h_b(2P) 10.2598 8.0E-5 0.0008 0 0 0 3 0
#
# the 2^3P_0 multiplet
#
create ThePEG::ParticleData chi_b0(2P)
setup chi_b0(2P) 110551 chi_b0(2P) 10.2325 0.000887 0.00887 0 0 0 1 0
#
# the 2^3P_1 multiplet
#
create ThePEG::ParticleData chi_b1(2P)
setup chi_b1(2P) 120553 chi_b1(2P) 10.25546 7.87E-5 0.000787 0 0 0 3 0
#
# the 2^3P_2 multiplet
#
create ThePEG::ParticleData chi_c2(2P)
setup chi_c2(2P) 100445 chi_c2(2P) 3.929 0.029 0.24 0 0 0 5 0
newdef chi_c2(2P):WidthLoCut 0.19
newdef chi_c2(2P):WidthUpCut 0.29
create ThePEG::ParticleData chi_b2(2P)
setup chi_b2(2P) 100555 chi_b2(2P) 10.26865 0.0001847 0.001847 0 0 0 5 0
#
# the 2^1D_2 multiplet
#
#
# the 2^3D_1 multiplet
#
create ThePEG::ParticleData Upsilon_1(2D)
setup Upsilon_1(2D) 130553 Upsilon_1(2D) 10.435 2.64E-5 0.000264 0 0 0 3 0
#
# the 2^3D_2 multiplet
#
create ThePEG::ParticleData Upsilon_2(2D)
setup Upsilon_2(2D) 120555 Upsilon_2(2D) 10.441 2.23E-5 0.000223 0 0 0 5 0
#
# the 2^3D_3 multiplet
#
create ThePEG::ParticleData Upsilon_3(2D)
setup Upsilon_3(2D) 100557 Upsilon_3(2D) 10.444 2.02E-5 0.000202 0 0 0 7 0
#
# the 3^1S_0 multiplet
#
create ThePEG::ParticleData eta_b(3S)
setup eta_b(3S) 200551 eta_b(3S) 10.337 0.0027 0.027 0 0 0 1 0
#
# the 3^3S_1 multiplet
#
create ThePEG::ParticleData Upsilon(3S)
setup Upsilon(3S) 200553 Upsilon(3S) 10.3552 2.032E-5 0.0002032 0 0 0 3 0
#
# the 3^1P_1 multiplet
#
#
# the 3^3P_0 multiplet
#
create ThePEG::ParticleData chi_b0(3P)
setup chi_b0(3P) 210551 chi_b0(3P) 10.5007 0.00071 0.0071 0 0 0 1 0
#
# the 3^3P_1 multiplet
#
create ThePEG::ParticleData chi_b1(3P)
setup chi_b1(3P) 220553 chi_b1(3P) 10.5134 6.6E-5 0.00066 0 0 0 3 0
#
# the 3^3P_2 multiplet
#
create ThePEG::ParticleData chi_b2(3P)
setup chi_b2(3P) 200555 chi_b2(3P) 10.524 0.000146 0.00146 0 0 0 5 0
#
# the 4^3S_1 multiplet
#
create ThePEG::ParticleData Upsilon(4S)
setup Upsilon(4S) 300553 Upsilon(4S) 10.5794 0.0205 0.11275 0 0 0 3 0
newdef Upsilon(4S):WidthLoCut 0.0205
newdef Upsilon(4S):WidthUpCut 0.205
#
# the 5^3S_1 multiplet
#
create ThePEG::ParticleData Upsilon(5S)
setup Upsilon(5S) 9000553 Upsilon(5S) 10.8852 0.037 0.185 0 0 0 3 0
#
# The meson which are not part of a multiplet
#
create ThePEG::ParticleData f_0(1500)
setup f_0(1500) 9030221 f_0(1500) 1.506 0.112 0.336 0 0 0 1 0
newdef f_0(1500):VariableRatio 1
create ThePEG::ParticleData sigma
setup sigma 9000221 sigma 0.86 0.88 1.155 0 0 0 1 0
newdef sigma:WidthLoCut 0.55
newdef sigma:WidthUpCut 1.76
create ThePEG::ParticleData f_0
setup f_0 9010221 f_0 0.965 0.1635 0.384 0 0 0 1 0
newdef f_0:VariableRatio 1
newdef f_0:WidthLoCut 0.256
newdef f_0:WidthUpCut 0.512
create ThePEG::ParticleData a_00
setup a_00 9000111 a_00 0.999 0.1778 0.4 0 0 0 1 0
newdef a_00:VariableRatio 1
newdef a_00:WidthLoCut 0.24
newdef a_00:WidthUpCut 0.56
create ThePEG::ParticleData a_0+
setup a_0+ 9000211 a_0+ 0.999 0.1778 0.4 0 3 0 1 0
newdef a_0+:VariableRatio 1
newdef a_0+:WidthLoCut 0.24
newdef a_0+:WidthUpCut 0.56
create ThePEG::ParticleData a_0-
setup a_0- -9000211 a_0- 0.999 0.1778 0.4 0 -3 0 1 0
makeanti a_0- a_0+
newdef a_0-:VariableRatio 1
newdef a_0-:WidthLoCut 0.24
newdef a_0-:WidthUpCut 0.56
create ThePEG::ParticleData eta(1405)
setup eta(1405) 9020221 eta(1405) 1.4098 0.0511 0.38325 0 0 0 1 0
newdef eta(1405):WidthLoCut 0.2555
newdef eta(1405):WidthUpCut 0.511
create ThePEG::ParticleData K_S0
setup K_S0 310 K_S0 0.497611 7.3514250055883E-15 0 26.842 0 0 1 0
create ThePEG::ParticleData K_L0
setup K_L0 130 K_L0 0.497611 1.2871947162427E-17 0 15330 0 0 1 1
create ThePEG::ParticleData kappa0
setup kappa0 9000311 kappa0 0.845 0.468 0.5895 0 0 0 1 0
newdef kappa0:WidthLoCut 0.207
newdef kappa0:WidthUpCut 0.972
create ThePEG::ParticleData kappabar0
setup kappabar0 -9000311 kappabar0 0.845 0.468 0.5895 0 0 0 1 0
makeanti kappabar0 kappa0
newdef kappabar0:WidthLoCut 0.207
newdef kappabar0:WidthUpCut 0.972
create ThePEG::ParticleData kappa-
setup kappa- -9000321 kappa- 0.845 0.468 0.5895 0 -3 0 1 0
newdef kappa-:WidthLoCut 0.207
newdef kappa-:WidthUpCut 0.972
create ThePEG::ParticleData kappa+
setup kappa+ 9000321 kappa+ 0.845 0.468 0.5895 0 3 0 1 0
makeanti kappa+ kappa-
newdef kappa+:WidthLoCut 0.207
newdef kappa+:WidthUpCut 0.972
diff --git a/src/snippets/Makefile.am b/src/snippets/Makefile.am
--- a/src/snippets/Makefile.am
+++ b/src/snippets/Makefile.am
@@ -1,50 +1,51 @@
BUILT_SOURCES = done-all-links
snippetsdir = ${pkgdatadir}/snippets
INPUTFILES = \
CellGridSampler.in \
Diffraction.in \
DipoleMerging.in \
DipoleShowerFiveFlavours.in \
DipoleShowerFourFlavours.in \
EECollider.in \
EPCollider.in \
PECollider.in \
HepMCFixedOrder.in \
HepMC.in \
Matchbox.in \
MB-DipoleShower.in \
MB.in \
MonacoSampler.in \
Particles-SetLonglivedParticlesStable.in \
PDF-CT10.in \
PDF-NNPDF30NLO.in \
+PionPDF.in \
PPCollider.in \
RivetFixedOrder.in \
Rivet.in \
YFS.in \
CMWinQtiledShower.in \
Dipole_AutoTunes_gss.in \
Tune-DotProduct.in \
Tune-DotProduct-Veto.in \
Tune-pT.in \
Tune-Q2.in \
EvolutionScheme-DotProduct.in \
EvolutionScheme-DotProduct-Veto.in \
EvolutionScheme-pT.in \
EvolutionScheme-Q2.in \
FixedTarget.in \
FixedTarget-PP.in \
H7Hadrons.in
dist_snippets_DATA = $(INPUTFILES)
CLEANFILES = done-all-links
done-all-links: $(INPUTFILES)
@echo "Linking input files"
@for i in $(INPUTFILES); do \
if test -f $(srcdir)/$$i -a ! -e $$i; then \
$(LN_S) -f $(srcdir)/$$i; fi; done
@touch done-all-links
diff --git a/src/snippets/PionPDF.in b/src/snippets/PionPDF.in
new file mode 100644
--- /dev/null
+++ b/src/snippets/PionPDF.in
@@ -0,0 +1,21 @@
+# -*- ThePEG-repository -*-
+
+#=============================================================================
+# Set up for charged pion beams
+#
+# The LHAPDF6 library must be available to use this snippet.
+#=============================================================================
+
+cd /Herwig/Partons
+create ThePEG::LHAPDF PionPDF ThePEGLHAPDF.so
+set PionPDF:PDFName xFitterPI_NLO_EIG
+cp HadronRemnants PionRemnants
+cp RemnantDecayer PionRemnantDecayer
+set PionRemnantDecayer:AllowTop Yes
+set PionRemnants:RemnantDecayer PionRemnantDecayer
+set PionPDF:RemnantHandler PionRemnants
+set /Herwig/Particles/pi+:PDF PionPDF
+set /Herwig/Particles/pi-:PDF PionPDF
+
+cd /
+