diff --git a/Analysis/NLORivetAnalysis.cc b/Analysis/NLORivetAnalysis.cc --- a/Analysis/NLORivetAnalysis.cc +++ b/Analysis/NLORivetAnalysis.cc @@ -1,310 +1,310 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the NLORivetAnalysis class. // #include #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 "Rivet/AnalysisHandler.hh" #include "Rivet/Tools/Logging.hh" using namespace ThePEG; NLORivetAnalysis::NLORivetAnalysis() - : _remnantId(82), _format(1),_unitchoice(), - _geneventPrecision(16), debug(false), _rivet(), _nevent(0) {} + : _remnantId(82),_unitchoice(), + debug(false), _rivet(), _nevent(0) {} namespace { // Special anonymous function for creating a genEvent. HepMC::GenEvent * makeEvent(tEventPtr event, tSubProPtr sub, long no, long remnantId, Energy eUnit, Length lUnit, CrossSection xsec, CrossSection xsecErr) { // generate beam particles const PPair& beam = event->incoming(); HepMC::GenParticlePtr b1 = HepMCTraits::newParticle(beam.first->momentum(),beam.first->id(), 1,eUnit); HepMC::GenParticlePtr b2 = HepMCTraits::newParticle(beam.second->momentum(),beam.second->id(), 1,eUnit); // generate remnants HepMC::GenParticlePtr r1 = HepMCTraits::newParticle(beam.first->momentum() - sub->incoming().first->momentum(), remnantId,1,eUnit); HepMC::GenParticlePtr 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::GenVertexPtr 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,_remnantId,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,_remnantId,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(make_pair(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 = 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/NLORivetAnalysis.h b/Analysis/NLORivetAnalysis.h --- a/Analysis/NLORivetAnalysis.h +++ b/Analysis/NLORivetAnalysis.h @@ -1,216 +1,206 @@ // -*- C++ -*- #ifndef THEPEG_NLORivetAnalysis_H #define THEPEG_NLORivetAnalysis_H // // This is the declaration of the NLORivetAnalysis class. // #include "ThePEG/Handlers/AnalysisHandler.h" #include "Rivet/AnalysisHandler.hh" namespace ThePEG { /** * Here is the documentation of the NLORivetAnalysis class. * * @see \ref NLORivetAnalysisInterfaces "The interfaces" * defined for NLORivetAnalysis. */ class NLORivetAnalysis: public ThePEG::AnalysisHandler { public: /** * The default constructor. */ NLORivetAnalysis(); public: /** @name Virtual functions required by the AnalysisHandler class. */ //@{ /** * Analyze a given Event. Note that a fully generated event * may be presented several times, if it has been manipulated in * between. The default version of this function will call transform * to make a lorentz transformation of the whole event, then extract * all final state particles and call analyze(tPVector) of this * analysis object and those of all associated analysis objects. The * default version will not, however, do anything on events which * have not been fully generated, or have been manipulated in any * way. * @param event pointer to the Event to be analyzed. * @param ieve the event number. * @param loop the number of times this event has been presented. * If negative the event is now fully generated. * @param state a number different from zero if the event has been * manipulated in some way since it was last presented. */ virtual void analyze(ThePEG::tEventPtr event, long ieve, int loop, int state); //@} 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(ThePEG::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(ThePEG::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 ThePEG::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 ThePEG::IBPtr fullclone() const; //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object. Called in the read phase. */ virtual void doinit(); /** * Initialize this object. Called in the run phase just before * a run begins. */ virtual void doinitrun(); /** * Finalize this object. Called in the run phase just after a * run has ended. Used eg. to write out statistics. */ virtual void dofinish(); //@} private: /** * The static object used to initialize the description of this class. * Indicates that this is a concrete class with persistent data. */ static ThePEG::ClassDescription initNLORivetAnalysis; /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ NLORivetAnalysis & operator=(const NLORivetAnalysis &) = delete; private: /** * The PDG ID to be used for remnants */ long _remnantId; /** - * The HepMC format - */ - int _format; - - /** * Selector for the choice of units */ int _unitchoice; /** - * Choice of output precision in GenEvent format - */ - unsigned int _geneventPrecision; - - /** * The Analyses to use */ vector _analyses; /** * The base name of the output file. */ string filename; /** * Enable debugging information from NLORivet */ bool debug; /** * The NLORivetAnalysisHandler */ Rivet::AnalysisHandler * _rivet; /** * Event count */ unsigned long _nevent; }; } #include "ThePEG/Utilities/ClassTraits.h" namespace ThePEG { /** @cond TRAITSPECIALIZATIONS */ /** This template specialization informs ThePEG about the * base classes of NLORivetAnalysis. */ template <> struct BaseClassTrait { /** Typedef of the first base class of NLORivetAnalysis. */ typedef AnalysisHandler NthBase; }; /** This template specialization informs ThePEG about the name of * the NLORivetAnalysis class and the shared object where it is defined. */ template <> struct ClassTraits : public ClassTraitsBase { /** Return a platform-independent class name */ static string className() { return "ThePEG::NLORivetAnalysis"; } /** * The name of a file containing the dynamic library where the class * NLORivetAnalysis is implemented. It may also include several, space-separated, * libraries if the class NLORivetAnalysis depends on other classes (base classes * excepted). In this case the listed libraries will be dynamically * linked in the order they are specified. */ static string library() { return "RivetAnalysis.so"; } }; /** @endcond */ } #endif /* THEPEG_NLORivetAnalysis_H */