6 #include "HepMC/PythiaWrapper.h"
7 #include "HepMC/IO_HEPEVT.h"
13 particleOrigin_(iConfig.getParameter<std::
string>(
"particleOrigin")),
14 forceTauPolarization_(iConfig.getParameter<std::
string>(
"forceTauPolarization")),
15 forceTauDecay_(iConfig.getParameter<std::
string>(
"forceTauDecay")),
16 generatorMode_(iConfig.getParameter<std::
string>(
"generatorMode")),
17 gunParticle_(iConfig.getParameter<int>(
"gunParticle")),
18 forceTauPlusHelicity_(iConfig.getParameter<int>(
"forceTauPlusHelicity")),
19 forceTauMinusHelicity_(iConfig.getParameter<int>(
"forceTauMinusHelicity")),
25 edm::LogInfo(
"MuonReplacement") <<
"[ParticleReplacer::ParticleReplacer] "
28 edm::LogInfo(
"MuonReplacement") <<
"[ParticleReplacer::ParticleReplacer] "
31 edm::LogInfo(
"MuonReplacement") <<
"[ParticleReplacer::ParticleReplacer] "
34 edm::LogInfo(
"MuonReplacement") <<
"[ParticleReplacer::ParticleReplacer] "
37 std::memset(
pol1_, 0, 4*
sizeof(
float));
38 std::memset(
pol2_, 0, 4*
sizeof(
float));
41 throw cms::Exception(
"Configuration") <<
"Generator mode other than Tauola is not supported" << std::endl;
43 throw cms::Exception(
"UnimplementedFeature") <<
"ParticleReplacerParticleGun is not usable yet." << std::endl;
70 throw cms::Exception(
"UnimplementedFeature") <<
"ParticleReplacerParticleGun does NOT support merging at HepMC level" << std::endl;
72 std::auto_ptr<HepMC::GenEvent> evt(0);
73 std::vector<HepMC::FourVector> muons_corrected;
74 muons_corrected.reserve(muons.size());
79 for(
unsigned int i=0;
i<muons_corrected.size(); ++
i) {
80 HepMC::FourVector&
muon = muons_corrected[
i];
87 std::cout <<
" ///////////////////// After py1ent, before pyhepc /////////////////////" << std::endl;
95 std::cout <<
" ///////////////////// After pyhepc, before vertex shift /////////////////////" << std::endl;
96 HepMC::HEPEVT_Wrapper::print_hepevt();
99 int nparticles = HepMC::HEPEVT_Wrapper::number_entries();
100 HepMC::ThreeVector
shift(0,0,0);
103 throw cms::Exception(
"LogicError") <<
"Particle origin set to primaryVertex, but pvtx is null!" << std::endl;
105 shift.setX(pvtx->
x()*10);
106 shift.setY(pvtx->
y()*10);
107 shift.setZ(pvtx->
z()*10);
109 for(
int i=1;
i <= nparticles; ++
i) {
111 throw cms::Exception(
"LogicError") <<
"Particle in HEPEVT is " << HepMC::HEPEVT_Wrapper::id(
i)
113 <<
" for index " <<
i << std::endl;
118 shift.setX(muon.
vx()*10);
119 shift.setY(muon.
vy()*10);
120 shift.setZ(muon.
vz()*10);
123 HepMC::HEPEVT_Wrapper::set_position(
i,
131 std::cout <<
" ///////////////////// After vertex shift, before pyhepc/tauola /////////////////////" << std::endl;
132 HepMC::HEPEVT_Wrapper::print_hepevt();
144 for(
unsigned int i=0;
i<muons_corrected.size(); ++
i) {
149 std::cout <<
" ///////////////////// After tauola, before pyhepc /////////////////////" << std::endl;
150 HepMC::HEPEVT_Wrapper::print_hepevt();
154 std::cout <<
" ///////////////////// After pyhepc, before vertex fix /////////////////////" << std::endl;
169 std::cout <<
" ///////////////////// After vertex fix, before pyexec /////////////////////" << std::endl;
176 std::cout <<
" ///////////////////// After pyhepc, before pyexec /////////////////////" << std::endl;
184 std::cout <<
" ///////////////////// After pyexec, before pyhepc /////////////////////" << std::endl;
190 HepMC::IO_HEPEVT
conv;
191 evt = std::auto_ptr<HepMC::GenEvent>(
new HepMC::GenEvent(*conv.read_next_event()));
196 std::cout << std::endl <<
"Vertices: " << std::endl;
197 for(HepMC::GenEvent::vertex_const_iterator iter = evt->vertices_begin(); iter != evt->vertices_end(); ++iter) {
198 std::cout <<
"Vertex " << (*iter)->id() <<
", barcode " << (*iter)->barcode() << std::endl;
200 HepMC::ThreeVector
point = (*iter)->point3d();
201 std::cout <<
" Point (" << point.x() <<
", " << point.y() <<
", " << point.z() <<
")" << std::endl;
204 for(HepMC::GenVertex::particles_in_const_iterator
pi = (*iter)->particles_in_const_begin();
pi != (*iter)->particles_in_const_end(); ++
pi) {
209 for(HepMC::GenVertex::particles_out_const_iterator
pi = (*iter)->particles_out_const_begin();
pi != (*iter)->particles_out_const_end(); ++
pi) {
222 for(reco::MuonCollection::const_iterator iMuon = muons.begin(); iMuon != muons.end(); ++iMuon) {
225 double E =
sqrt(
tauMass*
tauMass + vec.x()*vec.x() + vec.y()*vec.y() + vec.z()*vec.z());
227 corrected.push_back(HepMC::FourVector(vec.x(), vec.y(), vec.z(), E));
232 for(reco::MuonCollection::const_iterator iMuon = muons.begin(); iMuon != muons.end(); ++iMuon) {
233 corrected.push_back(iMuon->p4());
283 edm::LogError(
"MuonReplacement") <<
"[ParticleReplacerAlgoParticleGun::forceTauoladecays] "
284 <<
"Unknown value for forcing tau decays: " <<
forceTauDecay_ << std::endl;
288 if(
abs(pdg_id) != 15) {
289 edm::LogError(
"MuonReplacement") <<
"[ParticleReplacerAlgoParticleGun::tauola_forParticleGuns] "
290 <<
"Trying to decay something else than tau: pdg_id = " << pdg_id << std::endl;
333 else if (pdg_id == 15){
378 else return pol2_[2];
394 else return pol1_[2];
401 edm::LogError(
"MuonReplacement") <<
"[ParticleReplacerAlgoParticleGun::tauHelicity] "
402 <<
"tauHelicity undetermined, returning 0" << std::endl;
407 uint32_t randomNumber =
rand();
408 if(randomNumber%2 > 0.5)
return 1;
T getParameter(std::string const &) const
std::string generatorMode_
int forceTauMinusHelicity_
static HepMC::IO_HEPEVT conv
double y() const
y coordinate
void setPSet(const edm::ParameterSet &)
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_
gen::TauolaInterface * tauola_
void call_pylist(int mode)
float tauHelicity(int pdg_id)
std::string particleOrigin_
ParticleReplacerParticleGun(const edm::ParameterSet &, bool)
double z() const
y coordinate
virtual double vz() const
z coordinate of vertex position
double x() const
x coordinate
std::string forceTauDecay_
math::XYZTLorentzVector LorentzVector
Lorentz vector.
virtual ~ParticleReplacerParticleGun()
void tauola_(int *, int *)
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_
*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()