Inheritance diagram for Pythia8Hadronizer:

Public Member Functions
const char *	classname () const
bool	decay ()
bool	declareSpecialSettings (const std::vector< std::string >)
bool	declareStableParticles (const std::vector< int > &pdgIds)
void	finalizeEvent ()
bool	generatePartonsAndHadronize ()
bool	hadronize ()
bool	initializeForExternalPartons ()
bool	initializeForInternalPartons ()
	Pythia8Hadronizer (const edm::ParameterSet &params)
bool	readSettings (int)
bool	residualDecay ()
void	statistics ()
	~Pythia8Hadronizer ()
Private Types
enum	{ PP, PPbar, ElectronPositron }
Private Attributes
double	comEnergy
	Center-of-Mass energy.
std::auto_ptr< Pythia >	decayer
double	fBeam1PZ
double	fBeam2PZ
EmissionVetoHook *	fEmissionVetoHook
int	fInitialState
JetMatchingHook *	fJetMatchingHook
UserHooks *	fReweightUserHook
std::auto_ptr< LHAupLesHouches >	lhaUP
string	LHEInputFileName
unsigned int	maxEventsToPrint
	Events to print if verbosity.
ParameterCollector	parameters
std::auto_ptr< Pythia >	pythia
Event *	pythiaEvent
bool	pythiaHepMCVerbosity
	HepMC verbosity flag.
unsigned int	pythiaPylistVerbosity
	Pythia PYLIST Verbosity flag.
HepMC::I_Pythia8	toHepMC

Detailed Description

Definition at line 47 of file Pythia8Hadronizer.cc.

Member Enumeration Documentation

anonymous enum [private]

Enumerator:

PP
PPbar
ElectronPositron

Definition at line 90 of file Pythia8Hadronizer.cc.

{ PP, PPbar, ElectronPositron };

Constructor & Destructor Documentation

Pythia8Hadronizer::Pythia8Hadronizer ( const edm::ParameterSet & params )

Definition at line 111 of file Pythia8Hadronizer.cc.

References edm::errors::Configuration, gather_cfg::cout, decayer, ElectronPositron, Exception, edm::ParameterSet::exists(), fEmissionVetoHook, fInitialState, fJetMatchingHook, fReweightUserHook, gen::getEngineReference(), edm::ParameterSet::getUntrackedParameter(), PP, PPbar, pythia, and randomEngine.

                                                                  :
  BaseHadronizer(params),
  parameters(params.getParameter<edm::ParameterSet>("PythiaParameters")),
  comEnergy(params.getParameter<double>("comEnergy")),
  pythiaPylistVerbosity(params.getUntrackedParameter<int>("pythiaPylistVerbosity", 0)),
  pythiaHepMCVerbosity(params.getUntrackedParameter<bool>("pythiaHepMCVerbosity", false)),
  maxEventsToPrint(params.getUntrackedParameter<int>("maxEventsToPrint", 0)),
  LHEInputFileName(params.getUntrackedParameter<string>("LHEInputFileName","")),
  fInitialState(PP),
  fReweightUserHook(0),
  fJetMatchingHook(0),
  fEmissionVetoHook(0)
{
  randomEngine = &getEngineReference();

  //Old code that used Pythia8 own random engine
  //edm::Service<edm::RandomNumberGenerator> rng;
  //uint32_t seed = rng->mySeed();
  //Pythia8::Rndm::init(seed);

  RandomP8* RP8 = new RandomP8();

  // J.Y.: the following 3 parameters are hacked "for a reason"
  //
  if ( params.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 ( params.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 ( params.exists( "ElectronProtonInitialState" ) || params.exists( "PositronProtonInitialState" ) )
  {
    // throw on unknown initial state !
    throw edm::Exception(edm::errors::Configuration,"Pythia8Interface")
      <<" UNKNOWN INITIAL STATE. \n The allowed initial states are: PP, PPbar, ElectronPositron \n";
  }

  pythia.reset(new Pythia);
  decayer.reset(new Pythia);

  pythia->setRndmEnginePtr(RP8);
  decayer->setRndmEnginePtr(RP8);
    

  // Reweight user hook
  //
  if( params.exists( "reweightGen" ) )
    fReweightUserHook = new PtHatReweightUserHook();

  // PS matching prototype
  //
  if ( params.exists("jetMatching") )
  {
    edm::ParameterSet jmParams =
      params.getUntrackedParameter<edm::ParameterSet>("jetMatching");
    fJetMatchingHook = new JetMatchingHook( jmParams, &pythia->info );
  }

  // Emission veto
  //
  if ( params.exists("emissionVeto") )
  {   
    fEmissionVetoHook = new EmissionVetoHook(0);
    pythia->setUserHooksPtr( fEmissionVetoHook );
  }  

  int NHooks=0;
  if(fReweightUserHook) NHooks++;
  if(fJetMatchingHook) NHooks++;
  if(fEmissionVetoHook) NHooks++;
  if(NHooks > 1)
    throw edm::Exception(edm::errors::Configuration,"Pythia8Interface")
      <<" Too many User Hooks. \n Please choose one from: reweightGen, jetMatching, emissionVeto \n";

  if(fReweightUserHook) pythia->setUserHooksPtr(fReweightUserHook);
  if(fJetMatchingHook) pythia->setUserHooksPtr(fJetMatchingHook);
  if(fEmissionVetoHook) pythia->setUserHooksPtr(fEmissionVetoHook);
}

Pythia8Hadronizer::~Pythia8Hadronizer ( )

Definition at line 217 of file Pythia8Hadronizer.cc.

References fEmissionVetoHook.

{
// do we need to delete UserHooks/JetMatchingHook here ???

  if(fEmissionVetoHook) {delete fEmissionVetoHook; fEmissionVetoHook=0;}
}

Member Function Documentation

const char* Pythia8Hadronizer::classname ( ) const [inline]

Definition at line 67 of file Pythia8Hadronizer.cc.

{ return "Pythia8Hadronizer"; }

bool Pythia8Hadronizer::decay ( )

Definition at line 421 of file Pythia8Hadronizer.cc.

{
   return true;
}

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

Definition at line 358 of file Pythia8Hadronizer.cc.

References spr::find(), and pythia.

{
  for ( unsigned int iss=0; iss<settings.size(); iss++ )
  {
    if ( settings[iss].find("QED-brem-off") == std::string::npos ) continue;
    pythia->readString( "TimeShower:QEDshowerByL=off" );
  }

  return true;
}

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

Definition at line 330 of file Pythia8Hadronizer.cc.

References decayer, i, and pythia.

{
  for ( size_t i=0; i<pdgIds.size(); i++ )
  {
    // FIXME: need to double check if PID's are the same in Py6 & Py8,
    //        because the HepPDT translation tool is actually for **Py6** 
    // 
    // well, actually it looks like Py8 operates in PDT id's rather than Py6's
    //
    // int PyID = HepPID::translatePDTtoPythia( pdgIds[i] ); 
    int PyID = pdgIds[i]; 
    std::ostringstream pyCard ;
    pyCard << PyID <<":mayDecay=false";
    pythia->readString( pyCard.str() );
    // alternative:
    // set the 2nd input argument warn=false 
    // - this way Py8 will NOT print warnings about unknown particle code(s)
    // pythia->readString( pyCard.str(), false )
  }
      
  // init decayer
  decayer->readString("ProcessLevel:all = off"); // The trick!
  decayer->init();
   
  return true;
}

void Pythia8Hadronizer::finalizeEvent ( )

Definition at line 500 of file Pythia8Hadronizer.cc.

References abs, gather_cfg::cout, gen::BaseHadronizer::event(), gen::BaseHadronizer::eventInfo(), gen::BaseHadronizer::lheEvent(), gen::BaseHadronizer::lheRunInfo(), maxEventsToPrint, pythia, pythiaHepMCVerbosity, and pythiaPylistVerbosity.

{
  bool lhe = lheEvent() != 0;

  event()->set_signal_process_id(pythia->info.code());
  event()->set_event_scale(pythia->info.pTHat());       //FIXME

  //cout.precision(10);
  //cout << " pt = " << pythia->info.pTHat() << " weights = "
  //     << pythia->info.weight() << " "
  //     << fReweightUserHook->biasedSelectionWeight() << endl;

  if (event()->alphaQED() <= 0)
    event()->set_alphaQED( pythia->info.alphaEM() );
  if (event()->alphaQCD() <= 0)
    event()->set_alphaQCD( pythia->info.alphaS() );

  HepMC::GenCrossSection xsec;
  xsec.set_cross_section( pythia->info.sigmaGen() * 1e9,
                          pythia->info.sigmaErr() * 1e9);
  event()->set_cross_section(xsec);

  // Putting pdf info into the HepMC record
  // There is the overloaded pythia8 HepMCInterface method fill_next_event
  // that does this, but CMSSW GeneratorInterface does not fill HepMC
  // record according to the HepMC convention (stores f(x) instead of x*f(x)
  // and converts gluon PDG ID to zero). For this reason we use the
  // method fill_next_event (above) that does NOT this and fill pdf info here
  //
  int id1 = pythia->info.id1();
  int id2 = pythia->info.id2();
  if (id1 == 21) id1 = 0;
  if (id2 == 21) id2 = 0;
  double x1 = pythia->info.x1();
  double x2 = pythia->info.x2();
  //double Q = pythia->info.QRen();
  double Q = pythia->info.QFac();
  double pdf1 = pythia->info.pdf1() / pythia->info.x1();
  double pdf2 = pythia->info.pdf2() / pythia->info.x2();
  event()->set_pdf_info(HepMC::PdfInfo(id1,id2,x1,x2,Q,pdf1,pdf2));

  // Storing weights. Will be moved to pythia8 HepMCInterface
  //
  if (lhe && std::abs(lheRunInfo()->getHEPRUP()->IDWTUP) == 4)
    // translate mb to pb (CMS/Gen "convention" as of May 2009)
    event()->weights().push_back( pythia->info.weight() * 1.0e9 );
  else
    event()->weights().push_back( pythia->info.weight() );

  // now create the GenEventInfo product from the GenEvent and fill
  // the missing pieces
  eventInfo().reset( new GenEventInfoProduct( event().get() ) );

  // in pythia 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, pythia->info.pTHat()));
  }

  //******** Verbosity ********

  if (maxEventsToPrint > 0 &&
      (pythiaPylistVerbosity || pythiaHepMCVerbosity)) {
    maxEventsToPrint--;
    if (pythiaPylistVerbosity) {
      pythia->info.list(std::cout); 
      pythia->event.list(std::cout);
    } 

    if (pythiaHepMCVerbosity) {
      std::cout << "Event process = "
                << pythia->info.code() << "\n"
                << "----------------------" << std::endl;
      event()->print();
    }
  }
}

bool Pythia8Hadronizer::generatePartonsAndHadronize ( )

Definition at line 380 of file Pythia8Hadronizer.cc.

References gen::BaseHadronizer::event(), pythia, pythiaEvent, and toHepMC.

{
  if (!pythia->next()) return false;

  event().reset(new HepMC::GenEvent);
  toHepMC.fill_next_event(*pythiaEvent, event().get());

  return true;
}

bool Pythia8Hadronizer::hadronize ( )

Definition at line 391 of file Pythia8Hadronizer.cc.

References JetMatchingHook::beforeHadronization(), lhef::LHEEvent::count(), gen::BaseHadronizer::event(), fJetMatchingHook, lhef::LHERunInfo::kAccepted, lhef::LHERunInfo::kSelected, lhaUP, gen::BaseHadronizer::lheEvent(), LHEInputFileName, pythia, pythiaEvent, JetMatchingHook::resetMatchingStatus(), and toHepMC.

{
  if(LHEInputFileName == string()) lhaUP->loadEvent(lheEvent());

  if ( fJetMatchingHook ) 
  {
    fJetMatchingHook->resetMatchingStatus(); 
    fJetMatchingHook->beforeHadronization( lheEvent() );
  }

  bool py8next = pythia->next();
  // if (!pythia->next())
  if (!py8next)
  {
    lheEvent()->count( lhef::LHERunInfo::kSelected );
    event().reset();
    return false;
  }

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

  event().reset(new HepMC::GenEvent);
  toHepMC.fill_next_event(*pythiaEvent, event().get());

  return true;
}

bool Pythia8Hadronizer::initializeForExternalPartons ( )

Definition at line 280 of file Pythia8Hadronizer.cc.

References gather_cfg::cout, fJetMatchingHook, JetMatchingHook::init(), lhaUP, LHEInputFileName, gen::BaseHadronizer::lheRunInfo(), pythia, and pythiaEvent.

{

  std::cout << "Initializing for external partons" << std::endl;

  pythiaEvent = &pythia->event;
    
  if(LHEInputFileName != string()) {

    cout << endl;
    cout << "LHE Input File Name = " << LHEInputFileName << endl;
    cout << endl;
    pythia->init(LHEInputFileName);

  } else {

    lhaUP.reset(new LHAupLesHouches());
    lhaUP->loadRunInfo(lheRunInfo());
    pythia->init(lhaUP.get());

  }

  // PS matching prototype
  //
  if ( fJetMatchingHook ) 
  {
    // matcher will be init as well, inside init(...)
    //
    fJetMatchingHook->init ( lheRunInfo() );
  }

    return true;
}

bool Pythia8Hadronizer::initializeForInternalPartons ( )

Definition at line 250 of file Pythia8Hadronizer.cc.

References comEnergy, edm::errors::Configuration, ElectronPositron, Exception, fInitialState, PP, PPbar, pythia, and pythiaEvent.

{

  pythiaEvent = &pythia->event;

  if ( fInitialState == PP ) // default
  {
    pythia->init(2212, 2212, comEnergy);
  }
  else if ( fInitialState == PPbar )
  {
    pythia->init(2212, -2212, comEnergy);
  }
  else if ( fInitialState == ElectronPositron )
  {
    pythia->init(11, -11, comEnergy);
  }    
  else 
  {
    // throw on unknown initial state !
    throw edm::Exception(edm::errors::Configuration,"Pythia8Interface")
      <<" UNKNOWN INITIAL STATE. \n The allowed initial states are: PP, PPbar, ElectronPositron \n";
  }

  pythia->settings.listChanged();

  return true;
}

bool Pythia8Hadronizer::readSettings ( int )

Definition at line 225 of file Pythia8Hadronizer.cc.

References gen::ParameterCollector::begin(), gen::ParameterCollector::end(), geometryCSVtoXML::line, parameters, pythia, and pythiaPylistVerbosity.

{
  for ( ParameterCollector::const_iterator line = parameters.begin();
        line != parameters.end(); ++line ) {
    if (line->find("Random:") != std::string::npos)
      throw cms::Exception("PythiaError") << "Attempted to set random number "
        "using Pythia commands. Please use " "the RandomNumberGeneratorService."
        << std::endl;

    if (!pythia->readString(*line)) throw cms::Exception("PythiaError")
                                      << "Pythia 8 did not accept \""
                                      << *line << "\"." << std::endl;
  }

  if ( pythiaPylistVerbosity > 10 ) {
    if ( pythiaPylistVerbosity == 11 || pythiaPylistVerbosity == 13 )
      pythia->settings.listAll();
    if ( pythiaPylistVerbosity == 12 || pythiaPylistVerbosity == 13 )
      pythia->particleData.listAll();
  }

  return true;
}

bool Pythia8Hadronizer::residualDecay ( )

Definition at line 427 of file Pythia8Hadronizer.cc.

References decayer, gen::BaseHadronizer::event(), configurableAnalysis::GenParticle, and pythiaEvent.

{

  int NPartsBeforeDecays = pythiaEvent->size();
  int NPartsAfterDecays = event().get()->particles_size();
  int NewBarcode = NPartsAfterDecays;
   
  for ( int ipart=NPartsAfterDecays; ipart>NPartsBeforeDecays; ipart-- )
  {

    HepMC::GenParticle* part = event().get()->barcode_to_particle( ipart );

    if ( part->status() == 1 )
    {
      decayer->event.reset();
      Particle py8part(  part->pdg_id(), 93, 0, 0, 0, 0, 0, 0,
                         part->momentum().x(),
                         part->momentum().y(),
                         part->momentum().z(),
                         part->momentum().t(),
                         part->generated_mass() );
      HepMC::GenVertex* ProdVtx = part->production_vertex();
      py8part.vProd( ProdVtx->position().x(), ProdVtx->position().y(), 
                     ProdVtx->position().z(), ProdVtx->position().t() );
      py8part.tau( (decayer->particleData).tau0( part->pdg_id() ) );
      decayer->event.append( py8part );
      int nentries = decayer->event.size();
      if ( !decayer->event[nentries-1].mayDecay() ) continue;
      decayer->next();
      int nentries1 = decayer->event.size();
      // --> decayer->event.list(std::cout);
      if ( nentries1 <= nentries ) continue; //same number of particles, no decays...
            
      part->set_status(2);
            
      Particle& py8daughter = decayer->event[nentries]; // the 1st daughter
      HepMC::GenVertex* DecVtx = new HepMC::GenVertex( HepMC::FourVector(py8daughter.xProd(),
                                                       py8daughter.yProd(),
                                                       py8daughter.zProd(),
                                                       py8daughter.tProd()) );

      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( py8daughter.px(), py8daughter.py(), py8daughter.pz(), py8daughter.e() );
            
      HepMC::GenParticle* daughter =
                        new HepMC::GenParticle( pmom, py8daughter.id(), 1 );
            
      NewBarcode++;
      daughter->suggest_barcode( NewBarcode );
      DecVtx->add_particle_out( daughter );
                    
      for ( ipart=nentries+1; ipart<nentries1; ipart++ )
      {
        py8daughter = decayer->event[ipart];
        HepMC::FourVector pmomN( py8daughter.px(), py8daughter.py(), py8daughter.pz(), py8daughter.e() );           
        HepMC::GenParticle* daughterN =
                        new HepMC::GenParticle( pmomN, py8daughter.id(), 1 );
        NewBarcode++;
        daughterN->suggest_barcode( NewBarcode );
        DecVtx->add_particle_out( daughterN );
      }
            
      event().get()->add_vertex( DecVtx );

    }
 } 
   
 return true;
}

void Pythia8Hadronizer::statistics ( )

Definition at line 370 of file Pythia8Hadronizer.cc.

References pythia, gen::BaseHadronizer::runInfo(), and GenRunInfoProduct::setInternalXSec().

{
  pythia->statistics();

  double xsec = pythia->info.sigmaGen(); // cross section in mb
  xsec *= 1.0e9; // translate to pb (CMS/Gen "convention" as of May 2009)
  runInfo().setInternalXSec(xsec);
}