Py8GunBase.cc

Go to the documentation of this file.

#include "GeneratorInterface/Pythia8Interface/interface/Py8GunBase.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Concurrency/interface/SharedResourceNames.h"

// EvtGen plugin
//
#include "Pythia8Plugins/EvtGen.h"

using namespace Pythia8;

const std::vector<std::string> gen::Py8GunBase::p8SharedResources = {edm::SharedResourceNames::kPythia8};

namespace gen {

  Py8GunBase::Py8GunBase(edm::ParameterSet const& ps) : Py8InterfaceBase(ps) {
    // PGun specs
    //
    edm::ParameterSet pgun_params = ps.getParameter<edm::ParameterSet>("PGunParameters");

    // although there's the method ParameterSet::empty(),
    // it looks like it's NOT even necessary to check if it is,
    // before trying to extract parameters - if it is empty,
    // the default values seem to be taken
    //
    fPartIDs = pgun_params.getParameter<std::vector<int> >("ParticleID");
    fMinPhi = pgun_params.getParameter<double>("MinPhi");  // ,-3.14159265358979323846);
    fMaxPhi = pgun_params.getParameter<double>("MaxPhi");  // , 3.14159265358979323846);
  }

  // specific to Py8GunHad !!!
  //
  bool Py8GunBase::initializeForInternalPartons() {
    // NO MATTER what was this setting below, override it before init
    // - this is essencial for the PGun mode

    // Key requirement: switch off ProcessLevel, and thereby also PartonLevel.
    fMasterGen->readString("ProcessLevel:all = off");
    fMasterGen->readString("ProcessLevel::resonanceDecays=on");
    fMasterGen->init();

    // init decayer
    fDecayer->readString("ProcessLevel:all = off");  // Same trick!
    fDecayer->readString("ProcessLevel::resonanceDecays=on");
    fDecayer->init();

    if (useEvtGen) {
      edm::LogInfo("Pythia8Interface") << "Creating and initializing pythia8 EvtGen plugin";
      evtgenDecays.reset(new EvtGenDecays(fMasterGen.get(), evtgenDecFile, evtgenPdlFile));
      for (unsigned int i = 0; i < evtgenUserFiles.size(); i++)
        evtgenDecays->readDecayFile(evtgenUserFiles.at(i));
    }

    return true;
  }

  bool Py8GunBase::residualDecay() {
    Event* pythiaEvent = &(fMasterGen->event);

    int NPartsBeforeDecays = pythiaEvent->size() - 1;  // do NOT count the very 1st "system" particle
                                                       // in Pythia8::Event record; it does NOT even
                                                       // get translated by the HepMCInterface to the
                                                       // HepMC::GenEvent record !!!
    int NPartsAfterDecays = event().get()->particles_size();

    if (NPartsAfterDecays == NPartsBeforeDecays)
      return true;

    bool result = true;

    for (int ipart = NPartsAfterDecays; ipart > NPartsBeforeDecays; ipart--) {
      HepMC::GenParticle* part = event().get()->barcode_to_particle(ipart);

      if (part->status() == 1 && (fDecayer->particleData).canDecay(part->pdg_id())) {
        fDecayer->event.reset();
        Particle py8part(part->pdg_id(),
                         93,
                         0,
                         0,
                         0,
                         0,
                         0,
                         0,
                         part->momentum().x(),
                         part->momentum().y(),
                         part->momentum().z(),
                         part->momentum().t(),
                         part->generated_mass());
        HepMC::GenVertex* ProdVtx = part->production_vertex();
        py8part.vProd(
            ProdVtx->position().x(), ProdVtx->position().y(), ProdVtx->position().z(), ProdVtx->position().t());
        py8part.tau((fDecayer->particleData).tau0(part->pdg_id()));
        fDecayer->event.append(py8part);
        int nentries = fDecayer->event.size();
        if (!fDecayer->event[nentries - 1].mayDecay())
          continue;
        fDecayer->next();
        int nentries1 = fDecayer->event.size();
        if (nentries1 <= nentries)
          continue;  //same number of particles, no decays...

        part->set_status(2);

        result = toHepMC.fill_next_event(*(fDecayer.get()), event().get(), -1, true, part);
      }
    }

    return result;
  }

  void Py8GunBase::finalizeEvent() {
    //******** Verbosity ********

    if (maxEventsToPrint > 0 && (pythiaPylistVerbosity || pythiaHepMCVerbosity || pythiaHepMCVerbosityParticles)) {
      maxEventsToPrint--;
      if (pythiaPylistVerbosity) {
        fMasterGen->info.list();
        fMasterGen->event.list();
      }

      if (pythiaHepMCVerbosity) {
        std::cout << "Event process = " << fMasterGen->info.code() << "\n"
                  << "----------------------" << std::endl;
        event()->print();
      }
      if (pythiaHepMCVerbosityParticles) {
        std::cout << "Event process = " << fMasterGen->info.code() << "\n"
                  << "----------------------" << std::endl;
        ascii_io->write_event(event().get());
      }
    }

    return;
  }

  void Py8GunBase::statistics() {
    fMasterGen->stat();

    double xsec = fMasterGen->info.sigmaGen();  // cross section in mb
    xsec *= 1.0e9;                              // translate to pb (CMS/Gen "convention" as of May 2009)
    runInfo().setInternalXSec(xsec);
    return;
  }

  void Py8GunBase::evtGenDecay() {
    if (evtgenDecays.get())
      evtgenDecays->decay();
  }

}  // namespace gen