#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	residualDecay ()

void	statistics ()

	~Pythia6Hadronizer ()

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

edm::Event &	getEDMEvent () const

HepMC::GenEvent *	getGenEvent ()

GenEventInfoProduct *	getGenEventInfo ()

GenRunInfoProduct &	getGenRunInfo ()

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

void	resetEvent (HepMC::GenEvent *event)

void	resetEventInfo (GenEventInfoProduct *eventInfo)

void	setEDMEvent (edm::Event &event)

void	setLHEEvent (lhef::LHEEvent *event)

void	setLHERunInfo (lhef::LHERunInfo *runInfo)

	~BaseHadronizer ()

Static Public Member Functions
static JetMatching *	getJetMatching ()

Private Member Functions
void	fillTmpStorage ()

void	flushTmpStorage ()

void	imposeProperTime ()

Private Attributes
double	fCOMEnergy

bool	fConvertToPDG

bool	fDisplayPythiaBanner

bool	fDisplayPythiaCards

bool	fGluinoHadronsEnabled

bool	fHepMCVerbosity

bool	fImposeProperTime

unsigned int	fMaxEventsToPrint

Pythia6Service *	fPy6Service

unsigned int	fPythiaListVerbosity

bool	fStopHadronsEnabled

Static Private Attributes
static JetMatching *	fJetMatching = 0

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

std::auto_ptr < GenEventInfoProduct > &	eventInfo ()

lhef::LHEEvent *	lheEvent ()

lhef::LHERunInfo *	lheRunInfo ()

GenRunInfoProduct &	runInfo ()

Detailed Description

Definition at line 37 of file Pythia6Hadronizer.h.

Constructor & Destructor Documentation

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

Definition at line 92 of file Pythia6Hadronizer.cc.

References gen::call_pygive(), edm::errors::Configuration, gen::JetMatching::create(), edm::hlt::Exception, edm::ParameterSet::exists(), fConvertToPDG, fDisplayPythiaBanner, fDisplayPythiaCards, fGluinoHadronsEnabled, fImposeProperTime, fJetMatching, flushTmpStorage(), fStopHadronsEnabled, edm::ParameterSet::getParameter(), and edm::ParameterSet::getUntrackedParameter().

    : BaseHadronizer(ps),
      fPy6Service( new Pythia6ServiceWithCallback(ps) ), // this will store py6 params for further settings
      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 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();
    
 }

gen::Pythia6Hadronizer::~Pythia6Hadronizer ( )

Definition at line 163 of file Pythia6Hadronizer.cc.

References fJetMatching, and fPy6Service.

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

Member Function Documentation

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

Definition at line 865 of file Pythia6Hadronizer.cc.

 {
    return "gen::Pythia6Hadronizer";
 }

bool gen::Pythia6Hadronizer::decay ( )

Definition at line 446 of file Pythia6Hadronizer.cc.

 {
    return true;
 }

bool gen::Pythia6Hadronizer::declareSpecialSettings ( const std::vector< std::string > settings )

Definition at line 738 of file Pythia6Hadronizer.cc.

References gen::call_pygive(), edm::errors::Configuration, edm::hlt::Exception, spr::find(), gen::pycomp_(), gen::pydat1_, and relativeConstraints::value.

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

bool gen::Pythia6Hadronizer::declareStableParticles ( const std::vector< int > pdg )

Definition at line 716 of file Pythia6Hadronizer.cc.

References gen::call_pygive(), 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;
 }

void gen::Pythia6Hadronizer::fillTmpStorage ( )

private

Definition at line 187 of file Pythia6Hadronizer.cc.

References i, gen::pyjets_local, and mathSSE::return().

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

void gen::Pythia6Hadronizer::finalizeEvent ( )

Definition at line 206 of file Pythia6Hadronizer.cc.

References abs, gen::call_pylist(), gather_cfg::cout, gen::BaseHadronizer::event(), gen::BaseHadronizer::eventInfo(), fConvertToPDG, fHepMCVerbosity, lhef::LHEEvent::fillEventInfo(), lhef::LHEEvent::fillPdfInfo(), fImposeProperTime, fMaxEventsToPrint, fPythiaListVerbosity, imposeProperTime(), gen::BaseHadronizer::lheEvent(), gen::BaseHadronizer::lheRunInfo(), pydat1, pyint1, and pypars.

 {
 
    bool lhe = lheEvent() != 0;
 
    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 ) ;
 
    // 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().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
    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();
       }
    }
    
    return;
 }

void gen::Pythia6Hadronizer::flushTmpStorage ( )

private

Definition at line 169 of file Pythia6Hadronizer.cc.

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

bool gen::Pythia6Hadronizer::generatePartonsAndHadronize ( )

Definition at line 324 of file Pythia6Hadronizer.cc.

References gen::call_pygive(), conv, gen::BaseHadronizer::event(), fGluinoHadronsEnabled, fillTmpStorage(), flushTmpStorage(), fPy6Service, fStopHadronsEnabled, gen::FortranCallback::getInstance(), gen::pyglfr_(), pyint1, gen::pystfr_(), 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( conv.read_next_event() );
    
    // this is to deal with post-gen tools & residualDecay() that may reuse PYJETS
    //
    flushTmpStorage();
    fillTmpStorage();
       
    return true;
 }

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

inlinestatic

Definition at line 54 of file Pythia6Hadronizer.h.

References fJetMatching.

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

54 { return fJetMatching; }

gen::Pythia6Hadronizer::fJetMatching

static JetMatching * fJetMatching

Definition: Pythia6Hadronizer.h:82

bool gen::Pythia6Hadronizer::hadronize ( )

Definition at line 375 of file Pythia6Hadronizer.cc.

References gen::JetMatching::beforeHadronisation(), gen::call_pygive(), conv, 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_(), 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( conv.read_next_event() );
    
    // this is to deal with post-gen tools and/or residualDecay(), that may reuse PYJETS
    //
    flushTmpStorage();
    fillTmpStorage();
       
    return true;
 }

void gen::Pythia6Hadronizer::imposeProperTime ( )

private

Definition at line 781 of file Pythia6Hadronizer.cc.

References abs, gen::FortranInstance::call(), gen::BaseHadronizer::event(), fPy6Service, funct::log(), gen::pycomp_(), pydat1, gen::pyr_(), mathSSE::sqrt(), matplotRender::t, vbegin, vend, x, detailsBasic3DVector::y, and detailsBasic3DVector::z.

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

bool gen::Pythia6Hadronizer::initializeForExternalPartons ( )

Definition at line 631 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::Pythia6Service::setGeneralParams(), gen::FortranCallback::setLHERunInfo(), gen::Pythia6Service::setPYUPDAParams(), and gen::Pythia6Service::setSLHAFromHeader().

 {
    Pythia6Service::InstanceWrapper guard(fPy6Service);  // grab Py6 instance
 
    // note: CSA mode is NOT supposed to woirk 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;
 }

bool gen::Pythia6Hadronizer::initializeForInternalPartons ( )

Definition at line 684 of file Pythia6Hadronizer.cc.

References gen::call_pygive(), gen::Pythia6Service::closeSLHA(), fCOMEnergy, fGluinoHadronsEnabled, fPy6Service, fStopHadronsEnabled, gen::pyglrhad_(), gen::pystrhad_(), gen::Pythia6Service::setCSAParams(), gen::Pythia6Service::setGeneralParams(), gen::Pythia6Service::setPYUPDAParams(), and gen::Pythia6Service::setSLHAParams().

 {
    Pythia6Service::InstanceWrapper guard(fPy6Service);  // grab Py6 instance
     
    fPy6Service->setGeneralParams();   
    fPy6Service->setCSAParams();
    fPy6Service->setSLHAParams();
    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_();
    }
    
    call_pyinit("CMS", "p", "p", fCOMEnergy);
 
    fPy6Service->setPYUPDAParams(true);
    
    fPy6Service->closeSLHA();
    
    return true;
 }

bool gen::Pythia6Hadronizer::residualDecay ( )

Definition at line 451 of file Pythia6Hadronizer.cc.

References gen::BaseHadronizer::event(), fPy6Service, configurableAnalysis::GenParticle, i, dbtoconf::parent, gen::pycomp_(), gen::pydecy_(), gen::pyjets_local, and ntuplemaker::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;
 }

void gen::Pythia6Hadronizer::statistics ( )

Definition at line 847 of file Pythia6Hadronizer.cc.

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