Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F8309656
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
23 KB
Subscribers
None
View Options
diff --git a/Exsample2/GeneralSampler.cc b/Exsample2/GeneralSampler.cc
--- a/Exsample2/GeneralSampler.cc
+++ b/Exsample2/GeneralSampler.cc
@@ -1,445 +1,500 @@
// -*- C++ -*-
//
// GeneralSampler.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 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 GeneralSampler class.
//
#include "GeneralSampler.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Utilities/LoopGuard.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Handlers/StandardEventHandler.h"
#include "ThePEG/Handlers/StandardXComb.h"
#include <boost/progress.hpp>
using namespace Herwig;
GeneralSampler::GeneralSampler()
: theVerbose(false), theFlatSubprocesses(false),
isSampling(false), theUpdateAfter(1),
crossSectionCalls(0), gotCrossSections(false),
theIntegratedXSec(0.), theIntegratedXSecErr(0.),
- theSumWeights(0.), theSumWeights2(0.), norm(0.) {}
+ theSumWeights(0.), theSumWeights2(0.),
+ theAttempts(0), theAccepts(0),
+ norm(0.), runCombinationData(false) {}
GeneralSampler::~GeneralSampler() {}
IBPtr GeneralSampler::clone() const {
return new_ptr(*this);
}
IBPtr GeneralSampler::fullclone() const {
return new_ptr(*this);
}
double sign(double x) {
return x >= 0. ? 1. : -1.;
}
void GeneralSampler::initialize() {
if ( theBinSampler->isUnweighting() && eventHandler()->weighted() ) {
throw InitException() << "weighted events requested from unweighting bin sampler object.";
}
if ( !theBinSampler->isUnweighting() && !eventHandler()->weighted() ) {
throw InitException() << "unweighted events requested from weighted bin sampler object.";
}
if ( theBinSampler->isUnweighting() && theFlatSubprocesses ) {
throw InitException() << "cannot run flat subprocesses with unweighted sampling.";
}
if ( !samplers.empty() )
return;
boost::progress_display* progressBar = 0;
if ( !theVerbose ) {
cout << "integrating subprocesses";
progressBar = new boost::progress_display(eventHandler()->nBins(),cout);
}
for ( int b = 0; b < eventHandler()->nBins(); ++b ) {
Ptr<BinSampler>::ptr s = theBinSampler->cloneMe();
s->eventHandler(eventHandler());
s->bin(b);
const StandardEventHandler& eh = *eventHandler();
const StandardXComb& xc = *eh.xCombs()[b];
if ( eventHandler()->weighted() && theSamplingBias ) {
s->enhanceInitialPoints(theSamplingBias->enhanceInitialPoints(xc));
s->oversamplingFactor(theSamplingBias->oversamplingFactor(xc));
}
lastSampler = s;
s->initialize(theVerbose);
samplers[b] = s;
if ( !theVerbose )
++(*progressBar);
if ( s->nanPoints() && theVerbose ) {
cout << "warning: "
<< s->nanPoints() << " of "
<< s->allPoints() << " points with nan or inf weight.\n"
<< flush;
}
}
updateCrossSections(true);
if ( samplers.empty() ) {
throw Exception() << "No processes with non-zero cross section present."
<< Exception::abortnow;
}
if ( theVerbose )
cout << "estimated total cross section is ( "
<< integratedXSec()/nanobarn << " +/- "
<< integratedXSecErr()/nanobarn << " ) nb\n" << flush;
if ( progressBar )
delete progressBar;
}
double GeneralSampler::generate() {
long tries = 0;
long excptTries = 0;
gotCrossSections = false;
lastSampler = samplers.upper_bound(UseRandom::rnd())->second;
while ( true ) {
try {
lastSampler->generate(eventHandler()->weighted());
} catch (BinSampler::NewMaximum&) {
continue;
} catch(BinSampler::UpdateCrossSections) {
updateCrossSections();
lastSampler = samplers.upper_bound(UseRandom::rnd())->second;
tries = 0;
if ( ++excptTries == eventHandler()->maxLoop() )
break;
continue;
} catch (...) {
throw;
}
if ( isnan(lastSampler->lastWeight()) || isinf(lastSampler->lastWeight()) ) {
lastSampler = samplers.upper_bound(UseRandom::rnd())->second;
tries = 0;
if ( ++excptTries == eventHandler()->maxLoop() )
break;
continue;
}
+ ++theAttempts;
+
if ( eventHandler()->weighted() && lastSampler->lastWeight() == 0.0 ) {
lastSampler->accept();
lastSampler = samplers.upper_bound(UseRandom::rnd())->second;
tries = 0;
if ( ++excptTries == eventHandler()->maxLoop() )
break;
continue;
}
- // revisit unweighting at this point later
+ // revisit unweighting weighted samplers at this point later
break;
}
+ ++theAccepts;
+
if ( excptTries == eventHandler()->maxLoop() )
throw Exception()
<< "GeneralSampler::generate() : Maximum number of tries to re-run event "
<< "selection reached. Aborting now." << Exception::runerror;
lastPoint() = lastSampler->lastPoint();
lastSampler->accept();
if ( !eventHandler()->weighted() ) {
theSumWeights += sign(lastSampler->lastWeight());
theSumWeights2 += 1.0;
return sign(lastSampler->lastWeight());
} else {
- double w = lastSampler->lastWeight()/(norm*lastSampler->bias());
+ double w = lastSampler->lastWeight()/lastSampler->bias();
theSumWeights += w;
theSumWeights2 += sqr(w);
return w;
}
return 0.;
}
void GeneralSampler::rejectLast() {
lastSampler->reject();
if ( !eventHandler()->weighted() ) {
theSumWeights -= sign(lastSampler->lastWeight());
theSumWeights2 -= 1.0;
} else {
- double w = lastSampler->lastWeight()/(norm*lastSampler->bias());
+ double w = lastSampler->lastWeight()/lastSampler->bias();
theSumWeights -= w;
theSumWeights2 -= sqr(w);
}
+ --theAttempts;
+ --theAccepts;
}
void GeneralSampler::currentCrossSections() const {
if ( !isSampling )
return;
if ( gotCrossSections )
return;
if ( crossSectionCalls > 0 ) {
if ( ++crossSectionCalls == theUpdateAfter ) {
crossSectionCalls = 0;
} else return;
}
++crossSectionCalls;
double xsec = 0.;
double var = 0.;
for ( map<double,Ptr<BinSampler>::ptr>::const_iterator s = samplers.begin();
s != samplers.end(); ++s ) {
xsec += s->second->averageWeight();
var += s->second->averageWeightVariance();
}
theIntegratedXSec = xsec;
theIntegratedXSecErr = sqrt(var);
gotCrossSections = true;
}
void GeneralSampler::updateCrossSections(bool) {
double xsec = 0.;
double var = 0.;
double sumbias = 0.;
for ( map<double,Ptr<BinSampler>::ptr>::iterator s = samplers.begin();
s != samplers.end(); ++s ) {
if ( (isSampling && s->second == lastSampler) ||
!isSampling )
s->second->nextIteration();
if ( isSampling && s->second == lastSampler ) {
s->second->maxWeight(s->second->iterations().back().maxWeight());
s->second->minWeight(s->second->iterations().back().minWeight());
}
xsec += s->second->averageWeight();
var += s->second->averageWeightVariance();
if ( !theFlatSubprocesses ) {
sumbias += s->second->averageAbsWeight() * s->second->oversamplingFactor();
} else {
sumbias += 1.;
}
}
theIntegratedXSec = xsec;
theIntegratedXSecErr = sqrt(var);
norm = sumbias;
if ( sumbias == 0.0 ) {
samplers.clear();
theIntegratedXSec = ZERO;
theIntegratedXSecErr = ZERO;
return;
}
map<double,Ptr<BinSampler>::ptr> newSamplers;
double current = 0.;
for ( map<double,Ptr<BinSampler>::ptr>::iterator s = samplers.begin();
s != samplers.end(); ++s ) {
double abssw = s->second->averageAbsWeight() * s->second->oversamplingFactor();
if ( abssw == 0.0 )
continue;
if ( !theFlatSubprocesses ) {
s->second->bias(abssw/sumbias);
current += abssw;
} else {
- s->second->bias(1.);
+ s->second->bias(1./sumbias);
current += 1.;
}
newSamplers[current/sumbias] = s->second;
}
samplers = newSamplers;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void GeneralSampler::dofinish() {
set<string> compensating;
for ( map<double,Ptr<BinSampler>::ptr>::const_iterator s =
samplers.begin(); s != samplers.end(); ++s ) {
if ( s->second->compensating() ) {
compensating.insert(s->second->process());
}
if ( s->second->nanPoints() ) {
generator()->logWarning(Exception()
<< "warning: "
<< s->second->nanPoints() << " of "
<< s->second->allPoints() << " points with nan or inf weight\n"
<< "in " << s->second->process() << Exception::warning);
if ( theVerbose ) {
cout << "warning: "
<< s->second->nanPoints() << " of "
<< s->second->allPoints() << " points with nan or inf weight\n"
<< "in " << s->second->process() << "\n";
}
}
s->second->finalize(theVerbose);
}
updateCrossSections();
if ( theVerbose ) {
if ( !compensating.empty() ) {
cout << "warning: sampling for the following processes is still compensating:\n";
for ( set<string>::const_iterator c = compensating.begin();
c != compensating.end(); ++c )
cout << *c << "\n";
}
cout << "final integrated cross section is ( "
<< integratedXSec()/nanobarn << " +/- "
<< integratedXSecErr()/nanobarn << " ) nb\n" << flush;
}
if ( !compensating.empty() ) {
generator()->logWarning(Exception()
<< "Warning: Some samplers are still in compensating mode."
<< Exception::warning);
}
+ if ( runCombinationData ) {
+
+ ofstream data("sampling.dat");
+
+ double runXSec =
+ theSumWeights/theAttempts;
+ double runXSecErr =
+ (1./theAttempts)*(1./(theAttempts-1.))*
+ abs(theSumWeights2 - sqr(theSumWeights)/theAttempts);
+
+ data << setprecision(20);
+
+ data << "CrossSectionCombined "
+ << (integratedXSec()/nanobarn) << " +/- "
+ << (integratedXSecErr()/nanobarn) << "\n"
+ << "CrossSectionRun "
+ << runXSec << " +/- " << runXSecErr << "\n"
+ << "PointsAttempted " << theAttempts << "\n"
+ << "PointsAccepted " << theAccepts << "\n"
+ << "SumWeights " << theSumWeights << "\n"
+ << "SumWeights2 " << theSumWeights2 << "\n"
+ << flush;
+
+ }
+
SamplerBase::dofinish();
+
}
void GeneralSampler::doinitrun() {
SamplerBase::doinitrun();
for ( map<double,Ptr<BinSampler>::ptr>::iterator s = samplers.begin();
s != samplers.end(); ++s ) {
s->second->eventHandler(eventHandler());
s->second->initialize(false);
assert( !s->second->iterations().empty() );
s->second->maxWeight(s->second->iterations().back().maxWeight());
s->second->minWeight(s->second->iterations().back().minWeight());
}
isSampling = true;
}
void GeneralSampler::rebind(const TranslationMap & trans) {
for ( map<double,Ptr<BinSampler>::ptr>::iterator s =
samplers.begin(); s != samplers.end(); ++s )
s->second = trans.translate(s->second);
SamplerBase::rebind(trans);
}
IVector GeneralSampler::getReferences() {
IVector ret = SamplerBase::getReferences();
for ( map<double,Ptr<BinSampler>::ptr>::iterator s =
samplers.begin(); s != samplers.end(); ++s )
ret.push_back(s->second);
return ret;
}
void GeneralSampler::persistentOutput(PersistentOStream & os) const {
os << theBinSampler << theSamplingBias << theVerbose << theFlatSubprocesses
<< samplers << lastSampler << theUpdateAfter
<< theIntegratedXSec << theIntegratedXSecErr
- << theSumWeights << theSumWeights2 << norm;
+ << theSumWeights << theSumWeights2
+ << theAttempts << theAccepts
+ << norm << runCombinationData;
}
void GeneralSampler::persistentInput(PersistentIStream & is, int) {
is >> theBinSampler >> theSamplingBias >> theVerbose >> theFlatSubprocesses
>> samplers >> lastSampler >> theUpdateAfter
>> theIntegratedXSec >> theIntegratedXSecErr
- >> theSumWeights >> theSumWeights2 >> norm;
+ >> theSumWeights >> theSumWeights2
+ >> theAttempts >> theAccepts
+ >> norm >> runCombinationData;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<GeneralSampler,SamplerBase>
describeHerwigGeneralSampler("Herwig::GeneralSampler", "HwExsample2.so");
void GeneralSampler::Init() {
static ClassDocumentation<GeneralSampler> documentation
("A GeneralSampler class");
static Reference<GeneralSampler,BinSampler> interfaceBinSampler
("BinSampler",
"The bin sampler to be used.",
&GeneralSampler::theBinSampler, false, false, true, false, false);
static Reference<GeneralSampler,SamplingBias> interfaceSamplingBias
("SamplingBias",
"Set the sampling bias to apply.",
&GeneralSampler::theSamplingBias, false, false, true, true, false);
static Parameter<GeneralSampler,size_t> interfaceUpdateAfter
("UpdateAfter",
"Update cross sections every number of events.",
&GeneralSampler::theUpdateAfter, 1, 1, 0,
false, false, Interface::lowerlim);
static Switch<GeneralSampler,bool> interfaceVerbose
("Verbose",
"",
&GeneralSampler::theVerbose, false, false, false);
static SwitchOption interfaceVerboseOn
(interfaceVerbose,
"On",
"",
true);
static SwitchOption interfaceVerboseOff
(interfaceVerbose,
"Off",
"",
false);
static Switch<GeneralSampler,bool> interfaceFlatSubprocesses
("FlatSubprocesses",
"[debug] ",
&GeneralSampler::theFlatSubprocesses, false, false, false);
static SwitchOption interfaceFlatSubprocessesOn
(interfaceFlatSubprocesses,
"On",
"",
true);
static SwitchOption interfaceFlatSubprocessesOff
(interfaceFlatSubprocesses,
"Off",
"",
false);
interfaceFlatSubprocesses.rank(-1);
+ static Switch<GeneralSampler,bool> interfaceRunCombinationData
+ ("RunCombinationData",
+ "",
+ &GeneralSampler::runCombinationData, false, false, false);
+ static SwitchOption interfaceRunCombinationDataOn
+ (interfaceRunCombinationData,
+ "On",
+ "",
+ true);
+ static SwitchOption interfaceRunCombinationDataOff
+ (interfaceRunCombinationData,
+ "Off",
+ "",
+ false);
+
+ interfaceRunCombinationData.rank(-1);
+
}
diff --git a/Exsample2/GeneralSampler.h b/Exsample2/GeneralSampler.h
--- a/Exsample2/GeneralSampler.h
+++ b/Exsample2/GeneralSampler.h
@@ -1,313 +1,344 @@
// -*- C++ -*-
//
// GeneralSampler.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 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_GeneralSampler_H
#define Herwig_GeneralSampler_H
//
// This is the declaration of the GeneralSampler class.
//
#include "ThePEG/Handlers/SamplerBase.h"
#include "BinSampler.h"
#include "SamplingBias.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief A GeneralSampler class
*
* @see \ref GeneralSamplerInterfaces "The interfaces"
* defined for GeneralSampler.
*/
class GeneralSampler: public SamplerBase {
public:
/**
* Exception thrown in case too many attempts to unweight an event.
*/
struct MaxTryException : public ThePEG::Exception {};
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
GeneralSampler();
/**
* The destructor.
*/
virtual ~GeneralSampler();
//@}
public:
/** @name Virtual functions from SamplerBase. */
//@{
/**
* Initialize the the sampler, possibly doing presampling of the
* phase space.
*/
virtual void initialize();
/**
* Generarate a new phase space point and return a weight associated
* with it. This weight should preferably be 1.
*/
virtual double generate();
/**
* Reject the last chosen phase space point.
*/
virtual void rejectLast();
/**
* If the sampler is able to sample several different functions
* separately, this function should return the last chosen
* function. This default version always returns 0.
*/
virtual int lastBin() const { return lastSampler ? lastSampler->bin() : 0; }
/**
* Return the total integrated cross section determined from the
* Monte Carlo sampling so far.
*/
virtual CrossSection integratedXSec() const {
currentCrossSections();
return theIntegratedXSec * nanobarn;
}
/**
* Return the error on the total integrated cross section determined
* from the Monte Carlo sampling so far.
*/
virtual CrossSection integratedXSecErr() const {
currentCrossSections();
return theIntegratedXSecErr * nanobarn;
}
/**
* Return the overestimated integrated cross section.
*/
virtual CrossSection maxXSec() const {
currentCrossSections();
return abs(theIntegratedXSec) * nanobarn;
}
/**
* Return the sum of the weights returned by generate() so far (of
* the events that were not rejeted).
*/
virtual double sumWeights() const {
return theSumWeights;
}
/**
* Return the sum of the weights squaredreturned by generate() so far (of
* the events that were not rejeted).
*/
virtual double sumWeights2() const {
return theSumWeights2;
}
+
+ /**
+ * Return the number of attempts
+ */
+ unsigned long attempts() const {
+ return theAttempts;
+ }
+
+ /**
+ * Return the number of accepts
+ */
+ unsigned long accepts() const {
+ return theAccepts;
+ }
+
//@}
protected:
/**
* Calculate cross sections from samplers, assuming the start of a
* new iteration.
*/
void updateCrossSections(bool firstTime = false);
/**
* Calculate cross sections from samplers at current state.
*/
void currentCrossSections() const;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
protected:
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
/**
* Finalize this object. Called in the run phase just after a
* run has ended. Used eg. to write out statistics.
*/
virtual void dofinish();
/**
* Rebind pointer to other Interfaced objects. Called in the setup phase
* after all objects used in an EventGenerator has been cloned so that
* the pointers will refer to the cloned objects afterwards.
* @param trans a TranslationMap relating the original objects to
* their respective clones.
* @throws RebindException if no cloned object was found for a given
* pointer.
*/
virtual void rebind(const TranslationMap & trans);
/**
* Return a vector of all pointers to Interfaced objects used in this
* object.
* @return a vector of pointers.
*/
virtual IVector getReferences();
private:
/**
* The bin sampler to use.
*/
Ptr<BinSampler>::ptr theBinSampler;
/**
* The sampling bias to apply
*/
Ptr<SamplingBias>::ptr theSamplingBias;
/**
* Whether or not additional information should be printed to cout.
*/
bool theVerbose;
/**
* True, if subprocesses should be selected flat. This is a debug
* flag, cross section information and distributions will not be
* correct.
*/
bool theFlatSubprocesses;
/**
* True, if we are generating events.
*/
bool isSampling;
/**
* The selector map for the bin samplers.
*/
map<double,Ptr<BinSampler>::ptr> samplers;
/**
* The last selected bin sampler.
*/
Ptr<BinSampler>::tptr lastSampler;
/**
* The number of events after which cross sections should truly be
* updated. This is used to prevent exhaustive combination of
* statistics when HepMC events are written out.
*/
size_t theUpdateAfter;
/**
* The number of calls to currentCrossSections since the last
* update.
*/
mutable size_t crossSectionCalls;
/**
* True, if currentCrossSections has been called since the last call
* to generate.
*/
mutable bool gotCrossSections;
/**
* The integrated cross section in nanobarn
*/
mutable double theIntegratedXSec;
/**
* The integrated cross section error in nanobarn
*/
mutable double theIntegratedXSecErr;
/**
* The sum of weights
*/
double theSumWeights;
/**
* The sum of weights squared
*/
double theSumWeights2;
/**
+ * The number of attempts
+ */
+ unsigned long theAttempts;
+
+ /**
+ * The number of accepts
+ */
+ unsigned long theAccepts;
+
+ /**
* The sum of absolute cross section.
*/
double norm;
+ /**
+ * True, if information for combining unnormalized runs should be
+ * printed out
+ */
+ bool runCombinationData;
+
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
GeneralSampler & operator=(const GeneralSampler &);
};
}
#endif /* Herwig_GeneralSampler_H */
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Sat, Dec 21, 4:06 PM (1 d, 21 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023381
Default Alt Text
(23 KB)
Attached To
rHERWIGHG herwighg
Event Timeline
Log In to Comment