Herwig7Interface/interface/HepMC3Traits.h

0001 // -*- C++ -*-
0002 //
0003 // HepMCTraits.h is a part of ThePEG - Toolkit for HEP Event Generation
0004 // Copyright (C) 1999-2019 Leif Lonnblad
0005 //
0006 // ThePEG is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 #ifndef ThePEG_HepMC3Traits_H
0010 #define ThePEG_HepMC3Traits_H
0011
0012 #include <string>
0013 #include <map>
0014 #include <memory>
0015
0016 #include "HepMC3/GenEvent.h"
0017 #include "HepMC3/Units.h"
0018 #include <HepMC3/GenParticle.h>
0019 #include <HepMC3/GenPdfInfo.h>
0020
0021 #include "ThePEG/Config/ThePEG.h"
0022 #include <ThePEG/Repository/EventGenerator.h>
0023 #include <ThePEG/EventRecord/Event.h>
0024
0025 namespace HepMC3 {
0026   class GenVertex;
0027 }  // namespace HepMC3
0028
0029 namespace ThePEG {
0030
0031   /**
0032    * HepMCTraitsBase is a convenient base class for specializing the
0033    * HepMCTraits class to deal with different flavours of HepMC in the
0034    * HepMCConverter class. The default version will work for the CLHEP
0035    * implementation of HepMC. To use the HepMCConverter class for any
0036    * flavour of HepMC you have to specialize the HepMCTraits class
0037    * accordingly, possibly inheriting the functionality from the
0038    * HepMCTraitsBase class and only overriding the functions and
0039    * typedefs which are different. For the CLHEP flavour of HepMC you
0040    * only need to do <code>template&lt;&gt; struct
0041    * HepMCTraits&lt;HepMC::GenEvent&gt;: public
0042    * HepMCTraitsBase&lt;HepMC::GenEvent,HepMC::GenParticle,HepMC::GenVertex,
0043    * HepMC::Polarization&gt; {};</code> somewhere inside the ThePEG
0044    * namespace.  The boolean template argument determines whether the
0045    * HepMC implementation is specifying units or not.
0046    */
0047   template <typename HepMCEventT,
0048             typename HepMCParticleT,
0049             typename HepMCParticlePtrT,
0050             typename HepMCVertexT,
0051             typename HepMCVertexPtrT,
0052             typename HepMCPolarizationT,
0053             typename HepMCPdfInfoT>
0054
0055   struct HepMCTraitsBase {
0056     /** Typedef of the particle class. */
0057     typedef HepMCParticleT ParticleT;
0058
0059     /** Typedef of the event class. */
0060     typedef HepMCEventT EventT;
0061
0062     /** Typedef of the vertex class. */
0063     typedef HepMCVertexT VertexT;
0064
0065     /** Typedef of the polarization class. */
0066     typedef HepMCPolarizationT PolarizationT;
0067
0068     /** Typedef of the PdfInfo class. */
0069     typedef HepMCPdfInfoT PdfInfoT;
0070
0071     /** Typedef of a particle pointer */
0072     typedef HepMCParticlePtrT ParticlePtrT;
0073
0074     /** Typedef of a vertex pointer */
0075     typedef HepMCVertexPtrT VertexPtrT;
0076
0077     /** Create an event object with number \a evno and \a weight. */
0078     static EventT* newEvent(long evno, double weight, const map<string, double>& optionalWeights) {
0079       EventT* e = new EventT();
0080       e->set_event_number(evno);
0081       std::vector<std::string> wnames;
0082       std::vector<double> wvalues;
0083
0084       wnames.push_back("Default");
0085       wvalues.push_back(weight);
0086       for (map<string, double>::const_iterator w = optionalWeights.begin(); w != optionalWeights.end(); ++w) {
0087         wnames.push_back(w->first);
0088         wvalues.push_back(w->second);
0089       }
0090       e->run_info()->set_weight_names(wnames);
0091       e->weights() = wvalues;
0092
0093 #ifdef HEPMC_HAS_NAMED_WEIGHTS
0094       for (size_t i = 0; i < wnames.size(); i++)
0095         e->weights()[wnames[i]] = wvalues[i];
0096 #else
0097       e->weights() = wvalues;
0098 #endif
0099
0100       return e;
0101     }
0102
0103     /** Reset event weight and number of a re-used GenEvent. */
0104     static void resetEvent(EventT* e, long evno, double weight, const map<string, double>& optionalWeights) {
0105       e->set_event_number(evno);
0106       e->weights().clear();
0107       std::vector<std::string> wnames;
0108       std::vector<double> wvalues;
0109
0110       wnames.push_back("Default");
0111       wvalues.push_back(weight);
0112       for (map<string, double>::const_iterator w = optionalWeights.begin(); w != optionalWeights.end(); ++w) {
0113         wnames.push_back(w->first);
0114         wvalues.push_back(w->second);
0115       }
0116       e->run_info()->set_weight_names(wnames);
0117       e->weights() = wvalues;
0118
0119 #ifdef HEPMC_HAS_NAMED_WEIGHTS
0120       for (size_t i = 0; i < wnames.size(); i++)
0121         e->weights()[wnames[i]] = wvalues[i];
0122 #else
0123       e->weights() = wvalues;
0124 #endif
0125     }
0126
0127     /**
0128      * Return true if this version of HepMC accept user-defined units.
0129      */
0130     static bool hasUnits() {
0131 #ifdef HEPMC_HAS_UNITS
0132       return true;
0133 #else
0134       return false;
0135 #endif
0136     }
0137
0138     /**
0139      * Return the energy unit used in the installed version of HepMC.
0140      */
0141     static Energy defaultEnergyUnit() {
0142       //#ifndef HEPMC_HAS_UNITS
0143       return GeV;
0144       //#else
0145       //      return HepMC3::Units::default_momentum_unit() == HepMC3::Units::GEV ? GeV : MeV;
0146       //#endif
0147     }
0148
0149     /**
0150      * Return the length unit used in the installed version of HepMC.
0151      */
0152     static Length defaultLengthUnit() {
0153       //#ifndef HEPMC_HAS_UNITS
0154       return millimeter;
0155       //#else
0156       //      return HepMC3::Units::default_length_unit() == HepMC3::Units::MM ? millimeter : 10.0 * millimeter;
0157       //#endif
0158     }
0159
0160     /**
0161      * Return the momentum unit used by a given GenEvent object. If
0162      * HepMC does not support units this must return GeV.
0163      */
0164     static Energy momentumUnit(const EventT& e) {
0165 #ifdef HEPMC_HAS_UNITS
0166       return e.momentum_unit() == HepMC3::Units::MEV ? MeV : GeV;
0167 #else
0168       return GeV;
0169 #endif
0170     }
0171
0172     /**
0173      * Return the length unit used by a given GenEvent object. If
0174      * HepMC does not support units this must return millimeter.
0175      */
0176     static Length lengthUnit(const EventT& e) {
0177 #ifdef HEPMC_HAS_UNITS
0178       return e.length_unit() == HepMC3::Units::CM ? centimeter : millimeter;
0179 #else
0180       return millimeter;
0181 #endif
0182     }
0183
0184     /**
0185      * Set the units to be used by the given GenEvent object. If
0186      * HepMC does not support units this should be a no-op.
0187      */
0188 #ifdef HEPMC_HAS_UNITS
0189     static void setUnits(EventT& e, Energy momu, Length lenu) {
0190       e.use_units(momu == MeV ? HepMC3::Units::MEV : HepMC3::Units::GEV,
0191                   lenu == centimeter ? HepMC3::Units::CM : HepMC3::Units::MM);
0192     }
0193 #else
0194     static void setUnits(EventT&, Energy, Length) {}
0195 #endif
0196
0197     /** Set the \a scale, \f$\alpha_S\f$ (\a aS) and \f$\alpha_{EM}\f$
0198      (\a aEM) for the event \a e. The scale will be scaled with \a
0199      unit before given to the GenEvent. */
0200     static void setScaleAndAlphas(EventT& e, Energy2 scale, double aS, double aEM, Energy unit) {
0201       e.set_event_scale(sqrt(scale) / unit);
0202       e.set_alphaQCD(aS);
0203       e.set_alphaQED(aEM);
0204     }
0205
0206     /** Set the primary vertex, \a v, for the event \a e. */
0207     static void setSignalProcessVertex(EventT& e, VertexPtrT v) { e.set_signal_process_vertex(v); }
0208
0209     /** Set a vertex, \a v, for the event \a e. */
0210     static void addVertex(EventT& e, VertexPtrT v) { e.add_vertex(v); }
0211
0212     /** Create a new particle object with momentum \a p, PDG number \a
0213      id and status code \a status. The momentum will be scaled with
0214      \a unit which according to the HepMC documentation should be
0215      GeV. */
0216     static ParticlePtrT newParticle(const Lorentz5Momentum& p, long id, int status, Energy unit) {
0217       // Note that according to the documentation the momentum is stored in a
0218       // HepLorentzVector in GeV (event though the CLHEP standard is MeV).
0219       LorentzVector<double> p_scalar = p / unit;
0220       ParticlePtrT genp = new ParticleT(p_scalar, id, status);
0221       genp->setGeneratedMass(p.mass() / unit);
0222       return genp;
0223     }
0224
0225     /** Set the polarization directions, \a the and \a phi, for particle
0226      \a p. */
0227     static void setPolarization(ParticleT& genp, double the, double phi) {
0228       genp.set_polarization(PolarizationT(the, phi));
0229     }
0230
0231     /** Set the colour line (with index \a indx) to \a coline for
0232      particle \a p. */
0233     static void setColourLine(ParticleT& p, int indx, int coline) { p.set_flow(indx, coline); }
0234
0235     /** Create a new vertex. */
0236     static VertexPtrT newVertex() { return new VertexT(); }
0237
0238     /** Add an incoming particle, \a p, to the vertex, \a v. */
0239     static void addIncoming(VertexT& v, ParticlePtrT p) { v.add_particle_in(p); }
0240
0241     /** Add an outgoing particle, \a p, to the vertex, \a v. */
0242     static void addOutgoing(VertexT& v, ParticlePtrT p) { v.add_particle_out(p); }
0243
0244     /** Set the position \a p for the vertex, \a v. The length will be
0245      scaled with \a unit which normally should be millimeters. */
0246     static void setPosition(VertexT& v, const LorentzPoint& p, Length unit) {
0247       LorentzVector<double> p_scaled = p / unit;
0248       v.set_position(p_scaled);
0249     }
0250
0251     /** Set the beam particles for the event.*/
0252     static void setBeamParticles(EventT& e, ParticlePtrT p1, ParticlePtrT p2) {
0253       e.set_beam_particles(p1, p2);
0254       p1->set_status(4);
0255       p2->set_status(4);
0256     }
0257
0258     /** Set the PDF info for the event. */
0259 #ifdef HEPMC_HAS_PDF_INFO
0260     static void setPdfInfo(EventT& e, int id1, int id2, double x1, double x2, double scale, double xf1, double xf2) {
0261       HepMC3::GenPdfInfoPtr pdfinfo = std::make_shared<HepMC3::GenPdfInfo>();
0262       pdfinfo->set(id1, id2, x1, x2, scale, xf1, xf2);
0263       e.set_pdf_info(pdfinfo);
0264     }
0265 #else
0266     static void setPdfInfo(EventT&, int, int, double, double, double, double, double) {}
0267 #endif
0268
0269     /** Set the cross section info for the event. */
0270 #ifdef HEPMC_HAS_CROSS_SECTION
0271     static void setCrossSection(EventT& ev, double xs, double xserr) {
0272       std::shared_ptr<HepMC3::GenCrossSection> x = std::make_shared<HepMC3::GenCrossSection>();
0273       x->set_cross_section(xs, xserr);
0274       ev.set_cross_section(x);
0275     }
0276 #else
0277     static void setCrossSection(EventT&, double, double) {}
0278 #endif
0279   };
0280
0281   /**
0282    * The HepMCTraits class is used to deal with different flavours of
0283    * HepMC in the HepMCConverter class. To use the HepMCConverter class
0284    * for any flavour of HepMC you have to specialize the
0285    * HepMCTraits class accordingly, possibly inheriting the
0286    * functionality from the HepMCTraitsBase class and only overriding
0287    * the functions and typedefs which are different. For the CLHEP
0288    * flavour of HepMC you only need to do <code>template&lt;&gt; struct
0289    * HepMCTraits&lt;HepMC::GenEvent&gt;: public
0290    * HepMCTraitsBase&lt;HepMC::GenEvent,HepMC::GenParticle,HepMC::GenVertex,
0291    * HepMC::Polarization,HepMC::PdfInfo&gt; {};</code> somewhere inside the ThePEG
0292    * namespace.
0293    */
0294   template <typename HepMCEventT>
0295   struct HepMCTraits {};
0296
0297 }  // namespace ThePEG
0298
0299 #endif