diff --git a/GeneratorInterface/CascadeInterface/BuildFile.xml b/GeneratorInterface/CascadeInterface/BuildFile.xml deleted file mode 100644 index 8a2ad4f32ea62..0000000000000 --- a/GeneratorInterface/CascadeInterface/BuildFile.xml +++ /dev/null @@ -1,8 +0,0 @@ - - - - - - - - diff --git a/GeneratorInterface/CascadeInterface/plugins/BuildFile.xml b/GeneratorInterface/CascadeInterface/plugins/BuildFile.xml deleted file mode 100644 index 767782a893c5d..0000000000000 --- a/GeneratorInterface/CascadeInterface/plugins/BuildFile.xml +++ /dev/null @@ -1,13 +0,0 @@ - - - - - - - - - - - - - diff --git a/GeneratorInterface/CascadeInterface/plugins/Cascade2GeneratorFilter.cc b/GeneratorInterface/CascadeInterface/plugins/Cascade2GeneratorFilter.cc deleted file mode 100644 index 92a4d57c42d04..0000000000000 --- a/GeneratorInterface/CascadeInterface/plugins/Cascade2GeneratorFilter.cc +++ /dev/null @@ -1,10 +0,0 @@ -#include "GeneratorInterface/Core/interface/GeneratorFilter.h" -#include "GeneratorInterface/ExternalDecays/interface/ExternalDecayDriver.h" -#include "Cascade2Hadronizer.h" - -namespace gen { - typedef edm::GeneratorFilter Cascade2GeneratorFilter; -} - -using gen::Cascade2GeneratorFilter; -DEFINE_FWK_MODULE(Cascade2GeneratorFilter); diff --git a/GeneratorInterface/CascadeInterface/plugins/Cascade2Hadronizer.cc b/GeneratorInterface/CascadeInterface/plugins/Cascade2Hadronizer.cc deleted file mode 100644 index a547b2b47e18b..0000000000000 --- a/GeneratorInterface/CascadeInterface/plugins/Cascade2Hadronizer.cc +++ /dev/null @@ -1,812 +0,0 @@ -#include "GeneratorInterface/CascadeInterface/plugins/Cascade2Hadronizer.h" - -#include - -#include "HepMC/GenEvent.h" -#include "HepMC/PdfInfo.h" -#include "HepMC/PythiaWrapper6_4.h" -#include "GeneratorInterface/CascadeInterface/plugins/CascadeWrapper.h" - -#include "HepMC/HEPEVT_Wrapper.h" -#include "HepMC/IO_HEPEVT.h" -#include "HepMC/IO_GenEvent.h" - -#include "FWCore/Concurrency/interface/SharedResourceNames.h" -#include "FWCore/MessageLogger/interface/MessageLogger.h" -#include "GeneratorInterface/Core/interface/FortranCallback.h" -#include "GeneratorInterface/Core/interface/FortranInstance.h" - -HepMC::IO_HEPEVT hepevtio; - -#include "HepPID/ParticleIDTranslations.hh" - -#include "FWCore/Utilities/interface/EDMException.h" - -//-- Pythia6 routines and functionalities to pass around Pythia6 params - -#include "GeneratorInterface/Pythia6Interface/interface/Pythia6Service.h" -#include "GeneratorInterface/Pythia6Interface/interface/Pythia6Declarations.h" - -#include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h" - -#include "CLHEP/Random/RandomEngine.h" - -using namespace edm; -using namespace std; - -#define debug 0 - -static CLHEP::HepRandomEngine* cascade2RandomEngine; - -extern "C" { - -double dcasrn_(int* idummy) { - static int call = 0; - - double rdm_nb = cascade2RandomEngine->flat(); - - if (debug && ++call < 100) - cout << "dcasrn from c++, call: " << call << " random number: " << rdm_nb << endl; - - return rdm_nb; -} -} - -namespace gen { - - class Pythia6ServiceWithCallback : public Pythia6Service { - public: - Pythia6ServiceWithCallback(const edm::ParameterSet& pset) : Pythia6Service(pset) {} - - private: - void upInit() override { FortranCallback::getInstance()->fillHeader(); } - - void upEvnt() override { FortranCallback::getInstance()->fillEvent(); } - - bool upVeto() override { - bool veto = false; - if (!hepeup_.nup) - veto = true; //-- LHE Common Blocks - return (veto); - } - }; - - static struct { - int n, npad, k[5][pyjets_maxn]; - double p[5][pyjets_maxn], v[5][pyjets_maxn]; - } pyjets_local; - - const std::vector Cascade2Hadronizer::theSharedResources = {edm::SharedResourceNames::kPythia6, - gen::FortranInstance::kFortranInstance}; - - Cascade2Hadronizer::Cascade2Hadronizer(edm::ParameterSet const& pset) - : BaseHadronizer(pset), - fPy6Service(new Pythia6ServiceWithCallback(pset)), //-- this will store py6 parameters for further settings - - fComEnergy(pset.getParameter("comEnergy")), - fextCrossSection(pset.getUntrackedParameter("crossSection", -1.)), - fextCrossSectionError(pset.getUntrackedParameter("crossSectionError", -1.)), - fFilterEfficiency(pset.getUntrackedParameter("filterEfficiency", -1.)), - - fMaxEventsToPrint(pset.getUntrackedParameter("maxEventsToPrint", 0)), - fHepMCVerbosity(pset.getUntrackedParameter("pythiaHepMCVerbosity", false)), - fPythiaListVerbosity(pset.getUntrackedParameter("pythiaPylistVerbosity", 0)), - - fDisplayPythiaBanner(pset.getUntrackedParameter("displayPythiaBanner", false)), - fDisplayPythiaCards(pset.getUntrackedParameter("displayPythiaCards", false)) { - fParameters = pset.getParameter("Cascade2Parameters"); - - fConvertToPDG = false; - if (pset.exists("doPDGConvert")) - fConvertToPDG = pset.getParameter("doPDGConvert"); - - //-- silence Pythia6 banner printout, unless display requested - - if (!fDisplayPythiaBanner) { - if (!call_pygive("MSTU(12)=12345")) { - throw edm::Exception(edm::errors::Configuration, "PythiaError") << " Pythia did not accept MSTU(12)=12345"; - } - } - - //-- silence printouts from PYGIVE, unless display requested - - if (!fDisplayPythiaCards) { - if (!call_pygive("MSTU(13)=0")) { - throw edm::Exception(edm::errors::Configuration, "PythiaError") << " Pythia did not accept MSTU(13)=0"; - } - } - - //-- tmp stuff to deal with EvtGen corrupting pyjets - flushTmpStorage(); - } - - Cascade2Hadronizer::~Cascade2Hadronizer() { - if (fPy6Service != nullptr) - delete fPy6Service; - } - - void Cascade2Hadronizer::doSetRandomEngine(CLHEP::HepRandomEngine* v) { - cascade2RandomEngine = v; - fPy6Service->setRandomEngine(v); - } - - void Cascade2Hadronizer::flushTmpStorage() { - pyjets_local.n = 0; - pyjets_local.npad = 0; - - for (int ip = 0; ip < pyjets_maxn; ip++) { - for (int i = 0; i < 5; i++) { - pyjets_local.k[i][ip] = 0; - pyjets_local.p[i][ip] = 0.; - pyjets_local.v[i][ip] = 0.; - } - } - return; - } - - void Cascade2Hadronizer::fillTmpStorage() { - pyjets_local.n = pyjets.n; - pyjets_local.npad = pyjets.npad; - - for (int ip = 0; ip < pyjets_maxn; ip++) { - for (int i = 0; i < 5; i++) { - pyjets_local.k[i][ip] = pyjets.k[i][ip]; - pyjets_local.p[i][ip] = pyjets.p[i][ip]; - pyjets_local.v[i][ip] = pyjets.v[i][ip]; - } - } - return; - } - - void Cascade2Hadronizer::finalizeEvent() { - HepMC::PdfInfo pdf; - - //-- filling factorization "Q scale" now! pthat moved to binningValues() - - if (event()->signal_process_id() <= 0) - event()->set_signal_process_id(pypars.msti[0]); - if (event()->event_scale() <= 0) - event()->set_event_scale(pypars.pari[22]); - if (event()->alphaQED() <= 0) - event()->set_alphaQED(pyint1.vint[56]); - if (event()->alphaQCD() <= 0) - event()->set_alphaQCD(pyint1.vint[57]); - - //-- get pdf info directly from Pythia6 and set it up into HepMC::GenEvent - //-- S. Mrenna: Prefer vint block - - if (pdf.id1() <= 0) - pdf.set_id1(pyint1.mint[14] == 21 ? 0 : pyint1.mint[14]); - if (pdf.id2() <= 0) - pdf.set_id2(pyint1.mint[15] == 21 ? 0 : pyint1.mint[15]); - if (pdf.x1() <= 0) - pdf.set_x1(pyint1.vint[40]); - if (pdf.x2() <= 0) - pdf.set_x2(pyint1.vint[41]); - if (pdf.pdf1() <= 0) - pdf.set_pdf1(pyint1.vint[38] / pyint1.vint[40]); - if (pdf.pdf2() <= 0) - pdf.set_pdf2(pyint1.vint[39] / pyint1.vint[41]); - if (pdf.scalePDF() <= 0) - pdf.set_scalePDF(pyint1.vint[50]); - - event()->set_pdf_info(pdf); - - if (debug) { - cout << "standard Py6 event weight: pyint1.vint[96]: " << pyint1.vint[96] << endl; - cout << "event weight returned by PYEVWT: 1./(pyint1.vint[98]): " << 1. / (pyint1.vint[98]) << endl; - } - - //-- this is "standard" Py6 event weight (corresponds to PYINT1/VINT(97)) - // event()->weights().push_back(pyint1.vint[96]); - - //-- this is event weight as 1./VINT(99) (PYINT1/VINT(99) is returned by the PYEVWT) - // event()->weights().push_back(1./(pyint1.vint[98])); - - //-- all cascade events have weight = 1 - event()->weights().push_back(1.); - - //-- now create the GenEventInfo product from the GenEvent and fill the missing pieces - - eventInfo() = std::make_unique(event().get()); - - //-- in Pythia6 pthat is used to subdivide samples into different bins - //-- in LHE mode the binning is done by the external ME generator - //-- which is likely not pthat, so only filling it for Py6 internal mode - - eventInfo()->setBinningValues(vector(1, pypars.pari[16])); - - //-- here we treat long-lived particles - - if (pydat1.mstj[21] == 3 || pydat1.mstj[21] == 4) - imposeProperTime(); - - //-- convert particle IDs Py6 -> PDG (if requested) - - if (fConvertToPDG) { - for (HepMC::GenEvent::particle_iterator part = event()->particles_begin(); part != event()->particles_end(); - ++part) { - (*part)->set_pdg_id(HepPID::translatePythiatoPDT((*part)->pdg_id())); - } - } - - //-- service printouts, if requested - if (fMaxEventsToPrint > 0) { - fMaxEventsToPrint--; - if (fPythiaListVerbosity) - call_pylist(fPythiaListVerbosity); - if (fHepMCVerbosity) { - cout << "Event process = " << pypars.msti[0] << endl << "----------------------" << endl; - event()->print(); - } - } - - return; - } - - bool Cascade2Hadronizer::hadronize() { return false; } - - bool Cascade2Hadronizer::generatePartonsAndHadronize() { - //-- grab Py6 instance - Pythia6Service::InstanceWrapper guard(fPy6Service); - - FortranCallback::getInstance()->resetIterationsPerEvent(); - - //-- skip event if py6 considers it bad - if (pyint1.mint[50] != 0) { - event().reset(); - return false; - } - - //-- generation of the event with CASCADE - call_event(); - - //-- pythia pyhepc routine converts common PYJETS in common HEPEVT - call_pyhepc(1); - - //-- delete the created event from memory - event().reset(hepevtio.read_next_event()); - - //-- this is to deal with post-gen tools & residualDecay() that may reuse PYJETS - flushTmpStorage(); - fillTmpStorage(); - - return true; - } - - bool Cascade2Hadronizer::decay() { return true; } - - bool Cascade2Hadronizer::residualDecay() { - //-- grab Py6 instance - Pythia6Service::InstanceWrapper guard(fPy6Service); - - int NPartsBeforeDecays = pyjets_local.n; - int NPartsAfterDecays = event().get()->particles_size(); - int barcode = NPartsAfterDecays; - - // JVY: well, in principle, it's not a 100% fair to go up to NPartsBeforeDecays, - // because Photos will attach gamma's to existing vertices, i.e. in the middle - // of the event rather than at the end; but this will only shift pointers down, - // so we'll be going again over a few "original" particle... - // in the alternative, we may go all the way up to the beginning of the event - // and re-check if anything remains to decay, that's fine even if it'll take - // some extra CPU... - - for (int ipart = NPartsAfterDecays; ipart > NPartsBeforeDecays; ipart--) { - HepMC::GenParticle* part = event().get()->barcode_to_particle(ipart); - int status = part->status(); - if (status != 1) - continue; // check only "stable" particles, - // as some undecayed may still be marked as such - int pdgid = part->pdg_id(); - int py6ptr = pycomp_(pdgid); - - if (pydat3.mdcy[0][py6ptr - 1] != 1) - continue; // particle is not expected to decay - - int py6id = HepPID::translatePDTtoPythia(pdgid); - - // first, will need to zero out, then fill up PYJETS - // I better do it directly (by hands) rather than via py1ent - // - to avoid (additional) mass smearing - - if (part->momentum().t() <= part->generated_mass()) { - continue; // e == m -> 0-momentum, nothing to decay... - } - - else { - pyjets.n = 0; - - for (int i = 0; i < 5; i++) { - pyjets.k[i][0] = 0; - pyjets.p[i][0] = 0.; - pyjets.v[i][0] = 0.; - } - - pyjets.k[0][0] = 1; - pyjets.k[1][0] = py6id; - pyjets.p[4][0] = part->generated_mass(); - pyjets.p[3][0] = part->momentum().t(); - pyjets.p[0][0] = part->momentum().x(); - pyjets.p[1][0] = part->momentum().y(); - pyjets.p[2][0] = part->momentum().z(); - HepMC::GenVertex* prod_vtx = part->production_vertex(); - if (!prod_vtx) - continue; // in principle, should never happen but... - pyjets.v[0][0] = prod_vtx->position().x(); - pyjets.v[1][0] = prod_vtx->position().y(); - pyjets.v[2][0] = prod_vtx->position().z(); - pyjets.v[3][0] = prod_vtx->position().t(); - pyjets.v[4][0] = 0.; - pyjets.n = 1; - pyjets.npad = pyjets_local.npad; - } - - //-- now call Py6 decay routine - - int parent = 1; // since we always pass to Py6 a single particle - pydecy_(parent); - - //-- now attach decay products to mother - - for (int iprt1 = 1; iprt1 < pyjets.n; iprt1++) { - part->set_status(2); - - HepMC::GenVertex* DecVtx = new HepMC::GenVertex( - HepMC::FourVector(pyjets.v[0][iprt1], pyjets.v[1][iprt1], pyjets.v[2][iprt1], pyjets.v[3][iprt1])); - DecVtx->add_particle_in(part); // this will cleanup end_vertex if exists, replace with the new one - // I presume (vtx) barcode will be given automatically - - HepMC::FourVector pmom(pyjets.p[0][iprt1], pyjets.p[1][iprt1], pyjets.p[2][iprt1], pyjets.p[3][iprt1]); - - int dstatus = 0; - if (pyjets.k[0][iprt1] >= 1 && pyjets.k[0][iprt1] <= 10) - dstatus = 1; - else if (pyjets.k[0][iprt1] >= 11 && pyjets.k[0][iprt1] <= 20) - dstatus = 2; - else if (pyjets.k[0][iprt1] >= 21 && pyjets.k[0][iprt1] <= 30) - dstatus = 3; - else if (pyjets.k[0][iprt1] >= 31 && pyjets.k[0][iprt1] <= 100) - dstatus = pyjets.k[0][iprt1]; - - HepMC::GenParticle* daughter = - new HepMC::GenParticle(pmom, HepPID::translatePythiatoPDT(pyjets.k[1][iprt1]), dstatus); - barcode++; - daughter->suggest_barcode(barcode); - DecVtx->add_particle_out(daughter); - - int iprt2; - for (iprt2 = iprt1 + 1; iprt2 < pyjets.n; iprt2++) { // the pointer is shifted by -1, c++ style - - if (pyjets.k[2][iprt2] != parent) { - parent = pyjets.k[2][iprt2]; - break; // another parent particle; reset & break the loop - } - - HepMC::FourVector pmomN(pyjets.p[0][iprt2], pyjets.p[1][iprt2], pyjets.p[2][iprt2], pyjets.p[3][iprt2]); - - dstatus = 0; - if (pyjets.k[0][iprt2] >= 1 && pyjets.k[0][iprt2] <= 10) - dstatus = 1; - else if (pyjets.k[0][iprt2] >= 11 && pyjets.k[0][iprt2] <= 20) - dstatus = 2; - else if (pyjets.k[0][iprt2] >= 21 && pyjets.k[0][iprt2] <= 30) - dstatus = 3; - else if (pyjets.k[0][iprt2] >= 31 && pyjets.k[0][iprt2] <= 100) - dstatus = pyjets.k[0][iprt2]; - - HepMC::GenParticle* daughterN = - new HepMC::GenParticle(pmomN, HepPID::translatePythiatoPDT(pyjets.k[1][iprt2]), dstatus); - barcode++; - daughterN->suggest_barcode(barcode); - DecVtx->add_particle_out(daughterN); - } //-- end iprt2 loop - - iprt1 = iprt2 - 1; // reset counter such that it doesn't go over the same child more than once - // don't forget to offset back into c++ counting, as it's already +1 forward - - event().get()->add_vertex(DecVtx); - } //-- end iprt1 loop - } //-- end loop over decay products - - // now restore the very original Py6 event record - - if (pyjets_local.n != pyjets.n) { - // restore pyjets to its state as it was before external decays - - // might have been jammed by action above or by py1ent calls in EvtGen - - pyjets.n = pyjets_local.n; - pyjets.npad = pyjets_local.npad; - for (int ip = 0; ip < pyjets_local.n; ip++) { - for (int i = 0; i < 5; i++) { - pyjets.k[i][ip] = pyjets_local.k[i][ip]; - pyjets.p[i][ip] = pyjets_local.p[i][ip]; - pyjets.v[i][ip] = pyjets_local.v[i][ip]; - } - } - } - - return true; - } - - bool Cascade2Hadronizer::readSettings(int key) { - //-- grab Py6 instance - Pythia6Service::InstanceWrapper guard(fPy6Service); - - fPy6Service->setGeneralParams(); - - if (key == 0) - fPy6Service->setCSAParams(); - fPy6Service->setSLHAParams(); - - return true; - } - - bool Cascade2Hadronizer::initializeForExternalPartons() { return false; } - - bool Cascade2Hadronizer::initializeForInternalPartons() { - //-- grab Py6 instance - Pythia6Service::InstanceWrapper guard(fPy6Service); - - //-- change standard parameters of JETSET/PYTHIA - replace call_pytcha() - - fPy6Service->setGeneralParams(); - - fPy6Service->setCSAParams(); - fPy6Service->setSLHAParams(); - fPy6Service->setPYUPDAParams(false); - - pythia6PrintParameters(); - - //-- mstu(8) is set to NMXHEP in this dummy call (version >=6.404) - call_pyhepc(1); - - //-- initialise random number generator: has been changed to be CMSSW compliant - //-- dcasrn overloaded: call to rluxgo not needed anymore - //-- call_rluxgo(4,314159265,0,0); - - //-- initialise CASCADE parameters (default values) - call_casini(); - - //-- Read the parameters and pass them to the common blocks - //-- call_steer(); - - //-- initialise CASCADE parameters (user values) - - //-- retrieve all the different sets - vector AllSets = fParameters.getParameter >("parameterSets"); - - //-- loop over the different sets - for (unsigned i = 0; i < AllSets.size(); ++i) { - string Set = AllSets[i]; - vector Para_Set = fParameters.getParameter >(Set); - - //-- loop over all the parameters and stop in case of mistake - for (vector::const_iterator itPara = Para_Set.begin(); itPara != Para_Set.end(); ++itPara) { - if (!cascadeReadParameters(*itPara)) { - throw edm::Exception(edm::errors::Configuration, "CascadeError") - << " Cascade did not accept the following parameter: \"" << *itPara << "\"" << endl; - } - } //-- end loop over all the parameters - } //-- end loop over the different sets - - cainpu.plepin = -fComEnergy / 2; - cainpu.ppin = fComEnergy / 2; - - cascadePrintParameters(); - - //-- change standard parameters of CASCADE - call_cascha(); - - //-- change standard parameters of JETSET/PYTHIA - //-- call_pytcha(); - - //-- set up for running CASCADE (integration of the cross-section) - call_cascade(); - - //-- print cross-section result from integration - call_caend(1); - - fPy6Service->setPYUPDAParams(true); - - fPy6Service->closeSLHA(); - - return true; - } - - bool Cascade2Hadronizer::declareStableParticles(const vector& _pdg) { - vector pdg = _pdg; - for (size_t i = 0; i < pdg.size(); i++) { - int PyID = HepPID::translatePDTtoPythia(pdg[i]); - int pyCode = pycomp_(PyID); - - if (pyCode > 0) { - ostringstream pyCard; - pyCard << "MDCY(" << pyCode << ",1)=0"; - //-- cout << "pdg= " << pdg[i] << " " << pyCard.str() << endl; - - call_pygive(pyCard.str()); - } - } - - return true; - } - - void Cascade2Hadronizer::imposeProperTime() { - // this is practically a copy/paste of the original code by J.Alcaraz, - // taken directly from PythiaSource - - int dumm = 0; - HepMC::GenEvent::vertex_const_iterator vbegin = event()->vertices_begin(); - HepMC::GenEvent::vertex_const_iterator vend = event()->vertices_end(); - HepMC::GenEvent::vertex_const_iterator vitr = vbegin; - - for (; vitr != vend; ++vitr) { - HepMC::GenVertex::particle_iterator pbegin = (*vitr)->particles_begin(HepMC::children); - HepMC::GenVertex::particle_iterator pend = (*vitr)->particles_end(HepMC::children); - HepMC::GenVertex::particle_iterator pitr = pbegin; - - for (; pitr != pend; ++pitr) { - if ((*pitr)->end_vertex()) - continue; - if ((*pitr)->status() != 1) - continue; - - int pdgcode = abs((*pitr)->pdg_id()); - // Do nothing if the particle is not expected to decay - if (pydat3.mdcy[0][pycomp_(pdgcode) - 1] != 1) - continue; - - double ctau = pydat2.pmas[3][pycomp_(pdgcode) - 1]; - HepMC::FourVector mom = (*pitr)->momentum(); - HepMC::FourVector vin = (*vitr)->position(); - double x = 0.; - double y = 0.; - double z = 0.; - double t = 0.; - bool decayInRange = false; - while (!decayInRange) { - double unif_rand = fPy6Service->call(pyr_, &dumm); - // Value of 0 is excluded, so following line is OK - double proper_length = -ctau * log(unif_rand); - double factor = proper_length / mom.m(); - x = vin.x() + factor * mom.px(); - y = vin.y() + factor * mom.py(); - z = vin.z() + factor * mom.pz(); - t = vin.t() + factor * mom.e(); - - // Decay must be happen outside a cylindrical region - if (pydat1.mstj[21] == 4) { - if (sqrt(x * x + y * y) > pydat1.parj[72] || fabs(z) > pydat1.parj[73]) - decayInRange = true; - // Decay must be happen outside a given sphere - } else if (pydat1.mstj[21] == 3) { - if (sqrt(x * x + y * y + z * z) > pydat1.parj[71]) - decayInRange = true; - } - // Decay is always OK otherwise - else { - decayInRange = true; - } - } - - HepMC::GenVertex* vdec = new HepMC::GenVertex(HepMC::FourVector(x, y, z, t)); - event()->add_vertex(vdec); - vdec->add_particle_in((*pitr)); - } - } - - return; - } - - void Cascade2Hadronizer::statistics() { - runInfo().setExternalXSecLO(GenRunInfoProduct::XSec(fextCrossSection, fextCrossSectionError)); - - runInfo().setExternalXSecNLO(GenRunInfoProduct::XSec(-1., -1.)); - - if (!runInfo().internalXSec()) { - double sigma_CMS = 0.001 * caeffic.avgi; - double error_CMS = 0.001 * caeffic.sd; - runInfo().setInternalXSec(GenRunInfoProduct::XSec(sigma_CMS, error_CMS)); - } - - //-- print out generated event summary - call_caend(2); - - //-- write out some information from Pythia to the screen - call_pystat(1); - - return; - } - - const char* Cascade2Hadronizer::classname() const { return "gen::Cascade2Hadronizer"; } - - //-- Read the parameters and pass them to the common blocks - - bool Cascade2Hadronizer::cascadeReadParameters(const string& ParameterString) { - bool accepted = true; - - if (!strncmp(ParameterString.c_str(), "KE", 2)) - caluco.ke = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "IRES(1)", 7)) - capar6.ires[0] = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "KP", 2)) - caluco.kp = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "IRES(2)", 7)) - capar6.ires[1] = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "NFRAG", 5)) - cainpu.nfrag = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "IPST", 4)) - cashower.ipst = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "IPSIPOL", 7)) - jpsi.ipsipol = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - //-- from version 2.2.03 on - // else if(!strncmp(ParameterString.c_str(),"I23S",4)) - // jpsi.i23s = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]); - - else if (!strncmp(ParameterString.c_str(), "IFPS", 4)) - cainpu.ifps = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "ITIMSHR", 7)) - casshwr.itimshr = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "IRUNAEM", 7)) - capar1.irunaem = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "IRUNA", 5)) - capar1.iruna = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "IQ2", 3)) - capar1.iq2 = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "IPRO", 4)) - capar1.ipro = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "NFLAV", 5)) - caluco.nflav = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "INTER", 5)) - cainpu.inter = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "IHFLA", 5)) - cahflav.ihfla = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "IRPA", 4)) - cascol.irpa = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "IRPB", 4)) - cascol.irpb = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "IRPC", 4)) - cascol.irpc = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "ICCFM", 5)) - casshwr.iccfm = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "IFINAL", 6)) - cainpu.ifinal = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "IGLU", 4)) - cagluon.iglu = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "IRspl", 5)) - casprre.irspl = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "PT2CUT", 6)) - captcut.pt2cut[capar1.ipro - 1] = atof(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "NCB", 3)) - integr.ncb = atoi(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "ACC1", 4)) - integr.acc1 = atof(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "ACC2", 4)) - integr.acc2 = atof(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "SCALFAF", 7)) - scalf.scalfaf = atof(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else if (!strncmp(ParameterString.c_str(), "SCALFA", 6)) - scalf.scalfa = atof(&ParameterString[strcspn(ParameterString.c_str(), "=") + 1]); - - else - accepted = false; - - return accepted; - } - - void Cascade2Hadronizer::cascadePrintParameters() { - cout << "" << endl; - cout << "----------------------------------" << endl; - cout << "---- Cascade Parameters ----" << endl; - cout << "----------------------------------" << endl; - cout << "" << endl; - - cout << "computation switch: " << endl; - cout << "number of calls per iteration for bases: " << integr.ncb << endl; - cout << "relative precision for grid optimisation: " << integr.acc1 << " %" << endl; - cout << "relative precision for integration: " << integr.acc2 << " %" << endl; - cout << "" << endl; - - cout << "kinematics: " << endl; - cout << "flavour code of beam 1: " << caluco.ke << endl; - cout << "direct or resolved particle 1: " << capar6.ires[0] << endl; - cout << "pz of incoming beam 1: " << cainpu.plepin << " GeV" << endl; - cout << "flavour code of beam 2: " << caluco.kp << endl; - cout << "direct or resolved particle 2: " << capar6.ires[1] << endl; - cout << "pz of incoming beam 2: " << cainpu.ppin << " GeV" << endl; - cout << "number of active flavors: " << caluco.nflav << endl; - cout << "" << endl; - - cout << "hard subprocess selection: " << endl; - cout << "hard subprocess number (IPRO): " << capar1.ipro << endl; - cout << "IPRO = 10, switch to select QCD process g* g* -> q qbar: " << cascol.irpa << endl; - cout << "IPRO = 10, switch to select QCD process g* g -> g g: " << cascol.irpb << endl; - cout << "IPRO = 10, switch to select QCD process g* q -> g q: " << cascol.irpc << endl; - cout << "flavor of heavy quark produced (IPRO = 11, 504, 514): " << cahflav.ihfla << endl; - // cout<<"select vector meson state: "< -#include -#include - -#include "FWCore/ParameterSet/interface/ParameterSetfwd.h" -#include "GeneratorInterface/Core/interface/ParameterCollector.h" -#include "GeneratorInterface/Core/interface/BaseHadronizer.h" - -namespace HepMC { - class GenEvent; -} - -namespace CLHEP { - class HepRandomEngine; -} - -namespace gen { - - class Pythia6Service; - - class Cascade2Hadronizer : public BaseHadronizer { - public: - Cascade2Hadronizer(edm::ParameterSet const& ps); - ~Cascade2Hadronizer() override; - - bool readSettings(int); - bool initializeForExternalPartons(); //-- initializer for the LHE input - bool initializeForInternalPartons(); - - //-- Read the parameters and pass them to the common blocks - bool cascadeReadParameters(const std::string& ParameterString); - void cascadePrintParameters(); - void pythia6PrintParameters(); - bool declareStableParticles(const std::vector&); - bool declareSpecialSettings(const std::vector&) { return true; } - void statistics(); - - bool generatePartonsAndHadronize(); - bool hadronize(); //-- hadronizer for the LHE input - bool decay(); - bool residualDecay(); - void finalizeEvent(); - - const char* classname() const; - - private: - void doSetRandomEngine(CLHEP::HepRandomEngine* v) override; - std::vector const& doSharedResources() const override { return theSharedResources; } - - static const std::vector theSharedResources; - - //-- methods - - void flushTmpStorage(); - void fillTmpStorage(); - void imposeProperTime(); //-- to correctly treat particle decay - - //-- data members - - edm::ParameterSet fParameters; - - Pythia6Service* fPy6Service; - - double - fComEnergy; //-- irrelevant for setting py6 as hadronizer (or if anything, it should be picked up from LHERunInfoProduct) - double fextCrossSection; - double fextCrossSectionError; - double fFilterEfficiency; - - unsigned int fMaxEventsToPrint; - bool fHepMCVerbosity; - unsigned int fPythiaListVerbosity; //-- p6 specific - - bool fDisplayPythiaBanner; //-- p6 specific - bool fDisplayPythiaCards; //-- p6 specific - - bool fConvertToPDG; //-- conversion of Py6 PID's into PDG convention - }; -} // namespace gen - -#endif diff --git a/GeneratorInterface/CascadeInterface/plugins/CascadeWrapper.h b/GeneratorInterface/CascadeInterface/plugins/CascadeWrapper.h deleted file mode 100644 index 06e3189921059..0000000000000 --- a/GeneratorInterface/CascadeInterface/plugins/CascadeWrapper.h +++ /dev/null @@ -1,159 +0,0 @@ -#ifndef CASCADE_WRAPPER_H -#define CASCADE_WRAPPER_H - -//-- the pyhepc routine used by Pythia to fill the HEPEVT common block -//-- is using double precision and 4000 entries - -#include -#include - -//-- CASCADE Common Block Declarations - -extern "C" { -extern struct { int ke, kp, keb, kph, kgl, kpa, nflav; } caluco_; -} - -#define caluco caluco_ - -extern "C" { -extern struct { int lst[30], ires[2]; } capar6_; -} - -#define capar6 capar6_ - -extern "C" { -extern struct { - double plepin, ppin; - int nfrag, ilepto, ifps, ihf, inter, isemih, ifinal; -} cainpu_; -} - -#define cainpu cainpu_ - -extern "C" { -extern struct { int ipst; } cashower_; -} - -#define cashower cashower_ - -extern "C" { -extern struct { - int ipsipol, ipsiel1, ipsiel2; - //-- int i23s; //-- from version 2.2.03 on -} jpsi_; -} - -#define jpsi jpsi_ - -extern "C" { -extern struct { int iorder, itimshr, iccfm; } casshwr_; -} - -#define casshwr casshwr_ - -extern "C" { -extern struct { int ipro, iruna, iq2, irunaem; } capar1_; -} - -#define capar1 capar1_ - -extern "C" { -extern struct { int ihfla, kpsi, kchi; } cahflav_; -} - -#define cahflav cahflav_ - -extern "C" { -extern struct { int icolora, irespro, irpa, irpb, irpc, irpd, irpe, irpf, irpg; } cascol_; -} - -#define cascol cascol_ - -extern "C" { -extern struct { int iglu; } cagluon_; -} - -#define cagluon cagluon_ - -extern "C" { -extern struct { int irspl; } casprre_; -} - -#define casprre casprre_ - -extern "C" { -extern struct { double pt2cut[1000]; } captcut_; -} - -#define captcut captcut_ - -extern "C" { -extern struct { - double acc1, acc2; - int iint, ncb; -} integr_; -} - -#define integr integr_ - -extern "C" { -extern struct { double scalfa, scalfaf; } scalf_; -} - -#define scalf scalf_ - -extern "C" { -extern struct { char pdfpath[512]; } caspdf_; -} - -#define caspdf caspdf_ - -extern "C" { -extern struct { - double avgi, sd; - int nin, nout; -} caeffic_; -} - -#define caeffic caeffic_ - -//-- CASCADE routines declarations - -extern "C" { -void casini_(); -void steer_(); -void cascha_(); -void cascade_(); -void caend_(int* mode); -void event_(); -// void rluxgo_(int* mode1,int* mode2,int* mode3,int* mode4); -double dcasrn_(int* idummy); -} - -inline void call_casini() { casini_(); } -inline void call_steer() { steer_(); } -inline void call_cascha() { cascha_(); } -inline void call_cascade() { cascade_(); } -inline void call_caend(int mode) { caend_(&mode); } -inline void call_event() { event_(); } -// inline void call_rluxgo(int mode1, int mode2, int mode3, int mode4) { rluxgo_(&mode1, &mode2, &mode3, &mode4); } -// inline double call_dcasrn() { return(dcasrn_()); } - -//-- PYTHIA Common Block Declarations - -//-- PYTHIA routines declarations - -extern "C" { -void pytcha_(); -void pyedit_(int* mode); -} - -inline void call_pytcha() { pytcha_(); } -inline void call_pyedit(int mode) { pyedit_(&mode); } - -//inline void call_pyhepc(int mode) { pyhepc(&mode); } -//inline void call_pylist(int mode) { pylist(&mode); } -//inline void call_pystat(int mode) { pystat(&mode); } -//inline void call_pyevnt() { pyevnt(); } - -#endif // CASCADE_WRAPPER_H diff --git a/GeneratorInterface/CascadeInterface/python/Cascade2Parameters_cfi.py b/GeneratorInterface/CascadeInterface/python/Cascade2Parameters_cfi.py deleted file mode 100644 index e6073364ac2d6..0000000000000 --- a/GeneratorInterface/CascadeInterface/python/Cascade2Parameters_cfi.py +++ /dev/null @@ -1,35 +0,0 @@ -import FWCore.ParameterSet.Config as cms - - -Cascade2Parameters = cms.PSet( - parameterSets = cms.vstring('CascadeSettings'), - CascadeSettings = cms.vstring('NCB = 50000 ! number of calls per iteration for bases', - 'ACC1 = 1.0 ! relative precision for grid optimisation', - 'ACC2 = 0.5 ! relative precision for integration', - 'KE = 2212 ! flavour code of beam1', - 'IRES(1) = 1 ! direct or resolved particle 1', - 'KP = 2212 ! flavour code of beam2', - 'IRES(2) = 1 ! direct or resolved particle 2', - 'NFLAV = 5 ! number of active flavors', - 'IPRO = 10 ! hard subprocess number', - 'IRPA = 1 ! IPRO = 10, switch to select QCD process g* g* -> q qbar', - 'IRPB = 1 ! IPRO = 10, switch to select QCD process g* g -> g g', - 'IRPC = 1 ! IPRO = 10, switch to select QCD process g* q -> g q', - 'IHFLA = 4 ! flavor of heavy quark produced (IPRO = 11, 504, 514)', - # 'I23S = 0 ! select vector meson state (0 = 1S, 2 = 2S, 3 = 3S from version 2.2.03 on)', - 'IPSIPOL = 0 ! use matrix element including J/psi (Upsilon) polarisation (1 = on, 0 = off)', - 'PT2CUT = 0.0 ! pt2 cut in ME for massless partons' - 'NFRAG = 1 ! switch for fragmentation (1 = on, 0 = off)', - 'IFPS = 3 ! switch for parton shower: (0 = off - 1 = initial - 2 = final - 3 = initial & final)', - 'ITIMSHR = 1 ! switch for time like parton shower in intial state (0 = off,1 = on)', - 'ICCFM = 1 ! select CCFM or DGLAP evolution (CCFM = 1, DGLAP = 0)' - 'IFINAL = 1 ! scale switch for final state parton shower', - 'SCALFAF = 1.0 ! scale factor for final state parton shower', - 'IRspl = 4 ! switch for proton remnant treatment', - 'IPST = 0 ! keep track of intermediate state in parton shower ', - 'INTER = 0 ! mode of interaction for ep (photon exchange, Z-echange (not implemented))', - 'IRUNAEM = 1 ! switch for running alphaem (0 = off,1 = on)', - 'IRUNA = 1 ! switch for running alphas (0 = off,1 = on)', - 'IQ2 = 3 ! scale in alphas', - 'SCALFA = 1.0 ! scale factor for scale in alphas', - 'IGLU = 1010 ! select uPDF')) diff --git a/GeneratorInterface/CascadeInterface/python/Cascade2_example_test_cff.py b/GeneratorInterface/CascadeInterface/python/Cascade2_example_test_cff.py deleted file mode 100644 index da520f0f9abcf..0000000000000 --- a/GeneratorInterface/CascadeInterface/python/Cascade2_example_test_cff.py +++ /dev/null @@ -1,32 +0,0 @@ -import FWCore.ParameterSet.Config as cms - -from GeneratorInterface.CascadeInterface.Cascade2Parameters_cfi import Cascade2Parameters as Cascade2ParametersRef - -source = cms.Source("EmptySource") - -generator = cms.EDFilter('Cascade2GeneratorFilter', - PythiaParameters = cms.PSet( - processParameters = cms.vstring('PMAS(4,1) = 1.6 ! charm mass', - 'PMAS(5,1) = 4.75 ! bottom mass', - 'PMAS(25,1) = 125. ! higgs mass', - 'PARU(112) = 0.2 ! lambda QCD set A0', - 'MSTU(111) = 1 ! = 0 : alpha_s is fixed at the value PARU(111), = 1 : first-order running alpha_s, = 2 : second-order running alpha_s', - 'MSTU(112) = 4 ! nr of flavours wrt lambda_QCD', - 'MSTU(113) = 3 ! min nr of flavours for alphas', - 'MSTU(114) = 5 ! max nr of flavours for alphas', - 'MSTJ(48) = 1 ! (D = 0) 0 = no max. angle, 1 = max angle def. in PARJ(85)'), - parameterSets = cms.vstring('processParameters') - ), - - comEnergy = cms.double(7000.0), - crossSection = cms.untracked.double(-1), - crossSectionError = cms.untracked.double(-1), - filterEfficiency = cms.untracked.double(1), - maxEventsToPrint = cms.untracked.int32(2), - pythiaHepMCVerbosity = cms.untracked.bool(False), - pythiaPylistVerbosity = cms.untracked.int32(0), - - Cascade2Parameters = Cascade2ParametersRef - ) - -ProductionFilterSequence = cms.Sequence(generator)