CMS 3D CMS Logo

Public Member Functions | Private Member Functions | Private Attributes

AlpgenSource Class Reference

Inheritance diagram for AlpgenSource:
edm::ExternalInputSource edm::ConfigurableInputSource edm::InputSource edm::ProductRegistryHelper

List of all members.

Public Member Functions

 AlpgenSource (const edm::ParameterSet &params, const edm::InputSourceDescription &desc)
 Constructor.
virtual ~AlpgenSource ()
 Destructor.

Private Member Functions

virtual void beginRun (edm::Run &run)
template<typename T >
getParameter (AlpgenHeader::Parameter index, const T &defValue) const
 Function to get parameter by name from AlpgenHeader, w/ default.
template<typename T >
getParameter (AlpgenHeader::Parameter index) const
 Function to get parameter by name from AlpgenHeader.
unsigned int processID () const
virtual bool produce (edm::Event &event)
bool readAlpgenEvent (lhef::HEPEUP &hepeup)
 Read an event and put it into the HEPEUP.
std::string slhaMassLine (int pdgId, AlpgenHeader::Masses mass, const std::string &comment) const

Private Attributes

std::string extraHeaderFileName_
 Name of the extra header file.
std::string extraHeaderName_
 Name given to the extra header.
std::string fileName_
 Name of the input file.
AlpgenHeader header
 Alpgen _unw.par file as an AlpgenHeader.
std::auto_ptr< std::ifstream > inputFile_
 Pointer to the input file.
LHERunInfoProduct::Header lheAlpgenUnwParHeader
 Alpgen _unw.par file as a LHE header.
unsigned long nEvent_
 Number of events.
unsigned long skipEvents_
 Number of events to skip.
bool writeAlpgenParFile
bool writeAlpgenWgtFile
 configuration flags
bool writeExtraHeader

Detailed Description

Definition at line 29 of file AlpgenSource.cc.


Constructor & Destructor Documentation

AlpgenSource::AlpgenSource ( const edm::ParameterSet params,
const edm::InputSourceDescription desc 
)

Constructor.

Definition at line 91 of file AlpgenSource.cc.

References LHERunInfoProduct::Header::addLine(), LHERunInfoProduct::Header::begin(), LHERunInfoProduct::Header::end(), Exception, fileName_, edm::ExternalInputSource::fileNames(), header, lheAlpgenUnwParHeader, and AlpgenHeader::parse().

                                                                  :
  edm::ExternalInputSource(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::InRun>();
  produces<LHEEventProduct>();
}
AlpgenSource::~AlpgenSource ( ) [virtual]

Destructor.

Definition at line 144 of file AlpgenSource.cc.

{
}

Member Function Documentation

void AlpgenSource::beginRun ( edm::Run run) [private, virtual]

Reimplemented from edm::ConfigurableInputSource.

Definition at line 158 of file AlpgenSource.cc.

References LHERunInfoProduct::Header::addLine(), indexGen::comments, AlpgenHeader::ebeam, lhef::HEPRUP::EBMUP, Exception, extraHeaderFileName_, extraHeaderName_, fileName_, header, lhef::HEPRUP::IDBMUP, lhef::HEPRUP::IDWTUP, AlpgenHeader::ih2, inputFile_, lheAlpgenUnwParHeader, lhef::HEPRUP::LPRUP, AlpgenHeader::mb, AlpgenHeader::mc, AlpgenHeader::mh, AlpgenHeader::mt, AlpgenHeader::mw, AlpgenHeader::mz, lhef::HEPRUP::PDFGUP, lhef::HEPRUP::PDFSUP, processID(), edm::Run::put(), lhef::HEPRUP::resize(), slhaMassLine(), writeAlpgenParFile, writeAlpgenWgtFile, writeExtraHeader, lhef::HEPRUP::XERRUP, lhef::HEPRUP::XMAXUP, AlpgenHeader::xsec, AlpgenHeader::xsecErr, and lhef::HEPRUP::XSECUP.

{
  // 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_.c_str());
  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::auto_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(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 T &  defValue 
) const [private]

Function to get parameter by name from AlpgenHeader, w/ default.

Definition at line 291 of file AlpgenSource.cc.

References header, AlpgenHeader::params, and pos.

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

Function to get parameter by name from AlpgenHeader.

Definition at line 277 of file AlpgenSource.cc.

References header, AlpgenHeader::parameterName(), AlpgenHeader::params, and pos.

{
  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);
}
unsigned int AlpgenSource::processID ( ) const [private]

The Alpgen process ID. This is defined as processID() = 100*X + 10*Y + Z, where = ihrd, Y = ihvy, Z = njets

Definition at line 302 of file AlpgenSource.cc.

References header, AlpgenHeader::ihrd, AlpgenHeader::ihvy, and AlpgenHeader::njets.

Referenced by beginRun(), and readAlpgenEvent().

{
  // 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::produce ( edm::Event event) [private, virtual]

Implements edm::ConfigurableInputSource.

Definition at line 403 of file AlpgenSource.cc.

References Exception, alpgen::fixEventHiggsTTbar(), alpgen::fixEventMultiBoson(), alpgen::fixEventSingleTop(), alpgen::fixEventTTbar(), alpgen::fixEventWZ(), header, AlpgenHeader::ihrd, inputFile_, AlpgenHeader::itopprc, AlpgenHeader::masses, AlpgenHeader::mb, nEvent_, AlpgenHeader::params, readAlpgenEvent(), and skipEvents_.

{
  // The LHE Event Record
  lhef::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;
  }

  // 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.

  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::auto_ptr<LHEEventProduct> lheEvent(new LHEEventProduct(hepeup));
  event.put(lheEvent);

  return true;
}
bool AlpgenSource::readAlpgenEvent ( lhef::HEPEUP hepeup) [private]

Read an event and put it into the HEPEUP.

Definition at line 311 of file AlpgenSource.cc.

References relval_parameters_module::energy, reco::flavour(), i, lhef::HEPEUP::ICOLUP, lhef::HEPEUP::IDPRUP, lhef::HEPEUP::IDUP, inputFile_, lhef::HEPEUP::ISTUP, python::rootplot::utilities::ls(), lhef::HEPEUP::MOTHUP, processID(), lhef::HEPEUP::PUP, lhef::HEPEUP::resize(), lhef::HEPEUP::SCALUP, mathSSE::sqrt(), and lhef::HEPEUP::XWGTUP.

Referenced by produce().

{
  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;
}
std::string AlpgenSource::slhaMassLine ( int  pdgId,
AlpgenHeader::Masses  mass,
const std::string &  comment 
) const [private]

Converts the AlpgenHeader::Masses to a std::string formatted as a slhaMassLine to facilitate passing them to Alpgen.

Definition at line 148 of file AlpgenSource.cc.

References header, and AlpgenHeader::masses.

Referenced by beginRun().

{
  std::ostringstream ss;
  ss << std::setw(9) << pdgId << "     " << std::scientific
     << std::setprecision(9) << header.masses[mass] << "   # "
     << comment << std::endl;
  return ss.str();
}

Member Data Documentation

std::string AlpgenSource::extraHeaderFileName_ [private]

Name of the extra header file.

Definition at line 80 of file AlpgenSource.cc.

Referenced by beginRun().

std::string AlpgenSource::extraHeaderName_ [private]

Name given to the extra header.

Definition at line 83 of file AlpgenSource.cc.

Referenced by beginRun().

std::string AlpgenSource::fileName_ [private]

Name of the input file.

Definition at line 63 of file AlpgenSource.cc.

Referenced by AlpgenSource(), and beginRun().

Alpgen _unw.par file as an AlpgenHeader.

Definition at line 77 of file AlpgenSource.cc.

Referenced by AlpgenSource(), beginRun(), getParameter(), processID(), produce(), and slhaMassLine().

std::auto_ptr<std::ifstream> AlpgenSource::inputFile_ [private]

Pointer to the input file.

Definition at line 66 of file AlpgenSource.cc.

Referenced by beginRun(), produce(), and readAlpgenEvent().

Alpgen _unw.par file as a LHE header.

Definition at line 75 of file AlpgenSource.cc.

Referenced by AlpgenSource(), and beginRun().

unsigned long AlpgenSource::nEvent_ [private]

Number of events.

Definition at line 72 of file AlpgenSource.cc.

Referenced by produce().

unsigned long AlpgenSource::skipEvents_ [private]

Number of events to skip.

Definition at line 69 of file AlpgenSource.cc.

Referenced by produce().

Definition at line 87 of file AlpgenSource.cc.

Referenced by beginRun().

configuration flags

Definition at line 86 of file AlpgenSource.cc.

Referenced by beginRun().

Definition at line 88 of file AlpgenSource.cc.

Referenced by beginRun().