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 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::newParticle(beam.first->momentum(),beam.first->id(), 1,eUnit); HepMC::GenParticle * b2 = HepMCTraits::newParticle(beam.second->momentum(),beam.second->id(), 1,eUnit); // generate remnants HepMC::GenParticle * r1 = HepMCTraits::newParticle(beam.first->momentum() - sub->incoming().first->momentum(), _remnantId,1,eUnit); HepMC::GenParticle * r2 = HepMCTraits::newParticle(beam.second->momentum() - sub->incoming().second->momentum(), _remnantId,1,eUnit); // generate outgoing particles vector outgoing; for ( ParticleVector::const_iterator p = sub->outgoing().begin(); p != sub->outgoing().end(); ++p ) { outgoing.push_back(HepMCTraits::newParticle((**p).momentum(),(**p).id(), 1,eUnit)); } // generate one blob vertex HepMC::GenVertex * vertex = HepMCTraits::newVertex(); HepMCTraits::addIncoming(*vertex,b1); HepMCTraits::addIncoming(*vertex,b2); HepMCTraits::addOutgoing(*vertex,r1); HepMCTraits::addOutgoing(*vertex,r2); for ( vector::const_iterator p = outgoing.begin(); p != outgoing.end(); ++p ) HepMCTraits::addOutgoing(*vertex,*p); HepMC::GenEvent * ev = HepMCTraits::newEvent(no,event->weight()*sub->groupWeight(), event->optionalWeights()); HepMCTraits::setUnits(*ev,eUnit,lUnit); HepMCTraits::setBeamParticles(*ev,b1,b2); HepMCTraits::addVertex(*ev,vertex); HepMCTraits::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(event->primaryCollision()->handler()); assert(eh); CrossSection xsec = eh->integratedXSec(); CrossSection xsecErr = eh->integratedXSecErr(); tSubProPtr sub = event->primarySubProcess(); Ptr::tptr grp = dynamic_ptr_cast::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() << "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() << "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::initNLORivetAnalysis; // Definition of the static class description member. void NLORivetAnalysis::Init() { static ThePEG::ClassDocumentation documentation ("The NLORivetAnalysis class is a simple class to allow analyses" " from the Rivet library to be called from ThePEG"); static ThePEG::ParVector interfaceAnalyses ("Analyses", "The names of the Rivet analyses to use", &NLORivetAnalysis::_analyses, -1, "", "","" "", false, false, ThePEG::Interface::nolimits); static Parameter interfaceRemnantId ("RemnantId", "Set the PDG id to be used for remnants.", &NLORivetAnalysis::_remnantId, 82, 0, 0, false, false, Interface::nolimits); static Parameter 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 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 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::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() << "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::initRivetAnalysis; // Definition of the static class description member. void RivetAnalysis::Init() { static ThePEG::ClassDocumentation documentation ("The RivetAnalysis class is a simple class to allow analyses" " from the Rivet library to be called from ThePEG"); static ThePEG::ParVector interfaceAnalyses ("Analyses", "The names of the Rivet analyses to use", &RivetAnalysis::_analyses, -1, "", "","" "", false, false, ThePEG::Interface::nolimits); static ThePEG::ParVector interfacePaths ("Paths", "The directory paths where Rivet should look for analyses.", &RivetAnalysis::_paths, -1, "", "","" "", false, false, ThePEG::Interface::nolimits); static Parameter 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 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 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) ])