gen::Pythia6Hadronizer Class Reference

#include <Pythia6Hadronizer.h>

Inheritance diagram for gen::Pythia6Hadronizer:

Public Member Functions
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 ()
	Pythia6Hadronizer (edm::ParameterSet const &ps)
bool	readSettings (int)
bool	residualDecay ()
void	statistics ()
	~Pythia6Hadronizer () 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)

Static Public Member Functions
static JetMatching *	getJetMatching ()

Private Types
enum	{   PP, PPbar, ElectronPositron, ElectronProton,   PositronProton }

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	fBeam1PZ
double	fBeam2PZ
double	fCOMEnergy
bool	fConvertToPDG
bool	fDisplayPythiaBanner
bool	fDisplayPythiaCards
bool	fGluinoHadronsEnabled
bool	fHepMCVerbosity
bool	fImposeProperTime
int	fInitialState
unsigned int	fMaxEventsToPrint
Pythia6Service *	fPy6Service
unsigned int	fPythiaListVerbosity
bool	fStopHadronsEnabled

Static Private Attributes
static JetMatching *	fJetMatching = nullptr
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 37 of file Pythia6Hadronizer.h.

Member Enumeration Documentation

anonymous enum

anonymous enum

private

Enumerator
PP
PPbar
ElectronPositron
ElectronProton
PositronProton

Definition at line 82 of file Pythia6Hadronizer.h.

{ PP, PPbar, ElectronPositron, ElectronProton, PositronProton };

Constructor & Destructor Documentation

Pythia6Hadronizer()

gen::Pythia6Hadronizer::Pythia6Hadronizer

(

edm::ParameterSet const &

ps

)

Definition at line 97 of file Pythia6Hadronizer.cc.

References gen::call_pygive(), edm::errors::Configuration, gather_cfg::cout, gen::JetMatching::create(), ElectronPositron, ElectronProton, Exception, edm::ParameterSet::exists(), fBeam1PZ, fBeam2PZ, fConvertToPDG, fDisplayPythiaBanner, fDisplayPythiaCards, fGluinoHadronsEnabled, fImposeProperTime, fInitialState, fJetMatching, flushTmpStorage(), fStopHadronsEnabled, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), PositronProton, PP, and PPbar.

      : BaseHadronizer(ps),
        fPy6Service(new Pythia6ServiceWithCallback(ps)),  // this will store py6 params for further settings
        fInitialState(PP),
        fCOMEnergy(ps.getParameter<double>("comEnergy")),
        fHepMCVerbosity(ps.getUntrackedParameter<bool>("pythiaHepMCVerbosity", false)),
        fMaxEventsToPrint(ps.getUntrackedParameter<int>("maxEventsToPrint", 0)),
        fPythiaListVerbosity(ps.getUntrackedParameter<int>("pythiaPylistVerbosity", 0)),
        fDisplayPythiaBanner(ps.getUntrackedParameter<bool>("displayPythiaBanner", false)),
        fDisplayPythiaCards(ps.getUntrackedParameter<bool>("displayPythiaCards", false)) {
    // J.Y.: the following 3 parameters are hacked "for a reason"
    //
    if (ps.exists("PPbarInitialState")) {
      if (fInitialState == PP) {
        fInitialState = PPbar;
        edm::LogInfo("GeneratorInterface|Pythia6Interface")
            << "Pythia6 will be initialized for PROTON-ANTIPROTON INITIAL STATE. "
            << "This is a user-request change from the DEFAULT PROTON-PROTON initial state." << std::endl;
        std::cout << "Pythia6 will be initialized for PROTON-ANTIPROTON INITIAL STATE." << std::endl;
        std::cout << "This is a user-request change from the DEFAULT PROTON-PROTON initial state." << std::endl;
      } else {
        // probably need to throw on attempt to override ?
      }
    } else if (ps.exists("ElectronPositronInitialState")) {
      if (fInitialState == PP) {
        fInitialState = ElectronPositron;
        edm::LogInfo("GeneratorInterface|Pythia6Interface")
            << "Pythia6 will be initialized for ELECTRON-POSITRON INITIAL STATE. "
            << "This is a user-request change from the DEFAULT PROTON-PROTON initial state." << std::endl;
        std::cout << "Pythia6 will be initialized for ELECTRON-POSITRON INITIAL STATE." << std::endl;
        std::cout << "This is a user-request change from the DEFAULT PROTON-PROTON initial state." << std::endl;
      } else {
        // probably need to throw on attempt to override ?
      }
    } else if (ps.exists("ElectronProtonInitialState")) {
      if (fInitialState == PP) {
        fInitialState = ElectronProton;
        fBeam1PZ =
            (ps.getParameter<edm::ParameterSet>("ElectronProtonInitialState")).getParameter<double>("electronMomentum");
        fBeam2PZ =
            (ps.getParameter<edm::ParameterSet>("ElectronProtonInitialState")).getParameter<double>("protonMomentum");
        edm::LogInfo("GeneratorInterface|Pythia6Interface")
            << "Pythia6 will be initialized for ELECTRON-PROTON INITIAL STATE. "
            << "This is a user-request change from the DEFAULT PROTON-PROTON initial state." << std::endl;
        std::cout << "Pythia6 will be initialized for ELECTRON-PROTON INITIAL STATE." << std::endl;
        std::cout << "This is a user-request change from the DEFAULT PROTON-PROTON initial state." << std::endl;
      } else {
        // probably need to throw on attempt to override ?
      }
    } else if (ps.exists("PositronProtonInitialState")) {
      if (fInitialState == PP) {
        fInitialState = PositronProton;
        fBeam1PZ =
            (ps.getParameter<edm::ParameterSet>("ElectronProtonInitialState")).getParameter<double>("positronMomentum");
        fBeam2PZ =
            (ps.getParameter<edm::ParameterSet>("ElectronProtonInitialState")).getParameter<double>("protonMomentum");
        edm::LogInfo("GeneratorInterface|Pythia6Interface")
            << "Pythia6 will be initialized for POSITRON-PROTON INITIAL STATE. "
            << "This is a user-request change from the DEFAULT PROTON-PROTON initial state." << std::endl;
        std::cout << "Pythia6 will be initialized for POSITRON-PROTON INITIAL STATE." << std::endl;
        std::cout << "This is a user-request change from the DEFAULT PROTON-PROTON initial state." << std::endl;
      } else {
        // throw on unknown initial state !
        throw edm::Exception(edm::errors::Configuration, "Pythia6Interface")
            << " UNKNOWN INITIAL STATE. \n The allowed initial states are: PP, PPbar, ElectronPositron, "
               "ElectronProton, and PositronProton \n";
      }
    }

    // J.Y.: the following 4 params are "hacked", in the sense
    // that they're tracked but get in optionally;
    // this will be fixed once we update all applications
    //

    fStopHadronsEnabled = false;
    if (ps.exists("stopHadrons"))
      fStopHadronsEnabled = ps.getParameter<bool>("stopHadrons");

    fGluinoHadronsEnabled = false;
    if (ps.exists("gluinoHadrons"))
      fGluinoHadronsEnabled = ps.getParameter<bool>("gluinoHadrons");

    fImposeProperTime = false;
    if (ps.exists("imposeProperTime")) {
      fImposeProperTime = ps.getParameter<bool>("imposeProperTime");
    }

    fConvertToPDG = false;
    if (ps.exists("doPDGConvert"))
      fConvertToPDG = ps.getParameter<bool>("doPDGConvert");

    if (ps.exists("jetMatching")) {
      edm::ParameterSet jmParams = ps.getUntrackedParameter<edm::ParameterSet>("jetMatching");

      fJetMatching = JetMatching::create(jmParams).release();
    }

    // first of all, 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
    // NPartsBeforeDecays = 0;
    flushTmpStorage();
  }

~Pythia6Hadronizer()

gen::Pythia6Hadronizer::~Pythia6Hadronizer ( )

override

Definition at line 215 of file Pythia6Hadronizer.cc.

References fJetMatching, and fPy6Service.

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

Member Function Documentation

classname()

const char * gen::Pythia6Hadronizer::classname

(

)

const

Definition at line 910 of file Pythia6Hadronizer.cc.

{ return "gen::Pythia6Hadronizer"; }

decay()

bool gen::Pythia6Hadronizer::decay

(

)

Definition at line 509 of file Pythia6Hadronizer.cc.

{ return true; }

declareSpecialSettings()

bool gen::Pythia6Hadronizer::declareSpecialSettings

(

const std::vector< std::string > &

settings

)

Definition at line 796 of file Pythia6Hadronizer.cc.

References gen::call_pygive(), edm::errors::Configuration, Exception, spr::find(), contentValuesFiles::number, gen::pycomp_(), gen::pydat1_, and AlCaHLTBitMon_QueryRunRegistry::string.

                                                                                     {
    for (unsigned int iss = 0; iss < settings.size(); iss++) {
      if (settings[iss].find("QED-brem-off") == std::string::npos)
        continue;
      size_t fnd1 = settings[iss].find(':');
      if (fnd1 == std::string::npos)
        continue;

      std::string value = settings[iss].substr(fnd1 + 1);

      if (value == "all") {
        call_pygive("MSTJ(41)=3");
      } else {
        int number = atoi(value.c_str());
        int PyID = HepPID::translatePDTtoPythia(number);
        int pyCode = pycomp_(PyID);
        if (pyCode > 0) {
          // first of all, check if mstj(39) is 0 or if we're trying to override user's setting
          // if so, throw an exception and stop, because otherwise the user will get behaviour
          // that's different from what she/he expects !
          if (pydat1_.mstj[38] > 0 && pydat1_.mstj[38] != pyCode) {
            throw edm::Exception(edm::errors::Configuration, "Pythia6Interface")
                << " Fatal conflict: \n mandatory internal directive to set MSTJ(39)=" << pyCode
                << " overrides user setting MSTJ(39)=" << pydat1_.mstj[38]
                << " - user will not get expected behaviour \n";
          }
          std::ostringstream pyCard;
          pyCard << "MSTJ(39)=" << pyCode;
          call_pygive(pyCard.str());
        }
      }
    }

    return true;
  }

declareStableParticles()

bool gen::Pythia6Hadronizer::declareStableParticles

(

const std::vector< int > &

pdg

)

Definition at line 778 of file Pythia6Hadronizer.cc.

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

                                                                          {
    for (unsigned int i = 0; i < pdg.size(); i++) {
      int PyID = HepPID::translatePDTtoPythia(pdg[i]);
      // int PyID = pdg[i];
      int pyCode = pycomp_(PyID);
      if (pyCode > 0) {
        std::ostringstream pyCard;
        pyCard << "MDCY(" << pyCode << ",1)=0";
        /* this is a test printout... 
         std::cout << "pdg= " << pdg[i] << " " << pyCard.str() << std::endl; 
*/
        call_pygive(pyCard.str());
      }
    }

    return true;
  }

doSetRandomEngine()

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

overrideprivatevirtual

Reimplemented from gen::BaseHadronizer.

Definition at line 222 of file Pythia6Hadronizer.cc.

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

{ fPy6Service->setRandomEngine(v); }

doSharedResources()

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

inlineoverrideprivatevirtual

Reimplemented from gen::BaseHadronizer.

Definition at line 66 of file Pythia6Hadronizer.h.

References theSharedResources.

{ return theSharedResources; }

fillTmpStorage()

void gen::Pythia6Hadronizer::fillTmpStorage ( )

private

Definition at line 237 of file Pythia6Hadronizer.cc.

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

Referenced by generatePartonsAndHadronize(), and hadronize().

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

finalizeEvent()

void gen::Pythia6Hadronizer::finalizeEvent

(

)

Definition at line 251 of file Pythia6Hadronizer.cc.

References funct::abs(), gen::call_pylist(), gather_cfg::cout, callgraph::cs, gen::BaseHadronizer::event(), gen::BaseHadronizer::eventInfo(), fConvertToPDG, fDisplayPythiaCards, fHepMCVerbosity, lhef::LHEEvent::fillEventInfo(), lhef::LHEEvent::fillPdfInfo(), fImposeProperTime, alignBH_cfg::fixed, fMaxEventsToPrint, fPythiaListVerbosity, imposeProperTime(), gen::BaseHadronizer::lheEvent(), gen::BaseHadronizer::lheRunInfo(), pydat1, pyint1, pypars, and mathSSE::sqrt().

                                        {
    bool lhe = lheEvent() != nullptr;

    HepMC::PdfInfo pdf;

    // if we are in hadronizer mode, we can pass on information from
    // the LHE input
    if (lhe) {
      lheEvent()->fillEventInfo(event().get());
      lheEvent()->fillPdfInfo(&pdf);
    } else {
      // filling in 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]);
    }

    /* 9/9/2010 - JVY: This is the old piece of code - I can't remember why we implemented it this way.
                   However, it's causing problems with pdf1 & pdf2 when processing LHE samples,
           specifically, because both are set to -1, it tries to fill with Py6 numbers that
           are NOT valid/right at this point !
           In general, for LHE/ME event processing we should implement the correct calculation
           of the pdf's, rather than using py6 ones.

   // filling in factorization "Q scale" now! pthat moved to binningValues()

   if (!lhe || event()->signal_process_id() < 0) event()->set_signal_process_id( pypars.msti[0] );
   if (!lhe || event()->event_scale() < 0)       event()->set_event_scale( pypars.pari[22] );
   if (!lhe || event()->alphaQED() < 0)          event()->set_alphaQED( pyint1.vint[56] );
   if (!lhe || 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 (!lhe || pdf.id1() < 0)      pdf.set_id1( pyint1.mint[14] == 21 ? 0 : pyint1.mint[14] );
   if (!lhe || pdf.id2() < 0)      pdf.set_id2( pyint1.mint[15] == 21 ? 0 : pyint1.mint[15] );
   if (!lhe || pdf.x1() < 0)       pdf.set_x1( pyint1.vint[40] );
   if (!lhe || pdf.x2() < 0)       pdf.set_x2( pyint1.vint[41] );
   if (!lhe || pdf.pdf1() < 0)     pdf.set_pdf1( pyint1.vint[38] / pyint1.vint[40] );
   if (!lhe || pdf.pdf2() < 0)     pdf.set_pdf2( pyint1.vint[39] / pyint1.vint[41] );
   if (!lhe || pdf.scalePDF() < 0) pdf.set_scalePDF( pyint1.vint[50] );
*/

    event()->set_pdf_info(pdf);

    HepMC::GenCrossSection xsec;
    double cs = pypars.pari[0];  // cross section in mb
    cs *= 1.0e9;                 // translate to pb
    double cserr = cs / sqrt(pypars.msti[4]);
    xsec.set_cross_section(cs, cserr);
    event()->set_cross_section(xsec);

    // this is "standard" Py6 event weight (corresponds to PYINT1/VINT(97)
    //
    if (lhe && std::abs(lheRunInfo()->getHEPRUP()->IDWTUP) == 4)
      // translate mb to pb (CMS/Gen "convention" as of May 2009)
      event()->weights().push_back(pyint1.vint[96] * 1.0e9);
    else
      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]));

    // now create the GenEventInfo product from the GenEvent and fill
    // the missing pieces

    eventInfo() = std::make_unique<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
    if (!lhe) {
      eventInfo()->setBinningValues(std::vector<double>(1, pypars.pari[16]));
    }

    // here we treat long-lived particles
    //
    if (fImposeProperTime || 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) {
        std::cout << "Event process = " << pypars.msti[0] << std::endl << "----------------------" << std::endl;
        event()->print();
      }
    }

    // dump of all settings after all initializations
    if (fDisplayPythiaCards) {
      fDisplayPythiaCards = false;
      call_pylist(12);
      call_pylist(13);
      std::cout << "\n PYPARS \n" << std::endl;
      std::cout << std::setw(5) << std::fixed << "I" << std::setw(10) << std::fixed << "MSTP(I)" << std::setw(16)
                << std::fixed << "PARP(I)" << std::setw(10) << std::fixed << "MSTI(I)" << std::setw(16) << std::fixed
                << "PARI(I)" << std::endl;
      for (unsigned int ind = 0; ind < 200; ind++) {
        std::cout << std::setw(5) << std::fixed << ind + 1 << std::setw(10) << std::fixed << pypars.mstp[ind]
                  << std::setw(16) << std::fixed << pypars.parp[ind] << std::setw(10) << std::fixed << pypars.msti[ind]
                  << std::setw(16) << std::fixed << pypars.pari[ind] << std::endl;
      }
    }

    return;
  }

flushTmpStorage()

void gen::Pythia6Hadronizer::flushTmpStorage ( )

private

Definition at line 224 of file Pythia6Hadronizer.cc.

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

Referenced by generatePartonsAndHadronize(), hadronize(), and Pythia6Hadronizer().

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

generatePartonsAndHadronize()

bool gen::Pythia6Hadronizer::generatePartonsAndHadronize

(

)

Definition at line 396 of file Pythia6Hadronizer.cc.

References gen::call_pygive(), gen::BaseHadronizer::event(), fGluinoHadronsEnabled, fillTmpStorage(), flushTmpStorage(), fPy6Service, fStopHadronsEnabled, gen::FortranCallback::getInstance(), gen::pyglfr_(), pyint1, gen::pystfr_(), pythia6_conv, and gen::FortranCallback::resetIterationsPerEvent().

                                                      {
    Pythia6Service::InstanceWrapper guard(fPy6Service);  // grab Py6 instance

    FortranCallback::getInstance()->resetIterationsPerEvent();

    // generate event with Pythia6
    //

    if (fStopHadronsEnabled || fGluinoHadronsEnabled) {
      // call_pygive("MSTJ(1)=-1");
      call_pygive("MSTJ(14)=-1");
    }

    call_pyevnt();

    if (fStopHadronsEnabled || fGluinoHadronsEnabled) {
      // call_pygive("MSTJ(1)=1");
      call_pygive("MSTJ(14)=1");
      int ierr = 0;
      if (fStopHadronsEnabled) {
        pystfr_(ierr);
        if (ierr != 0)  // skip failed events
        {
          event().reset();
          return false;
        }
      }
      if (fGluinoHadronsEnabled)
        pyglfr_();
    }

    if (pyint1.mint[50] != 0)  // skip event if py6 considers it bad
    {
      event().reset();
      return false;
    }

    call_pyhepc(1);
    event().reset(pythia6_conv.read_next_event());

    // this is to deal with post-gen tools & residualDecay() that may reuse PYJETS
    //
    flushTmpStorage();
    fillTmpStorage();

    return true;
  }

getJetMatching()

static JetMatching* gen::Pythia6Hadronizer::getJetMatching ( )

inlinestatic

Definition at line 53 of file Pythia6Hadronizer.h.

References fJetMatching.

Referenced by gen::Pythia6ServiceWithCallback::upEvnt(), and gen::Pythia6ServiceWithCallback::upVeto().

{ return fJetMatching; }

hadronize()

bool gen::Pythia6Hadronizer::hadronize

(

)

Definition at line 444 of file Pythia6Hadronizer.cc.

References gen::JetMatching::beforeHadronisation(), gen::call_pygive(), lhef::LHEEvent::count(), gen::BaseHadronizer::event(), fGluinoHadronsEnabled, fillTmpStorage(), fJetMatching, flushTmpStorage(), fPy6Service, fStopHadronsEnabled, gen::FortranCallback::getInstance(), hepeup_, lhef::LHERunInfo::kAccepted, lhef::LHERunInfo::kSelected, gen::BaseHadronizer::lheEvent(), HEPEUP_::nup, gen::pyglfr_(), pyint1, pypars, gen::pystfr_(), pythia6_conv, gen::FortranCallback::resetIterationsPerEvent(), gen::JetMatching::resetMatchingStatus(), and gen::FortranCallback::setLHEEvent().

                                    {
    Pythia6Service::InstanceWrapper guard(fPy6Service);  // grab Py6 instance

    FortranCallback::getInstance()->setLHEEvent(lheEvent());
    FortranCallback::getInstance()->resetIterationsPerEvent();
    if (fJetMatching) {
      fJetMatching->resetMatchingStatus();
      fJetMatching->beforeHadronisation(lheEvent());
    }

    // generate event with Pythia6
    //
    if (fStopHadronsEnabled || fGluinoHadronsEnabled) {
      call_pygive("MSTJ(1)=-1");
      call_pygive("MSTJ(14)=-1");
    }

    call_pyevnt();

    if (FortranCallback::getInstance()->getIterationsPerEvent() > 1 || hepeup_.nup <= 0 || pypars.msti[0] == 1) {
      // update LHE matching statistics
      lheEvent()->count(lhef::LHERunInfo::kSelected);

      event().reset();
      return false;
    }

    // update LHE matching statistics
    //
    lheEvent()->count(lhef::LHERunInfo::kAccepted);

    if (fStopHadronsEnabled || fGluinoHadronsEnabled) {
      call_pygive("MSTJ(1)=1");
      call_pygive("MSTJ(14)=1");
      int ierr = 0;
      if (fStopHadronsEnabled) {
        pystfr_(ierr);
        if (ierr != 0)  // skip failed events
        {
          event().reset();
          return false;
        }
      }

      if (fGluinoHadronsEnabled)
        pyglfr_();
    }

    if (pyint1.mint[50] != 0)  // skip event if py6 considers it bad
    {
      event().reset();
      return false;
    }

    call_pyhepc(1);
    event().reset(pythia6_conv.read_next_event());

    // this is to deal with post-gen tools and/or residualDecay(), that may reuse PYJETS
    //
    flushTmpStorage();
    fillTmpStorage();

    return true;
  }

imposeProperTime()

void gen::Pythia6Hadronizer::imposeProperTime ( )

private

Definition at line 832 of file Pythia6Hadronizer.cc.

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(), submitPVValidationJobs::t, and x.

Referenced by finalizeEvent().

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

initializeForExternalPartons()

bool gen::Pythia6Hadronizer::initializeForExternalPartons

(

)

Definition at line 671 of file Pythia6Hadronizer.cc.

References gen::call_pygive(), gen::Pythia6Service::closeSLHA(), fGluinoHadronsEnabled, lhef::LHERunInfo::findHeader(), fJetMatching, fPy6Service, fStopHadronsEnabled, gen::FortranCallback::getInstance(), gen::JetMatching::init(), gen::BaseHadronizer::lheRunInfo(), gen::pyglrhad_(), gen::pystrhad_(), gen::FortranCallback::setLHERunInfo(), gen::Pythia6Service::setPYUPDAParams(), and gen::Pythia6Service::setSLHAFromHeader().

                                                       {
    Pythia6Service::InstanceWrapper guard(fPy6Service);  // grab Py6 instance

    // note: CSA mode is NOT supposed to work with external partons !!!

    // fPy6Service->setGeneralParams();

    fPy6Service->setPYUPDAParams(false);

    FortranCallback::getInstance()->setLHERunInfo(lheRunInfo());

    if (fStopHadronsEnabled) {
      // overwrite mstp(111), no matter what
      call_pygive("MSTP(111)=0");
      pystrhad_();
      call_pygive("MWID(302)=0");    // I don't know if this is specific to processing ME/LHE only,
      call_pygive("MDCY(302,1)=0");  // or this should also be the case for full event...
                                     // anyway, this comes from experience of processing MG events
    }

    if (fGluinoHadronsEnabled) {
      // overwrite mstp(111), no matter what
      call_pygive("MSTP(111)=0");
      pyglrhad_();
      //call_pygive("MWID(309)=0");
      //call_pygive("MDCY(309,1)=0");
    }

    call_pyinit("USER", "", "", 0.0);

    fPy6Service->setPYUPDAParams(true);

    std::vector<std::string> slha = lheRunInfo()->findHeader("slha");
    if (!slha.empty()) {
      edm::LogInfo("Generator|LHEInterface") << "Pythia6 hadronisation found an SLHA header, "
                                             << "will be passed on to Pythia." << std::endl;
      fPy6Service->setSLHAFromHeader(slha);
      fPy6Service->closeSLHA();
    }

    if (fJetMatching) {
      fJetMatching->init(lheRunInfo());
      // FIXME: the jet matching routine might not be interested in PS callback
      call_pygive("MSTP(143)=1");
    }

    return true;
  }

initializeForInternalPartons()

bool gen::Pythia6Hadronizer::initializeForInternalPartons

(

)

Definition at line 720 of file Pythia6Hadronizer.cc.

References gen::call_pygive(), gen::Pythia6Service::closeSLHA(), edm::errors::Configuration, ElectronPositron, ElectronProton, Exception, fBeam1PZ, fBeam2PZ, fCOMEnergy, fGluinoHadronsEnabled, fInitialState, fPy6Service, fStopHadronsEnabled, PositronProton, PP, PPbar, gen::pyglrhad_(), gen::pystrhad_(), and gen::Pythia6Service::setPYUPDAParams().

                                                       {
    Pythia6Service::InstanceWrapper guard(fPy6Service);  // grab Py6 instance

    fPy6Service->setPYUPDAParams(false);

    if (fStopHadronsEnabled) {
      // overwrite mstp(111), no matter what
      call_pygive("MSTP(111)=0");
      pystrhad_();
    }

    if (fGluinoHadronsEnabled) {
      // overwrite mstp(111), no matter what
      call_pygive("MSTP(111)=0");
      pyglrhad_();
    }

    if (fInitialState == PP)  // default
    {
      call_pyinit("CMS", "p", "p", fCOMEnergy);
    } else if (fInitialState == PPbar) {
      call_pyinit("CMS", "p", "pbar", fCOMEnergy);
    } else if (fInitialState == ElectronPositron) {
      call_pyinit("CMS", "e+", "e-", fCOMEnergy);
    } else if (fInitialState == ElectronProton) {
      // set p(1,i) & p(2,i) for the beams in pyjets !
      pyjets.p[0][0] = 0.;
      pyjets.p[1][0] = 0.;
      pyjets.p[2][0] = fBeam1PZ;
      pyjets.p[0][1] = 0.;
      pyjets.p[1][1] = 0.;
      pyjets.p[2][1] = fBeam2PZ;
      // call "3mon" frame & 0.0 win
      call_pyinit("3mom", "e-", "p", 0.0);
    } else if (fInitialState == PositronProton) {
      // set p(1,i) & p(2,i) for the beams in pyjets !
      pyjets.p[0][0] = 0.;
      pyjets.p[1][0] = 0.;
      pyjets.p[2][0] = fBeam1PZ;
      pyjets.p[0][1] = 0.;
      pyjets.p[1][1] = 0.;
      pyjets.p[2][1] = fBeam2PZ;
      // call "3mon" frame & 0.0 win
      call_pyinit("3mom", "e+", "p", 0.0);
    } else {
      // throw on unknown initial state !
      throw edm::Exception(edm::errors::Configuration, "Pythia6Interface")
          << " UNKNOWN INITIAL STATE. \n The allowed initial states are: PP, PPbar, ElectronPositron, ElectronProton, "
             "and PositronProton \n";
    }

    fPy6Service->setPYUPDAParams(true);

    fPy6Service->closeSLHA();

    return true;
  }

readSettings()

bool gen::Pythia6Hadronizer::readSettings

(

int

key

)

Definition at line 660 of file Pythia6Hadronizer.cc.

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

                                              {
    Pythia6Service::InstanceWrapper guard(fPy6Service);  // grab Py6 instance

    fPy6Service->setGeneralParams();
    if (key == 0)
      fPy6Service->setCSAParams();
    fPy6Service->setSLHAParams();

    return true;
  }

residualDecay()

bool gen::Pythia6Hadronizer::residualDecay

(

)

Definition at line 511 of file Pythia6Hadronizer.cc.

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.

                                        {
    // event().get()->print();

    Pythia6Service::InstanceWrapper guard(fPy6Service);  // grab Py6 instance

    // int nDocLines = pypars.msti[3];

    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 vertexes, i.e. in the middle
    // of the event rather than at the end; but this will only shift poiters 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);
        }

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

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

statistics()

void gen::Pythia6Hadronizer::statistics

(

)

Definition at line 896 of file Pythia6Hadronizer.cc.

References callgraph::cs, pypars, gen::BaseHadronizer::runInfo(), and GenRunInfoProduct::setInternalXSec().

                                     {
    if (!runInfo().internalXSec()) {
      // set xsec if not already done (e.g. from LHE cross section collector)
      double cs = pypars.pari[0];  // cross section in mb
      cs *= 1.0e9;                 // translate to pb (CMS/Gen "convention" as of May 2009)
      runInfo().setInternalXSec(cs);
      // FIXME: can we get the xsec statistical error somewhere?
    }

    call_pystat(1);

    return;
  }

Member Data Documentation

fBeam1PZ

double gen::Pythia6Hadronizer::fBeam1PZ

private

Definition at line 87 of file Pythia6Hadronizer.h.

Referenced by initializeForInternalPartons(), and Pythia6Hadronizer().

fBeam2PZ

double gen::Pythia6Hadronizer::fBeam2PZ

private

Definition at line 88 of file Pythia6Hadronizer.h.

Referenced by initializeForInternalPartons(), and Pythia6Hadronizer().

fCOMEnergy

double gen::Pythia6Hadronizer::fCOMEnergy

private

Definition at line 85 of file Pythia6Hadronizer.h.

Referenced by initializeForInternalPartons().

fConvertToPDG

bool gen::Pythia6Hadronizer::fConvertToPDG

private

Definition at line 112 of file Pythia6Hadronizer.h.

Referenced by finalizeEvent(), and Pythia6Hadronizer().

fDisplayPythiaBanner

bool gen::Pythia6Hadronizer::fDisplayPythiaBanner

private

Definition at line 98 of file Pythia6Hadronizer.h.

Referenced by Pythia6Hadronizer().

fDisplayPythiaCards

bool gen::Pythia6Hadronizer::fDisplayPythiaCards

private

Definition at line 99 of file Pythia6Hadronizer.h.

Referenced by finalizeEvent(), and Pythia6Hadronizer().

fGluinoHadronsEnabled

bool gen::Pythia6Hadronizer::fGluinoHadronsEnabled

private

Definition at line 106 of file Pythia6Hadronizer.h.

Referenced by generatePartonsAndHadronize(), hadronize(), initializeForExternalPartons(), initializeForInternalPartons(), and Pythia6Hadronizer().

fHepMCVerbosity

bool gen::Pythia6Hadronizer::fHepMCVerbosity

private

Definition at line 92 of file Pythia6Hadronizer.h.

Referenced by finalizeEvent().

fImposeProperTime

bool gen::Pythia6Hadronizer::fImposeProperTime

private

Definition at line 109 of file Pythia6Hadronizer.h.

Referenced by finalizeEvent(), and Pythia6Hadronizer().

fInitialState

int gen::Pythia6Hadronizer::fInitialState

private

Definition at line 84 of file Pythia6Hadronizer.h.

Referenced by initializeForInternalPartons(), and Pythia6Hadronizer().

fJetMatching

JetMatching * gen::Pythia6Hadronizer::fJetMatching = nullptr

staticprivate

Definition at line 90 of file Pythia6Hadronizer.h.

Referenced by getJetMatching(), hadronize(), initializeForExternalPartons(), Pythia6Hadronizer(), and ~Pythia6Hadronizer().

fMaxEventsToPrint

unsigned int gen::Pythia6Hadronizer::fMaxEventsToPrint

private

Definition at line 93 of file Pythia6Hadronizer.h.

Referenced by finalizeEvent().

fPy6Service

Pythia6Service* gen::Pythia6Hadronizer::fPy6Service

private

Definition at line 76 of file Pythia6Hadronizer.h.

Referenced by doSetRandomEngine(), generatePartonsAndHadronize(), hadronize(), imposeProperTime(), initializeForExternalPartons(), initializeForInternalPartons(), readSettings(), residualDecay(), and ~Pythia6Hadronizer().

fPythiaListVerbosity

unsigned int gen::Pythia6Hadronizer::fPythiaListVerbosity

private

Definition at line 97 of file Pythia6Hadronizer.h.

Referenced by finalizeEvent().

fStopHadronsEnabled

bool gen::Pythia6Hadronizer::fStopHadronsEnabled

private

Definition at line 105 of file Pythia6Hadronizer.h.

Referenced by generatePartonsAndHadronize(), hadronize(), initializeForExternalPartons(), initializeForInternalPartons(), and Pythia6Hadronizer().

theSharedResources

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

staticprivate

Initial value:

= {edm::SharedResourceNames::kPythia6,
                                                                          gen::FortranInstance::kFortranInstance}

Definition at line 117 of file Pythia6Hadronizer.h.

Referenced by doSharedResources().