diff --git a/MatrixElement/FxFx/FxFxEventHandler.cc b/MatrixElement/FxFx/FxFxEventHandler.cc
--- a/MatrixElement/FxFx/FxFxEventHandler.cc
+++ b/MatrixElement/FxFx/FxFxEventHandler.cc
@@ -1,652 +1,673 @@
// -*- C++ -*-
//
// Based on:
// FxFxEventHandler.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2019 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the FxFxEventHandler class.
//
#include "FxFxEventHandler.h"
#include "FxFxReader.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/Handlers/LuminosityFunction.h"
#include "ThePEG/Handlers/XComb.h"
#include "ThePEG/Handlers/CascadeHandler.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/PDF/PartonExtractor.h"
#include "ThePEG/Utilities/LoopGuard.h"
#include "ThePEG/EventRecord/Event.h"
#include "ThePEG/EventRecord/Collision.h"
#include "ThePEG/EventRecord/Step.h"
#include "ThePEG/EventRecord/SubProcess.h"
#include "ThePEG/Utilities/EnumIO.h"
#include "ThePEG/Utilities/Maths.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Utilities/Debug.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace ThePEG;
FxFxEventHandler::~FxFxEventHandler() {}
IBPtr FxFxEventHandler::clone() const {
return new_ptr(*this);
}
IBPtr FxFxEventHandler::fullclone() const {
return new_ptr(*this);
}
void FxFxEventHandler::doinit() {
EventHandler::doinit();
for ( int i = 0, N = readers().size(); i < N; ++i ) {
readers()[i]->init();
//FxFxReader & reader = *readers()[i];
//reader.initialize(*this);
//weightnames = reader.optWeightsNamesFunc();
}
/*
XSecStat* initxsecs = new XSecStat[weightnames.size()];
for(int ww = 0; ww < weightnames.size(); ww++){
initxsecs[ww].reset();
optstats.insert(std::make_pair<string,XSecStat>(weightnames[ww], initxsecs[ww]));
opthistStats.insert(std::make_pair<string,XSecStat>(weightnames[ww], initxsecs[ww]));
CrossSection initxs = 0.*picobarn;
optxs.insert(std::make_pair<string,CrossSection>(weightnames[ww], initxs));
}*/
ntries = 0;
}
void FxFxEventHandler::initialize() {
if ( lumiFnPtr() ) Repository::clog()
<< "The LuminosityFunction '" << lumiFnPtr()->name()
<< "' assigned to the FxFxEventHandler '" << name()
<< "' will not be active in this run. Instead the incoming "
<< "particles will be determined by the used FxFxReader objects.\n"
<< Exception::warning;
if ( readers().empty() )
throw FxFxInitError()
<< "No readers were defined for the FxFxEventHandler '"
<< name() << "'" << Exception::warning;
// Go through all the readers and collect information about cross
// sections and processes.
typedef map<int,tFxFxReaderPtr> ProcessMap;
ProcessMap processes;
PDPair incoming;
Energy MaxEA = ZERO;
Energy MaxEB = ZERO;
for ( int i = 0, N = readers().size(); i < N; ++i ) {
FxFxReader & reader = *readers()[i];
reader.initialize(*this);
weightnames = reader.optWeightsNamesFunc();
// Check that the incoming particles are consistent between the
// readers.
if ( !incoming.first ) {
incoming.first = getParticleData(reader.heprup.IDBMUP.first);
if ( !incoming.first )
Throw<FxFxInitError>()
<< "Unknown beam PID " << reader.heprup.IDBMUP.first
<< ". Have you created a matching BeamParticle object?"
<< Exception::runerror;
}
if ( !incoming.second ) {
incoming.second = getParticleData(reader.heprup.IDBMUP.second);
if ( !incoming.second )
Throw<FxFxInitError>()
<< "Unknown beam PID " << reader.heprup.IDBMUP.first
<< ". Have you created a matching BeamParticle object?"
<< Exception::runerror;
}
if ( incoming.first->id() != reader.heprup.IDBMUP.first ||
incoming.second->id() != reader.heprup.IDBMUP.second )
Repository::clog()
<< "The different FxFxReader objects in the "
<< "FxFxEventHandler '" << name() << "' have different "
<< "types of colliding particles." << Exception::warning;
MaxEA = max(MaxEA, reader.heprup.EBMUP.first*GeV);
MaxEB = max(MaxEB, reader.heprup.EBMUP.second*GeV);
// Check that the weighting of the events in the different readers
// is consistent with the ones requested for this event
// handler. Also collect the sum of the maximum weights.
if ( reader.negativeWeights() && weightOption() > 0 )
throw FxFxInitError()
<< "The reader '" << reader.name()
<< "' contains negatively weighted events, "
<< "which is not allowed for the FxFxEventHandler '"
<< name() << "'." << Exception::warning;
// Check that we do not have the same process numbers in different
// readers.
for ( int ip = 0; ip < reader.heprup.NPRUP; ++ip ) {
if ( reader.heprup.LPRUP[ip] ) {
ProcessMap::iterator pit = processes.find(reader.heprup.LPRUP[ip]);
if ( pit == processes.end() )
processes[reader.heprup.LPRUP[ip]] = readers()[i];
else if ( warnPNum ) {
Throw<FxFxPNumException>()
<< "In the FxFxEventHandler '"
<< name() << "', both the '" << pit->second->name() << "' and '"
<< reader.name() << "' contains sub-process number " << pit->first
<< ". This process may be double-counted in this run."
<< Exception::warning;
}
}
}
selector().insert(reader.stats.maxXSec(), i);
}
stats.maxXSec(selector().sum());
histStats.maxXSec(selector().sum());
for (map<string,XSecStat>::iterator it= optstats.begin(); it!=optstats.end(); ++it){
(it->second).maxXSec(selector().sum());
}
// Check that we have any cross section at all.
if ( stats.maxXSec() <= ZERO )
throw FxFxInitError()
<< "The sum of the cross sections of the readers in the "
<< "FxFxEventHandler '" << name()
<< "' was zero." << Exception::warning;
// We now create a LuminosityFunction object to inform others about
// the energy of the beam.
theIncoming = incoming;
lumiFn(new_ptr(LuminosityFunction(MaxEA, MaxEB)));
}
void FxFxEventHandler::doinitrun() {
EventHandler::doinitrun();
stats.reset();
histStats.reset();
for ( int i = 0, N = readers().size(); i < N; ++i ) {
readers()[i]->init();
FxFxReader & reader = *readers()[i];
reader.initialize(*this);
weightnames = reader.optWeightsNamesFunc();
}
// XSecStat initxsecs;
XSecStat* initxsecs = new XSecStat[weightnames.size()];
for(size_t ww = 0; ww < weightnames.size(); ww++){
initxsecs[ww].reset();
// optstats.insert(std::make_pair<string,XSecStat>(weightnames[ww], initxsecs[ww]));
// opthistStats.insert(std::make_pair<string,XSecStat>(weightnames[ww], initxsecs[ww]));
CrossSection initxs = 0.*picobarn;
//optxs.insert(std::make_pair<string,CrossSection>(weightnames[ww], initxs));
optstats.insert(std::make_pair(weightnames[ww], initxsecs[ww]));
opthistStats.insert(std::make_pair(weightnames[ww], initxsecs[ww]));
optxs.insert(std::make_pair(weightnames[ww], initxs));
}
ntries = 0;
}
EventPtr FxFxEventHandler::generateEvent() {
LoopGuard<EventLoopException,FxFxEventHandler>
loopGuard(*this, maxLoop());
while ( true ) {
loopGuard();
currentReader(readers()[selector().select(UseRandom::current())]);
skipEvents();
currentReader()->reset();
double weight = currentReader()->getEvent();
if ( weightOption() == unitweight && weight < 0.0 ) weight = 0.0;
if ( weightOption() == unitweight || weightOption() == unitnegweight ) {
CrossSection newmax = selector().reweight(weight);
if ( newmax > CrossSection() )
increaseMaxXSec(newmax);
}
select(weight/currentReader()->preweight);
histStats.select(weight);
if ( !weighted() ) {
if ( weightOption() == unitweight || weightOption() == unitnegweight ) {
if ( !rndbool(abs(weight)) ) continue;
weight = Math::sign(1.0, weight);
}
else if ( weight == 0.0 ) continue;
} else if ( weight == 0.0 ) continue;
accept();
// Divide by the bias introduced by the preweights in the reader.
weight /= currentReader()->preweight;
// fact for weight normalization
double fact = theNormWeight ? double(selector().sum()/picobarn) : 1.;
try {
theLastXComb = currentReader()->getXComb();
- currentEvent(new_ptr(Event(lastParticles(), this, generator()->runName(),
- generator()->currentEventNumber(), weight*fact)));
+ //whether to use the LHE event number or not for the event identification
+ if(UseLHEEvent==0 || currentReader()->LHEEventNum() == -1) {
+ currentEvent(new_ptr(Event(lastParticles(), this, generator()->runName(),
+ generator()->currentEventNumber(), weight*fact )));
+ }
+ else if(UseLHEEvent==1 && currentReader()->LHEEventNum() != -1) {
+ currentEvent(new_ptr(Event(lastParticles(), this, generator()->runName(),
+ currentReader()->LHEEventNum(), weight*fact )));
+ }
currentEvent()->optionalWeights() = currentReader()->optionalEventWeights();
// normalize the optional weights
for(map<string,double>::iterator it = currentEvent()->optionalWeights().begin();
it!=currentEvent()->optionalWeights().end();++it) {
if(it->first!="ecom"&& it->second!=-999 && it->second!=-111 && it->second!=-222 && it->second!=-333) { it->second *= fact; }
}
//print optional weights here
// cout << "event weight = " << weight << " fact = " << fact << endl;
/*for (map<string,double>::const_iterator it= currentReader()->optionalEventWeights().begin(); it!=currentReader()->optionalEventWeights().end(); ++it){
std::cout << it->first << " => " << it->second << '\n';
}
cout << endl;*/
//print npLO and npNLO
// cout << currentReader()->optionalEventnpLO() << "\t" << currentReader()->optionalEventnpNLO() << endl;
performCollision();
if ( !currentCollision() ) throw Veto();
return currentEvent();
}
catch (Veto) {
reject(weight);
}
catch (Stop) {
}
catch (Exception &) {
reject(weight);
throw;
}
}
}
void FxFxEventHandler::skipEvents() {
if ( weightOption() == 2 || weightOption() == -2 ) return; //does it make sense to skip events if we are using varying weights?
// Don't do this for readers which seem to generate events on the fly.
if ( currentReader()->active() || currentReader()->NEvents() <= 0 ) return;
// Estimate the fration of the total events available from
// currentReader() which will be requested.
double frac = currentReader()->stats.maxXSec()/stats.maxXSec();
if ( stats.accepted() > 0 )
frac *= double(stats.attempts())/double(stats.accepted());
else
frac *= double(stats.attempts() + 1);
double xscan = generator()->N()*frac/currentReader()->NEvents();
// Estimate the number of times we need to go through the events for
// the currentReader(), and how many events on average we need to
// skip for each attempted event to go through the file an integer
// number of times.
double nscan = ceil(xscan);
double meanskip = nscan/xscan - 1.0;
// Skip an average numer of steps with a Poissonian distribution.
currentReader()->
skip(UseRandom::rndPoisson(meanskip)%currentReader()->NEvents());
}
void FxFxEventHandler::select(double weight) {
stats.select(weight);
currentReader()->select(weight);
vector<double> w;
for (map<string,double>::const_iterator it= currentReader()->optionalEventWeights().begin(); it!=currentReader()->optionalEventWeights().end(); ++it){
w.push_back(it->second);
}
int ii = 0;
for (map<string,XSecStat>::iterator it= opthistStats.begin(); it!=opthistStats.end(); ++it){
(it->second).select(w[ii]);
ii++;
}
ii = 0;
for (map<string,XSecStat>::iterator it= optstats.begin(); it!=optstats.end(); ++it){
(it->second).select(w[ii]);
ii++;
}
}
tCollPtr FxFxEventHandler::performCollision() {
lastExtractor()->select(lastXCombPtr());
if ( CKKWHandler() ) CKKWHandler()->setXComb(lastXCombPtr());
currentCollision(new_ptr(Collision(lastParticles(), currentEvent(), this)));
if ( currentEvent() ) currentEvent()->addCollision(currentCollision());
currentStep(new_ptr(Step(currentCollision(), this)));
currentCollision()->addStep(currentStep());
currentStep()->addSubProcess(currentReader()->getSubProcess());
lastExtractor()->constructRemnants(lastXCombPtr()->partonBinInstances(),
subProcess(), currentStep());
if ( !currentReader()->cuts().passCuts(*currentCollision()) ) throw Veto();
initGroups();
if ( ThePEG_DEBUG_ITEM(1) ) {
if ( currentEvent() )
generator()->logfile() << *currentEvent();
else
generator()->logfile() << *currentCollision();
}
return continueCollision();
}
EventPtr FxFxEventHandler::continueEvent() {
try {
continueCollision();
}
catch (Veto) {
const double fact =
theNormWeight ?
double(selector().sum()/picobarn) : 1.;
reject(currentEvent()->weight()/fact);
}
catch (Stop) {
}
catch (Exception &) {
const double fact =
theNormWeight ?
double(selector().sum()/picobarn) : 1.;
reject(currentEvent()->weight()/fact);
throw;
}
return currentEvent();
}
void FxFxEventHandler::dofinish() {
EventHandler::dofinish();
if ( selector().compensating() ) generator()->log()
<< "Warning: The run was ended while the FxFxEventHandler '"
<< name() << "' was still trying to compensate for weights larger than 1. "
<< "The cross section estimates may therefore be statistically "
<< "inaccurate." << endl;
}
void FxFxEventHandler::statistics(ostream & os) const {
if ( statLevel() == 0 ) return;
string line = "======================================="
"=======================================\n";
if ( stats.accepted() <= 0 ) {
os << line << "No events generated by event handler '" << name() << "'."
<< endl;
return;
}
os << line << "Statistics for Les Houches event handler \'" << name() << "\':\n"
<< " "
<< "generated number of Cross-section\n"
<< " "
<< " events attempts (nb)\n";
os << line << "Total:" << setw(42) << stats.accepted() << setw(13)
<< stats.attempts() << setw(17)
<< ouniterr(stats.xSec(), stats.xSecErr(), nanobarn) << endl
<< line;
if ( statLevel() == 1 ) return;
if ( statLevel() == 2 ) {
os << "Per Les Houches Reader breakdown:\n";
for ( int i = 0, N = readers().size(); i < N; ++i ) {
FxFxReader & reader = *readers()[i];
string n = reader.name();
n.resize(37, ' ');
os << n << setw(11) << reader.stats.accepted() << setw(13)
<< reader.stats.attempts() << setw(17)
<< ouniterr(reader.stats.xSec(), reader.stats.xSecErr(), nanobarn)
<< endl;
}
os << line;
} else {
os << "Per Les Houches Reader (and process #) breakdown:\n";
for ( int i = 0, N = readers().size(); i < N; ++i ) {
FxFxReader & reader = *readers()[i];
string n = reader.name() + " (all)";
n.resize(37, ' ');
os << n << setw(11) << reader.stats.accepted() << setw(13)
<< reader.stats.attempts() << setw(17)
<< ouniterr(reader.stats.xSec(), reader.stats.xSecErr(), nanobarn)
<< endl;
CrossSection xsectot = reader.stats.xSec();
if ( xsectot != ZERO ) xsectot /= reader.stats.sumWeights();
typedef FxFxReader::StatMap::const_iterator const_iterator;
for ( const_iterator i = reader.statmap.begin();
i != reader.statmap.end(); ++i ) {
ostringstream ss;
ss << reader.name() << " (" << i->first << ")";
string n = ss.str();
n.resize(37, ' ');
os << n << setw(11) << i->second.accepted() << setw(13)
<< i->second.attempts() << setw(17)
<< ouniterr(i->second.sumWeights()*xsectot,
sqrt(i->second.sumWeights2())*xsectot, nanobarn) << endl;
}
os << line;
}
}
string warn = "Warning: Result may be statistically incorrect since\n"
" the following FxFxReaders were oversampled:\n";
for ( int i = 0, N = readers().size(); i < N; ++i ) {
FxFxReader & reader = *readers()[i];
if ( reader.NEvents() > 0 && reader.stats.attempts() > reader.NEvents() ) {
os << warn;
warn = "";
os << "'" << reader.name() << "' (by a factor "
<< double(reader.stats.attempts())/double(reader.NEvents())
<< ")" << endl;
}
}
}
void FxFxEventHandler::increaseMaxXSec(CrossSection maxxsec) {
stats.maxXSec(selector().sum());
histStats.maxXSec(selector().sum());
currentReader()->increaseMaxXSec(maxxsec);
}
void FxFxEventHandler::accept() {
ntries++;
stats.accept();
histStats.accept();
currentReader()->accept();
for (map<string,XSecStat>::iterator it= opthistStats.begin(); it!=opthistStats.end(); ++it){
(it->second).accept();
}
for (map<string,XSecStat>::iterator it= optstats.begin(); it!=optstats.end(); ++it){
(it->second).accept();
}
}
void FxFxEventHandler::reject(double w) {
ntries++;
stats.reject(w);
histStats.reject(w);
currentReader()->reject(w);
vector<double> wv;
for (map<string,double>::const_iterator it= currentReader()->optionalEventWeights().begin(); it!=currentReader()->optionalEventWeights().end(); ++it){
wv.push_back(it->second);
}
int ii = 0;
for (map<string,XSecStat>::iterator it= opthistStats.begin(); it!=opthistStats.end(); ++it){
(it->second).reject(wv[ii]);
ii++;
}
ii = 0;
for (map<string,XSecStat>::iterator it= optstats.begin(); it!=optstats.end(); ++it){
(it->second).reject(wv[ii]);
ii++;
}
}
map<string,CrossSection> FxFxEventHandler::optintegratedXSecMap() const {
map<string,CrossSection> result;
for (map<string,XSecStat>::const_iterator it= optstats.begin(); it!=optstats.end(); ++it){
result[it->first] = (it->second.sumWeights()/it->second.attempts()) * picobarn;
}
return result;
}
CrossSection FxFxEventHandler::histogramScale() const {
return histStats.xSec()/histStats.sumWeights();
}
CrossSection FxFxEventHandler::integratedXSec() const {
return histStats.xSec();
}
CrossSection FxFxEventHandler::integratedXSecErr() const {
return histStats.xSecErr();
}
int FxFxEventHandler::ntriesinternal() const {
return stats.attempts();
}
void FxFxEventHandler::persistentOutput(PersistentOStream & os) const {
os << stats << histStats << theReaders << theSelector
<< oenum(theWeightOption) << theUnitTolerance << theCurrentReader << warnPNum
- << theNormWeight;
+ << theNormWeight << UseLHEEvent;
}
void FxFxEventHandler::persistentInput(PersistentIStream & is, int) {
is >> stats >> histStats >> theReaders >> theSelector
>> ienum(theWeightOption) >> theUnitTolerance >> theCurrentReader >> warnPNum
- >> theNormWeight;
+ >> theNormWeight >> UseLHEEvent;;
}
ClassDescription<FxFxEventHandler>
FxFxEventHandler::initFxFxEventHandler;
// Definition of the static class description member.
void FxFxEventHandler::setUnitTolerance(double x) {
theUnitTolerance = x;
selector().tolerance(unitTolerance());
}
void FxFxEventHandler::Init() {
static ClassDocumentation<FxFxEventHandler> documentation
("This is the main class administrating the selection of hard "
"subprocesses from a set of ThePEG::FxFxReader objects.");
static RefVector<FxFxEventHandler,FxFxReader>
interfaceFxFxReaders
("FxFxReaders",
"Objects capable of reading events from an event file or an "
"external matrix element generator.",
&FxFxEventHandler::theReaders, -1, false, false, true, false, false);
static Switch<FxFxEventHandler,WeightOpt> interfaceWeightOption
("WeightOption",
"The different ways to weight events in the Les Houches event handler. "
"Whether weighted or not and whether or not negative weights are allowed.",
&FxFxEventHandler::theWeightOption, unitweight, true, false);
static SwitchOption interfaceWeightOptionUnitWeight
(interfaceWeightOption,
"UnitWeight",
"All events have unit weight.",
unitweight);
static SwitchOption interfaceWeightOptionNegUnitWeight
(interfaceWeightOption,
"NegUnitWeight",
"All events have weight +1 or maybe -1.",
unitnegweight);
static SwitchOption interfaceWeightOptionVarWeight
(interfaceWeightOption,
"VarWeight",
"Events may have varying but positive weights.",
varweight);
static SwitchOption interfaceWeightOptionVarNegWeight
(interfaceWeightOption,
"VarNegWeight",
"Events may have varying weights, both positive and negative.",
varnegweight);
static Switch<FxFxEventHandler,bool> interfaceWarnPNum
("WarnPNum",
"Warn if the same process number is used in more than one "
"FxFxReader.",
&FxFxEventHandler::warnPNum, true, true, false);
static SwitchOption interfaceWarnPNumWarning
(interfaceWarnPNum,
"Warning",
"Give a warning message.",
true);
static SwitchOption interfaceWarnPNumNoWarning
(interfaceWarnPNum,
"NoWarning",
"Don't give a warning message.",
false);
static Parameter<FxFxEventHandler,double> interfaceUnitTolerance
("UnitTolerance",
"If the <interface>WeightOption</interface> is set to unit weight, do not start compensating unless the a weight is found to be this much larger than unity.",
&FxFxEventHandler::theUnitTolerance, 1.0e-6, 0.0, 0,
true, false, Interface::lowerlim,
&FxFxEventHandler::setUnitTolerance,
(double(FxFxEventHandler::*)()const)(0),
(double(FxFxEventHandler::*)()const)(0),
(double(FxFxEventHandler::*)()const)(0),
(double(FxFxEventHandler::*)()const)(0));
static Switch<FxFxEventHandler,unsigned int> interfaceWeightNormalization
("WeightNormalization",
"How to normalize the output weights",
&FxFxEventHandler::theNormWeight, 0, false, false);
static SwitchOption interfaceWeightNormalizationUnit
(interfaceWeightNormalization,
"Normalized",
"Standard normalization, i.e. +/- for unweighted events",
0);
static SwitchOption interfaceWeightNormalizationCrossSection
(interfaceWeightNormalization,
"CrossSection",
"Normalize the weights to the max cross section in pb",
1);
+ static Switch<FxFxEventHandler,unsigned int> interfaceEventNumbering
+ ("EventNumbering",
+ "How to number the events",
+ &FxFxEventHandler::UseLHEEvent, 0, false, false);
+ static SwitchOption interfaceEventNumberingIncremental
+ (interfaceEventNumbering,
+ "Incremental",
+ "Standard incremental numbering (i.e. as they are generated)",
+ 0);
+ static SwitchOption interfaceEventNumberingLHE
+ (interfaceEventNumbering,
+ "LHE",
+ "Corresponding to the LHE event number",
+ 1);
interfaceFxFxReaders.rank(10);
interfaceWeightOption.rank(9);
}
diff --git a/MatrixElement/FxFx/FxFxEventHandler.h b/MatrixElement/FxFx/FxFxEventHandler.h
--- a/MatrixElement/FxFx/FxFxEventHandler.h
+++ b/MatrixElement/FxFx/FxFxEventHandler.h
@@ -1,447 +1,452 @@
// -*- C++ -*-
//
// FxFxEventHandler.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2019 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef THEPEG_FxFxEventHandler_H
#define THEPEG_FxFxEventHandler_H
//
// This is the declaration of the FxFxEventHandler class.
//
#include "ThePEG/Handlers/EventHandler.h"
#include "FxFxEventHandler.fh"
#include "FxFxReader.fh"
#include "ThePEG/Utilities/CompSelector.h"
#include "ThePEG/Utilities/XSecStat.h"
namespace ThePEG {
/**
* The FxFxEventHandler inherits from the general EventHandler
* class and administers the reading of events generated by external
* matrix element generator programs according to the Les Houches
* accord.
*
* The class has a list of <code>FxFxReader</code>s which
* typically are connected to files with event data produced by
* external matrix element generator programs. When an event is
* requested by FxFxEventHandler, one of the readers are chosen,
* an event is read in and then passed to the different
* <code>StepHandler</code> defined in the underlying
* EventHandler class.
*
* @see \ref FxFxEventHandlerInterfaces "The interfaces"
* defined for FxFxEventHandler.
*/
class FxFxEventHandler: public EventHandler {
public:
/**
* A vector of FxFxReader objects.
*/
typedef vector<FxFxReaderPtr> ReaderVector;
/**
* A selector of readers.
*/
typedef CompSelector<int,CrossSection> ReaderSelector;
/**
* Enumerate the weighting options.
*/
enum WeightOpt {
unitweight = 1, /**< All events have unit weight. */
unitnegweight = -1, /**< All events have wight +/- 1. */
varweight = 2, /**< Varying positive weights. */
varnegweight = -2 /**< Varying positive or negative weights. */
};
friend class FxFxHandler;
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
FxFxEventHandler()
- : theWeightOption(unitweight), theUnitTolerance(1.0e-6), warnPNum(true), theNormWeight(0)
+ : theWeightOption(unitweight), theUnitTolerance(1.0e-6), warnPNum(true), theNormWeight(0), UseLHEEvent(0)
{
selector().tolerance(unitTolerance());
}
/**
* The destructor.
*/
virtual ~FxFxEventHandler();
//@}
public:
/** @name Initialization and finalization functions. */
//@{
/**
* Initialize this event handler and all related objects needed to
* generate events.
*/
virtual void initialize();
/**
* Write out accumulated statistics about intergrated cross sections
* and stuff.
*/
virtual void statistics(ostream &) const;
/**
* Histogram scale. A histogram bin which has been filled with the
* weights associated with the Event objects should be scaled by
* this factor to give the correct cross section.
*/
virtual CrossSection histogramScale() const;
/**
* The estimated total integrated cross section of the processes
* generated in this run.
* @return 0 if no integrated cross section could be estimated.
*/
virtual CrossSection integratedXSec() const;
virtual int ntriesinternal() const;
/**
* The estimated error in the total integrated cross section of the
* processes generated in this run.
* @return 0 if no integrated cross section error could be estimated.
*/
virtual CrossSection integratedXSecErr() const;
virtual map<string,CrossSection> optintegratedXSecMap() const;
//@}
/** @name Functions used for the actual generation */
//@{
/**
* Generate an event.
*/
virtual EventPtr generateEvent();
/**
* Create the Event and Collision objects. Used by the
* generateEvent() function.
*/
virtual tCollPtr performCollision();
/**
* Continue generating an event if the generation has been stopped
* before finishing.
*/
virtual EventPtr continueEvent();
//@}
/** @name Functions to manipulate statistics. */
//@{
/**
* An event has been selected. Signal that an event has been
* selected with the given \a weight. If unit weights are requested,
* the event will be accepted with that weight. This also takes care
* of the statistics collection of the selected reader object.
*/
void select(double weight);
/**
* Accept the current event, taking care of the statistics
* collection of the corresponding reader objects.
*/
void accept();
/**
* Reject the current event, taking care of the statistics
* collection of the corresponding reader objects.
*/
void reject(double weight);
/**
* Increase the overestimated cross section for the selected reader.
*/
void increaseMaxXSec(CrossSection maxxsec);
/**
* Skip some events. To ensure a reader file is scanned an even
* number of times, skip a number of events for the selected reader.
*/
void skipEvents();
//@}
/** @name Simple access functions. */
//@{
/**
* The way weights are to be treated.
*/
WeightOpt weightOption() const { return theWeightOption; }
/**
* If the weight option is set to unit weight, do not start
* compensating unless the weight is this much larger than unity.
*/
double unitTolerance() const { return theUnitTolerance; }
/**
* Access the list of readers.
*/
const ReaderVector & readers() const { return theReaders; }
/**
* The selector to choose readers according to their overestimated
* cross section.
*/
const ReaderSelector & selector() const { return theSelector; }
/**
* The currently selected reader object.
*/
tFxFxReaderPtr currentReader() const { return theCurrentReader; }
/**
* Set the currently selected reader object.
*/
void currentReader(tFxFxReaderPtr x) { theCurrentReader = x; }
//@}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
/**
* The currently selected reader object.
*/
tFxFxReaderPtr theCurrentReader;
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
/**
* Finalize this object. Called in the run phase just after a
* run has ended. Used eg. to write out statistics.
*/
virtual void dofinish();
//@}
protected:
/**
* Access the list of readers.
*/
ReaderVector & readers() { return theReaders; }
/**
* The selector to choose readers according to their overestimated
* cross section.
*/
ReaderSelector & selector() { return theSelector; }
/**
* Helper function for the interface;
*/
void setUnitTolerance(double);
/**
* Collect statistics for this event handler.
*/
XSecStat stats;
map<string,XSecStat> optstats;
map<string,CrossSection> optxs;
int ntries;
map<string,XSecStat> OptStatsFunc() { return optstats; }
map<string,CrossSection> OptXsFunc() { return optxs; }
/**
* Collect statistics for this event handler. To be used for
* histogram scaling.
*/
XSecStat histStats;
map<string,XSecStat> opthistStats;
/*
* The weight identifiers for the events
*/
vector<string> weightnames;
private:
/**
* The list of readers.
*/
ReaderVector theReaders;
/**
* The selector to choose readers according to their overestimated
* cross section.
*/
ReaderSelector theSelector;
/**
* The way weights are to be treated.
*/
WeightOpt theWeightOption;
/**
* If the weight option is set to unit weight, do not start
* compensating unless the weight is this much larger than unity.
*/
double theUnitTolerance;
/**
* Warn if the same process number is used in more than one
* FxFxReader.
*/
bool warnPNum;
/**
* How to normalize the weights
*/
unsigned int theNormWeight;
+ /**
+ * How to number the events
+ */
+ unsigned int UseLHEEvent;
+
public:
/** @cond EXCEPTIONCLASSES */
/**
* Exception class used if no readers were assigned.
*/
class FxFxInitError: public InitException {};
/**
* Exception class used if the same process number is used by more
* than ne reader.
*/
class FxFxPNumException: public InitException {};
/** @endcond */
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<FxFxEventHandler> initFxFxEventHandler;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
FxFxEventHandler & operator=(const FxFxEventHandler &) = delete;
};
}
// CLASSDOC OFF
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of FxFxEventHandler. */
template <>
struct BaseClassTrait<FxFxEventHandler,1> {
/** Typedef of the first base class of FxFxEventHandler. */
typedef EventHandler NthBase;
};
/** This template specialization informs ThePEG about the name of
* the FxFxEventHandler class and the shared object where it is defined. */
template <>
struct ClassTraits<FxFxEventHandler>
: public ClassTraitsBase<FxFxEventHandler> {
/** Return a platform-independent class name */
static string className() { return "Herwig::FxFxEventHandler"; }
/** Return the name of the shared library be loaded to get access to
* the FxFxEventHandler class and every other class it uses
* (except the base class). */
static string library() { return "HwFxFx.so"; }
};
/** @endcond */
}
#endif /* THEPEG_FxFxEventHandler_H */
diff --git a/MatrixElement/FxFx/FxFxFileReader.cc b/MatrixElement/FxFx/FxFxFileReader.cc
--- a/MatrixElement/FxFx/FxFxFileReader.cc
+++ b/MatrixElement/FxFx/FxFxFileReader.cc
@@ -1,929 +1,946 @@
// -*- C++ -*-
//
// FxFxFileReader.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2019 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the FxFxFileReader class.
//
#include "FxFxFileReader.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 <sstream>
#include <iostream>
using namespace ThePEG;
FxFxFileReader::
FxFxFileReader(const FxFxFileReader & x)
: FxFxReader(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),
theIncludeFxFxTags(x.theIncludeFxFxTags),
theIncludeCentral(x.theIncludeCentral),
theDecayer(x.theDecayer) {}
FxFxFileReader::~FxFxFileReader() {}
IBPtr FxFxFileReader::clone() const {
return new_ptr(*this);
}
IBPtr FxFxFileReader::fullclone() const {
return new_ptr(*this);
}
bool FxFxFileReader::preInitialize() const {
return true;
}
void FxFxFileReader::doinit() {
FxFxReader::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;
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()
<< "FxFxFileReader::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()
<< "FxFxFileReader::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 FxFxFileReader::doinit()."
<< " Please check the Les Houches file.\n"
<< Exception::runerror;
tcPDPtr p = getParticleData(t);
if( !p )
throw SetupException()
<< "FxFxFileReader::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()
<< "FxFxFileReader::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 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>()
<< "FxFxFileReader::doinit() Decayer must be set using the "
<< "FxFxFileReader:Decayer"
<< " must be set to allow the creation of new"
<< " decay modes."
<< Exception::runerror;
if(!dm) {
dm = generator()->preinitCreateDecayMode(tag);
if(!dm)
Throw<SetupException>()
<< "FxFxFileReader::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, "Active", "set", "Yes");
if(dm->CC()) {
generator()->preinitInterface(dm->CC(), "BranchingRatio", "set", br.str());
generator()->preinitInterface(dm->CC(), "Active", "set", "Yes");
}
++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 FxFxFileReader::doinit()"
<< Exception::runerror;
}
void FxFxFileReader::initialize(FxFxEventHandler & eh) {
FxFxReader::initialize(eh);
if ( LHFVersion.empty() )
Throw<FxFxFileError>()
<< "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> FxFxFileReader::optWeightNamesFunc() { return optionalWeightsNames; }
vector<string> FxFxFileReader::optWeightsNamesFunc() { return optionalWeightsNames; }
void FxFxFileReader::open() {
if ( filename().empty() )
throw FxFxFileError()
<< "No Les Houches file name. "
<< "Use 'set " << name() << ":FileName'."
<< Exception::runerror;
cfile.open(filename());
if ( !cfile )
throw FxFxFileError()
<< "The FxFxFileReader '" << 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
*/
string str_weightgroup = "<weightgroup";
string str_arrow = ">";
string str_newline = "\n";
erase_substr(weightinfo, str_weightgroup);
erase_substr(weightinfo, str_arrow);
erase_substr(weightinfo, str_newline);
}
/* if we are reading a new weightgroup, go on
* until we find the end of it
*/
if(readingInitWeights_sc && !cfile.find("</weightgroup") && !cfile.find("<weightgroup")) {
hs = cfile.getline();
//cout << "hs=" << hs << endl;
//cout << "weightinfo= " << weightinfo << endl;
//fix for potential new lines:
if(!cfile.find("<weight") and !cfile.find("</weightgroup")) {
weightinfo = weightinfo + hs;
//cout << "weightinfo fixed= " << weightinfo << endl;
continue;
}
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) { string str_arrow = ">"; erase_substr(sub, str_arrow); 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
erase_substr(scinfo,stopDEL);
erase_substr(scinfo,startDEL);
scinfo = StringUtils::stripws(scinfo);
//cout << "scinfo = " << scinfo << endl;
/* fill in the map
* indicating the information to be appended to each scale
* i.e. scinfo for each scalname
*/
scalemap[scalename] = scinfo.c_str();
string str_id = "id=";
string str_prime = "'";
erase_substr(scalename, str_id);
erase_substr(scalename, str_prime);
optionalWeightsNames.push_back(scalename);
}
}
if ( cfile.find("<header") ) {
// 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";
if (theIncludeCentral) optionalWeightsNames.push_back(central);
// cout << "reading init finished" << endl;
if ( !cfile ) {
heprup.NPRUP = -42;
LHFVersion = "";
return;
}
}
bool FxFxFileReader::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);
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
*/
if(theIncludeFxFxTags) optionalWeights[npstrings.c_str()] = -999;
string ecomstring = "ecom";
optionalWeights[ecomstring.c_str()] = heprup.EBMUP.first+heprup.EBMUP.second;
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(std::isnan(hepeup.PUP[i][0])||std::isnan(hepeup.PUP[i][1])||
std::isnan(hepeup.PUP[i][2])||std::isnan(hepeup.PUP[i][3])||
std::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;
}
}
+ LHEeventnum = -1; // set the LHEeventnum to -1, this will be the default if the tag <event num> is not found
+
+
// 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) { string str_arrow = ">" ; erase_substr(sub, str_arrow); weightName = sub; }
if(wi==2) weightValue = atof(sub.c_str());
++wi;
} while (iss);
// store the optional weights found in the temporary map
//cout << "weightName, weightValue= " << weightName << ", " << weightValue << endl;
optionalWeightsTemp[weightName] = weightValue;
}
/* reading of aMCFast weights
*/
if(readingaMCFast) {
std::stringstream amcfstringstream;
amcfstringstream << "aMCFast " << cfile.getline();
std::string amcfstrings = amcfstringstream.str();
string str_newline = "\n";
erase_substr(amcfstrings,str_newline);
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
erase_substr(mg5clusinfo,stopDEL);
erase_substr(mg5clusinfo,startDEL);
string str_arrow = ">";
erase_substr(mg5clusinfo,str_arrow);
string str_quotation = "\"";
erase_substr(mg5clusinfo,str_quotation);
string str_newline= "\n";
erase_substr(mg5clusinfo,str_newline);
optionalWeights[mg5clusinfo.c_str()] = -222; //for the mg5 scale info weights we give them a weight -222 for future identification
}
+
+ //the event num tag
+ if(cfile.find("<event num")) {
+ string hs = cfile.getline();
+ string startDEL = "<event num"; //starting delimiter
+ string stopDEL = "</event num>"; //end delimiter
+ unsigned firstLim = hs.find(startDEL);
+ string LHEeventnum_str = hs.substr(firstLim);
+ erase_substr(LHEeventnum_str,stopDEL);
+ erase_substr(LHEeventnum_str,startDEL);
+ string str_arrow = ">";
+ erase_substr(LHEeventnum_str,str_arrow);
+ LHEeventnum = std::stol(LHEeventnum_str, nullptr, 10);
+ }
//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
erase_substr(mg5scinfo,stopDEL);
erase_substr(mg5scinfo,startDEL);
string str_arrow = ">";
erase_substr(mg5scinfo,str_arrow);
string str_newline= "\n";
erase_substr(mg5scinfo,str_newline);
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){
string itfirst = it->first;
erase_substr(itfirst, "'");
erase_substr(itfirst, "\"");
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
string it2first = it2->first;
erase_substr(it2first, "'");
erase_substr(it2first, "\"");
//cout << "itfirst, it2first = " << itfirst << "\t" << it2first << endl;
if(itfirst==it2first) {
string info = it2->second;
string str_newline = "\n";
erase_substr(info, str_newline);
//set the optional weights
optionalWeights[info] = it->second;
}
}
}
/* additionally, we set the "central" scale
* this is actually the default event weight
*/
string central = "central";
if (theIncludeCentral) optionalWeights[central] = hepeup.XWGTUP;
if ( !cfile ) return false;
return true;
}
void FxFxFileReader::close() {
cfile.close();
}
void FxFxFileReader::persistentOutput(PersistentOStream & os) const {
os << neve << LHFVersion << outsideBlock << headerBlock << initComments
<< initAttributes << eventComments << eventAttributes << theFileName
<< theQNumbers << theIncludeFxFxTags << theIncludeCentral << theDecayer;
}
void FxFxFileReader::persistentInput(PersistentIStream & is, int) {
is >> neve >> LHFVersion >> outsideBlock >> headerBlock >> initComments
>> initAttributes >> eventComments >> eventAttributes >> theFileName
>> theQNumbers >> theIncludeFxFxTags >> theIncludeCentral >> theDecayer;
ieve = 0;
}
ClassDescription<FxFxFileReader>
FxFxFileReader::initFxFxFileReader;
// Definition of the static class description member.
void FxFxFileReader::Init() {
static ClassDocumentation<FxFxFileReader> documentation
("ThePEG::FxFxFileReader 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<FxFxFileReader,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.",
&FxFxFileReader::theFileName, "", false, false);
interfaceFileName.fileType();
interfaceFileName.rank(11);
static Switch<FxFxFileReader,bool> interfaceQNumbers
("QNumbers",
"Whether or not to read search for and read a QNUMBERS"
" block in the header of the file.",
&FxFxFileReader::theQNumbers, false, false, false);
static SwitchOption interfaceQNumbersYes
(interfaceQNumbers,
"Yes",
"Use QNUMBERS",
true);
static SwitchOption interfaceQNumbersNo
(interfaceQNumbers,
"No",
"Don't use QNUMBERS",
false);
static Switch<FxFxFileReader,bool> interfaceIncludeFxFxTags
("IncludeFxFxTags",
"Include FxFx tags",
&FxFxFileReader::theIncludeFxFxTags, true, true, false);
static SwitchOption interfaceIncludeFxFxTagsYes
(interfaceIncludeFxFxTags,
"Yes",
"Use the FxFx tags",
true);
static SwitchOption interfaceIncludeFxFxTagsNo
(interfaceIncludeFxFxTags,
"No",
"Don't use the FxFx tags",
false);
static Switch<FxFxFileReader,bool> interfaceIncludeCentral
("IncludeCentral",
"Include definition of central weight",
&FxFxFileReader::theIncludeCentral, false, true, false);
static SwitchOption interfaceIncludeCentralYes
(interfaceIncludeCentral,
"Yes",
"include definition of central weight",
true);
static SwitchOption interfaceIncludeCentralNo
(interfaceIncludeCentral,
"No",
"Don't include definition of central weight",
false);
static Reference<FxFxFileReader,Decayer> interfaceDecayer
("Decayer",
"Decayer to use for any decays read from the QNUMBERS Blocks",
&FxFxFileReader::theDecayer, false, false, true, true, false);
}
void FxFxFileReader::erase_substr(std::string& subject, const std::string& search) {
size_t pos = 0;
while((pos = subject.find(search, pos)) != std::string::npos) {
subject.erase( pos, search.length() );
}
}
diff --git a/MatrixElement/FxFx/FxFxReader.cc b/MatrixElement/FxFx/FxFxReader.cc
--- a/MatrixElement/FxFx/FxFxReader.cc
+++ b/MatrixElement/FxFx/FxFxReader.cc
@@ -1,1601 +1,1605 @@
// -*- C++ -*-
//
// FxFxReader.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2019 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the FxFxReader class.
//
#include "FxFxReader.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 "FxFxEventHandler.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;
FxFxReader::FxFxReader(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) {}
FxFxReader::FxFxReader(const FxFxReader & 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) {}
FxFxReader::~FxFxReader() {}
void FxFxReader::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 FxFxReader::preInitialize() const {
if ( HandlerBase::preInitialize() ) return true;
if ( doInitPDFs && ! ( inPDF.first && inPDF.second ) ) return true;
return false;
}
void FxFxReader::doinit() {
HandlerBase::doinit();
open();
close();
if ( !heprup.IDBMUP.first || !heprup.IDBMUP.second )
Throw<FxFxInitError>()
<< "No information about incoming particles were found in "
<< "FxFxReader '" << 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<FxFxInitError>()
<< "No information about the energy of incoming particles were found in "
<< "FxFxReader '" << name() << "'." << Exception::warning;
if ( doInitPDFs && ! ( inPDF.first && inPDF.second ) ) {
initPDFs();
if ( ! ( inPDF.first && inPDF.second ) ) Throw<InitException>()
<< "FxFxReader '" << name()
<< "' could not create PDFBase objects in pre-initialization."
<< Exception::warning;
}
else if ( !inPDF.first || !inPDF.second ) Throw<FxFxInitError>()
<< "No information about the PDFs of incoming particles were found in "
<< "FxFxReader '" << name() << "'." << Exception::warning;
}
void FxFxReader::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>()
<< "FxFxReader '" << 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>()
<< "FxFxReader '" << 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>()
<< "FxFxReader '" << name()
<< "' could not find information about the PDFs used."
<< Exception::warning;
}
void FxFxReader::initialize(FxFxEventHandler & eh) {
Energy2 Smax = ZERO;
double Y = 0.0;
if ( !theCuts ) {
theCuts = eh.cuts();
if ( !theCuts ) Throw<FxFxInitError>()
<< "No Cuts object was assigned to the FxFxReader '"
<< name() << "' nor was one\nassigned to the controlling "
<< "FxFxEventHandler '" << 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<FxFxInitError>()
<< "No PartonExtractor object was assigned to the FxFxReader '"
<< name() << "' nor was one\nassigned to the controlling "
<< "FxFxEventHandler '" << 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<FxFxInitError>()
<< "The FxFxReader '" << name() << "' uses the same Cuts object "
<< "as another FxFxReader which has not got the same energies of "
<< "the colliding particles. For the generation to work properly "
<< "different FxFxReader 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));
close();
if ( !heprup.IDWTUP && useWeightWarnings )
Throw<FxFxInitError>()
<< "No information about the weighting scheme was found. The events "
<< "produced by FxFxReader " << name()
<< " may not be sampled correctly." << Exception::warning;
if ( abs(heprup.IDWTUP) > 1 && !eh.weighted() && useWeightWarnings )
Throw<FxFxInitError>()
<< "FxFxReader " << 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 FxFxReader 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<FxFxInitError>()
<< "FxFxReader " << name() << " has the IDWTUP flag set to "
<< heprup.IDWTUP
<< ", which does not correspond\nto the weight option "
<< eh.weightOption() << " set in "
<< "the FxFxEventHandler " << 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 FxFxReader::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;
vector<double> sumlprup;
vector<double> sumsqlprup;
vector<long> nscanned;
for ( int i = 0; ( maxScan() < 0 || i < maxScan() ) && readEvent(); ++i ) {
if ( !checkPartonBin() ) Throw<FxFxInitError>()
<< "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);
sumlprup.push_back(0.);
sumsqlprup.push_back(0.);
nscanned.push_back(0);
} else {
id = idit - lprup.begin();
}
++neve;
++oldeve[id];
oldsum += hepeup.XWGTUP;
sumlprup[id] += hepeup.XWGTUP;
sumsqlprup[id] += sqr(hepeup.XWGTUP);
++nscanned[id];
if ( cacheFile() ) {
if ( eventWeight() == 0.0 ) {
++cuteve;
continue;
}
cacheEvent();
}
++neweve[id];
newmax[id] = max(newmax[id], abs(eventWeight()));
if ( eventWeight() < 0.0 ) negw = true;
} //end of scanning events
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<FxFxInitError>()
<< "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];
}
}
if ( abs(heprup.IDWTUP) != 1 ) {
// Try to fix things if abs(heprup.IDWTUP) != 1.
double sumxsec = 0.0;
if(abs(heprup.IDWTUP)==3) {
for ( int id = 0; id < heprup.NPRUP; ++id ) sumxsec += heprup.XSECUP[id];
}
else {
for ( int id = 0; id < heprup.NPRUP; ++id ) {
//set the cross section directly from the event weights read
heprup.XSECUP[id] = sumlprup[id]/nscanned[id];
heprup.XERRUP[id] = (sumsqlprup[id]/nscanned[id] - sqr(sumlprup[id]/nscanned[id])) / nscanned[id];
if(fabs(heprup.XERRUP[id]) < 1e-10) { heprup.XERRUP[id] = 0; }
if(fabs(newmax[id]) < 1e-10) { newmax[id] = 0; }
if(heprup.XERRUP[id] < 0.) {
if( heprup.XERRUP[id]/(sumsqlprup[id]/nscanned[id])>-1e-10)
heprup.XERRUP[id] = 0.;
else {
Throw<FxFxInitError>()
<< "Negative error when scanning events in LesHouschesReader '" << name()
<< Exception::warning;
heprup.XERRUP[id] = 0.;
}
}
heprup.XERRUP[id] = sqrt( heprup.XERRUP[id] );
heprup.XMAXUP[id] = newmax[id];
//cout << "heprup.XMAXUP[id] = " << heprup.XMAXUP[id] << endl;
sumxsec += heprup.XSECUP[id];
}
}
//cout << "sumxsec = " << sumxsec << endl;
//weightScale = picobarn*neve*sumxsec/oldsum;
weightScale = 1.0*picobarn; // temporary fix?
// cout << "weightscale = " << weightScale/picobarn << endl;
}
}
if ( cacheFile() ) closeCacheFile();
if ( negw ) heprup.IDWTUP = min(-abs(heprup.IDWTUP), -1);
return neve;
}
void FxFxReader::setWeightScale(long) {}
void FxFxReader::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 FxFxReader::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 FxFxReader::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 FxFxReader::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 FxFxReader::fillEvent() {
if ( !particleIndex.empty() ) return;
particleIndex.clear();
colourIndex.clear();
colourIndex(0, tColinePtr());
createParticles();
createBeams();
}
void FxFxReader::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(FxFxReopenWarning()
<< "Reopening FxFxReader '" << name()
<< "' after accessing " << stats.attempts()
<< " events out of "
<< NEvents() << Exception::warning);
else
throw FxFxReopenWarning()
<< "More events requested than available in FxFxReader "
<< name() << Exception::runerror;
}
if ( cacheFile() ) {
closeCacheFile();
openReadCacheFile();
if ( !uncacheEvent() ) Throw<FxFxReopenError>()
<< "Could not reopen FxFxReader '" << name()
<< "'." << Exception::runerror;
} else {
close();
open();
if ( !readEvent() ) Throw<FxFxReopenError>()
<< "Could not reopen FxFxReader '" << name()
<< "'." << Exception::runerror;
}
}
void FxFxReader::reset() {
particleIndex.clear();
colourIndex.clear();
if ( theLastXComb ) theLastXComb->clean();
theLastXComb = tXCombPtr();
}
bool FxFxReader::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
// FxFxReader 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 FxFxReader::getEvent() {
if ( cacheFile() ) {
if ( !uncacheEvent() ) reopen();
} else {
if ( !readEvent() ) reopen();
}
++position;
double max = maxWeights[hepeup.IDPRUP]*maxFactor;
// cout << "maxFactor = " << maxFactor << " maxWeights[hepeup.IDPRUP] = " << maxWeights[hepeup.IDPRUP] << endl;
// normalize all the weights to the max weight
for(map<string,double>::iterator it=optionalWeights.begin();
it!=optionalWeights.end();++it) {
if(it->first!="ecom" && it->second!=-999 && it->second!=-111 && it->second!=-222 && it->second!=-333) {
it->second = (max != 0.0) ? it->second/max : 0.0;
}
}
//cout << "hepeup.XWGTUP = " << hepeup.XWGTUP << endl;
//cout << "max = " << max << " " << eventWeight() << endl;
return max != 0.0? eventWeight()/max: 0.0;
}
void FxFxReader::skip(long n) {
HoldFlag<> skipflag(skipping);
while ( n-- ) getEvent();
}
double FxFxReader::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 FxFxReader::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 FxFxReader::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<FxFxInconsistencyError>()
<< "Could not find appropriate PartonBin objects for event produced by "
<< "FxFxReader '" << 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<FxFxInconsistencyError>()
<< "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<FxFxInconsistencyError>()
<< "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 FxFxReader::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);
if(theMomentumTreatment == 1) mom.rescaleEnergy();
else if(theMomentumTreatment == 2) mom.rescaleMass();
PDPtr pd = getParticleData(hepeup.IDUP[i]);
if (!pd) {
Throw<FxFxInitError>()
<< "FxFxReader '" << 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<FxFxInconsistencyError>()
<< "FxFxReader " << 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<FxFxInconsistencyError>()
<< "To many incoming beam particles in the FxFxReader '"
<< 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<FxFxInconsistencyError>()
<< "To many incoming particles to hard subprocess in the "
<< "FxFxReader '" << 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<FxFxInconsistencyError>()
<< "Unknown status code (" << hepeup.ISTUP[i]
<< ") in the FxFxReader '" << 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<FxFxInconsistencyError>()
<< "Polarization must be between -1 and 1, not "
<< spinup << " as found in the "
<< "FxFx 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) {
+ if(iy==ix) continue;
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) {
+ 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(hepeup.ISTUP[particleIndex(external[iy])-1]<0) {
- if(external[iy]->colourInfo()->colourLines().empty()) continue;
+ if(external[iy]->colourInfo()->antiColourLines().empty()) continue;
anti2 = external[iy]->colourInfo()->antiColourLines();
}
else if(hepeup.ISTUP[particleIndex(external[iy])-1]>0) {
- if(external[iy]->colourInfo()->antiColourLines().empty()) continue;
+ if(external[iy]->colourInfo()->colourLines().empty()) continue;
anti2 = external[iy]->colourInfo()->colourLines();
}
for(unsigned int ic2=0;ic2<anti2.size();++ic2) {
- if(col[ic1]==anti2[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<FxFxInconsistencyError>()
<< "FxFxReader '" << 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<FxFxInconsistencyError>()
<< "FxFxReader '" << 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<FxFxInconsistencyError>()
<< "FxFxReader '" << 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<FxFxInconsistencyError>()
<< "FxFxReader '" << name() << "' found inconsistent colour "
<< "flow in Les Houches common block structure for decay of\n"
<< *external[0] << "\n"
<< hepeup << Exception::runerror;
}
}
}
}
void FxFxReader::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 FxFxReader::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 FxFxReader::openReadCacheFile() {
if ( cacheFile() ) closeCacheFile();
cacheFile().open(cacheFileName(), "r");
position = 0;
}
void FxFxReader::openWriteCacheFile() {
if ( cacheFile() ) closeCacheFile();
cacheFile().open(cacheFileName(), "w");
}
void FxFxReader::closeCacheFile() {
cacheFile().close();
}
void FxFxReader::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);
+ pos = mwrite(pos, LHEeventnum);
for(size_t ff = 0; ff < optionalWeightsNames.size(); ff++) {
pos = mwrite(pos, optionalWeightsNames[ff]);
}
/* for(int f = 0; f < optionalWeightsNames.size(); f++) {
cout << "optionalWeightsNames = " << optionalWeightsNames[f] << endl;
}*/
pos = mwrite(pos, optionalnpLO);
pos = mwrite(pos, optionalnpNLO);
pos = mwrite(pos, preweight);
cacheFile().write(&buff[0], buff.size(), 1);
}
bool FxFxReader::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);
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(size_t ff = 0; ff < optionalWeightsNames.size(); ff++) {
pos = mread(pos, optionalWeightsNames[ff]);
}
pos = mread(pos, optionalnpLO);
pos = mread(pos, optionalnpNLO);
pos = mread(pos, preweight);
+ pos = mread(pos, LHEeventnum);
// 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 FxFxReader::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
+ << theLastXComb << theMaxMultCKKW << theMinMultCKKW << lastweight << optionalWeights << optionalnpLO << optionalnpNLO << LHEeventnum
<< maxFactor << ounit(weightScale, picobarn) << xSecWeights << maxWeights
<< theMomentumTreatment << useWeightWarnings << theReOpenAllowed
<< theIncludeSpin;
}
void FxFxReader::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
+ >> theLastXComb >> theMaxMultCKKW >> theMinMultCKKW >> lastweight >> optionalWeights >> optionalnpLO >> optionalnpNLO >> LHEeventnum
>> maxFactor >> iunit(weightScale, picobarn) >> xSecWeights >> maxWeights
>> theMomentumTreatment >> useWeightWarnings >> theReOpenAllowed
>> theIncludeSpin;
}
AbstractClassDescription<FxFxReader>
FxFxReader::initFxFxReader;
// Definition of the static class description member.
void FxFxReader::setBeamA(long id) { heprup.IDBMUP.first = id; }
long FxFxReader::getBeamA() const { return heprup.IDBMUP.first; }
void FxFxReader::setBeamB(long id) { heprup.IDBMUP.second = id; }
long FxFxReader::getBeamB() const { return heprup.IDBMUP.second; }
void FxFxReader::setEBeamA(Energy e) { heprup.EBMUP.first = e/GeV; }
Energy FxFxReader::getEBeamA() const { return heprup.EBMUP.first*GeV; }
void FxFxReader::setEBeamB(Energy e) { heprup.EBMUP.second = e/GeV; }
Energy FxFxReader::getEBeamB() const { return heprup.EBMUP.second*GeV; }
void FxFxReader::setPDFA(PDFPtr pdf) { inPDF.first = pdf; }
PDFPtr FxFxReader::getPDFA() const { return inPDF.first; }
void FxFxReader::setPDFB(PDFPtr pdf) { inPDF.second = pdf; }
PDFPtr FxFxReader::getPDFB() const { return inPDF.second; }
void FxFxReader::Init() {
static ClassDocumentation<FxFxReader> documentation
("ThePEG::FxFxReader is an abstract base class to be used "
"for objects which reads event files or streams from matrix element "
"generators.");
static Parameter<FxFxReader,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,
&FxFxReader::setBeamA,
&FxFxReader::getBeamA, 0, 0, 0);
static Parameter<FxFxReader,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,
&FxFxReader::setBeamB,
&FxFxReader::getBeamB, 0, 0, 0);
static Parameter<FxFxReader,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,
&FxFxReader::setEBeamA, &FxFxReader::getEBeamA, 0, 0, 0);
static Parameter<FxFxReader,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,
&FxFxReader::setEBeamB, &FxFxReader::getEBeamB, 0, 0, 0);
static Reference<FxFxReader,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,
&FxFxReader::setPDFA, &FxFxReader::getPDFA, 0);
static Reference<FxFxReader,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,
&FxFxReader::setPDFB, &FxFxReader::getPDFB, 0);
static Parameter<FxFxReader,long> interfaceMaxScan
("MaxScan",
"The maximum number of events to scan to obtain information about "
"processes and cross section in the intialization.",
&FxFxReader::theMaxScan, -1, 0, 0,
true, false, false);
static Parameter<FxFxReader,string> interfaceCacheFileName
("CacheFileName",
"Name of file used to cache the events form the reader in a fast-readable "
"form. If empty, no cache file will be generated.",
&FxFxReader::theCacheFileName, "",
true, false);
interfaceCacheFileName.fileType();
static Switch<FxFxReader,bool> interfaceCutEarly
("CutEarly",
"Determines whether to apply cuts to events before converting to "
"ThePEG format.",
&FxFxReader::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<FxFxReader,PartonExtractor> interfacePartonExtractor
("PartonExtractor",
"The PartonExtractor object used to construct remnants. If no object is "
"provided the FxFxEventHandler object must provide one instead.",
&FxFxReader::thePartonExtractor, true, false, true, true, false);
static Reference<FxFxReader,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 FxFxEventHandler object must "
"provide one instead.",
&FxFxReader::theCuts, true, false, true, true, false);
static RefVector<FxFxReader,ReweightBase> interfaceReweights
("Reweights",
"A list of ThePEG::ReweightBase objects to modify this the weight of "
"this reader.",
&FxFxReader::reweights, 0, false, false, true, false);
static RefVector<FxFxReader,ReweightBase> interfacePreweights
("Preweights",
"A list of ThePEG::ReweightBase objects to bias the phase space for this "
"reader without influencing the actual cross section.",
&FxFxReader::preweights, 0, false, false, true, false);
static Switch<FxFxReader,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?",
&FxFxReader::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<FxFxReader,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.",
&FxFxReader::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<FxFxReader,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.",
&FxFxReader::theMaxMultCKKW, 0, 0, 0,
true, false, Interface::lowerlim);
static Parameter<FxFxReader,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.",
&FxFxReader::theMinMultCKKW, 0, 0, 0,
true, false, Interface::lowerlim);
static Switch<FxFxReader,unsigned int> interfaceMomentumTreatment
("MomentumTreatment",
"Treatment of the momenta supplied by the interface",
&FxFxReader::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<FxFxReader,bool> interfaceWeightWarnings
("WeightWarnings",
"Determines if warnings about possible weight incompatibilities should "
"be issued when this reader is initialized.",
&FxFxReader::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<FxFxReader,bool> interfaceAllowedTopReOpen
("AllowedToReOpen",
"Can the file be reopened if more events are requested than the file contains?",
&FxFxReader::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<FxFxReader,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",
&FxFxReader::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;
}
}
diff --git a/MatrixElement/FxFx/FxFxReader.h b/MatrixElement/FxFx/FxFxReader.h
--- a/MatrixElement/FxFx/FxFxReader.h
+++ b/MatrixElement/FxFx/FxFxReader.h
@@ -1,1006 +1,1016 @@
// -*- C++ -*-
//
// FxFxReader.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2019 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef THEPEG_FxFxReader_H
#define THEPEG_FxFxReader_H
// This is the declaration of the FxFxReader class.
#include "FxFx.h"
#include "ThePEG/Handlers/HandlerBase.h"
#include "ThePEG/Utilities/ObjectIndexer.h"
#include "ThePEG/Utilities/Exception.h"
#include "ThePEG/Utilities/XSecStat.h"
#include "ThePEG/PDF/PartonBinInstance.h"
#include "ThePEG/PDF/PartonBin.fh"
#include "ThePEG/MatrixElement/ReweightBase.h"
#include "FxFxEventHandler.fh"
#include "FxFxReader.fh"
#include "ThePEG/Utilities/CFile.h"
#include <cstdio>
#include <cstring>
namespace ThePEG {
/**
* FxFxReader is an abstract base class to be used for objects
* which reads event files or streams from matrix element
* generators. Derived classes must at least implement the open() and
* doReadEvent() methods to read in information about the whole run into
* the HEPRUP variable and next event into the HEPEUP variable
* respectively. Also the close() function to close the file or stream
* read must be implemented. Although these functions are named as if
* we are reading from event files, they could just as well implement
* the actual generation of events.
*
* After filling the HEPRUP and HEPEUP variables, which are protected
* and easily accesible from the sub-class, this base class will then
* be responsible for transforming this data to the ThePEG Event
* record in the getEvent() method. <code>FxFxReader</code>s can
* only be used inside FxFxEventHandler objects.
*
* In the initialization the virtual open() and scan() functions are
* called. Here the derived class must provide the information about
* the processes in the variables corresponding to the HEPRUP common
* block. Note that the IDWTUP is required to be +/- 1, and sub
* classes are required to change the information accordingly to
* ensure the correct corss section sampling. Note also that the
* controlling FxFxEventHandler may choose to generate weighted
* events even if IDWTUP is 1.
*
* Note that the information given per process in e.g. the XSECUP and
* XMAXUP vectors is not used by the FxFxEventHandler and by
* default the FxFxReader is not assumed to be able to actively
* choose between the sub-processes. Instead, the
* FxFxEventHandler can handle several FxFxReader objects
* and choose between them. However, a sub-class of FxFxReader
* may set the flag isActive, in which case it is assumed to be able
* to select between its sub-processes itself.
*
* The FxFxReader may be assigned a number ReweightBase objects
* which either completely reweights the events produced (in the
* reweights vector), or only biases the selection without influencing
* the cross section (in the preweights vector). Note that it is the
* responsibility of a sub-class to call the reweight() function and
* multiply the weight according to its return value (typically done
* in the readEvent() function).
*
* @see \ref FxFxReaderInterfaces "The interfaces"
* defined for FxFxReader.
* @see Event
* @see FxFxEventHandler
*/
class FxFxReader: public HandlerBase, public LastXCombInfo<> {
/**
* FxFxEventHandler should have access to our private parts.
*/
friend class FxFxEventHandler;
/**
* Map for accumulating statistics of cross sections per process
* number.
*/
typedef map<int,XSecStat> StatMap;
/**
* Map of XComb objects describing the incoming partons indexed by
* the corresponding PartonBin pair.
*/
typedef map<tcPBPair,XCombPtr> XCombMap;
/**
* A vector of pointers to ReweightBase objects.
*/
typedef vector<ReweightPtr> ReweightVector;
public:
/** @name Standard constructors and destructors. */
//@{
/**
* Default constructor. If the optional argument is true, the reader
* is assumed to be able to produce events on demand for a given
* process.
*/
FxFxReader(bool active = false);
/**
* Copy-constructor.
*/
FxFxReader(const FxFxReader &);
/**
* Destructor.
*/
virtual ~FxFxReader();
//@}
public:
/** @name Main virtual fuctions to be overridden in
* sub-classes. They are named as if we are reading from event
* files, but could equally well implement the actual generation of
* events. */
//@{
/**
* Open a file or stream with events and read in the run information
* into the heprup variable.
*/
virtual void open() = 0;
/**
* Read the next event from the file or stream into the
* corresponding protected variables. Return false if there is no
* more events.
*/
virtual bool doReadEvent() = 0;
/**
* Close the file or stream from which events have been read.
*/
virtual void close() = 0;
/**
* return the weight names
*/
// virtual vector<string> optWeightsNamesFunc();
virtual vector<string> optWeightsNamesFunc() = 0;
//virtual vector<string*> optWeightNamesFunc() = 0;
vector<string> optionalWeightsNames;
/**
* The ID (e.g. 100x, 2001) for the weight
*/
// vector<string> optionalWeightsNames;
//@}
/** @name Other important function which may be overridden in
* sub-classes which wants to bypass the basic HEPRUP or HEPEUP
* variables or otherwise facilitate the conversion to ThePEG
* objects. */
//@{
/**
* Initialize. This function is called by the FxFxEventHandler
* to which this object is assigned.
*/
virtual void initialize(FxFxEventHandler & eh);
/**
* Calls readEvent() or uncacheEvent() to read information into the
* FxFx common block variables. This function is called by the
* FxFxEventHandler if this reader has been selectod to
* produce an event.
*
* @return the weight asociated with this event. If negative weights
* are allowed it should be between -1 and 1, otherwise between 0
* and 1. If outside these limits the previously estimated maximum
* is violated. Note that the estimated maximum then should be
* updated from the outside.
*/
virtual double getEvent();
/**
* Calls doReadEvent() and performs pre-defined reweightings. A
* sub-class overrides this function it must make sure that the
* corresponding reweightings are done.
*/
virtual bool readEvent();
/**
* Skip \a n events. Used by FxFxEventHandler to make sure
* that a file is scanned an even number of times in case the events
* are not ramdomly distributed in the file.
*/
virtual void skip(long n);
/**
* Get an XComb object. Converts the information in the Les Houches
* common block variables to an XComb object describing the sub
* process. This is the way information is conveyed from the reader
* to the controlling FxFxEventHandler.
*/
tXCombPtr getXComb();
/**
* Get a SubProcess object corresponding to the information in the
* Les Houches common block variables.
*/
tSubProPtr getSubProcess();
/**
* Scan the file or stream to obtain information about cross section
* weights and particles etc. This function should fill the
* variables corresponding to the /HEPRUP/ common block. The
* function returns the number of events scanned.
*/
virtual long scan();
/**
* Take the information corresponding to the HEPRUP common block and
* initialize the statistics for this reader.
*/
virtual void initStat();
/**
* Reweights the current event using the reweights and preweights
* vectors. It is the responsibility of the sub-class to call this
* function after the HEPEUP information has been retrieved.
*/
double reweight();
/**
* Converts the information in the Les Houches common block
* variables into a Particle objects.
*/
virtual void fillEvent();
/**
* Removes the particles created in the last generated event,
* preparing to produce a new one.
*/
void reset();
/**
* Possibility for subclasses to recover from non-conformant
* settings of XMAXUP when an event file has been scanned with \a
* neve events. Should set weightScale so that the average XMAXUP
* times weightScale gives the cross section for a process. (This is
* needed for MadEvent).
*/
virtual void setWeightScale(long neve);
//@}
/** @name Access information about the current event. */
//@{
/**
* Return the size of this event in bytes. To be used for the cache
* file. \a npart is the number of particles. If \a npart is 0, the
* number is taken from NUP.
*/
static size_t eventSize(int N) {
return (N + 1)*sizeof(int) + // IDPRUP, ISTUP
(7*N + 4)*sizeof(double) + // XWGTUP, SCALUP, AQEDUP, AQCDUP, PUP,
// VTIMUP, SPINUP
N*sizeof(long) + // IDUP
2*N*sizeof(pair<int,int>) + // MOTHUP, ICOLUP
sizeof(pair<double,double>) + // XPDWUP.
2*sizeof(double); // lastweight and preweight
}
/**
* The current event weight given by XWGTUP times possible
* reweighting. Note that this is not necessarily the same as what
* is returned by getEvent(), which is scaled with the maximum
* weight.
*/
double eventWeight() const { return hepeup.XWGTUP*lastweight; }
/**
* Return the optional named weights associated to the current event.
*/
const map<string,double>& optionalEventWeights() const { return optionalWeights; }
+ /**
+ * Return the Les Houches event number associated with the current event
+ */
+ const long& LHEEventNum() const { return LHEeventnum; }
+
/**
* Return the optional npLO and npNLO
*/
const int& optionalEventnpLO() const { return optionalnpLO; }
const int& optionalEventnpNLO() const { return optionalnpNLO; }
/**
* The pair of PartonBinInstance objects describing the current
* incoming partons in the event.
*/
const PBIPair & partonBinInstances() const { return thePartonBinInstances; }
/**
* Return the instances of the beam particles for the current event.
*/
const PPair & beams() const { return theBeams; }
/**
* Return the instances of the incoming particles to the sub process
* for the current event.
*/
const PPair & incoming() const { return theIncoming; }
/**
* Return the instances of the outgoing particles from the sub process
* for the current event.
*/
const PVector & outgoing() const { return theOutgoing; }
/**
* Return the instances of the intermediate particles in the sub
* process for the current event.
*/
const PVector & intermediates() const { return theIntermediates; }
/**
* 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.
*/
int maxMultCKKW() const { return theMaxMultCKKW; }
/**
* 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.
*/
int minMultCKKW() const { return theMinMultCKKW; } //@}
/** @name Other inlined access functions. */
//@{
/**
* The number of events found in this reader. If less than zero the
* number of events are unlimited.
*/
long NEvents() const { return theNEvents; }
/**
* The number of events produced so far. Is reset to zero if an
* event file is reopened.
*/
long currentPosition() const { return position; }
/**
* The maximum number of events to scan to collect information about
* processes and cross sections. If less than 0, all events will be
* scanned.
*/
long maxScan() const { return theMaxScan; }
/**
* Return true if this reader is active.
*/
bool active() const { return isActive; }
/**
* True if negative weights may be produced.
*/
bool negativeWeights() const { return heprup.IDWTUP < 0; }
/**
* The collected cross section statistics for this reader.
*/
const XSecStat & xSecStats() const { return stats; }
/**
* Collected statistics about the individual processes.
*/
const StatMap & processStats() const { return statmap; }
/**
* Select the current event. It will later be rejected with a
* probability given by \a weight.
*/
void select(double weight) {
stats.select(weight);
statmap[hepeup.IDPRUP].select(weight);
}
/**
* Accept the current event assuming it was previously selcted.
*/
void accept() {
stats.accept();
statmap[hepeup.IDPRUP].accept();
}
/**
* Reject the current event assuming it was previously accepted.
*/
void reject(double w) {
stats.reject(w);
statmap[hepeup.IDPRUP].reject(w);
}
/**
* Increase the overestimated cross section for this reader.
*/
virtual void increaseMaxXSec(CrossSection maxxsec);
/**
* The PartonExtractor object used to construct remnants.
*/
tPExtrPtr partonExtractor() const { return thePartonExtractor; }
/**
* Return a possibly null pointer to a CascadeHandler to be used for
* CKKW-reweighting.
*/
tCascHdlPtr CKKWHandler() const { return theCKKW; }
/**
* The pairs of PartonBin objects describing the partons which can
* be extracted by the PartonExtractor object.
*/
const PartonPairVec & partonBins() const { return thePartonBins; }
/**
* The map of XComb objects indexed by the corresponding PartonBin
* pair.
*/
const XCombMap & xCombs() const { return theXCombs; }
/**
* The Cuts object to be used for this reader.
*/
const Cuts & cuts() const { return *theCuts; }
//@}
protected:
/** @name Functions for manipulating cache files. */
//@{
/**
* Name of file used to cache the events form the reader in a
* fast-readable form. If empty, no cache file will be generated.
*/
string cacheFileName() const { return theCacheFileName; }
/**
* Determines whether to apply cuts to events converting them to
* ThePEG format.
*/
bool cutEarly() const { return doCutEarly; }
/**
* File stream for the cache.
*/
CFile cacheFile() const { return theCacheFile;}
/**
* Open the cache file for reading.
*/
void openReadCacheFile();
/**
* Open the cache file for writing.
*/
void openWriteCacheFile();
/**
* Close the cache file;
*/
void closeCacheFile();
/**
* Write the current event to the cache file.
*/
void cacheEvent() const;
/**
* Read an event from the cache file. Return false if something went wrong.
*/
bool uncacheEvent();
/**
* Reopen a reader. If we have reached the end of an event file,
* reopen it and issue a warning if we have used up a large fraction
* of it.
*/
void reopen();
/**
* Helper function to write a variable to a memory location
*/
template <typename T>
static char * mwrite(char * pos, const T & t, size_t n = 1) {
std::memcpy(pos, &t, n*sizeof(T));
return pos + n*sizeof(T);
}
/**
* Helper function to read a variable from a memory location
*/
template <typename T>
static const char * mread(const char * pos, T & t, size_t n = 1) {
std::memcpy(&t, pos, n*sizeof(T));
return pos + n*sizeof(T);
}
//@}
/** @name Auxilliary virtual methods which may be verridden by sub-classes. */
//@{
/**
* Check the existence of a pair of PartonBin objects corresponding
* to the current event.
*
* @return false if no pair of suitable PartonBin objects was found.
*/
virtual bool checkPartonBin();
/**
* Create instances of all particles in the event and store them
* in particleIndex.
*/
virtual void createParticles();
/**
* Using the already created particles create a pair of
* PartonBinInstance objects corresponding to the incoming
* partons. Return the corresponding PartonBin objects.
*/
virtual tcPBPair createPartonBinInstances();
/**
* Create instances of the incoming beams in the event and store
* them in particleIndex. If no beam particles are included in the
* event they are created from the run info.
*/
virtual void createBeams();
/**
* Go through the mother indices and connect up the Particles.
*/
virtual void connectMothers();
//@}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* Standard Init function used to initialize the interfaces.
*/
static void Init();
protected:
/** @name Set functions for some variables not in the Les Houches accord. */
//@{
/**
* The number of events in this reader. If less than zero the number
* of events is unlimited.
*/
void NEvents(long x) { theNEvents = x; }
/**
* The map of XComb objects indexed by the corresponding PartonBin
* pair.
*/
XCombMap & xCombs() { return theXCombs; }
//@}
/** @name Standard (and non-standard) Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
/**
* Finalize this object. Called in the run phase just after a
* run has ended. Used eg. to write out statistics.
*/
virtual void dofinish() {
close();
HandlerBase::dofinish();
}
/**
* Return true if this object needs to be initialized before all
* other objects because it needs to extract PDFs from the event file.
*/
virtual bool preInitialize() const;
/**
* Called from doinit() to extract PDFs from the event file and add
* the corresponding objects to the current EventGenerator.
*/
virtual void initPDFs();
//@}
protected:
/**
* The HEPRUP common block.
*/
HEPRUP heprup;
/**
* The HEPEUP common block.
*/
HEPEUP hepeup;
/**
* The ParticleData objects corresponding to the incoming particles.
*/
tcPDPair inData;
/**
* The PDFBase objects which has been used for the beam particle
* when generating the events being read. Specified in the interface
* or derived from PDFGUP and PDFSUP.
*/
pair<PDFPtr,PDFPtr> inPDF;
/**
* The PDFBase object to be used in the subsequent generation.
*/
pair<cPDFPtr,cPDFPtr> outPDF;
/**
* The PartonExtractor object used to construct remnants.
*/
PExtrPtr thePartonExtractor;
/**
* A pointer to a CascadeHandler to be used for CKKW-reweighting.
*/
tCascHdlPtr theCKKW;
/**
* The pairs of PartonBin objects describing the partons which can
* be extracted by the PartonExtractor object.
*/
PartonPairVec thePartonBins;
/**
* The map of XComb objects indexed by the corresponding PartonBin
* pair.
*/
XCombMap theXCombs;
/**
* The Cuts object to be used for this reader.
*/
CutsPtr theCuts;
/**
* The number of events in this reader. If less than zero the number
* of events is unlimited.
*/
long theNEvents;
/**
* The number of events produced by this reader so far. Is reset
* every time an event file is reopened.
*/
long position;
/**
* The number of times this reader has been reopened.
*/
int reopened;
/**
* The maximum number of events to scan to collect information about
* processes and cross sections. If less than 0, all events will be
* scanned.
*/
long theMaxScan;
/**
* Flag to tell whether we are in the process of scanning.
*/
bool scanning;
/**
* True if this is an active reader.
*/
bool isActive;
/**
* Name of file used to cache the events form the reader in a
* fast-readable form. If empty, no cache file will be generated.
*/
string theCacheFileName;
/**
* Determines whether to apply cuts to events before converting them
* to ThePEG format.
*/
bool doCutEarly;
/**
* Collect statistics for this reader.
*/
XSecStat stats;
/**
* Collect statistics for each individual process.
*/
StatMap statmap;
/**
* The pair of PartonBinInstance objects describing the current
* incoming partons in the event.
*/
PBIPair thePartonBinInstances;
/**
* Association between ColourLines and colour indices in the current
* translation.
*/
ObjectIndexer<long,ColourLine> colourIndex;
/**
* Association between Particles and indices in the current
* translation.
*/
ObjectIndexer<long,Particle> particleIndex;
/**
* The instances of the beam particles for the current event.
*/
PPair theBeams;
/**
* The instances of the incoming particles to the sub process for
* the current event.
*/
PPair theIncoming;
/**
* The instances of the outgoing particles from the sub process for
* the current event.
*/
PVector theOutgoing;
/**
* The instances of the intermediate particles in the sub process for
* the current event.
*/
PVector theIntermediates;
/**
* File stream for the cache.
*/
CFile theCacheFile;
/**
* The reweight objects modifying the weights of this reader.
*/
ReweightVector reweights;
/**
* The preweight objects modifying the weights of this reader.
*/
ReweightVector preweights;
/**
* The factor with which this reader was last pre-weighted.
*/
double preweight;
/**
* Should the event be reweighted by PDFs used by the PartonExtractor?
*/
bool reweightPDF;
/**
* Should PDFBase objects be constructed from the information in the
* event file in the initialization?
*/
bool doInitPDFs;
/**
* 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.
*/
int theMaxMultCKKW;
/**
* 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.
*/
int theMinMultCKKW;
/**
* The weight multiplying the last read event due to PDF
* reweighting, CKKW reweighting or assigned reweight and preweight
* objects.
*/
double lastweight;
/**
* The optional weights associated to the last read events.
*/
map<string,double> optionalWeights;
+
+ /**
+ * The event number
+ */
+ long LHEeventnum;
/**
* If the maximum cross section of this reader has been increased
* with increaseMaxXSec(), this is the total factor with which it
* has been increased.
*/
double maxFactor;
/**
* npLO for FxFx merging
*/
int optionalnpLO;
/**
* npNLO for FxFx merging
*/
int optionalnpNLO;
/**
* The (reweighted) XWGTUP value should be scaled with this cross
* section when compared to the overestimated cross section.
*/
CrossSection weightScale;
/**
* Individual scales for different sub-processes if reweighted.
*/
vector<double> xSecWeights;
/**
* Individual maximum weights for individual (possibly reweighted)
* processes.
*/
map<int,double> maxWeights;
/**
* Is set to true when getEvent() is called from skip(int).
*/
bool skipping;
/**
* Option for the treatment of the momenta supplied
*/
unsigned int theMomentumTreatment;
/**
* Set to true if warnings about possible weight incompatibilities
* should be issued.
*/
bool useWeightWarnings;
/**
* Option to allow reopening of the file
*/
bool theReOpenAllowed;
/**
* Use the spin information
*/
bool theIncludeSpin;
private:
/** Access function for the interface. */
void setBeamA(long id);
/** Access function for the interface. */
long getBeamA() const;
/** Access function for the interface. */
void setBeamB(long id);
/** Access function for the interface. */
long getBeamB() const;
/** Access function for the interface. */
void setEBeamA(Energy e);
/** Access function for the interface. */
Energy getEBeamA() const;
/** Access function for the interface. */
void setEBeamB(Energy e);
/** Access function for the interface. */
Energy getEBeamB() const;
/** Access function for the interface. */
void setPDFA(PDFPtr);
/** Access function for the interface. */
PDFPtr getPDFA() const;
/** Access function for the interface. */
void setPDFB(PDFPtr);
/** Access function for the interface. */
PDFPtr getPDFB() const;
private:
/**
* Describe an abstract base class with persistent data.
*/
static AbstractClassDescription<FxFxReader> initFxFxReader;
/**
* Private and non-existent assignment operator.
*/
FxFxReader & operator=(const FxFxReader &) = delete;
public:
/** @cond EXCEPTIONCLASSES */
/** Exception class used by FxFxReader in case inconsistencies
* are encountered. */
class FxFxInconsistencyError: public Exception {};
/** Exception class used by FxFxReader in case more events
than available are requested. */
class FxFxReopenWarning: public Exception {};
/** Exception class used by FxFxReader in case reopening an
event file fails. */
class FxFxReopenError: public Exception {};
/** Exception class used by FxFxReader in case there is
information missing in the initialization phase. */
class FxFxInitError: public InitException {};
/** @endcond */
};
/// Stream output for HEPEUP
ostream & operator<<(ostream & os, const HEPEUP & h);
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/**
* This template specialization informs ThePEG about the
* base class of FxFxReader.
*/
template <>
struct BaseClassTrait<FxFxReader,1>: public ClassTraitsType {
/** Typedef of the base class of FxFxReader. */
typedef HandlerBase NthBase;
};
/**
* This template specialization informs ThePEG about the name of the
* FxFxReader class and the shared object where it is
* defined.
*/
template <>
struct ClassTraits<FxFxReader>
: public ClassTraitsBase<FxFxReader> {
/**
* Return the class name.
*/
static string className() { return "Herwig::FxFxReader"; }
/**
* Return the name of the shared library to be loaded to get access
* to the FxFxReader class and every other class it uses
* (except the base class).
*/
static string library() { return "HwFxFx.so"; }
};
/** @endcond */
}
#endif /* THEPEG_FxFxReader_H */
diff --git a/MatrixElement/Matchbox/Base/SubtractedME.h b/MatrixElement/Matchbox/Base/SubtractedME.h
--- a/MatrixElement/Matchbox/Base/SubtractedME.h
+++ b/MatrixElement/Matchbox/Base/SubtractedME.h
@@ -1,528 +1,536 @@
// -*- C++ -*-
//
// SubtractedME.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_SubtractedME_H
#define HERWIG_SubtractedME_H
//
// This is the declaration of the SubtractedME class.
//
#include "ThePEG/MatrixElement/MEGroup.h"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "Herwig/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
#include "Herwig/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h"
#include "Herwig/MatrixElement/Matchbox/MatchboxFactory.fh"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief SubtractedME represents a subtracted real emission matrix element.
*
* @see \ref SubtractedMEInterfaces "The interfaces"
* defined for SubtractedME.
*/
class SubtractedME:
public MEGroup,
public LastMatchboxXCombInfo {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
SubtractedME();
/**
* The destructor.
*/
virtual ~SubtractedME();
//@}
public:
/**
* Return the factory which produced this matrix element
*/
Ptr<MatchboxFactory>::tcptr factory() const;
/** @name Phasespace and subprocess information */
//@{
/**
* For the given event generation setup return an xcomb object
* appropriate to this matrix element.
*/
virtual StdXCombPtr makeXComb(Energy newMaxEnergy, const cPDPair & inc,
tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
const PBPair & newPartonBins, tCutsPtr newCuts,
const DiagramVector & newDiagrams, bool mir,
const PartonPairVec& allPBins,
tStdXCombPtr newHead = tStdXCombPtr(),
tMEPtr newME = tMEPtr());
/**
* Set the XComb object to be used in the next call to
* generateKinematics() and dSigHatDR().
*/
virtual void setXComb(tStdXCombPtr);
/**
* Return true, if the same additional random numbers
* should be presented to any of the dependent
* matrix elements.
*/
virtual bool uniformAdditional() const { return true; }
/**
* Return true, if the XComb steering this matrix element
* should keep track of the random numbers used to generate
* the last phase space point
*/
virtual bool keepRandomNumbers() const { return true; }
/**
* Given a process from the head matrix element,
* return a list of diagrams which should be considered for
* the given dependent matrix element.
*/
virtual MEBase::DiagramVector dependentDiagrams(const cPDVector& proc,
tMEPtr depME) const;
/**
* Return true, if SubProcessGroups should be
* setup from this MEGroup. If not, a single SubProcess
* is constructed from the data provided by the
* head matrix element.
*/
virtual bool subProcessGroups() const;
/**
* Return true, if one of the dependent subprocesses should be
* constructed in place of the one driven by the head matrix element
* or a full subprocess group.
*/
virtual bool selectDependentSubProcess() const { return false; }
/**
* Fill the projectors object of xcombs to choose subprocesses
* different than the one currently integrated.
*/
virtual void fillProjectors();
/**
* Return true, if projectors will be used
*/
virtual bool willProject() const {
return virtualShowerSubtraction() || loopSimSubtraction();
}
/**
* Return true, if this MEGroup will reweight the contributing cross
* sections.
*/
virtual bool groupReweighted() const {
return showerApproximation();
}
/**
* Reweight the head cross section
*/
virtual double reweightHead(const vector<tStdXCombPtr>&);
/**
* Reweight the dependent cross section
*/
virtual double reweightDependent(tStdXCombPtr, const vector<tStdXCombPtr>&);
/**
* Switch on or off that scales should be calculated from real emission kinematics
*/
void doRealEmissionScales();
//@}
/** @name Methods relevant to matching */
//@{
/**
* Inform this matrix element that a new phase space
* point is about to be generated, so all caches should
* be flushed.
*/
virtual void flushCaches() {
MEGroup::flushCaches();
if ( showerApproximation() )
showerApproximation()->resetBelowCutoff();
}
/**
* Return the shower approximation.
*/
Ptr<ShowerApproximation>::tptr showerApproximation() const;
/**
* Indicate that the shower real emission contribution should be subtracted.
*/
void doRealShowerSubtraction();
/**
* Return true, if the shower real emission contribution should be subtracted.
*/
bool realShowerSubtraction() const { return theRealShowerSubtraction; }
/**
* Indicate that the shower virtual contribution should be subtracted.
*/
void doVirtualShowerSubtraction();
/**
* Return true, if the shower virtual contribution should be subtracted.
*/
bool virtualShowerSubtraction() const { return theVirtualShowerSubtraction; }
/**
* Indicate that the loopsim matched virtual contribution should be subtracted.
*/
void doLoopSimSubtraction();
/**
* Return true, if the loopsim matched virtual contribution should be subtracted.
*/
bool loopSimSubtraction() const { return theLoopSimSubtraction; }
+ /**
+ * Return true, if this configuration of cross sections should not
+ * be included due to their relative magnitude. Arguments are head
+ * cross section and dependent cross section, including all
+ * reweights.
+ */
+ virtual bool discard(const CrossSection&, const CrossSection&) const { return false; }
+
//@}
/** @name Matrix element and dipole information */
//@{
/**
* Return the subtraction dipoles.
*/
vector<Ptr<SubtractionDipole>::ptr> dipoles();
/**
* Return the underlying born matrix elements.
*/
const vector<Ptr<MatchboxMEBase>::ptr>& borns() const;
/**
* Access the underlying born matrix elements,
* overriding the ones contained in the factory object.
*/
void setBorns(const vector<Ptr<MatchboxMEBase>::ptr>& newBorns) { theBorns = newBorns; }
/**
* Build up dipoles needed.
*/
void getDipoles();
/**
* Clone all dipoles.
*/
void cloneDipoles(const string& prefix = "");
/**
* Clone the real emission matrix element.
*/
void cloneRealME(const string& prefix = "");
/**
* Clone all dependencies.
*/
void cloneDependencies(const string& prefix = "") {
cloneDipoles(prefix);
cloneRealME(prefix);
}
/**
* Return all dipoles matching the given Born process
*/
vector<Ptr<SubtractionDipole>::ptr> splitDipoles(const cPDVector&);
//@}
/** @name Veto scale settings */
//@{
/**
* Set veto scales on the particles at the given
* SubProcess which has been generated using this
* matrix element.
*/
virtual void setVetoScales(tSubProPtr) const;
//@}
/** @name Diagnostic information */
//@{
/**
* Dump the setup to an ostream
*/
void print(ostream&) const;
/**
* Collect information on the last evaluated phasespace
* point for verification or debugging purposes. This
* only called, if the StdXCombGroup did accumulate
* a non-zero cross section from this ME group.
*/
virtual void lastEventStatistics();
/**
* Print debug information on the last event
*/
void printLastEvent(ostream&) const;
/**
* Check the subtraction for the last event
*/
void lastEventSubtraction();
/**
* Return true, if verbose
*/
bool verbose() const;
/**
* Return true, if verbose
*/
bool initVerbose() const;
//@}
/** @name Setup of Subtracted ME objects */
//@{
/**
* Return true if this object needs to be initialized before all
* other objects (except those for which this function also returns
* true). This default version always returns false, but subclasses
* may override it to return true.
*/
virtual bool preInitialize() const { return true; }
/**
* Simple envelope histogram to keep track of subtraction
*/
struct SubtractionHistogram {
/**
* The lower bound
*/
double lower;
/**
* The bins, indexed by upper bound.
*/
map<double,pair<double,double> > bins;
/**
* Constructor
*/
SubtractionHistogram(double low = 0.001,
double up = 10.,
unsigned int nbins = 100);
/**
* Book an event.
*/
void book(double inv, double diff) {
map<double,pair<double,double> >::iterator b =
bins.upper_bound(inv);
if ( b == bins.end() ) return;
b->second.first = min(b->second.first,diff);
b->second.second = max(b->second.second,diff);
}
/**
* Write to file given name and invariant.
*/
void dump(const std::string& prefix,
const int& plottype,
const bool& scatterplot,
const cPDVector& proc,
int i, int j) const;
/**
* Write to persistent ostream
*/
void persistentOutput(PersistentOStream&) const;
/**
* Read from persistent istream
*/
void persistentInput(PersistentIStream&);
};
//@}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
/**
* Finalize this object. Called in the run phase just after a
* run has ended. Used eg. to write out statistics.
*/
virtual void dofinish();
//@}
private:
/**
* The underlying born matrix elements, overriding the ones
* contained in the factory object.
*/
vector<Ptr<MatchboxMEBase>::ptr> theBorns;
/**
* Pointer to the head real emission ME casted to a MatchboxMEBase
* object.
*/
Ptr<MatchboxMEBase>::ptr theReal;
/**
* Define the key for the collinear subtraction data.
*/
typedef pair<cPDVector,pair<size_t, size_t> > CollinearSubtractionIndex;
/**
* subtraction data for collinear limits.
*/
map<CollinearSubtractionIndex,SubtractionHistogram> collinearHistograms;
/**
* names of files to which subtraction data is written for all phase space points individually
*/
map<CollinearSubtractionIndex,string> fnamesCollinearSubtraction;
/**
* Define the key for the soft subtraction data.
*/
typedef pair<cPDVector,size_t> SoftSubtractionIndex;
/**
* subtraction data for soft limits.
*/
map<SoftSubtractionIndex,SubtractionHistogram> softHistograms;
/**
* names of files to which subtraction data is written for all phase space points individually
*/
map<SoftSubtractionIndex,string> fnamesSoftSubtraction;
/**
* True, if the shower real emission contribution should be subtracted.
*/
bool theRealShowerSubtraction;
/**
* True, if the shower virtual contribution should be subtracted.
*/
bool theVirtualShowerSubtraction;
/**
* True, if the loopsim matched virtual contribution should be subtracted.
*/
bool theLoopSimSubtraction;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
SubtractedME & operator=(const SubtractedME &) = delete;
};
inline PersistentOStream& operator<<(PersistentOStream& os,
const SubtractedME::SubtractionHistogram& h) {
h.persistentOutput(os);
return os;
}
inline PersistentIStream& operator>>(PersistentIStream& is,
SubtractedME::SubtractionHistogram& h) {
h.persistentInput(is);
return is;
}
}
#endif /* HERWIG_SubtractedME_H */
diff --git a/MatrixElement/Matchbox/External/VBFNLO/VBFNLOPhasespace.cc b/MatrixElement/Matchbox/External/VBFNLO/VBFNLOPhasespace.cc
--- a/MatrixElement/Matchbox/External/VBFNLO/VBFNLOPhasespace.cc
+++ b/MatrixElement/Matchbox/External/VBFNLO/VBFNLOPhasespace.cc
@@ -1,275 +1,277 @@
// -*- C++ -*-
//
// VBFNLOPhasespace.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the VBFNLOPhasespace class.
//
#include "VBFNLOPhasespace.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "Herwig/Utilities/GSLBisection.h"
#include "ThePEG/Utilities/DynamicLoader.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/MatrixElement/Matchbox/Base/MatchboxAmplitude.h"
#include "VBFNLO/utilities/BLHAinterface.h"
#define DEFSTR(s) CPPSTR(s)
#define CPPSTR(s) #s
using namespace Herwig;
VBFNLOPhasespace::VBFNLOPhasespace() :
lastSqrtS(0*GeV), needToReshuffle(false), VBFNLOlib_(DEFSTR(VBFNLOLIB))
{}
void VBFNLOPhasespace::loadVBFNLO() {
if ( ! DynamicLoader::load(VBFNLOlib_+"/libVBFNLO.so") ) {
string error1 = DynamicLoader::lastErrorMessage;
if ( ! DynamicLoader::load(VBFNLOlib_+"/libVBFNLO.dylib") ) {
string error2 = DynamicLoader::lastErrorMessage;
if ( ! DynamicLoader::load("libVBFNLO.so") ) {
string error3 = DynamicLoader::lastErrorMessage;
if ( ! DynamicLoader::load("libVBFNLO.dylib") ) {
string error4 = DynamicLoader::lastErrorMessage;
throw Exception() << "VBFNLOPhasespace: failed to load libVBFNLO.so/dylib\n"
<< "Error messages are:\n\n"
<< "* " << VBFNLOlib_ << "/libVBFNLO.so:\n"
<< error1 << "\n"
<< "* " << VBFNLOlib_ << "/libVBFNLO.dylib:\n"
<< error2 << "\n"
<< "* libVBFNLO.so:\n"
<< error3 << "\n"
<< "* libVBFNLO.dylib:\n"
<< error4 << "\n"
<< Exception::runerror;
}
}
}
}
}
VBFNLOPhasespace::~VBFNLOPhasespace() {}
IBPtr VBFNLOPhasespace::clone() const {
return new_ptr(*this);
}
IBPtr VBFNLOPhasespace::fullclone() const {
return new_ptr(*this);
}
void VBFNLOPhasespace::setXComb(tStdXCombPtr xco) {
MatchboxPhasespace::setXComb(xco);
// test for resuffling
needToReshuffle = false;
if ( xco ) {
for ( cPDVector::const_iterator d = mePartonData().begin();
d != mePartonData().end(); ++d ) {
// Higgs is massive -> does not need reshuffling
if ( ( (**d).id() != ParticleID::h0 ) && ( (**d).hardProcessMass() != ZERO ) ) {
needToReshuffle = true;
break;
}
}
}
// set CMS energy
int pStatus = 0;
double zero = 0.0;
double value = sqrt(lastXCombPtr()->lastS())/GeV;
if (value && (value != lastSqrtS/GeV)) {
lastSqrtS = value*GeV;
string name = "sqrtS";
OLP_SetParameter(const_cast<char*>(name.c_str()),&value,&zero,&pStatus);
if ( !pStatus )
throw Exception() << "VBFNLOPhasespace::setXComb(): VBFNLO failed to set parameter '"
<< name << "' to " << value << "\n"
<< Exception::runerror;
}
}
double VBFNLOPhasespace::generateTwoToNKinematics(const double* random,
vector<Lorentz5Momentum>& momenta) {
double weight;
int id =
olpId()[ProcessType::oneLoopInterference] ?
olpId()[ProcessType::oneLoopInterference] :
olpId()[ProcessType::treeME2];
double* p = new double[4*momenta.size()];
OLP_PhaseSpacePoint(&id, const_cast<double*>(random), const_cast<double*>(random+1), p, &weight);
if (weight < 0) {
throw Exception() << "VBFNLOPhasespace::generateTwoToNKinematics(): Negative weight in VBFNLOPhaseSpace\n"
<< Exception::runerror;
}
if (weight == 0) {
delete[] p;
return 0;
}
for ( size_t i = 0; i < momenta.size(); ++i ) {
momenta[i].setT(p[4*i] *GeV);
momenta[i].setX(p[4*i+1]*GeV);
momenta[i].setY(p[4*i+2]*GeV);
momenta[i].setZ(p[4*i+3]*GeV);
momenta[i].rescaleMass();
}
delete[] p;
Energy beamenergy = sqrt(lastXCombPtr()->lastS())/2.;
double x1 = momenta[0].e()/beamenergy;
double x2 = momenta[1].e()/beamenergy;
Energy2 thisSHat = (momenta[0] + momenta[1]).m2();
// reshuffle so that particles have correct mass
if ( needToReshuffle ) {
// boost final-state into partonic CMS
Boost toCMS = (momenta[0]+momenta[1]).findBoostToCM();
for ( size_t i = 2; i < momenta.size(); ++i ) {
momenta[i].boost(toCMS);
}
// copied from MatchboxRambo phasespace
double xi;
ReshuffleEquation solve(sqrt(thisSHat),mePartonData().begin()+2,mePartonData().end(),
momenta.begin()+2,momenta.end());
GSLBisection solver(1e-10,1e-8,10000);
try {
xi = solver.value(solve,0.0,1.1);
} catch (GSLBisection::GSLerror) {
return 0.;
} catch (GSLBisection::IntervalError) {
return 0.;
}
weight *= pow(xi,3.*(momenta.size()-3.));
Energy num = ZERO;
Energy den = ZERO;
cPDVector::const_iterator d = mePartonData().begin()+2;
for ( vector<Lorentz5Momentum>::iterator k = momenta.begin()+2;
k != momenta.end(); ++k, ++d ) {
num += (*k).vect().mag2()/(*k).t();
Energy q = (*k).t();
(*k).setT(sqrt(sqr((**d).hardProcessMass())+xi*xi*sqr((*k).t())));
(*k).setVect(xi*(*k).vect());
weight *= q/(*k).t();
den += (*k).vect().mag2()/(*k).t();
(*k).setMass((**d).hardProcessMass());
}
// unboost
for ( size_t i = 2; i < momenta.size(); ++i ) {
momenta[i].boost(-toCMS);
}
}
if ( !matchConstraints(momenta) )
return 0.;
lastXCombPtr()->lastX1X2(make_pair(x1,x2));
lastXCombPtr()->lastSHat(thisSHat);
weight /= pow(thisSHat/GeV2,momenta.size()-4);
weight /= x1*x2;
+ fillDiagramWeights();
+
return weight;
}
int VBFNLOPhasespace::nDimPhasespace(int nFinal) const {
return 3*nFinal;
//get this from within VBFNLO
int pStatus = 0;
double value, zero;
string name = "PSdimension";
OLP_GetParameter(const_cast<char*>(name.c_str()),&value,&zero,&pStatus);
if ( pStatus != 1) {
throw Exception() << "VBFNLOPhasespace::nDimPhasespace(): Cannot get phasespace dimension in VBFNLOPhaseSpace\n"
<< "error code: " << pStatus << "\n"
<< Exception::runerror;
}
// one additional number (first) needed for channel selection
// one additional number (last) needed for global phi integration
return value+2;
}
Energy VBFNLOPhasespace::ReshuffleEquation::operator() (double xi) const {
cPDVector::const_iterator d = dataBegin;
vector<Lorentz5Momentum>::const_iterator p = momentaBegin;
Energy res = -w;
for ( ; d != dataEnd; ++d, ++p ) {
res += sqrt(sqr((**d).hardProcessMass()) +
xi*xi*sqr(p->t()));
}
return res;
}
void VBFNLOPhasespace::doinit() {
loadVBFNLO();
MatchboxPhasespace::doinit();
}
void VBFNLOPhasespace::doinitrun() {
loadVBFNLO();
MatchboxPhasespace::doinitrun();
}
void VBFNLOPhasespace::persistentOutput(PersistentOStream & os) const {
os << needToReshuffle << theLastXComb;
}
void VBFNLOPhasespace::persistentInput(PersistentIStream & is, int) {
is >> needToReshuffle >> theLastXComb;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<VBFNLOPhasespace,MatchboxPhasespace>
describeHerwigVBFNLOPhasespace("Herwig::VBFNLOPhasespace", "HwMatchboxVBFNLO.so");
void VBFNLOPhasespace::Init() {
static ClassDocumentation<VBFNLOPhasespace> documentation
("VBFNLOPhasespace is an interface to the internal phasespace generator "
"of VBFNLO. It uses the information passed via the BLHA interface to "
"obtain information on the required channels.");
}
diff --git a/m4/herwig.m4 b/m4/herwig.m4
--- a/m4/herwig.m4
+++ b/m4/herwig.m4
@@ -1,983 +1,1003 @@
dnl ##### THEPEG #####
AC_DEFUN([HERWIG_CHECK_THEPEG],
[
defaultlocation="${prefix}"
test "x$defaultlocation" = xNONE && defaultlocation="${ac_default_prefix}"
AC_MSG_CHECKING([for libThePEG in])
AC_ARG_WITH(thepeg,
AC_HELP_STRING([--with-thepeg=DIR],[location of ThePEG installation]),
[],
[with_thepeg="${defaultlocation}"])
AC_MSG_RESULT([$with_thepeg])
if test "x$with_thepeg" = "xno"; then
AC_MSG_ERROR([Cannot build Herwig without ThePEG. Please set --with-thepeg.])
fi
THEPEGLDFLAGS="-L${with_thepeg}/lib/ThePEG"
THEPEGHASLHAPDF="no"
if test -e ${with_thepeg}/lib/ThePEG/ThePEGLHAPDF.so ; then
THEPEGHASLHAPDF="yes"
fi
if test "${host_cpu}" == "x86_64" -a -e ${with_thepeg}/lib64/ThePEG/libThePEG.so ; then
THEPEGLDFLAGS="-L${with_thepeg}/lib64/ThePEG"
if test -e ${with_thepeg}/lib64/ThePEG/ThePEGLHAPDF.so ; then
THEPEGHASLHAPDF="yes"
fi
fi
if test "x$THEPEGHASLHAPDF" == "xno" ; then
AC_MSG_ERROR([Herwig requires ThePEG to be build with lhapdf.])
fi
THEPEGHASFASTJET="no"
if test -e ${with_thepeg}/lib/ThePEG/FastJetFinder.so ; then
THEPEGHASFASTJET="yes"
fi
if test "${host_cpu}" == "x86_64" -a -e ${with_thepeg}/lib64/ThePEG/libThePEG.so ; then
THEPEGLDFLAGS="-L${with_thepeg}/lib64/ThePEG"
if test -e ${with_thepeg}/lib64/ThePEG/FastJetFinder.so ; then
THEPEGHASFASTJET="yes"
fi
fi
if test "x$THEPEGHASFASTJET" == "xno" ; then
AC_MSG_ERROR([Herwig requires ThePEG to be build with FastJet.])
fi
THEPEGPATH="${with_thepeg}"
oldldflags="$LDFLAGS"
oldlibs="$LIBS"
LDFLAGS="$LDFLAGS $THEPEGLDFLAGS"
AC_CHECK_LIB([ThePEG],[debugThePEG],[],
[AC_MSG_ERROR([No ThePEG libraries in $THEPEGLDFLAGS. Please set --with-thepeg.])])
AC_SUBST([THEPEGLIB],[-lThePEG])
AC_SUBST(THEPEGLDFLAGS)
AC_SUBST(THEPEGPATH)
AC_SUBST(THEPEGHASLHAPDF)
AC_SUBST(THEPEGHASFASTJET)
LIBS="$oldlibs"
LDFLAGS="$oldldflags"
AC_MSG_CHECKING([for ThePEG headers in])
AC_ARG_WITH([thepeg-headers],
AC_HELP_STRING([--with-thepeg-headers=DIR],[location of ThePEG include directory]),
[],
[with_thepeg_headers="${with_thepeg}/include"])
AC_MSG_RESULT([$with_thepeg_headers])
if test "x$with_thepeg_headers" = "xno"; then
AC_MSG_ERROR([Cannot build Herwig without ThePEG headers. Please set --with-thepeg-headers.])
fi
THEPEGINCLUDE="-I$with_thepeg_headers"
oldcppflags="$CPPFLAGS"
CPPFLAGS="$CPPFLAGS $THEPEGINCLUDE"
AC_CHECK_HEADER([ThePEG/Config/ThePEG.h],[],
[AC_MSG_ERROR([No ThePEG headers in $with_thepeg_headers. Please set --with-thepeg-headers.])])
CPPFLAGS="$oldcppflags"
AC_SUBST(THEPEGINCLUDE)
AC_MSG_CHECKING([for HepMCAnalysis.so in ThePEG])
THEPEGHASHEPMC="no"
if test -e ${with_thepeg}/lib/ThePEG/HepMCAnalysis.so ; then
THEPEGHASHEPMC="yes"
fi
if test "${host_cpu}" == "x86_64" -a -e ${with_thepeg}/lib64/ThePEG/libThePEG.so ; then
THEPEGLDFLAGS="-L${with_thepeg}/lib64/ThePEG"
if test -e ${with_thepeg}/lib64/ThePEG/HepMCAnalysis.so ; then
THEPEGHASHEPMC="yes"
fi
fi
if test "x$THEPEGHASHEPMC" == "xno" ; then
CREATE_HEPMC="# create"
AC_MSG_RESULT([not found])
else
CREATE_HEPMC="create"
AC_MSG_RESULT([found])
fi
AC_SUBST([CREATE_HEPMC])
AC_MSG_CHECKING([for RivetAnalysis.so in ThePEG])
THEPEGHASRIVET="no"
if test -e ${with_thepeg}/lib/ThePEG/RivetAnalysis.so ; then
THEPEGHASRIVET="yes"
fi
if test "${host_cpu}" == "x86_64" -a -e ${with_thepeg}/lib64/ThePEG/libThePEG.so ; then
THEPEGLDFLAGS="-L${with_thepeg}/lib64/ThePEG"
if test -e ${with_thepeg}/lib64/ThePEG/RivetAnalysis.so ; then
THEPEGHASRIVET="yes"
fi
fi
if test "x$THEPEGHASRIVET" == "xno" ; then
CREATE_RIVET="# create"
AC_MSG_RESULT([not found])
else
CREATE_RIVET="create"
AC_MSG_RESULT([found])
fi
AM_CONDITIONAL(HAVE_RIVET,[test x$THEPEGHASRIVET = xyes])
AC_SUBST([CREATE_RIVET])
])
dnl ##### LOOPTOOLS #####
AC_DEFUN([HERWIG_LOOPTOOLS],
[
AC_REQUIRE([AC_PROG_FC])
AC_REQUIRE([AC_FC_LIBRARY_LDFLAGS])
AC_REQUIRE([AC_PROG_CC])
AC_REQUIRE([HERWIG_COMPILERFLAGS])
AC_MSG_CHECKING([if Looptools build works])
enable_looptools=yes
if test "x$GCC" = "xyes"; then
case "${host}" in
x86_64-*|*-darwin1*)
AM_FCFLAGS="$AM_FCFLAGS -fdefault-integer-8"
;;
esac
AC_LANG_PUSH([Fortran])
oldFCFLAGS="$FCFLAGS"
FCFLAGS="$AM_FCFLAGS"
AC_COMPILE_IFELSE(
AC_LANG_PROGRAM([],[ print *[,]"Hello"]),
[],
[AC_MSG_RESULT([no])
AC_MSG_ERROR([needs gfortran on 64bit machines])]
)
FCFLAGS="$oldFCFLAGS"
AC_LANG_POP([Fortran])
fi
AC_MSG_RESULT([$enable_looptools])
AC_LANG_PUSH([Fortran])
AC_MSG_CHECKING([checking if fortran compiler compiles argument mismatches])
AC_COMPILE_IFELSE(AC_LANG_SOURCE([[
program temp
call a(1.0D0)
end program
subroutine a(b)
double complex b
end]]),
[AC_MSG_RESULT([yes])],
[oldFCFLAGS="$FCFLAGS"
FCFLAGS="-std=legacy"
AC_MSG_CHECKING([checking if fortran compiler compiles argument mismatches with -std=legacy])
AC_COMPILE_IFELSE(AC_LANG_SOURCE([[
program temp
double precision b
double complex c
b = 1.0D0
c =1.0D0
call a(b)
call a(c)
end program
subroutine a(b)
double complex b
end]]),
[AC_MSG_RESULT([yes])
AM_FCFLAGS="$AM_FCFLAGS -std=legacy"],
[AC_MSG_RESULT([no])
AC_MSG_ERROR([fortran compiler won't compile LoopTools])])
FCFLAGS="$oldFCFLAGS -std=legacy"]
)
AC_MSG_CHECKING([checking if fortran compiler compiles long lines])
AC_COMPILE_IFELSE(AC_LANG_SOURCE([[
program temp
write (*,*) 'aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa'
end program]]),
[AC_MSG_RESULT([yes])],
[AC_MSG_RESULT([no])
oldFCFLAGS="$FCFLAGS"
FCFLAGS="-ffixed-line-length-none"
AC_MSG_CHECKING([checking if fortran compiler compiles long lines with -ffixed-line-length-none])
AC_COMPILE_IFELSE(AC_LANG_SOURCE([[
program temp
write (*,*) 'aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa'
end program]]),
[AC_MSG_RESULT([yes])
AM_FCFLAGS="$AM_FCFLAGS -ffixed-line-length-none"
FCFLAGS="$oldFCFLAGS -ffixed-line-length-none"],
[FCFLAGS="-132"
AC_MSG_CHECKING([checking if fortran compiler compiles long lines with -132])
AC_COMPILE_IFELSE(AC_LANG_SOURCE([[
program temp
write (*,*) 'aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa'
end program]]),
[AC_MSG_RESULT([yes])
AM_FCFLAGS="$AM_FCFLAGS -132"
FCFLAGS="$oldFCFLAGS -132"],
[AC_MSG_RESULT([no])
AC_MSG_ERROR([fortran compiler won't compile LoopTools])])])
]
)
AC_LANG_POP([Fortran])
AC_SUBST([F77],[$FC])
AC_SUBST([FFLAGS],[$FCFLAGS])
AC_SUBST([AM_FFLAGS],[$AM_FCFLAGS])
AC_SUBST([FLIBS],[$FCLIBS])
])
dnl ##### VBFNLO #####
AC_DEFUN([HERWIG_CHECK_VBFNLO],
[
AC_MSG_CHECKING([for VBFNLO])
AC_ARG_WITH([vbfnlo],
AS_HELP_STRING([--with-vbfnlo=DIR], [Installation path of VBFNLO]),
[],
[with_vbfnlo=no]
)
AC_MSG_RESULT([$with_vbfnlo])
AS_IF([test "x$with_vbfnlo" != "xno"],
[AC_CHECK_FILES(
${with_vbfnlo}/lib/VBFNLO/libVBFNLO.so,
[have_vbfnlo=lib], [have_vbfnlo=no])],
[have_vbfnlo=no])
AS_IF([test "x$with_vbfnlo" != "xno" -a "x$have_vbfnlo" = "xno" ],
[AC_CHECK_FILES(
${with_vbfnlo}/lib64/VBFNLO/libVBFNLO.so,
[have_vbfnlo=lib64], [have_vbfnlo=no])])
AS_IF([test "x$with_vbfnlo" != "xno" -a "x$have_vbfnlo" = "xno" ],
[AC_CHECK_FILES(
${with_vbfnlo}/lib/VBFNLO/libVBFNLO.dylib,
[have_vbfnlo=lib], [have_vbfnlo=no])])
AS_IF([test "x$with_vbfnlo" != "xno" -a "x$have_vbfnlo" = "xno" ],
[AC_CHECK_FILES(
${with_vbfnlo}/lib64/VBFNLO/libVBFNLO.dylib,
[have_vbfnlo=lib64], [have_vbfnlo=no])])
AS_IF([test "x$have_vbfnlo" = "xlib"],
[VBFNLOLIBS=${with_vbfnlo}/lib/VBFNLO
AC_SUBST(VBFNLOLIBS)
])
AS_IF([test "x$have_vbfnlo" = "xlib64"],
[VBFNLOLIBS=${with_vbfnlo}/lib64/VBFNLO
AC_SUBST(VBFNLOLIBS)
])
AS_IF([test "x$with_vbfnlo" != "xno" -a "x$have_vbfnlo" = "xno"],
[AC_MSG_ERROR([vbfnlo requested but not found])])
AM_CONDITIONAL(HAVE_VBFNLO,[test "x$have_vbfnlo" = "xlib" -o "x$have_vbfnlo" = "xlib64"])
if test "x$have_vbfnlo" = "xlib" -o "x$have_vbfnlo" = "xlib64" ; then
AC_REQUIRE([AC_PROG_SED])
VBFNLOINCLUDE=${with_vbfnlo}/include
AC_SUBST(VBFNLOINCLUDE)
VBFNLOLIB=$(echo ${with_vbfnlo}/${have_vbfnlo}/VBFNLO | $SED -e 's%/\+%/%g')
AC_SUBST(VBFNLOLIB)
LOAD_VBFNLO="library"
CREATE_VBFNLO="create"
INSERT_VBFNLO="insert"
SET_VBFNLO="set"
DO_VBFNLO="do"
MKDIR_VBFNLO="mkdir"
else
LOAD_VBFNLO="# library"
CREATE_VBFNLO="# create"
INSERT_VBFNLO="# insert"
SET_VBFNLO="# set"
DO_VBFNLO="# do"
MKDIR_VBFNLO="# mkdir"
fi
AC_SUBST([LOAD_VBFNLO])
AC_SUBST([CREATE_VBFNLO])
AC_SUBST([INSERT_VBFNLO])
AC_SUBST([SET_VBFNLO])
AC_SUBST([DO_VBFNLO])
AC_SUBST([MKDIR_VBFNLO])
])
dnl ##### njet #####
AC_DEFUN([HERWIG_CHECK_NJET],
[
AC_MSG_CHECKING([for njet])
AC_ARG_WITH([njet],
AS_HELP_STRING([--with-njet=DIR], [Installation path of njet]),
[],
[with_njet=no]
)
AC_MSG_RESULT([$with_njet])
AS_IF([test "x$with_njet" != "xno"],
[AC_CHECK_FILES(
${with_njet}/lib/libnjet2.so,
[have_njet=lib], [have_njet=no])],
[have_njet=no])
AS_IF([test "x$with_njet" != "xno" -a "x$have_njet" = "xno" ],
[AC_CHECK_FILES(
${with_njet}/lib64/libnjet2.so,
[have_njet=lib64], [have_njet=no])])
AS_IF([test "x$with_njet" != "xno" -a "x$have_njet" = "xno" ],
[AC_CHECK_FILES(
${with_njet}/lib/libnjet2.dylib,
[have_njet=lib], [have_njet=no])])
+AS_IF([test "x$with_njet" != "xno" -a "x$have_njet" = "xno" ],
+ [AC_CHECK_FILES(
+ ${with_njet}/lib/libnjet3.so,
+ [have_njet=lib], [have_njet=no])])
+
+AS_IF([test "x$with_njet" != "xno" -a "x$have_njet" = "xno" ],
+ [AC_CHECK_FILES(
+ ${with_njet}/lib64/libnjet3.so,
+ [have_njet=lib], [have_njet=no])])
+
+AS_IF([test "x$with_njet" != "xno" -a "x$have_njet" = "xno" ],
+ [AC_CHECK_FILES(
+ ${with_njet}/lib/libnjet3.dylib,
+ [have_njet=lib], [have_njet=no])])
+
+AS_IF([test "x$with_njet" != "xno" ],
+ [AC_CHECK_FILES(
+ ${with_njet}/include/njet.h,
+ [njet_include=include], [njet_include=include/njet])])
+
AS_IF([test "x$have_njet" = "xlib"],
[NJETLIBPATH=${with_njet}/lib
AC_SUBST(NJETLIBPATH)
- NJETINCLUDEPATH=${with_njet}/include
+ NJETINCLUDEPATH=${with_njet}/$njet_include
AC_SUBST(NJETINCLUDEPATH)
NJETPREFIX=${with_njet}
AC_SUBST(NJETPREFIX)
])
AS_IF([test "x$have_njet" = "xlib64"],
[NJETLIBPATH=${with_njet}/lib64
AC_SUBST(NJETLIBPATH)
- NJETINCLUDEPATH=${with_njet}/include
+ NJETINCLUDEPATH=${with_njet}/$njet_include
AC_SUBST(NJETINCLUDEPATH)
NJETPREFIX=${with_njet}
AC_SUBST(NJETPREFIX)
])
AS_IF([test "x$with_njet" != "xno" -a "x$have_njet" = "xno"],
[AC_MSG_ERROR([njet requested but not found])])
AM_CONDITIONAL(HAVE_NJET,[test "x$have_njet" = "xlib" -o "x$have_njet" = "xlib64"])
AS_IF([test "x$with_njet" != "xno"],[AC_MSG_CHECKING([for Njet version])
tmp_njetversion=[$(${with_njet}/bin/njet.py --version 2>&1 | grep -oE '[0-9]{4}$')]
AS_IF([test -z "$tmp_njetversion"],
[tmp_njetversion=1023])
AC_MSG_RESULT([${tmp_njetversion}])
NJET_VERSION=${tmp_njetversion}
AC_SUBST(NJET_VERSION)
])
if test "x$have_njet" = "xlib" -o "x$have_njet" = "xlib64" ; then
LOAD_NJET="library"
CREATE_NJET="create"
INSERT_NJET="insert"
DO_NJET="do"
else
LOAD_NJET="# library"
CREATE_NJET="# create"
INSERT_NJET="# insert"
DO_NJET="# do"
fi
AC_SUBST([LOAD_NJET])
AC_SUBST([CREATE_NJET])
AC_SUBST([INSERT_NJET])
AC_SUBST([DO_NJET])
])
dnl ##### gosam #####
AC_DEFUN([HERWIG_CHECK_GOSAM],
[
AC_MSG_CHECKING([for GoSam])
AC_ARG_WITH([gosam],
AS_HELP_STRING([--with-gosam=DIR], [Installation path of GoSam]),
[],
[with_gosam=no]
)
AC_MSG_RESULT([$with_gosam])
AS_IF([test "x$with_gosam" != "xno"],
[AC_CHECK_FILES(
${with_gosam}/bin/gosam.py,
[have_gosam=lib], [have_gosam=no])],
[have_gosam=no])
AS_IF([test "x$have_gosam" = "xlib"],
[GOSAMPREFIX=${with_gosam}
AC_SUBST(GOSAMPREFIX)
])
AS_IF([test "x$with_gosam" != "xno" -a "x$have_gosam" = "xno"],
[AC_MSG_ERROR([GoSam requested but not found])])
AS_IF([test "x$with_gosam" != "xno"],
[AC_MSG_CHECKING([for GoSam version >= 2.0.4])
tmp_gosamversion=[$(${with_gosam}/bin/gosam.py --version | grep 'GoSam.*rev' | cut -d' ' -f2)]
AX_COMPARE_VERSION([${tmp_gosamversion}],[lt],[2.0.4],
[AC_MSG_RESULT([no])
AC_MSG_ERROR([Herwig requires GoSam 2.0.4 or later, found ${tmp_gosamversion}])],
[AC_MSG_RESULT([yes])])])
AM_CONDITIONAL(HAVE_GOSAM,[test "x$have_gosam" = "xlib" ])
if test "x$have_gosam" = "xlib" ; then
LOAD_GOSAM="library"
CREATE_GOSAM="create"
INSERT_GOSAM="insert"
DO_GOSAM="do"
else
LOAD_GOSAM="# library"
CREATE_GOSAM="# create"
INSERT_GOSAM="# insert"
DO_GOSAM="# do"
fi
AC_SUBST([LOAD_GOSAM])
AC_SUBST([CREATE_GOSAM])
AC_SUBST([INSERT_GOSAM])
AC_SUBST([DO_GOSAM])
])
dnl ##### gosam-contrib #####
AC_DEFUN([HERWIG_CHECK_GOSAM_CONTRIB],
[
AC_MSG_CHECKING([for gosam-contrib])
AC_ARG_WITH([gosam-contrib],
AS_HELP_STRING([--with-gosam-contrib=DIR], [Installation path of gosam-contrib]),
[],
[with_gosam_contrib=no]
)
AC_MSG_RESULT([$with_gosam_contrib])
AS_IF([test "x$with_gosam_contrib" = "xno" -a "x$with_gosam" != "xno"],
[AC_CHECK_FILES(
${with_gosam}/lib/libsamurai.so,
[with_gosam_contrib=${with_gosam}], [])
])
AS_IF([test "x$with_gosam_contrib" = "xno" -a "x$with_gosam" != "xno"],
[AC_CHECK_FILES(
${with_gosam}/lib64/libsamurai.so,
[with_gosam_contrib=${with_gosam}], [])
])
AS_IF([test "x$with_gosam_contrib" = "xno" -a "x$with_gosam" != "xno"],
[AC_CHECK_FILES(
${with_gosam}/lib/libsamurai.dylib,
[with_gosam_contrib=${with_gosam}], [])
])
AS_IF([test "x$with_gosam_contrib" = "xno" -a "x$with_gosam" != "xno"],
[AC_CHECK_FILES(
${with_gosam}/lib64/libsamurai.dylib,
[with_gosam_contrib=${with_gosam}], [])
])
AS_IF([test "x$with_gosam_contrib" = "xno" -a "x$with_gosam" != "xno"],
[AC_MSG_ERROR([GoSam requested without requesting GoSam-Contrib])])
AS_IF([test "x$with_gosam_contrib" != "xno"],
[AC_CHECK_FILES(
${with_gosam_contrib}/lib/libsamurai.so,
[have_gosam_contrib=lib], [have_gosam_contrib=no])],
[have_gosam_contrib=no])
AS_IF([test "x$with_gosam_contrib" != "xno" -a "x$have_gosam_contrib" = "xno" ],
[AC_CHECK_FILES(
${with_gosam_contrib}/lib64/libsamurai.so,
[have_gosam_contrib=lib64], [have_gosam_contrib=no])])
AS_IF([test "x$with_gosam_contrib" != "xno" -a "x$have_gosam_contrib" = "xno" ],
[AC_CHECK_FILES(
${with_gosam_contrib}/lib/libsamurai.dylib,
[have_gosam_contrib=lib], [have_gosam_contrib=no])])
AS_IF([test "x$with_gosam_contrib" != "xno" -a "x$have_gosam_contrib" = "xno" ],
[AC_CHECK_FILES(
${with_gosam_contrib}/lib64/libsamurai.dylib,
[have_gosam_contrib=lib64], [have_gosam_contrib=no])])
AS_IF([test "x$have_gosam_contrib" != "xno"],
[GOSAMCONTRIBPREFIX=${with_gosam_contrib}
AC_SUBST(GOSAMCONTRIBPREFIX)
])
AS_IF([test "x$have_gosam_contrib" = "xlib"],
[GOSAMCONTRIBLIBS=${with_gosam_contrib}/lib
AC_SUBST(GOSAMCONTRIBLIBS)
])
AS_IF([test "x$have_gosam_contrib" = "xlib64"],
[GOSAMCONTRIBLIBS=${with_gosam_contrib}/lib64
AC_SUBST(GOSAMCONTRIBLIBS)
])
AS_IF([test "x$with_gosam_contrib" != "xno" -a "x$have_gosam_contrib" = "xno"],
[AC_MSG_ERROR([GoSam-Contrib requested but not found])])
AM_CONDITIONAL(HAVE_GOSAM_CONTRIB,[test "x$have_gosam_contrib" = "xlib" -o "x$have_gosam_contrib" = "xlib64"])
if test "x$have_gosam_contrib" = "xlib" -o "x$have_gosam_contrib" = "xlib64" ; then
LOAD_GOSAM_CONTRIB="library"
CREATE_GOSAM_CONTRIB="create"
INSERT_GOSAM_CONTRIB="insert"
else
LOAD_GOSAM_CONTRIB="# library"
CREATE_GOSAM_CONTRIB="# create"
INSERT_GOSAM_CONTRIB="# insert"
fi
AC_SUBST([LOAD_GOSAM_CONTRIB])
AC_SUBST([CREATE_GOSAM_CONTRIB])
AC_SUBST([INSERT_GOSAM_CONTRIB])
])
dnl ##### OpenLoops #####
AC_DEFUN([HERWIG_CHECK_OPENLOOPS],
[
AC_MSG_CHECKING([for OpenLoops])
AC_ARG_WITH([openloops],
AS_HELP_STRING([--with-openloops=DIR], [Installation path of OpenLoops]),
[],
[with_openloops=no]
)
AC_MSG_RESULT([$with_openloops])
AS_IF([test "x$with_openloops" != "xno"],
[AC_CHECK_FILES(
${with_openloops}/lib/libopenloops.so,
[have_openloops=lib], [have_openloops=no])],
[have_openloops=no])
AS_IF([test "x$with_openloops" != "xno" -a "x$have_openloops" = "xno" ],
[AC_CHECK_FILES(
${with_openloops}/lib/libopenloops.dylib,
[have_openloops=lib], [have_openloops=no])])
AS_IF([test "x$with_openloops" != "xno" -a "x$have_openloops" = "xno" ],
[AC_CHECK_FILES(
${with_openloops}/lib64/libopenloops.so,
[have_openloops=lib64], [have_openloops=no])])
AS_IF([test "x$with_openloops" != "xno" -a "x$have_openloops" = "xno" ],
[AC_CHECK_FILES(
${with_openloops}/lib64/libopenloops.dylib,
[have_openloops=lib64], [have_openloops=no])])
AS_IF([test "x$have_openloops" = "xlib"],
[OPENLOOPSLIBS=${with_openloops}/lib
AC_SUBST(OPENLOOPSLIBS)
])
AS_IF([test "x$have_openloops" = "xlib64"],
[OPENLOOPSLIBS=${with_openloops}/lib64
AC_SUBST(OPENLOOPSLIBS)
])
AS_IF([test "x$with_openloops" != "xno" -a "x$have_openloops" = "xno"],
[AC_MSG_ERROR([OpenLoops requested but not found])])
AM_CONDITIONAL(HAVE_OPENLOOPS,[test "x$have_openloops" = "xlib" -o "x$have_openloops" = "xlib64"])
if test "x$have_openloops" = "xlib" -o "x$have_openloops" = "xlib64" ; then
OPENLOOPSPREFIX=${with_openloops}
LOAD_OPENLOOPS="library"
CREATE_OPENLOOPS="create"
INSERT_OPENLOOPS="insert"
SET_OPENLOOPS="set"
DO_OPENLOOPS="do"
MKDIR_OPENLOOPS="mkdir"
else
LOAD_OPENLOOPS="# library"
CREATE_OPENLOOPS="# create"
INSERT_OPENLOOPS="# insert"
SET_OPENLOOPS="# set"
DO_OPENLOOPS="# do"
MKDIR_OPENLOOPS="# mkdir"
fi
AC_SUBST([OPENLOOPSPREFIX])
AC_SUBST([LOAD_OPENLOOPS])
AC_SUBST([CREATE_OPENLOOPS])
AC_SUBST([INSERT_OPENLOOPS])
AC_SUBST([SET_OPENLOOPS])
AC_SUBST([DO_OPENLOOPS])
AC_SUBST([MKDIR_OPENLOOPS])
])
#########################################
dnl ##### madgraph #####
AC_DEFUN([HERWIG_CHECK_MADGRAPH],
[
AC_MSG_CHECKING([for MadGraph])
AC_ARG_WITH([madgraph],
AS_HELP_STRING([--with-madgraph=DIR], [Installation path of MadGraph]),
[],
[with_madgraph=no]
)
AC_MSG_RESULT([$with_madgraph])
AS_IF([test "x$with_madgraph" != "xno"],
[AC_CHECK_FILES(
${with_madgraph}/bin/mg5_aMC,
[have_madgraph=yes], [have_madgraph=no])],
[have_madgraph=no])
AS_IF([test "x$have_madgraph" = "xyes"],
[MADGRAPHPREFIX=${with_madgraph}
AC_SUBST(MADGRAPHPREFIX)
])
AS_IF([test "x$with_madgraph" != "xno" -a "x$have_madgraph" = "xno"],
[AC_MSG_ERROR([MadGraph requested but not found])])
AM_CONDITIONAL(HAVE_MADGRAPH,[test "x$have_madgraph" = "xyes" ])
if test "x$have_madgraph" = "xyes" ; then
LOAD_MADGRAPH="library"
CREATE_MADGRAPH="create"
INSERT_MADGRAPH="insert"
SET_MADGRAPH="set"
DO_MADGRAPH="do"
else
LOAD_MADGRAPH="# library"
CREATE_MADGRAPH="# create"
INSERT_MADGRAPH="# insert"
SET_MADGRAPH="# set"
DO_MADGRAPH="# do"
fi
AC_SUBST([LOAD_MADGRAPH])
AC_SUBST([CREATE_MADGRAPH])
AC_SUBST([INSERT_MADGRAPH])
AC_SUBST([SET_MADGRAPH])
AC_SUBST([DO_MADGRAPH])
])
AC_DEFUN([HERWIG_CHECK_PYTHIA],
[
dnl check if a directory is specified for Pythia
AC_ARG_WITH(pythia,
[AC_HELP_STRING([--with-pythia=dir],
[Assume the given directory for Pythia])])
dnl search for the pythia-config script
if test "$with_pythia" = ""; then
AC_PATH_PROG(pythiaconfig, pythia8-config, no)
else
AC_PATH_PROG(pythiaconfig, pythia8-config, no, ${with_pythia}/bin)
fi
if test "${pythiaconfig}" = "no"; then
AC_MSG_CHECKING(Pythia)
AC_MSG_RESULT(no);
PYTHIA8LIB=
# $2
else
PYTHIA8DATA=`${pythiaconfig} --datadir`/xmldoc
PYTHIA8LIB="-L`${pythiaconfig} --libdir` -lpythia8"
fi
AC_SUBST(PYTHIA8DATA)
AC_SUBST(PYTHIA8LIB)
])
dnl CHECK PYTHIA END
dnl ##### EvtGen #####
AC_DEFUN([HERWIG_CHECK_EVTGEN],
[
AC_MSG_CHECKING([for evtgen])
AC_ARG_WITH([evtgen],
AS_HELP_STRING([--with-evtgen=DIR], [Installation path of EvtGen]),
[],
[with_evtgen=no]
)
AC_MSG_RESULT([$with_evtgen])
AS_IF([test "x$with_evtgen" != "xno"],
[AC_CHECK_FILES(
${with_evtgen}/lib/libEvtGenExternal.so,
[have_evtgen=lib], [have_evtgen=no])],
[have_evtgen=no])
AS_IF([test "x$with_evtgen" != "xno" -a "x$have_evtgen" = "xno"],
[AC_CHECK_FILES(
${with_evtgen}/lib64/libEvtGenExternal.so,
[have_evtgen=lib64], [have_evtgen=no])])
AS_IF([test "x$with_evtgen" != "xno" -a "x$have_evtgen" = "xno" ],
[AC_CHECK_FILES(
${with_evtgen}/lib/libEvtGenExternal.dylib,
[have_evtgen=lib], [have_evtgen=no])])
AS_IF([test "x$have_evtgen" = "xlib" -o "x$have_evtgen" = "xlib64" ],
[EVTGENPREFIX=${with_evtgen}
AC_SUBST(EVTGENPREFIX)
])
AS_IF([test "x$with_evtgen" != "xno" -a "x$have_evtgen" = "xno"],
[AC_MSG_ERROR([EvtGen requested but not found])])
AC_SUBST([EVTGENINCLUDE],[-I$EVTGENPREFIX/include])
AM_CONDITIONAL(HAVE_EVTGEN,[test "x$have_evtgen" = "xlib" -o "x$have_evtgen" = "xlib64"])
if test "x$have_evtgen" = "xlib" ; then
LOAD_EVTGEN_DECAYS="read EvtGenBDecays.in"
LOAD_EVTGEN_DECAYER="read EvtGenDecayer.in"
EVTGENLIBS="-L$with_evtgen/lib -lEvtGen -lEvtGenExternal"
elif test "x$have_evtgen" = "xlib64" ; then
LOAD_EVTGEN_DECAYS="read EvtGenBDecays.in"
LOAD_EVTGEN_DECAYER="read EvtGenDecayer.in"
EVTGENLIBS="-L$with_evtgen/lib64 -lEvtGen -lEvtGenExternal"
else
LOAD_EVTGEN_DECAYS="read HerwigBDecays.in"
LOAD_EVTGEN_DECAYER="#read EvtGenDecayer.in"
EVTGENLIBS=""
fi
AC_SUBST([LOAD_EVTGEN_DECAYS])
AC_SUBST([LOAD_EVTGEN_DECAYER])
AC_SUBST([EVTGENLIBS])
])
dnl ###### GSL ######
AC_DEFUN([HERWIG_CHECK_GSL],
[
AC_MSG_CHECKING([for gsl location])
GSLINCLUDE=""
GSLLIBS=""
AC_ARG_WITH(gsl,
AC_HELP_STRING([--with-gsl=DIR],[location of gsl installation @<:@default=system libs@:>@]),
[],
[with_gsl=system])
if test "x$with_gsl" = "xno"; then
AC_MSG_ERROR([libgsl is required. Please install the GNU scientific library and header files.])
fi
if test "x$with_gsl" = "xsystem"; then
AC_MSG_RESULT([in system libraries])
oldlibs="$LIBS"
AC_CHECK_LIB(m,main)
AC_CHECK_LIB(gslcblas,main)
AC_CHECK_LIB(gsl,main,[],
[
AC_MSG_ERROR([Cannot find libgsl. Please install the GNU scientific library and header files or use --with-gsl=.])
]
)
GSLLIBS="$LIBS"
GSLPATH="$with_gsl"
LIBS=$oldlibs
else
if test "`uname -m`" = "x86_64" -a -e "$with_gsl/lib64/libgsl.a" -a -d "$with_gsl/include/gsl"; then
AC_MSG_RESULT([found in $with_gsl])
GSLLIBS="-L$with_gsl/lib64 -R$with_gsl/lib64 -lgslcblas -lgsl"
GSLINCLUDE="-I$with_gsl/include"
GSLPATH="$with_gsl"
elif test -e "$with_gsl/lib/libgsl.a" -a -d "$with_gsl/include/gsl"; then
AC_MSG_RESULT([found in $with_gsl])
GSLLIBS="-L$with_gsl/lib -R$with_gsl/lib -lgslcblas -lgsl"
GSLINCLUDE="-I$with_gsl/include"
GSLPATH="$with_gsl"
else
AC_MSG_RESULT([not found])
AC_MSG_ERROR([Can't find $with_gsl/lib/libgsl.a or the headers in $with_gsl/include])
fi
fi
AC_SUBST(GSLINCLUDE)
AC_SUBST(GSLLIBS)
AC_SUBST(GSLPATH)
])
dnl ##### COMPILERFLAGS #####
AC_DEFUN([HERWIG_COMPILERFLAGS],
[
AC_REQUIRE([HERWIG_CHECK_GSL])
AC_REQUIRE([HERWIG_CHECK_THEPEG])
AC_REQUIRE([BOOST_REQUIRE])
AC_REQUIRE([AX_COMPILER_VENDOR])
AM_CPPFLAGS="-I\$(top_builddir)/include $THEPEGINCLUDE \$(GSLINCLUDE) \$(BOOST_CPPFLAGS)"
AC_MSG_CHECKING([for debugging mode])
AC_ARG_ENABLE(debug,
AC_HELP_STRING([--enable-debug],[debug mode, use --enable-debug=slow for additional options that slow down the run.]),
[],
[enable_debug=no]
)
AC_MSG_RESULT([$enable_debug])
if test "x$enable_debug" = "xno"; then
debugflags=""
else
debugflags="-g"
fi
dnl -Wfloat-equal -fvisibility-inlines-hidden -Wctor-dtor-privacy -Weffc++
if test -n "$GCC"; then
warnflags="-pedantic -Wall -Wextra -Wno-overloaded-virtual"
if test "x$enable_debug" = "xslow"; then
debugflags="$debugflags -fno-inline"
AM_CPPFLAGS="$AM_CPPFLAGS -D_GLIBCXX_DEBUG"
fi
fi
dnl do an actual capability check on ld instead of this workaround
case "${host}" in
*-darwin*)
;;
*)
AM_LDFLAGS="-Wl,--enable-new-dtags"
;;
esac
case "${ax_cv_cxx_compiler_vendor}" in
gnu)
AM_CXXFLAGS="-pedantic -Wall -W"
;;
clang)
AM_CXXFLAGS="-pedantic -Wall -Wno-overloaded-virtual -Wno-unused-function -Wno-unused-parameter"
dnl -Wno-unneeded-internal-declaration
;;
intel)
AM_CXXFLAGS="-strict-ansi -Wall -wd13000,1418,981,444,383,1599,1572,2259,980"
;;
esac
AC_SUBST(AM_CPPFLAGS)
AC_SUBST(AM_CFLAGS, ["$warnflags $debugflags"])
AC_SUBST(AM_CXXFLAGS,["$AM_CXXFLAGS $warnflags $debugflags"])
AC_SUBST(AM_FCFLAGS, ["$debugflags"])
AC_SUBST(AM_LDFLAGS)
])
AC_DEFUN([HERWIG_ENABLE_MODELS],
[
AC_MSG_CHECKING([if BSM models should be built])
AC_ARG_ENABLE(models,
AC_HELP_STRING([--disable-models],[Turn off compilation of BSM models.]),
[],
[enable_models=yes]
)
AC_MSG_RESULT([$enable_models])
LOAD_BSM=""
if test "$enable_models" = "yes"; then
LOAD_BSM="read BSMlibs.in"
fi
AC_SUBST(LOAD_BSM)
AM_CONDITIONAL(WANT_BSM,[test "$enable_models" = "yes"])
])
AC_DEFUN([HERWIG_OVERVIEW],
[
FCSTRING=`$FC --version | head -1`
CXXSTRING=`$CXX --version | head -1`
CCSTRING=`$CC --version | head -1`
if test "x$PYTHON" != "x:"
then
python_was_found="yes, using Python $PYTHON_VERSION"
else
python_was_found="no, requires Python >= 2.6"
fi
cat << _HW_EOF_ > config.herwig
*****************************************************
*** $PACKAGE_STRING configuration summary
*** Please include this information in bug reports!
***--------------------------------------------------
*** Prefix: $prefix
***
*** BSM models: $enable_models
*** UFO converter: ${python_was_found}
***
*** Herwig debug mode: $enable_debug
***
*** ThePEG: $with_thepeg
*** ThePEG headers: $with_thepeg_headers
***
*** GoSam: $with_gosam
*** GoSam-Contrib: $with_gosam_contrib
*** MadGraph: $with_madgraph
*** njet: $with_njet ${NJET_VERSION}
*** OpenLoops: $with_openloops
*** VBFNLO: $with_vbfnlo
***
*** EvtGen: $with_evtgen
*** GSL: $with_gsl
*** boost: ${BOOST_CPPFLAGS:-system}
*** Fastjet: ${fjconfig}
***
*** Host: $host
*** CC: $CCSTRING
*** CXX: $CXXSTRING
*** FC: $FCSTRING
***
*** CXXFLAGS: $CXXFLAGS
*** FCFLAGS: $FCFLAGS
*****************************************************
_HW_EOF_
])