diff --git a/Config/Containers.h b/Config/Containers.h --- a/Config/Containers.h +++ b/Config/Containers.h @@ -1,368 +1,369 @@ // -*- C++ -*- // // Containers.h is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2019 Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef ThePEG_Containers_H #define ThePEG_Containers_H /** \file * This file defines a number of containers. Some are just typedefs of * std containers, while others are wrappers around * std containers introduced in the hope of reducing the * amount of debugging and code duplication. * * Do not make changes in this file. If you need to modify any of the * standard containers used in ThePEG, edit a copy of this file and * include it in an alternative config file which can be included in * the main ThePEG.h config file * using the macro ThePEG_ALTERNATE_CONFIG. */ #include "ThePEG/Utilities/UnitIO.h" +#include "ThePEG/Config/Pointers.h" namespace ThePEG { /** A set of pointers to ParticleData objects. */ ThePEG_DECLARE_SET(PDPtr,ParticleDataSet); /** A vector of pointers to ParticleData objects. */ typedef vector PDVector; /** A vector of pointers to const ParticleData objects. */ typedef vector cPDVector; /** A vector of transient pointers to ParticleData objects. */ typedef vector tPDVector; /** A vector of transient pointers to const ParticleData objects. */ typedef vector tcPDVector; /** A set of pointers to MatcherBase objects. */ ThePEG_DECLARE_SET(PMPtr,MatcherSet); /** A set of pointers to DecayMode objects */ ThePEG_DECLARE_SET(DMPtr,DecayModeSet); /** A set of pointers to InterfacedBase objects */ ThePEG_DECLARE_SET(IBPtr,ObjectSet); /** A set of pointers to InterfacedBase objects */ ThePEG_DECLARE_SET(IBPtr,DependencySet); /** A map relating integers to ParticleData objects */ ThePEG_DECLARE_MAP(long,PDPtr,ParticleMap); /** A map relating character strings to InterfacedBase objects */ ThePEG_DECLARE_MAP(string,IBPtr,ObjectMap); /** A map relating InterfacedBase objects to corresponding * DependencySet containers */ ThePEG_DECLARE_MAP(IBPtr,DependencySet,DependencyMap); /** A vector of pointers to InterfacedBase objects. */ typedef vector IVector; /** A vector of pointers to const InterfacedBase objects. */ typedef vector CIVector; /** A vector of pointers to Particle objects. */ typedef vector ParticleVector; /** A vector of pointers to Particle objects. */ typedef vector PVector; /** A vector of pointers to const Particle objects. */ typedef vector cPVector; /** A vector of transient pointers to Particle objects. */ typedef vector tPVector; /** A vector of transient pointers to const Particle objects. */ typedef vector tcPVector; /** A list of pointers to Particle objects. */ typedef list ParticleList; /** A list of pointers to Particle objects. */ typedef list PList; /** A list of pointers to const Particle objects. */ typedef list cPList; /** A list of transient pointers to Particle objects. */ typedef list tPList; /** A list of transient pointers to const Particle objects. */ typedef list tcPList; /** A map relating character strings to bare pointers to InterfaceBase objects */ ThePEG_DECLARE_MAP(string,const InterfaceBase *,InterfaceMap); /** A rebinder for InterfacedBase objects. */ typedef Rebinder TranslationMap; /** A map relating character strings to EventGenerator objects */ ThePEG_DECLARE_MAP(string,EGPtr,GeneratorMap); /** A vector of pointers to AnalysisHandler objects. */ typedef vector AnalysisVector; /** A pair of pointers to ParticleData objects. */ typedef pair PDPair; /** A pair of pointers to const ParticleData objects. */ typedef pair cPDPair; /** A pair of transient pointers to ParticleData objects. */ typedef pair tPDPair; /** A pair of transient pointers to const ParticleData objects. */ typedef pair tcPDPair; /** A pair of pointers to Particle objects. */ typedef pair PPair; /** A pair of pointers to const Particle objects. */ typedef pair cPPair; /** A pair of transient pointers to const Particle objects. */ typedef pair tPPair; /** A pair of transient pointers to const Particle objects. */ typedef pair tcPPair; /** An Interval in scale. */ typedef Interval SInterval; /** A vector of intervals of scales. */ typedef vector SIntervalVector; /** A vector of pairs of transient pointers to PartonBins. */ typedef vector tPartonPairVec; /** A pair of transient pointers to ColourLine objects. */ typedef pair tColinePair; /** A set of pointers to DecayMode objects. */ ThePEG_DECLARE_SET(tDMPtr,DecaySet); /** A set oc character strings. */ ThePEG_DECLARE_SET(string,StringSet); /** A vector of energies. */ typedef vector EnergyVector; /** A vector of pointers to EventInfoBase objects. */ typedef vector EIVector; /** A vector of doubles. */ typedef vector DVector; /** A pair of doubles. */ typedef pair DPair; /** @name Global shift operators to simplify adding and removing * objects to containers. */ //@{ /** * Overload the left shift operator for vector to push_back objects to * a vector. * @param tv the vector being filled by push_back. * @param u the object being pushed back. * @return a referens to the vector. */ template vector & operator<<(vector & tv, const U & u) { tv.push_back(u); return tv; } /** * Overload the right shift operator for vector to pop objects from * a vector. * @param tv the vector being popped by pop_back. * @param u the object at the back of the vector before popping. * @return a referens to the vector. */ template vector & operator>>(vector & tv, U & u) { u = tv.back(); tv.pop_back(); return tv; } /** * Overload the left shift operator for stack to push objects to * a vector. * @param ts the stack being filled by push. * @param u the object being pushed. * @return a referens to the stack. */ template stack & operator<<(stack & ts, const U & u) { ts.push(u); return ts; } /** * Overload the right shift operator for stack to pop objects from * a stack. * @param ts the stack being popped. * @param u the object at the top of the stack before popping. * @return a referens to the stack. */ template stack & operator>>(stack & ts, U & u) { u = ts.top(); ts.pop(); return ts; } /** * Overload the left shift operator for deque to push_back objects to * a deque. * @param td the deque being filled by push_back. * @param u the object being pushed back. * @return a referens to the deque. */ template deque & operator<<(deque & td, const U & u) { td.push_back(u); return td; } /** * Overload the right shift operator for vector to pop objects from * a deque. * @param td the deque being popped by pop_front. * @param u the object at the front of the deque before popping. * @return a referens to the deque. */ template deque & operator>>(deque & td, U & u) { u = td.front(); td.pop_front(); return td; } /** * Overload the left shift operator for std::set to insert objects in * a set. * @param ts the set being filled by insert. * @param u the object being inserted. * @return a referens to the set. */ template set & operator<<(set & ts, const U & u) { ts.insert(u); return ts; } //@} /** @name Functions for I/O of containers of objects with unit. */ //@{ /** * Ouput a vector of objects with the specified unit. * @param os the stream used for output. * @param v the vector to be output. * @param u the unit to be used. */ template void ounitstream(OStream & os, const vector & v, UT & u) { os << v.size(); for ( typename vector::const_iterator i = v.begin(); i != v.end(); ++i ) os << ounit(*i, u); } /** * Input a vector of objects with the specified unit. * @param is the stream used for input. * @param v the vector to be input. * @param u the unit to be used. */ template void iunitstream(IStream & is, vector & v, UT & u) { typename vector::size_type l; is >> l; v.resize(l); for ( typename vector::iterator i = v.begin(); i != v.end(); ++i ) is >> iunit(*i, u); } /** * Ouput a set of objects with the specified unit. * @param os the stream used for output. * @param s the set to be output. * @param u the unit to be used. */ template void ounitstream(OStream & os, const set & s, UT & u) { os << s.size(); for ( typename set::const_iterator i = s.begin(); i != s.end(); ++i ) os << ounit(*i, u); } /** * Input a set of objects with the specified unit. * @param is the stream used for input. * @param s the set to be input. * @param u the unit to be used. */ template void iunitstream(IStream & is, set & s, UT & u) { s.clear(); typename set::size_type l; is >> l; T t; while ( l-- ) { is >> iunit(t, u); s.insert(t); } } /** * Ouput a map of keys and objects where the objects are output with * the specified unit. * @param os the stream used for output. * @param m the map to be output. * @param u the unit to be used for the mapped objects. */ template void ounitstream(OStream & os, const map & m, UT & u) { os << m.size(); for ( typename map::const_iterator i = m.begin(); i != m.end(); ++i ) os << i->first << ounit(i->second, u); } /** * Input a map of keys and objects where the objects are input with * the specified unit. * @param is the stream used for input. * @param m the map to be input. * @param u the unit to be used for the mapped objects. */ template void iunitstream(IStream & is, map & m, UT & u) { m.clear(); typename map::size_type l; is >> l; T t; K k; while ( l-- ) { is >> k >> iunit(t, u); m[k] = t; } } //@} } #endif /* ThePEG_Containers_H */ diff --git a/Config/HepMCHelper.h b/Config/HepMCHelper.h --- a/Config/HepMCHelper.h +++ b/Config/HepMCHelper.h @@ -1,180 +1,182 @@ // -*- C++ -*- // // HepMCHelper_HepMC.h is a part of ThePEG - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2019 The Herwig Collaboration // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is a helper header to implement HepMC conversions // #include "ThePEG/Vectors/HepMCTraits.h" #ifdef HAVE_HEPMC3 #include "HepMC3/GenEvent.h" #include "HepMC3/GenVertex.h" #include "HepMC3/GenParticle.h" #include "HepMC3/Version.h" #include "HepMC3/WriterAscii.h" #include "HepMC3/WriterHEPEVT.h" #include "HepMC3/WriterAsciiHepMC2.h" +#include "HepMC3/ReaderAscii.h" +#include "HepMC3/ReaderAsciiHepMC2.h" #ifdef HAVE_HEPMC3_WRITERROOT_H #include "HepMC3/WriterRoot.h" #endif #ifdef HAVE_HEPMC3_WRITERROOTTREE_H #include "HepMC3/WriterRootTree.h" #endif namespace HepMC3 { using PdfInfo=GenPdfInfo; using Polarization=std::pair; } namespace HepMC=HepMC3; #else #include "HepMC/GenEvent.h" #include "HepMC/GenVertex.h" #include "HepMC/GenParticle.h" #include "HepMC/Version.h" #include "HepMC/IO_BaseClass.h" #include "HepMC/IO_GenEvent.h" #include "HepMC/IO_AsciiParticles.h" namespace HepMC { #undef GenParticlePtr #undef GenVertexPtr typedef GenParticle * GenParticlePtr; typedef GenVertex * GenVertexPtr; } #endif namespace ThePEG { /** * Struct for HepMC conversion */ #ifndef HAVE_HEPMC3 template<> struct HepMCTraits : public HepMCTraitsBase { #else // This is version 3! template<> struct HepMCTraits : public HepMCTraitsBase { /** Create an event object with number \a evno and \a weight. */ static EventT * newEvent(long evno, double weight, const map& optionalWeights) { EventT * e = new EventT(HepMC::Units::GEV, HepMC::Units::MM); e->set_event_number(evno); e->set_event_number(evno); //std::vector wnames; std::vector wvalues; //wnames.push_back("Default"); wvalues.push_back(weight); for ( map::const_iterator w = optionalWeights.begin(); w != optionalWeights.end(); ++w ) { //wnames.push_back(w->first); wvalues.push_back(w->second); } //e->run_info()->set_weight_names(wnames); e->weights()=wvalues; return e; } /** Create a new vertex. */ static VertexPtrT newVertex() { return std::make_shared(VertexT()); } /** Set the \a scale, \f$\alpha_S\f$ (\a aS) and \f$\alpha_{EM}\f$ (\a aEM) for the event \a e. The scale will be scaled with \a unit before given to the GenEvent. */ static void setScaleAndAlphas(EventT & e, Energy2 scale, double aS, double aEM, Energy unit) { e.add_attribute("event_scale",std::make_shared(sqrt(scale)/unit)); e.add_attribute("mpi",std::make_shared(-1));//Please fix it later, once ThePEG authors respond e.add_attribute("signal_process_id",std::make_shared(0));//Please fix it later, once ThePEG authors respond e.add_attribute("alphaQCD",std::make_shared(aS)); e.add_attribute("alphaQED",std::make_shared(aEM)); } /** Set the colour line (with index \a indx) to \a coline for particle \a p. */ static void setColourLine(ParticleT & p, int indx, int coline) { p.add_attribute("flow"+std::to_string(indx),std::make_shared(coline)); } /** Add an incoming particle, \a p, to the vertex, \a v. */ static void addIncoming(VertexT & v, ParticlePtrT p) { v.add_particle_in(p); } /** Add an outgoing particle, \a p, to the vertex, \a v. */ static void addOutgoing(VertexT & v, ParticlePtrT p) { v.add_particle_out(p); } /** Set the primary vertex, \a v, for the event \a e. */ static void setSignalProcessVertex(EventT & e, VertexPtrT v) { e.add_vertex(v); e.add_attribute("signal_process_vertex", std::make_shared(v->id())); } /** Set a vertex, \a v, for the event \a e. */ static void addVertex(EventT & e, VertexPtrT v) { e.add_vertex(v); } /** Set the beam particles for the event.*/ static void setBeamParticles(EventT & e, ParticlePtrT p1, ParticlePtrT p2) { // e.set_beam_particles(p1,p2); p1->set_status(4); p2->set_status(4); e.set_beam_particles(p1, p2); } /** Create a new particle object with momentum \a p, PDG number \a id and status code \a status. The momentum will be scaled with \a unit which according to the HepMC documentation should be GeV. */ static ParticlePtrT newParticle(const Lorentz5Momentum & p, long id, int status, Energy unit) { // Note that according to the documentation the momentum is stored in a // HepLorentzVector in GeV (event though the CLHEP standard is MeV). HepMC::FourVector p_scalar(p.x()/unit, p.y()/unit, p.z()/unit, p.e()/unit); ParticlePtrT genp = std::make_shared(ParticleT(p_scalar, id, status)); genp->set_generated_mass(p.mass()/unit); return genp; } /** Set the polarization directions, \a the and \a phi, for particle \a p. */ static void setPolarization(ParticleT & genp, double the, double phi) { genp.add_attribute("theta",std::make_shared(the)); genp.add_attribute("phi",std::make_shared(phi)); } /** Set the position \a p for the vertex, \a v. The length will be scaled with \a unit which normally should be millimeters. */ static void setPosition(VertexT & v, const LorentzPoint & p, Length unit) { HepMC::FourVector v_scaled(p.x()/unit, p.y()/unit, p.z()/unit, p.e()/unit); v.set_position(v_scaled); } #endif }; } diff --git a/HepMC/HepMCEventHandler.cc b/HepMC/HepMCEventHandler.cc new file mode 100644 --- /dev/null +++ b/HepMC/HepMCEventHandler.cc @@ -0,0 +1,176 @@ +// -*- C++ -*- +// +// + + +#include "HepMCEventHandler.h" +#include "HepMCReader.h" +/* #include "ThePEG/" */ +#include "ThePEG/Interface/ClassDocumentation.h" +#include "ThePEG/Interface/Parameter.h" +#include "ThePEG/Interface/Reference.h" +#include "ThePEG/Utilities/Throw.h" +#include "ThePEG/Utilities/Debug.h" + +#include "ThePEG/Handlers/XComb.h" +#include "ThePEG/Handlers/CascadeHandler.h" +/* #include "ThePEG/Cuts/Cuts.h" */ +#include "ThePEG/PDF/PartonExtractor.h" + +#include "ThePEG/Persistency/PersistentOStream.h" +#include "ThePEG/Persistency/PersistentIStream.h" + + + +using namespace ThePEG; + +HepMCEventHandler::~HepMCEventHandler() {} + + +IBPtr HepMCEventHandler::clone() const { + return new_ptr(*this); +} + +IBPtr HepMCEventHandler::fullclone() const { + return new_ptr(*this); +} + + +void HepMCEventHandler::doinit() { + EventHandler::doinit(); + +} + +void HepMCEventHandler::doinitrun() { + // init reader + EventHandler::doinitrun(); + /* reader.initialize(*this); */ +} + + +void HepMCEventHandler::dofinish() { + EventHandler::dofinish(); +} + +void HepMCEventHandler::initialize() { + + if ( !reader() ) + Throw () + << "No reader were defined for the HepMCEventHandler '" + << name() << "'" << Exception::warning; + reader()->initialize(*this); + + PDPair incoming; + //TODO: just a work around + // set incoming particles fixed to protons + incoming.first = getParticleData(2212); + incoming.second = getParticleData(2212); + + theIncoming = incoming; +} + +EventPtr HepMCEventHandler::generateEvent() { + while( reader()->readEvent() ) { + /* Throw () */ + /* << "Couldn't read event from file " << reader()->filename() << "." */ + /* << Exception::runerror; */ + // TODO: Don't know if this works, in LesHouchesEventHandler LastXCombInfo::lastParticles + try { + theLastXComb = reader()->getXComb(); + // is used to access the incoming particles + /* PPair inParts = PPair(reader()->incoming()[0], reader()->incoming()[1]); */ + PPair inParts = PPair(reader()->incoming().first, reader()->incoming().second); + currentEvent(new_ptr(Event(inParts, this, generator()->runName(), + generator()->currentEventNumber()))); + + //TODO: not working right now + performCollision(); + //TODO: addParticles is a portected function + /* currentEvent()->addParticles(reader()->outgoing().begin(), reader()->outgoing().end()); */ + + return currentEvent(); + } + catch (Stop) { + } + catch (Exception &) { + throw; + } + } +} + +tCollPtr HepMCEventHandler::performCollision() { + lastExtractor()->select(lastXCombPtr()); + if ( CKKWHandler() ) CKKWHandler()->setXComb(lastXCombPtr()); + /* PPair inParts = PPair(reader()->incoming()[0], reader()->incoming()[1]); */ + /* PPair inParts = PPair(reader()->incoming().first, reader()->incoming().second); */ + currentCollision(new_ptr(Collision(lastParticles(), currentEvent(), this))); + + if ( currentEvent() ) currentEvent()->addCollision(currentCollision()); + currentStep(new_ptr(Step(currentCollision(), this))); + currentCollision()->addStep(currentStep()); + //TODO: need XComb with partons + currentStep()->addSubProcess(reader()->getSubProcess()); + lastExtractor()->constructRemnants(lastXCombPtr()->partonBinInstances(), + subProcess(), currentStep()); + /* if ( !currentReader()->cuts().passCuts(*currentCollision()) ) throw Veto(); */ + + initGroups(); + if ( ThePEG_DEBUG_ITEM(1) ) { + if ( currentEvent() ) + generator()->logfile() << *currentEvent(); + /* else */ + /* generator()->logfile() << *currentCollision(); */ + } + + return continueCollision(); +} + +EventPtr HepMCEventHandler::continueEvent() { + // TODO: Don't know at all if this will work + try { + continueCollision(); + } + catch (Veto) { + } + catch (Stop) { + } + catch (Exception &) { + throw; + } + return currentEvent(); +} + +void HepMCEventHandler::persistentOutput(PersistentOStream & os) const { + os << theReader << theIncoming; +} + +void HepMCEventHandler::persistentInput(PersistentIStream & is, int) { + is >> theReader >> theIncoming; +} + +ClassDescription HepMCEventHandler::initHepMCEventHandler; +// Definition of the static class description member. + +// init function +// with these paramters it is possible to set option in +// in file +void HepMCEventHandler::Init() { + static ClassDocumentation documentation + ("This is the main class administrating the event read " + "by ThePEG::HepMCReader object."); + + /* stores the reader Object */ + static Reference interfaceHepMCReader + ("HepMCReader", + "Object for reading HepMC files", + &HepMCEventHandler::theReader, + true, false, true, true, false); + /* TODO: more options */ + + /* static Reference interfacePartonExtractor */ + /* ("PartonExtractor", */ + /* "A PartonExtractor object to be used by subclasses which do not " */ + /* "provide their own. Note that this may be overridden by subclasses.", */ + /* &EventHandler::thePartonExtractor, true, false, true, true, false); */ + +} diff --git a/HepMC/HepMCEventHandler.fh b/HepMC/HepMCEventHandler.fh new file mode 100644 --- /dev/null +++ b/HepMC/HepMCEventHandler.fh @@ -0,0 +1,16 @@ +// -*- C++ -*- +// +// This is the forward declaration of the HepMCEventHandler class. +// +#ifndef ThePEG_HepMCEventHandler_FH +#define ThePEG_HepMCEventHandler_FH + +#include "ThePEG/Config/Pointers.h" + +namespace ThePEG { + +class HepMCEventHandler; +ThePEG_DECLARE_POINTERS(HepMCEventHandler,HepMCEventHandlerPtr); +} + +#endif diff --git a/HepMC/HepMCEventHandler.h b/HepMC/HepMCEventHandler.h new file mode 100644 --- /dev/null +++ b/HepMC/HepMCEventHandler.h @@ -0,0 +1,191 @@ + +#ifndef THEPEG_HepMCEventHandler_H +#define THEPEG_HepMCEventHandler_H + +#include "ThePEG/Handlers/EventHandler.h" +/* #include "HepMCReader.h" */ +#include "HepMCEventHandler.fh" +#include "HepMCReader.fh" + +namespace ThePEG { + + class HepMCEventHandler: public EventHandler { + public: + + /** + * + * The default constructor + */ + HepMCEventHandler() { + } + + /** + * The destructor + */ + virtual ~HepMCEventHandler(); + public: + + // With these functins variables of the class are + // stored to the disk + // + /** @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(); + + public: + /** @name Initialization and finalization functions. */ + //@{ + /** + * Initialize this event handler and all related objects needed to + * generate events. + */ + virtual void initialize(); + + + /** @name Functions used for the actual generation */ + //@{ + /** + * Generate an event. + */ + virtual EventPtr generateEvent(); + + /** + * Create the Event and Collision objects. Used by the + * generateEvent() function. + */ + virtual tCollPtr performCollision(); + + /** + * Continue generating an event if the generation has been stopped + * before finishing. + */ + virtual EventPtr continueEvent(); + //@} + /** + * Access the list of readers. + */ + const HepMCReaderPtr & reader() const { return theReader; } + + 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; + //@} + protected: + + /** + * Initialize this object. Called in the run phase just before + * a run begins. + */ + virtual void doinit(); + + virtual void doinitrun(); + + virtual void dofinish(); + + + private: + /** + * The reader + */ + HepMCReaderPtr theReader; + + public: + + /** @cond EXCEPTIONCLASSES */ + /** + * Exception class used if no readers were assigned. + */ + class HepMCInitError: public InitException {}; + + /** + * Exception class used if error while reading from file + */ + class HepMCOpenError: public Exception {}; + /** @endcond */ + + private: + + /** + * The static object used to initialize the description of this class. + * Indicates that this is a concrete class with persistent data. + */ + static ClassDescription initHepMCEventHandler; + + /** + * The assignment operator is private and must never be called. + * In fact, it should not even be implemented. + */ + HepMCEventHandler & operator=(const HepMCEventHandler &) = delete; + + + }; + + } + + // CLASSDOC OFF + // TODO: don't know what this is for + // It's a hook in the interface system of ThePEG + +#include "ThePEG/Utilities/ClassTraits.h" + + namespace ThePEG { + + /** @cond TRAITSPECIALIZATIONS */ + + /** This template specialization informs ThePEG about the + * base classes of HepMCEventHandler. */ + template <> + struct BaseClassTrait { + /** Typedef of the first base class of HepMCEventHandler. */ + typedef EventHandler NthBase; + }; + + /** This template specialization informs ThePEG about the name of + * the HepMCEventHandler 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::HepMCEventHandler"; } + /** Return the name of the shared library be loaded to get access to + * the HepMCEventHandler class and every other class it uses + * (except the base class). */ + static string library() { return "HepMCReader.so"; } + }; + + /** @endcond */ + + } +#endif /* THEPEG_HepMCEventHandler_H */ diff --git a/HepMC/HepMCReader.cc b/HepMC/HepMCReader.cc new file mode 100644 --- /dev/null +++ b/HepMC/HepMCReader.cc @@ -0,0 +1,622 @@ +// -*- C++ -*- +// +// HepMCReader.cc is a part of ThePEG - Toolkit for HEP Event Generation +// Copyright (C) 2020 Julian Lukwata +// +// ThePEG is licenced under version 3 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 HepMCReader class. +// + +/* #ifndef HAVE_HEPMC3 */ +/* #define HAVE_HEPMC3 1 */ +/* #endif */ + +#include "HepMCReader.h" +#include "HepMCEventHandler.h" +#ifdef HAVE_HEPMC3 +#include "HepMC3/ReaderAscii.h" +#include "HepMC3/ReaderAsciiHepMC2.h" +#include "HepMC3/GenEvent.h" +#include "HepMC3/GenParticle.h" +#include "HepMC3/GenVertex.h" +#include "HepMC3/FourVector.h" +#include "HepMC3/Attribute.h" +#else +#include "HepMC/IO_GenEvent.h" +#include "HepMC/GenEvent.h" +#include "HepMC/SimpeVector.h" +#endif + + +/* #include "ThePEG/Handlers/HandlerBase.h" */ + +#include "ThePEG/Interface/ClassDocumentation.h" +#include "ThePEG/Interface/Parameter.h" +#include "ThePEG/Interface/Switch.h" +#include "ThePEG/Interface/Reference.h" + +#include "ThePEG/Persistency/PersistentOStream.h" +#include "ThePEG/Persistency/PersistentIStream.h" + +/* #include "ThePEG/EventRecord/Event.h" */ +/* #include "ThePEG/PDF/NoPDF.h" */ +/* #include "ThePEG/PDF/BeamParticleData.h" */ +#include "ThePEG/PDF/RemnantHandler.h" +#include "ThePEG/PDF/PartonExtractor.h" +#include "ThePEG/PDF/PartonBin.h" +#include "ThePEG/Handlers/CascadeHandler.h" +#include "ThePEG/Cuts/Cuts.h" +#include "ThePEG/Handlers/XComb.h" +#include "ThePEG/Utilities/Throw.h" +/* #include "ThePEG/Utilities/HoldFlag.h" */ +/* #include "ThePEG/Utilities/Debug.h" */ +#include +#include + + +using namespace ThePEG; + +HepMCReader::HepMCReader() + : _eventNumber(1), _format(1), _filename(), + _unitchoice(), _event() {} + + // Cannot copy streams. + // Let doinitrun() take care of their initialization. +HepMCReader::HepMCReader(const HepMCReader & x) + : HandlerBase(x), LastXCombInfo<>(x), + _eventNumber(x._eventNumber), _format(x._format), _filename(x._filename), + _hepmcin(x._hepmcin), _unitchoice(x._unitchoice), _event(x._event), + thePartonExtractor(x.thePartonExtractor), theCKKW(x.theCKKW), + thePartonBins(x.thePartonBins), theXComb(x.theXComb), theCuts(x.theCuts), + theIncoming(x.theIncoming), theOutgoing(x.theOutgoing) {} + + +HepMCReader::~HepMCReader() {} + +IBPtr HepMCReader::fullclone() const { + return new_ptr(*this); +} + +IBPtr HepMCReader::clone() const { + return new_ptr(*this); +} + + +void HepMCReader::doinitrun() { + HandlerBase::doinitrun(); + // TODO: What does emax do + Energy emax = 100.0*GeV; + + /* PDPair inData = make_pair(getParticleData(2212), getParticleData(2212)); */ + + // partonBins() is not implemented + /* thePartonBins = partonExtractor()->getPartons(emax, inData, cuts()); */ + PDFCuts cuts1(cuts(), true, emax); + PBPair partonPair( + new_ptr(PartonBin(PDPtr(), PBPtr(), inData.first, PDFPtr(), cuts1)), + new_ptr(PartonBin(PDPtr(), PBPtr(), inData.second, PDFPtr(), cuts1))); + PBIPair partonBIPair( + new_ptr(PartonBinInstance(partonPair.first)), + new_ptr(PartonBinInstance(partonPair.second))); + xComb()->resetPartonBinInstances(partonBIPair); + + /* std::cout << "Size of PartonBins: " << partonBins().size() << std::endl; */ + /* theXComb = new_ptr(XComb(emax, inData, &eh, partonExtractor(), CKKWHandler(), */ + /* partonBins()[50], theCuts)); */ + + if ( _filename.empty() ) Throw() + << "No Filename for HepMC file was set." << Exception::runerror; +#ifdef HAVE_HEPMC3 + // create event to store stuff + /* _event = new HepMC3::GenEvent(); */ +#endif + + switch(_format){ + default: + case 1: { +#ifdef HAVE_HEPMC3 + _hepmcin = new HepMC3::ReaderAsciiHepMC2(_filename.c_str()); + // set reader option to correctly read flow + std::map< std::string, std::string > read_opt; + read_opt["particle_flows_are_separated"] = "1"; + _hepmcin->set_options(read_opt); +#else + _hepmcin = new HepMC::IO_GenEvent(_filename.c_str(), std::ios::in); +#endif + + } + break; +#ifdef HAVE_HEPMC3 + case 2: { + _hepmcin = new HepMC3::ReaderAscii(_filename.c_str()); + } +#endif + } + _eventNumber = 0; +} + +void HepMCReader::dofinish() { +#ifdef HAVE_HEPMC3 + _hepmcin->close(); + std::cout << "Closing AsciiReader.\n"; + /* delete _event; */ +#else + if (_hepmcin) { + delete _hepmcin; + _hepmcin = 0; + } +#endif + HandlerBase::dofinish(); + std::cout << "\nHepMCReader: read HepMC input.\n"; +} + +bool HepMCReader::doReadEvent() { +#ifdef HAVE_HEPMC3 + _event = std::make_shared(); + _hepmcin->read_event(*_event); +#else + _event = _hepmcin->read_next_event(); +#endif + // check if new evt was read + if ( !event() ) { + std::cout << "could not read event apparantly." << std::endl; + return false; + } + + return true; +} + +void HepMCReader::initialize(HepMCEventHandler & eh) { + + if ( !theCuts ) { + theCuts = eh.cuts(); + if ( !theCuts ) Throw() + << "No Cuts object was assigned to the HepMCReader '" + << name() << "' nor was one\nassigned to the controlling " + << "HepMCEventHandler '" << eh.name() << "'.\nAt least one of them " + << "needs to have a Cuts object." << Exception::runerror; + } + + theCKKW = eh.CKKWHandler(); + + if ( !partonExtractor() ) { + thePartonExtractor = eh.partonExtractor(); + if ( !partonExtractor() ) Throw() + << "No PartonExtractor object was assigned to the HepMCReader '" + << name() << "' nor was one\nassigned to the controlling " + << "HepMCEventHandler '" << eh.name() << "'.\nAt least one of them " + << "needs to have a PartonExtractor object." << Exception::runerror; + } + double E_double = 100.0; + Energy emax = E_double*GeV; + // TODO: just a work around + // set incoming particles fixed to protons + PDPair inData = make_pair(getParticleData(2212), getParticleData(2212)); + // partonBins() is not implemented + thePartonBins = partonExtractor()->getPartons(emax, inData, cuts()); + /* PDFCuts cuts1(cuts(), true, emax); */ + /* PBPair partonPair( */ + /* new_ptr(PartonBin(PDPtr(), PBPtr(), inData.first, PDFPtr(), cuts1)), */ + /* new_ptr(PartonBin(PDPtr(), PBPtr(), inData.second, PDFPtr(), cuts1))); */ + + /* for ( int i = 0, N = partonBins().size(); i < N; ++i ) { */ + /* theXCombs[partonBins()[i]] = */ + /* new_ptr(XComb(emax, inData, &eh, partonExtractor(), CKKWHandler(), */ + /* partonBins()[i], theCuts)); */ + /* partonExtractor()->nDims(partonBins()[i]); */ + /* } */ + /* outPDF = make_pair(partonExtractor()->getPDF(inData.first), */ + /* partonExtractor()->getPDF(inData.second)); */ + + /* std::cout << "Size of PartonBins: " << partonBins().size() << std::endl; */ + theXComb = new_ptr(XComb(emax, inData, &eh, partonExtractor(), CKKWHandler(), + partonBins()[1], + /* partonPair, */ + /* PBPair(new PartonBin(), new PartonBin()), */ + theCuts)); + /* theXComb = tXCombPtr(); */ + /* partonExtractor()->nDims(partonBins()[50]); */ + /* if(!xComb()->empty()) std::cout << "LastXComb is not empty." << std::endl; */ + /* std::cout << "the Max energy is : " << lastXCombPtr()->maxEnergy() << std::endl; */ + /* theXComb = new_ptr(XComb()); */ + /* _event = NULL; */ + +} + +bool HepMCReader::readEvent() { + + reset(); + + if ( !doReadEvent() ) return false; + + fillEvent(); + + return true; +} + +tXCombPtr HepMCReader::getXComb() { + if ( lastXCombPtr() ) return lastXCombPtr(); + fillEvent(); + connectMothers(); + //TODO: not implemented + /* tcPBPair sel = createPartonBinInstances(); */ + // set parton bin instances + thePartonBinInstances.first = + new_ptr(PartonBinInstance(incoming().first, partonBins()[0].first, + _e_scale)); + thePartonBinInstances.second = + new_ptr(PartonBinInstance(incoming().second, partonBins()[0].second, + _e_scale)); + /* PPair lastParts = PPair(incoming()[0], incoming()[1]); */ + tXCombPtr lastXC = xComb(); + // use default constructor + /* tXCombPtr lastXC = new_ptr(XComb()); */ + /* lastXC-lastParticles(lastParts); */ + // clean up the old XComb object before switching to a new one + if ( theLastXComb && theLastXComb != lastXC ) + theLastXComb->clean(); + theLastXComb = lastXC; + lastXCombPtr()->subProcess(SubProPtr()); + lastXCombPtr()->setPartonBinInstances(partonBinInstances(), + _e_scale); + /* if( !xComb()->empty()) std::cout << "xComb is not empyt." << std::endl; */ + +#ifdef HAVE_HEPMC3 + std::shared_ptr A_QCD = event()->attribute("alphaQCD"); + std::shared_ptr A_EM = event()->attribute("alphaEM"); + double a_QCD =A_QCD?(A_QCD->value()):0.0; + double a_EM = A_EM?(A_EM->value()):0.0; +#else + double a_QCD = event()->alphaQCD(); + double a_EM = event()->alphaQED(); +#endif + lastXCombPtr()->lastAlphaS(a_QCD); + lastXCombPtr()->lastAlphaEM(a_EM); + + return lastXCombPtr(); +} + +tSubProPtr HepMCReader::getSubProcess() { + getXComb(); + if ( subProcess() ) return subProcess(); + //TODO: Don't have lastPartons + lastXCombPtr()->subProcess(new_ptr(SubProcess(lastPartons(), tCollPtr(), this))); + subProcess()->setOutgoing(outgoing().begin(), outgoing().end()); + subProcess()->setIntermediates(intermediates().begin(), + intermediates().end()); + return subProcess(); +} + +void HepMCReader::fillEvent() { + if ( !particleIndex.empty() ) return; + // read eventNumber + _eventNumber = event()->event_number(); +#ifdef HAVE_HEPMC3 + std::shared_ptr E_scale = event()->attribute("event_scale"); + _e_scale = E_scale?(E_scale->value())*GeV2:0.0*GeV2; +#else + _e_scale = event()->event_scale()*GeV2; +#endif + //create particles + particleIndex.clear(); + colourIndex.clear(); + colourIndex(0, tColinePtr()); + // own functions + createParticles(); + //TODO: not implementes yet, don't know if needed + /* createBeams(); */ +} + +void HepMCReader::reset() { + particleIndex.clear(); + colourIndex.clear(); + if ( theLastXComb ) theLastXComb->clean(); + theLastXComb = tXCombPtr(); + /* delete _event; */ + if ( _eventNumber > 0 ) { + /* std::cout << "Clear event " << event()->event_number() */ + /* << "." << std::endl; */ +#ifdef HAVE_HEPMC3 + _event->clear(); +#else + delete _event; +#endif + _eventNumber = -1; + } +} + +void HepMCReader::createParticles() { + theBeams = PPair(); + /* PDPair incoming; */ + /* incoming = make_pair(getParticleData(2212), getParticleData(2212)); */ + theIncoming = PPair(); + /* theBeams = PVector(); */ + /* theIncoming = PVector(); */ + + /* incoming.first = getParticleData(2212); */ + /* incoming.second = getParticleData(2212); */ + + /* theIncoming = incoming; */ + theOutgoing = PVector(); + theIntermediates = PVector(); + + // to store the momentum + double x, y, z, t; + +#ifdef HAVE_HEPMC3 + for(auto part: event()->particles()) { + const HepMC3::FourVector partMom = part->momentum(); +#else + for( HepMC::GenEvent::particle_const_iterator pci = event()->particles_begin(); + pci != event()->particles_end(); ++pci ){ + HepMC::GenParticle * part = (*pci); + const HepMC::FourVector partMom = (*pci)->momentum(); +#endif + + x = partMom.px(); + y = partMom.py(); + z = partMom.pz(); + t = partMom.e(); + + Lorentz5Momentum mom(x*GeV, y*GeV, z*GeV, t*GeV); + + PDPtr pd = getParticleData(part->pdg_id()); + + PPtr p = pd->produceParticle(mom); + +#ifdef HAVE_HEPMC3 + tColinePtr c; + + std::shared_ptr Flow1(part->attribute("flow1")); + std::shared_ptr Flow2(part->attribute("flow2")); + if(Flow1) { + c = colourIndex(Flow1->value()); + c->addColoured(p); + } + if(Flow2){ + c = colourIndex(Flow2->value()); + c->addAntiColoured(p); + } + + // use id as identifier + particleIndex(part->id(), p); +#else + const HepMC::Flow * flow = &(*pci)->flow(); + + tColinePtr c = colourIndex(flow->icode(1)); + if ( c ) { + c->addColoured(p); + } + c = colourIndex(flow->icode(2)); + if ( c ) { + c->addAntiColoured(p); + } + // use barcode as identifier + particleIndex(part->barcode(), p); +#endif + + /* const ObjectIndexer & part_i = particleIndex; */ + // copy status + switch (part->status() ) { + // incoming beam particle + // TODO: use the first two particles as incoming particles + case 4: { +#ifdef HAVE_HEPMC3 + HepMC3::ConstGenVertexPtr ver = part->end_vertex(); + if(ver) { + x = ver->data().position.x(); + y = ver->data().position.y(); + z = ver->data().position.z(); + t = ver->data().position.t(); + + // hope it is in mm + p->setVertex(LorentzPoint(x*mm, y*mm, z*mm, t*mm)); + + } +#else + HepMC::GenVertex *ver = part->end_vertex(); + + if ( ver ) { + /* std::cout << "Particle status " << part->status() << std::endl; */ + /* std::cout << "x: " << std::scientific << ver->position().x() << std::endl; */ + /* std::cout << "y: " << ver->position().y() << std::endl; */ + /* std::cout << "z: " << ver->position().z() << std::endl; */ + x = ver->position().x(); + y = ver->position().y(); + z = ver->position().z(); + t = ver->position().t(); + // hope it is in mm + LorentzPoint pos(x*mm, y*mm, z*mm, t*mm); + // set vertex position + p->setVertex(pos); + } +#endif + if( !theIncoming.first ) { + theBeams.first = p; + theIncoming.first = p; + } + else if ( !theIncoming.second ) { + theBeams.second = p; + theIncoming.second = p; + } + else Throw() + << "To many incoming particles." << Exception::warning; + break; + } + //outgoing finalstate particle + case 1: + theOutgoing.push_back(p); + /* p->scale(sqr(hepeup.SCALUP*GeV)); */ + break; + // intermediate particle + case 11: + // decaying particles + case 2: + theIntermediates.push_back(p); + break; + default: + Throw() + << "Unknown status code (" << part->status() + << ") in the HepMCReader '" << name() << "'." + << Exception::runerror; + } + } +} + +/* void HepMCReader::createBeams() { */ +/* //TODO: Maybe not needed */ +/* } */ + +void HepMCReader::connectMothers() { + const ObjectIndexer & pi = particleIndex; + /* double test */ + double x, y, z, t; + +#ifdef HAVE_HEPMC3 + for(auto part: event()->particles()) { + /* std::cout << "Hallo." << std::endl; */ + HepMC3::ConstGenVertexPtr ver = part->production_vertex(); + // set position of production vertex + /* x = y = z = t = 0.0; */ + + if( ver ){ + x = ver->data().position.x(); + y = ver->data().position.y(); + z = ver->data().position.z(); + t = ver->data().position.t(); + + // hope it is in mm + LorentzPoint pos(x*mm, y*mm, z*mm, t*mm); + // set vertex position + if(part->status() != 4) + pi(part->id())->setVertex(pos); + + for(auto part_in: ver->particles_in()) { + pi(part_in->id())->addChild(pi(part->id())); + } + } +#else + for( HepMC::GenEvent::particle_const_iterator pci = event()->particles_begin(); + pci != event()->particles_end(); ++pci ){ + HepMC::GenVertex *ver = (*pci)->production_vertex(); + + if ( ver ) { + x = ver->position().x(); + y = ver->position().y(); + z = ver->position().z(); + t = ver->position().t(); + // hope it is in mm + LorentzPoint pos(x*mm, y*mm, z*mm, t*mm); + // set vertex position + if((*pci)->status() != 4) + pi((*pci)->barcode())->setVertex(pos); + + for(HepMC::GenVertex::particles_in_const_iterator pici = ver->particles_in_const_begin(); + pici != ver->particles_in_const_end(); ++pici ){ + pi((*pici)->barcode())->addChild(pi((*pci)->barcode())); + } + } +#endif + } +} + +void HepMCReader::persistentOutput(PersistentOStream & os) const { + os << _eventNumber << _format << _filename + << _unitchoice << theCuts + /* << _emax */ + << E_double + << theCKKW + << inData + << thePartonBins + << theXComb + << theIncoming << theOutgoing + /* << theLastXComb ; */ + << thePartonExtractor; +} + +void HepMCReader::persistentInput(PersistentIStream & is, int) { + is >> _eventNumber >> _format >> _filename + >> _unitchoice >> theCuts + /* >> _emax */ + >> E_double + >> theCKKW + >> inData + >> thePartonBins + >> theXComb + >> theIncoming >> theOutgoing + /* >> theLastXComb ; */ + >> thePartonExtractor; +} + +ClassDescription HepMCReader::initHepMCReader; +// Definition of the static class description member. + +void HepMCReader::Init() { + + static ClassDocumentation documentation + ("This event handler will read the event record from HepMC format."); + + static Parameter interfaceFilename + ("Filename", "Name of the input file", + &HepMCReader::_filename, ""); + + //TODO: Unit choice should be in input file + static Switch interfaceUnits + ("Units", + "Unit choice for energy and length", + &HepMCReader::_unitchoice, 0, false, false); + static SwitchOption interfaceUnitsGeV_mm + (interfaceUnits, + "GeV_mm", + "Use GeV and mm as units.", + 0); + static SwitchOption interfaceUnitsMeV_mm + (interfaceUnits, + "MeV_mm", + "Use MeV and mm as units.", + 1); + static SwitchOption interfaceUnitsGeV_cm + (interfaceUnits, + "GeV_cm", + "Use GeV and cm as units.", + 2); + static SwitchOption interfaceUnitsMeV_cm + (interfaceUnits, + "MeV_cm", + "Use MeV and cm as units.", + 3); + + static Reference interfaceCuts + ("Cuts", + "The Cuts object to be used for this reader. Note that these " + "must not be looser cuts than those used in the actual generation. " + "If no object is provided the HepMCEventHandler object must " + "provide one instead.", + &HepMCReader::theCuts, true, false, true, true, false); + + static Switch interfaceFormat + ("Format", +#ifdef HAVE_HEPMC3 + "Input format (1 = GenEvent, 2 = GenEventHepMC3", +#else + "Input format (1 = GenEvent", +#endif + &HepMCReader::_format, 1, false, false); + static SwitchOption interfaceFormatGenEvent + (interfaceFormat, + "GenEvent", + "IO_GenEvent format", + 1); + static SwitchOption interfaceFormatGenEventHepMC3 + (interfaceFormat, + "GenEventHepMC3", + "GenEvent in HepMC3", + 2); + + interfaceCuts.rank(8); + +} diff --git a/HepMC/HepMCReader.fh b/HepMC/HepMCReader.fh new file mode 100644 --- /dev/null +++ b/HepMC/HepMCReader.fh @@ -0,0 +1,16 @@ +// -*- C++ -*- +// +// This is the forward declaration of the HepMCReader class. +// +#ifndef ThePEG_HepMCReader_FH +#define ThePEG_HepMCReader_FH + +#include "ThePEG/Config/Pointers.h" + +namespace ThePEG { + +class HepMCReader; +ThePEG_DECLARE_CLASS_POINTERS(HepMCReader,HepMCReaderPtr); +} + +#endif diff --git a/HepMC/HepMCReader.h b/HepMC/HepMCReader.h new file mode 100644 --- /dev/null +++ b/HepMC/HepMCReader.h @@ -0,0 +1,795 @@ +// -*- C++ -*- +// +// HepMCReader.h is a part of ThePEG - Toolkit for HEP Event Generation +// Copyright (C) 2020 Julian Lukwata +// +// ThePEG is licenced under version 3 of the GPL, see COPYING for details. +// Please respect the MCnet academic guidelines, see GUIDELINES for details. +// +#ifndef THEPEG_HepMCReader_H +#define THEPEG_HepMCReader_H +// +// This is the declaration of the HepMCInFile class. +// +#include +#include +#include +#include "HepMCEventHandler.h" + + +#ifdef HAVE_HEPMC3 +#include "HepMC3/ReaderAscii.h" +#include "HepMC3/GenEvent.h" +#include "HepMC3/Units.h" +#else +#include "HepMC/IO_BaseClass.h" +#include "HepMC/GenEvent.h" +#endif + +#include "ThePEG/Handlers/HandlerBase.h" + +#include "ThePEG/Utilities/ObjectIndexer.h" +#include "ThePEG/Utilities/Exception.h" +/* #include "ThePEG/Config/Unitsystem.h" */ +/* #include "ThePEG/Utilities/XSecStat.h" */ +#include "ThePEG/PDF/PartonBinInstance.h" +#include "ThePEG/PDF/PartonBin.fh" +/* #include "ThePEG/MatrixElement/ReweightBase.h" */ +#include "HepMCReader.fh" +#include "HepMCEventHandler.fh" + + + +namespace ThePEG { + + /** + * HepMCReader is a class to be used for objects + * which reads event files . + * + * this base class will then be responsible for transforming this data to + * the ThePEG Event record in the getEvent() method. + * HepMCReaders can only be used inside HepMCEventHandler objects. + * + * + * The HepMCReader may be assigned a number ReweightBase objects + * which either completely reweights the events produced (in the + * reweights vector), or only biases the selection without influencing + * the cross section (in the preweights vector). Note that it is the + * responsibility of a sub-class to call the reweight() function and + * multiply the weight according to its return value (typically done + * in the readEvent() function). + * + * @see Event + * @see HepMCEventHandler + */ + class HepMCReader: public HandlerBase, public LastXCombInfo<> { + + /** + * HepMCEventHandler should have access to our private parts. + */ + friend class HepMCEventHandler; + + /** + * Map for accumulating statistics of cross sections per process + * number. + */ + /* typedef map StatMap; */ + + /** + * Map of XComb objects describing the incoming partons indexed by + * the corresponding PartonBin pair. + */ + /* typedef map XCombMap; */ + + + public: + + /** @name Standard constructors and destructors. */ + //@{ + /** + * Default constructor. + */ + HepMCReader(); + + /** + * Copy-constructor. + */ + HepMCReader(const HepMCReader &); + + /** + * Destructor. + */ + virtual ~HepMCReader(); + //@} + + public: + + /** @name Main virtual fuctions to be overridden in + * sub-classes. They are named as if we are reading from event + * files, but could equally well implement the actual generation of + * events. */ + //@{ + + /** + * Read the next event from the file or stream into the + * corresponding protected variables. Return false if there is no + * more events. + */ + virtual bool doReadEvent(); + + + /** + * return the weight names + */ + /* virtual vector optWeightsNamesFunc(); */ + + /** + * vector with the optional weights names + */ + + /* vector optionalWeightsNames; */ + + //@} + + /** @name Other important function which may be overridden in + * sub-classes which wants to bypass the basic HEPRUP or HEPEUP + * variables or otherwise facilitate the conversion to ThePEG + * objects. */ + //@{ + /** + * Initialize. This function is called by the HepMCEventHandler + * to which this object is assigned. + */ + virtual void initialize(HepMCEventHandler & eh); + + /** + * Calls readEvent() to read information into the + * variables. This function is called by the + * HepMCEventHandler if this reader has been selectod to + * produce an event. + * + * @return the weight asociated with this event. If negative weights + * are allowed it should be between -1 and 1, otherwise between 0 + * and 1. If outside these limits the previously estimated maximum + * is violated. Note that the estimated maximum then should be + * updated from the outside. + */ + /* virtual double getEvent(); */ + + /** + * Calls doReadEvent(). + */ + virtual bool readEvent(); + + /** + * Get an XComb object. Converts the information in the HepMC file + * to an XComb object describing the sub process. This is the way + * information is conveyed from the reader to the controlling HepMCEventHandler. + */ + tXCombPtr getXComb(); + + /** + * Get a SubProcess object corresponding to the information in the + */ + tSubProPtr getSubProcess(); + + /** + * Scan the file or stream to obtain information about cross section + * weights and particles etc. This function should fill the + * variables corresponding to the /HEPRUP/ common block. The + * function returns the number of events scanned. + */ + /* virtual long scan(); */ + + /** + * Take the information corresponding to the HEPRUP common block and + * initialize the statistics for this reader. + */ + /* virtual void initStat(); */ + + /** + * Reweights the current event using the reweights and preweights + * vectors. It is the responsibility of the sub-class to call this + * function after the HEPEUP information has been retrieved. + */ + /* double reweight(); */ + + /** + * Converts the information in the Les Houches common block + * variables into a Particle objects. + */ + virtual void fillEvent(); + + /** + * Removes the particles created in the last generated event, + * preparing to produce a new one. + */ + void reset(); + + //@} + + /** @name Access information about the current event. */ + //@{ + + /** + * The current event weight given by XWGTUP times possible + * reweighting. Note that this is not necessarily the same as what + * is returned by getEvent(), which is scaled with the maximum + * weight. + */ + /* double eventWeight() const { return hepeup.XWGTUP*lastweight; } */ + + /** + * Return the optional named weights associated to the current event. + */ + /* const map& optionalEventWeights() const { return optionalWeights; } */ + + /** + * Return the optional npLO and npNLO + */ + /* const int& optionalEventnpLO() const { return optionalnpLO; } */ + /* const int& optionalEventnpNLO() const { return optionalnpNLO; } */ + + /** + * The pair of PartonBinInstance objects describing the current + * incoming partons in the event. + */ + const PBIPair & partonBinInstances() const { return thePartonBinInstances; } + /** + * Return the instances of the beam particles for the current event. + */ + const PPair & beams() const { return theBeams; } + /* const PVector & beams() const { return theBeams; } */ + /** + * Return the instances of the incoming particles to the sub process + * for the current event. + */ + const PPair & incoming() const { return theIncoming; } + /* const PVector & incoming() const { return theIncoming; } */ + /** + * Return the instances of the outgoing particles from the sub process + * for the current event. + */ + const PVector & outgoing() const { return theOutgoing; } + /** + * Return the instances of the intermediate particles in the sub + * process for the current event. + */ + const PVector & intermediates() const { return theIntermediates; } + /** + * If this reader is to be used (possibly together with others) for + * CKKW reweighting and veto, this should give the multiplicity of + * outgoing particles in the highest multiplicity matrix element in + * the group. + */ + /* int maxMultCKKW() const { return theMaxMultCKKW; } */ + /** + * If this reader is to be used (possibly together with others) for + * CKKW reweighting and veto, this should give the multiplicity of + * outgoing particles in the lowest multiplicity matrix element in + * the group. + */ + /* int minMultCKKW() const { return theMinMultCKKW; } */ + //@} + + /** @name Other inlined access functions. */ + //@{ + /** + * The number of events found in this reader. If less than zero the + * number of events are unlimited. + */ + /* long NEvents() const { return theNEvents; } */ + + /** + * The number of events produced so far. Is reset to zero if an + * event file is reopened. + */ + /* long currentPosition() const { return position; } */ + + /** + * The maximum number of events to scan to collect information about + * processes and cross sections. If less than 0, all events will be + * scanned. + */ + /* long maxScan() const { return theMaxScan; } */ + + /** + * Return true if this reader is active. + */ + /* bool active() const { return isActive; } */ + + /** + * The collected cross section statistics for this reader. + */ + /* const XSecStat & xSecStats() const { return stats; } */ + + /** + * Collected statistics about the individual processes. + */ + /* const StatMap & processStats() const { return statmap; } */ + + /** + * Increase the overestimated cross section for this reader. + */ + /* virtual void increaseMaxXSec(CrossSection maxxsec); */ + + /** + * The PartonExtractor object used to construct remnants. + */ + tPExtrPtr partonExtractor() const { return thePartonExtractor; } + + /** + * Return a possibly null pointer to a CascadeHandler to be used for + * CKKW-reweighting. + */ + tCascHdlPtr CKKWHandler() const { return theCKKW; } + + /** + * The pairs of PartonBin objects describing the partons which can + * be extracted by the PartonExtractor object. + */ + const PartonPairVec & partonBins() const { return thePartonBins; } + + /** + * The map of XComb objects indexed by the corresponding PartonBin + * pair. + */ + /* const XCombMap & xComb() const { return theXComb; } */ + XCombPtr xComb() const { return theXComb; } + + /** + * The Cuts object to be used for this reader. + */ + /* CutsPtr cuts() const { return theCuts; } */ + const Cuts & cuts() const { return *theCuts; } + + /** + * + */ + const string filename() const { return _filename; } + + +#ifdef HAVE_HEPMC3 + std::shared_ptr event() { return std::const_pointer_cast(_event); } +#else + const HepMC::GenEvent * event() const { return _event; } +#endif + + //@} + + protected: + + + /** @name Auxilliary virtual methods which may be verridden by sub-classes. */ + //@{ + /** + * Check the existence of a pair of PartonBin objects corresponding + * to the current event. + * + * @return false if no pair of suitable PartonBin objects was found. + */ + /* virtual bool checkPartonBin(); */ + + /** + * Create instances of all particles in the event and store them + * in particleIndex. + */ + virtual void createParticles(); + + /** + * Using the already created particles create a pair of + * PartonBinInstance objects corresponding to the incoming + * partons. Return the corresponding PartonBin objects. + */ + /* virtual tcPBPair createPartonBinInstances(); */ + + /** + * Create instances of the incoming beams in the event and store + * them in particleIndex. If no beam particles are included in the + * event they are created from the run info. + */ + /* virtual void createBeams(); */ + + /** + * Go through the mother indices and connect up the Particles. + * TODO: Only needed if we consider intermediate steps :see_no_evil: + */ + virtual void connectMothers(); + //@} + + 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); + //@} + + /** + * Standard Init function used to initialize the interfaces. + */ + static void Init(); + + protected: + + /** @name Set functions for some variables not in the HepMC. */ + //@{ + /** + * The number of events in this reader. If less than zero the number + * of events is unlimited. + */ + /* void NEvents(long x) { theNEvents = x; } */ + + /** + * The map of XComb objects indexed by the corresponding PartuonBin + * pair. + */ + /* XCombMap & xCombs() { return theXCombs; } */ + //@} + + /** @name Standard (and non-standard) Interfaced functions. */ + //@{ + /** + * Initialize this object after the setup phase before saving an + * EventGenerator to disk. + * @throws InitException if object could not be initialized properly. + */ + /* 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(); + + /** + * Return true if this object needs to be initialized before all + * other objects because it needs to extract PDFs from the event file. + */ + /* virtual bool preInitialize() const; */ + + /** + * Called from doinit() to extract PDFs from the event file and add + * the corresponding objects to the current EventGenerator. + */ + /* virtual void initPDFs(); */ + //@} + + /** @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; + //@} + + protected: + + /** + * number of event + */ + int _eventNumber; + + /** + * The HepMC format + */ + int _format; + + /** + * The HepMC filename + */ + string _filename; + + /** + * The HepMC I/O handler + */ +#ifdef HAVE_HEPMC3 + HepMC3::Reader *_hepmcin; +#else + /* HepMC::IO_AsciiParticles *_hepmcio; */ + HepMC::IO_BaseClass *_hepmcin; +#endif + + /** + * Selector for the choice of units + */ + int _unitchoice; + + /** + * current HepMC event + */ +#ifdef HAVE_HEPMC3 + std::shared_ptr _event; // HepMC3::Units::GEV, HepMC3::Units:MM); +#else + HepMC::GenEvent *_event; +#endif + /** + * event scale + */ + Energy2 _e_scale; + + /** + * The ParticleData objects corresponding to the incoming particles. + */ + tcPDPair inData; + + /** + * The PDFBase objects which has been used for the beam particle + * when generating the events being read. Specified in the interface + * or derived from PDFGUP and PDFSUP. + */ + pair inPDF; + + /** + * The PDFBase object to be used in the subsequent generation. + */ + pair outPDF; + + /** + * The PartonExtractor object used to construct remnants. + */ + tPExtrPtr thePartonExtractor; + + /** + * A pointer to a CascadeHandler to be used for CKKW-reweighting. + */ + tCascHdlPtr theCKKW; + + /** + * The pairs of PartonBin objects describing the partons which can + * be extracted by the PartonExtractor object. + */ + PartonPairVec thePartonBins; + + /** + * The map of XComb objects indexed by the corresponding PartonBin + * pair. + */ + /* XCombMap theXCombs; */ + XCombPtr theXComb; + + /** + * The Cuts object to be used for this reader. + */ + CutsPtr theCuts; + + /** + * Determines whether to apply cuts to events before converting them + * to ThePEG format. + */ + bool doCutEarly; + + /** + * Collect statistics for this reader. + */ + /* XSecStat stats; */ + + /** + * Collect statistics for each individual process. + */ + /* StatMap statmap; */ + + /** + * The pair of PartonBinInstance objects describing the current + * incoming partons in the event. + */ + PBIPair thePartonBinInstances; + + /** + * Association between ColourLines and colour indices in the current + * translation. + */ + ObjectIndexer colourIndex; + + /** + * Association between Particles and indices in the current + * translation. + */ + ObjectIndexer particleIndex; + + /** + * The instances of the beam particles for the current event. + */ + PPair theBeams; + /* PVector theBeams; */ + + /** + * The instances of the incoming particles to the sub process for + * the current event. + */ + PPair theIncoming; + /* PVector theIncoming; */ + + /** + * The instances of the outgoing particles from the sub process for + * the current event. + */ + PVector theOutgoing; + + /** + * The instances of the intermediate particles in the sub process for + * the current event. + */ + PVector theIntermediates; + + /** + * Should PDFBase objects be constructed from the information in the + * event file in the initialization? + */ + bool doInitPDFs; + + /** + * If the maximum cross section of this reader has been increased + * with increaseMaxXSec(), this is the total factor with which it + * has been increased. + */ + /* double maxFactor; */ + + /** + * npLO for HepMC merging + */ + /* int optionalnpLO; */ + + /** + * npNLO for HepMC merging + */ + /* int optionalnpNLO; */ + + /** + * The (reweighted) XWGTUP value should be scaled with this cross + * section when compared to the overestimated cross section. + */ + /* CrossSection weightScale; */ + + /** + * Individual scales for different sub-processes if reweighted. + */ + /* vector xSecWeights; */ + + /** + * Individual maximum weights for individual (possibly reweighted) + * processes. + */ + /* map maxWeights; */ + + /** + * Option for the treatment of the momenta supplied + */ + /* unsigned int theMomentumTreatment; */ + + /** + * just for me + */ + + /* Energy _emax; */ + double E_double; + + private: + + /** Access function for the interface. */ + void setBeamA(long id); + /** Access function for the interface. */ + long getBeamA() const; + /** Access function for the interface. */ + void setBeamB(long id); + /** Access function for the interface. */ + long getBeamB() const; + /** Access function for the interface. */ + void setEBeamA(Energy e); + /** Access function for the interface. */ + Energy getEBeamA() const; + /** Access function for the interface. */ + void setEBeamB(Energy e); + /** Access function for the interface. */ + Energy getEBeamB() const; + /** Access function for the interface. */ + void setPDFA(PDFPtr); + /** Access function for the interface. */ + PDFPtr getPDFA() const; + /** Access function for the interface. */ + void setPDFB(PDFPtr); + /** Access function for the interface. */ + PDFPtr getPDFB() const; + + private: + + /** + * Describe an abstract base class with persistent data. + */ + static ClassDescription initHepMCReader; + + /** + * Private and non-existent assignment operator. + */ + HepMCReader & operator=(const HepMCReader &) = delete; + + public: + + /** @cond EXCEPTIONCLASSES */ + /** Exception class used by HepMCReader in case inconsistencies + * are encountered. */ + class HepMCInconsistencyError: public Exception {}; + + /** Exception class used by HepMCReader in case more events + than available are requested. */ + class HepMCReopenWarning: public Exception {}; + + /** Exception class used by HepMCReader in case there is + information missing in the initialization phase. */ + class HepMCInitError: public InitException {}; + /** @endcond */ + + }; + + /// Stream output for HEPEUP + /* ostream & operator<<(ostream & os, const HEPEUP & h); */ + +} + + +#include "ThePEG/Utilities/ClassTraits.h" + +namespace ThePEG { + + /** @cond TRAITSPECIALIZATIONS */ + + /** + * This template specialization informs ThePEG about the + * base class of HepMCReader. + */ + template <> + struct BaseClassTrait: public ClassTraitsType { + /** Typedef of the base class of HepMCReader. */ + typedef HandlerBase NthBase; + }; + + /** + * This template specialization informs ThePEG about the name of the + * HepMCReader class and the shared object where it is + * defined. + */ + template <> + struct ClassTraits + : public ClassTraitsBase { + /** + * Return the class name. + */ + static string className() { return "ThePEG::HepMCReader"; } + /** + * Return the name of the shared library to be loaded to get access + * to the HepMCReader class and every other class it uses + * (except the base class). + */ + static string library() { return "HepMCReader.so"; } + + }; + + /** @endcond */ + +} + + +#endif /* THEPEG_HepMCReader_H */ diff --git a/HepMC/Makefile.am b/HepMC/Makefile.am new file mode 100644 --- /dev/null +++ b/HepMC/Makefile.am @@ -0,0 +1,16 @@ +mySOURCES = HepMCReader.cc HepMCEventHandler.cc + +DOCFILES = HepMCEventHandler.h HepMCReader.h + +INCLUDEFILES = $(DOCFILES) HepMCEventHandler.fh \ + HepMCReader.fh + +pkglib_LTLIBRARIES = HepMCReader.la + +# Version info should be updated if any interface or persistent I/O +# function is changed +HepMCReader_la_LDFLAGS = $(AM_LDFLAGS) -module $(LIBTOOLVERSIONINFO) +HepMCReader_la_SOURCES = $(mySOURCES) $(INCLUDEFILES) + +include $(top_srcdir)/Config/Makefile.aminclude + diff --git a/Makefile.am b/Makefile.am --- a/Makefile.am +++ b/Makefile.am @@ -1,18 +1,19 @@ if JAVAGUI JAVADIR = java endif MAINDIRS = include Config Utilities Helicity Interface LesHouches Vectors \ PDT PDF Persistency Handlers MatrixElement Pointer ACDC \ - StandardModel Repository EventRecord Cuts Analysis lib src Doc + StandardModel Repository EventRecord Cuts Analysis lib src Doc \ + HepMC DIST_SUBDIRS = $(MAINDIRS) java SUBDIRS = $(MAINDIRS) $(JAVADIR) EXTRA_DIST = GUIDELINES ## DISTCHECK_CONFIGURE_FLAGS = --enable-unitchecks --without-lhapdf --with-gsl=$(GSLPATH) ACLOCAL_AMFLAGS = -I m4 DISTCLEANFILES = config.thepeg diff --git a/Vectors/HepMCConverter.tcc b/Vectors/HepMCConverter.tcc --- a/Vectors/HepMCConverter.tcc +++ b/Vectors/HepMCConverter.tcc @@ -1,337 +1,347 @@ // -*- C++ -*- // // HepMCConverter.tcc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2019 Leif Lonnblad // // ThePEG is licenced under version 3 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 HepMCConverter class. // #include "ThePEG/StandardModel/StandardModelBase.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/EventRecord/StandardSelectors.h" #include "ThePEG/EventRecord/Collision.h" #include "ThePEG/EventRecord/Step.h" #include "ThePEG/EventRecord/SubProcess.h" #include "ThePEG/Handlers/XComb.h" #include "ThePEG/Handlers/EventHandler.h" #include "ThePEG/PDF/PartonExtractor.h" #include "ThePEG/PDF/PDF.h" #include "ThePEG/PDT/StandardMatchers.h" #include "ThePEG/Utilities/Throw.h" namespace ThePEG { template typename HepMCConverter::GenEvent * HepMCConverter:: convert(const Event & ev, bool nocopies, Energy eunit, Length lunit) { HepMCConverter converter(ev, nocopies, eunit, lunit); return converter.geneve; } template void HepMCConverter:: convert(const Event & ev, GenEvent & gev, bool nocopies) { HepMCConverter converter(ev, gev, nocopies, Traits::momentumUnit(gev), Traits::lengthUnit(gev)); } template void HepMCConverter:: convert(const Event & ev, GenEvent & gev, bool nocopies, Energy eunit, Length lunit) { HepMCConverter converter(ev, gev, nocopies, eunit, lunit); } template HepMCConverter:: HepMCConverter(const Event & ev, bool nocopies, Energy eunit, Length lunit) : energyUnit(eunit), lengthUnit(lunit) { geneve = Traits::newEvent(ev.number(), ev.weight(), ev.optionalWeights()); init(ev, nocopies); } template HepMCConverter:: HepMCConverter(const Event & ev, GenEvent & gev, bool nocopies, Energy eunit, Length lunit) : energyUnit(eunit), lengthUnit(lunit) { geneve = &gev; Traits::resetEvent(geneve, ev.number(), ev.weight(), ev.optionalWeights()); init(ev, nocopies); } struct ParticleOrderNumberCmp { bool operator()(tcPPtr a, tcPPtr b) const { return a->number() < b->number(); } }; template void HepMCConverter::init(const Event & ev, bool nocopies) { if ( lengthUnit != millimeter && lengthUnit != centimeter ) Throw() << "Length unit used for HepMC::GenEvent was not MM nor CM." << Exception::runerror; if ( energyUnit != GeV && energyUnit != MeV ) Throw() << "Momentum unit used for HepMC::GenEvent was not GEV nor MEV." << Exception::runerror; Traits::setUnits(*geneve, energyUnit, lengthUnit); tcEHPtr eh; if ( ev.primaryCollision() && ( eh = dynamic_ptr_cast(ev.primaryCollision()->handler()) ) ) { // Get general event info if present. Traits::setScaleAndAlphas(*geneve, eh->lastScale(), eh->lastAlphaS(),eh->lastAlphaEM(), energyUnit); } // Extract all particles and order them. tcPVector all; ev.select(back_inserter(all), SelectAll()); stable_sort(all.begin(), all.end(), ParticleOrderNumberCmp()); vertices.reserve(all.size()*2); // Create GenParticle's and map them to the ThePEG particles. for ( int i = 0, N = all.size(); i < N; ++i ) { tcPPtr p = all[i]; if ( nocopies && p->next() ) continue; if ( pmap.find(p) != pmap.end() ) continue; pmap[p] = createParticle(p); if ( !p->children().empty() || p->next() ) { // If the particle has children it should have a decay vertex: vertices.push_back(Vertex()); decv[p] = &vertices.back(); vertices.back().in.insert(p); } if ( !p->parents().empty() || p->previous() || (p->children().empty() && !p->next()) ) { // If the particle has parents it should have a production // vertex. If neither parents or children it should still have a // dummy production vertex. vertices.push_back(Vertex()); prov[p] = &vertices.back(); vertices.back().out.insert(p); } } // Now go through the the particles again, and join the vertices. for ( int i = 0, N = all.size(); i < N; ++i ) { tcPPtr p = all[i]; if ( nocopies ) { if ( p->next() ) continue; for ( int i = 0, N = p->children().size(); i < N; ++i ) join(p, p->children()[i]->final()); tcPPtr pp = p; while ( pp->parents().empty() && pp->previous() ) pp = pp->previous(); for ( int i = 0, N = pp->parents().size(); i < N; ++i ) join(pp->parents()[i]->final(), p); } else { for ( int i = 0, N = p->children().size(); i < N; ++i ) join(p, p->children()[i]); if ( p->next() ) join(p, p->next()); for ( int i = 0, N = p->parents().size(); i < N; ++i ) join(p->parents()[i], p); if ( p->previous() ) join(p->previous(), p); } } // Time to create the GenVertex's for ( typename VertexMap::iterator it = prov.begin(); it != prov.end(); ++it ) if ( !member(vmap, it->second) ) vmap[it->second] = createVertex(it->second); for ( typename VertexMap::iterator it = decv.begin(); it != decv.end(); ++it ) if ( !member(vmap, it->second) ) vmap[it->second] = createVertex(it->second); // Now find the primary signal process vertex defined to be the // decay vertex of the first parton coming into the primary hard // sub-collision. + // TODO: for hepmc3 we need the right ordering +#ifdef HAVE_HEPMC3 + pmap[ev.incoming().first]->set_status(4); + pmap[ev.incoming().second]->set_status(4); + std::vector pvec; + for(auto const& it: pmap) + pvec.push_back(it.second); + geneve->add_tree(pvec); +#else tSubProPtr sub = ev.primarySubProcess(); if ( sub && sub->incoming().first ) { const Vertex * prim = decv[sub->incoming().first]; Traits::setSignalProcessVertex(*geneve, vmap[prim]); vmap.erase(prim); } // Then add the rest of the vertices. for ( typename GenVertexMap::iterator it = vmap.begin(); it != vmap.end(); ++it ) Traits::addVertex(*geneve, it->second); // and the incoming beam particles Traits::setBeamParticles(*geneve,pmap[ev.incoming().first], pmap[ev.incoming().second]); +#endif // and the PDF info setPdfInfo(ev); // and the cross section info Traits::setCrossSection(*geneve, eh->integratedXSec()/picobarn, eh->integratedXSecErr()/picobarn); for ( int i = 0, N = all.size(); i < N; ++i ) { tcPPtr p = all[i]; if ( pmap.find(p) == pmap.end() ) continue; GenParticlePtrT gp = pmap[p]; if ( p->hasColourInfo() ) { // Check if the particle is connected to colour lines, in which // case the lines are mapped to an integer and set in the // GenParticle's Flow info. tcColinePtr l; if ( (l = p->colourLine()) ) { if ( !member(flowmap, l) ) flowmap[l] = flowmap.size() + 500; Traits::setColourLine(*gp, 1, flowmap[l]); } if ( (l = p->antiColourLine()) ) { if ( !member(flowmap, l) ) flowmap[l] = flowmap.size() + 500; Traits::setColourLine(*gp, 2, flowmap[l]); } } if ( p->spinInfo() && p->spinInfo()->hasPolarization() ) { DPair pol = p->spinInfo()->polarization(); Traits::setPolarization(*gp, pol.first, pol.second); } } } template typename HepMCConverter::GenParticlePtrT HepMCConverter::createParticle(tcPPtr p) const { int status = 1; size_t nChildren = p->children().size(); if ( nChildren > 0 || p->next() ) status = 11; if ( nChildren > 1 ) { long id = p->data().id(); if ( BaryonMatcher::Check(id) || MesonMatcher::Check(id) || id == ParticleID::muminus || id == ParticleID::muplus || id == ParticleID::tauminus || id == ParticleID::tauplus ) { bool child = false; for(unsigned int ix=0;ixchildren()[ix]->id()==id) { child = true; break; } } if ( !child ) { if(p->data().widthCut()!=ZERO) { if(p->mass() <= p->data().massMax() && p->mass() >= p->data().massMin() ) status = 2; } else { status = 2; } } } } GenParticlePtrT gp = Traits::newParticle(p->momentum(), p->id(), p->status() ? p->status() : status, energyUnit); if ( p->spinInfo() && p->spinInfo()->hasPolarization() ) { DPair pol = p->spinInfo()->polarization(); Traits::setPolarization(*gp, pol.first, pol.second); } return gp; } template void HepMCConverter::join(tcPPtr parent, tcPPtr child) { Vertex * dec = decv[parent]; Vertex * pro = prov[child]; if ( !pro || !dec ) Throw() << "Found a reference to a ThePEG::Particle which was not in the Event." << Exception::eventerror; if ( pro == dec ) return; while ( !pro->in.empty() ) { dec->in.insert(*(pro->in.begin())); decv[*(pro->in.begin())] = dec; pro->in.erase(pro->in.begin()); } while ( !pro->out.empty() ) { dec->out.insert(*(pro->out.begin())); prov[*(pro->out.begin())] = dec; pro->out.erase(pro->out.begin()); } } template typename HepMCConverter::GenVertexPtrT HepMCConverter::createVertex(Vertex * v) { if ( !v ) Throw() << "Found internal null Vertex." << Exception::abortnow; GenVertexPtrT gv = Traits::newVertex(); // We assume that the vertex position is the average of the decay // vertices of all incoming and the creation vertices of all // outgoing particles in the lab. Note that this will probably not // be useful information for very small distances. LorentzPoint p; for ( tcParticleSet::iterator it = v->in.begin(); it != v->in.end(); ++it ) { p += (**it).labDecayVertex(); Traits::addIncoming(*gv, pmap[*it]); } for ( tcParticleSet::iterator it = v->out.begin(); it != v->out.end(); ++it ) { p += (**it).labVertex(); Traits::addOutgoing(*gv, pmap[*it]); } p /= double(v->in.size() + v->out.size()); Traits::setPosition(*gv, p, lengthUnit); return gv; } template void HepMCConverter::setPdfInfo(const Event & e) { // ids of the partons going into the primary sub process tSubProPtr sub = e.primarySubProcess(); int id1 = sub->incoming().first ->id(); int id2 = sub->incoming().second->id(); // get the event handler tcEHPtr eh = dynamic_ptr_cast(e.handler()); // get the values of x double x1 = eh->lastX1(); double x2 = eh->lastX2(); // get the pdfs pair pdfs; pdfs.first = eh->pdf(sub->incoming().first ); pdfs.second = eh->pdf(sub->incoming().second); // get the scale Energy2 scale = eh->lastScale(); // get the values of the pdfs double xf1 = pdfs.first.xfx(sub->incoming().first->dataPtr(), scale, x1); double xf2 = pdfs.second.xfx(sub->incoming().second->dataPtr(), scale, x2); Traits::setPdfInfo(*geneve, id1, id2, x1, x2, sqrt(scale/GeV2), xf1, xf2); } }