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