No OneTemporary
Actions

Size

65 KB

Subscribers

None

View Options

	diff --git a/Sampling/BinSampler.cc b/Sampling/BinSampler.cc
	--- a/Sampling/BinSampler.cc
	+++ b/Sampling/BinSampler.cc
	@@ -1,717 +1,728 @@
	// -- C++ --
	//
	// BinSampler.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 BinSampler class.
	//

	#include "BinSampler.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/Repository/Repository.h"

	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.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>

	#include "GeneralSampler.h"

	using namespace Herwig;

	BinSampler::BinSampler()
	: MultiIterationStatistics(),
	theBias(1.),
	theWeighted(false),
	theInitialPoints(1000000),
	theNIterations(1),
	theEnhancementFactor(1.0),
	theNonZeroInPresampling(false),
	theHalfPoints(false),
	theMaxNewMax(30),
	theReferenceWeight(1.0),
	theBin(-1),
	theInitialized(false),
	theRemapperPoints(0),
	theRemapChannelDimension(false),
	theLuminosityMapperBins(0),
	theGeneralMapperBins(0),
	theRemapperMinSelection(0.00001),
	theIntegrated(false),
	theRemappersFilled(false),
	- theHasGrids(false) {}
	+ theHasGrids(false),
	+ theKappa(1.){}

	BinSampler::~BinSampler() {}

	IBPtr BinSampler::clone() const {
	return new_ptr(*this);
	}

	IBPtr BinSampler::fullclone() const {
	return new_ptr(*this);
	}

	void BinSampler::sampler(Ptr<GeneralSampler>::tptr s) {
	theSampler = s;
	}

	Ptr<GeneralSampler>::tptr BinSampler::sampler() const {
	return theSampler;
	}

	string BinSampler::process() const {
	ostringstream os("");
	const StandardEventHandler& eh = *theEventHandler;
	const StandardXComb& xc = *eh.xCombs()[theBin];
	os << xc.matrixElement()->name() << " : ";
	os << xc.mePartonData()[0]->PDGName() << " "
	<< xc.mePartonData()[1]->PDGName() << " -> ";
	for ( cPDVector::const_iterator pid =
	xc.mePartonData().begin() + 2;
	pid != xc.mePartonData().end(); ++pid )
	os << (**pid).PDGName() << " ";
	return os.str();
	}

	string BinSampler::shortprocess() const {
	ostringstream os("");
	const StandardEventHandler& eh = *theEventHandler;
	const StandardXComb& xc = *eh.xCombs()[theBin];
	os << xc.mePartonData()[0]->id() << " "
	<< xc.mePartonData()[1]->id() << " : ";
	for ( cPDVector::const_iterator pid =
	xc.mePartonData().begin() + 2;
	pid != xc.mePartonData().end(); ++pid )
	os << (**pid).id() << " ";
	return os.str();
	}

	string BinSampler::id() const {
	ostringstream os("");
	const StandardEventHandler& eh = *theEventHandler;
	const StandardXComb& xc = *eh.xCombs()[theBin];
	string name = xc.matrixElement()->name();
	string::size_type i = name.find_first_of("[");
	string nameFirst = name.substr(0,i);
	i = name.find_first_of("]");
	string nameSecond = name.substr(i+1);
	os << nameFirst << nameSecond << ":";
	for ( cPDVector::const_iterator pid =
	xc.mePartonData().begin();
	pid != xc.mePartonData().end(); ++pid )
	os << (**pid).id() << (pid != (--xc.mePartonData().end()) ? "," : "");
	return os.str();
	}

	double BinSampler::evaluate(vector<double> p,
	bool remap) {
	double w = 1.0;
	if ( remap && !remappers.empty() ) {
	for ( size_t k = 0; k < p.size(); ++k ) {
	map<size_t,Remapper>::const_iterator r =
	remappers.find(k);
	if ( r != remappers.end() ) {
	pair<double,double> f = r->second.generate(p[k]);
	p[k] = f.first;
	w /= f.second;
	}
	}
	}
	try {
	w *= eventHandler()->dSigDR(p) / nanobarn;
	} catch (Veto&) {
	w = 0.0;
	} catch (...) {
	throw;
	}
	if (randomNumberString()!="")
	for ( size_t k = 0; k < p.size(); ++k ) {
	RandomNumberHistograms[RandomNumberIndex(id(),k)].first.book(p[k],w);
	RandomNumberHistograms[RandomNumberIndex(id(),k)].second+=w;
	}
	return w;
	}

	double BinSampler::generate() {
	double w = 1.;
	for ( size_t k = 0; k < lastPoint().size(); ++k ) {
	lastPoint()[k] = UseRandom::rnd();
	}
	try {
	w = evaluate(lastPoint());
	} catch (Veto&) {
	w = 0.0;
	} catch (...) {
	throw;
	}

	if ( !weighted() && initialized() ) {
	- double p = min(abs(w),referenceWeight())/referenceWeight();
	+ double p = min(abs(w),kappa()referenceWeight())/(kappa()referenceWeight());
	double sign = w >= 0. ? 1. : -1.;
	if ( p < 1 && UseRandom::rnd() > p )
	w = 0.;
	else
	- w = sign*max(abs(w),referenceWeight());
	+ w = signmax(abs(w),referenceWeight()kappa());
	}
	select(w);
	if ( w != 0.0 )
	accept();
	+ assert(kappa()==1.\|\|sampler()->almostUnweighted());
	return w;
	}

	void BinSampler::fillRemappers(bool progress) {

	if ( remappers.empty() )
	return;

	unsigned long nanPoints = 0;

	boost::progress_display* progressBar = 0;
	if ( progress ) {
	Repository::clog() << "warming up " << process();
	progressBar = new boost::progress_display(theRemapperPoints,Repository::clog());
	}

	unsigned long countzero =0;
	for ( unsigned long k = 0; k < theRemapperPoints; ++k,++countzero ) {
	if (countzero>=theRemapperPoints)break;
	double w = 1.;
	for ( size_t j = 0; j < lastPoint().size(); ++j ) {
	lastPoint()[j] = UseRandom::rnd();
	}
	try {
	w = evaluate(lastPoint(),false);
	} catch (Veto&) {
	w = 0.0;
	} catch (...) {
	throw;
	}

	if ( isnan(w) \|\| isinf(w) )
	++nanPoints;

	if ( theNonZeroInPresampling && w==0. ){
	k--;
	continue;
	}

	if ( w != 0.0 ) {
	countzero=0;
	for ( map<size_t,Remapper>::iterator r = remappers.begin();
	r != remappers.end(); ++r )
	r->second.fill(lastPoint()[r->first],w);
	}

	if ( progressBar )
	++(*progressBar);

	}

	if ( progressBar ) {
	delete progressBar;
	}

	if ( nanPoints ) {
	Repository::clog() << "Warning: " << nanPoints
	<< " out of " << theRemapperPoints << " points with nan or inf "
	<< "weight encountered while filling remappers.\n" << flush;
	}

	}

	void BinSampler::saveIntegrationData() const {

	XML::Element stats = MultiIterationStatistics::toXML();
	stats.appendAttribute("process",id());

	sampler()->grids().append(stats);

	}

	void BinSampler::readIntegrationData() {

	if ( theIntegrated )
	return;

	bool haveStats = false;

	list<XML::Element>::iterator sit = sampler()->grids().children().begin();
	for ( ; sit != sampler()->grids().children().end(); ++sit ) {
	if ( sit->type() != XML::ElementTypes::Element )
	continue;
	if ( sit->name() != "MultiIterationStatistics" )
	continue;
	string proc;
	sit->getFromAttribute("process",proc);
	if ( proc == id() ) {
	haveStats = true;
	break;
	}
	}

	if ( haveStats ) {
	MultiIterationStatistics::fromXML(*sit);
	sampler()->grids().erase(sit);
	theIntegrated = true;
	} else {
	throw Exception()
	<< "\n--------------------------------------------------------------------------------\n\n"
	<< "Expected integration data.\n\n"
	<< "* When using the build setup make sure the integrate command has been run.\n\n"
	<< "* Check the [EventGenerator].log file for further information.\n\n"
	<< "* Make sure that the Herwig folder can be found and that it contains a HerwigGrids.xml file.\n\n"
	<< "* If you have split the integration jobs, make sure that each integration job was finished.\n"
	<< " Afterwards delete the global HerwigGrids.xml file in the Herwig subfolder\n"
	<< " to automatically create an updated version of the global HerwigGrids.xml file.\n\n"
	<< "--------------------------------------------------------------------------------\n"
	<< Exception::abortnow;
	}

	}

	void BinSampler::saveRemappers() const {

	if ( remappers.empty() )
	return;

	XML::Element maps(XML::ElementTypes::Element,"Remappers");
	maps.appendAttribute("process",id());

	for ( map<size_t,Remapper>::const_iterator r = remappers.begin();
	r != remappers.end(); ++r ) {
	XML::Element rmap = r->second.toXML();
	rmap.appendAttribute("dimension",r->first);
	maps.append(rmap);
	}

	sampler()->grids().append(maps);

	}

	void BinSampler::setupRemappers(bool progress) {

	if ( !theRemapperPoints )
	return;

	if ( theRemappersFilled )
	return;

	lastPoint().resize(dimension());

	bool haveGrid = false;

	list<XML::Element>::iterator git = sampler()->grids().children().begin();
	for ( ; git != sampler()->grids().children().end(); ++git ) {
	if ( git->type() != XML::ElementTypes::Element )
	continue;
	if ( git->name() != "Remappers" )
	continue;
	string proc;
	git->getFromAttribute("process",proc);
	if ( proc == id() ) {
	haveGrid = true;
	break;
	}
	}

	if ( haveGrid ) {
	for ( list<XML::Element>::iterator cit = git->children().begin();
	cit != git->children().end(); ++cit ) {
	if ( cit->type() != XML::ElementTypes::Element )
	continue;
	if ( cit->name() != "Remapper" )
	continue;
	size_t dimension = 0;
	cit->getFromAttribute("dimension",dimension);
	remappers[dimension].fromXML(*cit);
	}
	sampler()->grids().erase(git);
	}

	if ( !haveGrid ) {

	const StandardEventHandler& eh = *eventHandler();
	const StandardXComb& xc = *eh.xCombs()[bin()];
	const pair<int,int>& pdims = xc.partonDimensions();

	set<int> remapped;

	if ( theRemapChannelDimension && xc.diagrams().size() > 1 &&
	dimension() > pdims.first + pdims.second ) {
	remappers[pdims.first] = Remapper(xc.diagrams().size(),theRemapperMinSelection,false);
	remapped.insert(pdims.first);
	}

	if ( theLuminosityMapperBins > 1 && dimension() >= pdims.first + pdims.second ) {
	for ( int n = 0; n < pdims.first; ++n ) {
	remappers[n] = Remapper(theLuminosityMapperBins,theRemapperMinSelection,true);
	remapped.insert(n);
	}
	for ( int n = dimension() - pdims.second; n < dimension(); ++n ) {
	remappers[n] = Remapper(theLuminosityMapperBins,theRemapperMinSelection,true);
	remapped.insert(n);
	}
	}

	if ( theGeneralMapperBins > 1 ) {
	for ( int n = 0; n < dimension(); n++ ) {
	if ( remapped.find(n) == remapped.end() ) {
	remappers[n] = Remapper(theGeneralMapperBins,theRemapperMinSelection,true);
	remapped.insert(n);
	}
	}
	}

	fillRemappers(progress);

	for ( map<size_t,Remapper>::iterator r = remappers.begin();
	r != remappers.end(); ++r ) {
	r->second.finalize();
	}

	}

	theRemappersFilled = true;

	}

	void BinSampler::runIteration(unsigned long points, bool progress) {

	boost::progress_display* progressBar = 0;
	if ( progress ) {
	Repository::clog() << "integrating " << process() << " , iteration "
	<< (iterations().size() + 1);
	progressBar = new boost::progress_display(points,Repository::clog());
	}

	double w=0.;
	double maxweight=0;
	int numlastmax=0;
	unsigned long countzero =0;
	int newmax=0;
	for ( unsigned long k = 0; k < points; ++k,++countzero ) {
	if (countzero>=points)break;
	w=abs(generate());

	if(theNonZeroInPresampling && w==0.0){
	k--;
	continue;
	}
	if (w!=0.0)
	countzero =0;
	numlastmax++;
	if (theHalfPoints&&maxweight<w&&
	numlastmax<(int)(points/2.)){
	if(++newmax>theMaxNewMax){
	throw Exception()
	<< "\n--------------------------------------------------------------------------------\n\n"
	<< "To many new Maxima.\n\n"
	<< "* With the option:\n\n"
	<< "* set Sampler:BinSampler:HalfPoints Yes\n\n"
	<< "* for every new maximum weight found until the half of the persampling points\n"
	<< "* the counter is set to zero. We count the number of new maxima.\n"
	<< "* You have reached: "<<newmax<<"\n"
	<< "* Did you apply reasonable cuts to the process?\n"
	<< "* You can set the maximum allowed new maxima by:"
	<< "* set Sampler:BinSampler:MaxNewMax N\n\n"
	<< "--------------------------------------------------------------------------------\n"
	<< Exception::abortnow;
	}


	maxweight=w;
	k=0;
	numlastmax=0;
	}

	if ( progress ) {
	++(*progressBar);
	}

	}

	if ( progress ) {
	Repository::clog() << "integrated ( "
	<< averageWeight() << " +/- " << sqrt(averageWeightVariance())
	<< " ) nb\nepsilon = "
	<< (abs(maxWeight()) != 0. ? averageAbsWeight()/abs(maxWeight()) : 0.);
	if ( !iterations().empty() )
	Repository::clog() << " chi2 = " << chi2();
	Repository::clog() << "\n";
	Repository::clog() << "--------------------------------------------------------------------------------\n";
	}

	if ( progressBar )
	delete progressBar;

	}

	void BinSampler::initialize(bool progress) {
	lastPoint().resize(dimension());
	if (randomNumberString()!="")
	for(size_t i=0;i<lastPoint().size();i++){
	RandomNumberHistograms[RandomNumberIndex(id(),i)] = make_pair( RandomNumberHistogram(),0.);
	}
	if ( initialized() )
	return;
	if ( !sampler()->grids().children().empty() ) {
	nIterations(1);
	}
	if ( !integrated() ) {
	unsigned long points = initialPoints();
	for ( unsigned long k = 0; k < nIterations(); ++k ) {
	runIteration(points,progress);
	if ( k < nIterations() - 1 ) {
	points = (unsigned long)(points*enhancementFactor());
	adapt();
	nextIteration();
	}
	}
	}
	isInitialized();
	}


	void BinSampler::finalize(bool){
	if (theRandomNumbers!="")
	for ( map<RandomNumberIndex,pair<RandomNumberHistogram,double> >::
	const_iterator b = RandomNumberHistograms.begin();
	b != RandomNumberHistograms.end(); ++b ) {
	b->second.first.dump(randomNumberString(), b->first.first,shortprocess(),b->first.second);
	}

	}




	BinSampler::RandomNumberHistogram::
	RandomNumberHistogram(double low,
	double up,
	unsigned int nbins)
	: lower(low) {
	nbins = nbins + 1;

	double c = up / (nbins-1.);

	for ( unsigned int k = 1; k < nbins; ++k ) {
	bins[low+c*k] = 0.;
	binsw1[low+c*k] = 0.;
	}

	}


	void BinSampler::RandomNumberHistogram::
	dump(const std::string& folder,const std::string& prefix, const std::string& process,
	const int NR) const {
	ostringstream fname("");
	std::string prefix2;
	std::string prefix3=prefix;
	std::remove_copy(prefix.begin(), prefix.end(), std::back_inserter(prefix2), '.');
	prefix3=prefix2;prefix2.clear();
	std::remove_copy(prefix3.begin(), prefix3.end(), std::back_inserter(prefix2), ':');
	prefix3=prefix2;prefix2.clear();
	std::remove_copy(prefix3.begin(), prefix3.end(), std::back_inserter(prefix2), ',');
	fname << "RN-"<< NR ;
	ofstream out((folder+"/"+prefix2+fname.str()+".dat").c_str());
	double sumofweights=0.;
	for ( map<double,double >::const_iterator b = bins.begin();b != bins.end(); ++b )
	sumofweights+=b->second;
	double sumofweights2=0.;
	for ( map<double,double >::const_iterator b = binsw1.begin();b != binsw1.end(); ++b )
	sumofweights2+=b->second;
	map<double,double >::const_iterator b2 = binsw1.begin();
	if ( sumofweights == 0 ) {
	cerr << "Not enough statistic accumulated for "
	<< process << " skipping random number diagnostic.\n"
	<< flush;
	return;
	}

	for ( map<double,double >::const_iterator b = bins.begin();
	b != bins.end(); ++b, ++b2) {
	out << " " << b->first
	<< " " << b->second/sumofweights*100.
	<< " " << b2->second/sumofweights2*100.
	<< "\n" << flush;
	}
	double xmin = -0.01;
	double xmax = 1.01;
	ofstream gpout((folder+"/"+prefix2+fname.str()+".gp").c_str());
	gpout << "set terminal epslatex color solid\n"
	<< "set output '" << prefix2+fname.str() << "-plot.tex'\n"
	<< "set xrange [" << xmin << ":" << xmax << "]\n";
	gpout << "set xlabel 'rn "<<NR <<"' \n";
	gpout << "set size 0.5,0.6\n";
	gpout << "plot '" << prefix2+fname.str()
	<< ".dat' u ($1):($2) w boxes lc rgbcolor \"blue\" t '{\\tiny "<<process <<" }',";
	gpout << " '" << prefix2+fname.str();
	gpout << ".dat' u ($1):($3) w boxes lc rgbcolor \"red\" t '';";
	gpout << "reset\n";
	}








	// If needed, insert default implementations of virtual function defined
	// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).


	void BinSampler::persistentOutput(PersistentOStream & os) const {
	MultiIterationStatistics::put(os);
	os << theBias << theWeighted << theInitialPoints << theNIterations
	<< theEnhancementFactor << theNonZeroInPresampling << theHalfPoints
	<< theMaxNewMax << theReferenceWeight
	<< theBin << theInitialized << theLastPoint
	<< theEventHandler << theSampler << theRandomNumbers
	<< theRemapperPoints << theRemapChannelDimension
	- << theLuminosityMapperBins << theGeneralMapperBins;
	+ << theLuminosityMapperBins << theGeneralMapperBins << theKappa;
	}

	void BinSampler::persistentInput(PersistentIStream & is, int) {
	MultiIterationStatistics::get(is);
	is >> theBias >> theWeighted >> theInitialPoints >> theNIterations
	>> theEnhancementFactor >> theNonZeroInPresampling >> theHalfPoints
	>> theMaxNewMax >> theReferenceWeight
	>> theBin >> theInitialized >> theLastPoint
	>> theEventHandler >> theSampler >> theRandomNumbers
	>> theRemapperPoints >> theRemapChannelDimension
	- >> theLuminosityMapperBins >> theGeneralMapperBins;
	+ >> theLuminosityMapperBins >> theGeneralMapperBins >> theKappa;
	}


	// * 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<BinSampler,MultiIterationStatistics>
	describeHerwigBinSampler("Herwig::BinSampler", "HwSampling.so");

	void BinSampler::Init() {

	static ClassDocumentation<BinSampler> documentation
	("BinSampler samples XCombs bins. This default implementation performs flat MC integration.");

	static Parameter<BinSampler,unsigned long> interfaceInitialPoints
	("InitialPoints",
	"The number of points to use for initial integration.",
	&BinSampler::theInitialPoints, 1000000, 1, 0,
	false, false, Interface::lowerlim);

	static Parameter<BinSampler,size_t> interfaceNIterations
	("NIterations",
	"The number of iterations to perform initially.",
	&BinSampler::theNIterations, 1, 1, 0,
	false, false, Interface::lowerlim);

	static Parameter<BinSampler,double> interfaceEnhancementFactor
	("EnhancementFactor",
	"The enhancement factor for the number of points in the next iteration.",
	&BinSampler::theEnhancementFactor, 2.0, 1.0, 0,
	false, false, Interface::lowerlim);

	static Switch<BinSampler,bool> interfaceNonZeroInPresampling
	("NonZeroInPresampling",
	"Switch on to count only non zero weights in presampling.",
	&BinSampler::theNonZeroInPresampling, true, false, false);
	static SwitchOption interfaceNonZeroInPresamplingYes
	(interfaceNonZeroInPresampling,
	"Yes",
	"",
	true);
	static SwitchOption interfaceNonZeroInPresamplingNo
	(interfaceNonZeroInPresampling,
	"No",
	"",
	false);

	static Switch<BinSampler,bool> interfaceHalfPoints
	("HalfPoints",
	"Switch on to reset the counter of points if new maximumis was found in the first 1/2 points.",
	&BinSampler::theHalfPoints, true, false, false);
	static SwitchOption interfaceHalfPointsYes
	(interfaceHalfPoints,
	"Yes",
	"",
	true);
	static SwitchOption interfaceHalfPointsNo
	(interfaceHalfPoints,
	"No",
	"",
	false);

	static Parameter<BinSampler,int> interfaceMaxNewMax
	("MaxNewMax",
	"The maximum number of allowed new maxima in combination with the HalfPoints option.",
	&BinSampler::theMaxNewMax, 30, 1, 0,
	false, false, Interface::lowerlim);


	static Parameter<BinSampler,string> interfaceRandomNumbers
	("RandomNumbers",
	"Prefix for distributions of the random numbers.",
	&BinSampler::theRandomNumbers, "",
	false, false);

	static Parameter<BinSampler,unsigned long> interfaceRemapperPoints
	("RemapperPoints",
	"The number of points to be used for filling remappers.",
	&BinSampler::theRemapperPoints, 10000, 0, 0,
	false, false, Interface::lowerlim);

	static Switch<BinSampler,bool> interfaceRemapChannelDimension
	("RemapChannelDimension",
	"Switch on remapping of the channel dimension.",
	&BinSampler::theRemapChannelDimension, true, false, false);
	static SwitchOption interfaceRemapChannelDimensionYes
	(interfaceRemapChannelDimension,
	"Yes",
	"",
	true);
	static SwitchOption interfaceRemapChannelDimensionNo
	(interfaceRemapChannelDimension,
	"No",
	"",
	false);

	static Parameter<BinSampler,unsigned long> interfaceLuminosityMapperBins
	("LuminosityMapperBins",
	"The number of bins to be used for remapping parton luminosities.",
	&BinSampler::theLuminosityMapperBins, 0, 0, 0,
	false, false, Interface::lowerlim);

	static Parameter<BinSampler,unsigned long> interfaceGeneralMapperBins
	("GeneralMapperBins",
	"The number of bins to be used for remapping other phase space dimensions.",
	&BinSampler::theGeneralMapperBins, 0, 0, 0,
	false, false, Interface::lowerlim);

	static Parameter<BinSampler,double> interfaceRemapperMinSelection
	("RemapperMinSelection",
	"The minimum bin selection probability for remappers.",
	&BinSampler::theRemapperMinSelection, 0.00001, 0.0, 1.0,
	false, false, Interface::limited);

	+
	+ static Parameter<BinSampler,double> interfaceKappa
	+ ("Kappa",
	+ "In AllmostUnweighted mode unweight to Kappa ReferenceWeight.",
	+ &BinSampler::theKappa, 1., 0.000001, 1.0,
	+ false, false, Interface::limited);
	+
	+
	+
	}

	diff --git a/Sampling/BinSampler.h b/Sampling/BinSampler.h
	--- a/Sampling/BinSampler.h
	+++ b/Sampling/BinSampler.h
	@@ -1,567 +1,587 @@
	// -- C++ --
	//
	// BinSampler.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_BinSampler_H
	#define Herwig_BinSampler_H
	//
	// This is the declaration of the BinSampler class.
	//

	#include "ThePEG/Handlers/StandardEventHandler.h"
	#include "ThePEG/Utilities/Exception.h"
	#include "ThePEG/Repository/UseRandom.h"

	#include "MultiIterationStatistics.h"
	#include "Remapper.h"

	namespace Herwig {

	using namespace ThePEG;

	class GeneralSampler;

	/**
	* \ingroup Matchbox
	* \author Simon Platzer
	*
	* \brief BinSampler samples XCombs bins. This default implementation
	* performs flat MC integration.
	*
	* @see \ref BinSamplerInterfaces "The interfaces"
	* defined for BinSampler.
	*/
	class BinSampler: public Herwig::MultiIterationStatistics {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	BinSampler();

	/**
	* The destructor.
	*/
	virtual ~BinSampler();
	//@}

	public:

	/**
	* Clone this object.
	*/
	Ptr<BinSampler>::ptr cloneMe() const {
	return dynamic_ptr_cast<Ptr<BinSampler>::ptr>(clone());
	}

	public:

	/**
	* Evaluate the cross section
	*/
	double evaluate(vector<double> p,
	bool remap = true);

	/**
	* Return the bias with which this sampler is selected. The sampler
	* needs to divide out this bias in its weight calculation.
	*/
	double bias() const { return theBias; }

	/**
	* Set the bias with which this sampler is selected.
	*/
	void bias(double b) { theBias = b; }

	/**
	* Set the event handler
	*/
	void eventHandler(tStdEHPtr eh) { theEventHandler = eh; }

	/**
	* Return the event handler
	*/
	tStdEHPtr eventHandler() const { return theEventHandler; }

	/**
	* Set the containing sampler
	*/
	void sampler(Ptr<GeneralSampler>::tptr);

	/**
	* Get the containing sampler
	*/
	Ptr<GeneralSampler>::tptr sampler() const;

	/**
	* Return the bin
	*/
	int bin() const { return theBin; }

	/**
	* Set the bin
	*/
	void bin(int b) { theBin = b; }

	/**
	* Return a string describing the process handled by this sampler.
	*/
	string process() const;

	/**
	* Return a short string describing the process handled by this sampler.
	*/
	string shortprocess() const;

	/**
	* Return a string identifying the process handled by this sampler.
	*/
	string id() const;

	/**
	* Return the last generated point.
	*/
	const vector<double>& lastPoint() const { return theLastPoint; }

	/**
	* Access the last generated point.
	*/
	vector<double>& lastPoint() { return theLastPoint; }

	/**
	* Return the reference weight to be used
	*/
	double referenceWeight() const { return theReferenceWeight; }

	/**
	* Set the reference weight to be used
	*/
	void referenceWeight(double w) { theReferenceWeight = w; }

	/**
	* Return true, if this sampler can provide unweighted events; if
	* the proposal density is not an overestimate, weights larger than
	* one can be generated, the handling of these points being subject
	* to the GeneralSampler class.
	*/
	virtual bool canUnweight() const { return true; }

	/**
	* Return true, if this sampler adapts on the fly while generating
	* events. Cross sections in the GeneralSampler class are calculated
	* from adding up the cross sections quoted by individual samplers.
	*/
	virtual bool adaptsOnTheFly() const { return false; }

	/**
	* If this sampler features a compensation algorithm, return true if
	* more events need to be generated to finish the compensation.
	*/
	virtual bool compensating() const { return false; }

	/**
	* Return true, if weighted events should be generated
	*/
	bool weighted() const { return theWeighted; }

	/**
	* Indicate that weighted events should be generated
	*/
	void doWeighted(bool yes = true) { theWeighted = yes; }

	/**
	* Exception to be thrown if cross section information should be updated.
	*/
	struct NextIteration {};

	/**
	* Generate the next point and return its weight; store the point in
	* lastPoint().
	*/
	virtual double generate();

	/**
	* Fill and finalize the remappers present
	*/
	void fillRemappers(bool progress);

	/**
	* Write remappers to grid file
	*/
	void saveRemappers() const;

	/**
	* Write integration data to grid files
	*/
	void saveIntegrationData() const;

	/**
	* Save grid data
	*/
	virtual void saveGrid() const {}

	/**
	* Read integration data from grid files
	*/
	void readIntegrationData();

	/**
	* Read remappers from grid file
	*/
	void setupRemappers(bool progress);

	/**
	* Run a single iteration of n points, optionally printing a
	* progress bar to cout. Calls generate n times.
	*/
	void runIteration(unsigned long n, bool progress);

	/**
	* Adapt this sampler after an iteration has been run
	*/
	virtual void adapt() {}

	/**
	* Initialize this bin sampler. This default version calls runIteration.
	*/
	virtual void initialize(bool progress);

	/**
	* Return true, if this sampler has already been initialized.
	*/
	bool initialized() const { return theInitialized; }

	/**
	* Indicate that this sampler has already been initialized.
	*/
	void isInitialized() { theInitialized = true; }

	/**
	* Return true, if integration has already been performed
	*/
	bool integrated() const { return theIntegrated; }

	/**
	* Return true, if remappers have been set up
	*/
	bool remappersFilled() const { return theRemappersFilled; }

	/**
	* Return true, if this sampler has already read grid data.
	*/
	bool hasGrids() const { return theHasGrids; }

	/**
	* Indicate that this sampler has already read grid data.
	*/
	void didReadGrids() { theHasGrids = true; }

	/**
	* Finalize this sampler.
	*/
	virtual void finalize(bool);

	/**
	* Return the total integrated cross section determined from the
	* Monte Carlo sampling so far.
	*/
	virtual CrossSection integratedXSec() const {
	return averageWeight()*nanobarn;
	}

	/**
	* Return the error on the total integrated cross section determined
	* from the Monte Carlo sampling so far.
	*/
	virtual CrossSection integratedXSecErr() const {
	return sqrt(abs(averageWeightVariance()))*nanobarn;
	}

	/**
	* Define the key for the collinear subtraction data.
	*/


	struct RandomNumberHistogram {

	/**
	* The lower bound
	*/
	double lower;

	/**
	* The bins, indexed by upper bound.
	*/
	map<double,double > bins;

	map<double,double > binsw1;
	/**
	* Constructor
	*/
	RandomNumberHistogram(double low = 0.0,
	double up = 1.,
	unsigned int nbins = 20);

	/**
	* Book an event.
	*/
	void book(double inv, double weight) {
	map<double,double>::iterator b = bins.upper_bound(inv);
	if ( b == bins.end() ) return;
	b->second = b->second+weight;
	map<double,double>::iterator b2 = binsw1.upper_bound(inv);
	if ( b2 == binsw1.end() ) return;
	b2->second = b2->second+1.;

	}

	/**
	* Write to file given name and invariant.
	*/
	void dump(const std::string& folder,const std::string& prefix, const std::string& process,const int NR)const;


	};

	typedef pair<string,size_t > RandomNumberIndex;

	map<RandomNumberIndex,pair<RandomNumberHistogram,double> > RandomNumberHistograms;

	public:

	/**
	* Return the dimension.
	*/
	int dimension() const { return theEventHandler->nDim(bin()); }

	/**
	* Return the number of points to be used for initial integration.
	*/
	unsigned long initialPoints() const { return theInitialPoints; }

	/**
	* Set the number of points to be used for initial integration.
	*/
	void initialPoints(unsigned long n) { theInitialPoints = n; }

	/**
	* Return the number of iterations to be considered for initialization.
	*/
	size_t nIterations() const { return theNIterations; }

	/**
	* Set the number of iterations to be considered for initialization.
	*/
	void nIterations(size_t n) { theNIterations = n; }

	/**
	* Set the factor to enhance the number of points for the next
	* iteration.
	*/
	void enhancementFactor(double f) { theEnhancementFactor = f; }

	/**
	* Return the factor to enhance the number of points for the next
	* iteration.
	*/
	double enhancementFactor() const { return theEnhancementFactor; }

	/**
	* Return the folder for the random number plots.
	*/
	string randomNumberString() const {return theRandomNumbers;}

	+ /**
	+ * In the AlmostUnweighted mode we do not need to unweight
	+ * the events to the reference weight.
	+ * Kappa reduces effectivly the reference weight.
	+ * This can be used for processes, where unweighting
	+ * is hardly feasable.
	+ */
	+ double kappa() const {return theKappa;}
	+
	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:

	/**
	* The bias with which this sampler is selected.
	*/
	double theBias;

	/**
	* True, if weighted events should be generated
	*/
	bool theWeighted;

	/**
	* The number of points to use for initial integration.
	*/
	unsigned long theInitialPoints;

	/**
	* The number of iterations to be considered for initialization.
	*/
	size_t theNIterations;

	/**
	* Factor to enhance the number of points for the next iteration.
	*/
	double theEnhancementFactor;

	/**
	* Switch to count only non zero weights in presampling.
	*/

	bool theNonZeroInPresampling;

	/**
	* Switch to require that we get half of the points
	* in each iteration below the maximum weight of the iteration.
	*/

	bool theHalfPoints;

	/**
	* The maximum number of allowed new maxima,
	* in combination with HalfPoints, in order to prevent unstable
	* processes.
	*/
	int theMaxNewMax;

	/**
	* The reference weight to be used
	*/
	double theReferenceWeight;

	/**
	* The bin to be sampled.
	*/
	int theBin;

	/**
	* Wether or not this sampler has already been initialized.
	*/
	bool theInitialized;

	/**
	* The last generated point.
	*/
	vector<double> theLastPoint;

	/**
	* The event handler to be used.
	*/
	tStdEHPtr theEventHandler;

	/**
	* The containing sampler
	*/
	Ptr<GeneralSampler>::tptr theSampler;

	/**
	* Folder for the random number plots.
	*/
	string theRandomNumbers;

	/**
	* Remapper objects indexed by dimension
	*/
	map<size_t,Remapper> remappers;

	/**
	* The number of points to be used for initial filling of the remappers
	*/
	unsigned long theRemapperPoints;

	/**
	* True if channels should get a remapper
	*/
	bool theRemapChannelDimension;

	/**
	* The number of bins to be used for luminosity dimensions
	*/
	unsigned long theLuminosityMapperBins;

	/**
	* The number of bins to be used for any other dimension
	*/
	unsigned long theGeneralMapperBins;

	/**
	* The minimum selection probability for remapper bins
	*/
	double theRemapperMinSelection;

	/**
	* True, if integration has already be performed
	*/
	bool theIntegrated;

	/**
	* True, if remappers have been set up
	*/
	bool theRemappersFilled;

	/**
	* True, if this sampler has already read grid data.
	*/
	bool theHasGrids;

	+
	+
	+ /**
	+ * In the AlmostUnweighted mode we do not need to unweight
	+ * the events to the reference weight.
	+ * Kappa reduces effectivly the reference weight.
	+ * This can be used for processes, where unweighting
	+ * is hardly feasable.
	+ */
	+ double theKappa;
	+
	private:

	/**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	BinSampler & operator=(const BinSampler &);

	};

	}

	#endif /* Herwig_BinSampler_H */
	diff --git a/Sampling/CellGrids/CellGridSampler.cc b/Sampling/CellGrids/CellGridSampler.cc
	--- a/Sampling/CellGrids/CellGridSampler.cc
	+++ b/Sampling/CellGrids/CellGridSampler.cc
	@@ -1,344 +1,345 @@
	// -- C++ --
	//
	// CellGridSampler.cpp 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 CellGridSampler class.
	//

	#include "CellGridSampler.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Repository/Repository.h"
	#include "ThePEG/Utilities/DescribeClass.h"

	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/ParVector.h"
	#include "ThePEG/Interface/Switch.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>

	#include "CellGridSampler.h"
	#include "Herwig/Sampling/GeneralSampler.h"

	using namespace Herwig;
	using namespace ExSample;

	CellGridSampler::CellGridSampler()
	: BinSampler(), SimpleCellGrid(),
	theExplorationPoints(1000), theExplorationSteps(8),
	theGain(0.3), theEpsilon(0.01),
	theMinimumSelection(0.0001), theLuminositySplits(0),
	theChannelSplits(0), theAllChannelSplits(false),
	theUnweightCells(true) {}

	CellGridSampler::~CellGridSampler() {}

	IBPtr CellGridSampler::clone() const {
	return new_ptr(*this);
	}

	IBPtr CellGridSampler::fullclone() const {
	return new_ptr(*this);
	}

	double CellGridSampler::generate() {
	UseRandom rnd;
	double w = SimpleCellGrid::sample(rnd,*this,lastPoint(),
	!weighted() && initialized() && theUnweightCells,
	!initialized());
	if ( !weighted() && initialized() ) {
	- double p = min(abs(w),referenceWeight())/referenceWeight();
	+ double p = min(abs(w),kappa()referenceWeight())/(kappa()referenceWeight());
	double sign = w >= 0. ? 1. : -1.;
	if ( p < 1 && UseRandom::rnd() > p )
	w = 0.;
	else
	- w = sign*max(abs(w),referenceWeight());
	+ w = signmax(abs(w),referenceWeight()kappa());
	}
	select(w);
	if ( w != 0.0 )
	accept();
	+ assert(kappa()==1.\|\|sampler()->almostUnweighted());
	return w;
	}

	void CellGridSampler::adapt() {
	UseRandom rnd;
	set<SimpleCellGrid*> newCells;
	SimpleCellGrid::adapt(theGain,theEpsilon,newCells);
	SimpleCellGrid::explore(theExplorationPoints,rnd,*this,newCells,Repository::clog());
	SimpleCellGrid::setWeights();
	SimpleCellGrid::updateIntegral();
	SimpleCellGrid::minimumSelection(theMinimumSelection);
	}

	void CellGridSampler::saveGrid() const {
	XML::Element grid = SimpleCellGrid::toXML();
	grid.appendAttribute("process",id());
	sampler()->grids().append(grid);
	}

	void CellGridSampler::initialize(bool progress) {

	bool haveGrid = false;
	list<XML::Element>::iterator git = sampler()->grids().children().begin();
	for ( ; git != sampler()->grids().children().end(); ++git ) {
	if ( git->type() != XML::ElementTypes::Element )
	continue;
	if ( git->name() != "CellGrid" )
	continue;
	string proc;
	git->getFromAttribute("process",proc);
	if ( proc == id() ) {
	haveGrid = true;
	break;
	}
	}

	if ( haveGrid ) {
	SimpleCellGrid::fromXML(*git);
	sampler()->grids().erase(git);
	didReadGrids();
	}

	lastPoint().resize(dimension());
	if (randomNumberString()!="")
	for(size_t i=0;i<lastPoint().size();i++){
	RandomNumberHistograms[RandomNumberIndex(id(),i)] = make_pair( RandomNumberHistogram(),0.);
	}

	if ( initialized() ) {
	if ( !hasGrids() )
	throw Exception() << "CellGridSampler: Require existing grid when starting to run.\n"
	<< "Did you miss setting --setupfile?"
	<< Exception::abortnow;
	return;
	}

	if ( haveGrid ) {
	if ( !integrated() )
	runIteration(initialPoints(),progress);
	isInitialized();
	return;
	}

	SimpleCellGrid::boundaries(vector<double>(dimension(),0.0),vector<double>(dimension(),1.0));
	SimpleCellGrid::weightInformation().resize(dimension());

	UseRandom rnd;

	boost::progress_display* progressBar = 0;
	if ( progress ) {
	Repository::clog() << "exploring " << process();
	progressBar = new boost::progress_display(theExplorationSteps,Repository::clog());
	}
	std::set<SimpleCellGrid*> newCells;

	if ( pre_adaption_splits().empty() &&
	(theLuminositySplits \|\| theChannelSplits \|\| theAllChannelSplits) ) {
	const StandardEventHandler& eh = *eventHandler();
	const StandardXComb& xc = *eh.xCombs()[bin()];
	the_pre_adaption_splits.resize(dimension(),0);
	const pair<int,int>& pdims = xc.partonDimensions();
	if ( theLuminositySplits && dimension() >= pdims.first + pdims.second ) {
	for ( int n = 0; n < pdims.first; ++n )
	the_pre_adaption_splits[n] = theLuminositySplits;
	for ( int n = dimension() - pdims.second; n < dimension(); ++n )
	the_pre_adaption_splits[n] = theLuminositySplits;
	}
	if ( theChannelSplits && xc.diagrams().size() &&
	dimension() > pdims.first + pdims.second ) {
	the_pre_adaption_splits[pdims.first] = theChannelSplits;
	}
	if ( theAllChannelSplits && xc.diagrams().size() > 1 &&
	dimension() > pdims.first + pdims.second ) {
	the_pre_adaption_splits[pdims.first] = xc.diagrams().size() - 1;
	}
	}

	for(int splitdim=0; splitdim<min(dimension(),(int)pre_adaption_splits().size());splitdim++)
	SimpleCellGrid::splitter(splitdim,pre_adaption_splits()[splitdim]);

	SimpleCellGrid::explore(theExplorationPoints,rnd,*this,newCells,Repository::clog());


	bool notAll = false;
	for ( std::size_t step = 1; step < theExplorationSteps; ++step ) {
	newCells.clear();
	SimpleCellGrid::adapt(theGain,theEpsilon,newCells);
	if ( progressBar )
	++(*progressBar);
	if ( newCells.empty() ) {
	notAll = true;
	break;
	}
	SimpleCellGrid::explore(theExplorationPoints,rnd,*this,newCells,Repository::clog());
	}

	if ( progressBar )
	++(*progressBar);

	SimpleCellGrid::setWeights();
	SimpleCellGrid::updateIntegral();
	SimpleCellGrid::minimumSelection(theMinimumSelection);

	if ( progressBar ) {
	if ( notAll )
	cout << "\n" << flush;
	delete progressBar;
	}

	unsigned long points = initialPoints();
	for ( unsigned long k = 0; k < nIterations(); ++k ) {
	runIteration(points,progress);
	if ( k < nIterations() - 1 ) {
	points = (unsigned long)(points*enhancementFactor());
	adapt();
	nextIteration();
	}
	}

	didReadGrids();
	isInitialized();

	}

	void CellGridSampler::finalize(bool) {
	XML::Element grid = SimpleCellGrid::toXML();
	grid.appendAttribute("process",id());
	sampler()->grids().append(grid);
	if (randomNumberString()!="")
	for ( map<RandomNumberIndex,pair<RandomNumberHistogram,double> >::
	const_iterator b = RandomNumberHistograms.begin();
	b != RandomNumberHistograms.end(); ++b ) {
	b->second.first.dump(randomNumberString(), b->first.first,shortprocess(),b->first.second);
	}


	}

	// If needed, insert default implementations of virtual function defined
	// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).

	void CellGridSampler::persistentOutput(PersistentOStream & os) const {
	os << theExplorationPoints << theExplorationSteps
	<< theGain << theEpsilon << theMinimumSelection
	<< the_pre_adaption_splits
	<< theLuminositySplits << theChannelSplits
	<< theAllChannelSplits << theUnweightCells;
	}

	void CellGridSampler::persistentInput(PersistentIStream & is, int) {
	is >> theExplorationPoints >> theExplorationSteps
	>> theGain >> theEpsilon >> theMinimumSelection
	>> the_pre_adaption_splits
	>> theLuminositySplits >> theChannelSplits
	>> theAllChannelSplits >> theUnweightCells;
	}

	// * 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<CellGridSampler,BinSampler>
	describeHerwigCellGridSampler("Herwig::CellGridSampler", "HwSampling.so");

	void CellGridSampler::Init() {

	static ClassDocumentation<CellGridSampler> documentation
	("CellGridSampler samples XCombs bins using CellGrids.");

	static Parameter<CellGridSampler,size_t> interfaceExplorationPoints
	("ExplorationPoints",
	"The number of points to use for cell exploration.",
	&CellGridSampler::theExplorationPoints, 1000, 1, 0,
	false, false, Interface::lowerlim);

	static Parameter<CellGridSampler,size_t> interfaceExplorationSteps
	("ExplorationSteps",
	"The number of exploration steps to perform.",
	&CellGridSampler::theExplorationSteps, 8, 1, 0,
	false, false, Interface::lowerlim);

	static Parameter<CellGridSampler,double> interfaceGain
	("Gain",
	"The gain factor used for adaption.",
	&CellGridSampler::theGain, 0.3, 0.0, 1.0,
	false, false, Interface::limited);

	static Parameter<CellGridSampler,double> interfaceEpsilon
	("Epsilon",
	"The efficieny threshold used for adaption.",
	&CellGridSampler::theEpsilon, 0.01, 0.0, 1.0,
	false, false, Interface::limited);

	static Parameter<CellGridSampler,double> interfaceMinimumSelection
	("MinimumSelection",
	"The minimum cell selection probability.",
	&CellGridSampler::theMinimumSelection, 0.0001, 0.0, 1.0,
	false, false, Interface::limited);

	static ParVector<CellGridSampler,int> interfacethe_pre_adaption_splits
	("preadaptionsplit",
	"The splittings for each dimension befor adaption.",
	&CellGridSampler::the_pre_adaption_splits, 1., -1, 0.0, 0.0, 0,
	false, false, Interface::lowerlim);

	static Parameter<CellGridSampler,int> interfaceLuminositySplits
	("LuminositySplits",
	"",
	&CellGridSampler::theLuminositySplits, 0, 0, 0,
	false, false, Interface::lowerlim);

	static Parameter<CellGridSampler,int> interfaceChannelSplits
	("ChannelSplits",
	"",
	&CellGridSampler::theChannelSplits, 0, 0, 0,
	false, false, Interface::lowerlim);

	static Switch<CellGridSampler,bool> interfaceAllChannelSplits
	("AllChannelSplits",
	"",
	&CellGridSampler::theAllChannelSplits, false, false, false);
	static SwitchOption interfaceAllChannelSplitsOn
	(interfaceAllChannelSplits,
	"On",
	"",
	true);
	static SwitchOption interfaceAllChannelSplitsOff
	(interfaceAllChannelSplits,
	"Off",
	"",
	false);

	static Switch<CellGridSampler,bool> interfaceUnweightCells
	("UnweightCells",
	"",
	&CellGridSampler::theUnweightCells, true, false, false);
	static SwitchOption interfaceUnweightCellsYes
	(interfaceUnweightCells,
	"Yes",
	"",
	true);
	static SwitchOption interfaceUnweightCellsNo
	(interfaceUnweightCells,
	"No",
	"",
	false);


	}

	diff --git a/Sampling/MonacoSampler.cc b/Sampling/MonacoSampler.cc
	--- a/Sampling/MonacoSampler.cc
	+++ b/Sampling/MonacoSampler.cc
	@@ -1,406 +1,399 @@
	// -- C++ --
	//
	// MonacoSampler.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 MonacoSampler class.
	//

	#include "MonacoSampler.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/Repository/Repository.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>

	#include "MonacoSampler.h"
	#include "Herwig/Sampling/GeneralSampler.h"

	using namespace Herwig;

	MonacoSampler::MonacoSampler()
	: BinSampler(),
	theAlpha(0.875),
	- theKappa(1.),
	theGridDivisions(48),
	theIterationPoints(0) {}

	MonacoSampler::~MonacoSampler() {}

	IBPtr MonacoSampler::clone() const {
	return new_ptr(*this);
	}

	IBPtr MonacoSampler::fullclone() const {
	return new_ptr(*this);
	}

	double MonacoSampler::generate() {
	double w = 1.;
	// cout<<"\npoint: ";
	std::valarray<int> upperb(dimension());
	for ( int k = 0; k < dimension(); ++k ) {
	double div = (1 - UseRandom::rnd()) * theGridDivisions;
	upperb[k] = static_cast<int>(div);
	double gupper, glower;
	if ( upperb[k] <= 0 ) {
	upperb[k] = 0;
	glower = 0.;
	gupper = theGrid(k,0);
	} else if (upperb[k] >= static_cast<int>(theGridDivisions)) {
	upperb[k] = theGridDivisions-1;
	glower = theGrid(k,theGridDivisions-2);
	gupper = theGrid(k,theGridDivisions-1);
	} else {
	glower = theGrid(k,upperb[k]-1);
	gupper = theGrid(k,upperb[k]);
	}
	double gdiff = gupper - glower;
	lastPoint()[k] = glower + (div-upperb[k])*gdiff;
	w = gdiff theGridDivisions;
	}
	// cout<<lastPoint()[k]<<" ";

	try {
	w *= eventHandler()->dSigDR(lastPoint()) / nanobarn;
	} catch (Veto&) {
	w = 0.0;
	} catch (...) {
	throw;
	}
	// only store numbers
	double wgt = w;
	if ( isnan(wgt) \|\| isinf(wgt) ) wgt = 0;
	// save results for later grid optimization
	theIterationPoints++;
	for ( int k = 0; k < dimension(); ++k ) {
	theGridData(k,upperb[k]) += wgt*wgt;
	}

	if (randomNumberString()!="")
	for ( size_t k = 0; k < lastPoint().size(); ++k ) {
	RandomNumberHistograms[RandomNumberIndex(id(),k)].first.book(lastPoint()[k],wgt);
	RandomNumberHistograms[RandomNumberIndex(id(),k)].second+=wgt;
	}

	if ( !weighted() && initialized() ) {
	- double p = min(abs(w),theKappareferenceWeight())/(theKappareferenceWeight());
	+ double p = min(abs(w),kappa()referenceWeight())/(kappa()referenceWeight());
	double sign = w >= 0. ? 1. : -1.;
	if ( p < 1 && UseRandom::rnd() > p )
	w = 0.;
	else
	- w = signmax(abs(w),theKappareferenceWeight());
	+ w = signmax(abs(w),kappa()referenceWeight());
	}
	select(w);
	- assert(theKappa==1.\|\|sampler()->almostUnweighted());
	+ assert(kappa()==1.\|\|sampler()->almostUnweighted());
	if ( w != 0.0 )
	accept();
	return w;
	}

	void MonacoSampler::saveGrid() const {
	XML::Element grid = toXML();
	grid.appendAttribute("process",id());
	sampler()->grids().append(grid);
	}

	void MonacoSampler::initialize(bool progress) {

	//read in grid
	bool haveGrid = false;
	list<XML::Element>::iterator git = sampler()->grids().children().begin();
	for ( ; git != sampler()->grids().children().end(); ++git ) {
	if ( git->type() != XML::ElementTypes::Element )
	continue;
	if ( git->name() != "Monaco" )
	continue;
	string proc;
	git->getFromAttribute("process",proc);
	if ( proc == id() ) {
	haveGrid = true;
	break;
	}
	}

	if ( haveGrid ) {
	fromXML(*git);
	sampler()->grids().erase(git);
	didReadGrids();
	} else {
	// flat grid
	theGrid.resize(dimension(),theGridDivisions);
	for (int k = 0; k < dimension(); k++)
	for (size_t l = 0; l < theGridDivisions; l++)
	theGrid(k,l) = (l+1)/static_cast<double>(theGridDivisions);
	theGridData = boost::numeric::ublas::zero_matrix<double>(dimension(),theGridDivisions);
	theIterationPoints = 0;
	}

	lastPoint().resize(dimension());
	if (randomNumberString()!="")
	for(size_t i=0;i<lastPoint().size();i++){
	RandomNumberHistograms[RandomNumberIndex(id(),i)] = make_pair( RandomNumberHistogram(),0.);
	}

	if ( initialized() ) {
	if ( !hasGrids() )
	throw Exception() << "MonacoSampler: Require existing grid when starting to run.\n"
	<< "Did you miss setting --setupfile?"
	<< Exception::abortnow;
	return;
	}

	if ( haveGrid ) {
	if ( !integrated() ) {
	runIteration(initialPoints(),progress);
	adapt();
	}
	isInitialized();
	return;
	}

	// if ( !sampler()->grids().children().empty() ) {
	// nIterations(1);
	// }
	unsigned long points = initialPoints();
	for ( unsigned long k = 0; k < nIterations(); ++k ) {
	runIteration(points,progress);
	if ( k < nIterations() - 1 ) {
	points = (unsigned long)(points*enhancementFactor());
	adapt();
	nextIteration();
	}
	}
	adapt();
	didReadGrids();
	isInitialized();
	}

	void MonacoSampler::adapt() {

	int dim = dimension();

	// refine grid
	std::valarray<double> gridcumul(dim);
	for (int k=0; k<dim; ++k) {
	double gridold = theGridData(k,0);
	double gridnew = theGridData(k,1);
	theGridData(k,0) = (gridold + gridnew) / 2.0;
	gridcumul[k] = theGridData(k,0);
	for (size_t l=1; l<theGridDivisions-1; ++l) {
	theGridData(k,l) = gridold + gridnew;
	gridold = gridnew;
	gridnew = theGridData(k,l+1);
	theGridData(k,l) = (theGridData(k,l) + gridnew) / 3.0;
	gridcumul[k] += theGridData(k,l);
	}
	theGridData(k,theGridDivisions-1) = (gridnew + gridold) / 2.0;
	gridcumul[k] += theGridData(k,theGridDivisions-1);
	}

	for (int k=0; k<dim; ++k) {
	double rc = 0.;
	std::valarray<double> ri(theGridDivisions);
	for (size_t l=0; l<theGridDivisions; ++l) {
	ri[l] = 0.;
	if ((theGridData(k,l) >= 0) && (gridcumul[k] != 0)) {
	theGridData(k,l) = max( 1.0e-30, theGridData(k,l) );
	double gpart = gridcumul[k] / theGridData(k,l);
	ri[l] = pow( (gpart - 1.0) / (gpart * log( gpart )), theAlpha);
	} else {
	ri[l] = pow( 1. / log( 1e30 ), theAlpha);
	}
	rc += ri[l];
	}
	rc /= theGridDivisions;

	double gridold = 0, gridnew = 0.;
	double deltar = 0.;
	unsigned int m = 0;
	std::valarray<double> theGridRowNew(theGridDivisions);

	for (size_t l = 0; l < theGridDivisions; ++l) {
	deltar += ri[l];
	gridold = gridnew;
	gridnew = theGrid(k,l);
	for (; deltar > rc; m++) {
	deltar -= rc;
	theGridRowNew[m] = gridnew - (gridnew - gridold) * deltar / ri[l];
	}
	}
	for (size_t l = 0; l < theGridDivisions-1; ++l) {
	theGrid(k,l) = theGridRowNew[l];
	}
	theGrid(k,theGridDivisions-1) = 1.0;
	}

	theGridData = boost::numeric::ublas::zero_matrix<double>(dimension(),theGridDivisions);
	theIterationPoints = 0;
	}

	void MonacoSampler::finalize(bool) {
	// save grid
	adapt();
	XML::Element grid = MonacoSampler::toXML();
	grid.appendAttribute("process",id());
	sampler()->grids().append(grid);

	if (randomNumberString()!="")
	for ( map<RandomNumberIndex,pair<RandomNumberHistogram,double> >::
	const_iterator b = RandomNumberHistograms.begin();
	b != RandomNumberHistograms.end(); ++b ) {
	b->second.first.dump(randomNumberString(), b->first.first,shortprocess(),b->first.second);
	}

	}

	void MonacoSampler::fromXML(const XML::Element& grid) {
	int dim = 0;
	grid.getFromAttribute("Dimension",dim);
	if ( dim != dimension() ) {
	throw std::runtime_error("[MonacoSampler] Number of dimensions in grid file does not match expectation.");
	}

	size_t griddivisions = 0;
	grid.getFromAttribute("GridDivisions",griddivisions);
	boost::numeric::ublas::matrix<double> tmpgrid(dim,griddivisions);

	pair<multimap<pair<int,string>,list<XML::Element>::iterator>::const_iterator,multimap<pair<int,string>,list<XML::Element>::iterator>::const_iterator> cit;
	cit = grid.findAll(XML::ElementTypes::Element,"GridVector");

	if ( cit.first->second == grid.children().end() )
	throw std::runtime_error("[MonacoSampler] Expected a GridVector element.");

	for (multimap<pair<int,string>,list<XML::Element>::iterator>::const_iterator iit=cit.first; iit!=cit.second; ++iit) {
	const XML::Element& gridvector = *iit->second;


	int k = 0;
	gridvector.getFromAttribute("Index",k);
	if ( k >= dim ) {
	throw std::runtime_error("[MonacoSampler] Index of grid dimension larger than grid size.");
	} else {
	list<XML::Element>::const_iterator git;
	git = gridvector.findFirst(XML::ElementTypes::ParsedCharacterData,"");
	if ( git == gridvector.children().end() )
	throw std::runtime_error("[MonacoSampler] Expected grid data.");
	istringstream bdata(git->content());

	for ( size_t l = 0; l < griddivisions; ++l ) {
	bdata >> tmpgrid(k,l);
	}
	}

	}

	// store back into main variable
	// if griddivisions do not match, rebin preserving bin density

	theGrid.resize(dim,theGridDivisions);
	theIterationPoints = 0;
	double divratio = griddivisions / static_cast<double>(theGridDivisions);
	for (int k = 0; k < dim; k++) {
	double xold = 0, xnew = 0, deltar = 0;
	size_t l = 0;
	for (size_t m = 0; m < griddivisions; m++) {
	deltar += 1;
	xold = xnew;
	xnew = tmpgrid(k,m);
	for (; deltar > divratio; l++) {
	deltar -= divratio;
	theGrid(k,l) = xnew - (xnew - xold) * deltar;
	}
	}
	theGrid(k,theGridDivisions-1) = 1.0;
	}

	theGridData = boost::numeric::ublas::zero_matrix<double>(dimension(),theGridDivisions);

	}

	XML::Element MonacoSampler::toXML() const {

	XML::Element grid(XML::ElementTypes::Element,"Monaco");
	grid.appendAttribute("Dimension",dimension());
	grid.appendAttribute("GridDivisions",theGridDivisions);

	for ( int k = 0; k < dimension(); ++k ) {
	XML::Element gridvector(XML::ElementTypes::Element,"GridVector");
	gridvector.appendAttribute("Index",k);
	ostringstream bdata;
	bdata << setprecision(17);
	for ( size_t l = 0; l < theGridDivisions; ++l )
	bdata << theGrid(k,l) << " ";

	XML::Element belem(XML::ElementTypes::ParsedCharacterData,bdata.str());
	gridvector.append(belem);
	grid.append(gridvector);
	}

	return grid;

	}

	// If needed, insert default implementations of virtual function defined
	// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).


	void MonacoSampler::persistentOutput(PersistentOStream & os) const {
	BinSampler::put(os);
	- os << theAlpha << theKappa << theGridDivisions;
	+ os << theAlpha << theGridDivisions;
	}

	void MonacoSampler::persistentInput(PersistentIStream & is, int) {
	BinSampler::get(is);
	- is >> theAlpha >> theKappa >> theGridDivisions;
	+ is >> theAlpha >> theGridDivisions;
	}


	// * 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<MonacoSampler,BinSampler>
	describeHerwigMonacoSampler("Herwig::MonacoSampler", "HwSampling.so");

	void MonacoSampler::Init() {

	static ClassDocumentation<MonacoSampler> documentation
	("MonacoSampler samples XCombs bins. This implementation performs weighted MC integration using Monaco, an adapted Vegas algorithm.");

	static Parameter<MonacoSampler,double> interfaceAlpha
	("Alpha",
	"Rate of grid modification (0 for no modification).",
	&MonacoSampler::theAlpha, 0.875, 0.0, 0,
	false, false, Interface::lowerlim);

	- static Parameter<MonacoSampler,double> interfaceKappa
	- ("Kappa",
	- "In AllmostUnweighted mode unweight to Kappa ReferenceWeight.",
	- &MonacoSampler::theKappa, 1., 0.0001, 1.0,
	- false, false, Interface::limited);
	-
	static Parameter<MonacoSampler,size_t> interfaceGridDivisions
	("GridDivisions",
	"The number of divisions per grid dimension.",
	&MonacoSampler::theGridDivisions, 48, 1, 0,
	false, false, Interface::lowerlim);

	}
	diff --git a/Sampling/MonacoSampler.h b/Sampling/MonacoSampler.h
	--- a/Sampling/MonacoSampler.h
	+++ b/Sampling/MonacoSampler.h
	@@ -1,199 +1,190 @@
	// -- C++ --
	//
	// MonacoSampler.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_MonacoSampler_H
	#define Herwig_MonacoSampler_H
	//
	// This is the declaration of the MonacoSampler class.
	//

	#include "Herwig/Sampling/BinSampler.h"

	#include "Herwig/Utilities/XML/Element.h"

	#include <boost/numeric/ublas/matrix.hpp>

	namespace Herwig {

	using namespace ThePEG;

	/**
	* \ingroup Matchbox
	* \author Michael Rauch
	*
	* \brief MonacoSampler samples XCombs bins using using Monaco
	*
	* @see \ref MonacoSamplerInterfaces "The interfaces"
	* defined for MonacoSampler.
	*/
	class MonacoSampler: public Herwig::BinSampler {

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	MonacoSampler();

	/**
	* The destructor.
	*/
	virtual ~MonacoSampler();
	//@}

	public:

	/**
	* Clone this object.
	*/
	Ptr<MonacoSampler>::ptr cloneMe() const {
	return dynamic_ptr_cast<Ptr<MonacoSampler>::ptr>(clone());
	}

	public:


	/**
	* Generate the next point and return its weight; store the point in
	* lastPoint().
	*/
	virtual double generate();

	/**
	* Adapt this sampler after an iteration has been run
	*/
	virtual void adapt();

	/**
	* Save grid data
	*/
	virtual void saveGrid() const;

	/**
	* Initialize this bin sampler. This default version calls runIteration.
	*/
	virtual void initialize(bool progress);

	/**
	* Finalize this sampler.
	*/
	virtual void finalize(bool);

	/**
	* Fill Monaco grid data from an XML element
	*/
	virtual void fromXML(const XML::Element&);

	/**
	* Return an XML element for the data of the Monaco grid
	*/
	virtual XML::Element toXML() 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).

	private:

	/**
	* Rate of grid modification (0 for no modification)
	*/
	double theAlpha;

	/**
	- * In the AlmostUnweighted mode we do not need to unweight
	- * the events to the reference weight.
	- * Kappa reduces effectivly the reference weight.
	- * This can be used for processes, where unweighting
	- * is hardly feasable.
	- */
	- double theKappa;
	-
	- /**
	* Number of grid divisions per dimension
	*/
	size_t theGridDivisions;

	/**
	* Grid boundaries
	* (first index: dimension of random numbers,
	* second index: dimension of partitions per random number)
	*/
	boost::numeric::ublas::matrix<double> theGrid;

	/**
	* Collected value per grid bin
	*/
	boost::numeric::ublas::matrix<double> theGridData;

	private:

	/**
	* Number of points collected in iteration so far
	*/
	size_t theIterationPoints;

	private:

	/**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	MonacoSampler & operator=(const MonacoSampler &);

	};

	}

	#endif /* Herwig_MonacoSampler_H */

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 4:14 PM (1 d, 15 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3796085
Default Alt Text: (65 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions