Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

get rid of option to use pythia6 to decay particles in FastSim #7155

Merged
merged 1 commit into from Jan 13, 2015
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
49 changes: 0 additions & 49 deletions FastSimulation/ParticleDecay/doc/ParticleDecay.doc

This file was deleted.

34 changes: 0 additions & 34 deletions FastSimulation/ParticleDecay/interface/Pythia6jets.h

This file was deleted.

37 changes: 12 additions & 25 deletions FastSimulation/ParticleDecay/interface/PythiaDecays.h
Expand Up @@ -4,51 +4,38 @@
#include <memory>
#include <vector>

// Needed for Pythia8
#include <HepMC/GenEvent.h>
#include <HepMC/GenParticle.h>
#include <Pythia8/Pythia.h>
#include "Pythia8Plugins/HepMC2.h"
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
#include "GeneratorInterface/ExternalDecays/interface/ExternalDecayDriver.h"
// TrajectoryManager does not compile when forward declaring P8RndmEngine
#include "GeneratorInterface/Pythia8Interface/interface/P8RndmEngine.h"

class ParticlePropagator;
class Pythia6jets;
class RawParticle;

namespace gen {
class P8RndmEngine;
}

namespace CLHEP {
class HepRandomEngine;
}

namespace gen {
class P8RndmEngine;
namespace Pythia8 {
class Pythia;
}

typedef std::vector<RawParticle> DaughterParticleList;
typedef DaughterParticleList::const_iterator DaughterParticleIterator;

namespace gen { class Pythia6Service; } // remove?
//using namespace gen;

class PythiaDecays
{
public:
PythiaDecays(std::string program);
~PythiaDecays();

const DaughterParticleList&
particleDaughtersPy6(ParticlePropagator& particle, CLHEP::HepRandomEngine*);
const DaughterParticleList&
particleDaughtersPy8(ParticlePropagator& particle, CLHEP::HepRandomEngine*);
PythiaDecays();
~PythiaDecays(){};
const DaughterParticleList & particleDaughters(ParticlePropagator& particle, CLHEP::HepRandomEngine*);

private:

DaughterParticleList theList;
std::string program_;
// for Pythia6:
gen::Pythia6Service *pyservice;
Pythia6jets* pyjets;
// for Pythia8:
std::auto_ptr<Pythia8::Pythia> decayer;
std::unique_ptr<gen::P8RndmEngine> p8RndmEngine;
};
Expand Down
1 change: 0 additions & 1 deletion FastSimulation/ParticleDecay/interface/RandomP8.h

This file was deleted.

83 changes: 0 additions & 83 deletions FastSimulation/ParticleDecay/src/Pythia6jets.cc

This file was deleted.

95 changes: 6 additions & 89 deletions FastSimulation/ParticleDecay/src/PythiaDecays.cc
@@ -1,37 +1,12 @@
// HepMC Headers
#include "HepMC/PythiaWrapper6_4.h"

// Pythia6 framework integration service Headers
#include "GeneratorInterface/Pythia6Interface/interface/Pythia6Service.h"

// FAMOS Headers
#include "FastSimulation/ParticlePropagator/interface/ParticlePropagator.h"
#include "FastSimulation/ParticleDecay/interface/PythiaDecays.h"
#include "FastSimulation/ParticleDecay/interface/Pythia6jets.h"

#include "GeneratorInterface/Pythia8Interface/interface/P8RndmEngine.h"

#include "FWCore/ServiceRegistry/interface/RandomEngineSentry.h"

// Needed for Pythia6
#define PYTHIA6PYDECY pythia6pydecy_
#include <Pythia8/Pythia.h>
#include "Pythia8Plugins/HepMC2.h"

extern "C" {
void PYTHIA6PYDECY(int *ip);
}

PythiaDecays::PythiaDecays(std::string program)
PythiaDecays::PythiaDecays()
{
program_=program;
if (program_ == "pythia6") {
//// Pythia6:
pyjets = new Pythia6jets();
pyservice = new gen::Pythia6Service();
// The PYTHIA decay tables will be initialized later
} else if (program_ == "pythia8") {

//// Pythia8:

// inspired by method Pythia8Hadronizer::residualDecay() in GeneratorInterface/Pythia8Interface/src/Py8GunBase.cc
decayer.reset(new Pythia8::Pythia);
p8RndmEngine.reset(new gen::P8RndmEngine);
Expand All @@ -41,29 +16,18 @@ PythiaDecays::PythiaDecays(std::string program)
decayer->settings.flag("ProcessLevel:resonanceDecays",false);
decayer->init();

// forbid all decays
// forbid all decays
// (decays are allowed selectively in the particleDaughters function)
Pythia8::ParticleData & pdt = decayer->particleData;
int pid = 1;
while(pdt.nextId(pid) > pid){
pid = pdt.nextId(pid);
pdt.mayDecay(pid,false);
}

} else {
std::cout << "WARNING: you are requesting an option which is not available in PythiaDecays::PythiaDecays " << std::endl;
}

}

PythiaDecays::~PythiaDecays() {
if (program_ == "pythia6") {
delete pyjets;
delete pyservice;
}
}

const DaughterParticleList&
PythiaDecays::particleDaughtersPy8(ParticlePropagator& particle, CLHEP::HepRandomEngine* engine)
PythiaDecays::particleDaughters(ParticlePropagator& particle, CLHEP::HepRandomEngine* engine)
{
edm::RandomEngineSentry<gen::P8RndmEngine> sentry(p8RndmEngine.get(), engine);

Expand Down Expand Up @@ -106,50 +70,3 @@ PythiaDecays::particleDaughtersPy8(ParticlePropagator& particle, CLHEP::HepRando

return theList;
}

const DaughterParticleList&
PythiaDecays::particleDaughtersPy6(ParticlePropagator& particle, CLHEP::HepRandomEngine* engine)
{
edm::RandomEngineSentry<gen::Pythia6Service> sentry(pyservice, engine);

gen::Pythia6Service::InstanceWrapper guard(pyservice); // grab Py6 context

// Pythia6jets pyjets;
int ip;

pyjets->k(1,1) = 1;
pyjets->k(1,2) = particle.pid();
pyjets->p(1,1) = particle.Px();
pyjets->p(1,2) = particle.Py();
pyjets->p(1,3) = particle.Pz();
pyjets->p(1,4) = std::max(particle.mass(),particle.e());
pyjets->p(1,5) = particle.mass();
pyjets->v(1,1) = particle.X();
pyjets->v(1,2) = particle.Y();
pyjets->v(1,3) = particle.Z();
pyjets->v(1,4) = particle.T();
pyjets->n() = 1;

ip = 1;
PYTHIA6PYDECY(&ip);

// Fill the list of daughters
theList.clear();
if ( pyjets->n()==1 ) return theList;

theList.resize(pyjets->n()-1,RawParticle());

for (int i=2;i<=pyjets->n();++i) {

theList[i-2].SetXYZT(pyjets->p(i,1),pyjets->p(i,2),
pyjets->p(i,3),pyjets->p(i,4));
theList[i-2].setVertex(pyjets->v(i,1),pyjets->v(i,2),
pyjets->v(i,3),pyjets->v(i,4));
theList[i-2].setID(pyjets->k(i,2));
theList[i-2].setMass(pyjets->p(i,5));

}

return theList;

}