#include <Cascade2Hadronizer.h>

Inheritance diagram for gen::Cascade2Hadronizer:

Public Member Functions
	Cascade2Hadronizer (edm::ParameterSet const &ps)

void	cascadePrintParameters ()

bool	cascadeReadParameters (const std::string &ParameterString)

const char *	classname () const

bool	decay ()

bool	declareSpecialSettings (const std::vector< std::string > &)

bool	declareStableParticles (const std::vector< int > &)

void	finalizeEvent ()

bool	generatePartonsAndHadronize ()

bool	hadronize ()

bool	initializeForExternalPartons ()

bool	initializeForInternalPartons ()

void	pythia6PrintParameters ()

bool	readSettings (int)

bool	residualDecay ()

void	statistics ()

	~Cascade2Hadronizer () override

Public Member Functions inherited from gen::BaseHadronizer
	BaseHadronizer (edm::ParameterSet const &ps)

void	cleanLHE ()

void	generateLHE (edm::LuminosityBlock const &lumi, CLHEP::HepRandomEngine *rengine, unsigned int ncpu)

edm::Event &	getEDMEvent () const

std::unique_ptr< HepMC::GenEvent >	getGenEvent ()

std::unique_ptr< GenEventInfoProduct >	getGenEventInfo ()

virtual std::unique_ptr< GenLumiInfoHeader >	getGenLumiInfoHeader () const

GenRunInfoProduct &	getGenRunInfo ()

std::unique_ptr< lhef::LHEEvent >	getLHEEvent ()

const std::shared_ptr< lhef::LHERunInfo > &	getLHERunInfo () const

const std::string &	gridpackPath () const

int	randomIndex () const

const std::string &	randomInitConfigDescription () const

void	randomizeIndex (edm::LuminosityBlock const &lumi, CLHEP::HepRandomEngine *rengine)

void	resetEvent (std::unique_ptr< HepMC::GenEvent > event)

void	resetEventInfo (std::unique_ptr< GenEventInfoProduct > eventInfo)

virtual bool	select (HepMC::GenEvent *) const

void	setEDMEvent (edm::Event &event)

void	setLHEEvent (std::unique_ptr< lhef::LHEEvent > event)

void	setLHERunInfo (std::unique_ptr< lhef::LHERunInfo > runInfo)

void	setRandomEngine (CLHEP::HepRandomEngine *v)

std::vector< std::string > const &	sharedResources () const

virtual	~BaseHadronizer () noexcept(false)

Private Member Functions
void	doSetRandomEngine (CLHEP::HepRandomEngine *v) override

std::vector< std::string > const &	doSharedResources () const override

void	fillTmpStorage ()

void	flushTmpStorage ()

void	imposeProperTime ()

Private Attributes
double	fComEnergy

bool	fConvertToPDG

bool	fDisplayPythiaBanner

bool	fDisplayPythiaCards

double	fextCrossSection

double	fextCrossSectionError

double	fFilterEfficiency

bool	fHepMCVerbosity

unsigned int	fMaxEventsToPrint

edm::ParameterSet	fParameters

Pythia6Service *	fPy6Service

unsigned int	fPythiaListVerbosity

Static Private Attributes
static const std::vector< std::string >	theSharedResources

Additional Inherited Members
Protected Member Functions inherited from gen::BaseHadronizer
std::unique_ptr< HepMC::GenEvent > &	event ()

std::unique_ptr< GenEventInfoProduct > &	eventInfo ()

lhef::LHEEvent *	lheEvent ()

lhef::LHERunInfo *	lheRunInfo ()

GenRunInfoProduct &	runInfo ()

Protected Attributes inherited from gen::BaseHadronizer
std::string	lheFile_

int	randomIndex_

Detailed Description

Definition at line 24 of file Cascade2Hadronizer.h.

Constructor & Destructor Documentation

◆ Cascade2Hadronizer()

gen::Cascade2Hadronizer::Cascade2Hadronizer ( edm::ParameterSet const & ps )

Definition at line 80 of file Cascade2Hadronizer.cc.

       : BaseHadronizer(pset),
         fPy6Service(new Pythia6ServiceWithCallback(pset)),  //-- this will store py6 parameters for further settings
  
         fComEnergy(pset.getParameter<double>("comEnergy")),
         fextCrossSection(pset.getUntrackedParameter<double>("crossSection", -1.)),
         fextCrossSectionError(pset.getUntrackedParameter<double>("crossSectionError", -1.)),
         fFilterEfficiency(pset.getUntrackedParameter<double>("filterEfficiency", -1.)),
  
         fMaxEventsToPrint(pset.getUntrackedParameter<int>("maxEventsToPrint", 0)),
         fHepMCVerbosity(pset.getUntrackedParameter<bool>("pythiaHepMCVerbosity", false)),
         fPythiaListVerbosity(pset.getUntrackedParameter<int>("pythiaPylistVerbosity", 0)),
  
         fDisplayPythiaBanner(pset.getUntrackedParameter<bool>("displayPythiaBanner", false)),
         fDisplayPythiaCards(pset.getUntrackedParameter<bool>("displayPythiaCards", false)) {
     fParameters = pset.getParameter<ParameterSet>("Cascade2Parameters");
  
     fConvertToPDG = false;
     if (pset.exists("doPDGConvert"))
       fConvertToPDG = pset.getParameter<bool>("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();
   }

References gen::call_pygive(), edm::errors::Configuration, Exception, fConvertToPDG, fDisplayPythiaBanner, fDisplayPythiaCards, flushTmpStorage(), fParameters, and muonDTDigis_cfi::pset.

◆ ~Cascade2Hadronizer()

gen::Cascade2Hadronizer::~Cascade2Hadronizer ( )

override

Definition at line 121 of file Cascade2Hadronizer.cc.

                                           {
     if (fPy6Service != nullptr)
       delete fPy6Service;
   }

References fPy6Service.

Member Function Documentation

◆ cascadePrintParameters()

void gen::Cascade2Hadronizer::cascadePrintParameters ( )

Definition at line 724 of file Cascade2Hadronizer.cc.

                                                   {
     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: "<<jpsi.i23s<<endl;
     cout << "use matrix element including J/psi (Upsilon) polarisation: " << jpsi.ipsipol << endl;
     cout << "pt2 cut in matrix element for massless partons: " << captcut.pt2cut[capar1.ipro - 1] << " GeV^2" << endl;
     cout << "" << endl;
  
     cout << "parton shower and fragmentation: " << endl;
     cout << "switch for fragmentation: " << cainpu.nfrag << endl;
     cout << "switch for parton shower (0 off - 1 initial state - 2 final state - 3 initial & final state): "
          << cainpu.ifps << endl;
     cout << "switch for time like parton shower in initial state: " << casshwr.itimshr << endl;
     cout << "select CCFM (1) or DGLAP (0) evolution: " << casshwr.iccfm << endl;
     cout << "scale switch for final state parton shower: " << cainpu.ifinal << endl;
     cout << "scale factor for final state parton shower: " << scalf.scalfaf << endl;
     cout << "switch for proton remnant treatment: " << casprre.irspl << endl;
     cout << "keep track of intermediate state in parton shower: " << cashower.ipst << endl;
     cout << "mode of interaction for e p: " << cainpu.inter << endl;
     cout << "" << endl;
  
     cout << "pdfs, couplings and scales: " << endl;
     cout << "switch for running alphaem: " << capar1.irunaem << endl;
     cout << "switch for running alphas: " << capar1.iruna << endl;
     cout << "scale in alphas: " << capar1.iq2 << endl;
     cout << "scale factor for scale in alphas: " << scalf.scalfa << endl;
     cout << "unintegrated pdf: " << cagluon.iglu << endl;
     cout << "path where updf are stored: " << caspdf.pdfpath << endl;
     cout << "" << endl;
  
     cout << "center of mass energy, cross section, filter efficiency: " << endl;
     cout << "center of mass energy: " << fComEnergy << " GeV" << endl;
     cout << "external LO cross section: " << fextCrossSection << " +/- " << fextCrossSectionError << " pb" << endl;
     cout << "filter efficiency: " << fFilterEfficiency << endl;
     cout << "" << endl;
   }

References cagluon, cahflav, cainpu, caluco, capar1, capar6, captcut, cascol, cashower, caspdf, casprre, casshwr, gather_cfg::cout, fComEnergy, fextCrossSection, fextCrossSectionError, fFilterEfficiency, integr, jpsi, and scalf.

Referenced by initializeForInternalPartons().

◆ cascadeReadParameters()

bool gen::Cascade2Hadronizer::cascadeReadParameters ( const std::string & ParameterString )

Definition at line 624 of file Cascade2Hadronizer.cc.

                                                                               {
     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;
   }

References cms::dd::accepted(), cagluon, cahflav, cainpu, caluco, capar1, capar6, captcut, cascol, cashower, casprre, casshwr, integr, jpsi, and scalf.

Referenced by initializeForInternalPartons().

◆ classname()

const char * gen::Cascade2Hadronizer::classname ( ) const

Definition at line 620 of file Cascade2Hadronizer.cc.

620 { return "gen::Cascade2Hadronizer"; }

◆ decay()

bool gen::Cascade2Hadronizer::decay ( )

Definition at line 275 of file Cascade2Hadronizer.cc.

275 { return true; }

◆ declareSpecialSettings()

bool gen::Cascade2Hadronizer::declareSpecialSettings ( const std::vector< std::string > & )

inline

Definition at line 38 of file Cascade2Hadronizer.h.

38 { return true; }

◆ declareStableParticles()

bool gen::Cascade2Hadronizer::declareStableParticles ( const std::vector< int > & _pdg )

Definition at line 515 of file Cascade2Hadronizer.cc.

                                                                          {
     vector<int> 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;
   }

References gen::call_pygive(), mps_fire::i, and gen::pycomp_().

◆ doSetRandomEngine()

void gen::Cascade2Hadronizer::doSetRandomEngine ( CLHEP::HepRandomEngine * v )

overrideprivatevirtual

Reimplemented from gen::BaseHadronizer.

Definition at line 126 of file Cascade2Hadronizer.cc.

                                                                     {
     cascade2RandomEngine = v;
     fPy6Service->setRandomEngine(v);
   }

References cascade2RandomEngine, fPy6Service, gen::Pythia6Service::setRandomEngine(), and gen::v.

◆ doSharedResources()

std::vector<std::string> const& gen::Cascade2Hadronizer::doSharedResources ( ) const

inlineoverrideprivatevirtual

Reimplemented from gen::BaseHadronizer.

Definition at line 51 of file Cascade2Hadronizer.h.

51 { return theSharedResources; }

References theSharedResources.

◆ fillTmpStorage()

void gen::Cascade2Hadronizer::fillTmpStorage ( )

private

Definition at line 145 of file Cascade2Hadronizer.cc.

                                           {
     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;
   }

References mps_fire::i, and gen::pyjets_local.

Referenced by generatePartonsAndHadronize().

◆ finalizeEvent()

void gen::Cascade2Hadronizer::finalizeEvent ( )

Definition at line 159 of file Cascade2Hadronizer.cc.

                                          {
     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().reset(new GenEventInfoProduct(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<double>(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;
   }

References gen::call_pylist(), gather_cfg::cout, debug, gen::BaseHadronizer::event(), gen::BaseHadronizer::eventInfo(), fConvertToPDG, fHepMCVerbosity, fMaxEventsToPrint, fPythiaListVerbosity, get, imposeProperTime(), pydat1, pyint1, and pypars.

◆ flushTmpStorage()

void gen::Cascade2Hadronizer::flushTmpStorage ( )

private

Definition at line 131 of file Cascade2Hadronizer.cc.

                                            {
     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;
   }

References mps_fire::i, and gen::pyjets_local.

Referenced by Cascade2Hadronizer(), and generatePartonsAndHadronize().

◆ generatePartonsAndHadronize()

bool gen::Cascade2Hadronizer::generatePartonsAndHadronize ( )

Definition at line 247 of file Cascade2Hadronizer.cc.

                                                        {
     //-- 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;
   }

References call_event(), gen::BaseHadronizer::event(), fillTmpStorage(), flushTmpStorage(), fPy6Service, gen::FortranCallback::getInstance(), hepevtio, pyint1, and gen::FortranCallback::resetIterationsPerEvent().

◆ hadronize()

bool gen::Cascade2Hadronizer::hadronize ( )

Definition at line 245 of file Cascade2Hadronizer.cc.

245 { return false; }

◆ imposeProperTime()

void gen::Cascade2Hadronizer::imposeProperTime ( )

private

Definition at line 533 of file Cascade2Hadronizer.cc.

                                             {
     // 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;
   }

References funct::abs(), gen::FortranInstance::call(), class-composition::children, gen::BaseHadronizer::event(), DQMScaleToClient_cfi::factor, fPy6Service, dqm-mbProfile::log, gen::pycomp_(), pydat1, gen::pyr_(), mathSSE::sqrt(), OrderedSet::t, vbegin, and vend.

Referenced by finalizeEvent().

◆ initializeForExternalPartons()

bool gen::Cascade2Hadronizer::initializeForExternalPartons ( )

Definition at line 443 of file Cascade2Hadronizer.cc.

443 { return false; }

◆ initializeForInternalPartons()

bool gen::Cascade2Hadronizer::initializeForInternalPartons ( )

Definition at line 445 of file Cascade2Hadronizer.cc.

                                                         {
     //-- 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<string> AllSets = fParameters.getParameter<vector<string> >("parameterSets");
  
     //-- loop over the different sets
     for (unsigned i = 0; i < AllSets.size(); ++i) {
       string Set = AllSets[i];
       vector<string> Para_Set = fParameters.getParameter<vector<string> >(Set);
  
       //-- loop over all the parameters and stop in case of mistake
       for (vector<string>::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;
   }

References cainpu, call_caend(), call_cascade(), call_cascha(), call_casini(), cascadePrintParameters(), cascadeReadParameters(), gen::Pythia6Service::closeSLHA(), edm::errors::Configuration, Exception, fComEnergy, fParameters, fPy6Service, edm::ParameterSet::getParameter(), mps_fire::i, pythia6PrintParameters(), gen::Pythia6Service::setCSAParams(), gen::Pythia6Service::setGeneralParams(), gen::Pythia6Service::setPYUPDAParams(), and gen::Pythia6Service::setSLHAParams().

◆ pythia6PrintParameters()

void gen::Cascade2Hadronizer::pythia6PrintParameters ( )

Definition at line 787 of file Cascade2Hadronizer.cc.

                                                   {
     cout << "" << endl;
     cout << "----------------------------------" << endl;
     cout << "----    Pythia6 Parameters    ----" << endl;
     cout << "----------------------------------" << endl;
     cout << "" << endl;
  
     cout << "charm mass: " << pydat2.pmas[0][3] << " GeV (default value =  1.5 GeV)" << endl;
     cout << "bottom mass: " << pydat2.pmas[0][4] << " GeV (default value = 4.8 GeV)" << endl;
     cout << "Higgs mass: " << pydat2.pmas[0][24] << " GeV (default value = 115 GeV)" << endl;
  
     cout << "lambda QCD: " << pydat1.paru[111] << " GeV (default value = 0.25 GeV)" << endl;
  
     cout << "alphas behaviour: " << pydat1.mstu[110] << " (default value = 1)" << endl;
     cout << "nr of flavours wrt lambda QCD: " << pydat1.mstu[111] << " (default value = 5)" << endl;
     cout << "min nr of flavours for alphas: " << pydat1.mstu[112] << " (default value = 3)" << endl;
     cout << "max nr of flavours for alphas: " << pydat1.mstu[113] << " (default value = 5)" << endl;
  
     cout << "maximum angle settings: " << pydat1.mstj[47] << " (default value = 0)" << endl;
  
     cout << "" << endl;
   }

References gather_cfg::cout, and pydat1.

Referenced by initializeForInternalPartons().

◆ readSettings()

bool gen::Cascade2Hadronizer::readSettings ( int key )

Definition at line 430 of file Cascade2Hadronizer.cc.

                                                {
     //-- grab Py6 instance
     Pythia6Service::InstanceWrapper guard(fPy6Service);
  
     fPy6Service->setGeneralParams();
  
     if (key == 0)
       fPy6Service->setCSAParams();
     fPy6Service->setSLHAParams();
  
     return true;
   }

References fPy6Service, crabWrapper::key, gen::Pythia6Service::setCSAParams(), gen::Pythia6Service::setGeneralParams(), and gen::Pythia6Service::setSLHAParams().

◆ residualDecay()

bool gen::Cascade2Hadronizer::residualDecay ( )

Definition at line 277 of file Cascade2Hadronizer.cc.

                                          {
     //-- 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;
   }

References gen::BaseHadronizer::event(), fPy6Service, GenParticle::GenParticle, mps_fire::i, class-composition::parent, EgammaValidation_cff::pdgid, gen::pycomp_(), gen::pydecy_(), gen::pyjets_local, and mps_update::status.

◆ statistics()

void gen::Cascade2Hadronizer::statistics ( )

Definition at line 600 of file Cascade2Hadronizer.cc.

                                       {
     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;
   }

References caeffic, call_caend(), fextCrossSection, fextCrossSectionError, gen::BaseHadronizer::runInfo(), GenRunInfoProduct::setExternalXSecLO(), GenRunInfoProduct::setExternalXSecNLO(), and GenRunInfoProduct::setInternalXSec().

Member Data Documentation

◆ fComEnergy

double gen::Cascade2Hadronizer::fComEnergy

private

Definition at line 68 of file Cascade2Hadronizer.h.

Referenced by cascadePrintParameters(), and initializeForInternalPartons().

◆ fConvertToPDG

bool gen::Cascade2Hadronizer::fConvertToPDG

private

Definition at line 80 of file Cascade2Hadronizer.h.

Referenced by Cascade2Hadronizer(), and finalizeEvent().

◆ fDisplayPythiaBanner

bool gen::Cascade2Hadronizer::fDisplayPythiaBanner

private

Definition at line 77 of file Cascade2Hadronizer.h.

Referenced by Cascade2Hadronizer().

◆ fDisplayPythiaCards

bool gen::Cascade2Hadronizer::fDisplayPythiaCards

private

Definition at line 78 of file Cascade2Hadronizer.h.

Referenced by Cascade2Hadronizer().

◆ fextCrossSection

double gen::Cascade2Hadronizer::fextCrossSection

private

Definition at line 69 of file Cascade2Hadronizer.h.

Referenced by cascadePrintParameters(), and statistics().

◆ fextCrossSectionError

double gen::Cascade2Hadronizer::fextCrossSectionError

private

Definition at line 70 of file Cascade2Hadronizer.h.

Referenced by cascadePrintParameters(), and statistics().

◆ fFilterEfficiency

double gen::Cascade2Hadronizer::fFilterEfficiency

private

Definition at line 71 of file Cascade2Hadronizer.h.

Referenced by cascadePrintParameters().

◆ fHepMCVerbosity

bool gen::Cascade2Hadronizer::fHepMCVerbosity

private

Definition at line 74 of file Cascade2Hadronizer.h.

Referenced by finalizeEvent().

◆ fMaxEventsToPrint

unsigned int gen::Cascade2Hadronizer::fMaxEventsToPrint

private

Definition at line 73 of file Cascade2Hadronizer.h.

Referenced by finalizeEvent().

◆ fParameters

edm::ParameterSet gen::Cascade2Hadronizer::fParameters

private

Definition at line 63 of file Cascade2Hadronizer.h.

Referenced by Cascade2Hadronizer(), and initializeForInternalPartons().

◆ fPy6Service

Pythia6Service* gen::Cascade2Hadronizer::fPy6Service

private

Definition at line 65 of file Cascade2Hadronizer.h.

Referenced by doSetRandomEngine(), generatePartonsAndHadronize(), imposeProperTime(), initializeForInternalPartons(), readSettings(), residualDecay(), and ~Cascade2Hadronizer().

◆ fPythiaListVerbosity

unsigned int gen::Cascade2Hadronizer::fPythiaListVerbosity

private

Definition at line 75 of file Cascade2Hadronizer.h.

Referenced by finalizeEvent().

◆ theSharedResources

const std::vector< std::string > gen::Cascade2Hadronizer::theSharedResources

staticprivate

Initial value:

= {edm::SharedResourceNames::kPythia6,

gen::FortranInstance::kFortranInstance}

Definition at line 53 of file Cascade2Hadronizer.h.

Referenced by doSharedResources().

Public Member Functions

Private Member Functions

Private Attributes

Static Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ Cascade2Hadronizer()

◆ ~Cascade2Hadronizer()

Member Function Documentation

◆ cascadePrintParameters()

◆ cascadeReadParameters()

◆ classname()

◆ decay()

◆ declareSpecialSettings()

◆ declareStableParticles()

◆ doSetRandomEngine()

◆ doSharedResources()

◆ fillTmpStorage()

◆ finalizeEvent()

◆ flushTmpStorage()

◆ generatePartonsAndHadronize()

◆ hadronize()

◆ imposeProperTime()

◆ initializeForExternalPartons()

◆ initializeForInternalPartons()

◆ pythia6PrintParameters()

◆ readSettings()

◆ residualDecay()

◆ statistics()

Member Data Documentation

◆ fComEnergy

◆ fConvertToPDG

◆ fDisplayPythiaBanner

◆ fDisplayPythiaCards

◆ fextCrossSection

◆ fextCrossSectionError

◆ fFilterEfficiency

◆ fHepMCVerbosity

◆ fMaxEventsToPrint

◆ fParameters

◆ fPy6Service

◆ fPythiaListVerbosity

◆ theSharedResources