Page MenuHomeHEPForge
Files F8723572

No OneTemporary
Actions

View Options

diff --git a/Models/General/BSMModel.cc b/Models/General/BSMModel.cc
--- a/Models/General/BSMModel.cc
+++ b/Models/General/BSMModel.cc
@@ -1,341 +1,347 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the BSMModel class.
//
#include "BSMModel.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/Utilities/StringUtils.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/MassGenerator.h"
#include "ThePEG/PDT/WidthGenerator.h"
#include "ThePEG/PDT/DecayMode.h"
using namespace Herwig;
BSMModel::BSMModel() : decayFile_(), topModesFromFile_(false),
tolerance_(1e-6)
{}
void BSMModel::persistentOutput(PersistentOStream & os) const {
os << decayFile_ << topModesFromFile_ << tolerance_ ;
}
void BSMModel::persistentInput(PersistentIStream & is, int) {
is >> decayFile_ >> topModesFromFile_ >> tolerance_ ;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeAbstractClass<BSMModel,Herwig::StandardModel>
describeHerwigBSMModel("Herwig::BSMModel", "BSMModel.so");
void BSMModel::Init() {
static ClassDocumentation<BSMModel> documentation
("The BSMModel class provides a base class for BSM models including the"
" features to read decays in the SLHA format");
static Parameter<BSMModel,string> interfaceDecayFileName
("DecayFileName",
"Name of the file from which to read decays in the SLHA format",
&BSMModel::decayFile_, "",
false, false);
static Switch<BSMModel,bool> interfaceTopModes
("TopModes",
"Whether ro use the Herwig++ SM top decays or those from the SLHA file",
&BSMModel::topModesFromFile_, false, false, false);
static SwitchOption interfaceTopModesFile
(interfaceTopModes,
"File",
"Take the modes from the files",
true);
static SwitchOption interfaceTopModesHerwig
(interfaceTopModes,
"Herwig",
"Use the SM ones", false);
static Parameter<BSMModel,double> interfaceBRTolerance
("BRTolerance",
"Tolerance for the sum of branching ratios to be difference from one.",
&BSMModel::tolerance_, 1e-6, 1e-8, 0.01,
false, false, Interface::limited);
}
void BSMModel::doinit() {
StandardModel::doinit();
// check if need to read decays
if(decayFile()=="") return;
// read decays
CFileLineReader cfile;
cfile.open(decayFile_);
if( !cfile ) throw SetupException()
<< "BSMModel::doinit - An error occurred in opening the "
<< "decay file \"" << decayFile_ << "\"."
<< Exception::runerror;
//Before reading the spectrum/decay files the SM higgs
//decay modes, mass and width generators need to be turned off.
PDPtr h0 = getParticleData(ParticleID::h0);
h0->widthGenerator(WidthGeneratorPtr());
h0->massGenerator(MassGenPtr());
h0->width(ZERO);
DecaySet::const_iterator dit = h0->decayModes().begin();
DecaySet::const_iterator dend = h0->decayModes().end();
for( ; dit != dend; ++dit ) {
generator()->preinitInterface(*dit, "BranchingRatio", "set", "0.");
generator()->preinitInterface(*dit, "OnOff", "set", "Off");
}
// if taking the top modes from the file
// delete the SM stuff
if(topModesFromFile_) {
PDPtr top = getParticleData(ParticleID::t);
top->widthGenerator(WidthGeneratorPtr());
top->massGenerator(MassGenPtr());
DecaySet::const_iterator dit = top->decayModes().begin();
DecaySet::const_iterator dend = top->decayModes().end();
for( ; dit != dend; ++dit ) {
generator()->preinitInterface(*dit, "BranchingRatio", "set", "0.");
generator()->preinitInterface(*dit, "OnOff", "set", "Off");
}
}
// read first line and check if this is a Les Houches event file
cfile.readline();
bool lesHouches = cfile.find("<LesHouchesEvents");
bool reading = !lesHouches;
if(lesHouches) cfile.readline();
// function pointer for putting all characters to lower case.
int (*pf)(int) = tolower;
while (true) {
string line = cfile.getline();
// check for start of slha block in SLHA files
if(lesHouches && !reading) {
if(line.find("<slha")==0) reading = true;
if(!cfile.readline()) break;
continue;
}
// ignore comment lines
if(line[0] == '#') {
if(!cfile.readline()) break;
continue;
}
// make everything lower case
transform(line.begin(), line.end(), line.begin(), pf);
// start of a block
if(line.find("decay") == 0) {
readDecay(cfile, line);
+ continue;
}
else if( lesHouches && line.find("</slha") == 0 ) {
reading = false;
break;
}
if(!cfile.readline()) break;
}
}
void BSMModel::readDecay(CFileLineReader & cfile,
string decay) const{
// extract parent PDG code and width
long parent(0);
Energy width(ZERO);
istringstream iss(decay);
string dummy;
iss >> dummy >> parent >> iunit(width, GeV);
PDPtr inpart = getParticleData(parent);
if(!topModesFromFile_&&abs(parent)==ParticleID::t) {
cfile.readline();
return;
}
if(!inpart) throw SetupException()
<< "BSMModel::readDecay() - "
<< "A ParticleData object with the PDG code "
<< parent << " does not exist. "
<< Exception::runerror;
inpart->width(width);
if( width > ZERO ) inpart->cTau(hbarc/width);
inpart->widthCut(5.*width);
Energy inMass = inpart->mass();
- string prefix("decaymode " + inpart->name() + "->");
+ string prefix(inpart->name() + "->");
double brsum(0.);
unsigned int nmode = 0;
while(cfile.readline()) {
string line = cfile.getline();
// skip comments
if(line[0] == '#') continue;
// reached the end
if( line[0] == 'B' || line[0] == 'b' ||
line[0] == 'D' || line[0] == 'd' ||
line[0] == '<' ) {
cfile.resetline();
break;
}
// read the mode
// get the branching ratio and no of decay products
istringstream is(line);
double brat(0.);
unsigned int nda(0),npr(0);
is >> brat >> nda;
vector<tcPDPtr> products,bosons;
Energy mout(ZERO),moutnoWZ(ZERO);
string tag = prefix;
while( true ) {
long t;
is >> t;
if( is.fail() ) break;
if( t == abs(parent) ) {
throw SetupException()
<< "An error occurred while read a decay of the "
<< inpart->PDGName() << ". One of its products has the same PDG code "
<< "as the parent particle. Please check the SLHA file.\n"
<< Exception::runerror;
}
tcPDPtr p = getParticleData(t);
if( !p ) {
throw SetupException()
<< "BSMModel::readDecay() - An unknown PDG code has been encounterd "
<< "while reading a decay mode. ID: " << t
<< Exception::runerror;
}
++npr;
tag += p->name() + ",";
Energy mass = p->mass();
mout += mass;
if(abs(p->id())==ParticleID::Wplus||p->id()==ParticleID::Z0) {
bosons.push_back(p);
}
else {
products.push_back(p);
moutnoWZ += mass;
}
}
if( npr != nda ) {
throw SetupException()
<< "BSMModel::readDecay - While reading a decay of the "
<< inpart->PDGName() << " from an SLHA file, an inconsistency "
<< "between the number of decay products and the value in "
<< "the 'NDA' column was found. Please check if the spectrum "
<< "file is correct.\n"
<< Exception::warning;
}
if( npr > 1 ) {
tag.replace(tag.size() - 1, 1, ";");
// normal option
if(mout<=inMass) {
inpart->stable(false);
brsum += brat;
createDecayMode(tag, brat);
}
// no possible off-shell gauge bosons throw it away
else if(bosons.empty() || bosons.size()>2 ||
moutnoWZ>=inMass) {
cerr << "BSMModel::readDecay() "
<< "The decay " << tag << " cannot proceed for on-shell "
<< "particles, skipping it.\n";
}
else {
Energy maxMass = inMass - moutnoWZ;
string newTag = prefix;
for(unsigned int ix=0;ix<products.size();++ix)
newTag += products[ix]->name() + ",";
if(bosons.size()==1) {
cerr << "BSMModel::readDecay() "
<< "The decay " << tag << " cannot proceed for on-shell\n"
<< "particles, replacing gauge boson with its decay products\n";
vector<pair<double,string> > modes =
createWZDecayModes(newTag,brat,bosons[0],maxMass);
for(unsigned int ix=0;ix<modes.size();++ix) {
modes[ix].second.replace(modes[ix].second.size() - 1, 1, ";");
createDecayMode(modes[ix].second,modes[ix].first);
brsum += modes[ix].first;
}
}
else if(bosons.size()==2) {
bool identical = bosons[0]->id()==bosons[1]->id();
if(maxMass>bosons[0]->mass()&&maxMass>bosons[1]->mass()) {
cerr << "BSMModel::readDecay() "
<< "The decay " << tag << " cannot proceed for on-shell\n"
<< "particles, replacing one of the gauge bosons"
<< " with its decay products\n";
unsigned int imax = identical ? 1 : 2;
if(imax==2) brat *= 0.5;
for(unsigned int ix=0;ix<imax;++ix) {
string newTag2 = newTag+bosons[ix]->name()+',';
unsigned int iother = ix==0 ? 1 : 0;
vector<pair<double,string> > modes =
createWZDecayModes(newTag2,brat,bosons[iother],maxMass);
for(unsigned int ix=0;ix<modes.size();++ix) {
modes[ix].second.replace(modes[ix].second.size() - 1, 1, ";");
createDecayMode(modes[ix].second,modes[ix].first);
brsum += modes[ix].first;
}
}
}
else {
cerr << "BSMModel::readDecay() "
<< "The decay " << tag << " cannot proceed for on-shell\n"
<< "particles, and has too many off-shell gauge bosons,"
<< " skipping it.\n";
}
}
else {
cerr << "BSMModel::readDecay() "
<< "The decay " << tag << " cannot proceed for on-shell\n"
<< "particles, and has too many outgoing gauge bosons skipping it.\n";
}
}
}
}
if( abs(brsum - 1.) > tolerance_ && nmode!=0 ) {
cerr << "Warning: The total branching ratio for " << inpart->PDGName()
<< " from the spectrum file does not sum to 1. The branching fractions"
<< " will be rescaled.\n";
cerr << setprecision(13) << abs(brsum - 1.) << "\n";
}
}
void BSMModel::createDecayMode(string tag, double brat) const {
- ostringstream cmd;
- cmd << tag << string(" ")
- << setprecision(13) << brat << string(" 1 /Herwig/Decays/Mambo");
- Repository::exec(cmd.str(), cerr);
+ tDMPtr dm = generator()->findDecayMode(tag);
+ if(!dm) {
+ dm = generator()->preinitCreateDecayMode(tag);
+ }
+ generator()->preinitInterface(dm, "OnOff", "set", "On");
+ generator()->preinitInterface(dm, "Decayer", "set","/Herwig/Decays/Mambo");
+ ostringstream brf;
+ brf << setprecision(13)<< brat;
+ generator()->preinitInterface(dm, "BranchingRatio","set", brf.str());
}
vector<pair<double,string> >
BSMModel::createWZDecayModes(string tag, double brat,
tcPDPtr boson, Energy maxMass) const {
vector<pair<double,string> > modes;
double sum(0.);
for(DecaySet::const_iterator dit=boson->decayModes().begin();
dit!=boson->decayModes().end();++dit) {
tcDMPtr mode = *dit;
if(!mode->on()) continue;
string extra;
Energy outMass(ZERO);
for(ParticleMSet::const_iterator pit=mode->products().begin();
pit!=mode->products().end();++pit) {
extra += (**pit).name() + ",";
outMass += (**pit).mass();
}
if(outMass<maxMass) {
sum += mode->brat();
modes.push_back(make_pair(mode->brat(),tag+extra));
}
}
for(unsigned int ix=0;ix<modes.size();++ix)
modes[ix].first *= brat/sum;
return modes;
}
diff --git a/Models/General/BSMModel.h b/Models/General/BSMModel.h
--- a/Models/General/BSMModel.h
+++ b/Models/General/BSMModel.h
@@ -1,135 +1,142 @@
// -*- C++ -*-
#ifndef Herwig_BSMModel_H
#define Herwig_BSMModel_H
//
// This is the declaration of the BSMModel class.
//
#include "Herwig++/Models/StandardModel/StandardModel.h"
#include "ThePEG/Utilities/CFileLineReader.h"
namespace Herwig {
using namespace ThePEG;
/**
* Here is the documentation of the BSMModel class.
*
* @see \ref BSMModelInterfaces "The interfaces"
* defined for BSMModel.
*/
class BSMModel: public Herwig::StandardModel {
public:
/**
* The default constructor.
*/
BSMModel();
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/**
* Get name of SLHA decay file
*/
const string & decayFile() const {return decayFile_;}
/**
* Set name of SLHA decay file
*/
void decayFile(string in) {decayFile_ = in;}
/**
* Read decaymodes from LHA file
* @param ifs input stream containg data
* @param decay string containing name of parent and value of total width
*/
void readDecay(CFileLineReader & ifs, string decay) const;
/**
* Create a DecayMode object in the repository
* @param tag The tag identifying the decay mode including the prefix
* 'decaymode'
* @param brat Branching ratio of this mode
*/
void createDecayMode(string tag, double brat) const;
/**
* Create a DecayMode object in the repository
* @param tag The tag identifying the decay mode including the prefix
* 'decaymode'
* @param brat Branching ratio of this mode
*/
vector<pair<double,string> > createWZDecayModes(string tag, double brat,
tcPDPtr boson,
Energy maxMass) const;
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
+ /**
+ * Overloaded function from Interfaced
+ */
+ virtual bool preInitialize() const {
+ return true;
+ }
+
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
BSMModel & operator=(const BSMModel &);
private:
/**
* Name of the decay file
*/
string decayFile_;
/**
* Whether or not to replace the top decay modes with those from
* the SLHA files
*/
bool topModesFromFile_;
/**
* Tolerance for branching ratios
*/
double tolerance_;
};
}
#endif /* Herwig_BSMModel_H */
diff --git a/PDT/SMHiggsMassGenerator.cc b/PDT/SMHiggsMassGenerator.cc
--- a/PDT/SMHiggsMassGenerator.cc
+++ b/PDT/SMHiggsMassGenerator.cc
@@ -1,84 +1,86 @@
// -*- C++ -*-
//
// SMHiggsMassGenerator.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SMHiggsMassGenerator class.
//
#include "SMHiggsMassGenerator.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
void SMHiggsMassGenerator::persistentOutput(PersistentOStream & os) const {
os << _shape << _hwidth;
}
void SMHiggsMassGenerator::persistentInput(PersistentIStream & is, int) {
is >> _shape >> _hwidth;
}
ClassDescription<SMHiggsMassGenerator> SMHiggsMassGenerator::initSMHiggsMassGenerator;
// Definition of the static class description member.
void SMHiggsMassGenerator::Init() {
static ClassDocumentation<SMHiggsMassGenerator> documentation
("The SMHiggsMassGenerator class implements the mass distribution for the"
" Higgs boson as in hep-ph/9505211.",
"The Higgs mass was distributed as in \\cite{Seymour:1995qg}.",
"%\\cite{Seymour:1995qg}\n"
"\\bibitem{Seymour:1995qg}\n"
" M.~H.~Seymour,\n"
" %``The Higgs boson line shape and perturbative unitarity,''\n"
" Phys.\\ Lett.\\ B {\\bf 354}, 409 (1995)\n"
" [arXiv:hep-ph/9505211].\n"
" %%CITATION = PHLTA,B354,409;%%\n"
);
static Switch<SMHiggsMassGenerator,unsigned int> interfaceHiggsShape
("HiggsShape",
"The distribution for the Higgs mass",
&SMHiggsMassGenerator::_shape, 1, false, false);
static SwitchOption interfaceHiggsShapeNormal
(interfaceHiggsShape,
"Normal",
"The standard running width distribution",
0);
static SwitchOption interfaceHiggsShapeImproved
(interfaceHiggsShape,
"Improved",
"The improved shape of hep-ph/9505211",
1);
}
bool SMHiggsMassGenerator::accept(const ParticleData & part) const {
if(part.id()!=ParticleID::h0) return false;
return GenericMassGenerator::accept(part);
}
void SMHiggsMassGenerator::doinit() {
- if(particle()->widthGenerator()) {
- _hwidth=dynamic_ptr_cast<GenericWidthGeneratorPtr>(particle()->widthGenerator());
+ if(particle()->massGenerator()==this) {
+ if(particle()->widthGenerator()) {
+ _hwidth=dynamic_ptr_cast<GenericWidthGeneratorPtr>(particle()->widthGenerator());
+ }
+ if(!_hwidth) throw InitException()
+ << "Must be using the Herwig::GenericWidthGenerator in "
+ << "SMHiggsMassGenerator::doinit()" << Exception::runerror;
}
- if(!_hwidth) throw InitException()
- << "Must be using the Herwig::GenericWidthGenerator in "
- << "SMHiggsMassGenerator::doinit()" << Exception::runerror;
GenericMassGenerator::doinit();
}
void SMHiggsMassGenerator::dataBaseOutput(ofstream & output,bool header) {
if(header) output << "update Mass_Generators set parameters=\"";
output << "newdef " << fullName() << ":BreitWignerShape " << _shape << "\n";
if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl;
}

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 8:42 PM (22 h, 38 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242328
Default Alt Text
(19 KB)

Event Timeline