Skip to content

Commit

Permalink
Merge pull request cms-sw#7255 from mkirsano/updateresidualdecay
Browse files Browse the repository at this point in the history 
Update residual decays and verbosity
  • Loading branch information
@cmsbuild @bendavid
cmsbuild authored and bendavid committed Jan 18, 2015
1 parent 279d01f commit de3c7cb
Show file tree
Hide file tree
Showing 5 changed files with 96 additions and 91 deletions.
4 changes: 4 additions & 0 deletions GeneratorInterface/Pythia8Interface/interface/Py8InterfaceBase.h
View file
Expand Up @@ -7,6 +7,8 @@
#include "GeneratorInterface/Core/interface/ParameterCollector.h"
#include "GeneratorInterface/Pythia8Interface/interface/P8RndmEngine.h"

#include "HepMC/IO_AsciiParticles.h"

#include <Pythia8/Pythia.h>
#include <Pythia8Plugins/HepMC2.h>

Expand Down Expand Up @@ -45,7 +47,9 @@ namespace gen {

unsigned int pythiaPylistVerbosity;
bool pythiaHepMCVerbosity;
bool pythiaHepMCVerbosityParticles;
unsigned int maxEventsToPrint;
HepMC::IO_AsciiParticles* ascii_io;

private:

Expand Down
18 changes: 14 additions & 4 deletions GeneratorInterface/Pythia8Interface/plugins/HepMCA2.h
View file
Expand Up @@ -34,7 +34,8 @@ class Pythia8ToHepMCA : public IO_BaseClass {
virtual ~Pythia8ToHepMCA() {;}

// Alternative method to convert Pythia events into HepMC ones.
bool append_event( Pythia8::Event& pyev, GenEvent* evt, GenParticle* rootpart);
bool append_event( Pythia8::Event& pyev, GenEvent* evt, GenParticle* rootpart,
int ibarcode = -1, Pythia8::Info* pyinfo = 0, Pythia8::Settings* pyset = 0);

// Read out values for some switches.
bool print_inconsistency() const {return m_print_inconsistency;}
Expand Down Expand Up @@ -70,7 +71,8 @@ class Pythia8ToHepMCA : public IO_BaseClass {
// Read one event from Pythia8, append GenEvent
// and return T/F = success/failure.

inline bool Pythia8ToHepMCA::append_event( Pythia8::Event& pyev, GenEvent* evt, GenParticle* rootpart) {
inline bool Pythia8ToHepMCA::append_event( Pythia8::Event& pyev, GenEvent* evt, GenParticle* rootpart,
int ibarcode, Pythia8::Info* pyinfo, Pythia8::Settings* pyset) {

// 1. Error if no event passed.
if (!evt) {
Expand All @@ -85,6 +87,9 @@ inline bool Pythia8ToHepMCA::append_event( Pythia8::Event& pyev, GenEvent* evt,
double lenFac = HepMC::Units::conversion_factor(HepMC::Units::MM,
evt->length_unit());

int NewBarcode = evt->particles_size();
if (ibarcode > -1) NewBarcode = ibarcode;

GenVertex* prod_vtx0 = new GenVertex();
prod_vtx0->add_particle_in( rootpart );
evt->add_vertex( prod_vtx0 );
Expand All @@ -100,7 +105,8 @@ inline bool Pythia8ToHepMCA::append_event( Pythia8::Event& pyev, GenEvent* evt,
FourVector( momFac * pyev[i].px(), momFac * pyev[i].py(),
momFac * pyev[i].pz(), momFac * pyev[i].e() ),
pyev[i].id(), pyev[i].statusHepMC() );
hepevt_particles[i]->suggest_barcode(0);
if (ibarcode !=0) NewBarcode++;
hepevt_particles[i]->suggest_barcode(NewBarcode);
hepevt_particles[i]->set_generated_mass( momFac * pyev[i].m() );

// Colour flow uses index 1 and 2.
Expand Down Expand Up @@ -181,6 +187,10 @@ inline bool Pythia8ToHepMCA::append_event( Pythia8::Event& pyev, GenEvent* evt,
}
}

// If hadronization switched on then no final coloured particles.
bool doHadr = (pyset == 0) ? m_free_parton_warnings
: pyset->flag("HadronLevel:all") && pyset->flag("HadronLevel:Hadronize");

// 4. Check for particles which come from nowhere, i.e. are without
// mothers or daughters. These need to be attached to a vertex, or else
// they will never become part of the event.
Expand All @@ -194,7 +204,7 @@ inline bool Pythia8ToHepMCA::append_event( Pythia8::Event& pyev, GenEvent* evt,
}

// Also check for free partons (= gluons and quarks; not diquarks?).
if ( m_free_parton_warnings ) {
if ( doHadr && m_free_parton_warnings ) {
if ( hepevt_particles[i]->pdg_id() == 21 &&
!hepevt_particles[i]->end_vertex() ) {
std::cerr << "gluon without end vertex " << i << std::endl;
Expand Down
11 changes: 9 additions & 2 deletions GeneratorInterface/Pythia8Interface/plugins/Pythia8Hadronizer.cc
View file
Expand Up @@ -583,7 +583,7 @@ bool Pythia8Hadronizer::residualDecay()

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

if ( part->status() == 1 )
if ( part->status() == 1 && (fDecayer->particleData).canDecay(part->pdg_id()) )
{
fDecayer->event.reset();
Particle py8part( part->pdg_id(), 93, 0, 0, 0, 0, 0, 0,
Expand Down Expand Up @@ -637,7 +637,8 @@ void Pythia8Hadronizer::finalizeEvent()
//******** Verbosity ********

if (maxEventsToPrint > 0 &&
(pythiaPylistVerbosity || pythiaHepMCVerbosity)) {
(pythiaPylistVerbosity || pythiaHepMCVerbosity ||
pythiaHepMCVerbosityParticles) ) {
maxEventsToPrint--;
if (pythiaPylistVerbosity) {
fMasterGen->info.list(std::cout);
Expand All @@ -650,6 +651,12 @@ void Pythia8Hadronizer::finalizeEvent()
<< "----------------------" << std::endl;
event()->print();
}
if (pythiaHepMCVerbosityParticles) {
std::cout << "Event process = "
<< fMasterGen->info.code() << "\n"
<< "----------------------" << std::endl;
ascii_io->write_event(event().get());
}
}
}

Expand Down
151 changes: 66 additions & 85 deletions GeneratorInterface/Pythia8Interface/src/Py8GunBase.cc
View file
@@ -1,6 +1,8 @@
#include "GeneratorInterface/Pythia8Interface/interface/Py8GunBase.h"
#include "FWCore/Concurrency/interface/SharedResourceNames.h"

#include "GeneratorInterface/Pythia8Interface/plugins/HepMCA2.h"

using namespace Pythia8;

const std::vector<std::string> gen::Py8GunBase::p8SharedResources = { edm::SharedResourceNames::kPythia8 };
Expand Down Expand Up @@ -60,119 +62,98 @@ bool Py8GunBase::initializeForInternalPartons()
bool Py8GunBase::residualDecay()
{

Event* pythiaEvent = &(fMasterGen->event);
Event* pythiaEvent = &(fMasterGen->event);

int NPartsBeforeDecays = pythiaEvent->size()-1; // do NOT count the very 1st "system" particle
// in Pythia8::Event record; it does NOT even
// get translated by the HepMCInterface to the
// HepMC::GenEvent record !!!
int NPartsAfterDecays = event().get()->particles_size();
int NewBarcode = NPartsAfterDecays;

for ( int ipart=NPartsAfterDecays; ipart>NPartsBeforeDecays; ipart-- )
{
int NPartsBeforeDecays = pythiaEvent->size()-1; // do NOT count the very 1st "system" particle
// in Pythia8::Event record; it does NOT even
// get translated by the HepMCInterface to the
// HepMC::GenEvent record !!!
int NPartsAfterDecays = event().get()->particles_size();

if(NPartsAfterDecays == NPartsBeforeDecays) return true;

HepMC::Pythia8ToHepMCA toHepMCA;
bool result = true;

for ( int ipart=NPartsAfterDecays; ipart>NPartsBeforeDecays; ipart-- )
{

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

if ( part->status() == 1 )
{
fDecayer->event.reset();
Particle py8part( part->pdg_id(), 93, 0, 0, 0, 0, 0, 0,
if ( part->status() == 1 && (fDecayer->particleData).canDecay(part->pdg_id()) )
{
fDecayer->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(),
HepMC::GenVertex* ProdVtx = part->production_vertex();
py8part.vProd( ProdVtx->position().x(), ProdVtx->position().y(),
ProdVtx->position().z(), ProdVtx->position().t() );
py8part.tau( (fDecayer->particleData).tau0( part->pdg_id() ) );
fDecayer->event.append( py8part );
int nentries = fDecayer->event.size();
if ( !fDecayer->event[nentries-1].mayDecay() ) continue;
fDecayer->next();
int nentries1 = fDecayer->event.size();
if ( nentries1 <= nentries ) continue; //same number of particles, no decays...

part->set_status(2);
py8part.tau( (fDecayer->particleData).tau0( part->pdg_id() ) );
fDecayer->event.append( py8part );
int nentries = fDecayer->event.size();
if ( !fDecayer->event[nentries-1].mayDecay() ) continue;
fDecayer->next();
int nentries1 = fDecayer->event.size();
if ( nentries1 <= nentries ) continue; //same number of particles, no decays...

Particle& py8daughter = fDecayer->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
part->set_status(2);

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 ( int ipart1=nentries+1; ipart1<nentries1; ipart1++ )
{
py8daughter = fDecayer->event[ipart1];
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 );
result = toHepMCA.append_event( fDecayer->event, event().get(), part);

}
}
}
}

return true;
return result;

}

void Py8GunBase::finalizeEvent()
{


// FIXME !!!

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

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

if (pythiaHepMCVerbosity)
{
std::cout << "Event process = "
<< fMasterGen->info.code() << "\n"
<< "----------------------" << std::endl;
event()->print();
}
}
if (maxEventsToPrint > 0 &&
(pythiaPylistVerbosity || pythiaHepMCVerbosity ||
pythiaHepMCVerbosityParticles) )
{
maxEventsToPrint--;
if (pythiaPylistVerbosity)
{
fMasterGen->info.list(std::cout);
fMasterGen->event.list(std::cout);
}

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

return;
return;
}

void Py8GunBase::statistics()
{

fMasterGen->stat();
fMasterGen->stat();

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

}

Expand Down
3 changes: 3 additions & 0 deletions GeneratorInterface/Pythia8Interface/src/Py8InterfaceBase.cc
View file
Expand Up @@ -22,8 +22,11 @@ Py8InterfaceBase::Py8InterfaceBase( edm::ParameterSet const& ps )

pythiaPylistVerbosity = ps.getUntrackedParameter<int>("pythiaPylistVerbosity", 0);
pythiaHepMCVerbosity = ps.getUntrackedParameter<bool>("pythiaHepMCVerbosity", false);
pythiaHepMCVerbosityParticles = ps.getUntrackedParameter<bool>("pythiaHepMCVerbosityParticles", false);
maxEventsToPrint = ps.getUntrackedParameter<int>("maxEventsToPrint", 0);

if(pythiaHepMCVerbosityParticles)
ascii_io = new HepMC::IO_AsciiParticles("cout", std::ios::out);
}

bool Py8InterfaceBase::readSettings( int )
Expand Down

0 comments on commit de3c7cb

Please sign in to comment.