Pythia8Interface/src/Py8GunBase.cc

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