No OneTemporary
Actions

Size

19 KB

Subscribers

None

View Options

	diff --git a/Analysis/NLORivetAnalysis.cc b/Analysis/NLORivetAnalysis.cc
	--- a/Analysis/NLORivetAnalysis.cc
	+++ b/Analysis/NLORivetAnalysis.cc
	@@ -1,304 +1,309 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the NLORivetAnalysis class.
	//

	#include "NLORivetAnalysis.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/ParVector.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Vectors/HepMCConverter.h"
	#include "ThePEG/Config/HepMCHelper.h"
	#include "ThePEG/EventRecord/Event.h"
	#include "ThePEG/EventRecord/SubProcess.h"
	#include "ThePEG/EventRecord/SubProcessGroup.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Repository/CurrentGenerator.h"
	#include "HepMC/GenEvent.h"
	#include "Rivet/AnalysisHandler.hh"
	#include "Rivet/Tools/Logging.hh"
	#include "HepMC/IO_AsciiParticles.h"
	#include "HepMC/IO_GenEvent.h"
	#include <config.h>

	using namespace ThePEG;

	NLORivetAnalysis::NLORivetAnalysis()
	: _remnantId(82), _format(1),_unitchoice(),
	_geneventPrecision(16), debug(false), _rivet(), _nevent(0) {}

	HepMC::GenEvent * NLORivetAnalysis::makeEvent(tEventPtr event, tSubProPtr sub, long no,
	Energy eUnit, Length lUnit,
	CrossSection xsec, CrossSection xsecErr) const {

	// generate beam particles
	const PPair& beam = event->incoming();
	HepMC::GenParticle * b1 =
	HepMCTraits<HepMC::GenEvent>::newParticle(beam.first->momentum(),beam.first->id(),
	1,eUnit);
	HepMC::GenParticle * b2 =
	HepMCTraits<HepMC::GenEvent>::newParticle(beam.second->momentum(),beam.second->id(),
	1,eUnit);

	// generate remnants
	HepMC::GenParticle * r1 =
	HepMCTraits<HepMC::GenEvent>::newParticle(beam.first->momentum() -
	sub->incoming().first->momentum(),
	_remnantId,1,eUnit);
	HepMC::GenParticle * r2 =
	HepMCTraits<HepMC::GenEvent>::newParticle(beam.second->momentum() -
	sub->incoming().second->momentum(),
	_remnantId,1,eUnit);

	// generate outgoing particles
	vector<HepMC::GenParticle*> outgoing;
	for ( ParticleVector::const_iterator p = sub->outgoing().begin();
	p != sub->outgoing().end(); ++p ) {
	outgoing.push_back(HepMCTraits<HepMC::GenEvent>::newParticle((p).momentum(),(p).id(),
	1,eUnit));
	}

	// generate one blob vertex
	HepMC::GenVertex * vertex = HepMCTraits<HepMC::GenEvent>::newVertex();

	HepMCTraits<HepMC::GenEvent>::addIncoming(*vertex,b1);
	HepMCTraits<HepMC::GenEvent>::addIncoming(*vertex,b2);

	HepMCTraits<HepMC::GenEvent>::addOutgoing(*vertex,r1);
	HepMCTraits<HepMC::GenEvent>::addOutgoing(*vertex,r2);

	for ( vector<HepMC::GenParticle*>::const_iterator p = outgoing.begin();
	p != outgoing.end(); ++p )
	HepMCTraits<HepMC::GenEvent>::addOutgoing(vertex,p);

	HepMC::GenEvent * ev =
	HepMCTraits<HepMC::GenEvent>::newEvent(no,event->weight()*sub->groupWeight(),
	event->optionalWeights());
	HepMCTraits<HepMC::GenEvent>::setUnits(*ev,eUnit,lUnit);
	HepMCTraits<HepMC::GenEvent>::setBeamParticles(*ev,b1,b2);

	HepMCTraits<HepMC::GenEvent>::addVertex(*ev,vertex);

	HepMCTraits<HepMC::GenEvent>::setCrossSection(*ev,xsec/picobarn,
	xsecErr/picobarn);

	return ev;

	}

	void NLORivetAnalysis::analyze(ThePEG::tEventPtr event, long ieve, int loop, int state) {
	Energy eUnit;
	Length lUnit;
	switch (_unitchoice) {
	default: eUnit = GeV; lUnit = millimeter; break;
	case 1: eUnit = MeV; lUnit = millimeter; break;
	case 2: eUnit = GeV; lUnit = centimeter; break;
	case 3: eUnit = MeV; lUnit = centimeter; break;
	}

	tcEHPtr eh = dynamic_ptr_cast<tcEHPtr>(event->primaryCollision()->handler());
	assert(eh);

	CrossSection xsec = eh->integratedXSec();
	CrossSection xsecErr = eh->integratedXSecErr();

	tSubProPtr sub = event->primarySubProcess();
	Ptr<SubProcessGroup>::tptr grp =
	dynamic_ptr_cast<Ptr<SubProcessGroup>::tptr>(sub);

	AnalysisHandler::analyze(event, ieve, loop, state);
	// Rotate to CMS, extract final state particles and call analyze(particles).
	// convert to hepmc

	HepMC::GenEvent * hepmc =
	makeEvent(event,sub,_nevent,eUnit,lUnit,xsec,xsecErr);


	CurrentGenerator::Redirect stdout(cout);
	if(_rivet){
	#if ThePEG_RIVET_VERSION == 1
	_rivet->analyze(*hepmc);
	#elif ThePEG_RIVET_VERSION > 1
	try {
	_rivet->analyze(*hepmc);
	} catch (const YODA::Exception & e) {
	Throw<Exception>() << "Warning: Rivet/Yoda got the exception: "<< e.what()<<"\n"
	<< Exception::warning;
	}
	#else
	#error "Unknown ThePEG_RIVET_VERSION"
	#endif
	}
	// delete hepmc event
	delete hepmc;

	if ( grp ) {

	for ( SubProcessVector::const_iterator s = grp->dependent().begin();
	s != grp->dependent().end(); ++s ) {

	hepmc = makeEvent(event,*s,_nevent,eUnit,lUnit,xsec,xsecErr);

	if ( _rivet ){
	#if ThePEG_RIVET_VERSION == 1
	_rivet->analyze(*hepmc);
	#elif ThePEG_RIVET_VERSION > 1
	try {
	_rivet->analyze(*hepmc);
	} catch (const YODA::Exception & e) {
	Throw<Exception>() << "Warning: Rivet/Yoda got the exception: "<< e.what()<<"\n"
	<< Exception::warning;
	}
	#else
	#error "Unknown ThePEG_RIVET_VERSION"
	#endif
	}
	// delete hepmc event
	delete hepmc;


	}

	}

	++_nevent;

	}


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

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

	void NLORivetAnalysis::persistentOutput(ThePEG::PersistentOStream & os) const {
	os << _analyses << filename << debug;
	}

	void NLORivetAnalysis::persistentInput(ThePEG::PersistentIStream & is, int) {
	is >> _analyses >> filename >> debug;
	}

	ThePEG::ClassDescription<NLORivetAnalysis> NLORivetAnalysis::initNLORivetAnalysis;
	// Definition of the static class description member.

	void NLORivetAnalysis::Init() {

	static ThePEG::ClassDocumentation<NLORivetAnalysis> documentation
	("The NLORivetAnalysis class is a simple class to allow analyses"
	" from the Rivet library to be called from ThePEG");

	static ThePEG::ParVector<NLORivetAnalysis,string> interfaceAnalyses
	("Analyses",
	"The names of the Rivet analyses to use",
	&NLORivetAnalysis::_analyses, -1, "", "","" "",
	false, false, ThePEG::Interface::nolimits);

	static Parameter<NLORivetAnalysis,long> interfaceRemnantId
	("RemnantId",
	"Set the PDG id to be used for remnants.",
	&NLORivetAnalysis::_remnantId, 82, 0, 0,
	false, false, Interface::nolimits);

	static Parameter<NLORivetAnalysis,string> interfaceFilename
	("Filename",
	#if ThePEG_RIVET_VERSION == 1
	"The name of the file where the AIDA histograms are put. If empty, "
	"the run name will be used instead. '.aida' will in any case be "
	"appended to the file name.",
	#elif ThePEG_RIVET_VERSION > 1
	"The name of the file where the YODA histograms are put. If empty, "
	"the run name will be used instead. '.yoda' will in any case be "
	"appended to the file name.",
	#else
	#error "Unknown ThePEG_RIVET_VERSION"
	#endif
	&NLORivetAnalysis::filename, "", true, false);


	static Switch<NLORivetAnalysis,bool> interfaceDebug
	("Debug",
	"Enable debug information from Rivet",
	&NLORivetAnalysis::debug, false, true, false);
	static SwitchOption interfaceDebugNo
	(interfaceDebug,
	"No",
	"Disable debug information.",
	false);
	static SwitchOption interfaceDebugYes
	(interfaceDebug,
	"Yes",
	"Enable debug information from Rivet.",
	true);


	interfaceAnalyses.rank(10);

	}

	void NLORivetAnalysis::dofinish() {
	AnalysisHandler::dofinish();
	if( _nevent > 0 && _rivet ) {
	CurrentGenerator::Redirect stdout(cout);
	+#if ThePEG_RIVET_VERSION > 2
	+ _rivet->setCrossSection(generator()->integratedXSec()/picobarn,
	+ generator()->integratedXSecErr()/picobarn);
	+#else
	_rivet->setCrossSection(generator()->integratedXSec()/picobarn);
	+#endif
	_rivet->finalize();

	string fname = filename;
	#if ThePEG_RIVET_VERSION == 1
	if ( fname.empty() )
	fname = generator()->path() + "/" + generator()->runName() + ".aida";
	#elif ThePEG_RIVET_VERSION > 1
	if ( fname.empty() )
	fname = generator()->path() + "/" + generator()->runName() + ".yoda";
	#else
	#error "Unknown ThePEG_RIVET_VERSION"
	#endif
	_rivet->writeData(fname);
	}
	delete _rivet;
	- _rivet = 0;
	+ _rivet = nullptr;
	}

	void NLORivetAnalysis::doinit() {
	AnalysisHandler::doinit();
	if(_analyses.empty())
	throw ThePEG::Exception() << "Must have at least one analysis loaded in "
	<< "NLORivetAnalysis::doinitrun()"
	<< ThePEG::Exception::runerror;

	// check that analysis list is available
	_rivet = new Rivet::AnalysisHandler; //(fname);
	_rivet->addAnalyses(_analyses);
	if ( _rivet->analysisNames().size() != _analyses.size() ) {
	throw ThePEG::Exception()
	<< "Rivet could not find all requested analyses.\n"
	<< "Use 'rivet --list-analyses' to check availability.\n"
	<< ThePEG::Exception::runerror;
	}
	delete _rivet;
	_rivet = 0;
	}

	void NLORivetAnalysis::doinitrun() {
	AnalysisHandler::doinitrun();
	// create NLORivet analysis handler
	CurrentGenerator::Redirect stdout(cout);
	_rivet = new Rivet::AnalysisHandler; //(fname);
	_rivet->addAnalyses(_analyses);
	// check that analysis list is still available
	if ( _rivet->analysisNames().size() != _analyses.size() ) {
	throw ThePEG::Exception()
	<< "Rivet could not find all requested analyses.\n"
	<< "Use 'rivet --list-analyses' to check availability.\n"
	<< ThePEG::Exception::runerror;
	}
	if ( debug )
	Rivet::Log::setLevel("Rivet",Rivet::Log::DEBUG);
	}
	diff --git a/Analysis/RivetAnalysis.cc b/Analysis/RivetAnalysis.cc
	--- a/Analysis/RivetAnalysis.cc
	+++ b/Analysis/RivetAnalysis.cc
	@@ -1,189 +1,194 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the RivetAnalysis class.
	//

	#include "ThePEG/Interface/Interfaced.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/ParVector.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Vectors/HepMCConverter.h"
	#include "ThePEG/Config/HepMCHelper.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Repository/CurrentGenerator.h"
	#include "HepMC/GenEvent.h"
	#include "RivetAnalysis.h"
	#include "Rivet/AnalysisHandler.hh"
	#include "Rivet/Tools/RivetPaths.hh"
	#include "Rivet/Tools/Logging.hh"
	#include <config.h>

	using namespace ThePEG;

	RivetAnalysis::RivetAnalysis() : debug(false), _rivet(), _nevent(0)
	{}

	void RivetAnalysis::analyze(ThePEG::tEventPtr event, long ieve, int loop, int state) {
	++_nevent;
	AnalysisHandler::analyze(event, ieve, loop, state);
	// Rotate to CMS, extract final state particles and call analyze(particles).
	// convert to hepmc
	HepMC::GenEvent * hepmc = ThePEG::HepMCConverter<HepMC::GenEvent>::convert(*event);
	// analyse the event
	CurrentGenerator::Redirect stdout(cout);
	if ( _rivet ){
	#if ThePEG_RIVET_VERSION == 1
	_rivet->analyze(*hepmc);
	#elif ThePEG_RIVET_VERSION > 1
	try {
	_rivet->analyze(*hepmc);
	} catch (const YODA::Exception & e) {
	Throw<Exception>() << "Warning: Rivet/Yoda got the exception: "<< e.what()<<"\n"
	<< Exception::warning;
	}
	#else
	#error "Unknown ThePEG_RIVET_VERSION"
	#endif
	}
	// delete hepmc event
	delete hepmc;
	}

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

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

	void RivetAnalysis::persistentOutput(ThePEG::PersistentOStream & os) const {
	os << _analyses << _paths << filename << debug;
	}

	void RivetAnalysis::persistentInput(ThePEG::PersistentIStream & is, int) {
	is >> _analyses >> _paths >> filename >> debug;
	}

	ThePEG::ClassDescription<RivetAnalysis> RivetAnalysis::initRivetAnalysis;
	// Definition of the static class description member.

	void RivetAnalysis::Init() {

	static ThePEG::ClassDocumentation<RivetAnalysis> documentation
	("The RivetAnalysis class is a simple class to allow analyses"
	" from the Rivet library to be called from ThePEG");

	static ThePEG::ParVector<RivetAnalysis,string> interfaceAnalyses
	("Analyses",
	"The names of the Rivet analyses to use",
	&RivetAnalysis::_analyses, -1, "", "","" "",
	false, false, ThePEG::Interface::nolimits);

	static ThePEG::ParVector<RivetAnalysis,string> interfacePaths
	("Paths",
	"The directory paths where Rivet should look for analyses.",
	&RivetAnalysis::_paths, -1, "", "","" "",
	false, false, ThePEG::Interface::nolimits);

	static Parameter<RivetAnalysis,string> interfaceFilename
	("Filename",
	#if ThePEG_RIVET_VERSION == 1
	"The name of the file where the AIDA histograms are put. If empty, "
	"the run name will be used instead. '.aida' will in any case be "
	"appended to the file name.",
	#elif ThePEG_RIVET_VERSION > 1
	"The name of the file where the YODA histograms are put. If empty, "
	"the run name will be used instead. '.yoda' will in any case be "
	"appended to the file name.",
	#else
	#error "Unknown ThePEG_RIVET_VERSION"
	#endif
	&RivetAnalysis::filename, "", true, false);


	static Switch<RivetAnalysis,bool> interfaceDebug
	("Debug",
	"Enable debug information from Rivet",
	&RivetAnalysis::debug, false, true, false);
	static SwitchOption interfaceDebugNo
	(interfaceDebug,
	"No",
	"Disable debug information.",
	false);
	static SwitchOption interfaceDebugYes
	(interfaceDebug,
	"Yes",
	"Enable debug information from Rivet.",
	true);


	interfaceAnalyses.rank(10);

	}

	void RivetAnalysis::dofinish() {
	AnalysisHandler::dofinish();
	if( _nevent > 0 && _rivet ) {
	CurrentGenerator::Redirect stdout(cout);
	+#if ThePEG_RIVET_VERSION > 2
	+ _rivet->setCrossSection(generator()->integratedXSec()/picobarn,
	+ generator()->integratedXSecErr()/picobarn);
	+#else
	_rivet->setCrossSection(generator()->integratedXSec()/picobarn);
	+#endif
	_rivet->finalize();

	string fname = filename;
	#if ThePEG_RIVET_VERSION == 1
	if ( fname.empty() )
	fname = generator()->path() + "/" + generator()->runName() + ".aida";
	#elif ThePEG_RIVET_VERSION > 1
	if ( fname.empty() )
	fname = generator()->path() + "/" + generator()->runName() + ".yoda";
	#else
	#error "Unknown ThePEG_RIVET_VERSION"
	#endif
	_rivet->writeData(fname);
	}
	delete _rivet;
	- _rivet = 0;
	+ _rivet = nullptr;
	}

	void RivetAnalysis::doinit() {
	AnalysisHandler::doinit();
	if(_analyses.empty())
	throw ThePEG::Exception() << "Must have at least one analysis loaded in "
	<< "RivetAnalysis::doinitrun()"
	<< ThePEG::Exception::runerror;

	// check that analysis list is available
	_rivet = new Rivet::AnalysisHandler; //(fname);
	for ( int i = 0, N = _paths.size(); i < N; ++i ) Rivet::addAnalysisLibPath(_paths[i]);
	_rivet->addAnalyses(_analyses);
	if ( _rivet->analysisNames().size() != _analyses.size() ) {
	throw ThePEG::Exception()
	<< "Rivet could not find all requested analyses.\n"
	<< "Use 'rivet --list-analyses' to check availability.\n"
	<< ThePEG::Exception::runerror;
	}
	delete _rivet;
	_rivet = 0;
	}

	void RivetAnalysis::doinitrun() {
	AnalysisHandler::doinitrun();
	// create Rivet analysis handler
	CurrentGenerator::Redirect stdout(cout);
	_rivet = new Rivet::AnalysisHandler; //(fname);
	for ( int i = 0, N = _paths.size(); i < N; ++i ) Rivet::addAnalysisLibPath(_paths[i]);
	_rivet->addAnalyses(_analyses);
	// check that analysis list is still available
	if ( _rivet->analysisNames().size() != _analyses.size() ) {
	throw ThePEG::Exception()
	<< "Rivet could not find all requested analyses.\n"
	<< "Use 'rivet --list-analyses' to check availability.\n"
	<< ThePEG::Exception::runerror;
	}
	if ( debug )
	Rivet::Log::setLevel("Rivet",Rivet::Log::DEBUG);
	}
	diff --git a/m4/rivet.m4 b/m4/rivet.m4
	--- a/m4/rivet.m4
	+++ b/m4/rivet.m4
	@@ -1,79 +1,81 @@
	dnl ##### RIVET #####
	AC_DEFUN([THEPEG_CHECK_RIVET],
	[
	AC_REQUIRE([THEPEG_CHECK_HEPMC])
	AC_REQUIRE([THEPEG_CHECK_GSL])
	AC_MSG_CHECKING([for Rivet location])
	RIVETINCLUDE=""
	LOAD_RIVET=""
	RIVETLIBS="-lRivet"

	AC_ARG_WITH(rivet,
	AC_HELP_STRING([--with-rivet=DIR],[Location of Rivet installation @<:@default=system libs@:>@]),
	[],
	[with_rivet=system])

	if test "x$with_hepmc" = "xno"; then
	with_rivet=no
	fi

	if test "x$with_rivet" = "xno"; then
	AC_MSG_RESULT([Rivet support disabled.])
	elif test "x$with_rivet" = "xsystem"; then
	AC_MSG_RESULT([in system libraries])
	oldlibs="$LIBS"
	LIBS="$LIBS $HEPMCLIBS"
	AC_CHECK_LIB(Rivet,main,
	[],
	[with_rivet=no
	AC_MSG_WARN([Rivet >= 1.3 not found in system libraries])
	])
	RIVETLIBS="$LIBS"
	LIBS=$oldlibs
	else
	AC_MSG_RESULT([$with_rivet])
	RIVETINCLUDE="$( $with_rivet/bin/rivet-config --cppflags )"
	RIVETLIBS="-L$with_rivet/lib -R$with_rivet/lib -lRivet"
	if test "${host_cpu}" == "x86_64" -a -e $with_rivet/lib64/libRivet.so ; then
	RIVETLIBS="-L$with_rivet/lib64 -R$with_rivet/lib64 -lRivet"
	fi
	fi

	if test "x$with_rivet" != "xno"; then
	LOAD_RIVET="library RivetAnalysis.so"
	# Now lets see if the libraries work properly
	oldLIBS="$LIBS"
	oldLDFLAGS="$LDFLAGS"
	oldCPPFLAGS="$CPPFLAGS"
	LIBS="$LIBS `echo $HEPMCLIBS \| sed -e 's! -R.* ! !'` `echo $RIVETLIBS \| sed -e 's! -R.* ! !'` `echo $GSLLIBS \| sed -e 's! -R.* ! !'`"
	LDFLAGS="$LDFLAGS"
	CPPFLAGS="$CPPFLAGS $HEPMCINCLUDE $RIVETINCLUDE $GSLINCLUDE"

	# check Rivet
	AC_MSG_CHECKING([that Rivet works])
	AC_LINK_IFELSE([AC_LANG_PROGRAM([[#include <Rivet/AnalysisHandler.hh>
	]],[[Rivet::AnalysisHandler foo; foo.writeData("foo");]])],[AC_MSG_RESULT([yes])],[AC_MSG_RESULT([no])
	AC_MSG_RESULT([No. Use '--with-rivet=' to set a path to Rivet >= 1.3'.])
	with_rivet="no"
	LOAD_RIVET=""
	])

	LIBS="$oldLIBS"
	LDFLAGS="$oldLDFLAGS"
	CPPFLAGS="$oldCPPFLAGS"
	fi

	rivetversion=1
	if test "x$with_rivet" = "xsystem"; then
	- echo $( rivet-config --version ) \| grep -q '^1\.' \|\| rivetversion=2
	+ echo $( rivet-config --version ) \| grep -q '^2\.' && rivetversion=2
	+ echo $( rivet-config --version ) \| grep -q '^3\.' && rivetversion=3
	elif test "x$with_rivet" != "xno"; then
	- echo $( "$with_rivet/bin/rivet-config" --version ) \| grep -q '^1\.' \|\| rivetversion=2
	+ echo $( "$with_rivet/bin/rivet-config" --version ) \| grep -q '^2\.' && rivetversion=2
	+ echo $( "$with_rivet/bin/rivet-config" --version ) \| grep -q '^3\.' && rivetversion=3
	fi

	-AC_DEFINE_UNQUOTED([ThePEG_RIVET_VERSION], [$rivetversion], [Rivet histogram variant (1=AIDA or 2=YODA)])
	+AC_DEFINE_UNQUOTED([ThePEG_RIVET_VERSION], [$rivetversion], [Rivet major version (1,2,3)])

	AM_CONDITIONAL(HAVE_RIVET,[test "x$with_rivet" != "xno"])
	AC_SUBST(RIVETINCLUDE)
	AC_SUBST(RIVETLIBS)
	AC_SUBST(LOAD_RIVET)
	])

File Metadata

Mime Type: text/x-diff
Expires: Sat, May 3, 6:25 AM (1 d, 12 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4976992
Default Alt Text: (19 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions