Pythia8Interface/plugins/LHAupLesHouches.cc

0001 #include <algorithm>
0002 #include <iostream>
0003 #include <iterator>
0004 #include <sstream>
0005 #include <string>
0006 #include <memory>
0007 #include <cassert>
0008
0009 #include "GeneratorInterface/Pythia8Interface/plugins/LHAupLesHouches.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011
0012 using namespace Pythia8;
0013
0014 bool LHAupLesHouches::setInit() {
0015   if (!runInfo)
0016     return false;
0017   const lhef::HEPRUP &heprup = *runInfo->getHEPRUP();
0018
0019   setBeamA(heprup.IDBMUP.first, heprup.EBMUP.first, heprup.PDFGUP.first, heprup.PDFSUP.first);
0020   setBeamB(heprup.IDBMUP.second, heprup.EBMUP.second, heprup.PDFGUP.second, heprup.PDFSUP.second);
0021   setStrategy(heprup.IDWTUP);
0022
0023   for (int i = 0; i < heprup.NPRUP; i++)
0024     addProcess(heprup.LPRUP[i], heprup.XSECUP[i], heprup.XERRUP[i], heprup.XMAXUP[i]);
0025
0026   //hadronisation->onInit().emit();
0027
0028   //runInfo.reset();
0029
0030   //fill SLHA header information if available
0031   std::vector<std::string> slha = runInfo->findHeader("slha");
0032   if (!slha.empty()) {
0033     std::string slhaheader;
0034     for (std::vector<std::string>::const_iterator iter = slha.begin(); iter != slha.end(); ++iter) {
0035       slhaheader.append(*iter);
0036     }
0037     infoPtr->setHeader("slha", slhaheader);
0038   }
0039
0040   //will be used to work around missing initialization inside pythia8
0041   if (!fEvAttributes) {
0042     fEvAttributes = new std::map<std::string, std::string>;
0043   } else {
0044     fEvAttributes->clear();
0045   }
0046
0047   return true;
0048 }
0049
0050 bool LHAupLesHouches::setEvent(int inProcId) {
0051   if (!event)
0052     return false;
0053
0054   if (event->getReadAttempts() > 0)
0055     return false;  // record already tried
0056
0057   const lhef::HEPEUP &hepeup = *event->getHEPEUP();
0058
0059   if (!hepeup.NUP)
0060     return false;
0061
0062   setProcess(hepeup.IDPRUP, hepeup.XWGTUP, hepeup.SCALUP, hepeup.AQEDUP, hepeup.AQCDUP);
0063
0064   const std::vector<float> &scales = event->scales();
0065
0066   unsigned int iscale = 0;
0067   for (int i = 0; i < hepeup.NUP; i++) {
0068     //retrieve scale corresponding to each particle
0069     double scalein = -1.;
0070
0071     //handle clustering scales if present,
0072     //applies to outgoing partons only
0073     if (setScalesFromLHEF_ && !scales.empty() && hepeup.ISTUP[i] == 1) {
0074       if (iscale >= scales.size()) {
0075         edm::LogError("InvalidLHEInput") << "Pythia8 requires"
0076                                          << "cluster scales for all outgoing partons or for none" << std::endl;
0077       }
0078       scalein = scales[iscale];
0079       ++iscale;
0080     }
0081
0082     addParticle(hepeup.IDUP[i],
0083                 hepeup.ISTUP[i],
0084                 hepeup.MOTHUP[i].first,
0085                 hepeup.MOTHUP[i].second,
0086                 hepeup.ICOLUP[i].first,
0087                 hepeup.ICOLUP[i].second,
0088                 hepeup.PUP[i][0],
0089                 hepeup.PUP[i][1],
0090                 hepeup.PUP[i][2],
0091                 hepeup.PUP[i][3],
0092                 hepeup.PUP[i][4],
0093                 hepeup.VTIMUP[i],
0094                 hepeup.SPINUP[i],
0095                 scalein);
0096   }
0097
0098   if (!infoPtr->eventAttributes) {
0099     fEvAttributes->clear();
0100     infoPtr->eventAttributes = fEvAttributes;
0101   } else {
0102     infoPtr->eventAttributes = fEvAttributes;  // make sure still there
0103     infoPtr->eventAttributes->clear();
0104   }
0105
0106   //fill parton multiplicities if available
0107   int npLO = event->npLO();
0108   int npNLO = event->npNLO();
0109
0110   //default value of -99 indicates tags were not present in the original LHE file
0111   //don't pass to pythia in this case to emulate pythia internal lhe reader behaviour
0112   if (npLO != -99) {
0113     char buffer[100];
0114     snprintf(buffer, 100, "%i", npLO);
0115     (*infoPtr->eventAttributes)["npLO"] = buffer;
0116   }
0117   if (npNLO != -99) {
0118     char buffer[100];
0119     snprintf(buffer, 100, "%i", npNLO);
0120     (*infoPtr->eventAttributes)["npNLO"] = buffer;
0121   }
0122
0123   //add #rwgt info from comments
0124   const std::vector<std::string> &comments = event->getComments();
0125   for (unsigned i = 0; i < comments.size(); i++) {
0126     if (comments[i].rfind("#rwgt", 0) == 0)
0127       (*infoPtr->eventAttributes)["#rwgt"] = comments[i];
0128   }
0129
0130   const lhef::LHEEvent::PDF *pdf = event->getPDF();
0131   if (pdf) {
0132     this->setPdf(pdf->id.first,
0133                  pdf->id.second,
0134                  pdf->x.first,
0135                  pdf->x.second,
0136                  pdf->scalePDF,
0137                  pdf->xPDF.first,
0138                  pdf->xPDF.second,
0139                  true);
0140   } else {
0141     this->setPdf(hepeup.IDUP[0],
0142                  hepeup.IDUP[1],
0143                  hepeup.PUP[0][3] / runInfo->getHEPRUP()->EBMUP.first,
0144                  hepeup.PUP[1][3] / runInfo->getHEPRUP()->EBMUP.second,
0145                  0.,
0146                  0.,
0147                  0.,
0148                  false);
0149   }
0150
0151   event->attempted();
0152
0153   return true;
0154 }