10 #include "HepMC/PythiaWrapper.h"
11 #include "HepMC/IO_HEPEVT.h"
15 namespace ParticleReplacerGunVar{
23 particleOrigin_(iConfig.getParameter<std::
string>(
"particleOrigin")),
24 forceTauPolarization_(iConfig.getParameter<std::
string>(
"forceTauPolarization")),
25 forceTauDecay_(iConfig.getParameter<std::
string>(
"forceTauDecay")),
26 generatorMode_(iConfig.getParameter<std::
string>(
"generatorMode")),
27 gunParticle_(iConfig.getParameter<int>(
"gunParticle")),
28 forceTauPlusHelicity_(iConfig.getParameter<int>(
"forceTauPlusHelicity")),
29 forceTauMinusHelicity_(iConfig.getParameter<int>(
"forceTauMinusHelicity")),
37 <<
"The RandomNumberProducer module requires the RandomNumberGeneratorService\n"
38 "which appears to be absent. Please add that service to your configuration\n"
39 "or remove the modules that require it." << std::endl;
46 edm::LogInfo(
"MuonReplacement") <<
"[ParticleReplacer::ParticleReplacer] "
49 edm::LogInfo(
"MuonReplacement") <<
"[ParticleReplacer::ParticleReplacer] "
52 edm::LogInfo(
"MuonReplacement") <<
"[ParticleReplacer::ParticleReplacer] "
55 edm::LogInfo(
"MuonReplacement") <<
"[ParticleReplacer::ParticleReplacer] "
58 std::memset(
pol1_, 0, 4*
sizeof(
float));
59 std::memset(
pol2_, 0, 4*
sizeof(
float));
62 throw cms::Exception(
"Configuration") <<
"Generator mode other than Tauola is not supported" << std::endl;
64 throw cms::Exception(
"UnimplementedFeature") <<
"ParticleReplacerParticleGun is not usable yet." << std::endl;
91 throw cms::Exception(
"UnimplementedFeature") <<
"ParticleReplacerParticleGun does NOT support merging at HepMC level" << std::endl;
93 std::auto_ptr<HepMC::GenEvent> evt(0);
94 std::vector<HepMC::FourVector> muons_corrected;
95 muons_corrected.reserve(muons.size());
100 for(
unsigned int i=0;
i<muons_corrected.size(); ++
i) {
101 HepMC::FourVector&
muon = muons_corrected[
i];
108 std::cout <<
" ///////////////////// After py1ent, before pyhepc /////////////////////" << std::endl;
116 std::cout <<
" ///////////////////// After pyhepc, before vertex shift /////////////////////" << std::endl;
117 HepMC::HEPEVT_Wrapper::print_hepevt();
120 int nparticles = HepMC::HEPEVT_Wrapper::number_entries();
121 HepMC::ThreeVector
shift(0,0,0);
124 throw cms::Exception(
"LogicError") <<
"Particle origin set to primaryVertex, but pvtx is null!" << std::endl;
126 shift.setX(pvtx->
x()*10);
127 shift.setY(pvtx->
y()*10);
128 shift.setZ(pvtx->
z()*10);
130 for(
int i=1;
i <= nparticles; ++
i) {
132 throw cms::Exception(
"LogicError") <<
"Particle in HEPEVT is " << HepMC::HEPEVT_Wrapper::id(
i)
134 <<
" for index " <<
i << std::endl;
139 shift.setX(muon.
vx()*10);
140 shift.setY(muon.
vy()*10);
141 shift.setZ(muon.
vz()*10);
144 HepMC::HEPEVT_Wrapper::set_position(
i,
152 std::cout <<
" ///////////////////// After vertex shift, before pyhepc/tauola /////////////////////" << std::endl;
153 HepMC::HEPEVT_Wrapper::print_hepevt();
165 for(
unsigned int i=0;
i<muons_corrected.size(); ++
i) {
170 std::cout <<
" ///////////////////// After tauola, before pyhepc /////////////////////" << std::endl;
171 HepMC::HEPEVT_Wrapper::print_hepevt();
175 std::cout <<
" ///////////////////// After pyhepc, before vertex fix /////////////////////" << std::endl;
190 std::cout <<
" ///////////////////// After vertex fix, before pyexec /////////////////////" << std::endl;
197 std::cout <<
" ///////////////////// After pyhepc, before pyexec /////////////////////" << std::endl;
205 std::cout <<
" ///////////////////// After pyexec, before pyhepc /////////////////////" << std::endl;
211 HepMC::IO_HEPEVT
conv;
212 evt = std::auto_ptr<HepMC::GenEvent>(
new HepMC::GenEvent(*conv.read_next_event()));
217 std::cout << std::endl <<
"Vertices: " << std::endl;
218 for(HepMC::GenEvent::vertex_const_iterator iter = evt->vertices_begin(); iter != evt->vertices_end(); ++iter) {
219 std::cout <<
"Vertex " << (*iter)->id() <<
", barcode " << (*iter)->barcode() << std::endl;
221 HepMC::ThreeVector
point = (*iter)->point3d();
222 std::cout <<
" Point (" << point.x() <<
", " << point.y() <<
", " << point.z() <<
")" << std::endl;
225 for(HepMC::GenVertex::particles_in_const_iterator
pi = (*iter)->particles_in_const_begin();
pi != (*iter)->particles_in_const_end(); ++
pi) {
230 for(HepMC::GenVertex::particles_out_const_iterator
pi = (*iter)->particles_out_const_begin();
pi != (*iter)->particles_out_const_end(); ++
pi) {
243 for(reco::MuonCollection::const_iterator iMuon = muons.begin(); iMuon != muons.end(); ++iMuon) {
246 double E =
sqrt(
tauMass*
tauMass + vec.x()*vec.x() + vec.y()*vec.y() + vec.z()*vec.z());
248 corrected.push_back(HepMC::FourVector(vec.x(), vec.y(), vec.z(), E));
253 for(reco::MuonCollection::const_iterator iMuon = muons.begin(); iMuon != muons.end(); ++iMuon) {
254 corrected.push_back(iMuon->p4());
304 edm::LogError(
"MuonReplacement") <<
"[ParticleReplacerAlgoParticleGun::forceTauoladecays] "
305 <<
"Unknown value for forcing tau decays: " <<
forceTauDecay_ << std::endl;
309 if(
abs(pdg_id) != 15) {
310 edm::LogError(
"MuonReplacement") <<
"[ParticleReplacerAlgoParticleGun::tauola_forParticleGuns] "
311 <<
"Trying to decay something else than tau: pdg_id = " << pdg_id << std::endl;
354 else if (pdg_id == 15){
399 else return pol2_[2];
415 else return pol1_[2];
422 edm::LogError(
"MuonReplacement") <<
"[ParticleReplacerAlgoParticleGun::tauHelicity] "
423 <<
"tauHelicity undetermined, returning 0" << std::endl;
428 uint32_t randomNumber =
rand();
429 if(randomNumber%2 > 0.5)
return 1;
T getParameter(std::string const &) const
std::string generatorMode_
int forceTauMinusHelicity_
virtual void SetDecayRandomEngine(CLHEP::HepRandomEngine *decayRandomEngine)
CLHEP::HepRandomEngine * decayRandomEngine
static HepMC::IO_HEPEVT conv
double y() const
y coordinate
void correctTauMass(const reco::MuonCollection &muons, std::vector< HepMC::FourVector > &corrected)
gen::Pythia6Service pythia_
float randomPolarization()
virtual double vy() const
y coordinate of vertex position
void tauola_forParticleGun(int tau_idx, int pdg_id, const HepMC::FourVector &particle_momentum)
std::vector< Muon > MuonCollection
collection of Muon objects
std::string forceTauPolarization_
void call_pylist(int mode)
float tauHelicity(int pdg_id)
std::string particleOrigin_
ParticleReplacerParticleGun(const edm::ParameterSet &, bool)
double z() const
y coordinate
virtual CLHEP::HepRandomEngine & getEngine() const =0
Use this to get the random number engine, this is the only function most users should call...
virtual double vz() const
z coordinate of vertex position
double x() const
x coordinate
CLHEP::HepRandomEngine * decayRandomEngine
std::string forceTauDecay_
math::XYZTLorentzVector LorentzVector
Lorentz vector.
virtual ~ParticleReplacerParticleGun()
std::auto_ptr< HepMC::GenEvent > produce(const reco::MuonCollection &muons, const reco::Vertex *pvtx=0, const HepMC::GenEvent *genEvt=0)
static unsigned int const shift
virtual double vx() const
x coordinate of vertex position
int forceTauPlusHelicity_
gen::TauolaInterfaceBase * tauola_
T get(const Candidate &c)
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
void forceTauolaTauDecays()