Page MenuHomeHEPForge
Files F8309029

No OneTemporary
Actions

View Options

diff --git a/LesHouches/LesHouchesFileReader.cc b/LesHouches/LesHouchesFileReader.cc
--- a/LesHouches/LesHouchesFileReader.cc
+++ b/LesHouches/LesHouchesFileReader.cc
@@ -1,873 +1,874 @@
// -*- 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"
#include <boost/algorithm/string.hpp>
#include <boost/lexical_cast.hpp>
#include <sstream>
#include <iostream>
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();
//cout << "theDecayer->fullName() = " << theDecayer->fullName() << endl;
//if(theDecayer) {
// cout << "DECAYER RETURNS TRUE" << endl;
// }
// 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;
unsigned int newNumber(0);
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;
if(split.size()>=2) {
name = StringUtils::stripws(split[1]);
}
else {
++newNumber;
ostringstream tname;
tname << "NP" << newNumber;
name = tname.str();
}
// 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) {
// Throw<SetupException>() << "Standard model mass for PID "
// << id
// << " will not be changed."
// << Exception::warning;
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;
}
if ( abs(inpart->id()) == 6 ||
abs(inpart->id()) == 15 ||
abs(inpart->id()) == 23 ||
abs(inpart->id()) == 24 ||
abs(inpart->id()) == 25 ) {
Throw<SetupException>() << "\n"
"************************************************************************\n"
"* Your LHE file changes the width of " << inpart->PDGName() << ".\n"
"* This can cause serious problems in the event generation!\n"
"************************************************************************\n"
"\n" << Exception::warning;
}
else if (inpart->width() > ZERO && width <= ZERO) {
Throw<SetupException>() << "\n"
"************************************************************************\n"
"* Your LHE file zeroes the non-zero width of " << inpart->PDGName() << ".\n"
"* If " << inpart->PDGName() << " is a decaying SM particle,\n"
"* this can cause serious problems in the event generation!\n"
"************************************************************************\n"
"\n" << Exception::warning;
}
// 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[0] != '<' && 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) {
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;
}
//vector<string> LesHouchesFileReader::optWeightNamesFunc() { return optionalWeightsNames; }
vector<string> LesHouchesFileReader::optWeightsNamesFunc() { return optionalWeightsNames; }
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"];
//cout << LHFVersion << endl;
if ( LHFVersion.empty() ) return;
bool readingHeader = false;
bool readingInit = false;
headerBlock = "";
char (cwgtinfo_weights_info[250][15]);
string hs;
int cwgtinfo_nn(0);
// Loop over all lines until we hit the </init> tag.
bool readingInitWeights = false, readingInitWeights_sc = false;
string weightinfo;
while ( cfile.readline() ) {
// found the init block for multiple weights
if(cfile.find("<initrwgt")) { /*cout << "reading weights" << endl;*/ readingInitWeights = true; }
// found end of init block for multiple weights: end the while loop
if(cfile.find("</initrwgt")) { readingInitWeights = false; readingInitWeights_sc = false; continue;}
// found the end of init block
if(cfile.find("</init")) { readingInit = false; break; }
/* read the weight information block
* optionalWeightNames will contain information about the weights
* this will be appended to the weight information
*/
if(readingInitWeights) {
string scalename = "";
if(cfile.find("<weightgroup")) {
/* we found a weight group
* start reading the information
* within it
*/
readingInitWeights_sc = true;
weightinfo = cfile.getline();
/* to make it shorter, erase some stuff
*/
boost::erase_all(weightinfo, "<weightgroup");
boost::erase_all(weightinfo, ">");
boost::erase_all(weightinfo, "\n");
}
/* if we are reading a new weightgroup, go on
* until we find the end of it
*/
if(readingInitWeights_sc && !cfile.find("</weightgroup")) {
// cout << "weightinfo = " << weightinfo << endl;
hs = cfile.getline();
istringstream isc(hs);
int ws = 0;
/* get the name that will be used to identify the scale
*/
do {
string sub; isc >> sub;
if(ws==1) { boost::erase_all(sub, ">"); scalename = sub; }
++ws;
} while (isc);
/* now get the relevant information
* e.g. scales or PDF sets used
*/
string startDEL = "'>"; //starting delimiter
string stopDEL = "</weight>"; //end delimiter
unsigned firstLim = hs.find(startDEL); //find start of delimiter
unsigned lastLim = hs.find(stopDEL); //find end of delimitr
string scinfo = hs.substr(firstLim); //define the information for the scale
boost::erase_all(scinfo,stopDEL);
boost::erase_all(scinfo,startDEL);
scinfo = weightinfo + scinfo;
/* fill in the map
* indicating the information to be appended to each scale
* i.e. scinfo for each scalname
*/
scalemap[scalename] = scinfo.c_str();
boost::erase_all(scalename, "id=");
boost::erase_all(scalename, "'");
optionalWeightsNames.push_back(scalename);
}
}
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";
}
if ( (cfile.find("<init ") && !cfile.find("<initrwgt")) || cfile.find("<init>") ) {
//cout << "found init block" << endl;
// 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;
}
}
}
if ( cfile.find("</header") ) {
readingHeader = false;
headerBlock += cfile.getline() + "\n";
}
if ( readingHeader ) {
/* We are in the process of reading the header block. Dump the
line to headerBlock.*/
headerBlock += cfile.getline() + "\n";
}
if ( readingInit ) {
// Here we found a comment line. Dump it to initComments.
initComments += cfile.getline() + "\n";
}
}
string central = "central";
optionalWeightsNames.push_back(central);
// cout << "reading init finished" << endl;
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;
optionalWeights.clear();
optionalWeightsTemp.clear();
// 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());
/* information necessary for FxFx merging:
* the npLO and npNLO tags
*/
istringstream ievat(cfile.getline());
- int we(0), npLO(-10), npNLO(-10);
+ int we(0), npLO(-99), npNLO(-99);
do {
string sub; ievat >> sub;
if(we==2) { npLO = atoi(sub.c_str()); }
if(we==5) { npNLO = atoi(sub.c_str()); }
++we;
} while (ievat);
optionalnpLO = npLO;
optionalnpNLO = npNLO;
std::stringstream npstringstream;
npstringstream << "np " << npLO << " " << npNLO;
std::string npstrings = npstringstream.str();
/* the FxFx merging information
* becomes part of the optionalWeights, labelled -999
* for future reference
*/
optionalWeights[npstrings.c_str()] = -999;
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;
//print momenta for debugging
// cout << hepeup.PUP[i][0] << " " << hepeup.PUP[i][1] << " " << hepeup.PUP[i][2] << " " << hepeup.PUP[i][3] << " " << hepeup.PUP[i][4] << endl;
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.
// read until the end of rwgt is found
bool readingWeights = false, readingaMCFast = false, readingMG5ClusInfo = false;
while ( cfile.readline() && !cfile.find("</event>")) {
if(cfile.find("</applgrid")) { readingaMCFast=false; } //aMCFast weights end
if(cfile.find("</rwgt")) { readingWeights=false; } //optional weights end
if(cfile.find("</clustering")) { readingMG5ClusInfo=false; } // mg5 mclustering scale info end
/* reading of optional weights
*/
if(readingWeights) {
if(!cfile.find("<wgt")) { continue; }
istringstream iss(cfile.getline());
int wi = 0;
double weightValue(0);
string weightName = "";
// we need to put the actual weight value into a double
do {
string sub; iss >> sub;
if(wi==1) { boost::erase_all(sub, ">"); weightName = sub; }
if(wi==2) weightValue = atof(sub.c_str());
++wi;
} while (iss);
// store the optional weights found in the temporary map
optionalWeightsTemp[weightName] = weightValue;
}
/* reading of aMCFast weights
*/
if(readingaMCFast) {
std::stringstream amcfstringstream;
amcfstringstream << "aMCFast " << cfile.getline();
std::string amcfstrings = amcfstringstream.str();
+ boost::erase_all(amcfstrings,"\n");
optionalWeights[amcfstrings.c_str()] = -111; //for the aMCFast we give them a weight -111 for future identification
}
/* read additional MG5 Clustering information
* used in LO merging
*/
if(readingMG5ClusInfo) {
string hs = cfile.getline();
string startDEL = "<clus scale="; //starting delimiter
string stopDEL = "</clus>"; //end delimiter
unsigned firstLim = hs.find(startDEL); //find start of delimiter
// unsigned lastLim = hs.find(stopDEL); //find end of delimitr
string mg5clusinfo = hs.substr(firstLim); //define the information for the scale
boost::erase_all(mg5clusinfo,stopDEL);
boost::erase_all(mg5clusinfo,startDEL);
boost::erase_all(mg5clusinfo,">");
boost::erase_all(mg5clusinfo,"\"");
optionalWeights[mg5clusinfo.c_str()] = -222; //for the mg5 scale info weights we give them a weight -222 for future identification
}
//store MG5 clustering information
if(cfile.find("<scales")) {
string hs = cfile.getline();
string startDEL = "<scales"; //starting delimiter
string stopDEL = "</scales>"; //end delimiter
unsigned firstLim = hs.find(startDEL); //find start of delimiter
// unsigned lastLim = hs.find(stopDEL); //find end of delimitr
string mg5scinfo = hs.substr(firstLim); //define the information for the scale
boost::erase_all(mg5scinfo,stopDEL);
boost::erase_all(mg5scinfo,startDEL);
boost::erase_all(mg5scinfo,">");
optionalWeights[mg5scinfo.c_str()] = -333; //for the mg5 scale info weights we give them a weight -333 for future identification
}
//determine start of aMCFast weights
if(cfile.find("<applgrid")) { readingaMCFast = true;}
//determine start of optional weights
if(cfile.find("<rwgt")) { readingWeights = true; }
//determine start of MG5 clustering scale information
if(cfile.find("<clustering")) { readingMG5ClusInfo = true;}
}
// loop over the optional weights and add the extra information as found in the init
for (map<string,double>::const_iterator it=optionalWeightsTemp.begin(); it!=optionalWeightsTemp.end(); ++it){
for (map<string,string>::const_iterator it2=scalemap.begin(); it2!=scalemap.end(); ++it2){
//find the scale id in the scale information and add this information
if(it->first==it2->first) {
string info = it2->second + " " + it->first;
boost::erase_all(info, "\n");
//set the optional weights
optionalWeights[info] = it->second;
/*cout << "info = " << info << endl;
std::cout << it->first << " => " << it->second << '\n';*/
}
}
}
/* additionally, we set the "central" scale
* this is actually the default event weight
*/
string central = "central";
optionalWeights[central] = hepeup.XWGTUP;
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;
+ << 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);
}
diff --git a/LesHouches/LesHouchesReader.cc b/LesHouches/LesHouchesReader.cc
--- a/LesHouches/LesHouchesReader.cc
+++ b/LesHouches/LesHouchesReader.cc
@@ -1,1566 +1,1566 @@
// -*- C++ -*-
//
// LesHouchesReader.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 LesHouchesReader class.
//
#include "LesHouchesReader.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Command.h"
#include "config.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDF/PartonExtractor.h"
#include "ThePEG/PDF/NoPDF.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/EventRecord/TmpTransform.h"
#include "ThePEG/Utilities/UtilityBase.h"
#include "ThePEG/Handlers/XComb.h"
#include "ThePEG/Handlers/CascadeHandler.h"
#include "LesHouchesEventHandler.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Utilities/HoldFlag.h"
#include "ThePEG/Utilities/Debug.h"
#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
using namespace ThePEG;
LesHouchesReader::LesHouchesReader(bool active)
: theNEvents(0), position(0), reopened(0), theMaxScan(-1), scanning(false),
isActive(active), theCacheFileName(""), doCutEarly(true),
preweight(1.0), reweightPDF(false), doInitPDFs(false),
theMaxMultCKKW(0), theMinMultCKKW(0), lastweight(1.0), maxFactor(1.0), optionalnpLO(0), optionalnpNLO(0),
weightScale(1.0*picobarn), skipping(false), theMomentumTreatment(0),
useWeightWarnings(true),theReOpenAllowed(true), theIncludeSpin(true) {}
LesHouchesReader::LesHouchesReader(const LesHouchesReader & x)
: HandlerBase(x), LastXCombInfo<>(x), heprup(x.heprup), hepeup(x.hepeup),
inData(x.inData), inPDF(x.inPDF), outPDF(x.outPDF),
thePartonExtractor(x.thePartonExtractor), thePartonBins(x.thePartonBins),
theXCombs(x.theXCombs), theCuts(x.theCuts),
theNEvents(x.theNEvents), position(x.position), reopened(x.reopened),
theMaxScan(x.theMaxScan), scanning(false),
isActive(x.isActive),
theCacheFileName(x.theCacheFileName), doCutEarly(x.doCutEarly),
stats(x.stats), statmap(x.statmap),
thePartonBinInstances(x.thePartonBinInstances),
reweights(x.reweights), preweights(x.preweights),
preweight(x.preweight), reweightPDF(x.reweightPDF),
doInitPDFs(x.doInitPDFs),
theMaxMultCKKW(x.theMaxMultCKKW), theMinMultCKKW(x.theMinMultCKKW),
lastweight(x.lastweight), maxFactor(x.maxFactor),
weightScale(x.weightScale), xSecWeights(x.xSecWeights),
maxWeights(x.maxWeights), skipping(x.skipping),
theMomentumTreatment(x.theMomentumTreatment),
useWeightWarnings(x.useWeightWarnings),
theReOpenAllowed(x.theReOpenAllowed),
theIncludeSpin(x.theIncludeSpin) {}
LesHouchesReader::~LesHouchesReader() {}
void LesHouchesReader::doinitrun() {
HandlerBase::doinitrun();
stats.reset();
for ( StatMap::iterator i = statmap.begin(); i != statmap.end(); ++i )
i->second.reset();
open();
if ( cacheFileName().length() ) openReadCacheFile();
position = 0;
reopened = 0;
}
bool LesHouchesReader::preInitialize() const {
if ( HandlerBase::preInitialize() ) return true;
if ( doInitPDFs && ! ( inPDF.first && inPDF.second ) ) return true;
return false;
}
void LesHouchesReader::doinit() {
HandlerBase::doinit();
open();
close();
if ( !heprup.IDBMUP.first || !heprup.IDBMUP.second )
Throw<LesHouchesInitError>()
<< "No information about incoming particles were found in "
<< "LesHouchesReader '" << name() << "'." << Exception::warning;
inData = make_pair(getParticleData(heprup.IDBMUP.first),
getParticleData(heprup.IDBMUP.second));
if ( heprup.EBMUP.first <= 0.0 || heprup.EBMUP.second <= 0.0 )
Throw<LesHouchesInitError>()
<< "No information about the energy of incoming particles were found in "
<< "LesHouchesReader '" << name() << "'." << Exception::warning;
if ( doInitPDFs && ! ( inPDF.first && inPDF.second ) ) {
initPDFs();
if ( ! ( inPDF.first && inPDF.second ) ) Throw<InitException>()
<< "LesHouchesReader '" << name()
<< "' could not create PDFBase objects in pre-initialization."
<< Exception::warning;
}
else if ( !inPDF.first || !inPDF.second ) Throw<LesHouchesInitError>()
<< "No information about the PDFs of incoming particles were found in "
<< "LesHouchesReader '" << name() << "'." << Exception::warning;
}
void LesHouchesReader::initPDFs() {
if ( inPDF.first && inPDF.second ) return;
string remhname;
if ( heprup.PDFSUP.first && !inPDF.first) {
inPDF.first = dynamic_ptr_cast<PDFPtr>
(generator()->preinitCreate("ThePEG::LHAPDF", fullName() + "/PDFA",
"ThePEGLHAPDF.so"));
if ( !inPDF.first ) {
Throw<InitException>()
<< "LesHouchesReader '" << name() << "' could not use information "
<< "about the PDFs used because the LHAPDF library was not properly "
"defined." << Exception::warning;
return;
}
remhname = fullName() + "/DummyRemH";
generator()->preinitCreate("ThePEG::NoRemnants", remhname);
generator()->preinitInterface(inPDF.first, "RemnantHandler",
"set", remhname);
if ( heprup.PDFGUP.first > 0 && heprup.PDFGUP.first < 10 ) {
ostringstream os;
os << heprup.PDFGUP.first << " " << heprup.PDFSUP.first;
generator()->preinitInterface(inPDF.first, "PDFLIBNumbers",
"set", os.str());
} else {
ostringstream os;
os << heprup.PDFGUP.first*1000 + heprup.PDFSUP.first;
generator()->preinitInterface(inPDF.first, "PDFNumber",
"set", os.str());
}
generator()->preinitInterface(inPDF.first, "RangeException",
"newdef", "Freeze");
}
if ( heprup.PDFSUP.second && !inPDF.second) {
inPDF.second = dynamic_ptr_cast<PDFPtr>
(generator()->preinitCreate("ThePEG::LHAPDF", fullName() + "/PDFB",
"ThePEGLHAPDF.so"));
if ( !inPDF.second ) {
Throw<InitException>()
<< "LesHouchesReader '" << name() << "' could not use information "
<< "about the PDFs used because the LHAPDF library was not properly "
"defined." << Exception::warning;
return;
}
if ( remhname == "" ) {
remhname = fullName() + "/DummyRemH";
generator()->preinitCreate("ThePEG::NoRemnants", remhname);
}
generator()->preinitInterface(inPDF.second, "RemnantHandler",
"set", remhname);
if ( heprup.PDFGUP.second > 0 && heprup.PDFGUP.second < 10 ) {
ostringstream os;
os << heprup.PDFGUP.second << " " << heprup.PDFSUP.second;
generator()->preinitInterface(inPDF.second, "PDFLIBNumbers",
"set", os.str());
} else {
ostringstream os;
os << heprup.PDFGUP.second*1000 + heprup.PDFSUP.second;
generator()->preinitInterface(inPDF.second, "PDFNumber",
"set", os.str());
}
generator()->preinitInterface(inPDF.second, "RangeException",
"newdef", "Freeze");
}
if ( ! ( inPDF.first && inPDF.second ) ) Throw<InitException>()
<< "LesHouchesReader '" << name()
<< "' could not find information about the PDFs used."
<< Exception::warning;
}
void LesHouchesReader::initialize(LesHouchesEventHandler & eh) {
Energy2 Smax = ZERO;
double Y = 0.0;
if ( !theCuts ) {
theCuts = eh.cuts();
if ( !theCuts ) Throw<LesHouchesInitError>()
<< "No Cuts object was assigned to the LesHouchesReader '"
<< name() << "' nor was one\nassigned to the controlling "
<< "LesHouchesEventHandler '" << eh.name() << "'.\nAt least one of them "
<< "needs to have a Cuts object." << Exception::runerror;
Smax = cuts().SMax();
Y = cuts().Y();
}
theCKKW = eh.CKKWHandler();
if ( !partonExtractor() ) {
thePartonExtractor = eh.partonExtractor();
if ( !partonExtractor() ) Throw<LesHouchesInitError>()
<< "No PartonExtractor object was assigned to the LesHouchesReader '"
<< name() << "' nor was one\nassigned to the controlling "
<< "LesHouchesEventHandler '" << eh.name() << "'.\nAt least one of them "
<< "needs to have a PartonExtractor object." << Exception::runerror;
}
open();
Energy emax = 2.0*sqrt(heprup.EBMUP.first*heprup.EBMUP.second)*GeV;
theCuts->initialize(sqr(emax),
0.5*log(heprup.EBMUP.first/heprup.EBMUP.second));
if ( Smax > ZERO && ( Smax != cuts().SMax() || Y != cuts().Y() ) )
Throw<LesHouchesInitError>()
<< "The LesHouchesReader '" << name() << "' uses the same Cuts object "
<< "as another LesHouchesReader which has not got the same energies of "
<< "the colliding particles. For the generation to work properly "
<< "different LesHouchesReader object with different colliding particles "
<< "must be assigned different (although possibly identical) Cuts "
<< "objects." << Exception::warning;
thePartonBins = partonExtractor()->getPartons(emax, inData, cuts());
for ( int i = 0, N = partonBins().size(); i < N; ++i ) {
theXCombs[partonBins()[i]] =
new_ptr(XComb(emax, inData, &eh, partonExtractor(), CKKWHandler(),
partonBins()[i], theCuts));
partonExtractor()->nDims(partonBins()[i]);
}
outPDF = make_pair(partonExtractor()->getPDF(inData.first),
partonExtractor()->getPDF(inData.second));
// SP: re-interpret 3/4 -> 1/2; see discussion w/ Leif and Keith at LH 2013
if ( abs(heprup.IDWTUP) == 3 )
heprup.IDWTUP = heprup.IDWTUP < 0 ? -1 : 1;
if ( abs(heprup.IDWTUP) == 4 )
heprup.IDWTUP = heprup.IDWTUP < 0 ? -2 : 2;
close();
if ( !heprup.IDWTUP && useWeightWarnings )
Throw<LesHouchesInitError>()
<< "No information about the weighting scheme was found. The events "
<< "produced by LesHouchesReader " << name()
<< " may not be sampled correctly." << Exception::warning;
if ( abs(heprup.IDWTUP) > 1 && !eh.weighted() && useWeightWarnings )
Throw<LesHouchesInitError>()
<< "LesHouchesReader " << name() << " has the IDWTUP flag set to "
<< heprup.IDWTUP << " which is not supported by this reader, the "
<< "produced events may not be sampled correctly. It is up to "
<< "sub-classes of LesHouchesReader to correctly convert to match IDWTUP "
<< "+/- 1. Will try to make intelligent guesses to get "
<< "correct statistics.\nIn most cases this should be sufficient. "
<< "Unset <interface>WeightWarnings</interface> to avoid this message,"
<< "or set <interface>Weighted</interface> to on."
<< Exception::warning;
if ( heprup.IDWTUP != eh.weightOption() && abs(heprup.IDWTUP) < 3 &&
useWeightWarnings )
Throw<LesHouchesInitError>()
<< "LesHouchesReader " << name() << " has the IDWTUP flag set to "
<< heprup.IDWTUP
<< ", which does not correspond\nto the weight option "
<< eh.weightOption() << " set in "
<< "the LesHouchesEventHandler " << eh.name() << ".\n\n"
<< "Use the following handler setting instead:\n"
<< " set " << eh.name() << ":WeightOption " << heprup.IDWTUP
<< "\nWill try to make intelligent guesses to get "
<< "correct statistics. In most cases this should be sufficient. "
<< "Unset <interface>WeightWarnings</interface> to avoid this message"
<< Exception::warning;
scan();
initStat();
}
long LesHouchesReader::scan() {
open();
// Shall we write the events to a cache file for fast reading? If so
// we write to a temporary file if the caches events should be
// randomized.
if ( cacheFileName().length() ) openWriteCacheFile();
// Keep track of the number of events scanned.
long neve = 0;
long cuteve = 0;
bool negw = false;
// If the open() has not already gotten information about subprocesses
// and cross sections we have to scan through the events.
if ( !heprup.NPRUP || cacheFile() || abs(heprup.IDWTUP) != 1 ) { // why scan if IDWTUP != 1?
HoldFlag<> isScanning(scanning);
double oldsum = 0.0;
vector<int> lprup;
vector<double> newmax;
vector<long> oldeve;
vector<long> neweve;
for ( int i = 0; ( maxScan() < 0 || i < maxScan() ) && readEvent(); ++i ) {
if ( !checkPartonBin() ) Throw<LesHouchesInitError>()
<< "Found event in LesHouchesReader '" << name()
<< "' which cannot be handeled by the assigned PartonExtractor '"
<< partonExtractor()->name() << "'." << Exception::runerror;
vector<int>::iterator idit =
find(lprup.begin(), lprup.end(), hepeup.IDPRUP);
int id = lprup.size();
if ( idit == lprup.end() ) {
lprup.push_back(hepeup.IDPRUP);
newmax.push_back(0.0);
neweve.push_back(0);
oldeve.push_back(0);
} else {
id = idit - lprup.begin();
}
++neve;
++oldeve[id];
oldsum += hepeup.XWGTUP;
if ( cacheFile() ) {
if ( eventWeight() == 0.0 ) {
++cuteve;
continue;
}
cacheEvent();
}
++neweve[id];
newmax[id] = max(newmax[id], abs(eventWeight()));
if ( eventWeight() < 0.0 ) negw = true;
}
// std::cout << "eventWeight() = " << eventWeight() << endl;
xSecWeights.resize(oldeve.size(), 1.0);
for ( int i = 0, N = oldeve.size(); i < N; ++i )
if ( oldeve[i] ) xSecWeights[i] = double(neweve[i])/double(oldeve[i]);
if ( maxScan() < 0 || neve > NEvents() ) NEvents(neve - cuteve);
if ( lprup.size() == heprup.LPRUP.size() ) {
for ( int id = 0, N = lprup.size(); id < N; ++id ) {
vector<int>::iterator idit =
find(heprup.LPRUP.begin(), heprup.LPRUP.end(), hepeup.IDPRUP);
if ( idit == heprup.LPRUP.end() ) {
Throw<LesHouchesInitError>()
<< "When scanning events, the LesHouschesReader '" << name()
<< "' found undeclared processes." << Exception::warning;
heprup.NPRUP = 0;
break;
}
int idh = idit - heprup.LPRUP.begin();
heprup.XMAXUP[idh] = newmax[id];
}
}
if ( heprup.NPRUP == 0 ) {
// No heprup block was supplied or something went wrong.
heprup.NPRUP = lprup.size();
heprup.LPRUP.resize(lprup.size());
heprup.XMAXUP.resize(lprup.size());
for ( int id = 0, N = lprup.size(); id < N; ++id ) {
heprup.LPRUP[id] = lprup[id];
heprup.XMAXUP[id] = newmax[id];
}
} else if ( abs(heprup.IDWTUP) != 1 ) {
// Try to fix things if abs(heprup.IDWTUP) != 1.
double sumxsec = 0.0;
for ( int id = 0; id < heprup.NPRUP; ++id ) sumxsec += heprup.XSECUP[id];
weightScale = picobarn*neve*sumxsec/oldsum;
}
}
if ( cacheFile() ) closeCacheFile();
if ( negw ) heprup.IDWTUP = min(-abs(heprup.IDWTUP), -1);
return neve;
}
void LesHouchesReader::setWeightScale(long) {}
void LesHouchesReader::initStat() {
stats.reset();
statmap.clear();
if ( heprup.NPRUP <= 0 ) return;
double sumx = 0.0;
xSecWeights.resize(heprup.NPRUP, 1.0);
maxWeights.clear();
for ( int ip = 0; ip < heprup.NPRUP; ++ip ) {
sumx = max(heprup.XMAXUP[ip]*xSecWeights[ip], sumx);
statmap[heprup.LPRUP[ip]] =
XSecStat(heprup.XMAXUP[ip]*weightScale*xSecWeights[ip]);
maxWeights[heprup.LPRUP[ip]] = heprup.XMAXUP[ip];
}
stats.maxXSec(sumx*weightScale);
maxFactor = 1.0;
}
void LesHouchesReader::increaseMaxXSec(CrossSection maxxsec) {
for ( int i = 0; i < heprup.NPRUP; ++i )
statmap[heprup.LPRUP[i]].maxXSec(statmap[heprup.LPRUP[i]].maxXSec()*
maxxsec/stats.maxXSec());
maxFactor *= maxxsec/stats.maxXSec();
stats.maxXSec(maxxsec);
}
tXCombPtr LesHouchesReader::getXComb() {
if ( lastXCombPtr() ) return lastXCombPtr();
fillEvent();
connectMothers();
tcPBPair sel = createPartonBinInstances();
tXCombPtr lastXC = xCombs()[sel];
// clean up the old XComb object before switching to a new one
if ( theLastXComb && theLastXComb != lastXC )
theLastXComb->clean();
theLastXComb = lastXC;
lastXCombPtr()->subProcess(SubProPtr());
lastXCombPtr()->setPartonBinInstances(partonBinInstances(),
sqr(hepeup.SCALUP)*GeV2);
lastXCombPtr()->lastAlphaS(hepeup.AQCDUP);
lastXCombPtr()->lastAlphaEM(hepeup.AQEDUP);
return lastXCombPtr();
}
tSubProPtr LesHouchesReader::getSubProcess() {
getXComb();
if ( subProcess() ) return subProcess();
lastXCombPtr()->subProcess(new_ptr(SubProcess(lastPartons(), tCollPtr(), this)));
subProcess()->setOutgoing(outgoing().begin(), outgoing().end());
subProcess()->setIntermediates(intermediates().begin(),
intermediates().end());
return subProcess();
}
void LesHouchesReader::fillEvent() {
if ( !particleIndex.empty() ) return;
particleIndex.clear();
colourIndex.clear();
colourIndex(0, tColinePtr());
createParticles();
createBeams();
}
void LesHouchesReader::reopen() {
// If we didn't know how many events there were, we know now.
if ( NEvents() <= 0 ) NEvents(position);
++reopened;
// How large fraction of the events have we actually used? And how
// large will we have used if we go through the file again?
double frac = double(stats.attempts())/double(NEvents());
if ( frac*double(reopened + 1)/double(reopened) > 1.0 &&
NEvents() - stats.attempts() <
generator()->N() - generator()->currentEventNumber() ) {
if(theReOpenAllowed)
generator()->logWarning(LesHouchesReopenWarning()
<< "Reopening LesHouchesReader '" << name()
<< "' after accessing " << stats.attempts()
<< " events out of "
<< NEvents() << Exception::warning);
else
throw LesHouchesReopenWarning()
<< "More events requested than available in LesHouchesReader "
<< name() << Exception::runerror;
}
if ( cacheFile() ) {
closeCacheFile();
openReadCacheFile();
if ( !uncacheEvent() ) Throw<LesHouchesReopenError>()
<< "Could not reopen LesHouchesReader '" << name()
<< "'." << Exception::runerror;
} else {
close();
open();
if ( !readEvent() ) Throw<LesHouchesReopenError>()
<< "Could not reopen LesHouchesReader '" << name()
<< "'." << Exception::runerror;
}
}
void LesHouchesReader::reset() {
particleIndex.clear();
colourIndex.clear();
if ( theLastXComb ) theLastXComb->clean();
theLastXComb = tXCombPtr();
}
bool LesHouchesReader::readEvent() {
reset();
if ( !doReadEvent() ) return false;
// If we are just skipping event we do not need to reweight or do
// anything fancy.
if ( skipping ) { return true; }
if ( cacheFile() && !scanning ) { return true; }
// Reweight according to the re- and pre-weights objects in the
// LesHouchesReader base class.
lastweight = reweight();
if ( !reweightPDF && !cutEarly() ) return true;
// We should try to reweight the PDFs or make early cuts here.
fillEvent();
double x1 = incoming().first->momentum().plus()/
beams().first->momentum().plus();
if ( reweightPDF &&
inPDF.first && outPDF.first && inPDF.first != outPDF.first ) {
if ( hepeup.XPDWUP.first <= 0.0 )
hepeup.XPDWUP.first =
inPDF.first->xfx(inData.first, incoming().first->dataPtr(),
sqr(hepeup.SCALUP*GeV), x1);
double xf = outPDF.first->xfx(inData.first, incoming().first->dataPtr(),
sqr(hepeup.SCALUP*GeV), x1);
lastweight *= xf/hepeup.XPDWUP.first;
hepeup.XPDWUP.first = xf;
}
double x2 = incoming().second->momentum().minus()/
beams().second->momentum().minus();
if ( reweightPDF &&
inPDF.second && outPDF.second && inPDF.second != outPDF.second ) {
if ( hepeup.XPDWUP.second <= 0.0 )
hepeup.XPDWUP.second =
inPDF.second->xfx(inData.second, incoming().second->dataPtr(),
sqr(hepeup.SCALUP*GeV), x2);
double xf =
outPDF.second->xfx(inData.second, incoming().second->dataPtr(),
sqr(hepeup.SCALUP*GeV), x2);
lastweight *= xf/hepeup.XPDWUP.second;
hepeup.XPDWUP.second = xf;
}
if ( cutEarly() ) {
if ( !cuts().initSubProcess((incoming().first->momentum() +
incoming().second->momentum()).m2(),
0.5*log(x1/x2)) ) lastweight = 0.0;
tSubProPtr sub = getSubProcess();
TmpTransform<tSubProPtr> tmp(sub, Utilities::getBoostToCM(sub->incoming()));
if ( !cuts().passCuts(*sub) ) lastweight = 0.0;
}
return true;
}
double LesHouchesReader::getEvent() {
if ( cacheFile() ) {
if ( !uncacheEvent() ) reopen();
} else {
if ( !readEvent() ) reopen();
}
++position;
double max = maxWeights[hepeup.IDPRUP]*maxFactor;
return max != 0.0? eventWeight()/max: 0.0;
}
void LesHouchesReader::skip(long n) {
HoldFlag<> skipflag(skipping);
while ( n-- ) getEvent();
}
double LesHouchesReader::reweight() {
preweight = 1.0;
if ( reweights.empty() && preweights.empty() &&
!( CKKWHandler() && maxMultCKKW() > 0 && maxMultCKKW() > minMultCKKW() ) )
return 1.0;
fillEvent();
getSubProcess();
for ( int i = 0, N = preweights.size(); i < N; ++i ) {
preweights[i]->setXComb(lastXCombPtr());
preweight *= preweights[i]->weight();
}
double weight = preweight;
for ( int i = 0, N = reweights.size(); i < N; ++i ) {
reweights[i]->setXComb(lastXCombPtr());
weight *= reweights[i]->weight();
}
// If we are caching events we do not want to do CKKW reweighting.
if ( cacheFile() ) return weight;
if ( CKKWHandler() && maxMultCKKW() > 0 && maxMultCKKW() > minMultCKKW() ) {
CKKWHandler()->setXComb(lastXCombPtr());
weight *= CKKWHandler()->reweightCKKW(minMultCKKW(), maxMultCKKW());
}
return weight;
}
bool LesHouchesReader::checkPartonBin() {
// First find the positions of the incoming partons.
pair< vector<int>, vector<int> > inc;
for ( int i = 0; i < hepeup.NUP; ++i ) {
if ( hepeup.ISTUP[i] == -9 ) {
if ( inc.first.empty() ) inc.first.push_back(i);
else if ( inc.second.empty() ) inc.second.push_back(i);
}
else if ( hepeup.ISTUP[i] == -1 ) {
if ( inc.first.size() &&
hepeup.MOTHUP[i].first == inc.first.back() + 1 )
inc.first.push_back(i);
else if ( inc.second.size() &&
hepeup.MOTHUP[i].first == inc.second.back() + 1 )
inc.second.push_back(i);
else if ( inc.first.empty() ) {
inc.first.push_back(-1);
inc.first.push_back(i);
}
else if ( inc.second.empty() ) {
inc.second.push_back(-1);
inc.second.push_back(i);
}
else if ( inc.first.size() <= inc.second.size() )
inc.first.push_back(i);
else
inc.second.push_back(i);
}
}
// Now store the corresponding id numbers
pair< vector<long>, vector<long> > ids;
ids.first.push_back(inc.first[0] < 0? heprup.IDBMUP.first:
hepeup.IDUP[inc.first[0]]);
for ( int i = 1, N = inc.first.size(); i < N; ++i )
ids.first.push_back(hepeup.IDUP[inc.first[i]]);
ids.second.push_back(inc.second[0] < 0? heprup.IDBMUP.second:
hepeup.IDUP[inc.second[0]]);
for ( int i = 1, N = inc.second.size(); i < N; ++i )
ids.second.push_back(hepeup.IDUP[inc.second[i]]);
// Find the correct pair of parton bins.
PBPair pbp;
for ( int i = 0, N = partonBins().size(); i < N; ++i ) {
tcPBPtr curr = partonBins()[i].first;
int icurr = inc.first.size() - 1;
while ( curr && icurr >= 0 ) {
if ( curr->parton()->id () != ids.first[icurr] ) break;
curr = curr->incoming();
--icurr;
}
if(!(!partonBins()[i].first->incoming() &&
!partonBins()[i].first->particle() &&
partonBins()[i].first->parton()->id () == ids.first[0] &&
( inc.first.size()==1 ||
(inc.first.size()==2 && ids.first[0]==ids.first[1]))) &&
( curr || icurr >= 0 ) ) continue;
curr = partonBins()[i].second;
icurr = inc.second.size() - 1;
while ( curr && icurr >= 0 ) {
if ( curr->parton()->id () != ids.second[icurr] ) break;
curr = curr->incoming();
--icurr;
}
if(!(!partonBins()[i].second->incoming() &&
!partonBins()[i].second->particle() &&
partonBins()[i].second->parton()->id () == ids.second[0] &&
( inc.second.size()==1 ||
(inc.second.size()==2 && ids.second[0]==ids.second[1]))) &&
( curr || icurr >= 0 ) ) continue;
pbp = partonBins()[i];
}
// If we are only checking we return here.
return ( pbp.first && pbp.second );
}
namespace {
bool recursionNotNull(tcPBPtr bin, tcPPtr p) {
while ( bin && p ) {
if ( p->dataPtr() != bin->parton() ) break;
bin = bin->incoming();
p = p->parents().size()? p->parents()[0]: tPPtr();
}
return bin || p;
}
}
tcPBPair LesHouchesReader::createPartonBinInstances() {
tcPBPair sel;
for ( int i = 0, N = partonBins().size(); i < N; ++i ) {
tcPBPtr bin = partonBins()[i].first;
tcPPtr p = incoming().first;
if ( recursionNotNull(bin,p) ) continue;
bin = partonBins()[i].second;
p = incoming().second;
if ( recursionNotNull(bin,p) ) continue;
sel = partonBins()[i];
break;
}
if ( !sel.first || !sel.second ) Throw<LesHouchesInconsistencyError>()
<< "Could not find appropriate PartonBin objects for event produced by "
<< "LesHouchesReader '" << name() << "'." << Exception::runerror;
Direction<0> dir(true);
thePartonBinInstances.first =
new_ptr(PartonBinInstance(incoming().first, sel.first,
-sqr(hepeup.SCALUP*GeV)));
if ( thePartonBinInstances.first->xi() > 1.00001 ) {
Throw<LesHouchesInconsistencyError>()
<< "Found an event with momentum fraction larger than unity (x1="
<< thePartonBinInstances.first->xi()
<< "). The event will be skipped." << Exception::warning;
throw Veto();
}
dir.reverse();
thePartonBinInstances.second =
new_ptr(PartonBinInstance(incoming().second, sel.second,
-sqr(hepeup.SCALUP*GeV)));
if ( thePartonBinInstances.second->xi() > 1.00001 ) {
Throw<LesHouchesInconsistencyError>()
<< "Found an event with momentum fraction larger than unity (x2="
<< thePartonBinInstances.second->xi()
<< "). The event will be skipped." << Exception::warning;
throw Veto();
}
return sel;
}
void LesHouchesReader::createParticles() {
theBeams = PPair();
theIncoming = PPair();
theOutgoing = PVector();
theIntermediates = PVector();
for ( int i = 0, N = hepeup.IDUP.size(); i < N; ++i ) {
if ( !hepeup.IDUP[i] ) continue;
Lorentz5Momentum mom(hepeup.PUP[i][0]*GeV, hepeup.PUP[i][1]*GeV,
hepeup.PUP[i][2]*GeV, hepeup.PUP[i][3]*GeV,
hepeup.PUP[i][4]*GeV);
// cout << hepeup.PUP[i][0] << " " << hepeup.PUP[i][1] << " " << hepeup.PUP[i][2] << " " << hepeup.PUP[i][3] << " " << hepeup.PUP[i][4] << endl;
if(theMomentumTreatment == 1) mom.rescaleEnergy();
else if(theMomentumTreatment == 2) mom.rescaleMass();
// cout << hepeup.PUP[i][0] << " " << hepeup.PUP[i][1] << " " << hepeup.PUP[i][2] << " " << hepeup.PUP[i][3] << " " << hepeup.PUP[i][4] << endl;
PDPtr pd = getParticleData(hepeup.IDUP[i]);
if (!pd) {
Throw<LesHouchesInitError>()
<< "LesHouchesReader '" << name() << "' found unknown particle ID "
<< hepeup.IDUP[i]
<< " in Les Houches common block structure.\n"
<< "You need to define the new particle in an input file.\n"
<< Exception::runerror;
}
if ( ! pd->coloured()
&& ( hepeup.ICOLUP[i].first != 0 || hepeup.ICOLUP[i].second != 0 ) ) {
Throw<LesHouchesInconsistencyError>()
<< "LesHouchesReader " << name() << ": " << pd->PDGName()
<< " is not a coloured particle.\nIt should not have "
<< "(anti-)colour lines " << hepeup.ICOLUP[i].first
<< ' ' << hepeup.ICOLUP[i].second
<< " set; the event file needs to be fixed."
<< Exception::runerror;
}
PPtr p = pd->produceParticle(mom);
if(hepeup.ICOLUP[i].first>=0 && hepeup.ICOLUP[i].second >=0) {
tColinePtr c = colourIndex(hepeup.ICOLUP[i].first);
if ( c ) {
c->addColoured(p);
}
c = colourIndex(hepeup.ICOLUP[i].second);
if ( c ) c->addAntiColoured(p);
}
else {
tColinePtr c1 = colourIndex(abs(hepeup.ICOLUP[i].first ));
tColinePtr c2 = colourIndex(abs(hepeup.ICOLUP[i].second));
if(pd->hasColour()) {
c1->addColouredIndexed(p,1);
c2->addColouredIndexed(p,2);
}
else {
c1->addAntiColouredIndexed(p,1);
c2->addAntiColouredIndexed(p,2);
}
}
particleIndex(i + 1, p);
switch ( hepeup.ISTUP[i] ) {
case -9:
if ( !theBeams.first ) theBeams.first = p;
else if ( !theBeams.second ) theBeams.second = p;
else Throw<LesHouchesInconsistencyError>()
<< "To many incoming beam particles in the LesHouchesReader '"
<< name() << "'." << Exception::runerror;
break;
case -1:
if ( !theIncoming.first ) theIncoming.first = p;
else if ( !theIncoming.second ) theIncoming.second = p;
else if ( particleIndex(theIncoming.first) == hepeup.MOTHUP[i].first )
theIncoming.first = p;
else if ( particleIndex(theIncoming.second) == hepeup.MOTHUP[i].first )
theIncoming.second = p;
else Throw<LesHouchesInconsistencyError>()
<< "To many incoming particles to hard subprocess in the "
<< "LesHouchesReader '" << name() << "'." << Exception::runerror;
p->scale(sqr(hepeup.SCALUP*GeV));
break;
case 1:
theOutgoing.push_back(p);
p->scale(sqr(hepeup.SCALUP*GeV));
break;
case -2:
case 2:
case 3:
theIntermediates.push_back(p);
break;
default:
Throw<LesHouchesInconsistencyError>()
<< "Unknown status code (" << hepeup.ISTUP[i]
<< ") in the LesHouchesReader '" << name() << "'."
<< Exception::runerror;
}
// value 9 is defined as "Unknown or unpolarized particles"
double spinup = hepeup.SPINUP[i];
if ( abs(spinup - 9) < 1.0e-3 )
spinup = 0.;
if ( spinup < -1. || spinup > 1. ) {
Throw<LesHouchesInconsistencyError>()
<< "Polarization must be between -1 and 1, not "
<< spinup << " as found in the "
<< "LesHouches event file.\nThe event file needs to be fixed."
<< Exception::runerror;
}
if( theIncludeSpin
&& abs(pd->id()) == ParticleID::tauminus
&& spinup !=0) {
if(pd->iSpin() == PDT::Spin1Half ) {
vector<Helicity::SpinorWaveFunction> wave;
Helicity::SpinorWaveFunction(wave,p,Helicity::outgoing,true);
RhoDMatrix rho(pd->iSpin(),true);
rho(0,0) = 0.5*(1.-spinup);
rho(1,1) = 0.5*(1.+spinup);
p->spinInfo()->rhoMatrix() = rho;
p->spinInfo()-> DMatrix() = rho;
}
}
}
// check the colour flows, and if necessary create any sources/sinks
// hard process
// get the particles in the hard process
PVector external;
for ( int i = 0, N = hepeup.IDUP.size(); i < N; ++i ) {
unsigned int moth;
switch ( hepeup.ISTUP[i] ) {
case -1:
external.push_back(particleIndex.find(i+1));
break;
case 1: case 2: case 3:
moth = hepeup.MOTHUP[i].first;
if(moth!=0 && (hepeup.ISTUP[moth]==-1||hepeup.ISTUP[moth]==-2||
hepeup.ISTUP[moth]==-9))
external.push_back(particleIndex.find(i+1));
moth = hepeup.MOTHUP[i].second;
if(moth!=0 && (hepeup.ISTUP[moth]==-1||hepeup.ISTUP[moth]==-2||
hepeup.ISTUP[moth]==-9))
external.push_back(particleIndex.find(i+1));
break;
case -2: case -9: default:
break;
}
}
// check the incoming/outgoing lines match
vector<tColinePtr> unMatchedColour,unMatchedAntiColour;
for(unsigned int ix=0;ix<external.size();++ix) {
vector<tcColinePtr>
col = external[ix]->colourInfo()-> colourLines();
vector<tcColinePtr>
anti = external[ix]->colourInfo()->antiColourLines();
if(hepeup.ISTUP[particleIndex(external[ix])-1]<0)
swap(col,anti);
if(!col.empty()) {
for(unsigned int ic1=0;ic1<col.size();++ic1) {
bool matched=false;
for(unsigned int iy=0;iy<external.size();++iy) {
vector<tcColinePtr> col2;
if(hepeup.ISTUP[particleIndex(external[iy])-1]<0) {
if(external[iy]->colourInfo()->colourLines().empty()) continue;
col2 = external[iy]->colourInfo()->colourLines();
}
else if(hepeup.ISTUP[particleIndex(external[iy])-1]>0) {
if(external[iy]->colourInfo()->antiColourLines().empty()) continue;
col2 = external[iy]->colourInfo()->antiColourLines();
}
for(unsigned int ic2=0;ic2<col2.size();++ic2) {
if(col[ic1]==col2[ic2]) {
matched=true;
break;
}
}
if(matched) break;
}
if(!matched) unMatchedColour.push_back(const_ptr_cast<tColinePtr>(col[ic1]));
}
}
if(!anti.empty()) {
for(unsigned int ic1=0;ic1<col.size();++ic1) {
bool matched=false;
for(unsigned int iy=0;iy<external.size();++iy) {
vector<tcColinePtr> anti2;
if(hepeup.ISTUP[particleIndex(external[iy])-1]<0) {
if(external[iy]->colourInfo()->colourLines().empty()) continue;
anti2 = external[iy]->colourInfo()->antiColourLines();
}
else if(hepeup.ISTUP[particleIndex(external[iy])-1]>0) {
if(external[iy]->colourInfo()->antiColourLines().empty()) continue;
anti2 = external[iy]->colourInfo()->colourLines();
}
for(unsigned int ic2=0;ic2<anti2.size();++ic2) {
if(col[ic1]==anti2[ic2]) {
matched=true;
break;
}
}
if(matched) break;
}
if(!matched) unMatchedAntiColour.push_back(const_ptr_cast<tColinePtr>(anti[ic1]));
}
}
}
// might have source/sink
if( unMatchedColour.size() + unMatchedAntiColour.size() != 0) {
if(unMatchedColour.size() == 3 ) {
unMatchedColour[0]->setSourceNeighbours(unMatchedColour[1],
unMatchedColour[2]);
}
else if(unMatchedColour.size() != 0 && ThePEG_DEBUG_LEVEL) {
Throw<LesHouchesInconsistencyError>()
<< "LesHouchesReader '" << name() << "' found inconsistent colour "
<< "flow in Les Houches common block structure for hard process.\n"
<< hepeup << Exception::runerror;
}
if(unMatchedAntiColour.size() == 3 ) {
unMatchedAntiColour[0]->setSinkNeighbours(unMatchedAntiColour[1],
unMatchedAntiColour[2]);
}
else if(unMatchedAntiColour.size() != 0 && ThePEG_DEBUG_LEVEL) {
Throw<LesHouchesInconsistencyError>()
<< "LesHouchesReader '" << name() << "' found inconsistent colour "
<< "flow in Les Houches common block structure for hard process.\n"
<< hepeup << Exception::runerror;
}
}
// any subsequent decays
for ( int i = 0, N = hepeup.IDUP.size(); i < N; ++i ) {
if(hepeup.ISTUP[i] !=2 && hepeup.ISTUP[i] !=3) continue;
PVector external;
external.push_back(particleIndex.find(i+1));
for ( int j = 0; j < N; ++j ) {
if(hepeup.MOTHUP[j].first==i+1|| hepeup.MOTHUP[j].second==i+1)
external.push_back(particleIndex.find(j+1));
}
// check the incoming/outgoing lines match
vector<tColinePtr> unMatchedColour,unMatchedAntiColour;
for(unsigned int ix=0;ix<external.size();++ix) {
vector<tcColinePtr>
col = external[ix]->colourInfo()-> colourLines();
vector<tcColinePtr>
anti = external[ix]->colourInfo()->antiColourLines();
if(ix==0) swap(col,anti);
if(!col.empty()) {
for(unsigned int ic1=0;ic1<col.size();++ic1) {
bool matched=false;
for(unsigned int iy=0;iy<external.size();++iy) {
if(iy==ix) continue;
vector<tcColinePtr> col2;
if(iy==0) {
if(external[iy]->colourInfo()->colourLines().empty()) continue;
col2 = external[iy]->colourInfo()->colourLines();
}
else {
if(external[iy]->colourInfo()->antiColourLines().empty()) continue;
col2 = external[iy]->colourInfo()->antiColourLines();
}
for(unsigned int ic2=0;ic2<col2.size();++ic2) {
if(col[ic1]==col2[ic2]) {
matched=true;
break;
}
}
if(matched) break;
}
if(!matched) unMatchedColour.push_back(const_ptr_cast<tColinePtr>(col[ic1]));
}
}
if(!anti.empty()) {
for(unsigned int ic1=0;ic1<anti.size();++ic1) {
bool matched=false;
for(unsigned int iy=0;iy<external.size();++iy) {
if(iy==ix) continue;
vector<tcColinePtr> anti2;
if(iy==0) {
if(external[iy]->colourInfo()->antiColourLines().empty()) continue;
anti2 = external[iy]->colourInfo()->antiColourLines();
}
else {
if(external[iy]->colourInfo()->colourLines().empty()) continue;
anti2 = external[iy]->colourInfo()->colourLines();
}
for(unsigned int ic2=0;ic2<anti2.size();++ic2) {
if(anti[ic1]==anti2[ic2]) {
matched=true;
break;
}
}
if(matched) break;
}
if(!matched) unMatchedAntiColour.push_back(const_ptr_cast<tColinePtr>(anti[ic1]));
}
}
}
// might have source/sink
if( unMatchedColour.size() + unMatchedAntiColour.size() != 0) {
if(unMatchedColour.size() == 3 ) {
unMatchedColour[0]->setSourceNeighbours(unMatchedColour[1],
unMatchedColour[2]);
}
else if(unMatchedColour.size() != 0 && ThePEG_DEBUG_LEVEL) {
Throw<LesHouchesInconsistencyError>()
<< "LesHouchesReader '" << name() << "' found inconsistent colour "
<< "flow in Les Houches common block structure for decay of \n"
<< *external[0] << "\n"
<< hepeup << Exception::runerror;
}
if(unMatchedAntiColour.size() == 3 ) {
unMatchedAntiColour[0]->setSinkNeighbours(unMatchedAntiColour[1],
unMatchedAntiColour[2]);
}
else if(unMatchedAntiColour.size() != 0 && ThePEG_DEBUG_LEVEL) {
Throw<LesHouchesInconsistencyError>()
<< "LesHouchesReader '" << name() << "' found inconsistent colour "
<< "flow in Les Houches common block structure for decay of\n"
<< *external[0] << "\n"
<< hepeup << Exception::runerror;
}
}
}
}
void LesHouchesReader::createBeams() {
if ( !theBeams.first && dynamic_ptr_cast<Ptr<NoPDF>::tcp>(inPDF.first) ) {
theBeams.first = theIncoming.first;
}
else if ( !theBeams.first ) {
theBeams.first = getParticleData(heprup.IDBMUP.first)->produceParticle();
double m = theBeams.first->mass()/GeV;
theBeams.first->set5Momentum
(Lorentz5Momentum(ZERO, ZERO,
sqrt(sqr(heprup.EBMUP.first) - sqr(m))*GeV,
heprup.EBMUP.first*GeV, m*GeV));
hepeup.IDUP.push_back(heprup.IDBMUP.first);
hepeup.ISTUP.push_back(-9);
hepeup.MOTHUP.push_back(make_pair(0, 0));
hepeup.ICOLUP.push_back(make_pair(0, 0));
hepeup.VTIMUP.push_back(0.0);
hepeup.SPINUP.push_back(0.0);
particleIndex(hepeup.IDUP.size(), theBeams.first);
hepeup.MOTHUP[particleIndex(theIncoming.first) - 1].first =
hepeup.IDUP.size();
}
if ( !theBeams.second && dynamic_ptr_cast<Ptr<NoPDF>::tcp>(inPDF.second) ) {
theBeams.second = theIncoming.second;
}
else if ( !theBeams.second ) {
theBeams.second = getParticleData(heprup.IDBMUP.second)->produceParticle();
double m = theBeams.second->mass()/GeV;
theBeams.second->set5Momentum
(Lorentz5Momentum(ZERO, ZERO,
-sqrt(sqr(heprup.EBMUP.second) - sqr(m))*GeV,
heprup.EBMUP.second*GeV, m*GeV));
hepeup.IDUP.push_back(heprup.IDBMUP.second);
hepeup.ISTUP.push_back(-9);
hepeup.MOTHUP.push_back(make_pair(0, 0));
hepeup.ICOLUP.push_back(make_pair(0, 0));
hepeup.VTIMUP.push_back(0.0);
hepeup.SPINUP.push_back(0.0);
particleIndex(hepeup.IDUP.size(), theBeams.second);
hepeup.MOTHUP[particleIndex(theIncoming.second) - 1].first =
hepeup.IDUP.size();
}
}
void LesHouchesReader::connectMothers() {
const ObjectIndexer<long,Particle> & pi = particleIndex;
for ( int i = 0, N = hepeup.IDUP.size(); i < N; ++i ) {
if ( pi(hepeup.MOTHUP[i].first) )
pi(hepeup.MOTHUP[i].first)->addChild(pi(i + 1));
if ( pi(hepeup.MOTHUP[i].second)
&& hepeup.MOTHUP[i].second != hepeup.MOTHUP[i].first )
pi(hepeup.MOTHUP[i].second)->addChild(pi(i + 1));
}
}
void LesHouchesReader::openReadCacheFile() {
if ( cacheFile() ) closeCacheFile();
cacheFile().open(cacheFileName(), "r");
position = 0;
}
void LesHouchesReader::openWriteCacheFile() {
if ( cacheFile() ) closeCacheFile();
cacheFile().open(cacheFileName(), "w");
}
void LesHouchesReader::closeCacheFile() {
cacheFile().close();
}
void LesHouchesReader::cacheEvent() const {
static vector<char> buff;
cacheFile().write(&hepeup.NUP, sizeof(hepeup.NUP));
buff.resize(eventSize(hepeup.NUP));
char * pos = &buff[0];
pos = mwrite(pos, hepeup.IDPRUP);
pos = mwrite(pos, hepeup.XWGTUP);
pos = mwrite(pos, hepeup.XPDWUP);
pos = mwrite(pos, hepeup.SCALUP);
pos = mwrite(pos, hepeup.AQEDUP);
pos = mwrite(pos, hepeup.AQCDUP);
pos = mwrite(pos, hepeup.IDUP[0], hepeup.NUP);
pos = mwrite(pos, hepeup.ISTUP[0], hepeup.NUP);
pos = mwrite(pos, hepeup.MOTHUP[0], hepeup.NUP);
pos = mwrite(pos, hepeup.ICOLUP[0], hepeup.NUP);
for ( int i = 0; i < hepeup.NUP; ++i )
pos = mwrite(pos, hepeup.PUP[i][0], 5);
pos = mwrite(pos, hepeup.VTIMUP[0], hepeup.NUP);
pos = mwrite(pos, hepeup.SPINUP[0], hepeup.NUP);
pos = mwrite(pos, lastweight);
pos = mwrite(pos, optionalWeights);
for(int ff = 0; ff < optionalWeightsNames.size(); ff++) {
pos = mwrite(pos, optionalWeightsNames[ff]);
}
pos = mwrite(pos, optionalnpLO);
pos = mwrite(pos, optionalnpNLO);
pos = mwrite(pos, preweight);
cacheFile().write(&buff[0], buff.size(), 1);
}
bool LesHouchesReader::uncacheEvent() {
reset();
static vector<char> buff;
if ( cacheFile().read(&hepeup.NUP, sizeof(hepeup.NUP)) != 1 )
return false;
buff.resize(eventSize(hepeup.NUP));
if ( cacheFile().read(&buff[0], buff.size()) != 1 ) return false;
const char * pos = &buff[0];
pos = mread(pos, hepeup.IDPRUP);
pos = mread(pos, hepeup.XWGTUP);
pos = mread(pos, hepeup.XPDWUP);
pos = mread(pos, hepeup.SCALUP);
pos = mread(pos, hepeup.AQEDUP);
pos = mread(pos, hepeup.AQCDUP);
hepeup.IDUP.resize(hepeup.NUP);
pos = mread(pos, hepeup.IDUP[0], hepeup.NUP);
hepeup.ISTUP.resize(hepeup.NUP);
pos = mread(pos, hepeup.ISTUP[0], hepeup.NUP);
hepeup.MOTHUP.resize(hepeup.NUP);
pos = mread(pos, hepeup.MOTHUP[0], hepeup.NUP);
hepeup.ICOLUP.resize(hepeup.NUP);
pos = mread(pos, hepeup.ICOLUP[0], hepeup.NUP);
hepeup.PUP.resize(hepeup.NUP, vector<double>(5));
for ( int i = 0; i < hepeup.NUP; ++i )
pos = mread(pos, hepeup.PUP[i][0], 5);
hepeup.VTIMUP.resize(hepeup.NUP);
pos = mread(pos, hepeup.VTIMUP[0], hepeup.NUP);
hepeup.SPINUP.resize(hepeup.NUP);
pos = mread(pos, hepeup.SPINUP[0], hepeup.NUP);
pos = mread(pos, lastweight);
pos = mread(pos, optionalWeights);
for(int ff = 0; ff < optionalWeightsNames.size(); ff++) {
pos = mread(pos, optionalWeightsNames[ff]);
}
pos = mread(pos, optionalnpLO);
pos = mread(pos, optionalnpNLO);
pos = mread(pos, preweight);
// If we are skipping, we do not have to do anything else.
if ( skipping ) return true;
if ( CKKWHandler() && maxMultCKKW() > 0 && maxMultCKKW() > minMultCKKW() ) {
// The cached event has not been submitted to CKKW reweighting, so
// we do that now.
fillEvent();
getSubProcess();
CKKWHandler()->setXComb(lastXCombPtr());
lastweight *= CKKWHandler()->reweightCKKW(minMultCKKW(), maxMultCKKW());
}
return true;
}
void LesHouchesReader::persistentOutput(PersistentOStream & os) const {
os << heprup.IDBMUP << heprup.EBMUP << heprup.PDFGUP << heprup.PDFSUP
<< heprup.IDWTUP << heprup.NPRUP << heprup.XSECUP << heprup.XERRUP
<< heprup.XMAXUP << heprup.LPRUP << hepeup.NUP << hepeup.IDPRUP
<< hepeup.XWGTUP << hepeup.XPDWUP << hepeup.SCALUP << hepeup.AQEDUP
<< hepeup.AQCDUP << hepeup.IDUP << hepeup.ISTUP << hepeup.MOTHUP
<< hepeup.ICOLUP << hepeup.PUP << hepeup.VTIMUP << hepeup.SPINUP
<< inData << inPDF << outPDF << thePartonExtractor << theCKKW
<< thePartonBins << theXCombs << theCuts << theNEvents << position
<< reopened << theMaxScan << isActive
<< theCacheFileName << doCutEarly << stats << statmap
<< thePartonBinInstances
<< theBeams << theIncoming << theOutgoing << theIntermediates
<< reweights << preweights << preweight << reweightPDF << doInitPDFs
<< theLastXComb << theMaxMultCKKW << theMinMultCKKW << lastweight << optionalWeights << optionalnpLO << optionalnpNLO
<< maxFactor << ounit(weightScale, picobarn) << xSecWeights << maxWeights
<< theMomentumTreatment << useWeightWarnings << theReOpenAllowed
<< theIncludeSpin;
}
void LesHouchesReader::persistentInput(PersistentIStream & is, int) {
if ( cacheFile() ) closeCacheFile();
is >> heprup.IDBMUP >> heprup.EBMUP >> heprup.PDFGUP >> heprup.PDFSUP
>> heprup.IDWTUP >> heprup.NPRUP >> heprup.XSECUP >> heprup.XERRUP
>> heprup.XMAXUP >> heprup.LPRUP >> hepeup.NUP >> hepeup.IDPRUP
>> hepeup.XWGTUP >> hepeup.XPDWUP >> hepeup.SCALUP >> hepeup.AQEDUP
>> hepeup.AQCDUP >> hepeup.IDUP >> hepeup.ISTUP >> hepeup.MOTHUP
>> hepeup.ICOLUP >> hepeup.PUP >> hepeup.VTIMUP >> hepeup.SPINUP
>> inData >> inPDF >> outPDF >> thePartonExtractor >> theCKKW
>> thePartonBins >> theXCombs >> theCuts >> theNEvents >> position
>> reopened >> theMaxScan >> isActive
>> theCacheFileName >> doCutEarly >> stats >> statmap
>> thePartonBinInstances
>> theBeams >> theIncoming >> theOutgoing >> theIntermediates
>> reweights >> preweights >> preweight >> reweightPDF >> doInitPDFs
>> theLastXComb >> theMaxMultCKKW >> theMinMultCKKW >> lastweight >> optionalWeights >> optionalnpLO >> optionalnpNLO
>> maxFactor >> iunit(weightScale, picobarn) >> xSecWeights >> maxWeights
>> theMomentumTreatment >> useWeightWarnings >> theReOpenAllowed
>> theIncludeSpin;
}
AbstractClassDescription<LesHouchesReader>
LesHouchesReader::initLesHouchesReader;
// Definition of the static class description member.
void LesHouchesReader::setBeamA(long id) { heprup.IDBMUP.first = id; }
long LesHouchesReader::getBeamA() const { return heprup.IDBMUP.first; }
void LesHouchesReader::setBeamB(long id) { heprup.IDBMUP.second = id; }
long LesHouchesReader::getBeamB() const { return heprup.IDBMUP.second; }
void LesHouchesReader::setEBeamA(Energy e) { heprup.EBMUP.first = e/GeV; }
Energy LesHouchesReader::getEBeamA() const { return heprup.EBMUP.first*GeV; }
void LesHouchesReader::setEBeamB(Energy e) { heprup.EBMUP.second = e/GeV; }
Energy LesHouchesReader::getEBeamB() const { return heprup.EBMUP.second*GeV; }
void LesHouchesReader::setPDFA(PDFPtr pdf) { inPDF.first = pdf; }
PDFPtr LesHouchesReader::getPDFA() const { return inPDF.first; }
void LesHouchesReader::setPDFB(PDFPtr pdf) { inPDF.second = pdf; }
PDFPtr LesHouchesReader::getPDFB() const { return inPDF.second; }
void LesHouchesReader::Init() {
static ClassDocumentation<LesHouchesReader> documentation
("ThePEG::LesHouchesReader is an abstract base class to be used "
"for objects which reads event files or streams from matrix element "
"generators.");
static Parameter<LesHouchesReader,long> interfaceBeamA
("BeamA",
"The PDG id of the incoming particle along the positive z-axis. "
"If zero the corresponding information is to be deduced from the "
"event stream/file.",
0, 0, 0, 0,
true, false, false,
&LesHouchesReader::setBeamA,
&LesHouchesReader::getBeamA, 0, 0, 0);
static Parameter<LesHouchesReader,long> interfaceBeamB
("BeamB",
"The PDG id of the incoming particle along the negative z-axis. "
"If zero the corresponding information is to be deduced from the "
"event stream/file.",
0, 0, 0, 0,
true, false, false,
&LesHouchesReader::setBeamB,
&LesHouchesReader::getBeamB, 0, 0, 0);
static Parameter<LesHouchesReader,Energy> interfaceEBeamA
("EBeamA",
"The energy of the incoming particle along the positive z-axis. "
"If zero the corresponding information is to be deduced from the "
"event stream/file.",
0, GeV, ZERO, ZERO, 1000000000.0*GeV,
true, false, true,
&LesHouchesReader::setEBeamA, &LesHouchesReader::getEBeamA, 0, 0, 0);
static Parameter<LesHouchesReader,Energy> interfaceEBeamB
("EBeamB",
"The energy of the incoming particle along the negative z-axis. "
"If zero the corresponding information is to be deduced from the "
"event stream/file.",
0, GeV, ZERO, ZERO, 1000000000.0*GeV,
true, false, true,
&LesHouchesReader::setEBeamB, &LesHouchesReader::getEBeamB, 0, 0, 0);
static Reference<LesHouchesReader,PDFBase> interfacePDFA
("PDFA",
"The PDF used for incoming particle along the positive z-axis. "
"If null the corresponding information is to be deduced from the "
"event stream/file.",
0, true, false, true, true, false,
&LesHouchesReader::setPDFA, &LesHouchesReader::getPDFA, 0);
static Reference<LesHouchesReader,PDFBase> interfacePDFB
("PDFB",
"The PDF used for incoming particle along the negative z-axis. "
"If null the corresponding information is to be deduced from the "
"event stream/file.",
0, true, false, true, true, false,
&LesHouchesReader::setPDFB, &LesHouchesReader::getPDFB, 0);
static Parameter<LesHouchesReader,long> interfaceMaxScan
("MaxScan",
"The maximum number of events to scan to obtain information about "
"processes and cross section in the intialization.",
&LesHouchesReader::theMaxScan, -1, 0, 0,
true, false, false);
static Parameter<LesHouchesReader,string> interfaceCacheFileName
("CacheFileName",
- "Name of file used to cache the events form the reader in a fast-readable "
+ "Name of file used to cache the events from the reader in a fast-readable "
"form. If empty, no cache file will be generated.",
&LesHouchesReader::theCacheFileName, "",
true, false);
interfaceCacheFileName.fileType();
static Switch<LesHouchesReader,bool> interfaceCutEarly
("CutEarly",
"Determines whether to apply cuts to events before converting to "
"ThePEG format.",
&LesHouchesReader::doCutEarly, true, true, false);
static SwitchOption interfaceCutEarlyYes
(interfaceCutEarly,
"Yes",
"Event are cut before converted.",
true);
static SwitchOption interfaceCutEarlyNo
(interfaceCutEarly,
"No",
"Events are not cut before converted.",
false);
static Reference<LesHouchesReader,PartonExtractor> interfacePartonExtractor
("PartonExtractor",
"The PartonExtractor object used to construct remnants. If no object is "
"provided the LesHouchesEventHandler object must provide one instead.",
&LesHouchesReader::thePartonExtractor, true, false, true, true, false);
static Reference<LesHouchesReader,Cuts> interfaceCuts
("Cuts",
"The Cuts object to be used for this reader. Note that these "
"must not be looser cuts than those used in the actual generation. "
"If no object is provided the LesHouchesEventHandler object must "
"provide one instead.",
&LesHouchesReader::theCuts, true, false, true, true, false);
static RefVector<LesHouchesReader,ReweightBase> interfaceReweights
("Reweights",
"A list of ThePEG::ReweightBase objects to modify this the weight of "
"this reader.",
&LesHouchesReader::reweights, 0, false, false, true, false);
static RefVector<LesHouchesReader,ReweightBase> interfacePreweights
("Preweights",
"A list of ThePEG::ReweightBase objects to bias the phase space for this "
"reader without influencing the actual cross section.",
&LesHouchesReader::preweights, 0, false, false, true, false);
static Switch<LesHouchesReader,bool> interfaceReweightPDF
("ReweightPDF",
"If the PDFs used in the generation for this reader is different "
"from the ones assumed by the associated PartonExtractor object, "
"should the events be reweighted to fit the latter?",
&LesHouchesReader::reweightPDF, false, true, false);
static SwitchOption interfaceReweightPDFNo
(interfaceReweightPDF, "No", "The event weights are kept as they are.",
false);
static SwitchOption interfaceReweightPDFYes
(interfaceReweightPDF,
"Yes", "The events are reweighted.", true);
static Switch<LesHouchesReader,bool> interfaceInitPDFs
("InitPDFs",
"If no PDFs were specified in <interface>PDFA</interface> or "
"<interface>PDFB</interface>for this reader, try to extract the "
"information from the event file and assign the relevant PDFBase"
"objects when the reader is initialized.",
&LesHouchesReader::doInitPDFs, false, true, false);
static SwitchOption interfaceInitPDFsYes
(interfaceInitPDFs,
"Yes",
"Extract PDFs during initialization.",
true);
static SwitchOption interfaceInitPDFsNo
(interfaceInitPDFs,
"No",
"Do not extract PDFs during initialization.",
false);
static Parameter<LesHouchesReader,int> interfaceMaxMultCKKW
("MaxMultCKKW",
"If this reader is to be used (possibly together with others) for CKKW-"
"reweighting and veto, this should give the multiplicity of outgoing "
"particles in the highest multiplicity matrix element in the group. "
"If set to zero, no CKKW procedure should be applied.",
&LesHouchesReader::theMaxMultCKKW, 0, 0, 0,
true, false, Interface::lowerlim);
static Parameter<LesHouchesReader,int> interfaceMinMultCKKW
("MinMultCKKW",
"If this reader is to be used (possibly together with others) for CKKW-"
"reweighting and veto, this should give the multiplicity of outgoing "
"particles in the lowest multiplicity matrix element in the group. If "
"larger or equal to <interface>MaxMultCKKW</interface>, no CKKW "
"procedure should be applied.",
&LesHouchesReader::theMinMultCKKW, 0, 0, 0,
true, false, Interface::lowerlim);
static Switch<LesHouchesReader,unsigned int> interfaceMomentumTreatment
("MomentumTreatment",
"Treatment of the momenta supplied by the interface",
&LesHouchesReader::theMomentumTreatment, 0, false, false);
static SwitchOption interfaceMomentumTreatmentAccept
(interfaceMomentumTreatment,
"Accept",
"Just accept the momenta given",
0);
static SwitchOption interfaceMomentumTreatmentRescaleEnergy
(interfaceMomentumTreatment,
"RescaleEnergy",
"Rescale the energy supplied so it is consistent with the mass",
1);
static SwitchOption interfaceMomentumTreatmentRescaleMass
(interfaceMomentumTreatment,
"RescaleMass",
"Rescale the mass supplied so it is consistent with the"
" energy and momentum",
2);
static Switch<LesHouchesReader,bool> interfaceWeightWarnings
("WeightWarnings",
"Determines if warnings about possible weight incompatibilities should "
"be issued when this reader is initialized.",
&LesHouchesReader::useWeightWarnings, true, true, false);
static SwitchOption interfaceWeightWarningsWarnAboutWeights
(interfaceWeightWarnings,
"WarnAboutWeights",
"Warn about possible incompatibilities with the weight option in the "
"Les Houches common block and the requested weight treatment.",
true);
static SwitchOption interfaceWeightWarningsDontWarnAboutWeights
(interfaceWeightWarnings,
"DontWarnAboutWeights",
"Do not warn about possible incompatibilities with the weight option "
"in the Les Houches common block and the requested weight treatment.",
false);
static Switch<LesHouchesReader,bool> interfaceAllowedTopReOpen
("AllowedToReOpen",
"Can the file be reopened if more events are requested than the file contains?",
&LesHouchesReader::theReOpenAllowed, true, false, false);
static SwitchOption interfaceAllowedTopReOpenYes
(interfaceAllowedTopReOpen,
"Yes",
"Allowed to reopen the file",
true);
static SwitchOption interfaceAllowedTopReOpenNo
(interfaceAllowedTopReOpen,
"No",
"Not allowed to reopen the file",
false);
static Switch<LesHouchesReader,bool> interfaceIncludeSpin
("IncludeSpin",
"Use the spin information present in the event file, for tau leptons"
" only as this is the only case which makes any sense",
&LesHouchesReader::theIncludeSpin, true, false, false);
static SwitchOption interfaceIncludeSpinYes
(interfaceIncludeSpin,
"Yes",
"Use the spin information",
true);
static SwitchOption interfaceIncludeSpinNo
(interfaceIncludeSpin,
"No",
"Don't use the spin information",
false);
interfaceCuts.rank(8);
interfacePartonExtractor.rank(7);
interfaceBeamA.rank(5);
interfaceBeamB.rank(4);
interfaceEBeamA.rank(3);
interfaceEBeamB.rank(2);
interfaceMaxMultCKKW.rank(1.5);
interfaceMinMultCKKW.rank(1.0);
interfaceBeamA.setHasDefault(false);
interfaceBeamB.setHasDefault(false);
interfaceEBeamA.setHasDefault(false);
interfaceEBeamB.setHasDefault(false);
interfaceMaxMultCKKW.setHasDefault(false);
interfaceMinMultCKKW.setHasDefault(false);
}
namespace ThePEG {
ostream & operator<<(ostream & os, const HEPEUP & h) {
os << "<event>\n"
<< " " << setw(4) << h.NUP
<< " " << setw(6) << h.IDPRUP
<< " " << setw(14) << h.XWGTUP
<< " " << setw(14) << h.SCALUP
<< " " << setw(14) << h.AQEDUP
<< " " << setw(14) << h.AQCDUP << "\n";
for ( int i = 0; i < h.NUP; ++i )
os << " " << setw(8) << h.IDUP[i]
<< " " << setw(2) << h.ISTUP[i]
<< " " << setw(4) << h.MOTHUP[i].first
<< " " << setw(4) << h.MOTHUP[i].second
<< " " << setw(4) << h.ICOLUP[i].first
<< " " << setw(4) << h.ICOLUP[i].second
<< " " << setw(14) << h.PUP[i][0]
<< " " << setw(14) << h.PUP[i][1]
<< " " << setw(14) << h.PUP[i][2]
<< " " << setw(14) << h.PUP[i][3]
<< " " << setw(14) << h.PUP[i][4]
<< " " << setw(1) << h.VTIMUP[i]
<< " " << setw(1) << h.SPINUP[i] << std::endl;
os << "</event>" << std::endl;
return os;
}
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 2:00 PM (11 h, 43 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023051
Default Alt Text
(88 KB)

Event Timeline