HepMCTraits.h

// -*- C++ -*-
//
// HepMCTraits.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_HepMCTraits_H
#define ThePEG_HepMCTraits_H

#ifdef HAVE_HEPMC3
#include "HepMC3/GenEvent.h"
namespace HepMC3 {
class GenEvent;
class GenParticle;
class GenVertex;
class GenPdfInfo;
}
namespace HepMC3 {
using PdfInfo=GenPdfInfo;
using Polarization=std::pair<double,double>;
}
namespace HepMC=HepMC3;
#else
#include "HepMC/GenEvent.h"
namespace HepMC {

class GenEvent;
class GenParticle;
class GenVertex;
class Polarization;
#ifndef HEPMC_GENPDFINFO_H
class PdfInfo;
#endif
}
#endif

 
namespace ThePEG {

template <typename HepMCEventT,
  typename HepMCParticleT, typename HepMCParticlePtrT,
  typename HepMCVertexT, typename HepMCVertexPtrT,
          typename HepMCPolarizationT, typename HepMCPdfInfoT>

struct HepMCTraitsBase {

  typedef HepMCParticleT ParticleT;

  typedef HepMCEventT EventT;

  typedef HepMCVertexT VertexT;

  typedef HepMCPolarizationT PolarizationT;

  typedef HepMCPdfInfoT PdfInfoT;

  typedef HepMCParticlePtrT ParticlePtrT;

  typedef HepMCVertexPtrT VertexPtrT;

  static EventT * newEvent(long evno, double weight,
                           const map<string,double>& optionalWeights) {
    EventT * e = new EventT();
    e->set_event_number(evno);
    std::vector<std::string> wnames;
    std::vector<double> wvalues;
    
    wnames.push_back("Default");
    wvalues.push_back(weight);
    for ( map<string,double>::const_iterator w = optionalWeights.begin();
          w != optionalWeights.end(); ++w ) {
    wnames.push_back(w->first);
    wvalues.push_back(w->second);
    }


#ifdef HAVE_HEPMC3
  e->run_info()->set_weight_names(wnames);
  e->weights()=wvalues;

#else 
#ifdef HEPMC_HAS_NAMED_WEIGHTS
      for (size_t i=0;i<wnames.size();i++) e->weights()[wnames[i]] = wvalues[i]; 
#else
      e->weights()=wvalues;
#endif
#endif

    
    
    return e;
  }

  static void resetEvent(EventT * e, long evno, double weight,
                         const map<string,double>& optionalWeights) {
    e->set_event_number(evno);
    e->weights().clear();
    std::vector<std::string> wnames;
    std::vector<double> wvalues;
    
    wnames.push_back("Default");
    wvalues.push_back(weight);
    for ( map<string,double>::const_iterator w = optionalWeights.begin();
          w != optionalWeights.end(); ++w ) {
    wnames.push_back(w->first);
    wvalues.push_back(w->second);
    }


#ifdef HAVE_HEPMC3
  e->run_info()->set_weight_names(wnames);
  e->weights()=wvalues;

#else 
#ifdef HEPMC_HAS_NAMED_WEIGHTS
      for (size_t i=0;i<wnames.size();i++) e->weights()[wnames[i]] = wvalues[i]; 
#else
      e->weights()=wvalues;
#endif
#endif
  }

  static bool hasUnits() {
#ifdef HEPMC_HAS_UNITS   
    return true;
#else
    return false;
#endif
  }

  static Energy defaultEnergyUnit() {
#ifndef HEPMC_HAS_UNITS
    return GeV;
#else
    return HepMC::Units::default_momentum_unit() == HepMC::Units::GEV? GeV: MeV;
#endif
  }

  static Length defaultLengthUnit() {
#ifndef HEPMC_HAS_UNITS
    return millimeter;
#else
    return HepMC::Units::default_length_unit() == HepMC::Units::MM?
      millimeter: 10.0*millimeter;
#endif
  }

  static Energy momentumUnit(const EventT & e) {
#ifdef HEPMC_HAS_UNITS   
    return e.momentum_unit() == HepMC::Units::MEV? MeV: GeV;
#else
    return GeV;
#endif
  }

  static Length lengthUnit(const EventT & e) {
#ifdef HEPMC_HAS_UNITS   
    return e.length_unit() == HepMC::Units::CM? centimeter: millimeter;
#else
    return millimeter;
#endif
  }

#ifdef HEPMC_HAS_UNITS   
  static void setUnits(EventT & e, Energy momu, Length lenu) {
    e.use_units(momu == MeV? HepMC::Units::MEV: HepMC::Units::GEV,
                lenu == centimeter? HepMC::Units::CM: HepMC::Units::MM);
  }
#else
  static void setUnits(EventT &, Energy, Length) {}
#endif


  static void setScaleAndAlphas(EventT & e, Energy2 scale,
                                double aS,  double aEM, Energy unit) {
    e.set_event_scale(sqrt(scale)/unit);
    e.set_alphaQCD(aS);
    e.set_alphaQED(aEM);
  }

  static void setSignalProcessVertex(EventT & e, VertexPtrT v) {
    e.set_signal_process_vertex(v);
  }

  static void addVertex(EventT & e, VertexPtrT v) {
    e.add_vertex(v);
  }

  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).
    LorentzVector<double> p_scalar = p/unit;
    ParticlePtrT genp =
      new ParticleT(p_scalar, id, status);
    genp->setGeneratedMass(p.mass()/unit);
    return genp;
  }

  static void setPolarization(ParticleT & genp, double the, double phi) {
    genp.set_polarization(PolarizationT(the, phi));
  }

  static void setColourLine(ParticleT & p, int indx, int coline) {
    p.set_flow(indx, coline);
  }

  static VertexPtrT newVertex() {
    return new VertexT();
  }

  static void addIncoming(VertexT & v, ParticlePtrT p) {
    v.add_particle_in(p);
  }

  static void addOutgoing(VertexT & v, ParticlePtrT p) {
    v.add_particle_out(p);
  }

  static void setPosition(VertexT & v, const LorentzPoint & p, Length unit) {
    LorentzVector<double> p_scaled = p/unit;
    v.set_position(p_scaled);
  }

  static void setBeamParticles(EventT & e, ParticlePtrT p1, ParticlePtrT p2) {
    e.set_beam_particles(p1,p2);
    p1->set_status(4);
    p2->set_status(4);
  }

#ifdef HEPMC_HAS_PDF_INFO
  static void setPdfInfo(EventT & e, int id1, int id2, double x1, double x2,
                         double scale, double xf1, double xf2) {
#ifdef HAVE_HEPMC3
    
    HepMC::GenPdfInfoPtr pdfinfo = std::make_shared<HepMC::GenPdfInfo>();
    pdfinfo->set(id1, id2, x1, x2, scale, xf1, xf2);
    e.set_pdf_info(pdfinfo);
#else
    e.set_pdf_info(PdfInfoT(id1, id2, x1, x2, scale, xf1, xf2));
#endif
  }
#else
  static void setPdfInfo(EventT &, int, int, double, double,
                         double, double, double) {}
#endif

#ifdef HEPMC_HAS_CROSS_SECTION
  static void setCrossSection(EventT & ev, double xs, double xserr) {
#ifdef HAVE_HEPMC3
     std::shared_ptr<HepMC::GenCrossSection> x =std::make_shared<HepMC::GenCrossSection>();
    x->set_cross_section(xs,xserr);
    ev.set_cross_section(x);
#else
    HepMC::GenCrossSection x;
    x.set_cross_section(xs, xserr);
    ev.set_cross_section(x);

#endif
  }
#else
  static void setCrossSection(EventT &, double, double) {}
#endif

};

template <typename HepMCEventT>
struct HepMCTraits {};
}

#endif