Page MenuHomeHEPForge

No OneTemporary

diff --git a/Handlers/StandardEventHandler.cc b/Handlers/StandardEventHandler.cc
--- a/Handlers/StandardEventHandler.cc
+++ b/Handlers/StandardEventHandler.cc
@@ -1,803 +1,804 @@
// -*- C++ -*-
//
// StandardEventHandler.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2011 Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the StandardEventHandler class.
//
#include "StandardEventHandler.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include "ThePEG/Handlers/StdXCombGroup.h"
#include "ThePEG/Handlers/SubProcessHandler.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Utilities/SimplePhaseSpace.h"
#include "ThePEG/Utilities/LoopGuard.h"
#include "ThePEG/Utilities/Debug.h"
#include "ThePEG/PDF/PartonExtractor.h"
#include "ThePEG/MatrixElement/MEBase.h"
#include "ThePEG/MatrixElement/MEGroup.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Handlers/LuminosityFunction.h"
#include "ThePEG/Handlers/SamplerBase.h"
#include "ThePEG/Handlers/CascadeHandler.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/Config/algorithm.h"
#include <iomanip>
#include <sstream>
using namespace ThePEG;
StandardEventHandler::StandardEventHandler()
: EventHandler(false), collisionCuts(true), theBinStrategy(2), weightedEvents(false),
theLumiDim(0) {
setupGroups();
}
StandardEventHandler::~StandardEventHandler() {}
void StandardEventHandler::reject(double w) {
tStdXCombPtr last = dynamic_ptr_cast<tStdXCombPtr>(lastXCombPtr());
if ( !last ) return;
last->reject(w);
}
void StandardEventHandler::reweight(double factor) const {
tStdXCombPtr last = dynamic_ptr_cast<tStdXCombPtr>(lastXCombPtr());
if ( !currentEvent() || !last )
return;
double weight = currentEvent()->weight();
last->reweight(weight,factor*weight);
currentEvent()->weight(factor*weight);
for ( map<string,double>::iterator w =
currentEvent()->optionalWeights().begin();
w != currentEvent()->optionalWeights().end(); ++w )
w->second *= factor;
}
void StandardEventHandler::doupdate() {
EventHandler::doupdate();
bool redo = touched();
UpdateChecker::check(theIncomingA, redo);
UpdateChecker::check(theIncomingB, redo);
for_each(subProcesses(), UpdateChecker(redo));
if ( !redo ) return;
theIncoming.first = theIncomingA;
theIncoming.second = theIncomingB;
for ( SubHandlerList::iterator sit = subProcesses().begin();
sit != subProcesses().end(); ++sit )
if ( !(**sit).pExtractor()->canHandle(incoming()) )
throw StandardEventHandlerUpdateException()
<< "Cannot use the parton extractor '" << (**sit).pExtractor()->name()
<< "' in the SubProcessHandler '" << (**sit).name() << "' in the "
<< "StandardEventHandler '" << name() << "' since it cannot handle "
<< "the requested types of incoming particles ("
<< theIncomingA->name() << "," << theIncomingB->name() << ").";
}
void StandardEventHandler::doinit() {
EventHandler::doinit();
if ( !lumiFnPtr() ) throw StandardEventHandlerUpdateException()
<< "The StandardEventHandler '" << name() << "' does not have any "
<< "LuminosityFunction object assigned to it, which it needs to be "
<< "able to generate events." << Exception::warning;
}
IBPtr StandardEventHandler::clone() const {
return new_ptr(*this);
}
IBPtr StandardEventHandler::fullclone() const {
return new_ptr(*this);
}
void StandardEventHandler::
addME(Energy maxEnergy, tSubHdlPtr sub, tPExtrPtr extractor, tCutsPtr cuts,
tCascHdlPtr ckkw, tMEPtr me, const PBPair & pBins,
const PartonPairVec& allPBins) {
typedef MEBase::DiagramVector DiagramVector;
typedef map<string,DiagramVector> DiagramMap;
cPDPair pin(pBins.first->parton(), pBins.second->parton());
DiagramVector diag = me->diagrams();
DiagramMap tdiag;
DiagramMap tmdiag;
for ( int i = 0, N = diag.size(); i < N; ++i ) {
cPDPair din(diag[i]->partons()[0], diag[i]->partons()[1]);
if (!me->noMirror())
if ( din.first->id() < din.second->id() ) swap(din.first, din.second);
if ( din == pin ) tdiag[diag[i]->getTag()].push_back(diag[i]);
if (!me->noMirror())
if ( din.first == pin.second && din.second == pin.first )
tmdiag[diag[i]->getTag()].push_back(diag[i]);
}
if ( tdiag.empty() ) tdiag = tmdiag;
for ( DiagramMap::iterator dit = tdiag.begin(); dit != tdiag.end(); ++dit ) {
cPDPair din(dit->second.back()->partons()[0],
dit->second.back()->partons()[1]);
// check
assert(me->noMirror() ? din == pin : true);
StdXCombPtr xcomb = me->makeXComb(maxEnergy, incoming(), this, sub, extractor,
ckkw, pBins, cuts, dit->second, din != pin,
allPBins);
if ( xcomb->checkInit() ) xCombs().push_back(xcomb);
else generator()->logWarning(
StandardEventHandlerInitError() << "The matrix element '"
<< xcomb->matrixElement()->name()
<< "' cannot generate the diagram '"
<< dit->first
<< "' when used together with the parton extractor '"
<< xcomb->pExtractor()->name()
<< "'.\nThe corresponding diagram is switched off, "
<< "check the collision energy and/or the cuts."
<< Exception::warning);
}
}
void StandardEventHandler::initGroups() {
tStdXCombPtr lastXC = dynamic_ptr_cast<tStdXCombPtr>(lastXCombPtr());
if ( lastXC ) optGroups = lastXC->subProcessHandler()->groups();
EventHandler::initGroups();
}
tCollPtr StandardEventHandler::performCollision() {
tStdXCombPtr lastXC = dynamic_ptr_cast<tStdXCombPtr>(lastXCombPtr());
if ( CKKWHandler() ) CKKWHandler()->setXComb(lastXCombPtr());
lastExtractor()->select(lastXC);
currentCollision(new_ptr(Collision(lastParticles(), currentEvent(), this)));
if ( currentEvent() ) currentEvent()->addCollision(currentCollision());
currentStep(new_ptr(Step(currentCollision())));
currentCollision()->addStep(currentStep());
currentStep()->addSubProcess(lastXC->construct());
lastExtractor()->construct(lastXC->partonBinInstances(), currentStep());
if ( collisionCuts )
if ( !lastCuts().passCuts(*currentCollision()) ) throw Veto();
initGroups();
if ( ThePEG_DEBUG_ITEM(1) ) {
if ( currentEvent() )
generator()->logfile() << *currentEvent();
else
generator()->logfile() << *currentCollision();
}
return continueCollision();
}
void StandardEventHandler::setScale(Energy2 scale) {
lastXCombPtr()->lastScale(scale);
}
void StandardEventHandler::initialize() {
theLumiDim = lumiFn().nDim(incoming());
Energy maxEnergy = lumiFn().maximumCMEnergy();
xCombs().clear();
xSecs().clear();
cuts()->initialize(sqr(maxEnergy), lumiFn().Y());
for ( SubHandlerList::const_iterator sit = subProcesses().begin();
sit != subProcesses().end(); ++sit ) {
CutsPtr kincuts = (**sit).cuts()? (**sit).cuts(): cuts();
if ( (**sit).cuts() ) kincuts->initialize(sqr(maxEnergy), lumiFn().Y());
PExtrPtr pextract = (**sit).pExtractor();
tCascHdlPtr ckkw = (**sit).CKKWHandler();
if ( !ckkw ) ckkw = CKKWHandler();
PartonPairVec vpc = pextract->getPartons(maxEnergy, incoming(), *kincuts);
// The last parton bin pair was in fact the bins corresponding to
// the incoming particles, so we remove them, but save them to
// keep them referenced.
PBPair orig = vpc.back();
vpc.pop_back();
for ( PartonPairVec::iterator ppit = vpc.begin();
ppit != vpc.end(); ++ppit )
for ( MEVector::const_iterator meit = (**sit).MEs().begin();
meit != (**sit).MEs().end(); ++meit ) {
addME(maxEnergy, *sit, pextract, kincuts, ckkw, *meit, *ppit,vpc);
}
}
xSecs().resize(xCombs().size());
theMaxDims.clear();
switch ( binStrategy() ) {
case 0: {
theMaxDims.push_back(0);
for ( int i = 0, N = xCombs().size(); i < N; ++i )
theMaxDims[0] = max(theMaxDims[0], xCombs()[i]->nDim());
break;
}
case 1: {
for ( int i = 0, N = xCombs().size(); i < N; ++i )
theMEXMap[xCombs()[i]->matrixElement()].push_back(xCombs()[i]);
MEXMap::const_iterator mei = theMEXMap.begin();
for ( int i = 0, N = theMEXMap.size(); i < N; ++i, ++mei) {
theMaxDims.push_back(0);
for ( int j = 0, M = mei->second.size(); j < M; ++j )
theMaxDims[i] = max(theMaxDims[i], mei->second[j]->nDim());
}
break;
}
case 2: {
for ( int i = 0, N = xCombs().size(); i < N; ++i )
theMaxDims.push_back(xCombs()[i]->nDim());
break;
}
}
sampler()->setEventHandler(this);
sampler()->initialize();
}
CrossSection StandardEventHandler::
dSigDR(const pair<double,double> ll, Energy2 maxS,
int ibin, int nr, const double * r) {
PPair inc = make_pair(incoming().first->produceParticle(),
incoming().second->produceParticle());
SimplePhaseSpace::CMS(inc, maxS);
XVector xv;
switch ( binStrategy() ) {
case 0:
xv = xCombs();
break;
case 1: {
MEXMap::iterator mei = theMEXMap.begin();
for ( int i = 0; i < ibin; ++i) ++mei;
xv = mei->second;
break;
}
case 2:
xv = XVector(1, xCombs()[ibin]);
break;
}
xSecs().resize(xv.size());
for ( int i = 0, N = xv.size(); i < N; ++i ) xv[i]->prepare(inc);
CrossSection sum = ZERO;
for ( int i = 0, N = xv.size(); i < N; ++i )
xSecs()[i] = ( sum += xv[i]->dSigDR(ll, nr, r) );
return sum;
}
tStdXCombPtr StandardEventHandler::select(int bin, double & weight) {
int i = 0;
if ( binStrategy() != 2 )
i = upper_bound(xSecs().begin(), xSecs().end(), rnd()*xSecs().back())
- xSecs().begin();
tStdXCombPtr lastXC;
switch ( binStrategy() ) {
case 0:
lastXC = xCombs()[i];
break;
case 1: {
MEXMap::iterator mei = theMEXMap.begin();
for ( int j = 0; j < bin; ++j) ++mei;
lastXC = mei->second[i];
break;
}
case 2:
lastXC = xCombs()[bin];
break;
}
// clean up the old XComb object before switching to a new one
if ( theLastXComb && theLastXComb != lastXC ) theLastXComb->clean();
theLastXComb = lastXC;
lastXC->matrixElement()->setXComb(lastXC);
weight /= lastXC->matrixElement()->preWeight();
lastXC->select(weight);
lastXC->accept();
return lastXC;
}
int StandardEventHandler::nBins() const {
switch ( binStrategy() ) {
case 0: return 1;
case 1: return theMEXMap.size();
case 2: return xCombs().size();
}
return -1;
}
struct Stat {
Stat() : attempted(0), accepted(0), sumw(0.0), sumw2(),
maxXSec(CrossSection()), totsum(0.0), totrerr(0.0) {}
Stat(long att, long acc, double w, double w2, CrossSection x,
double sumw, double genrerr)
: attempted(att), accepted(acc), sumw(w), sumw2(w2), maxXSec(x),
totsum(sumw), totrerr(genrerr) {}
inline CrossSection xSec() const {
return totsum != 0.0? maxXSec*sumw/totsum: maxXSec;
}
inline CrossSection xSecErr() const {
if ( totsum == 0.0 ) return maxXSec;
if ( sumw == 0.0 ) return xSec();
return abs(xSec())*sqrt(sqr(totrerr) + sumw2/sqr(sumw));
}
long attempted;
long accepted;
double sumw;
double sumw2;
CrossSection maxXSec;
double totsum;
double totrerr;
const Stat & operator+= (const Stat & s) {
attempted += s.attempted;
accepted += s.accepted;
sumw += s.sumw;
sumw2 += s.sumw2;
// replaced by taking whatever is in the one to be added, see
// comments below; SP 2013/03/14
/*
totsum = max(totsum, s.totsum);
if ( totsum != 0.0 )
maxXSec = max(maxXSec, s.maxXSec);
else
maxXSec += s.maxXSec;
*/
// whereever Stat is used nothing happens with the above code,
// because those quantities are always the same for all Stat
// objects involved, so long as cross sections and sum of weights
// are positive; but it fails, if we have negative cross sections
// (which can happen when running real-subtraction of a NLO
// without the Born and virtual contributions). To this extent
// it's commented out and replaced by just taking whatever is in
// the statistics to be added. There are more issues in the
// statistics treatment going on here which urgently need to be
// looked at, and the Stat object itself needs to be made much
// more transparent and cleaned up from misleading constructions
// like this and the old one below the next two lines.
totsum = s.totsum;
maxXSec = s.maxXSec;
totrerr = s.totrerr;
return *this;
}
};
void StandardEventHandler::statistics(ostream & os) const {
if ( statLevel() == 0 ) return;
map<cPDPair, Stat> partonMap;
map<MEPtr, Stat> meMap;
map<PExtrPtr, Stat> extractMap;
Stat tot;
double genrerr = sampler()->integratedXSecErr()/sampler()->integratedXSec();
for ( int i = 0, N = xCombs().size(); i < N; ++i ) {
const StandardXComb & x = *xCombs()[i];
Stat s(x.stats().attempts(), x.stats().accepted(),
x.stats().sumWeights(), x.stats().sumWeights2(),
sampler()->integratedXSec(), sampler()->sumWeights(), genrerr);
partonMap[x.partons()] += s;
meMap[x.matrixElement()] += s;
extractMap[x.pExtractor()] += s;
tot += s;
}
string line = string(78, '=') + "\n";
if ( tot.accepted <= 0 ) {
os << line << "No events generated by event handler '" << name() << "'."
<< endl;
return;
}
os << line << "Statistics for event handler \'" << name() << "\':\n"
<< " "
<< "generated number of Cross-section\n"
<< " "
<< " events attempts (nb)\n";
os << line << "Total (from weighted events): including vetoed events" << setw(23)
<< ouniterr(sampler()->integratedXSec(),
sampler()->integratedXSecErr(), nanobarn)
<< endl;
- os << line << "Total (from unweighted events):" << setw(17) << tot.accepted << setw(13)
+ os << line << "Total (from "
+ << (weighted() ? " weighted" : "unweighted") << " events):"
+ << setw(17) << tot.accepted << setw(13)
<< tot.attempted << setw(17)
<< ouniterr(tot.xSec(),tot.xSecErr() , nanobarn)
<< endl << line;
if ( statLevel() == 1 ) return;
os << "Per matrix element breakdown:\n";
for ( map<MEPtr, Stat>::iterator i = meMap.begin();
i != meMap.end(); ++i ) {
string n = i->first->name();
n.resize(37, ' ');
os << n << setw(11) << i->second.accepted << setw(13)
<< i->second.attempted << setw(17)
<< ouniterr(i->second.xSec(), i->second.xSecErr(), nanobarn)
<< endl;
}
os << line;
if ( statLevel() == 2 ) return;
os << "Per parton extractor breakdown:\n";
for ( map<PExtrPtr, Stat>::iterator i = extractMap.begin();
i != extractMap.end(); ++i ) {
string n = i->first->name();
n.resize(37, ' ');
os << n << setw(11) << i->second.accepted << setw(13)
<< i->second.attempted << setw(17)
<< ouniterr(i->second.xSec(), i->second.xSecErr(), nanobarn)
<< endl;
}
os << line;
os << "Per incoming partons breakdown:\n";
for ( map<cPDPair, Stat>::iterator i = partonMap.begin();
i != partonMap.end(); ++i ) {
string n = i->first.first->PDGName() + " " + i->first.second->PDGName();
n.resize(37, ' ');
os << n << setw(11) << i->second.accepted << setw(13)
<< i->second.attempted << setw(17)
<< ouniterr(i->second.xSec(), i->second.xSecErr(), nanobarn)
<< endl;
}
os << line;
if ( statLevel() == 3 ) return;
os << "Detailed breakdown:\n";
CrossSection xsectot = sampler()->integratedXSec()/sampler()->sumWeights();
for ( int i = 0, N = xCombs().size(); i < N; ++i ) {
const StandardXComb & x = *xCombs()[i];
os << "(" << x.pExtractor()->name() << ") "
<< x.partons().first->PDGName() << " "
<< x.partons().second->PDGName()
<< " (" << x.matrixElement()->name() << " "
<< x.lastDiagram()->getTag() << ") " << endl
<< setw(48) << x.stats().accepted() << setw(13) << x.stats().attempts()
<< setw(17)
<< ouniterr(x.stats().sumWeights()*xsectot,
sqrt(x.stats().sumWeights2())*xsectot, nanobarn) << endl;
}
os << line;
}
CrossSection StandardEventHandler::histogramScale() const {
Stat tot;
for ( int i = 0, N = xCombs().size(); i < N; ++i ) {
const StandardXComb & x = *xCombs()[i];
Stat s;
s = Stat(x.stats().attempts(), x.stats().accepted(),
x.stats().sumWeights(), x.stats().sumWeights2(),
sampler()->integratedXSec(), sampler()->sumWeights(), 1.0);
tot += s;
}
return tot.xSec()/tot.sumw;
}
CrossSection StandardEventHandler::integratedXSec() const {
if ( sampler()->integratedXSec() == ZERO )
return sampler()->maxXSec();
Stat tot;
for ( int i = 0, N = xCombs().size(); i < N; ++i ) {
const StandardXComb & x = *xCombs()[i];
Stat s;
s = Stat(x.stats().attempts(), x.stats().accepted(),
x.stats().sumWeights(), x.stats().sumWeights2(),
sampler()->integratedXSec(), sampler()->sumWeights(), 1.0);
tot += s;
}
return tot.xSec();
-
}
CrossSection StandardEventHandler::integratedXSecErr() const {
if ( sampler()->integratedXSec() == ZERO )
return sampler()->maxXSec();
Stat tot;
double genrerr = sampler()->integratedXSecErr()/sampler()->integratedXSec();
for ( int i = 0, N = xCombs().size(); i < N; ++i ) {
const StandardXComb & x = *xCombs()[i];
Stat s;
s = Stat(x.stats().attempts(), x.stats().accepted(),
x.stats().sumWeights(), x.stats().sumWeights2(),
sampler()->integratedXSec(), sampler()->sumWeights(), genrerr);
tot += s;
}
return tot.xSecErr();
}
CrossSection StandardEventHandler::integratedXSecNoReweight() const {
if ( sampler()->integratedXSec() == ZERO )
return sampler()->maxXSec();
Stat tot;
for ( int i = 0, N = xCombs().size(); i < N; ++i ) {
const StandardXComb & x = *xCombs()[i];
Stat s;
s = Stat(x.stats().attempts(), x.stats().accepted(),
x.stats().sumWeightsNoReweight(), x.stats().sumWeights2NoReweight(),
sampler()->integratedXSec(), sampler()->sumWeights(), 1.0);
tot += s;
}
return tot.xSec();
}
CrossSection StandardEventHandler::integratedXSecErrNoReweight() const {
if ( sampler()->integratedXSec() == ZERO )
return sampler()->maxXSec();
Stat tot;
double genrerr = sampler()->integratedXSecErr()/sampler()->integratedXSec();
for ( int i = 0, N = xCombs().size(); i < N; ++i ) {
const StandardXComb & x = *xCombs()[i];
Stat s;
s = Stat(x.stats().attempts(), x.stats().accepted(),
x.stats().sumWeightsNoReweight(), x.stats().sumWeights2NoReweight(),
sampler()->integratedXSec(), sampler()->sumWeights(), genrerr);
tot += s;
}
return tot.xSecErr();
}
void StandardEventHandler::doinitrun() {
EventHandler::doinitrun();
for ( SubHandlerList::iterator sit = subProcesses().begin();
sit != subProcesses().end(); ++sit )
(**sit).initrun();
sampler()->initrun();
for ( int i = 0, N = xCombs().size(); i < N; ++i )
xCombs()[i]->reset();
}
CrossSection StandardEventHandler::dSigDR(const vector<double> & r) {
double jac = 1.0;
pair<double,double> ll = lumiFn().generateLL(&r[0], jac);
Energy2 maxS = sqr(lumiFn().maximumCMEnergy())/exp(ll.first + ll.second);
int bin = sampler()->lastBin();
CrossSection x = jac*lumiFn().value(incoming(), ll.first, ll.second)
*dSigDR(ll, maxS, bin, nDim(bin) - lumiDim(), &r[lumiDim()]);
return x;
}
EventPtr StandardEventHandler::generateEvent() {
LoopGuard<EventLoopException,StandardEventHandler>
loopGuard(*this, maxLoop());
while (1) {
loopGuard();
EventHandler::clean();
double weight = sampler()->generate();
tStdXCombPtr lastXC = select(sampler()->lastBin(), weight);
try {
lumiFn().select(lastXC);
currentEventBoost() = lumiFn().getBoost();
currentEvent(new_ptr(Event(lastParticles(), this, generator()->runName(),
generator()->currentEventNumber(), weight)));
performCollision();
if ( !currentCollision() ) throw Veto();
currentEvent()->transform(currentEventBoost());
return currentEvent();
}
catch (Veto) {
reject(currentEvent()->weight());
}
catch (Stop) {
break;
}
catch (Exception &) {
reject(currentEvent()->weight());
throw;
}
}
return currentEvent();
}
EventPtr StandardEventHandler::continueEvent() {
if ( !generator() ) throw StandardEventHandlerInitError()
<< "The event handler '" << name() << "' had not been isolated "
<< "from the setup phase before it was used." << Exception::maybeabort;
try {
continueCollision();
}
catch (Veto) {
reject(currentEvent()->weight());
}
catch (Stop) {
}
catch (Exception &) {
reject(currentEvent()->weight());
throw;
}
return currentEvent();
}
void StandardEventHandler::select(tXCombPtr newXComb) {
EventHandler::select(newXComb);
lastExtractor()->select(newXComb);
}
void StandardEventHandler::clean() {
if ( theLastXComb ) theLastXComb->clean();
for (size_t i=0; i < theXCombs.size(); ++i )
if ( theXCombs[i] ) theXCombs[i]->clean();
EventHandler::clean();
}
void StandardEventHandler::dofinish() {
clean();
EventHandler::dofinish();
}
ClassDescription<StandardEventHandler>
StandardEventHandler::initStandardEventHandler; //
void StandardEventHandler::setIncomingA(PDPtr p) {
theIncomingA = p;
theIncoming.first = p;
}
void StandardEventHandler::setIncomingB(PDPtr p) {
theIncomingB = p;
theIncoming.second = p;
}
void StandardEventHandler::Init() {
static ClassDocumentation<StandardEventHandler> documentation
("This is the standard event handler to generate hard sub-processes "
"within ThePEG. It must specify a pair of incoming particle beams "
"in <interface>BeamA</interface> and <interface>BeamB</interface> "
"and a suiteable <interface>LuminosityFunction</interface>. In "
"addition at least one object describing the sub-processes to be "
"generated must be specified in "
"<interface>SubProcessHandlers</interface>.");
static Reference<StandardEventHandler,ParticleData> interfaceIncomingA
("BeamA",
"The type of particles in first beam",
&StandardEventHandler::theIncomingA, false, false, true, false,
&StandardEventHandler::setIncomingA);
static Reference<StandardEventHandler,ParticleData> interfaceIncomingB
("BeamB",
"The type of particles in second beam",
&StandardEventHandler::theIncomingB, false, false, true, false,
&StandardEventHandler::setIncomingB);
static RefVector<StandardEventHandler,SubProcessHandler> interfaceSubhandlers
("SubProcessHandlers",
"The list of sub-process handlers used in this StandardEventHandler. ",
&StandardEventHandler::theSubProcesses, 0, false, false, true, false);
static Reference<StandardEventHandler,Cuts> interfaceCuts
("Cuts",
"Common kinematical cuts for this StandardEventHandler. These cuts "
"may be overidden in individual sub-process handlers.",
&StandardEventHandler::theCuts, false, false, true, false);
static Switch<StandardEventHandler,bool> interfaceCollisionCuts
("CollisionCuts",
"Switch on or off cuts on collision objects",
&StandardEventHandler::collisionCuts, true, false, false);
static SwitchOption interfaceCollisionCutsOn
(interfaceCollisionCuts,
"On",
"Switch on cuts on collision objects",
true);
static SwitchOption interfaceCollisionCutsOff
(interfaceCollisionCuts,
"Off",
"Switch off cuts on collision cuts",
false);
static Switch<StandardEventHandler,int> interfaceBinStrategy
("BinStrategy",
"The strategy to be used when sampling different ThePEG::XComb "
"objects. An ThePEG::XComb objet represents a pair of incoming "
"parton types as defined by a THePEG::PartonExtractor and a "
"matrix element.",
&StandardEventHandler::theBinStrategy, 2, false, false);
static SwitchOption interfaceBinStrategy0
(interfaceBinStrategy,
"AllAtOnce",
"All bins are sampled together.",
0);
static SwitchOption interfaceBinStrategy1
(interfaceBinStrategy,
"PerME",
"All bins which have the same matrix element object are sampled together.",
1);
static SwitchOption interfaceBinStrategy2
(interfaceBinStrategy,
"Individual",
"All bins are sampled individually.",
2);
static Switch<StandardEventHandler,bool> interfaceWeighted
("Weighted",
"If switched on, this event Handler will produce weighted events",
&StandardEventHandler::weightedEvents, false, false);
static SwitchOption interfaceWeightedTrue
(interfaceWeighted, "On",
"This EventHandler produces weighted events.", true);
static SwitchOption interfaceWeightedFalse
(interfaceWeighted, "Off",
"This EventHandler produces unweighted events.", false);
static Reference<StandardEventHandler,SamplerBase> interfaceSampler
("Sampler",
"The phase space sampler responsible for generating phase space"
"points according to the cross section given by this event handler",
&StandardEventHandler::theSampler, false, false, true, true);
interfaceSubhandlers.rank(11);
interfaceIncomingA.rank(3);
interfaceIncomingB.rank(2);
}
void StandardEventHandler::persistentOutput(PersistentOStream & os) const {
os << theIncomingA << theIncomingB << theSubProcesses << theCuts << collisionCuts
<< theXCombs << ounit(theXSecs, nanobarn)
<< theBinStrategy << theMaxDims << theMEXMap
<< weightedEvents
<< theSampler << theLumiDim;
}
void StandardEventHandler::persistentInput(PersistentIStream & is, int) {
is >> theIncomingA >> theIncomingB >> theSubProcesses >> theCuts >> collisionCuts
>> theXCombs >> iunit(theXSecs, nanobarn)
>> theBinStrategy >> theMaxDims>> theMEXMap
>> weightedEvents
>> theSampler >> theLumiDim;
}
diff --git a/Utilities/UnitIO.h b/Utilities/UnitIO.h
--- a/Utilities/UnitIO.h
+++ b/Utilities/UnitIO.h
@@ -1,297 +1,305 @@
// -*- C++ -*-
//
// UnitIO.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2011 Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_UnitIO_H
#define ThePEG_UnitIO_H
// This is the declaration of the IUnit and OUnit classes and
// associated templated functions.
#include <complex>
#include <iomanip>
#include <sstream>
#include <cstdlib>
#include <cmath>
// Workarounds for OS X
#if defined __APPLE__ && defined __MACH__
extern "C" int isnan(double) throw();
extern "C" int isinf(double) throw();
#endif
namespace ThePEG {
/**
* The OUnit< class is used to
* facilitate output of unitful numbers to a
* persistent stream. An Energy can hence be written like
* this:<BR> <code>os
* << ounit(x, GeV);</code><BR> Also containers of unitful
* numbers may be written like this, as well as LorentzVector and
* ThreeVector.
*
* @see PersistentOStream
* @see PersistentIStream
*
*/
template <typename T, typename UT>
struct OUnit {
/** Constructor given an object to be written assuming the given
* unit. */
OUnit(const T & t, const UT & u): theX(t), theUnit(u) {}
/** Copy constructor */
OUnit(const OUnit<T,UT> & iu): theX(iu.theX), theUnit(iu.theUnit) {}
/** Reference to the object to be written. */
const T & theX;
/** The unit assumed when writing the object. */
const UT & theUnit;
};
/**
* The IUnit class is used to facilitate input of unitful numbers from
* and to a persistent stream. An Energy can hence be read like
* this:<BR> <code>is >> iunit(x, GeV);</code><BR> Also containers of
* unitful numbers may be read like this, as well as LorentzVector and
* ThreeVector.
*
* @see PersistentOStream
* @see PersistentIStream
*
*/
template <typename T, typename UT>
struct IUnit {
/** Constructor given an object to be read assuming the given
* unit. */
IUnit(T & t, const UT & u): theX(t), theUnit(u) {}
/** Copy constructor */
IUnit(const IUnit<T,UT> & iu): theX(iu.theX), theUnit(iu.theUnit) {}
/** Reference to the object to be read. */
T & theX;
/** The unit assumed when reading the object. */
const UT & theUnit;
};
/** Helper function creating a OUnit object given an object and a
* unit. */
template <typename T, typename UT>
inline OUnit<T,UT> ounit(const T & t, const UT & ut) {
return OUnit<T,UT>(t, ut);
}
/** Helper function creating a IUnit object given an object and a
* unit. */
template <typename T, typename UT>
inline IUnit<T,UT> iunit(T & t, const UT & ut) {
return IUnit<T,UT>(t, ut);
}
/** Helper function writing out an object with a given unit to an
* output stream. */
template <typename OStream, typename T, typename UT>
void ounitstream(OStream & os, const T & t, UT & u) {
os << t/u;
}
/** Helper function reading an object with a given unit from an
* input stream. */
template <typename IStream, typename T, typename UT>
void iunitstream(IStream & is, T & t, UT & u) {
double d;
is >> d;
t = d*u;;
}
/** Helper function reading a complex object with a given unit from an
* input stream. */
template <typename IStream, typename T, typename UT>
void iunitstream(IStream & is, std::complex<T> & t, UT & u) {
std::complex<double> d;
is >> d;
t = d*u;;
}
/** Output an OUnit object to a stream. */
template <typename OStream, typename T, typename UT>
OStream & operator<<(OStream & os, const OUnit<T,UT> & u) {
ounitstream(os, u.theX, u.theUnit);
return os;
}
/** Input an IUnit object from a stream. */
template <typename IStream, typename T, typename UT>
IStream & operator>>(IStream & is, const IUnit<T,UT> & u) {
iunitstream(is, u.theX, u.theUnit);
return is;
}
/**
* OUnitErr is used to write out unitful numbers with an error
* estimate on a standard ostream. using the helper function ouniterr
* an energy <code>e</code> with an error estimate <code>de</code> can
* be written out as eg. <code>cout << ouniterr(e, de,
* GeV);</code>. The result will be presented in scientific format
* (with the exponent divisible by three) with the relevant number of
* significant digits with a single digit in parenthesis indicating
* the error in the least significant digit,
* eg. <code>1.23(2)e+03</code>.
*/
template <typename T, typename UT>
struct OUnitErr {
/** Constructor given an object to be written assuming the given
* unit. */
OUnitErr(const T & t, const T & dt, const UT & u): x(t/u), dx(dt/u) {}
/** The number to be written. */
double x;
/** The estimated error of the number to be written. */
double dx;
};
/** Helper function creating a OUnitErr object. */
template <typename T, typename UT>
inline OUnitErr<T,UT> ouniterr(const T & t, const T & dt, const UT & ut) {
return OUnitErr<T,UT>(t, dt, ut);
}
/** Helper function creating a OUnitErr object. */
inline OUnitErr<double,double> ouniterr(double t, double dt) {
return OUnitErr<double,double>(t, dt, 1.0);
}
/** Output an OUnitErr object to a stream. */
template <typename OStream, typename T, typename UT>
OStream & operator<<(OStream & os, const OUnitErr<T,UT> & u) {
if ( isnan(u.x) || isinf(u.x) ) return os << u.x;
- if ( isnan(u.dx) || isinf(u.dx) ) return os << u.x << '(' << u.dx << ')';
+ if ( isnan(u.dx) || isinf(u.dx) ) {
+ ostringstream out;
+ out << u.x << '(' << u.dx << ')';
+ return os << out.str();
+ }
double dx = min(u.dx, abs(u.x));
if ( dx <= 0.0 ) return os << u.x;
+ double x = abs(u.x);
ostringstream osse;
osse << std::scientific << setprecision(0) << dx;
string sse = osse.str();
string::size_type ee = sse.find('e');
- long m = static_cast<long>(round(abs(u.x)/std::pow(10.0,std::atoi(sse.substr(ee + 1).c_str()))));
+ long m = static_cast<long>(round(abs(x)/std::pow(10.0,std::atoi(sse.substr(ee + 1).c_str()))));
int powx = m <= 0? os.precision(): int(log10(double(m)));
if ( m <= 0 || powx > os.precision() ) sse[0]='0';
ostringstream oss;
- oss << std::scientific << setprecision(powx) << u.x;
+ oss << std::scientific << setprecision(powx) << x;
string ss = oss.str();
string::size_type e = ss.find('e');
ostringstream out;
int pp = std::atoi(ss.substr(e + 1).c_str());
if ( pp%3 == 0 )
out << ss.substr(0, e) << "(" << sse[0] << ")" << ss.substr(e);
else if ( (pp - 1)%3 == 0 ) {
ostringstream oss;
- oss << std::scientific << setprecision(powx) << u.x/10.0;
+ oss << std::scientific << setprecision(powx) << x/10.0;
string ss = oss.str();
string::size_type e = ss.find('e');
if ( powx == 0 )
out << ss.substr(0, e) << "0(" << sse[0] << "0)" << ss.substr(e);
else if ( powx == 1 )
out << ss.substr(0, ss.find('.'))
<< ss.substr(ss.find('.') + 1, e - ss.find('.') - 1)
<< "(" << sse[0] << ")" << ss.substr(e);
else {
swap(ss[ss.find('.')], ss[ss.find('.') + 1]);
out << ss.substr(0, e) << "(" << sse[0] << ")" << ss.substr(e);
}
}
else {
ostringstream oss;
- oss << std::scientific << setprecision(powx) << u.x*10.0;
+ oss << std::scientific << setprecision(powx) << x*10.0;
string ss = oss.str();
string::size_type e = ss.find('e');
if ( powx == 0 )
out << "0." << ss.substr(0, e) << "(" << sse[0] << ")" << ss.substr(e);
else {
swap(ss[ss.find('.')], ss[ss.find('.') - 1]);
out << ss.substr(0, ss.find('.')) << "0" << ss.substr(ss.find('.'), e)
<< "(" << sse[0] << ")" << ss.substr(e);
}
}
- return os << out.str();
+ string res = out.str();
+ if ( u.x < 0.0 )
+ res = "-" + res;
+ return os << res;
}
/**
* The IUnitErr class is used to facilitate input of unitful numbers
* with error estimates written out using the OUnitErr class.
*
*/
template <typename T, typename UT>
struct IUnitErr {
/** Constructor given an object to be read assuming the given
* unit. */
IUnitErr(T & t, T & dt, const UT & u): x(t), dx(dt), ut(u) {}
/** Reference to the object to be read. */
T & x;
/** The estimated error of the number to be read. */
T & dx;
/** The unit assumed when reading the object. */
UT ut;
};
/** Helper function creating a IUnitErr object. */
template <typename T, typename UT>
inline IUnitErr<T,UT> iuniterr(T & t, T & dt, const UT & ut) {
return IUnitErr<T,UT>(t, dt, ut);
}
/** Helper function creating a OUnitErr object. */
inline IUnitErr<double,double> iuniterr(double & t, double & dt) {
return IUnitErr<double,double>(t, dt, 1.0);
}
/** Input an IUnit object from a stream. */
template <typename IStream, typename T, typename UT>
IStream & operator>>(IStream & is, const IUnitErr<T,UT> & u) {
string s;
double x = 0.0;
double dx = 0.0;
double ex = 1.0;
is >> s;
string::size_type open = s.find('(');
string::size_type close = s.find(')');
string se = "0";
string sp = "1";
double pe = 1.0;
if ( open != string::npos && close != string::npos ) {
se = s.substr(open + 1);
sp += s.substr(close + 1);
string::size_type dot = s.find('.');
if ( dot != string::npos && dot < open ) pe = std::pow(10.0, 1.0 - (open - dot));
}
istringstream(s) >> x;
istringstream(se) >> dx;
istringstream(sp) >> ex;
u.x = x*ex*u.ut;
u.dx = dx*ex*pe*u.ut;
return is;
}
}
#endif /* ThePEG_UnitIO_H */

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 8:02 PM (1 d, 5 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805956
Default Alt Text
(35 KB)

Event Timeline