Herwig6Instance.cc

Go to the documentation of this file.

#include <iostream>
#include <sstream>
#include <cstdlib>
#include <cstring>
#include <cctype>
#include <string>
#include <cassert>

#ifdef _POSIX_C_SOURCE
#include <sys/time.h>
#include <csignal>
#include <csetjmp>
#endif

#include <CLHEP/Random/RandomEngine.h>

#include <HepMC/HerwigWrapper.h>

#include "FWCore/ParameterSet/interface/FileInPath.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "GeneratorInterface/Herwig6Interface/interface/herwig.h"
#include "GeneratorInterface/Herwig6Interface/interface/Herwig6Instance.h"

#include "params.inc"

extern "C" void jminit_();
__attribute__((visibility("hidden"))) void dummy() {
  int dummyInt = 0;
  jminit_();
  jimmin_();
  jminit_();
  hwmsct_(&dummyInt);
  jmefin_();
}

// Added for reading in the paremeter files (orig. L. Sonnenschein)
extern "C" void lunread_(const char filename[], const int length);

using namespace gen;

// implementation for the Fortran callbacks from Herwig

// this gets a random number from the currently running instance
// so, if the have separate instances for let's say decayers based on
// Herwig than the Herwig pp collision instance, we get independent
// random numbers.  Also true for FastSim, since FastSim uses its own
// random numbers.  This means FastSim needs take a Herwig6Instance
// instance of its own instead of calling into Herwig directly.
double gen::hwrgen_(int *idummy) {
  Herwig6Instance *instance = FortranInstance::getInstance<Herwig6Instance>();
  assert(instance != nullptr);
  assert(instance->randomEngine != nullptr);
  return instance->randomEngine->flat();
}

void gen::cms_hwwarn_(char fn[6], int *code, int *exit) {
  std::string function(fn, 6);
  *exit = FortranInstance::getInstance<Herwig6Instance>()->hwwarn(function, *code);
}

extern "C" {
void hwaend_() {}

void cmsending_(int *ecode) {
  throw cms::Exception("Herwig6Error") << "Herwig6 stopped run with error code " << *ecode << "." << std::endl;
}
}

struct gen::TimeoutHolder {
  sigjmp_buf buf;
};

// Herwig6Instance methods

Herwig6Instance::Herwig6Instance() : randomEngine(nullptr), timeoutPrivate() {}

Herwig6Instance::~Herwig6Instance() {}

// timeout tool

#ifdef _POSIX_C_SOURCE
// some deep POSIX hackery to catch HERWIG sometimes (O(10k events) with
// complicated topologies) getting caught in and endless loop :-(

void Herwig6Instance::_timeout_sighandler(int signr) {
  Herwig6Instance *instance = FortranInstance::getInstance<Herwig6Instance>();
  assert(instance != nullptr);
  assert(instance->timeoutPrivate);
  siglongjmp(instance->timeoutPrivate->buf, 1);
}

bool Herwig6Instance::timeout(unsigned int secs, void (*fn)()) {
  if (timeoutPrivate)
    throw cms::Exception("ReentrancyProblem") << "Herwig6Instance::timeout() called recursively." << std::endl;
  struct sigaction saOld;
  std::memset(&saOld, 0, sizeof saOld);

  struct itimerval itv;
  timerclear(&itv.it_value);
  timerclear(&itv.it_interval);
  itv.it_value.tv_sec = 0;
  itv.it_interval.tv_sec = 0;
  setitimer(ITIMER_VIRTUAL, &itv, nullptr);

  sigset_t ss;
  sigemptyset(&ss);
  sigaddset(&ss, SIGVTALRM);

  sigprocmask(SIG_UNBLOCK, &ss, nullptr);
  sigprocmask(SIG_BLOCK, &ss, nullptr);

  timeoutPrivate = std::make_unique<TimeoutHolder>();
  if (sigsetjmp(timeoutPrivate->buf, 1)) {
    timeoutPrivate.reset();

    itv.it_value.tv_sec = 0;
    itv.it_interval.tv_sec = 0;
    setitimer(ITIMER_VIRTUAL, &itv, nullptr);
    sigprocmask(SIG_UNBLOCK, &ss, nullptr);
    return true;
  }

  itv.it_value.tv_sec = secs;
  itv.it_interval.tv_sec = secs;
  setitimer(ITIMER_VIRTUAL, &itv, nullptr);

  struct sigaction sa;
  std::memset(&sa, 0, sizeof sa);
  sa.sa_handler = &Herwig6Instance::_timeout_sighandler;
  sa.sa_flags = SA_ONESHOT;
  sigemptyset(&sa.sa_mask);

  sigaction(SIGVTALRM, &sa, &saOld);
  sigprocmask(SIG_UNBLOCK, &ss, nullptr);

  try {
    fn();
  } catch (...) {
    timeoutPrivate.reset();

    itv.it_value.tv_sec = 0;
    itv.it_interval.tv_sec = 0;
    setitimer(ITIMER_VIRTUAL, &itv, nullptr);

    sigaction(SIGVTALRM, &saOld, nullptr);

    throw;
  }

  timeoutPrivate.reset();

  itv.it_value.tv_sec = 0;
  itv.it_interval.tv_sec = 0;
  setitimer(ITIMER_VIRTUAL, &itv, nullptr);

  sigaction(SIGVTALRM, &saOld, nullptr);

  return false;
}
#else
bool Herwig6Instance::timeout(unsigned int secs, void (*fn)()) {
  fn();
  return false;
}
#endif

bool Herwig6Instance::hwwarn(const std::string &fn, int code) { return false; }

// regular Herwig6Instance methods

bool Herwig6Instance::give(const std::string &line) {
  typedef std::istringstream::traits_type traits;

  const char *p = line.c_str(), *q;
  p += std::strspn(p, " \t\r\n");

  for (q = p; std::isalnum(*q); q++)
    ;
  std::string name(p, q - p);

  const ConfigParam *param;
  for (param = configParams; param->name; param++)
    if (name == param->name)
      break;
  if (!param->name)
    return false;

  p = q + std::strspn(q, " \t\r\n");

  std::size_t pos = 0;
  std::size_t mult = 1;
  for (unsigned int i = 0; i < 3; i++) {
    if (!param->dim[i].size)
      break;

    if (*p++ != (i ? ',' : '('))
      return false;

    p += std::strspn(p, " \t\r\n");

    for (q = p; std::isdigit(*q); q++)
      ;
    std::istringstream ss(std::string(p, q - p));
    std::size_t index;
    ss >> index;
    if (ss.bad() || ss.peek() != traits::eof())
      return false;

    if (index < param->dim[i].offset)
      return false;
    index -= param->dim[i].offset;
    if (index >= param->dim[i].size)
      return false;

    p = q + std::strspn(q, " \t\r\n");

    pos += mult * index;
    mult *= param->dim[i].size;
  }

  if (param->dim[0].size) {
    if (*p++ != ')')
      return false;
    p += std::strspn(p, " \t\r\n");
  }

  if (*p++ != '=')
    return false;
  p += std::strspn(p, " \t\r\n");

  for (q = p; *q && (std::isalnum(*q) || std::strchr(".-+", *q)); q++)
    ;
  std::istringstream ss(std::string(p, q - p));

  p = q + std::strspn(q, " \t\r\n");
  if (*p && *p != '!')
    return false;

  switch (param->type) {
    case kInt: {
      int value;
      ss >> value;
      if (ss.bad() || ss.peek() != traits::eof())
        return false;

      ((int *)param->ptr)[pos] = value;
      break;
    }
    case kDouble: {
      double value;
      ss >> value;
      if (ss.bad() || ss.peek() != traits::eof())
        return false;

      ((double *)param->ptr)[pos] = value;
      break;
    }
    case kLogical: {
      std::string value_;
      ss >> value_;
      if (ss.bad() || ss.peek() != traits::eof())
        return false;

      for (std::string::iterator iter = value_.begin(); iter != value_.end(); ++iter)
        *iter = std::tolower(*iter);
      bool value;
      if (value_ == "yes" || value_ == "true" || value_ == "1")
        value = true;
      else if (value_ == "no" || value_ == "false" || value_ == "0")
        value = false;
      else
        return false;

      ((int *)param->ptr)[pos] = value;
      break;
    }
  }

  return true;
}

void Herwig6Instance::openParticleSpecFile(const std::string fileName) {
  edm::FileInPath fileAndPath(fileName);
  // WARING : This will call HWWARN if file does not exist.
  lunread_(fileAndPath.fullPath().c_str(), strlen(fileAndPath.fullPath().c_str()));

  return;
}