No OneTemporary
Actions

Size

22 KB

Subscribers

None

View Options

	diff --git a/LesHouches/LesHouchesFileReader.cc b/LesHouches/LesHouchesFileReader.cc
	--- a/LesHouches/LesHouchesFileReader.cc
	+++ b/LesHouches/LesHouchesFileReader.cc
	@@ -1,665 +1,666 @@
	// -- C++ --
	//
	// LesHouchesFileReader.cc is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2011 Leif Lonnblad
	//
	// ThePEG 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 LesHouchesFileReader class.
	//

	#include "LesHouchesFileReader.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Utilities/Throw.h"
	#include "ThePEG/PDT/DecayMode.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"

	using namespace ThePEG;

	LesHouchesFileReader::
	LesHouchesFileReader(const LesHouchesFileReader & x)
	: LesHouchesReader(x), neve(x.neve), ieve(0),
	LHFVersion(x.LHFVersion), outsideBlock(x.outsideBlock),
	headerBlock(x.headerBlock), initComments(x.initComments),
	initAttributes(x.initAttributes), eventComments(x.eventComments),
	eventAttributes(x.eventAttributes),
	theFileName(x.theFileName), theQNumbers(x.theQNumbers),
	theDecayer(x.theDecayer) {}

	LesHouchesFileReader::~LesHouchesFileReader() {}

	IBPtr LesHouchesFileReader::clone() const {
	return new_ptr(*this);
	}

	IBPtr LesHouchesFileReader::fullclone() const {
	return new_ptr(*this);
	}

	bool LesHouchesFileReader::preInitialize() const {
	return true;
	}

	void LesHouchesFileReader::doinit() {
	LesHouchesReader::doinit();
	// are we using QNUMBERS
	if(!theQNumbers) return;
	// parse the header block and create
	// any new particles needed in QNUMBERS blocks
	string block = headerBlock;
	string line = "";
	bool readingSLHA = false;
	int (*pf)(int) = tolower;
	do {
	line = StringUtils::car(block,"\r\n");
	block = StringUtils::cdr(block,"\r\n");
	if(line[0]=='#') continue;
	// are we reading the SLHA block
	if(readingSLHA) {
	// reached the end of slha block ?
	if(line.find("</slha") != string::npos) {
	readingSLHA = false;
	break;
	}
	// remove trailing comment from line
	vector<string> split = StringUtils::split(line,"#");
	// check for a qnumbers block
	transform(split[0].begin(), split[0].end(), split[0].begin(), pf);
	// if not contine
	if(split[0].find("block qnumbers")==string::npos)
	continue;
	// get name from comment
	string name = StringUtils::stripws(split[1]);
	// extract the PDG code
	split = StringUtils::split(split[0]," ");
	istringstream is(split[2]);
	long PDGCode(0);
	is >> PDGCode;
	// get the charge, spin, colour and whether an antiparticle
	int charge(0),spin(0),colour(0),anti(0);
	for(unsigned int ix=0;ix<4;++ix) {
	line = StringUtils::car(block,"\r\n");
	block = StringUtils::cdr(block,"\r\n");
	int dummy[2];
	istringstream is(line);
	is >> dummy[0] >> dummy[1];
	switch (dummy[0]) {
	case 1:
	charge = dummy[1];
	break;
	case 2:
	spin = dummy[1];
	break;
	case 3:
	colour = dummy[1];
	break;
	case 4:
	anti = dummy[1];
	break;
	default:
	assert(false);
	}
	}
	// check if particles already exist
	PDPair newParticle;
	newParticle.first = getParticleData(PDGCode);
	if(newParticle.first) Throw<SetupException>()
	<< "Particle with PDG code " << PDGCode
	<< " whose creation was requested in a QNUMBERS Block"
	<< " already exists. Retaining the original particle"
	<< Exception::warning;
	if(anti) {
	newParticle.second = getParticleData(-PDGCode);
	if(newParticle.second) Throw<SetupException>()
	<< "Anti-particle with PDG code " << -PDGCode
	<< " whose creation was requested in a QNUMBERS Block"
	<< " already exists. Retaining the original particle"
	<< Exception::warning;
	if(( newParticle.first && !newParticle.second ) \|\|
	( newParticle.second && !newParticle.first ) )
	Throw<SetupException>()
	<< "Either particle or anti-particle with PDG code " << PDGCode
	<< " whose creation was requested in a QNUMBERS Block"
	<< " already exists, but not both the particle and antiparticle. "
	<< " Something dodgy here stopping"
	<< Exception::runerror;
	}
	// already exists continue
	if(newParticle.first) continue;
	// create the particles
	// particle with no anti particle
	if( anti == 0 ) {
	// construct the name
	if(name=="") {
	ostringstream temp;
	temp << PDGCode;
	name = temp.str();
	}
	// create the ParticleData object
	newParticle.first = ParticleData::Create(PDGCode,name);
	}
	// particle anti-particle pair
	else {
	// construct the names
	string nameAnti;
	if(name=="") {
	ostringstream temp;
	temp << PDGCode;
	name = temp.str();
	ostringstream temp2;
	temp << -PDGCode;
	nameAnti = temp2.str();
	}
	else {
	nameAnti=name;
	for(string::iterator it=nameAnti.begin();it!=nameAnti.end();++it) {
	if(*it=='+') nameAnti.replace(it,it+1,"-");
	else if(*it=='-') nameAnti.replace(it,it+1,"+");
	}
	if(nameAnti==name) nameAnti += "bar";
	}
	// create the ParticleData objects
	newParticle = ParticleData::Create(PDGCode,name,nameAnti);
	}
	// set the particle properties
	if(colour==1) colour = 0;
	newParticle.first->iColour(PDT::Colour(colour));
	newParticle.first->iSpin (PDT::Spin (spin ));
	newParticle.first->iCharge(PDT::Charge(charge));
	// register it
	generator()->preinitRegister(newParticle.first,
	"/Herwig/Particles/"+newParticle.first->PDGName());
	// set the antiparticle properties
	if(newParticle.second) {
	if(colour==3\|\|colour==6) colour *= -1;
	charge = -charge;
	newParticle.second->iColour(PDT::Colour(colour));
	newParticle.second->iSpin (PDT::Spin (spin ));
	newParticle.second->iCharge(PDT::Charge(charge));
	// register it
	generator()->preinitRegister(newParticle.second,
	"/Herwig/Particles/"+newParticle.second->PDGName());
	}
	}
	// start of SLHA block ?
	else if(line.find("<slha") != string::npos) {
	readingSLHA = true;
	}
	}
	while(line!="");
	// now set any masses/decay modes
	block = headerBlock;
	line="";
	readingSLHA=false;
	bool ok=true;
	do {
	line = StringUtils::car(block,"\r\n");
	block = StringUtils::cdr(block,"\r\n");
	// are we reading the SLHA block
	if(readingSLHA) {
	// reached the end?
	if(line.find("</slha") == 0 ) {
	readingSLHA = false;
	break;
	}
	// make lower case
	transform(line.begin(),line.end(),line.begin(), pf);
	// found the mass block ?
	if(line.find("block mass")!=string::npos) {
	// read it
	line = StringUtils::car(block,"\r\n");
	// check not at end
	while(line[0] != 'D' && line[0] != 'B' &&
	line[0] != 'd' && line[0] != 'b' &&
	line != "") {
	// skip comment lines
	if(line[0] == '#') {
	block = StringUtils::cdr(block,"\r\n");
	line = StringUtils::car(block,"\r\n");
	continue;
	}
	// get the mass and PGD code
	istringstream temp(line);
	long id;
	double mass;
	temp >> id >> mass;
	// skip resetting masses on SM particles
	// as it can cause problems later on in event generation
	if(abs(id)<=6 \|\| (abs(id)>=11 && abs(id)<=16) \|\|
	abs(id)==23 \|\| abs(id)==24) {
	block = StringUtils::cdr(block,"\r\n");
	line = StringUtils::car(block,"\r\n");
	continue;
	}
	// magnitude of mass for susy models
	mass = abs(mass);
	// set the mass
	tPDPtr particle = getParticleData(id);
	if(!particle) throw SetupException()
	<< "LesHouchesFileReader::doinit() - Particle with PDG code not"
	<< id << " not found." << Exception::runerror;
	const InterfaceBase * ifb = BaseRepository::FindInterface(particle,
	"NominalMass");
	ostringstream os;
	os << mass;
	ifb->exec(*particle, "set", os.str());
	// read the next line
	block = StringUtils::cdr(block,"\r\n");
	line = StringUtils::car(block,"\r\n");
	};
	}
	// found a decay block
	else if(line.find("decay") == 0) {
	// get PGD code and width
	istringstream iss(line);
	string dummy;
	long parent(0);
	Energy width(ZERO);
	iss >> dummy >> parent >> iunit(width, GeV);
	// get the ParticleData object
	PDPtr inpart = getParticleData(parent);
	if(!inpart) {
	throw SetupException()
	<< "LesHouchesFileReader::doinit() - A ParticleData object with the PDG code "
	<< parent << " does not exist. " << Exception::runerror;
	return;
	}
	// set the width
	inpart->width(width);
	if( width > ZERO ) {
	inpart->cTau(hbarc/width);
	inpart->widthCut(5.*width);
	inpart->stable(false);
	}
	// construct prefix for DecayModes
	string prefix(inpart->name() + "->"), tag(prefix),line("");
	unsigned int nmode(0);
	// read any decay modes
	line = StringUtils::car(block,"\r\n");
	while(line[0] != 'D' && line[0] != 'B' &&
	line[0] != 'd' && line[0] != 'b' &&
	line != "") {
	// skip comments
	if(line[0] == '#') {
	block = StringUtils::cdr(block,"\r\n");
	line = StringUtils::car(block,"\r\n");
	continue;
	}
	// read decay mode and construct the tag
	istringstream is(line);
	double brat(0.);
	unsigned int nda(0),npr(0);
	is >> brat >> nda;
	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 in LesHouchesFileReader::doinit()."
	<< " Please check the Les Houches file.\n"
	<< Exception::runerror;
	tcPDPtr p = getParticleData(t);
	if( !p )
	throw SetupException()
	<< "LesHouchesFileReader::doinit() -"
	<< " An unknown PDG code has been encounterd "
	<< "while reading a decay mode. ID: " << t
	<< Exception::runerror;
	++npr;
	tag += p->name() + ",";
	}
	if( npr != nda )
	throw SetupException()
	<< "LesHouchesFileReader::doinit() - 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;
	// create the DecayMode
	if( npr > 1 ) {
	if( nmode==0 ) {
	generator()->preinitInterface(inpart, "VariableRatio" , "set","false");
	if(inpart->massGenerator()) {
	ok = false;
	Throw<SetupException>()
	<< inpart->PDGName() << " already has a MassGenerator set"
	<< " this is incompatible with using QNUMBERS "
	<< "Use\n"
	<< "set " << inpart->fullName() << ":Mass_generator NULL\n"
	<< "to fix this." << Exception::warning;
	}
	if(inpart->widthGenerator()) {
	ok = false;
	Throw<SetupException>()
	<< inpart->PDGName() << " already has a WidthGenerator set"
	<< " this is incompatible with using QNUMBERS "
	<< "Use\n"
	<< "set " << inpart->fullName() << ":Width_generator NULL\n"
	<< "to fix this." << Exception::warning;
	}
	unsigned int ntemp=0;
	for(DecaySet::const_iterator dit = inpart->decayModes().begin();
	dit != inpart->decayModes().end(); ++dit ) {
	if((**dit).on()) ++ntemp;
	}
	if(ntemp!=0) {
	ok = false;
	Throw<SetupException>()
	<< inpart->PDGName() << " already has DecayModes"
	<< " this is incompatible with using QNUMBERS "
	<< "Use\n"
	<< "do " << inpart->fullName() << ":SelectDecayModes none\n"
	<< " to fix this." << Exception::warning;
	}
	}
	inpart->stable(false);
	tag.replace(tag.size() - 1, 1, ";");
	DMPtr dm = generator()->findDecayMode(tag);
	+ if(!theDecayer) Throw<SetupException>()
	+ << "LesHouchesFileReader::doinit() Decayer must be set using the "
	+ << "LesHouchesFileReader:Decayer"
	+ << " must be set to allow the creation of new"
	+ << " decay modes."
	+ << Exception::runerror;
	if(!dm) {
	- if(!theDecayer) Throw<SetupException>()
	- << "LesHouchesFileReader::doinit() Decayer"
	- << " must be set to allow the creation of new"
	- << " decay modes"
	- << Exception::runerror;
	dm = generator()->preinitCreateDecayMode(tag);
	if(!dm)
	Throw<SetupException>()
	<< "LesHouchesFileReader::doinit() - Needed to create "
	<< "new decaymode but one could not be created for the tag "
	<< tag << Exception::warning;
	}
	generator()->preinitInterface(dm, "Decayer", "set",
	theDecayer->fullName());
	ostringstream br;
	br << setprecision(13) << brat;
	generator()->preinitInterface(dm, "BranchingRatio", "set", br.str());
	generator()->preinitInterface(dm, "OnOff", "set", "On");
	if(dm->CC()) {
	generator()->preinitInterface(dm->CC(), "BranchingRatio", "set", br.str());
	generator()->preinitInterface(dm->CC(), "OnOff", "set", "On");
	}
	++nmode;
	}
	tag=prefix;
	// read the next line
	block = StringUtils::cdr(block,"\r\n");
	line = StringUtils::car(block,"\r\n");
	};
	if(nmode>0) {
	inpart->update();
	if(inpart->CC())
	inpart->CC()->update();
	}
	}
	}
	// start of SLHA block ?
	else if(line.find("<slha") != string::npos) {
	readingSLHA = true;
	}
	}
	while(line!="");
	if(!ok)
	throw SetupException() << "Problem reading QNUMBERS blocks in LesHouchesFileReader::doinit()"
	<< Exception::runerror;
	}

	void LesHouchesFileReader::initialize(LesHouchesEventHandler & eh) {
	LesHouchesReader::initialize(eh);
	if ( LHFVersion.empty() )
	Throw<LesHouchesFileError>()
	<< "The file associated with '" << name() << "' does not contain a "
	<< "proper formatted Les Houches event file. The events may not be "
	<< "properly sampled." << Exception::warning;
	}

	void LesHouchesFileReader::open() {
	if ( filename().empty() )
	throw LesHouchesFileError()
	<< "No Les Houches file name. "
	<< "Use 'set " << name() << ":FileName'."
	<< Exception::runerror;
	cfile.open(filename());
	if ( !cfile )
	throw LesHouchesFileError()
	<< "The LesHouchesFileReader '" << name() << "' could not open the "
	<< "event file called '" << theFileName << "'."
	<< Exception::runerror;

	cfile.readline();
	if ( !cfile.find("<LesHouchesEvents") ) return;
	map<string,string> attributes =
	StringUtils::xmlAttributes("LesHouchesEvents", cfile.getline());
	LHFVersion = attributes["version"];
	if ( LHFVersion.empty() ) return;

	bool readingHeader = false;
	bool readingInit = false;
	headerBlock = "";

	// Loop over all lines until we hit the </init> tag.
	while ( cfile.readline() && !cfile.find("</init>") ) {
	if ( cfile.find("<header") ) {
	// We have hit the header block, so we should dump this and all
	// following lines to headerBlock until we hit the end of it.
	readingHeader = true;
	headerBlock = cfile.getline() + "\n";
	}
	else if ( cfile.find("<init") ) {
	// We have hit the init block, so we should expect to find the
	// standard information in the following. But first check for
	// attributes.
	initAttributes = StringUtils::xmlAttributes("init", cfile.getline());
	readingInit = true;
	cfile.readline();
	if ( !( cfile >> heprup.IDBMUP.first >> heprup.IDBMUP.second
	>> heprup.EBMUP.first >> heprup.EBMUP.second
	>> heprup.PDFGUP.first >> heprup.PDFGUP.second
	>> heprup.PDFSUP.first >> heprup.PDFSUP.second
	>> heprup.IDWTUP >> heprup.NPRUP ) ) {
	heprup.NPRUP = -42;
	LHFVersion = "";
	return;
	}
	heprup.resize();

	for ( int i = 0; i < heprup.NPRUP; ++i ) {
	cfile.readline();
	if ( !( cfile >> heprup.XSECUP[i] >> heprup.XERRUP[i]
	>> heprup.XMAXUP[i] >> heprup.LPRUP[i] ) ) {
	heprup.NPRUP = -42;
	LHFVersion = "";
	return;
	}
	}
	}
	else if ( cfile.find("</header") ) {
	readingHeader = false;
	headerBlock += cfile.getline() + "\n";
	}
	else if ( readingHeader ) {
	// We are in the process of reading the header block. Dump the
	// line to headerBlock.
	headerBlock += cfile.getline() + "\n";
	}
	else if ( readingInit ) {
	// Here we found a comment line. Dump it to initComments.
	initComments += cfile.getline() + "\n";
	}
	else {
	// We found some other stuff outside the standard tags.
	outsideBlock += cfile.getline() + "\n";
	}
	}
	if ( !cfile ) {
	heprup.NPRUP = -42;
	LHFVersion = "";
	return;
	}

	}

	bool LesHouchesFileReader::doReadEvent() {
	if ( !cfile ) return false;
	if ( LHFVersion.empty() ) return false;
	if ( heprup.NPRUP < 0 ) return false;
	eventComments = "";
	outsideBlock = "";
	hepeup.NUP = 0;
	hepeup.XPDWUP.first = hepeup.XPDWUP.second = 0.0;

	// Keep reading lines until we hit the next event or the end of
	// the event block. Save any inbetween lines. Exit if we didn't
	// find an event.
	while ( cfile.readline() && !cfile.find("<event") )
	outsideBlock += cfile.getline() + "\n";

	// We found an event. First scan for attributes.
	eventAttributes = StringUtils::xmlAttributes("event", cfile.getline());
	if ( !cfile.readline() ) return false;

	// The first line determines how many subsequent particle lines we
	// have.
	if ( !( cfile >> hepeup.NUP >> hepeup.IDPRUP >> hepeup.XWGTUP
	>> hepeup.SCALUP >> hepeup.AQEDUP >> hepeup.AQCDUP ) )
	return false;
	hepeup.resize();
	// Read all particle lines.
	for ( int i = 0; i < hepeup.NUP; ++i ) {
	if ( !cfile.readline() ) return false;
	if ( !( cfile >> hepeup.IDUP[i] >> hepeup.ISTUP[i]
	>> hepeup.MOTHUP[i].first >> hepeup.MOTHUP[i].second
	>> hepeup.ICOLUP[i].first >> hepeup.ICOLUP[i].second
	>> hepeup.PUP[i][0] >> hepeup.PUP[i][1] >> hepeup.PUP[i][2]
	>> hepeup.PUP[i][3] >> hepeup.PUP[i][4]
	>> hepeup.VTIMUP[i] >> hepeup.SPINUP[i] ) )
	return false;
	if(isnan(hepeup.PUP[i][0])\|\|isnan(hepeup.PUP[i][1])\|\|
	isnan(hepeup.PUP[i][2])\|\|isnan(hepeup.PUP[i][3])\|\|
	isnan(hepeup.PUP[i][4]))
	throw Exception()
	<< "nan's as momenta in Les Houches file "
	<< Exception::eventerror;
	if(hepeup.MOTHUP[i].first -1==i \|\|
	hepeup.MOTHUP[i].second-1==i) {
	throw Exception()
	<< "Particle has itself as a mother in Les Houches "
	<< "file, this is not allowed\n"
	<< Exception::eventerror;
	}
	}

	// Now read any additional comments and named weights.
	optionalWeights.clear();
	while ( cfile.readline() && !cfile.find("</event>") ) {
	// parse weightinfo elements
	if ( cfile.find("<weightinfo") ) {
	map<string,string> weightAttributes =
	StringUtils::xmlAttributes("weightinfo",cfile.getline());
	if ( !cfile.readline() ) return false;
	double weightValue; cfile >> weightValue;
	if ( !cfile.readline() ) return false;
	if ( !cfile.find("</weightinfo>") ) {
	throw LesHouchesFileError()
	<< "error while parsing weight information in "
	<< theFileName << Exception::eventerror;
	}
	// ignore the single unnamed weight, but warn if multiple of the
	// same name
	map<string,string>::const_iterator weightName = weightAttributes.find("name");
	if ( weightName != weightAttributes.end() ) {
	map<string,double>::iterator wit =
	optionalWeights.find(weightName->second);
	if ( wit == optionalWeights.end() ) {
	optionalWeights[weightName->second] = weightValue;
	} else {
	throw LesHouchesFileError()
	<< "multiple weights of the same name in "
	<< theFileName << Exception::eventerror;
	}
	}
	continue;
	}
	eventComments += cfile.getline() + "\n";
	}

	if ( !cfile ) return false;
	return true;

	}

	void LesHouchesFileReader::close() {
	cfile.close();
	}

	void LesHouchesFileReader::persistentOutput(PersistentOStream & os) const {
	os << neve << LHFVersion << outsideBlock << headerBlock << initComments
	<< initAttributes << eventComments << eventAttributes << theFileName
	<< theQNumbers << theDecayer;
	}

	void LesHouchesFileReader::persistentInput(PersistentIStream & is, int) {
	is >> neve >> LHFVersion >> outsideBlock >> headerBlock >> initComments
	>> initAttributes >> eventComments >> eventAttributes >> theFileName
	>> theQNumbers >> theDecayer;
	ieve = 0;
	}

	ClassDescription<LesHouchesFileReader>
	LesHouchesFileReader::initLesHouchesFileReader;
	// Definition of the static class description member.

	void LesHouchesFileReader::Init() {

	static ClassDocumentation<LesHouchesFileReader> documentation
	("ThePEG::LesHouchesFileReader is an base class to be used for objects "
	"which reads event files from matrix element generators. This class is "
	"able to read plain event files conforming to the Les Houches Event File "
	"accord.");

	static Parameter<LesHouchesFileReader,string> interfaceFileName
	("FileName",
	"The name of a file containing events conforming to the Les Houches "
	"protocol to be read into ThePEG. A file name ending in "
	"<code>.gz</code> will be read from a pipe which uses "
	"<code>zcat</code>. If a file name ends in <code>\|</code> the "
	"preceeding string is interpreted as a command, the output of which "
	"will be read through a pipe.",
	&LesHouchesFileReader::theFileName, "", false, false);

	interfaceFileName.fileType();
	interfaceFileName.rank(11);

	static Switch<LesHouchesFileReader,bool> interfaceQNumbers
	("QNumbers",
	"Whether or not to read search for and read a QNUMBERS"
	" block in the header of the file.",
	&LesHouchesFileReader::theQNumbers, false, false, false);
	static SwitchOption interfaceQNumbersYes
	(interfaceQNumbers,
	"Yes",
	"Use QNUMBERS",
	true);
	static SwitchOption interfaceQNumbersNo
	(interfaceQNumbers,
	"No",
	"Don't use QNUMBERS",
	false);

	static Reference<LesHouchesFileReader,Decayer> interfaceDecayer
	("Decayer",
	"Decayer to use for any decays read from the QNUMBERS Blocks",
	&LesHouchesFileReader::theDecayer, false, false, true, true, false);

	}

File Metadata

Mime Type: text/x-diff
Expires: Mon, Jan 20, 9:13 PM (1 d, 1 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4214351
Default Alt Text: (22 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions