Page MenuHomeHEPForge
Files F8309052

No OneTemporary
Actions

View Options

diff --git a/DipoleShower/Base/DipoleSplittingGenerator.cc b/DipoleShower/Base/DipoleSplittingGenerator.cc
--- a/DipoleShower/Base/DipoleSplittingGenerator.cc
+++ b/DipoleShower/Base/DipoleSplittingGenerator.cc
@@ -1,610 +1,611 @@
// -*- C++ -*-
//
// DipoleSplittingGenerator.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the DipoleSplittingGenerator class.
//
#include <config.h>
#include "DipoleSplittingGenerator.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/DipoleShower/DipoleShowerHandler.h"
using namespace Herwig;
DipoleSplittingGenerator::DipoleSplittingGenerator()
: HandlerBase(),
theExponentialGenerator(0), prepared(false), presampling(false),
theDoCompensate(false) {
if ( ShowerHandler::currentHandler() )
setGenerator(ShowerHandler::currentHandler()->generator());
}
DipoleSplittingGenerator::~DipoleSplittingGenerator() {
if ( theExponentialGenerator ) {
delete theExponentialGenerator;
theExponentialGenerator = 0;
}
}
IBPtr DipoleSplittingGenerator::clone() const {
return new_ptr(*this);
}
IBPtr DipoleSplittingGenerator::fullclone() const {
return new_ptr(*this);
}
void DipoleSplittingGenerator::wrap(Ptr<DipoleSplittingGenerator>::ptr other) {
assert(!prepared);
theOtherGenerator = other;
}
void DipoleSplittingGenerator::resetVariations() {
for ( map<string,double>::iterator w = currentWeights.begin();
w != currentWeights.end(); ++w )
w->second = 1.;
}
void DipoleSplittingGenerator::prepare(const DipoleSplittingInfo& sp) {
generatedSplitting = sp;
generatedSplitting.splittingKinematics(splittingKernel()->splittingKinematics());
generatedSplitting.splittingParameters().resize(splittingKernel()->nDimAdditional());
if ( wrapping() ) {
generatedSplitting.emitterData(theSplittingKernel->emitter(generatedSplitting.index()));
generatedSplitting.spectatorData(theSplittingKernel->spectator(generatedSplitting.index()));
generatedSplitting.emissionData(theSplittingKernel->emission(generatedSplitting.index()));
parameters.resize(theOtherGenerator->nDim());
prepared = true;
return;
}
generatedSplitting.emitterData(splittingKernel()->emitter(generatedSplitting.index()));
generatedSplitting.spectatorData(splittingKernel()->spectator(generatedSplitting.index()));
generatedSplitting.emissionData(splittingKernel()->emission(generatedSplitting.index()));
presampledSplitting = generatedSplitting;
prepared = true;
parameters.resize(nDim());
theExponentialGenerator =
new exsample::exponential_generator<DipoleSplittingGenerator,UseRandom>();
theExponentialGenerator->sampling_parameters().maxtry = maxtry();
theExponentialGenerator->sampling_parameters().presampling_points = presamplingPoints();
theExponentialGenerator->sampling_parameters().freeze_grid = freezeGrid();
+ theExponentialGenerator->detuning(detuning());
theExponentialGenerator->docompensate(theDoCompensate);
theExponentialGenerator->function(this);
theExponentialGenerator->initialize();
}
void DipoleSplittingGenerator::fixParameters(const DipoleSplittingInfo& sp,
Energy optHardPt) {
assert(generator());
assert(!presampling);
assert(prepared);
assert(sp.index() == generatedSplitting.index());
generatedSplitting.scale(sp.scale());
parameters[3] = sp.scale()/generator()->maximumCMEnergy();
generatedSplitting.hardPt(sp.hardPt());
parameters[0] = splittingKinematics()->ptToRandom(optHardPt == ZERO ?
generatedSplitting.hardPt() :
min(generatedSplitting.hardPt(),optHardPt),
sp.scale(),
sp.emitterX(), sp.spectatorX(),
generatedSplitting.index(),
*splittingKernel());
size_t shift = 4;
if ( generatedSplitting.index().emitterPDF().pdf() &&
generatedSplitting.index().spectatorPDF().pdf() ) {
generatedSplitting.emitterX(sp.emitterX());
generatedSplitting.spectatorX(sp.spectatorX());
parameters[4] = sp.emitterX();
parameters[5] = sp.spectatorX();
shift += 2;
}
if ( generatedSplitting.index().emitterPDF().pdf() &&
!generatedSplitting.index().spectatorPDF().pdf() ) {
generatedSplitting.emitterX(sp.emitterX());
parameters[4] = sp.emitterX();
++shift;
}
if ( !generatedSplitting.index().emitterPDF().pdf() &&
generatedSplitting.index().spectatorPDF().pdf() ) {
generatedSplitting.spectatorX(sp.spectatorX());
parameters[4] = sp.spectatorX();
++shift;
}
if ( splittingReweight() ) {
parameters[shift] = splittingReweight()->evaluate(sp);
++shift;
}
if ( splittingKernel()->nDimAdditional() )
copy(sp.lastSplittingParameters().begin(),sp.lastSplittingParameters().end(),parameters.begin()+shift);
if ( sp.emitter() )
generatedSplitting.emitter(sp.emitter());
if ( sp.spectator() )
generatedSplitting.spectator(sp.spectator());
}
int DipoleSplittingGenerator::nDim() const {
assert(!wrapping());
assert(prepared);
int ret = 4; // 0 pt, 1 z, 2 phi, 3 scale, 4/5 xs + parameters
if ( generatedSplitting.index().emitterPDF().pdf() ) {
++ret;
}
if ( generatedSplitting.index().spectatorPDF().pdf() ) {
++ret;
}
if ( splittingReweight() )
++ret;
ret += splittingKernel()->nDimAdditional();
return ret;
}
const vector<bool>& DipoleSplittingGenerator::sampleFlags() {
assert(!wrapping());
if ( !theFlags.empty() )
return theFlags;
theFlags.resize(nDim(),false);
theFlags[0] = true; theFlags[1] = true; theFlags[2] = true; // 0 pt, 1 z, 2 phi
return theFlags;
}
const pair<vector<double>,vector<double> >& DipoleSplittingGenerator::support() {
assert(!wrapping());
if ( !theSupport.first.empty() )
return theSupport;
vector<double> lower(nDim(),0.);
vector<double> upper(nDim(),1.);
pair<double,double> kSupport =
generatedSplitting.splittingKinematics()->kappaSupport(generatedSplitting);
pair<double,double> xSupport =
generatedSplitting.splittingKinematics()->xiSupport(generatedSplitting);
lower[0] = kSupport.first;
lower[1] = xSupport.first;
upper[0] = kSupport.second;
upper[1] = xSupport.second;
if ( splittingReweight() ) {
pair<double,double> bounds = splittingReweight()->reweightBounds(generatedSplitting.index());
int pos = 4;
if ( generatedSplitting.index().emitterPDF().pdf() ) {
++pos;
}
if ( generatedSplitting.index().spectatorPDF().pdf() ) {
++pos;
}
lower[pos] = bounds.first;
upper[pos] = bounds.second;
}
theSupport.first = lower;
theSupport.second = upper;
return theSupport;
}
void DipoleSplittingGenerator::startPresampling() {
assert(!wrapping());
splittingKernel()->startPresampling(generatedSplitting.index());
presampling = true;
}
void DipoleSplittingGenerator::stopPresampling() {
assert(!wrapping());
splittingKernel()->stopPresampling(generatedSplitting.index());
presampling = false;
}
bool DipoleSplittingGenerator::haveOverestimate() const {
assert(!wrapping());
assert(prepared);
return
generatedSplitting.splittingKinematics()->haveOverestimate() &&
splittingKernel()->haveOverestimate(generatedSplitting);
}
bool DipoleSplittingGenerator::overestimate(const vector<double>& point) {
assert(!wrapping());
assert(prepared);
assert(!presampling);
assert(haveOverestimate());
if ( ! generatedSplitting.splittingKinematics()->generateSplitting(point[0],point[1],point[2],
generatedSplitting,
*splittingKernel()) )
return 0.;
generatedSplitting.splittingKinematics()->prepareSplitting(generatedSplitting);
return
( generatedSplitting.splittingKinematics()->jacobianOverestimate() *
splittingKernel()->overestimate(generatedSplitting) *
(splittingReweight() ? splittingReweight()->evaluate(generatedSplitting) : 1.) );
}
double DipoleSplittingGenerator::invertOverestimateIntegral(double value) const {
assert(!wrapping());
assert(prepared);
assert(!presampling);
assert(haveOverestimate());
return
splittingKernel()->invertOverestimateIntegral(generatedSplitting,value);
}
double DipoleSplittingGenerator::evaluate(const vector<double>& point) {
assert(!wrapping());
assert(prepared);
assert(generator());
DipoleSplittingInfo& split =
( !presampling ? generatedSplitting : presampledSplitting );
split.continuesEvolving();
size_t shift = 4;
if ( presampling ) {
split.scale(point[3] * generator()->maximumCMEnergy());
if ( split.index().emitterPDF().pdf() &&
split.index().spectatorPDF().pdf() ) {
split.emitterX(point[4]);
split.spectatorX(point[5]);
shift += 2;
}
if ( split.index().emitterPDF().pdf() &&
!split.index().spectatorPDF().pdf() ) {
split.emitterX(point[4]);
++shift;
}
if ( !split.index().emitterPDF().pdf() &&
split.index().spectatorPDF().pdf() ) {
split.spectatorX(point[4]);
++shift;
}
if ( splittingReweight() )
++shift;
if ( splittingKernel()->nDimAdditional() )
copy(point.begin()+shift,point.end(),split.splittingParameters().begin());
split.hardPt(split.splittingKinematics()->ptMax(split.scale(),
split.emitterX(),
split.spectatorX(),
split.index(),
*splittingKernel()));
}
if ( ! split.splittingKinematics()->generateSplitting(point[0],point[1],point[2],split,*splittingKernel()) ) {
split.lastValue(0.);
return 0.;
}
split.splittingKinematics()->prepareSplitting(split);
if ( split.stoppedEvolving() ) {
split.lastValue(0.);
return 0.;
}
if ( !presampling )
splittingKernel()->clearAlphaPDFCache();
double kernel = splittingKernel()->evaluate(split);
if ( splittingReweight() ) {
if ( !presampling )
kernel *= splittingReweight()->evaluate(split);
else
kernel *= point[shift-1];
}
double jac = split.splittingKinematics()->jacobian();
// multiply in the profile scales when relevant
assert(ShowerHandler::currentHandler());
if ( ShowerHandler::currentHandler()->firstInteraction() &&
ShowerHandler::currentHandler()->profileScales() &&
!presampling ) {
Energy hard = ShowerHandler::currentHandler()->hardScale();
if ( hard > ZERO )
kernel *= ShowerHandler::currentHandler()->profileScales()->hardScaleProfile(hard,split.lastPt());
}
split.lastValue( abs(jac) * kernel );
if ( isnan(split.lastValue()) || isinf(split.lastValue()) ) {
generator()->log() << "DipoleSplittingGenerator:evaluate(): problematic splitting kernel encountered for "
<< splittingKernel()->name() << "\n" << flush;
split.lastValue(0.0);
}
if ( kernel < 0. )
return 0.;
return split.lastValue();
}
void DipoleSplittingGenerator::doGenerate(map<string,double>& variations,
Energy optCutoff) {
assert(!wrapping());
double res = 0.;
Energy startPt = generatedSplitting.hardPt();
double optKappaCutoff = 0.0;
if ( optCutoff > splittingKinematics()->IRCutoff() ) {
optKappaCutoff = splittingKinematics()->ptToRandom(optCutoff,
generatedSplitting.scale(),
generatedSplitting.emitterX(),
generatedSplitting.spectatorX(),
generatedSplitting.index(),
*splittingKernel());
}
resetVariations();
while (true) {
try {
if ( optKappaCutoff == 0.0 ) {
res = theExponentialGenerator->generate();
} else {
res = theExponentialGenerator->generate(optKappaCutoff);
}
} catch (exsample::exponential_regenerate&) {
resetVariations();
generatedSplitting.hardPt(startPt);
continue;
} catch (exsample::hit_and_miss_maxtry&) {
throw DipoleShowerHandler::RedoShower();
} catch (exsample::selection_maxtry&) {
throw DipoleShowerHandler::RedoShower();
}
break;
}
for ( map<string,double>::const_iterator w = currentWeights.begin();
w != currentWeights.end(); ++w ) {
map<string,double>::iterator v = variations.find(w->first);
if ( v != variations.end() )
v->second *= w->second;
else
variations[w->first] = w->second;
}
if ( res == 0. ) {
generatedSplitting.lastPt(0.0*GeV);
generatedSplitting.didStopEvolving();
} else {
generatedSplitting.continuesEvolving();
if ( theMCCheck )
theMCCheck->book(generatedSplitting.emitterX(),
generatedSplitting.spectatorX(),
generatedSplitting.scale(),
startPt,
generatedSplitting.lastPt(),
generatedSplitting.lastZ(),
1.);
}
}
Energy DipoleSplittingGenerator::generate(const DipoleSplittingInfo& split,
map<string,double>& variations,
Energy optHardPt,
Energy optCutoff) {
fixParameters(split,optHardPt);
if ( wrapping() ) {
return theOtherGenerator->generateWrapped(generatedSplitting,variations,optHardPt,optCutoff);
}
doGenerate(variations,optCutoff);
return generatedSplitting.lastPt();
}
Energy DipoleSplittingGenerator::generateWrapped(DipoleSplittingInfo& split,
map<string,double>& variations,
Energy optHardPt,
Energy optCutoff) {
assert(!wrapping());
DipoleSplittingInfo backup = generatedSplitting;
generatedSplitting = split;
fixParameters(split,optHardPt);
try {
doGenerate(variations,optCutoff);
} catch (...) {
split = generatedSplitting;
generatedSplitting = backup;
throw;
}
Energy pt = generatedSplitting.lastPt();
split = generatedSplitting;
generatedSplitting = backup;
return pt;
}
void DipoleSplittingGenerator::completeSplitting(DipoleSplittingInfo& sp) const {
pair<bool,bool> conf = sp.configuration();
sp = generatedSplitting;
sp.configuration(conf);
}
Ptr<DipoleSplittingKernel>::tptr DipoleSplittingGenerator::splittingKernel() const {
if ( wrapping() )
return theOtherGenerator->splittingKernel();
return theSplittingKernel;
}
Ptr<DipoleSplittingReweight>::tptr DipoleSplittingGenerator::splittingReweight() const {
if ( wrapping() )
return theOtherGenerator->splittingReweight();
return theSplittingReweight;
}
Ptr<DipoleSplittingKinematics>::tptr DipoleSplittingGenerator::splittingKinematics() const {
if ( wrapping() )
return theOtherGenerator->splittingKinematics();
return theSplittingKernel->splittingKinematics();
}
void DipoleSplittingGenerator::splittingKernel(Ptr<DipoleSplittingKernel>::tptr sp) {
theSplittingKernel = sp;
if ( theSplittingKernel->mcCheck() )
theMCCheck = theSplittingKernel->mcCheck();
}
void DipoleSplittingGenerator::splittingReweight(Ptr<DipoleSplittingReweight>::tptr sp) {
theSplittingReweight = sp;
}
void DipoleSplittingGenerator::debugGenerator(ostream& os) const {
os << "--- DipoleSplittingGenerator ---------------------------------------------------\n";
os << " generating splittings using\n"
<< " splittingKernel = " << splittingKernel()->name()
<< " splittingKinematics = " << generatedSplitting.splittingKinematics()->name() << "\n"
<< " to sample splittings of type:\n";
os << generatedSplitting;
os << "--------------------------------------------------------------------------------\n";
}
void DipoleSplittingGenerator::debugLastEvent(ostream& os) const {
os << "--- DipoleSplittingGenerator ---------------------------------------------------\n";
os << " last generated event:\n";
os << generatedSplitting;
os << "--------------------------------------------------------------------------------\n";
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void DipoleSplittingGenerator::persistentOutput(PersistentOStream & os) const {
os << theOtherGenerator << theSplittingKernel << theSplittingReweight << theMCCheck << theDoCompensate;
}
void DipoleSplittingGenerator::persistentInput(PersistentIStream & is, int) {
is >> theOtherGenerator >> theSplittingKernel >> theSplittingReweight >> theMCCheck >> theDoCompensate;
}
ClassDescription<DipoleSplittingGenerator> DipoleSplittingGenerator::initDipoleSplittingGenerator;
// Definition of the static class description member.
void DipoleSplittingGenerator::Init() {
static ClassDocumentation<DipoleSplittingGenerator> documentation
("DipoleSplittingGenerator is used by the dipole shower "
"to sample splittings from a given dipole splitting kernel.");
static Reference<DipoleSplittingGenerator,DipoleSplittingKernel> interfaceSplittingKernel
("SplittingKernel",
"Set the splitting kernel to sample from.",
&DipoleSplittingGenerator::theSplittingKernel, false, false, true, false, false);
static Reference<DipoleSplittingGenerator,DipoleSplittingReweight> interfaceSplittingReweight
("SplittingReweight",
"Set the splitting reweight.",
&DipoleSplittingGenerator::theSplittingReweight, false, false, true, true, false);
static Reference<DipoleSplittingGenerator,DipoleMCCheck> interfaceMCCheck
("MCCheck",
"[debug option] MCCheck",
&DipoleSplittingGenerator::theMCCheck, false, false, true, true, false);
interfaceMCCheck.rank(-1);
}
diff --git a/DipoleShower/Base/DipoleSplittingGenerator.h b/DipoleShower/Base/DipoleSplittingGenerator.h
--- a/DipoleShower/Base/DipoleSplittingGenerator.h
+++ b/DipoleShower/Base/DipoleSplittingGenerator.h
@@ -1,476 +1,481 @@
// -*- C++ -*-
//
// DipoleSplittingGenerator.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_DipoleSplittingGenerator_H
#define HERWIG_DipoleSplittingGenerator_H
//
// This is the declaration of the DipoleSplittingGenerator class.
//
#include "ThePEG/Handlers/HandlerBase.h"
#include "Herwig/DipoleShower/Kernels/DipoleSplittingKernel.h"
#include "DipoleSplittingReweight.h"
#include "Herwig/DipoleShower/Utility/DipoleMCCheck.h"
#include "Herwig/Sampling/exsample/exponential_generator.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup DipoleShower
* \author Simon Platzer
*
* \brief DipoleSplittingGenerator is used by the dipole shower
* to sample splittings from a given dipole splitting kernel.
*
* @see \ref DipoleSplittingGeneratorInterfaces "The interfaces"
* defined for DipoleSplittingGenerator.
*/
class DipoleSplittingGenerator: public HandlerBase {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
DipoleSplittingGenerator();
/**
* The destructor.
*/
virtual ~DipoleSplittingGenerator();
//@}
public:
/**
* Return the dipole splitting kernel.
*/
Ptr<DipoleSplittingKernel>::tptr splittingKernel() const;
/**
* Return the dipole splitting reweight.
*/
Ptr<DipoleSplittingReweight>::tptr splittingReweight() const;
/**
* Return the dipole splitting kinematics.
*/
Ptr<DipoleSplittingKinematics>::tptr splittingKinematics() const;
/**
* Set the dipole splitting kernel.
*/
void splittingKernel(Ptr<DipoleSplittingKernel>::tptr sp);
/**
* Set the dipole splitting reweight.
*/
void splittingReweight(Ptr<DipoleSplittingReweight>::tptr sp);
/**
* Make a wrapper around another generator.
*/
void wrap(Ptr<DipoleSplittingGenerator>::ptr other);
/**
* Return true, if this is actually a wrapper around
* another splitting generator.
*/
bool wrapping() const { return theOtherGenerator; }
public:
/**
* Reset the current variations to one
*/
void resetVariations();
/**
* Prepare to fill the given splitting.
*/
void prepare(const DipoleSplittingInfo&);
/**
* Fix parameters from the fiven DipoleSplittingInfo
* and generate the next splitting. Return the
* pt selected for the next splitting.
*/
Energy generate(const DipoleSplittingInfo&,
map<string,double>& variations,
Energy optHardPt = ZERO,
Energy optCutoff = ZERO);
/**
* Fix parameters from the fiven DipoleSplittingInfo
* and generate the next splitting. Return the
* pt selected for the next splitting when called
* from a wrapping generator.
*/
Energy generateWrapped(DipoleSplittingInfo&,
map<string,double>& variations,
Energy optHardPt = ZERO,
Energy optCutoff = ZERO);
/**
* Complete the given splitting.
*/
void completeSplitting(DipoleSplittingInfo&) const;
/**
* Return the last generated splitting
*/
const DipoleSplittingInfo& lastSplitting() const { return generatedSplitting; }
public:
/**
* Print debug information on the splitting
* handled.
*/
void debugGenerator(ostream&) const;
/**
* Print debug information on the last
* generated event.
*/
void debugLastEvent(ostream&) const;
protected:
/**
* Update parameters given a splitting.
*/
void fixParameters(const DipoleSplittingInfo&,
Energy optHardPt = ZERO);
/**
* With the parameters previuosly supplied
* through fixParameters generate the next
* splitting.
*/
void doGenerate(map<string,double>& variations,
Energy optCutoff = ZERO);
public:
/**
* Return the number of random numbers
* needed to sample this kernel.
*/
int nDim() const;
/**
* Flag, which variables are free variables.
*/
const vector<bool>& sampleFlags();
/**
* Return the support of the splitting kernel.
* The lower bound on the first variable is
* assumed to correspond to the cutoff on the
* evolution variable.
*/
const pair<vector<double>,vector<double> >& support();
/**
* Return the parameter point associated to the splitting
* previously supplied through fixParameters.
*/
const vector<double>& parameterPoint() const { return parameters; }
/**
* Indicate that presampling of this kernel
* will be performed in the next calls to
* evaluate until stopPresampling() is called.
*/
void startPresampling();
/**
* Indicate that presampling of this kernel
* is done until startPresampling() is called.
*/
void stopPresampling();
/**
* Return the number of points to presample this
* splitting generator.
*/
unsigned long presamplingPoints() const { return splittingKernel()->presamplingPoints(); }
/**
* Return the maximum number of trials
* to generate a splitting.
*/
unsigned long maxtry() const { return splittingKernel()->maxtry(); }
/**
* Return the number of accepted points after which the grid should
* be frozen
*/
unsigned long freezeGrid() const { return splittingKernel()->freezeGrid(); }
/**
+ * Return the detuning factor applied to the sampling overestimate kernel
+ */
+ double detuning() const { return splittingKernel()->detuning(); }
+
+ /**
* Return true, if this splitting generator
* is able to deliver an overestimate to the sampled
* kernel.
*/
bool haveOverestimate() const;
/**
* Return an overestimate to the sampled kernel.
*/
bool overestimate(const vector<double>&);
/**
* Invert the integral over the overestimate to equal
* the given value.
*/
double invertOverestimateIntegral(double) const;
/**
* Evalute the splitting kernel.
*/
double evaluate(const vector<double>&);
/**
* Indicate that a veto with the given kernel value and overestimate has occured.
*/
void veto(const vector<double>&, double p, double r) {
splittingKernel()->veto(generatedSplitting, p, r, currentWeights);
}
/**
* Indicate that an accept with the given kernel value and overestimate has occured.
*/
void accept(const vector<double>&, double p, double r) {
splittingKernel()->accept(generatedSplitting, p, r, currentWeights);
}
/**
* True, if sampler should apply compensation
*/
void doCompensate(bool yes = true) { theDoCompensate = yes; }
public:
/**@name Wrap to the exsample2 interface until this is finally cleaned up. */
//@{
inline const vector<bool>& variable_flags () {
return sampleFlags();
}
inline size_t evolution_variable () const { return 0; }
inline double evolution_cutoff () { return support().first[0]; }
inline const vector<double>& parameter_point () const {
return parameterPoint();
}
inline void start_presampling () {
startPresampling();
}
inline void stop_presampling () {
stopPresampling();
}
inline size_t dimension () const {
return nDim();
}
inline unsigned long presampling_points () const {
return presamplingPoints();
}
//@}
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).
private:
/**
* Pointer to another generator to wrap around.
*/
Ptr<DipoleSplittingGenerator>::ptr theOtherGenerator;
/**
* The dipole splitting kernel to sample
* splitting from.
*/
Ptr<DipoleSplittingKernel>::ptr theSplittingKernel;
/**
* The dipole splitting reweight.
*/
Ptr<DipoleSplittingReweight>::ptr theSplittingReweight;
/**
* Pointer to the exponential generator
*/
exsample::exponential_generator<DipoleSplittingGenerator,UseRandom>*
theExponentialGenerator;
/**
* The dipole splitting to be completed.
*/
DipoleSplittingInfo generatedSplitting;
/**
* A backup of the dipole splitting to be
* completed, if this generator is presampled.
*/
DipoleSplittingInfo presampledSplitting;
/**
* True, if prepared to sample splittings
* of a given kind.
*/
bool prepared;
/**
* Wether or not the kernel is currently
* being presampled.
*/
bool presampling;
/**
* The parameter point.
*/
vector<double> parameters;
/**
* The sampling flags
*/
vector<bool> theFlags;
/**
* The support.
*/
pair<vector<double>,vector<double> > theSupport;
/**
* Pointer to a check histogram object
*/
Ptr<DipoleMCCheck>::ptr theMCCheck;
/**
* True, if sampler should apply compensation
*/
bool theDoCompensate;
/**
* The currently used weight map
*/
map<string,double> currentWeights;
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<DipoleSplittingGenerator> initDipoleSplittingGenerator;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
DipoleSplittingGenerator & operator=(const DipoleSplittingGenerator &);
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of DipoleSplittingGenerator. */
template <>
struct BaseClassTrait<Herwig::DipoleSplittingGenerator,1> {
/** Typedef of the first base class of DipoleSplittingGenerator. */
typedef HandlerBase NthBase;
};
/** This template specialization informs ThePEG about the name of
* the DipoleSplittingGenerator class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::DipoleSplittingGenerator>
: public ClassTraitsBase<Herwig::DipoleSplittingGenerator> {
/** Return a platform-independent class name */
static string className() { return "Herwig::DipoleSplittingGenerator"; }
/**
* The name of a file containing the dynamic library where the class
* DipoleSplittingGenerator is implemented. It may also include several, space-separated,
* libraries if the class DipoleSplittingGenerator depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwDipoleShower.so"; }
};
/** @endcond */
}
#endif /* HERWIG_DipoleSplittingGenerator_H */
diff --git a/DipoleShower/Kernels/DipoleSplittingKernel.cc b/DipoleShower/Kernels/DipoleSplittingKernel.cc
--- a/DipoleShower/Kernels/DipoleSplittingKernel.cc
+++ b/DipoleShower/Kernels/DipoleSplittingKernel.cc
@@ -1,362 +1,369 @@
// -*- C++ -*-
//
// DipoleSplittingKernel.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the DipoleSplittingKernel class.
//
#include "DipoleSplittingKernel.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Shower/ShowerHandler.h"
using namespace Herwig;
DipoleSplittingKernel::DipoleSplittingKernel()
: HandlerBase(), theScreeningScale(0.0*GeV),
thePresamplingPoints(2000), theMaxtry(100000),
theFreezeGrid(500000),
+ theDetuning(1.0),
theStrictLargeN(false),
theFactorizationScaleFactor(1.0),
theRenormalizationScaleFactor(1.0),
theRenormalizationScaleFreeze(1.*GeV),
theFactorizationScaleFreeze(1.*GeV),
theVirtualitySplittingScale(false),
presampling(false) {}
DipoleSplittingKernel::~DipoleSplittingKernel() {}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void DipoleSplittingKernel::persistentOutput(PersistentOStream & os) const {
os << theAlphaS << ounit(theScreeningScale,GeV) << theSplittingKinematics << thePDFRatio
- << thePresamplingPoints << theMaxtry << theFreezeGrid
+ << thePresamplingPoints << theMaxtry << theFreezeGrid << theDetuning
<< theFlavour << theMCCheck << theStrictLargeN
<< theFactorizationScaleFactor
<< theRenormalizationScaleFactor
<< ounit(theRenormalizationScaleFreeze,GeV)
<< ounit(theFactorizationScaleFreeze,GeV)
<< theVirtualitySplittingScale;
}
void DipoleSplittingKernel::persistentInput(PersistentIStream & is, int) {
is >> theAlphaS >> iunit(theScreeningScale,GeV) >> theSplittingKinematics >> thePDFRatio
- >> thePresamplingPoints >> theMaxtry >> theFreezeGrid
+ >> thePresamplingPoints >> theMaxtry >> theFreezeGrid >> theDetuning
>> theFlavour >> theMCCheck >> theStrictLargeN
>> theFactorizationScaleFactor
>> theRenormalizationScaleFactor
>> iunit(theRenormalizationScaleFreeze,GeV)
>> iunit(theFactorizationScaleFreeze,GeV)
>> theVirtualitySplittingScale;
}
double DipoleSplittingKernel::alphaPDF(const DipoleSplittingInfo& split,
Energy optScale,
double rScaleFactor,
double fScaleFactor) const {
Energy pt = optScale == ZERO ? split.lastPt() : optScale;
Energy2 scale = ZERO;
if ( !virtualitySplittingScale() ) {
scale = sqr(pt) + sqr(theScreeningScale);
} else {
scale = sqr(splittingKinematics()->QFromPt(pt,split)) + sqr(theScreeningScale);
}
Energy2 rScale = sqr(theRenormalizationScaleFactor*rScaleFactor)*scale;
rScale = rScale > sqr(renormalizationScaleFreeze()) ? rScale : sqr(renormalizationScaleFreeze());
Energy2 fScale = sqr(theFactorizationScaleFactor*fScaleFactor)*scale;
fScale = fScale > sqr(factorizationScaleFreeze()) ? fScale : sqr(factorizationScaleFreeze());
double alphas = 1.0;
double pdf = 1.0;
// check if we are potentially reweighting and cache evaluations
bool evaluatePDF = true;
bool evaluateAlphaS = true;
bool variations =
!ShowerHandler::currentHandler()->showerVariations().empty() &&
!presampling;
if ( variations ) {
/*
cerr << "--------------------------------------------------------------------------------\n"
<< flush;
cerr << "variations ... central scale/GeV = "
<< sqrt(scale/GeV2) << " r = "
<< rScaleFactor << " f = " << fScaleFactor << "\n"
<< flush;
*/
map<double,double>::const_iterator pit = thePDFCache.find(fScaleFactor);
evaluatePDF = (pit == thePDFCache.end());
if ( !evaluatePDF ) {
//cerr << "PDF is cached: ";
pdf = pit->second;
//cerr << pdf << "\n" << flush;
}
map<double,double>::const_iterator ait = theAlphaSCache.find(rScaleFactor);
evaluateAlphaS = (ait == theAlphaSCache.end());
if ( !evaluateAlphaS ) {
//cerr << "alphas is cached: ";
alphas = ait->second;
//cerr << alphas << "\n" << flush;
}
}
if ( evaluateAlphaS )
alphas = alphaS()->value(rScale);
if ( evaluatePDF ) {
if ( split.index().initialStateEmitter() ) {
assert(pdfRatio());
pdf *=
split.lastEmitterZ() *
(*pdfRatio())(split.index().emitterPDF(), fScale,
split.index().emitterData(),split.emitterData(),
split.emitterX(),split.lastEmitterZ());
}
if ( split.index().initialStateSpectator() ) {
assert(pdfRatio());
pdf *=
split.lastSpectatorZ() *
(*pdfRatio())(split.index().spectatorPDF(), fScale,
split.index().spectatorData(),split.spectatorData(),
split.spectatorX(),split.lastSpectatorZ());
}
}
if ( evaluatePDF && variations ) {
//cerr << "caching PDF = " << pdf << "\n" << flush;
thePDFCache[fScaleFactor] = pdf;
}
if ( evaluateAlphaS && variations ) {
//cerr << "caching alphas = " << alphas << "\n" << flush;
theAlphaSCache[rScaleFactor] = alphas;
}
double ret = alphas*pdf / (2.*Constants::pi);
if ( ret < 0. )
ret = 0.;
return ret;
}
void DipoleSplittingKernel::accept(const DipoleSplittingInfo& split,
double, double,
map<string,double>& weights) const {
if ( ShowerHandler::currentHandler()->showerVariations().empty() )
return;
double reference = alphaPDF(split);
assert(reference > 0.);
for ( map<string,ShowerHandler::ShowerVariation>::const_iterator var =
ShowerHandler::currentHandler()->showerVariations().begin();
var != ShowerHandler::currentHandler()->showerVariations().end(); ++var ) {
if ( ( ShowerHandler::currentHandler()->firstInteraction() && var->second.firstInteraction ) ||
( !ShowerHandler::currentHandler()->firstInteraction() && var->second.secondaryInteractions ) ) {
/*
cerr << "reweighting in accept for: "
<< var->first
<< " using "
<< var->second.renormalizationScaleFactor << " "
<< var->second.factorizationScaleFactor << " "
<< var->second.firstInteraction << " "
<< var->second.secondaryInteractions << "\n" << flush;
*/
double varied = alphaPDF(split,ZERO,
var->second.renormalizationScaleFactor,
var->second.factorizationScaleFactor);
if ( varied != reference ) {
map<string,double>::iterator wi = weights.find(var->first);
if ( wi != weights.end() )
wi->second *= varied/reference;
else
weights[var->first] = varied/reference;
}
}
}
}
void DipoleSplittingKernel::veto(const DipoleSplittingInfo& split,
double p, double r,
map<string,double>& weights) const {
if ( ShowerHandler::currentHandler()->showerVariations().empty() )
return;
double reference = alphaPDF(split);
// this is dangerous, but we have no other choice currently -- need to
// carefully check for the effects; the assumption is that if the central
// one ius zero, then so will be the variations.
if ( reference == 0.0 )
return;
for ( map<string,ShowerHandler::ShowerVariation>::const_iterator var =
ShowerHandler::currentHandler()->showerVariations().begin();
var != ShowerHandler::currentHandler()->showerVariations().end(); ++var ) {
if ( ( ShowerHandler::currentHandler()->firstInteraction() && var->second.firstInteraction ) ||
( !ShowerHandler::currentHandler()->firstInteraction() && var->second.secondaryInteractions ) ) {
/*
cerr << "reweighting in veto for: "
<< var->first
<< " using "
<< var->second.renormalizationScaleFactor << " "
<< var->second.factorizationScaleFactor << " "
<< var->second.firstInteraction << " "
<< var->second.secondaryInteractions << "\n" << flush;
*/
double varied = alphaPDF(split,ZERO,
var->second.renormalizationScaleFactor,
var->second.factorizationScaleFactor);
if ( varied != reference ) {
map<string,double>::iterator wi = weights.find(var->first);
if ( wi != weights.end() )
wi->second *= (r - varied*p/reference) / (r-p);
else
weights[var->first] = (r - varied*p/reference) / (r-p);
}
}
}
}
AbstractClassDescription<DipoleSplittingKernel> DipoleSplittingKernel::initDipoleSplittingKernel;
// Definition of the static class description member.
void DipoleSplittingKernel::Init() {
static ClassDocumentation<DipoleSplittingKernel> documentation
("DipoleSplittingKernel is the base class for all kernels "
"used within the dipole shower.");
static Reference<DipoleSplittingKernel,AlphaSBase> interfaceAlphaS
("AlphaS",
"The strong coupling to be used by this splitting kernel.",
&DipoleSplittingKernel::theAlphaS, false, false, true, true, false);
static Parameter<DipoleSplittingKernel,Energy> interfaceScreeningScale
("ScreeningScale",
"A colour screening scale",
&DipoleSplittingKernel::theScreeningScale, GeV, 0.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Reference<DipoleSplittingKernel,DipoleSplittingKinematics> interfaceSplittingKinematics
("SplittingKinematics",
"The splitting kinematics to be used by this splitting kernel.",
&DipoleSplittingKernel::theSplittingKinematics, false, false, true, false, false);
static Reference<DipoleSplittingKernel,PDFRatio> interfacePDFRatio
("PDFRatio",
"Set the optional PDF ratio object to evaluate this kernel",
&DipoleSplittingKernel::thePDFRatio, false, false, true, true, false);
static Parameter<DipoleSplittingKernel,unsigned long> interfacePresamplingPoints
("PresamplingPoints",
"The number of points used to presample this kernel.",
&DipoleSplittingKernel::thePresamplingPoints, 2000, 1, 0,
false, false, Interface::lowerlim);
static Parameter<DipoleSplittingKernel,unsigned long> interfaceMaxtry
("Maxtry",
"The maximum number of attempts to generate a splitting.",
&DipoleSplittingKernel::theMaxtry, 10000, 1, 0,
false, false, Interface::lowerlim);
static Parameter<DipoleSplittingKernel,unsigned long> interfaceFreezeGrid
("FreezeGrid",
"",
&DipoleSplittingKernel::theFreezeGrid, 500000, 1, 0,
false, false, Interface::lowerlim);
static Reference<DipoleSplittingKernel,ParticleData> interfaceFlavour
("Flavour",
"Set the flavour to be produced if ambiguous.",
&DipoleSplittingKernel::theFlavour, false, false, true, true, false);
static Reference<DipoleSplittingKernel,DipoleMCCheck> interfaceMCCheck
("MCCheck",
"[debug option] MCCheck",
&DipoleSplittingKernel::theMCCheck, false, false, true, true, false);
interfaceMCCheck.rank(-1);
static Switch<DipoleSplittingKernel,bool> interfaceStrictLargeN
("StrictLargeN",
"Work in a strict large-N limit.",
&DipoleSplittingKernel::theStrictLargeN, false, false, false);
static SwitchOption interfaceStrictLargeNOn
(interfaceStrictLargeN,
"On",
"Replace C_F -> C_A/2 where present",
true);
static SwitchOption interfaceStrictLargeNOff
(interfaceStrictLargeN,
"Off",
"Keep C_F=4/3",
false);
interfaceStrictLargeN.rank(-2);
static Parameter<DipoleSplittingKernel,double> interfaceFactorizationScaleFactor
("FactorizationScaleFactor",
"The factorization scale factor.",
&DipoleSplittingKernel::theFactorizationScaleFactor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
interfaceFactorizationScaleFactor.rank(-2);
static Parameter<DipoleSplittingKernel,double> interfaceRenormalizationScaleFactor
("RenormalizationScaleFactor",
"The renormalization scale factor.",
&DipoleSplittingKernel::theRenormalizationScaleFactor, 1.0, 0.0, 0,
false, false, Interface::lowerlim);
interfaceRenormalizationScaleFactor.rank(-2);
static Parameter<DipoleSplittingKernel,Energy> interfaceRenormalizationScaleFreeze
("RenormalizationScaleFreeze",
"The freezing scale for the renormalization scale.",
&DipoleSplittingKernel::theRenormalizationScaleFreeze, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Parameter<DipoleSplittingKernel,Energy> interfaceFactorizationScaleFreeze
("FactorizationScaleFreeze",
"The freezing scale for the factorization scale.",
&DipoleSplittingKernel::theFactorizationScaleFreeze, GeV, 1.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Switch<DipoleSplittingKernel,bool> interfaceVirtualitySplittingScale
("VirtualitySplittingScale",
"Use the virtuality as the splitting scale.",
&DipoleSplittingKernel::theVirtualitySplittingScale, false, false, false);
static SwitchOption interfaceVirtualitySplittingScaleYes
(interfaceVirtualitySplittingScale,
"Yes",
"Use vrituality.",
true);
static SwitchOption interfaceVirtualitySplittingScaleNo
(interfaceVirtualitySplittingScale,
"No",
"Use transverse momentum.",
false);
+ static Parameter<DipoleSplittingKernel,double> interfaceDetuning
+ ("Detuning",
+ "A value to detune the overestimate kernel.",
+ &DipoleSplittingKernel::theDetuning, 1.0, 1.0, 0,
+ false, false, Interface::lowerlim);
+
}
diff --git a/DipoleShower/Kernels/DipoleSplittingKernel.h b/DipoleShower/Kernels/DipoleSplittingKernel.h
--- a/DipoleShower/Kernels/DipoleSplittingKernel.h
+++ b/DipoleShower/Kernels/DipoleSplittingKernel.h
@@ -1,480 +1,495 @@
// -*- C++ -*-
//
// DipoleSplittingKernel.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_DipoleSplittingKernel_H
#define HERWIG_DipoleSplittingKernel_H
//
// This is the declaration of the DipoleSplittingKernel class.
//
#include "ThePEG/Handlers/HandlerBase.h"
#include "ThePEG/StandardModel/AlphaSBase.h"
#include "ThePEG/PDF/PDF.h"
#include "Herwig/DipoleShower/Utility/PDFRatio.h"
#include "Herwig/DipoleShower/Base/DipoleSplittingInfo.h"
#include "Herwig/DipoleShower/Kinematics/DipoleSplittingKinematics.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup DipoleShower
* \author Simon Platzer
*
* \brief DipoleSplittingKernel is the base class for all kernels
* used within the dipole shower.
*
* @see \ref DipoleSplittingKernelInterfaces "The interfaces"
* defined for DipoleSplittingKernel.
*/
class DipoleSplittingKernel: public HandlerBase {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
DipoleSplittingKernel();
/**
* The destructor.
*/
virtual ~DipoleSplittingKernel();
//@}
public:
/**
* Return the alpha_s to be used
*/
Ptr<AlphaSBase>::tptr alphaS() const { return theAlphaS; }
/**
* Set the alpha_s to be used
*/
void alphaS(Ptr<AlphaSBase>::tptr ap) { theAlphaS = ap; }
/**
* Return the splitting kinematics object
*/
Ptr<DipoleSplittingKinematics>::tptr splittingKinematics() const { return theSplittingKinematics; }
/**
* Return the mc check object
*/
Ptr<DipoleMCCheck>::ptr mcCheck() const { return theMCCheck; }
/**
* Set the splitting kinematics object
*/
void splittingKinematics(Ptr<DipoleSplittingKinematics>::tptr sp) { theSplittingKinematics = sp; }
/**
* Return the PDFRatio object
*/
Ptr<PDFRatio>::tptr pdfRatio() const { return thePDFRatio; }
/**
* Set the PDFRatio object
*/
void pdfRatio(Ptr<PDFRatio>::tptr sp) { thePDFRatio = sp; }
/**
* Return the number of additional parameter
* random variables needed to evaluate this kernel
* except the momentum fractions of incoming partons.
* These will be accessible through the
* lastSplittingParameters() container of the splitting
* info object.
*/
virtual int nDimAdditional() const { return 0; }
/**
* Set the freezing value for the renormalization scale
*/
void renormalizationScaleFreeze(Energy s) { theRenormalizationScaleFreeze = s; }
/**
* Set the freezing value for the factorization scale
*/
void factorizationScaleFreeze(Energy s) { theFactorizationScaleFreeze = s; }
/**
* Get the freezing value for the renormalization scale
*/
Energy renormalizationScaleFreeze() const { return theRenormalizationScaleFreeze; }
/**
* Get the freezing value for the factorization scale
*/
Energy factorizationScaleFreeze() const { return theFactorizationScaleFreeze; }
public:
/**
* Return true, if this splitting kernel
* applies to the given dipole index.
*/
virtual bool canHandle(const DipoleIndex&) const = 0;
/**
* Return true, if this splitting kernel is
* the same for the given index a, as the given
* splitting kernel for index b.
*/
virtual bool canHandleEquivalent(const DipoleIndex& a,
const DipoleSplittingKernel& sk,
const DipoleIndex& b) const = 0;
/**
* Return the emitter data after splitting, given
* a dipole index.
*/
virtual tcPDPtr emitter(const DipoleIndex&) const = 0;
/**
* Return the emission data after splitting, given
* a dipole index.
*/
virtual tcPDPtr emission(const DipoleIndex&) const = 0;
/**
* Return the spectator data after splitting, given
* a dipole index.
*/
virtual tcPDPtr spectator(const DipoleIndex&) const = 0;
/**
* Return the flavour produced, if this cannot
* be determined from the dipole.
*/
PDPtr flavour() const { return theFlavour; }
/**
* Return true, if this splitting kernel is supposed to work in a
* strict large-N limit, i.e. replacing C_F by C_A/2
*/
bool strictLargeN() const { return theStrictLargeN; }
public:
/**
* Inform this splitting kernel, that it is being
* presampled until a call to stopPresampling
*/
virtual void startPresampling(const DipoleIndex&) {
presampling = true;
}
/**
* Inform this splitting kernel, that it is not being
* presampled until a call to startPresampling
*/
virtual void stopPresampling(const DipoleIndex&) {
presampling = false;
}
/**
* Return the number of points to presample this
* splitting generator.
*/
unsigned long presamplingPoints() const { return thePresamplingPoints; }
/**
* Return the maximum number of trials
* to generate a splitting.
*/
unsigned long maxtry() const { return theMaxtry; }
/**
* Return the number of accepted points after which the grid should
* be frozen
*/
unsigned long freezeGrid() const { return theFreezeGrid; }
/**
* Set the number of accepted points after which the grid should
* be frozen
*/
void freezeGrid(unsigned long n) { theFreezeGrid = n; }
/**
+ * Set a detuning factor to be applied to the sampling overestimate kernel
+ */
+ void detuning(double d) { theDetuning = d; }
+
+ /**
+ * Return the detuning factor applied to the sampling overestimate kernel
+ */
+ double detuning() const { return theDetuning; }
+
+ /**
* Evaluate this splitting kernel for the given
* dipole splitting.
*/
virtual double evaluate(const DipoleSplittingInfo&) const = 0;
/**
* Clear the alphaPDF cache
*/
void clearAlphaPDFCache() const {
theAlphaSCache.clear();
thePDFCache.clear();
}
/**
* Update the variations vector at the given splitting using the indicated
* kernel and overestimate values.
*/
virtual void accept(const DipoleSplittingInfo&, double, double, map<string,double>&) const;
/**
* Update the variations vector at the given splitting using the indicated
* kernel and overestimate values.
*/
virtual void veto(const DipoleSplittingInfo&, double, double, map<string,double>&) const;
/**
* Return true, if this kernel is capable of
* delivering an overestimate to the kernel, and
* of inverting the integral over the overestimate
* w.r.t. the phasepsace provided by the given
* DipoleSplittingInfo object.
*/
virtual bool haveOverestimate(const DipoleSplittingInfo&) const { return false; }
/**
* Return the overestimate to this splitting kernel
* for the given dipole splitting.
*/
virtual double overestimate(const DipoleSplittingInfo&) const { return -1.; }
/**
* Invert the integral over the overestimate
* w.r.t. the phasepsace provided by the given
* DipoleSplittingInfo object to equal
* the given value.
*/
virtual double invertOverestimateIntegral(const DipoleSplittingInfo&, double) const { return -1.; }
public:
/**
* Get the factorization scale factor
*/
double factorizationScaleFactor() const { return theFactorizationScaleFactor; }
/**
* Set the factorization scale factor
*/
void factorizationScaleFactor(double f) { theFactorizationScaleFactor = f; }
/**
* Get the renormalization scale factor
*/
double renormalizationScaleFactor() const { return theRenormalizationScaleFactor; }
/**
* Set the renormalization scale factor
*/
void renormalizationScaleFactor(double f) { theRenormalizationScaleFactor = f; }
protected:
/**
* Return the common factor of (alphas/2pi)*(pdf ratio)
*/
double alphaPDF(const DipoleSplittingInfo&,
Energy optScale = ZERO,
double rScaleFactor = 1.0,
double fScaleFactor = 1.0) const;
/**
* Return true, if the virtuality of the splitting should be used as the
* argument of alphas rather than the pt
*/
bool virtualitySplittingScale() const { return theVirtualitySplittingScale; }
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();
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The alpha_s to be used.
*/
Ptr<AlphaSBase>::ptr theAlphaS;
/**
* An optional 'colour screening' scale
* for alternative intrinsic pt generation.
*/
Energy theScreeningScale;
/**
* The splitting kinematics to be used.
*/
Ptr<DipoleSplittingKinematics>::ptr theSplittingKinematics;
/**
* An optional PDF ratio object to be used
* when evaluating this kernel.
*/
Ptr<PDFRatio>::ptr thePDFRatio;
/**
* The number of points to presample this
* splitting generator.
*/
unsigned long thePresamplingPoints;
/**
* The maximum number of trials
* to generate a splitting.
*/
unsigned long theMaxtry;
/**
* Return the number of accepted points after which the grid should
* be frozen
*/
unsigned long theFreezeGrid;
/**
+ * The detuning factor applied to the sampling overestimate kernel
+ */
+ double theDetuning;
+
+ /**
* The flavour produced, if this cannot
* be determined from the dipole.
*/
PDPtr theFlavour;
/**
* Pointer to a check histogram object
*/
Ptr<DipoleMCCheck>::ptr theMCCheck;
/**
* True, if this splitting kernel is supposed to work in a
* strict large-N limit, i.e. replacing C_F by C_A/2
*/
bool theStrictLargeN;
/**
* The factorization scale factor.
*/
double theFactorizationScaleFactor;
/**
* The renormalization scale factor.
*/
double theRenormalizationScaleFactor;
/**
* A freezing value for the renormalization scale
*/
Energy theRenormalizationScaleFreeze;
/**
* A freezing value for the factorization scale
*/
Energy theFactorizationScaleFreeze;
/**
* True, if the virtuality of the splitting should be used as the
* argument of alphas rather than the pt
*/
bool theVirtualitySplittingScale;
/**
* Cache for alphas evaluations
*/
mutable map<double,double> theAlphaSCache;
/**
* Cache for PDF evaluations
*/
mutable map<double,double> thePDFCache;
/**
* True, if we are presampling
*/
bool presampling;
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is an abstract class with persistent data.
*/
static AbstractClassDescription<DipoleSplittingKernel> initDipoleSplittingKernel;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
DipoleSplittingKernel & operator=(const DipoleSplittingKernel &);
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of DipoleSplittingKernel. */
template <>
struct BaseClassTrait<Herwig::DipoleSplittingKernel,1> {
/** Typedef of the first base class of DipoleSplittingKernel. */
typedef HandlerBase NthBase;
};
/** This template specialization informs ThePEG about the name of
* the DipoleSplittingKernel class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::DipoleSplittingKernel>
: public ClassTraitsBase<Herwig::DipoleSplittingKernel> {
/** Return a platform-independent class name */
static string className() { return "Herwig::DipoleSplittingKernel"; }
/**
* The name of a file containing the dynamic library where the class
* DipoleSplittingKernel is implemented. It may also include several, space-separated,
* libraries if the class DipoleSplittingKernel depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwDipoleShower.so"; }
};
/** @endcond */
}
#endif /* HERWIG_DipoleSplittingKernel_H */
diff --git a/Sampling/exsample/cell.h b/Sampling/exsample/cell.h
--- a/Sampling/exsample/cell.h
+++ b/Sampling/exsample/cell.h
@@ -1,304 +1,306 @@
// -*- C++ -*-
//
// cell.h is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
//
// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
//
// ExSample is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
#ifndef EXSAMPLE_cell_h_included
#define EXSAMPLE_cell_h_included
#include "utility.h"
#include "adaption_info.h"
#include "statistics.h"
namespace exsample {
/// \brief Information contained in a leaf cell
class cell_info {
public:
/// the default constructor
cell_info();
/// construct from boundaries and adaption info
cell_info(const std::vector<double>& ll,
const std::vector<double>& ur,
const adaption_info& ainfo);
/// construct from boundaries, flags for variables to be sampled,
/// and adaption info
cell_info(const std::vector<double>& ll,
const std::vector<double>& ur,
const std::vector<bool>& sampled_variables,
const adaption_info& ainfo);
public:
/// generate a flat trial point in this cell
template<class Random>
void select(Random&,
std::vector<double>&);
/// generate a flat trial point in this cell
/// only for the variables falgged as true
template<class Random>
void select(Random&,
std::vector<double>&,
const std::vector<bool>&);
/// indicate a function value for the given point
void selected(const std::vector<double>&,
double,
const adaption_info&);
/// indicate that a point has been
/// accepted in this cell
void accept() { ++accepted_; }
/// reject a previously accepted event
void reject() { --accepted_; }
public:
/// return true, if below efficiency threshold
bool bad(const adaption_info& ainfo) const {
return ((static_cast<double>(accepted_)/static_cast<double>(attempted_)) <
ainfo.efficiency_threshold);
}
/// suggest a split and indicate wether it is worth
/// to be performed
std::pair<std::size_t,double> get_split(const adaption_info&,
bool&) const;
/// explore this cell performing a flat sampling,
/// updating the given statistics object and pre-filling
/// the efficiency histogram by a trial unweighting
template<class Random, class Function, class SlaveStatistics>
void explore(Random&, const adaption_info&, Function*, statistics*,
- SlaveStatistics& opt);
+ SlaveStatistics& opt, double detuning = 1.0);
/// explore this cell in a more refined way, which
/// is however not suited for already calculating integrals
/// and stuff
template<class Random, class Function>
- void explore(Random&, const adaption_info&, Function*);
+ void explore(Random&, const adaption_info&, Function*,
+ double detuning = 1.0);
public:
/// get the current overestimate
double overestimate() const { return overestimate_; }
/// return the position of the last maximum
const std::vector<double>& last_max_position() const { return last_max_position_; }
/// set the current overestimate and maximum position
void overestimate(double v, const std::vector<double>& pos) {
overestimate_ = v;
last_max_position_ = pos;
}
/// get the volume
double volume() const { return volume_; }
/// get the lower left corner
const std::vector<double>& lower_left() const { return lower_left_; }
/// get the upper right corner
const std::vector<double>& upper_right() const { return upper_right_; }
/// get the number of attempted events
unsigned long attempted() const { return attempted_; }
/// get the number of accepted events
unsigned long accepted() const { return accepted_; }
public:
/// return the number of missing events
/// for the given parameter bin id
int parametric_missing(const bit_container<parameter_hash_bits>& id) const;
/// set the number of missing events
/// for the given parameter bin id
void parametric_missing(const bit_container<parameter_hash_bits>& id, int n);
/// increase to the number of missing events
/// for the given parameter bin id
void increase_parametric_missing(const bit_container<parameter_hash_bits>& id);
/// decrease to the number of missing events
/// for the given parameter bin id
void decrease_parametric_missing(const bit_container<parameter_hash_bits>& id);
/// return true, if the cell is compensating in
/// at least one parameter bin
bool parametric_compensating() const {
return !parametric_missing_map_.empty();
}
/// return true, if the cell contains the
/// indicated parameter point
bool contains_parameter(const std::vector<double>& point,
const std::vector<bool>& sampled) const;
public:
/// put to ostream
template<class OStream>
void put(OStream& os) const;
/// get from istream
template<class IStream>
void get(IStream& is);
private:
/// the value of the overestimate in this cell
double overestimate_;
/// the volume of this cell
double volume_;
/// the lower left corner of this cell
std::vector<double> lower_left_;
/// the upper right corner of this cell
std::vector<double> upper_right_;
/// midpoint of this cell
std::vector<double> mid_point_;
/// the position of the last encountered
/// maximum in this cell
std::vector<double> last_max_position_;
/// left-right statistics of average weight
std::vector<std::pair<double,double> > avg_weight_;
/// the number of attempts in this cell
unsigned long attempted_;
/// the number of accepted events in this cell
unsigned long accepted_;
/// an optional map of parameter bin ids
/// to the number of missing events
std::map<bit_container<parameter_hash_bits>,int> parametric_missing_map_;
};
/// \brief the general cell class
class cell {
public:
/// default constructor
cell();
/// construct from boundaries and adaption info
cell(const std::vector<double>& ll,
const std::vector<double>& ur,
const adaption_info& ainfo);
/// construct from boundaries, flags for variables to be sampled,
/// and adaption info
cell(const std::vector<double>& ll,
const std::vector<double>& ur,
const std::vector<bool>& sampled_variables,
const adaption_info& ainfo);
/// copy constructor
cell(const cell& x);
/// assignment
cell& operator=(const cell& x);
public:
/// split this cell, exploring the
/// child not containing the current overestimate
template<class Random, class Function>
std::pair<cell,cell> split(std::pair<std::size_t,double> split_d,
Random& rnd_gen,
Function* f,
const adaption_info& ainfo,
const std::vector<bool>& sampled =
- std::vector<bool>());
+ std::vector<bool>(),
+ double detuning = 1.0);
public:
/// return the split dimension
std::size_t split_dimension() const { return split_dimension_; }
/// return the split value
double split_point() const { return split_point_; }
/// return the integral
double integral() const { return integral_; }
/// access the integral
double& integral() { return integral_; }
/// set the integral
void integral(double v) { integral_ = v; }
/// access the number of missing events
int& missing_events() { return missing_events_; }
/// return the number of missing events
int missing_events() const { return missing_events_; }
/// set the number of missing events
void missing_events(int n) { missing_events_ = n; }
/// access the cell_info object
cell_info& info() { assert(cell_info_); return *cell_info_; }
/// return the cell_info object
const cell_info& info() const { assert(cell_info_); return *cell_info_; }
public:
/// put to ostream
template<class OStream>
void put(OStream& os) const;
/// get from istream
template<class IStream>
void get(IStream& is);
private:
/// the dimension along this cell
/// was split
std::size_t split_dimension_;
/// the value, where this cell was split
double split_point_;
/// the integral of the absolute value
/// of the overestimate over all the
/// children cells
double integral_;
/// the number of missing events in this cell
int missing_events_;
/// a pointer to the cell info object,
/// if this is a leaf cell
boost::scoped_ptr<cell_info> cell_info_;
};
}
#include "cell.icc"
#endif // EXSAMPLE_cell_h_included
diff --git a/Sampling/exsample/cell.icc b/Sampling/exsample/cell.icc
--- a/Sampling/exsample/cell.icc
+++ b/Sampling/exsample/cell.icc
@@ -1,487 +1,490 @@
// -*- C++ -*-
//
// cell.icc is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
//
// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
//
// ExSample is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
namespace exsample {
template<class Random>
void cell_info::select (Random& rnd_gen,
std::vector<double>& p) {
std::transform(lower_left_.begin(),lower_left_.end(),
upper_right_.begin(),p.begin(),
rnd_gen);
++attempted_;
}
template<class Random>
void cell_info::select (Random& rnd_gen,
std::vector<double>& p,
const std::vector<bool>& sample) {
conditional_transform(lower_left_.begin(),lower_left_.end(),
upper_right_.begin(),sample.begin(),
p.begin(),rnd_gen);
++attempted_;
}
template<class Random, class Function, class SlaveStatistics>
void cell_info::explore(Random& rnd_gen,
const adaption_info& ainfo,
Function* function, statistics* stats,
- SlaveStatistics& opt) {
+ SlaveStatistics& opt,
+ double detuning) {
function->start_presampling();
unsigned long n_sampled = 0;
std::vector<double> ll = lower_left_;
std::vector<double> ur = upper_right_;
double val = 0.;
std::vector<double> pos (ll.size());
std::vector< std::pair<double,std::vector<double> > > vals;
unsigned long ivalnonzero = 0;
while (n_sampled < ainfo.presampling_points) {
std::transform(ll.begin(),ll.end(),
ur.begin(),pos.begin(),
rnd_gen);
- val = function->evaluate(pos);
+ val = function->evaluate(pos) * detuning;
vals.push_back( std::pair<double,std::vector<double> > (val,pos) );
if ( val != 0 ) ivalnonzero++;
++n_sampled;
}
while ( ivalnonzero > 0 ) {
double avg = 0;
double err = 0;
double maxval = 0;
std::vector<double> maxpos (ll.size());
unsigned long imax;
for ( unsigned long ival=0; ival < vals.size(); ival++ ) {
val = std::abs(vals[ival].first);
if ( val == 0 ) continue;
avg += val;
err += sqr(val);
if ( val > maxval ) {
maxval = val;
maxpos = vals[ival].second;
imax = ival;
}
}
avg /= ivalnonzero;
err /= ivalnonzero;
err = sqrt(err-sqr(avg));
if ( maxval <= avg+sqrt(ivalnonzero/2.)*err ) {
overestimate_ = maxval;
last_max_position_ = maxpos;
break;
}
vals.erase(vals.begin()+imax);
ivalnonzero--;
}
for ( unsigned long ival=0; ival < vals.size(); ival++ ) {
val = vals[ival].first;
stats->presampled(val);
opt.select(val);
selected(pos,std::abs(val),ainfo);
}
function->stop_presampling();
}
template<class Random, class Function>
void cell_info::explore (Random& rnd_gen,
- const adaption_info& ainfo, Function* function) {
+ const adaption_info& ainfo, Function* function,
+ double detuning) {
function->start_presampling();
unsigned long n_sampled = 0;
std::vector<double> ll = lower_left_;
std::vector<double> ur = upper_right_;
double val = 0.;
std::vector<double> pos (ll.size());
std::vector< std::pair<double,std::vector<double> > > vals;
while (n_sampled < ainfo.presampling_points) {
std::transform(ll.begin(),ll.end(),
ur.begin(),pos.begin(),
rnd_gen);
- val = function->evaluate(pos);
+ val = function->evaluate(pos) * detuning;
if ( std::abs(val) > 0 )
vals.push_back( std::pair<double,std::vector<double> > (std::abs(val),pos) );
++n_sampled;
}
while ( vals.size() > 0 ) {
double avg = 0;
double err = 0;
double maxval = 0;
std::vector<double> maxpos (ll.size());
unsigned long imax;
for ( unsigned long ival=0; ival < vals.size(); ival++ ) {
double thisval = vals[ival].first;
avg += thisval;
err += sqr(thisval);
if ( thisval > maxval ) {
maxval = thisval;
maxpos = vals[ival].second;
imax = ival;
}
}
avg /= vals.size();
err /= vals.size();
err = sqrt(err-sqr(avg));
if ( maxval <= avg+sqrt(vals.size()/2.)*err ) {
overestimate_ = maxval;
last_max_position_ = maxpos;
break;
}
vals.erase(vals.begin()+imax);
}
function->stop_presampling();
}
template<class OStream>
void cell_info::put (OStream& os) const {
os << overestimate_;
ostream_traits<OStream>::separator(os);
os << volume_;
ostream_traits<OStream>::separator(os);
os << lower_left_.size();
ostream_traits<OStream>::separator(os);
for (std::size_t k = 0; k < lower_left_.size(); ++k) {
os << lower_left_[k];
ostream_traits<OStream>::separator(os);
}
for (std::size_t k = 0; k < upper_right_.size(); ++k) {
os << upper_right_[k];
ostream_traits<OStream>::separator(os);
}
for (std::size_t k = 0; k < mid_point_.size(); ++k) {
os << mid_point_[k];
ostream_traits<OStream>::separator(os);
}
for (std::size_t k = 0; k < last_max_position_.size(); ++k) {
os << last_max_position_[k];
ostream_traits<OStream>::separator(os);
}
for (std::size_t k = 0; k < avg_weight_.size(); ++k) {
os << avg_weight_[k].first;
ostream_traits<OStream>::separator(os);
os << avg_weight_[k].second;
ostream_traits<OStream>::separator(os);
}
os << attempted_;
ostream_traits<OStream>::separator(os);
os << accepted_;
ostream_traits<OStream>::separator(os);
os << parametric_missing_map_.size();
ostream_traits<OStream>::separator(os);
for ( std::map<bit_container<parameter_hash_bits>,int>::const_iterator p =
parametric_missing_map_.begin(); p != parametric_missing_map_.end(); ++p ) {
p->first.put(os);
os << p->second;
ostream_traits<OStream>::separator(os);
}
}
template<class IStream>
void cell_info::get (IStream& is) {
std::size_t dim;
is >> overestimate_ >> volume_ >> dim;
lower_left_.resize(dim);
for (std::size_t k = 0; k < lower_left_.size(); ++k) {
is >> lower_left_[k];
}
upper_right_.resize(dim);
for (std::size_t k = 0; k < upper_right_.size(); ++k) {
is >> upper_right_[k];
}
mid_point_.resize(dim);
for (std::size_t k = 0; k < mid_point_.size(); ++k) {
is >> mid_point_[k];
}
last_max_position_.resize(dim);
for (std::size_t k = 0; k < last_max_position_.size(); ++k) {
is >> last_max_position_[k];
}
avg_weight_.resize(dim);
for (std::size_t k = 0; k < avg_weight_.size(); ++k) {
is >> avg_weight_[k].first >> avg_weight_[k].second;
}
is >> attempted_ >> accepted_ >> dim;
for ( size_t k = 0; k < dim; ++k ) {
bit_container<parameter_hash_bits> in;
in.get(is);
is >> parametric_missing_map_[in];
}
}
template<class Random, class Function>
std::pair<cell,cell >
cell::split (std::pair<std::size_t,double> split_d,
Random& rnd_gen,
Function* function,
const adaption_info& ainfo,
- const std::vector<bool>& sampled) {
+ const std::vector<bool>& sampled,
+ double detuning) {
assert(!missing_events() && !info().parametric_compensating());
split_dimension_ = split_d.first;
split_point_ = split_d.second;
std::vector<double> lower_left1 = info().lower_left();
std::vector<double> upper_right1 = info().upper_right();
std::vector<double> lower_left2 = info().lower_left();
std::vector<double> upper_right2 = info().upper_right();
upper_right1[split_dimension_] = split_point_;
lower_left2[split_dimension_] = split_point_;
std::pair<cell,cell> children;
if (sampled.empty())
children = std::pair<cell,cell>(cell(lower_left1,upper_right1,ainfo),
cell(lower_left2,upper_right2,ainfo));
else
children = std::pair<cell,cell> (cell(lower_left1,upper_right1,sampled,ainfo),
cell(lower_left2,upper_right2,sampled,ainfo));
if (info().last_max_position()[split_dimension_] <= split_point_) {
children.first.info().overestimate(info().overestimate(),info().last_max_position());
- children.second.info().explore(rnd_gen,ainfo,function);
+ children.second.info().explore(rnd_gen,ainfo,function,detuning);
} else {
children.second.info().overestimate(info().overestimate(),info().last_max_position());
- children.first.info().explore(rnd_gen,ainfo,function);
+ children.first.info().explore(rnd_gen,ainfo,function,detuning);
}
cell_info_.reset(0);
children.first.integral(children.first.info().volume() * children.first.info().overestimate());
children.second.integral(children.second.info().volume() * children.second.info().overestimate());
return children;
}
template<class OStream>
void cell::put (OStream& os) const {
os << split_dimension_;
ostream_traits<OStream>::separator(os);
os << split_point_;
ostream_traits<OStream>::separator(os);
os << integral_;
ostream_traits<OStream>::separator(os);
os << missing_events_;
ostream_traits<OStream>::separator(os);
if (cell_info_) {
os << "has_cell_info";
ostream_traits<OStream>::separator(os);
cell_info_->put(os);
} else {
os << "has_no_cell_info";
ostream_traits<OStream>::separator(os);
}
}
template<class IStream>
void cell::get (IStream& is) {
std::string info_tag;
is >> split_dimension_ >> split_point_
>> integral_ >> missing_events_
>> info_tag;
if (info_tag == "has_cell_info") {
cell_info_.reset(new cell_info());
cell_info_->get(is);
}
}
inline cell_info::cell_info()
: overestimate_(0.), volume_(0.),
lower_left_(), upper_right_(), mid_point_(),
last_max_position_(), avg_weight_(),
attempted_(0), accepted_(0) {}
inline cell_info::cell_info(const std::vector<double>& ll,
const std::vector<double>& ur,
const adaption_info& ainfo)
: overestimate_(0.), volume_(),
lower_left_(ll), upper_right_(ur), mid_point_(),
last_max_position_(),
avg_weight_(std::vector<std::pair<double,double> >
(ainfo.dimension,std::make_pair(0.,0.))),
attempted_(0), accepted_(0) {
std::vector<double> delta;
std::transform(ur.begin(),ur.end(),
ll.begin(),std::back_inserter(delta),
std::minus<double>());
volume_ =
std::accumulate(delta.begin(),delta.end(),1.,std::multiplies<double>());
std::transform(ur.begin(),ur.end(),
ll.begin(),std::back_inserter(mid_point_),
std::plus<double>());
for (std::size_t k = 0; k < ainfo.dimension; ++k)
mid_point_[k] /= 2.;
}
inline cell_info::cell_info(const std::vector<double>& ll,
const std::vector<double>& ur,
const std::vector<bool>& sampled_variables,
const adaption_info& ainfo)
: overestimate_(0.), volume_(),
lower_left_(ll), upper_right_(ur), mid_point_(),
last_max_position_(),
avg_weight_(std::vector<std::pair<double,double> >
(ainfo.dimension,std::make_pair(0.,0.))),
attempted_(0), accepted_(0) {
std::vector<double> delta;
conditional_transform(ur.begin(),ur.end(),
ll.begin(),sampled_variables.begin(),
std::back_inserter(delta),
std::minus<double>());
volume_ =
std::accumulate(delta.begin(),delta.end(),1.,std::multiplies<double>());
std::transform(ur.begin(),ur.end(),
ll.begin(),std::back_inserter(mid_point_),
std::plus<double>());
for (std::size_t k = 0; k < ainfo.dimension; ++k)
mid_point_[k] /= 2.;
}
inline int cell_info::parametric_missing(const bit_container<parameter_hash_bits>& id) const {
std::map<bit_container<parameter_hash_bits>,int>::const_iterator mit
= parametric_missing_map_.find(id);
if (mit == parametric_missing_map_.end())
return 0;
return mit->second;
}
inline void cell_info::parametric_missing(const bit_container<parameter_hash_bits>& id, int n) {
if (n == 0) {
std::map<bit_container<parameter_hash_bits>,int>::iterator mit
= parametric_missing_map_.find(id);
if (mit != parametric_missing_map_.end())
parametric_missing_map_.erase(mit);
return;
}
parametric_missing_map_[id] = n;
}
inline void cell_info::increase_parametric_missing(const bit_container<parameter_hash_bits>& id) {
std::map<bit_container<parameter_hash_bits>,int>::iterator mit
= parametric_missing_map_.find(id);
if (mit != parametric_missing_map_.end()) {
mit->second += 1;
if (mit->second == 0) parametric_missing_map_.erase(mit);
} else parametric_missing_map_[id] = 1;
}
inline void cell_info::decrease_parametric_missing(const bit_container<parameter_hash_bits>& id) {
std::map<bit_container<parameter_hash_bits>,int>::iterator mit
= parametric_missing_map_.find(id);
if (mit != parametric_missing_map_.end()) {
mit->second -= 1;
if (mit->second == 0) parametric_missing_map_.erase(mit);
} else assert(false);
}
inline void cell_info::selected(const std::vector<double>& p,
double weight,
const adaption_info& ainfo) {
for (std::size_t k = 0; k < p.size(); ++k) {
if (ainfo.adapt[k]) {
if (p[k] < mid_point_[k])
avg_weight_[k].first += weight;
else
avg_weight_[k].second += weight;
}
}
}
inline std::pair<std::size_t,double> cell_info::get_split (const adaption_info& ainfo,
bool& worth) const {
std::size_t split_d = 0;
double gain = 0.;
for (std::size_t k = 0; k < ainfo.dimension; ++k) {
double xgain = 0.;
double left = avg_weight_[k].first;
double right = avg_weight_[k].second;
if (left+right > 0.) {
xgain = std::abs(left-right)/(left+right);
}
if (xgain > gain) {
gain = xgain;
split_d = k;
}
}
worth = (gain >= ainfo.gain_threshold);
return std::make_pair(split_d,mid_point_[split_d]);
}
inline bool cell_info::contains_parameter (const std::vector<double>& point,
const std::vector<bool>& sampled) const {
std::vector<double>::const_iterator p = point.begin();
std::vector<double>::const_iterator l = lower_left_.begin();
std::vector<double>::const_iterator u = upper_right_.begin();
std::vector<bool>::const_iterator f = sampled.begin();
for (; p < point.end(); ++p, ++f, ++l, ++u)
if (!(*f)) {
if (((*l) > (*p)) ||
((*u) < (*p)))
return false;
}
return true;
}
inline cell::cell()
: split_dimension_(0), split_point_(0.),
integral_(0.), missing_events_(0),
cell_info_(0) {}
inline cell::cell(const std::vector<double>& ll,
const std::vector<double>& ur,
const adaption_info& ainfo)
: split_dimension_(0), split_point_(0.),
integral_(0.), missing_events_(0),
cell_info_(new cell_info(ll,ur,ainfo)) {}
inline cell::cell(const std::vector<double>& ll,
const std::vector<double>& ur,
const std::vector<bool>& sampled_variables,
const adaption_info& ainfo)
: split_dimension_(0), split_point_(0.),
integral_(0.), missing_events_(0),
cell_info_(new cell_info(ll,ur,sampled_variables,ainfo)) {}
inline cell::cell(const cell& x)
: split_dimension_(x.split_dimension_),
split_point_(x.split_point_),
integral_(x.integral_),
missing_events_(x.missing_events_),
cell_info_(0) {
if (x.cell_info_)
cell_info_.reset(new cell_info(*x.cell_info_));
}
inline cell& cell::operator=(const cell& x) {
if (this == &x)
return *this;
split_dimension_ = x.split_dimension_;
split_point_ = x.split_point_;
integral_ = x.integral_;
missing_events_ = x.missing_events_;
if (x.cell_info_)
cell_info_.reset(new cell_info(*x.cell_info_));
return *this;
}
}
diff --git a/Sampling/exsample/exponential_generator.h b/Sampling/exsample/exponential_generator.h
--- a/Sampling/exsample/exponential_generator.h
+++ b/Sampling/exsample/exponential_generator.h
@@ -1,243 +1,249 @@
// -*- C++ -*-
//
// exponential_generator.h is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
//
// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
//
// ExSample is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
#ifndef EXSAMPLE_exponential_generator_h_included
#define EXSAMPLE_exponential_generator_h_included
#include "cell.h"
#include "selectors.h"
#include "statistics.h"
#include "linear_interpolator.h"
#include "binary_tree.h"
namespace exsample {
/// \brief Exception thrown, if the exponential_generator has just changed its
/// state. The attempt of generating an event should be repeated.
struct exponential_regenerate{};
/// \brief The generator for sudakov-type distributions.
template<class Function, class Random>
class exponential_generator {
public:
/// default constructor
exponential_generator()
: function_(0), check_events_(0), adaption_info_(), root_cell_(),
rnd_gen_(), did_split_(false), initialized_(false),
evolution_variable_(0), evolution_cutoff_(0.),
sample_variables_(), sample_other_variables_(),
parameter_splits_(),
last_cell_(), last_point_(), last_value_(0.),
last_parameter_bin_(), exponents_(),
last_exponent_integrand_(),
last_exponent_(), compensating_(false),
integral_accessor_(), missing_accessor_(),
parametric_selector_(), exponent_selector_(),
parametric_sampler_(), attempts_(0), accepts_(0),
- splits_(0), docompensate_(false) {}
+ splits_(0), docompensate_(false), detuning_(1.0) {}
public:
/// initialize this generator
void initialize();
/// finalize this generator
void finalize() {}
/// generate an event, returning
/// the sign of the weight or zero
/// for an event below the evolution cutoff
double generate();
/// generate an event, returning
/// the sign of the weight or zero
/// for an event below the evolution cutoff
double generate(double cutoff) {
double oldcut = evolution_cutoff_;
evolution_cutoff_ = cutoff;
double w = 0.0;
try {
w = generate();
} catch(...) {
evolution_cutoff_ = oldcut;
throw;
}
evolution_cutoff_ = oldcut;
return w;
}
/// return the last sampled phase space point
const std::vector<double>& last_point() const { return last_point_; }
/// return the last evaluated function
double last_value() const { return last_value_; }
/// indicate that the last generated point has been rejected
void reject() {
last_cell_->info().reject();
}
public:
/// return true, if this generator has been initialized
bool initialized() const { return initialized_; }
/// return true, if at least one split has been performed
bool did_split() const { return did_split_; }
/// access the function
Function& function() { return *function_; }
/// set the function
void function(Function * f) { function_ = f; }
/// access the adaption_info object
adaption_info& sampling_parameters() { return adaption_info_; }
/// indicate, if compensation should be applied
void docompensate(bool yes = true) { docompensate_ = yes; }
+ /// set the detuning parameter
+ void detuning(double val) { detuning_ = val; }
+
public:
/// put to ostream
template<class OStream>
void put(OStream& os) const;
/// get from istream
template<class IStream>
void get(IStream& is);
private:
/// check for and possibly split
/// the last selected cell
bool split();
/// get the projection of the density integrating over every
/// variable to be sampled, except the evolution variable for the
/// indicated parameter point. the k'th entry in
/// last_exponent_integrand_ is the value in the evolution
/// variable bin from evolution_splits_[k] to
/// evolution_splits_[k+1]
void get_exponent();
/// compensate
void compensate();
/// get all parameter points to build
/// all possible sub tree hashes
std::set<std::vector<double> > parameter_points();
/// get all parameter points to build
/// all possible sub tree hashes
void recursive_parameter_points(std::set<std::vector<double> >&,
std::vector<double>&,
size_t);
/// function to be sampled
Function * function_;
/// the number of events after which
/// a cell is checked for splits
unsigned long check_events_;
/// the adaption info object
adaption_info adaption_info_;
/// the root cell
binary_tree<cell> root_cell_;
/// the random number generator to be used
rnd_generator<Random> rnd_gen_;
/// wether a split has already been performed
bool did_split_;
/// wether this generator has been initialized
bool initialized_;
/// the position of the evolution variable
std::size_t evolution_variable_;
/// the cutoff on the evolution variable
double evolution_cutoff_;
/// flags of variables to be sampled
/// including the evolution variable
std::vector<bool> sample_variables_;
/// flags of variables to be sampled
/// excluding the evolution variable
std::vector<bool> sample_other_variables_;
/// the splits in any parameter done so far
/// (including the evolution variable)
std::map<std::size_t,std::vector<double> > parameter_splits_;
/// the last selected cell
binary_tree<cell>::iterator last_cell_;
/// the last sampled phasespace point
std::vector<double> last_point_;
/// the last function value
double last_value_;
/// the last parameter bin id
bit_container<parameter_hash_bits> last_parameter_bin_;
/// map parameter bin ids to exponent interpolations
std::map<bit_container<parameter_hash_bits>,linear_interpolator > exponents_;
/// the last exponent integrand
std::vector<double> last_exponent_integrand_;
/// the last exponent
std::map<bit_container<parameter_hash_bits>,linear_interpolator >::iterator last_exponent_;
/// wether or not we are compensating
bool compensating_;
/// the integral accessor to be used
integral_accessor integral_accessor_;
/// the missing events accessor to be used
parametric_missing_accessor missing_accessor_;
/// the parametric selector to be used
parametric_selector parametric_selector_;
/// the parametric selector to be used for parameter bins
parametric_selector exponent_selector_;
/// the parametric sampler to be used
parametric_sampling_selector<rnd_generator<Random> > parametric_sampler_;
/// the number of trials in the veto loo so far
unsigned long attempts_;
/// the number of accepted events so far
unsigned long accepts_;
/// number of splits done
unsigned long splits_;
/// true, if compensation should be applied
bool docompensate_;
+ /// a detuning factor to be applied to the overestimate
+ double detuning_;
+
};
}
#include "exponential_generator.icc"
#endif // EXSAMPLE_exponential_generator_h_included
diff --git a/Sampling/exsample/exponential_generator.icc b/Sampling/exsample/exponential_generator.icc
--- a/Sampling/exsample/exponential_generator.icc
+++ b/Sampling/exsample/exponential_generator.icc
@@ -1,385 +1,386 @@
// -*- C++ -*-
//
// exponential_generator.icc is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
//
// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
//
// ExSample is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
namespace exsample {
template<class Function, class Random>
void exponential_generator<Function,Random>::initialize() {
adaption_info_.dimension = function_->dimension();
adaption_info_.lower_left = function_->support().first;
adaption_info_.upper_right = function_->support().second;
if (adaption_info_.adapt.empty())
adaption_info_.adapt = std::vector<bool>(adaption_info_.dimension,true);
evolution_variable_ = function_->evolution_variable();
evolution_cutoff_ = function_->evolution_cutoff();
sample_variables_ = function_->variable_flags();
sample_other_variables_ = sample_variables_;
sample_other_variables_[evolution_variable_] = false;
last_point_.resize(adaption_info_.dimension);
parametric_selector_ = parametric_selector(&last_point_,sample_other_variables_);
exponent_selector_ = parametric_selector(&last_point_,sample_variables_);
missing_accessor_ = parametric_missing_accessor(&last_parameter_bin_);
parametric_sampler_ = parametric_sampling_selector<rnd_generator<Random> >
(&last_point_,&last_parameter_bin_,sample_other_variables_,rnd_gen_);
if (initialized_) return;
splits_ = 0;
for ( std::size_t k = 0; k < adaption_info_.dimension; ++k ) {
if ( sample_other_variables_[k] )
continue;
parameter_splits_[k].push_back(adaption_info_.lower_left[k]);
parameter_splits_[k].push_back(adaption_info_.upper_right[k]);
}
root_cell_ =
binary_tree<cell>(cell(adaption_info_.lower_left,
adaption_info_.upper_right,
sample_other_variables_,
adaption_info_));
- root_cell_.value().info().explore(rnd_gen_,adaption_info_,function_);
+ root_cell_.value().info().explore(rnd_gen_,adaption_info_,function_,detuning_);
root_cell_.value().integral(root_cell_.value().info().volume() * root_cell_.value().info().overestimate());
last_exponent_integrand_.resize(1);
check_events_ = adaption_info_.presampling_points;
initialized_ = true;
}
template<class Function, class Random>
bool exponential_generator<Function,Random>::split () {
if (adaption_info_.freeze_grid <= accepts_)
return false;
if (compensating_)
return false;
if (!(*last_cell_).info().bad(adaption_info_)) return false;
bool dosplit = false;
std::pair<std::size_t,double> sp =
(*last_cell_).info().get_split(adaption_info_,dosplit);
if (!dosplit) return false;
if (!adaption_info_.adapt[sp.first]) return false;
if (splits_ == parameter_hash_bits/2)
return false;
++splits_;
- last_cell_.node().split((*last_cell_).split(sp,rnd_gen_,function_,adaption_info_,sample_other_variables_));
+ last_cell_.node().split((*last_cell_).split(sp,rnd_gen_,function_,adaption_info_,
+ sample_other_variables_,detuning_));
if ( !sample_other_variables_[sp.first] ) {
if ( std::find(parameter_splits_[sp.first].begin(),parameter_splits_[sp.first].end(),sp.second)
== parameter_splits_[sp.first].end() ) {
parameter_splits_[sp.first].push_back(sp.second);
std::sort(parameter_splits_[sp.first].begin(),parameter_splits_[sp.first].end());
if ( sp.first == evolution_variable_ ) {
last_exponent_integrand_.push_back(0.);
}
}
}
did_split_ = true;
last_point_ = function_->parameter_point();
root_cell_.tree_accumulate(parametric_selector_,integral_accessor_,std::plus<double>());
exponents_.clear();
get_exponent();
return true;
}
template<class Function, class Random>
void exponential_generator<Function,Random>::get_exponent () {
last_parameter_bin_.reset();
root_cell_.subtree_hash (exponent_selector_,last_parameter_bin_);
last_exponent_ = exponents_.find(last_parameter_bin_);
if (last_exponent_ != exponents_.end())
return;
exponents_[last_parameter_bin_] = linear_interpolator();
last_exponent_ = exponents_.find(last_parameter_bin_);
double old_evo = last_point_[evolution_variable_];
std::vector<double>::iterator exp_it = last_exponent_integrand_.begin();
for (std::vector<double>::iterator esp = parameter_splits_[evolution_variable_].begin();
esp < boost::prior(parameter_splits_[evolution_variable_].end()); ++esp, ++exp_it) {
last_point_[evolution_variable_] = (*esp + *boost::next(esp))/2.;
*exp_it = root_cell_.accumulate(parametric_selector_,integral_accessor_,std::plus<double>());
}
exp_it = boost::prior(last_exponent_integrand_.end());
double total = 0.;
for (std::vector<double>::iterator esp = boost::prior(parameter_splits_[evolution_variable_].end());
esp > parameter_splits_[evolution_variable_].begin(); --esp, --exp_it) {
last_exponent_->second.set_interpolation(*esp,total);
total += (*exp_it) * ((*esp) - (*boost::prior(esp)));
}
last_exponent_->second.set_interpolation(parameter_splits_[evolution_variable_].front(),total);
last_point_[evolution_variable_] = old_evo;
}
template<class Function, class Random>
std::set<std::vector<double> >
exponential_generator<Function,Random>::parameter_points() {
std::set<std::vector<double> > res;
std::vector<double> pt(adaption_info_.dimension,0.);
recursive_parameter_points(res,pt,0);
return res;
}
template<class Function, class Random>
void exponential_generator<Function,Random>::
recursive_parameter_points(std::set<std::vector<double> >& res,
std::vector<double>& pt,
size_t current) {
if ( current == adaption_info_.dimension ) {
res.insert(pt);
return;
}
if ( sample_variables_[current] ) {
recursive_parameter_points(res,pt,current+1);
return;
}
for ( std::vector<double>::const_iterator sp =
parameter_splits_[current].begin();
sp != boost::prior(parameter_splits_[current].end()); ++sp ) {
pt[current] = (*sp + *boost::next(sp))/2.;
recursive_parameter_points(res,pt,current+1);
}
}
template<class Function, class Random>
void exponential_generator<Function,Random>::compensate() {
if (!did_split_ || !docompensate_) {
assert(did_split_ || last_cell_ == root_cell_.begin());
exponents_.clear();
last_cell_->info().overestimate(last_value_,last_point_);
last_cell_->integral(last_cell_->info().volume() * last_cell_->info().overestimate());
last_point_ = function_->parameter_point();
get_exponent();
return;
}
std::vector<double> themaxpoint = last_point_;
std::set<std::vector<double> > id_points
= parameter_points();
for ( std::set<std::vector<double> >::const_iterator id =
id_points.begin(); id != id_points.end(); ++id ) {
last_point_ = *id;
get_exponent();
}
std::map<bit_container<parameter_hash_bits>,linear_interpolator >
old_exponents = exponents_;
double old_oe = last_cell_->info().overestimate();
last_cell_->info().overestimate(last_value_,themaxpoint);
last_cell_->integral(last_cell_->info().volume() * last_cell_->info().overestimate());
exponents_.clear();
for ( std::set<std::vector<double> >::const_iterator id =
id_points.begin(); id != id_points.end(); ++id ) {
last_point_ = *id;
get_exponent();
std::map<bit_container<parameter_hash_bits>,linear_interpolator >::iterator
old_exp = old_exponents.find(last_parameter_bin_);
std::map<bit_container<parameter_hash_bits>,linear_interpolator >::iterator
new_exp = exponents_.find(last_parameter_bin_);
assert(old_exp != old_exponents.end() && new_exp != exponents_.end());
double old_norm = 1. - std::exp(-(old_exp->second)(adaption_info_.lower_left[evolution_variable_]));
double new_norm = 1. - std::exp(-(new_exp->second)(adaption_info_.lower_left[evolution_variable_]));
for (binary_tree<cell>::iterator it = root_cell_.begin();
it != root_cell_.end(); ++it) {
if ( !it->info().contains_parameter(last_point_,sample_variables_) )
continue;
double old_int = 0.;
double new_int = 0.;
for ( std::vector<double>::const_iterator sp = parameter_splits_[evolution_variable_].begin();
sp != boost::prior(parameter_splits_[evolution_variable_].end()); ++sp ) {
if ( *sp >= it->info().lower_left()[evolution_variable_] &&
*sp < it->info().upper_right()[evolution_variable_] ) {
double xl = *sp;
double xxl = *boost::next(sp);
double old_al =
(old_exp->second.interpolation()[xxl] - old_exp->second.interpolation()[xl]) /
(xxl-xl);
double old_bl =
(xxl * old_exp->second.interpolation()[xl] -
xl * old_exp->second.interpolation()[xxl]) /
(xxl-xl);
double new_al =
(new_exp->second.interpolation()[xxl] - new_exp->second.interpolation()[xl]) /
(xxl-xl);
double new_bl =
(xxl * new_exp->second.interpolation()[xl] -
xl * new_exp->second.interpolation()[xxl]) /
(xxl-xl);
if ( std::abs(old_al) > std::numeric_limits<double>::epsilon() ) {
old_int += (exp(-(old_al*xl+old_bl)) - exp(-(old_al*xxl+old_bl)))/old_al;
} else {
old_int += (xxl-xl)*exp(-old_bl);
}
if ( std::abs(new_al) > std::numeric_limits<double>::epsilon() ) {
new_int += (exp(-(new_al*xl+new_bl)) - exp(-(new_al*xxl+new_bl)))/new_al;
} else {
new_int += (xxl-xl)*exp(-new_bl);
}
}
}
double scaling;
if (it != last_cell_) {
if (old_int > std::numeric_limits<double>::epsilon() &&
new_int > std::numeric_limits<double>::epsilon())
scaling = ((old_norm * new_int) /
(new_norm * old_int)) - 1.;
else
scaling = 0.;
} else {
if (old_int > std::numeric_limits<double>::epsilon() &&
new_int > std::numeric_limits<double>::epsilon())
scaling = ((last_value_ * old_norm * new_int) /
(old_oe * new_norm * old_int)) - 1.;
else
scaling = 0.;
}
it->info().parametric_missing(last_parameter_bin_,
it->info().parametric_missing(last_parameter_bin_) +
static_cast<int>(round(scaling * it->info().attempted())));
if (it->info().parametric_missing(last_parameter_bin_) != 0) {
compensating_ = true;
}
}
}
last_point_ = function_->parameter_point();
}
template<class Function, class Random>
double exponential_generator<Function,Random>::generate () {
if (compensating_) {
compensating_ = false;
for (binary_tree<cell>::iterator it = root_cell_.begin();
it != root_cell_.end(); ++it)
if (it->info().parametric_compensating()) {
compensating_ = true;
break;
}
parametric_sampler_.compensate(compensating_);
}
last_point_ = function_->parameter_point();
if (last_point_[evolution_variable_] < evolution_cutoff_) {
return 0.;
}
unsigned long n_hit_miss = 0;
unsigned long n_select = 0;
double minus_log_r;
root_cell_.tree_accumulate(parametric_selector_,integral_accessor_,std::plus<double>());
get_exponent();
while (true) {
n_select = 0;
minus_log_r = -std::log(rnd_gen_()) + last_exponent_->second(last_point_[evolution_variable_]);
if (!last_exponent_->second.invertible(minus_log_r)) {
return 0.;
}
try {
last_point_[evolution_variable_] = last_exponent_->second.unique_inverse(minus_log_r);
} catch (constant_interpolation& c) {
last_point_[evolution_variable_] = rnd_gen_(c.range.first,c.range.second);
}
assert(!std::isnan(last_point_[evolution_variable_]) &&
!std::isinf(last_point_[evolution_variable_]));
if (last_point_[evolution_variable_] < evolution_cutoff_) {
return 0.;
}
++attempts_;
if (compensating_) {
root_cell_.tree_accumulate(missing_accessor_,std::plus<int>());
}
if (parameter_splits_[evolution_variable_].size() > 2)
root_cell_.tree_accumulate(parametric_selector_,integral_accessor_,std::plus<double>());
if (did_split_)
while ((last_cell_ = root_cell_.select(parametric_sampler_)) == root_cell_.end()) {
root_cell_.tree_accumulate(missing_accessor_,std::plus<int>());
if(++n_select > adaption_info_.maxtry)
throw selection_maxtry();
}
else
last_cell_ = root_cell_.begin();
last_cell_->info().select(rnd_gen_,last_point_,sample_other_variables_);
last_value_ = function_->evaluate(last_point_);
assert(last_value_ >= 0.);
last_cell_->info().selected(last_point_,last_value_,adaption_info_);
if (last_value_ > last_cell_->info().overestimate()) {
if ( std::abs(last_value_)/last_cell_->info().overestimate() > 2. ) {
last_value_ =
last_cell_->info().overestimate()*
(1.+exp(2.*(2.-std::abs(last_value_)/last_cell_->info().overestimate())));
}
compensate();
throw exponential_regenerate();
}
if (last_cell_->info().attempted() % check_events_ == 0) {
if (split()) {
throw exponential_regenerate();
}
}
if (last_value_/last_cell_->info().overestimate() > rnd_gen_()) {
function_->accept(last_point_,last_value_,last_cell_->info().overestimate());
break;
}
if ( last_value_ != 0.0 ) {
function_->veto(last_point_,last_value_,last_cell_->info().overestimate());
}
if(++n_hit_miss > adaption_info_.maxtry)
throw hit_and_miss_maxtry();
}
if (last_value_ == 0.)
return 0.;
++accepts_;
++check_events_;
last_cell_->info().accept();
return 1.;
}
template<class Function, class Random>
template<class OStream>
void exponential_generator<Function,Random>::put (OStream& os) const {
os << check_events_; ostream_traits<OStream>::separator(os);
adaption_info_.put(os);
root_cell_.put(os);
os << did_split_; ostream_traits<OStream>::separator(os);
os << initialized_; ostream_traits<OStream>::separator(os);
os << evolution_variable_; ostream_traits<OStream>::separator(os);
os << evolution_cutoff_; ostream_traits<OStream>::separator(os);
os << sample_variables_; ostream_traits<OStream>::separator(os);
os << sample_other_variables_; ostream_traits<OStream>::separator(os);
os << parameter_splits_; ostream_traits<OStream>::separator(os);
// last_cell_ is selected new so we ignore it here
os << last_point_; ostream_traits<OStream>::separator(os);
os << last_value_; ostream_traits<OStream>::separator(os);
last_parameter_bin_.put(os);
os << exponents_.size(); ostream_traits<OStream>::separator(os);
for ( std::map<bit_container<parameter_hash_bits>,linear_interpolator >::const_iterator
ex = exponents_.begin(); ex != exponents_.end() ; ++ex ) {
ex->first.put(os);
ex->second.put(os);
}
os << last_exponent_integrand_; ostream_traits<OStream>::separator(os);
os << compensating_; ostream_traits<OStream>::separator(os);
os << attempts_; ostream_traits<OStream>::separator(os);
os << accepts_; ostream_traits<OStream>::separator(os);
os << splits_; ostream_traits<OStream>::separator(os);
os << docompensate_; ostream_traits<OStream>::separator(os);
}
template<class Function, class Random>
template<class IStream>
void exponential_generator<Function,Random>::get (IStream& is) {
is >> check_events_;
adaption_info_.get(is);
root_cell_.get(is);
is >> did_split_ >> initialized_ >> evolution_variable_
>> evolution_cutoff_ >> sample_variables_ >> sample_other_variables_
>> parameter_splits_;
// last_cell_ is selected new so we ignore it here
is >> last_point_ >> last_value_;
last_parameter_bin_.get(is);
size_t dim; is >> dim;
for ( size_t k = 0; k < dim ; ++k ) {
bit_container<parameter_hash_bits> key;
key.get(is);
exponents_[key].get(is);
}
is >> last_exponent_integrand_;
last_exponent_ = exponents_.find(last_parameter_bin_);
is >> compensating_ >> attempts_ >> accepts_ >> splits_ >> docompensate_;
}
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 2:06 PM (16 h, 2 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023063
Default Alt Text
(106 KB)

Event Timeline