AlpgenSource.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>
#include <sstream>
#include <fstream>
#include <string>
#include <memory>
#include <cmath>

#include <boost/algorithm/string/classification.hpp>
#include <boost/algorithm/string/split.hpp>

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/Run.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/InputSourceMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Sources/interface/ProducerSourceFromFiles.h"

#include "SimDataFormats/GeneratorProducts/interface/LesHouches.h"
#include "SimDataFormats/GeneratorProducts/interface/LHECommonBlocks.h"
#include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"

#include "GeneratorInterface/AlpgenInterface/interface/AlpgenHeader.h"
#include "GeneratorInterface/AlpgenInterface/interface/AlpgenEventRecordFixes.h"

class AlpgenSource : public edm::ProducerSourceFromFiles {
public:
  AlpgenSource(const edm::ParameterSet &params, const edm::InputSourceDescription &desc);

  ~AlpgenSource() override;

private:
  bool setRunAndEventInfo(edm::EventID &, edm::TimeValue_t &, edm::EventAuxiliary::ExperimentType &) override;
  void produce(edm::Event &event) override;
  void beginRun(edm::Run &run) override;

  template <typename T>
  T getParameter(AlpgenHeader::Parameter index) const;
  template <typename T>
  T getParameter(AlpgenHeader::Parameter index, const T &defValue) const;

  std::string slhaMassLine(int pdgId, AlpgenHeader::Masses mass, const std::string &comment) const;

  unsigned int processID() const;

  bool readAlpgenEvent(lhef::HEPEUP &hepeup);

  std::string fileName_;

  std::unique_ptr<std::ifstream> inputFile_;

  unsigned long skipEvents_;

  unsigned long nEvent_;

  LHERunInfoProduct::Header lheAlpgenUnwParHeader;
  AlpgenHeader header;

  std::unique_ptr<lhef::HEPEUP> hepeup_;

  std::string extraHeaderFileName_;

  std::string extraHeaderName_;

  bool writeAlpgenWgtFile;
  bool writeAlpgenParFile;
  bool writeExtraHeader;
};

AlpgenSource::AlpgenSource(const edm::ParameterSet &params, const edm::InputSourceDescription &desc)
    : edm::ProducerSourceFromFiles(params, desc, false),
      skipEvents_(params.getUntrackedParameter<unsigned int>("skipEvents", 0)),
      nEvent_(0),
      lheAlpgenUnwParHeader("AlpgenUnwParFile"),
      extraHeaderFileName_(params.getUntrackedParameter<std::string>("extraHeaderFileName", "")),
      extraHeaderName_(params.getUntrackedParameter<std::string>("extraHeaderName", "")),
      writeAlpgenWgtFile(params.getUntrackedParameter<bool>("writeAlpgenWgtFile", true)),
      writeAlpgenParFile(params.getUntrackedParameter<bool>("writeAlpgenParFile", true)),
      writeExtraHeader(params.getUntrackedParameter<bool>("writeExtraHeader", false)) {
  std::vector<std::string> allFileNames = fileNames();

  // Only one filename
  if (allFileNames.size() != 1)
    throw cms::Exception("Generator|AlpgenInterface") << "AlpgenSource needs exactly one file specified "
                                                         "for now."
                                                      << std::endl;

  fileName_ = allFileNames[0];

  // Strip the "file:" prefix
  if (fileName_.find("file:") != 0)
    throw cms::Exception("Generator|AlpgenInterface") << "AlpgenSource only supports the file: scheme "
                                                         "for now."
                                                      << std::endl;
  fileName_.erase(0, 5);

  // Open the _unw.par file to store additional
  // informations in the LHERunInfoProduct
  std::ifstream reader((fileName_ + "_unw.par").c_str());
  if (!reader.good())
    throw cms::Exception("Generator|AlpgenInterface")
        << "AlpgenSource was unable to open the file \"" << fileName_ << "_unw.par\"." << std::endl;

  // A full copy of the _unw.par file in the LHE header.
  char buffer[256];
  while (reader.getline(buffer, sizeof buffer))
    lheAlpgenUnwParHeader.addLine(std::string(buffer) + "\n");

  // Parse that header to setup an Alpgen header,
  // which will be used in the production itself.
  if (!header.parse(lheAlpgenUnwParHeader.begin(), lheAlpgenUnwParHeader.end()))
    throw cms::Exception("Generator|AlpgenInterface") << "AlpgenSource was unable to parse the Alpgen "
                                                      << "unweighted parameter file." << std::endl;

  // Declare the products.
  produces<LHERunInfoProduct, edm::Transition::BeginRun>();
  produces<LHEEventProduct>();
}

AlpgenSource::~AlpgenSource() {}

std::string AlpgenSource::slhaMassLine(int pdgId, AlpgenHeader::Masses mass, const std::string &comment) const {
  std::ostringstream ss;
  ss << std::setw(9) << pdgId << "     " << std::scientific << std::setprecision(9) << header.masses[mass] << "   # "
     << comment << std::endl;
  return ss.str();
}

void AlpgenSource::beginRun(edm::Run &run) {
  // At this point, the lheUnwParHeader has the full contents of the _unw.par
  // file. So we can get the HEPRUP information from the LHE header itself.
  lhef::HEPRUP heprup;

  // Get basic run information.
  // Beam identity.
  heprup.IDBMUP.first = 2212;
  switch (getParameter<int>(AlpgenHeader::ih2)) {
    case 1:
      heprup.IDBMUP.second = 2212;
      break;
    case -1:
      heprup.IDBMUP.second = -2212;
      break;
    default:
      throw cms::Exception("Generator|AlpgenInterface") << "AlpgenSource was unable to understand the ih2 "
                                                        << "parameter." << std::endl;
  }

  // Beam energy.
  heprup.EBMUP.second = heprup.EBMUP.first = getParameter<double>(AlpgenHeader::ebeam);

  // PDF info. Initially, Alpgen doesn't fill it.
  heprup.PDFGUP.first = -1;
  heprup.PDFGUP.second = -1;
  heprup.PDFSUP.first = -1;
  heprup.PDFSUP.second = -1;

  // Unweighted events.
  heprup.IDWTUP = 3;

  // Only one process.
  heprup.resize(1);

  // Cross section and error.
  heprup.XSECUP[0] = header.xsec;
  heprup.XERRUP[0] = header.xsecErr;

  // Maximum weight.
  heprup.XMAXUP[0] = header.xsec;

  // Process code for Pythia.
  heprup.LPRUP[0] = processID();

  // Comments on top.
  LHERunInfoProduct::Header comments;
  comments.addLine("\n");
  comments.addLine("\tExtracted by AlpgenInterface\n");

  // Add SLHA header containing particle masses from Alpgen.
  // Pythia6Hadronisation will feed the masses to Pythia automatically.
  LHERunInfoProduct::Header slha("slha");
  slha.addLine("\n# SLHA header containing masses from Alpgen\n");
  slha.addLine("Block MASS      #  Mass spectrum (kinematic masses)\n");
  slha.addLine("#       PDG   Mass\n");
  slha.addLine(slhaMassLine(4, AlpgenHeader::mc, "charm    pole mass"));
  slha.addLine(slhaMassLine(5, AlpgenHeader::mb, "bottom   pole mass"));
  slha.addLine(slhaMassLine(6, AlpgenHeader::mt, "top      pole mass"));
  slha.addLine(slhaMassLine(23, AlpgenHeader::mz, "Z        mass"));
  slha.addLine(slhaMassLine(24, AlpgenHeader::mw, "W        mass"));
  slha.addLine(slhaMassLine(25, AlpgenHeader::mh, "H        mass"));

  char buffer[512];

  // We also add the information on weighted events.
  LHERunInfoProduct::Header lheAlpgenWgtHeader("AlpgenWgtFile");
  if (writeAlpgenWgtFile) {
    std::ifstream wgtascii((fileName_ + ".wgt").c_str());
    while (wgtascii.getline(buffer, 512)) {
      lheAlpgenWgtHeader.addLine(std::string(buffer) + "\n");
    }
  }

  LHERunInfoProduct::Header lheAlpgenParHeader("AlpgenParFile");
  if (writeAlpgenParFile) {
    std::ifstream parascii((fileName_ + ".par").c_str());
    while (parascii.getline(buffer, 512)) {
      lheAlpgenParHeader.addLine(std::string(buffer) + "\n");
    }
  }

  // If requested by the user, we also add any specific header provided.
  // Nota bene: the header is put in the LHERunInfoProduct AS IT WAS GIVEN.
  // That means NO CROSS-CHECKS WHATSOEVER. Use with care.
  LHERunInfoProduct::Header extraHeader(extraHeaderName_);
  if (writeExtraHeader) {
    std::ifstream extraascii(extraHeaderFileName_.c_str());
    while (extraascii.getline(buffer, 512)) {
      extraHeader.addLine(std::string(buffer) + "\n");
    }
  }

  // Build the final Run info object. Backwards-compatible order.
  std::unique_ptr<LHERunInfoProduct> runInfo(new LHERunInfoProduct(heprup));
  runInfo->addHeader(comments);
  runInfo->addHeader(lheAlpgenUnwParHeader);
  if (writeAlpgenWgtFile)
    runInfo->addHeader(lheAlpgenWgtHeader);
  if (writeAlpgenParFile)
    runInfo->addHeader(lheAlpgenParHeader);
  runInfo->addHeader(slha);
  if (writeExtraHeader)
    runInfo->addHeader(extraHeader);
  run.put(std::move(runInfo));

  // Open the .unw file in the heap, and set the global pointer to it.
  inputFile_.reset(new std::ifstream((fileName_ + ".unw").c_str()));
  if (!inputFile_->good())
    throw cms::Exception("Generator|AlpgenInterface")
        << "AlpgenSource was unable to open the file \"" << fileName_ << ".unw\"." << std::endl;
}

template <typename T>
T AlpgenSource::getParameter(AlpgenHeader::Parameter index) const {
  std::map<AlpgenHeader::Parameter, double>::const_iterator pos = header.params.find(index);
  if (pos == header.params.end())
    throw cms::Exception("Generator|AlpgenInterface")
        << "Requested Alpgen parameter \"" << AlpgenHeader::parameterName(index)
        << "\" "
           "not found in Alpgen parameter file."
        << std::endl;

  return T(pos->second);
}

template <typename T>
T AlpgenSource::getParameter(AlpgenHeader::Parameter index, const T &defValue) const {
  std::map<AlpgenHeader::Parameter, double>::const_iterator pos = header.params.find(index);
  if (pos == header.params.end())
    return defValue;
  else
    return T(pos->second);
}

unsigned int AlpgenSource::processID() const {
  // return 661; // The original, old thing.
  // digits #XYZ: X = ihrd, Y = ihvy, Z = njets
  return header.ihrd * 100 + getParameter<unsigned int>(AlpgenHeader::ihvy, 0) * 10 +
         getParameter<unsigned int>(AlpgenHeader::njets, 0);
}

bool AlpgenSource::readAlpgenEvent(lhef::HEPEUP &hepeup) {
  char buffer[512];
  double dummy;
  int nPart;
  double sWgtRes;
  double sQ;

  inputFile_->getline(buffer, sizeof buffer);
  if (!inputFile_->good())
    return false;

  std::istringstream ls(buffer);

  // Event number and process don't matter (or do they?)
  ls >> dummy >> dummy;

  // Number of particles in the record, sample's average weight and Q scale
  ls >> nPart >> sWgtRes >> sQ;

  if (ls.bad() || nPart < 1 || nPart > 1000)
    return false;

  // Make room for the particles listed in the Alpgen file
  hepeup.resize(nPart);

  // Scales, weights and process ID.
  hepeup.SCALUP = sQ;
  hepeup.XWGTUP = sWgtRes;
  hepeup.IDPRUP = processID();

  // Incoming lines
  for (int i = 0; i != 2; ++i) {
    inputFile_->getline(buffer, sizeof buffer);
    std::istringstream ls(buffer);
    int flavour;
    ls >> flavour;
    int colour1;
    ls >> colour1;
    int colour2;
    ls >> colour2;
    double zMomentum;
    ls >> zMomentum;

    if (inputFile_->bad())
      return false;

    // Setting the HEPEUP of the incoming lines.
    hepeup.IDUP[i] = flavour;
    hepeup.ISTUP[i] = -1;
    hepeup.MOTHUP[i].first = 0;
    hepeup.MOTHUP[i].second = 0;
    hepeup.PUP[i][0] = 0.;
    hepeup.PUP[i][1] = 0.;
    hepeup.PUP[i][2] = zMomentum;
    hepeup.PUP[i][3] = std::fabs(zMomentum);
    hepeup.PUP[i][4] = 0.;
    if (colour1)
      colour1 += 500;
    if (colour2)
      colour2 += 500;
    hepeup.ICOLUP[i].first = colour1;
    hepeup.ICOLUP[i].second = colour2;
  }

  // Outgoing lines
  for (int i = 2; i != nPart; ++i) {
    inputFile_->getline(buffer, sizeof buffer);
    std::istringstream ls(buffer);
    int flavour;
    ls >> flavour;
    int colour1;
    ls >> colour1;
    int colour2;
    ls >> colour2;
    double px, py, pz, mass;
    ls >> px >> py >> pz >> mass;
    double energy = std::sqrt(px * px + py * py + pz * pz + mass * mass);

    if (inputFile_->bad())
      return false;

    // Setting the HEPEUP of the outgoing lines.
    hepeup.IDUP[i] = flavour;
    hepeup.ISTUP[i] = 1;
    hepeup.MOTHUP[i].first = 1;
    hepeup.MOTHUP[i].second = 2;
    hepeup.PUP[i][0] = px;
    hepeup.PUP[i][1] = py;
    hepeup.PUP[i][2] = pz;
    hepeup.PUP[i][3] = energy;
    hepeup.PUP[i][4] = mass;
    if (colour1)
      colour1 += 500;
    if (colour2)
      colour2 += 500;
    hepeup.ICOLUP[i].first = colour1;
    hepeup.ICOLUP[i].second = colour2;
  }

  return true;
}

bool AlpgenSource::setRunAndEventInfo(edm::EventID &, edm::TimeValue_t &, edm::EventAuxiliary::ExperimentType &) {
  // The LHE Event Record
  hepeup_.reset(new lhef::HEPEUP);

  lhef::HEPEUP &hepeup = *hepeup_;

  // Read the .unw file until it is over.
  for (;;) {
    if (!readAlpgenEvent(hepeup)) {
      if (inputFile_->eof())
        return false;

      throw cms::Exception("Generator|AlpgenInterface")
          << "AlpgenSource is not able read event no. " << nEvent_ << std::endl;
    }

    nEvent_++;
    if (skipEvents_ > 0)
      skipEvents_--;
    else
      break;
  }
  return true;
}

void AlpgenSource::produce(edm::Event &event) {
  // Here are the Alpgen routines for filling up the rest
  // of the LHE Event Record. The .unw file has the information
  // in a compressed way, e.g. it doesn't list the W boson -
  // one has to reconstruct it from the e nu pair.
  lhef::HEPEUP &hepeup = *hepeup_;

  switch (header.ihrd) {
    case 1:
    case 2:
    case 3:
    case 4:
    case 10:
    case 14:
    case 15:
      alpgen::fixEventWZ(hepeup);
      break;
    case 5:
      alpgen::fixEventMultiBoson(hepeup);
      break;
    case 6:
      alpgen::fixEventTTbar(hepeup);
      break;
    case 8:
      alpgen::fixEventHiggsTTbar(hepeup);
      break;
    case 13:
      alpgen::fixEventSingleTop(hepeup, header.masses[AlpgenHeader::mb], int(header.params[AlpgenHeader::itopprc]));
      break;
    case 7:
    case 9:
    case 11:
    case 12:
    case 16:
      // No fixes needed.
      break;

    default:
      throw cms::Exception("Generator|AlpgenInterface") << "Unrecognized IHRD process code" << std::endl;
  }

  // Create the LHEEventProduct and put it into the Event.
  std::unique_ptr<LHEEventProduct> lheEvent(new LHEEventProduct(hepeup));
  event.put(std::move(lheEvent));

  hepeup_.reset();
}

DEFINE_FWK_INPUT_SOURCE(AlpgenSource);