AlpgenInterface/plugins/AlpgenSource.cc

0001 #include <iomanip>
0002 #include <iostream>
0003 #include <memory>
0004
0005 #include <cmath>
0006 #include <fstream>
0007 #include <memory>
0008 #include <sstream>
0009 #include <string>
0010
0011 #include <boost/algorithm/string/classification.hpp>
0012 #include <boost/algorithm/string/split.hpp>
0013
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/Run.h"
0016 #include "FWCore/Framework/interface/MakerMacros.h"
0017 #include "FWCore/Framework/interface/InputSourceMacros.h"
0018 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0019 #include "FWCore/Utilities/interface/Exception.h"
0020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0021 #include "FWCore/Sources/interface/ProducerSourceFromFiles.h"
0022
0023 #include "SimDataFormats/GeneratorProducts/interface/LesHouches.h"
0024 #include "SimDataFormats/GeneratorProducts/interface/LHECommonBlocks.h"
0025 #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h"
0026 #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
0027
0028 #include "GeneratorInterface/AlpgenInterface/interface/AlpgenHeader.h"
0029 #include "GeneratorInterface/AlpgenInterface/interface/AlpgenEventRecordFixes.h"
0030
0031 class AlpgenSource : public edm::ProducerSourceFromFiles {
0032 public:
0033   /// Constructor
0034   AlpgenSource(const edm::ParameterSet &params, const edm::InputSourceDescription &desc);
0035
0036   /// Destructor
0037   ~AlpgenSource() override;
0038
0039 private:
0040   bool setRunAndEventInfo(edm::EventID &, edm::TimeValue_t &, edm::EventAuxiliary::ExperimentType &) override;
0041   void produce(edm::Event &event) override;
0042   void beginRun(edm::Run &run) override;
0043
0044   /// Function to get parameter by name from AlpgenHeader.
0045   template <typename T>
0046   T getParameter(AlpgenHeader::Parameter index) const;
0047   /// Function to get parameter by name from AlpgenHeader, w/ default.
0048   template <typename T>
0049   T getParameter(AlpgenHeader::Parameter index, const T &defValue) const;
0050
0051   /// Converts the AlpgenHeader::Masses to a std::string
0052   /// formatted as a slhaMassLine to facilitate passing them to Alpgen.
0053   std::string slhaMassLine(int pdgId, AlpgenHeader::Masses mass, const std::string &comment) const;
0054
0055   /// The Alpgen process ID. This is defined as
0056   /// processID() = 100*X + 10*Y + Z,
0057   /// where = ihrd, Y = ihvy, Z = njets
0058   unsigned int processID() const;
0059
0060   /// Read an event and put it into the HEPEUP.
0061   bool readAlpgenEvent(lhef::HEPEUP &hepeup);
0062
0063   /// Name of the input file
0064   std::string fileName_;
0065
0066   /// Pointer to the input file
0067   std::unique_ptr<std::ifstream> inputFile_;
0068
0069   /// Number of events to skip
0070   unsigned long skipEvents_;
0071
0072   /// Number of events
0073   unsigned long nEvent_;
0074
0075   /// Alpgen _unw.par file as a LHE header
0076   LHERunInfoProduct::Header lheAlpgenUnwParHeader;
0077   /// Alpgen _unw.par file as an AlpgenHeader
0078   AlpgenHeader header;
0079
0080   std::unique_ptr<lhef::HEPEUP> hepeup_;
0081
0082   /// Name of the extra header file
0083   std::string extraHeaderFileName_;
0084
0085   /// Name given to the extra header
0086   std::string extraHeaderName_;
0087
0088   /// configuration flags
0089   bool writeAlpgenWgtFile;
0090   bool writeAlpgenParFile;
0091   bool writeExtraHeader;
0092 };
0093
0094 AlpgenSource::AlpgenSource(const edm::ParameterSet &params, const edm::InputSourceDescription &desc)
0095     : edm::ProducerSourceFromFiles(params, desc, false),
0096       skipEvents_(params.getUntrackedParameter<unsigned int>("skipEvents", 0)),
0097       nEvent_(0),
0098       lheAlpgenUnwParHeader("AlpgenUnwParFile"),
0099       extraHeaderFileName_(params.getUntrackedParameter<std::string>("extraHeaderFileName", "")),
0100       extraHeaderName_(params.getUntrackedParameter<std::string>("extraHeaderName", "")),
0101       writeAlpgenWgtFile(params.getUntrackedParameter<bool>("writeAlpgenWgtFile", true)),
0102       writeAlpgenParFile(params.getUntrackedParameter<bool>("writeAlpgenParFile", true)),
0103       writeExtraHeader(params.getUntrackedParameter<bool>("writeExtraHeader", false)) {
0104   std::vector<std::string> allFileNames = fileNames(0);
0105
0106   // Only one filename
0107   if (allFileNames.size() != 1)
0108     throw cms::Exception("Generator|AlpgenInterface") << "AlpgenSource needs exactly one file specified "
0109                                                          "for now."
0110                                                       << std::endl;
0111
0112   fileName_ = allFileNames[0];
0113
0114   // Strip the "file:" prefix
0115   if (fileName_.find("file:") != 0)
0116     throw cms::Exception("Generator|AlpgenInterface") << "AlpgenSource only supports the file: scheme "
0117                                                          "for now."
0118                                                       << std::endl;
0119   fileName_.erase(0, 5);
0120
0121   // Open the _unw.par file to store additional
0122   // informations in the LHERunInfoProduct
0123   std::ifstream reader((fileName_ + "_unw.par").c_str());
0124   if (!reader.good())
0125     throw cms::Exception("Generator|AlpgenInterface")
0126         << "AlpgenSource was unable to open the file \"" << fileName_ << "_unw.par\"." << std::endl;
0127
0128   // A full copy of the _unw.par file in the LHE header.
0129   char buffer[256];
0130   while (reader.getline(buffer, sizeof buffer))
0131     lheAlpgenUnwParHeader.addLine(std::string(buffer) + "\n");
0132
0133   // Parse that header to setup an Alpgen header,
0134   // which will be used in the production itself.
0135   if (!header.parse(lheAlpgenUnwParHeader.begin(), lheAlpgenUnwParHeader.end()))
0136     throw cms::Exception("Generator|AlpgenInterface") << "AlpgenSource was unable to parse the Alpgen "
0137                                                       << "unweighted parameter file." << std::endl;
0138
0139   // Declare the products.
0140   produces<LHERunInfoProduct, edm::Transition::BeginRun>();
0141   produces<LHEEventProduct>();
0142 }
0143
0144 AlpgenSource::~AlpgenSource() {}
0145
0146 std::string AlpgenSource::slhaMassLine(int pdgId, AlpgenHeader::Masses mass, const std::string &comment) const {
0147   std::ostringstream ss;
0148   ss << std::setw(9) << pdgId << "     " << std::scientific << std::setprecision(9) << header.masses[mass] << "   # "
0149      << comment << std::endl;
0150   return ss.str();
0151 }
0152
0153 void AlpgenSource::beginRun(edm::Run &run) {
0154   // At this point, the lheUnwParHeader has the full contents of the _unw.par
0155   // file. So we can get the HEPRUP information from the LHE header itself.
0156   lhef::HEPRUP heprup;
0157
0158   // Get basic run information.
0159   // Beam identity.
0160   heprup.IDBMUP.first = 2212;
0161   switch (getParameter<int>(AlpgenHeader::ih2)) {
0162     case 1:
0163       heprup.IDBMUP.second = 2212;
0164       break;
0165     case -1:
0166       heprup.IDBMUP.second = -2212;
0167       break;
0168     default:
0169       throw cms::Exception("Generator|AlpgenInterface") << "AlpgenSource was unable to understand the ih2 "
0170                                                         << "parameter." << std::endl;
0171   }
0172
0173   // Beam energy.
0174   heprup.EBMUP.second = heprup.EBMUP.first = getParameter<double>(AlpgenHeader::ebeam);
0175
0176   // PDF info. Initially, Alpgen doesn't fill it.
0177   heprup.PDFGUP.first = -1;
0178   heprup.PDFGUP.second = -1;
0179   heprup.PDFSUP.first = -1;
0180   heprup.PDFSUP.second = -1;
0181
0182   // Unweighted events.
0183   heprup.IDWTUP = 3;
0184
0185   // Only one process.
0186   heprup.resize(1);
0187
0188   // Cross section and error.
0189   heprup.XSECUP[0] = header.xsec;
0190   heprup.XERRUP[0] = header.xsecErr;
0191
0192   // Maximum weight.
0193   heprup.XMAXUP[0] = header.xsec;
0194
0195   // Process code for Pythia.
0196   heprup.LPRUP[0] = processID();
0197
0198   // Comments on top.
0199   LHERunInfoProduct::Header comments;
0200   comments.addLine("\n");
0201   comments.addLine("\tExtracted by AlpgenInterface\n");
0202
0203   // Add SLHA header containing particle masses from Alpgen.
0204   // Pythia6Hadronisation will feed the masses to Pythia automatically.
0205   LHERunInfoProduct::Header slha("slha");
0206   slha.addLine("\n# SLHA header containing masses from Alpgen\n");
0207   slha.addLine("Block MASS      #  Mass spectrum (kinematic masses)\n");
0208   slha.addLine("#       PDG   Mass\n");
0209   slha.addLine(slhaMassLine(4, AlpgenHeader::mc, "charm    pole mass"));
0210   slha.addLine(slhaMassLine(5, AlpgenHeader::mb, "bottom   pole mass"));
0211   slha.addLine(slhaMassLine(6, AlpgenHeader::mt, "top      pole mass"));
0212   slha.addLine(slhaMassLine(23, AlpgenHeader::mz, "Z        mass"));
0213   slha.addLine(slhaMassLine(24, AlpgenHeader::mw, "W        mass"));
0214   slha.addLine(slhaMassLine(25, AlpgenHeader::mh, "H        mass"));
0215
0216   char buffer[512];
0217
0218   // We also add the information on weighted events.
0219   LHERunInfoProduct::Header lheAlpgenWgtHeader("AlpgenWgtFile");
0220   if (writeAlpgenWgtFile) {
0221     std::ifstream wgtascii((fileName_ + ".wgt").c_str());
0222     while (wgtascii.getline(buffer, 512)) {
0223       lheAlpgenWgtHeader.addLine(std::string(buffer) + "\n");
0224     }
0225   }
0226
0227   LHERunInfoProduct::Header lheAlpgenParHeader("AlpgenParFile");
0228   if (writeAlpgenParFile) {
0229     std::ifstream parascii((fileName_ + ".par").c_str());
0230     while (parascii.getline(buffer, 512)) {
0231       lheAlpgenParHeader.addLine(std::string(buffer) + "\n");
0232     }
0233   }
0234
0235   // If requested by the user, we also add any specific header provided.
0236   // Nota bene: the header is put in the LHERunInfoProduct AS IT WAS GIVEN.
0237   // That means NO CROSS-CHECKS WHATSOEVER. Use with care.
0238   LHERunInfoProduct::Header extraHeader(extraHeaderName_);
0239   if (writeExtraHeader) {
0240     std::ifstream extraascii(extraHeaderFileName_.c_str());
0241     while (extraascii.getline(buffer, 512)) {
0242       extraHeader.addLine(std::string(buffer) + "\n");
0243     }
0244   }
0245
0246   // Build the final Run info object. Backwards-compatible order.
0247   std::unique_ptr<LHERunInfoProduct> runInfo(new LHERunInfoProduct(heprup));
0248   runInfo->addHeader(comments);
0249   runInfo->addHeader(lheAlpgenUnwParHeader);
0250   if (writeAlpgenWgtFile)
0251     runInfo->addHeader(lheAlpgenWgtHeader);
0252   if (writeAlpgenParFile)
0253     runInfo->addHeader(lheAlpgenParHeader);
0254   runInfo->addHeader(slha);
0255   if (writeExtraHeader)
0256     runInfo->addHeader(extraHeader);
0257   run.put(std::move(runInfo));
0258
0259   // Open the .unw file in the heap, and set the global pointer to it.
0260   inputFile_ = std::make_unique<std::ifstream>((fileName_ + ".unw").c_str());
0261   if (!inputFile_->good())
0262     throw cms::Exception("Generator|AlpgenInterface")
0263         << "AlpgenSource was unable to open the file \"" << fileName_ << ".unw\"." << std::endl;
0264 }
0265
0266 template <typename T>
0267 T AlpgenSource::getParameter(AlpgenHeader::Parameter index) const {
0268   std::map<AlpgenHeader::Parameter, double>::const_iterator pos = header.params.find(index);
0269   if (pos == header.params.end())
0270     throw cms::Exception("Generator|AlpgenInterface")
0271         << "Requested Alpgen parameter \"" << AlpgenHeader::parameterName(index)
0272         << "\" "
0273            "not found in Alpgen parameter file."
0274         << std::endl;
0275
0276   return T(pos->second);
0277 }
0278
0279 template <typename T>
0280 T AlpgenSource::getParameter(AlpgenHeader::Parameter index, const T &defValue) const {
0281   std::map<AlpgenHeader::Parameter, double>::const_iterator pos = header.params.find(index);
0282   if (pos == header.params.end())
0283     return defValue;
0284   else
0285     return T(pos->second);
0286 }
0287
0288 unsigned int AlpgenSource::processID() const {
0289   // return 661; // The original, old thing.
0290   // digits #XYZ: X = ihrd, Y = ihvy, Z = njets
0291   return header.ihrd * 100 + getParameter<unsigned int>(AlpgenHeader::ihvy, 0) * 10 +
0292          getParameter<unsigned int>(AlpgenHeader::njets, 0);
0293 }
0294
0295 bool AlpgenSource::readAlpgenEvent(lhef::HEPEUP &hepeup) {
0296   char buffer[512];
0297   double dummy;
0298   int nPart;
0299   double sWgtRes;
0300   double sQ;
0301
0302   inputFile_->getline(buffer, sizeof buffer);
0303   if (!inputFile_->good())
0304     return false;
0305
0306   std::istringstream ls(buffer);
0307
0308   // Event number and process don't matter (or do they?)
0309   ls >> dummy >> dummy;
0310
0311   // Number of particles in the record, sample's average weight and Q scale
0312   ls >> nPart >> sWgtRes >> sQ;
0313
0314   if (ls.bad() || nPart < 1 || nPart > 1000)
0315     return false;
0316
0317   // Make room for the particles listed in the Alpgen file
0318   hepeup.resize(nPart);
0319
0320   // Scales, weights and process ID.
0321   hepeup.SCALUP = sQ;
0322   hepeup.XWGTUP = sWgtRes;
0323   hepeup.IDPRUP = processID();
0324
0325   // Incoming lines
0326   for (int i = 0; i != 2; ++i) {
0327     inputFile_->getline(buffer, sizeof buffer);
0328     std::istringstream ls(buffer);
0329     int flavour;
0330     ls >> flavour;
0331     int colour1;
0332     ls >> colour1;
0333     int colour2;
0334     ls >> colour2;
0335     double zMomentum;
0336     ls >> zMomentum;
0337
0338     if (inputFile_->bad())
0339       return false;
0340
0341     // Setting the HEPEUP of the incoming lines.
0342     hepeup.IDUP[i] = flavour;
0343     hepeup.ISTUP[i] = -1;
0344     hepeup.MOTHUP[i].first = 0;
0345     hepeup.MOTHUP[i].second = 0;
0346     hepeup.PUP[i][0] = 0.;
0347     hepeup.PUP[i][1] = 0.;
0348     hepeup.PUP[i][2] = zMomentum;
0349     hepeup.PUP[i][3] = std::fabs(zMomentum);
0350     hepeup.PUP[i][4] = 0.;
0351     if (colour1)
0352       colour1 += 500;
0353     if (colour2)
0354       colour2 += 500;
0355     hepeup.ICOLUP[i].first = colour1;
0356     hepeup.ICOLUP[i].second = colour2;
0357   }
0358
0359   // Outgoing lines
0360   for (int i = 2; i != nPart; ++i) {
0361     inputFile_->getline(buffer, sizeof buffer);
0362     std::istringstream ls(buffer);
0363     int flavour;
0364     ls >> flavour;
0365     int colour1;
0366     ls >> colour1;
0367     int colour2;
0368     ls >> colour2;
0369     double px, py, pz, mass;
0370     ls >> px >> py >> pz >> mass;
0371     double energy = std::sqrt(px * px + py * py + pz * pz + mass * mass);
0372
0373     if (inputFile_->bad())
0374       return false;
0375
0376     // Setting the HEPEUP of the outgoing lines.
0377     hepeup.IDUP[i] = flavour;
0378     hepeup.ISTUP[i] = 1;
0379     hepeup.MOTHUP[i].first = 1;
0380     hepeup.MOTHUP[i].second = 2;
0381     hepeup.PUP[i][0] = px;
0382     hepeup.PUP[i][1] = py;
0383     hepeup.PUP[i][2] = pz;
0384     hepeup.PUP[i][3] = energy;
0385     hepeup.PUP[i][4] = mass;
0386     if (colour1)
0387       colour1 += 500;
0388     if (colour2)
0389       colour2 += 500;
0390     hepeup.ICOLUP[i].first = colour1;
0391     hepeup.ICOLUP[i].second = colour2;
0392   }
0393
0394   return true;
0395 }
0396
0397 bool AlpgenSource::setRunAndEventInfo(edm::EventID &, edm::TimeValue_t &, edm::EventAuxiliary::ExperimentType &) {
0398   // The LHE Event Record
0399   hepeup_ = std::make_unique<lhef::HEPEUP>();
0400
0401   lhef::HEPEUP &hepeup = *hepeup_;
0402
0403   // Read the .unw file until it is over.
0404   for (;;) {
0405     if (!readAlpgenEvent(hepeup)) {
0406       if (inputFile_->eof())
0407         return false;
0408
0409       throw cms::Exception("Generator|AlpgenInterface")
0410           << "AlpgenSource is not able read event no. " << nEvent_ << std::endl;
0411     }
0412
0413     nEvent_++;
0414     if (skipEvents_ > 0)
0415       skipEvents_--;
0416     else
0417       break;
0418   }
0419   return true;
0420 }
0421
0422 void AlpgenSource::produce(edm::Event &event) {
0423   // Here are the Alpgen routines for filling up the rest
0424   // of the LHE Event Record. The .unw file has the information
0425   // in a compressed way, e.g. it doesn't list the W boson -
0426   // one has to reconstruct it from the e nu pair.
0427   lhef::HEPEUP &hepeup = *hepeup_;
0428
0429   switch (header.ihrd) {
0430     case 1:
0431     case 2:
0432     case 3:
0433     case 4:
0434     case 10:
0435     case 14:
0436     case 15:
0437       alpgen::fixEventWZ(hepeup);
0438       break;
0439     case 5:
0440       alpgen::fixEventMultiBoson(hepeup);
0441       break;
0442     case 6:
0443       alpgen::fixEventTTbar(hepeup);
0444       break;
0445     case 8:
0446       alpgen::fixEventHiggsTTbar(hepeup);
0447       break;
0448     case 13:
0449       alpgen::fixEventSingleTop(hepeup, header.masses[AlpgenHeader::mb], int(header.params[AlpgenHeader::itopprc]));
0450       break;
0451     case 7:
0452     case 9:
0453     case 11:
0454     case 12:
0455     case 16:
0456       // No fixes needed.
0457       break;
0458
0459     default:
0460       throw cms::Exception("Generator|AlpgenInterface") << "Unrecognized IHRD process code" << std::endl;
0461   }
0462
0463   // Create the LHEEventProduct and put it into the Event.
0464   std::unique_ptr<LHEEventProduct> lheEvent(new LHEEventProduct(hepeup));
0465   event.put(std::move(lheEvent));
0466
0467   hepeup_.reset();
0468 }
0469
0470 DEFINE_FWK_INPUT_SOURCE(AlpgenSource);