Pythia8Interface/src/Py8GunBase.cc

0001 #include <memory>
0002
0003 #include "GeneratorInterface/Pythia8Interface/interface/Py8GunBase.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/Concurrency/interface/SharedResourceNames.h"
0006
0007 // EvtGen plugin
0008 //
0009 #include "Pythia8Plugins/EvtGen.h"
0010
0011 using namespace Pythia8;
0012
0013 const std::vector<std::string> gen::Py8GunBase::p8SharedResources = {edm::SharedResourceNames::kPythia8};
0014
0015 namespace gen {
0016
0017   Py8GunBase::Py8GunBase(edm::ParameterSet const& ps) : Py8InterfaceBase(ps) {
0018     // PGun specs
0019     //
0020     edm::ParameterSet pgun_params = ps.getParameter<edm::ParameterSet>("PGunParameters");
0021
0022     // although there's the method ParameterSet::empty(),
0023     // it looks like it's NOT even necessary to check if it is,
0024     // before trying to extract parameters - if it is empty,
0025     // the default values seem to be taken
0026     //
0027     fPartIDs = pgun_params.getParameter<std::vector<int> >("ParticleID");
0028     fMinPhi = pgun_params.getParameter<double>("MinPhi");  // ,-3.14159265358979323846);
0029     fMaxPhi = pgun_params.getParameter<double>("MaxPhi");  // , 3.14159265358979323846);
0030   }
0031
0032   // specific to Py8GunHad !!!
0033   //
0034   bool Py8GunBase::initializeForInternalPartons() {
0035     // NO MATTER what was this setting below, override it before init
0036     // - this is essencial for the PGun mode
0037
0038     // Key requirement: switch off ProcessLevel, and thereby also PartonLevel.
0039     fMasterGen->readString("ProcessLevel:all = off");
0040     fMasterGen->readString("ProcessLevel::resonanceDecays=on");
0041     fMasterGen->init();
0042
0043     // init decayer
0044     fDecayer->readString("ProcessLevel:all = off");  // Same trick!
0045     fDecayer->readString("ProcessLevel::resonanceDecays=on");
0046     fDecayer->init();
0047
0048     if (useEvtGen) {
0049       edm::LogInfo("Pythia8Interface") << "Creating and initializing pythia8 EvtGen plugin";
0050       evtgenDecays = std::make_shared<EvtGenDecays>(fMasterGen.get(), evtgenDecFile, evtgenPdlFile);
0051       for (unsigned int i = 0; i < evtgenUserFiles.size(); i++)
0052         evtgenDecays->readDecayFile(evtgenUserFiles.at(i));
0053     }
0054
0055     return true;
0056   }
0057
0058   bool Py8GunBase::residualDecay() {
0059     Event* pythiaEvent = &(fMasterGen->event);
0060
0061     int NPartsBeforeDecays = pythiaEvent->size() - 1;  // do NOT count the very 1st "system" particle
0062                                                        // in Pythia8::Event record; it does NOT even
0063                                                        // get translated by the HepMCInterface to the
0064                                                        // HepMC::GenEvent record !!!
0065     int NPartsAfterDecays = event().get()->particles_size();
0066
0067     if (NPartsAfterDecays == NPartsBeforeDecays)
0068       return true;
0069
0070     bool result = true;
0071
0072     for (int ipart = NPartsAfterDecays; ipart > NPartsBeforeDecays; ipart--) {
0073       HepMC::GenParticle* part = event().get()->barcode_to_particle(ipart);
0074
0075       if (part->status() == 1 && (fDecayer->particleData).canDecay(part->pdg_id())) {
0076         fDecayer->event.reset();
0077         Particle py8part(part->pdg_id(),
0078                          93,
0079                          0,
0080                          0,
0081                          0,
0082                          0,
0083                          0,
0084                          0,
0085                          part->momentum().x(),
0086                          part->momentum().y(),
0087                          part->momentum().z(),
0088                          part->momentum().t(),
0089                          part->generated_mass());
0090         HepMC::GenVertex* ProdVtx = part->production_vertex();
0091         py8part.vProd(
0092             ProdVtx->position().x(), ProdVtx->position().y(), ProdVtx->position().z(), ProdVtx->position().t());
0093         py8part.tau((fDecayer->particleData).tau0(part->pdg_id()));
0094         fDecayer->event.append(py8part);
0095         int nentries = fDecayer->event.size();
0096         if (!fDecayer->event[nentries - 1].mayDecay())
0097           continue;
0098         fDecayer->next();
0099         int nentries1 = fDecayer->event.size();
0100         if (nentries1 <= nentries)
0101           continue;  //same number of particles, no decays...
0102
0103         part->set_status(2);
0104
0105         result = toHepMC.fill_next_event(*(fDecayer.get()), event().get(), -1, true, part);
0106       }
0107     }
0108
0109     return result;
0110   }
0111
0112   void Py8GunBase::finalizeEvent() {
0113     //******** Verbosity ********
0114
0115     if (maxEventsToPrint > 0 && (pythiaPylistVerbosity || pythiaHepMCVerbosity || pythiaHepMCVerbosityParticles)) {
0116       maxEventsToPrint--;
0117       if (pythiaPylistVerbosity) {
0118         fMasterGen->info.list();
0119         fMasterGen->event.list();
0120       }
0121
0122       if (pythiaHepMCVerbosity) {
0123         std::cout << "Event process = " << fMasterGen->info.code() << "\n"
0124                   << "----------------------" << std::endl;
0125         event()->print();
0126       }
0127       if (pythiaHepMCVerbosityParticles) {
0128         std::cout << "Event process = " << fMasterGen->info.code() << "\n"
0129                   << "----------------------" << std::endl;
0130         ascii_io->write_event(event().get());
0131       }
0132     }
0133
0134     return;
0135   }
0136
0137   void Py8GunBase::statistics() {
0138     fMasterGen->stat();
0139
0140     double xsec = fMasterGen->info.sigmaGen();  // cross section in mb
0141     xsec *= 1.0e9;                              // translate to pb (CMS/Gen "convention" as of May 2009)
0142     runInfo().setInternalXSec(xsec);
0143     return;
0144   }
0145
0146   void Py8GunBase::evtGenDecay() {
0147     if (evtgenDecays.get())
0148       evtgenDecays->decay();
0149   }
0150
0151 }  // namespace gen