9 #include "HepMC/IO_HEPEVT.h"
12 #define TXGIVE txgive_
14 void TXGIVE(
const char*,
int length);
19 #define TXGIVE_INIT txgive_init_
27 generatorMode_(pset.getParameter<std::string>(
"generatorMode")),
30 maxNumberOfAttempts_(pset.getUntrackedParameter<int>(
"maxNumberOfAttempts", 1000))
49 LogInfo(
"Replacer") <<
"won't do any transformation with the given mumu";
54 LogInfo(
"Replacer") <<
"will transform mumu into tautau";
59 LogInfo(
"Replacer") <<
"will transform mumu into taunu (as coming from a W boson)";
64 LogInfo(
"Replacer") <<
"Will transform mu-nu into tau-nu. No mass correction will be made.";
69 throw cms::Exception(
"ParticleReplacerClass") <<
"Unknown transformation mode!\n";
85 std::string minVisibleTransverseMomentumLine = pset.
getUntrackedParameter<std::string>(
"minVisibleTransverseMomentum",
"");
88 const char* startptr = minVisibleTransverseMomentumLine.c_str();
90 double d = strtod(startptr, &endptr);
91 if(*endptr ==
'\0' && endptr != startptr)
95 cuts[0].
pt_ = cuts[1].
pt_ = d;
97 minVisPtCuts_.push_back(std::vector<MinVisPtCut>(cuts, cuts+2));
104 pos = minVisibleTransverseMomentumLine.find(
';', prev);
105 if(
pos == std::string::npos)
pos = minVisibleTransverseMomentumLine.length();
107 std::string sub = minVisibleTransverseMomentumLine.substr(prev,
pos - prev);
108 std::vector<MinVisPtCut>
cuts;
109 const char* sub_c = sub.c_str();
110 while(*sub_c !=
'\0')
112 const char* sep = std::strchr(sub_c,
'_');
113 if(sep ==
NULL)
throw cms::Exception(
"Configuration") <<
"Minimum transverse parameter string must contain an underscore to separate type from pt threshold" << std::endl;
114 std::string
type(sub_c, sep);
125 else throw cms::Exception(
"Configuration") <<
"'" << type <<
"' is not a valid type. Allowed values are elec1,mu1,had1,tau1,elec2,mu2,had2,tau2" << std::endl;
128 cut.
pt_ = strtod(sep + 1, &endptr);
129 if(endptr == sep + 1)
throw cms::Exception(
"Configuration") <<
"No pt threshold given" << std::endl;
139 if(fileService_.isAvailable()) {
140 outTree = fileService_->make<TTree>(
"event_generation",
"This tree stores information about the event generation");
145 if(!rng.isAvailable()) {
147 <<
"The RandomNumberProducer module requires the RandomNumberGeneratorService\n"
148 "which appears to be absent. Please add that service to your configuration\n"
149 "or remove the modules that require it." << std::endl;
173 using namespace HepMC;
176 throw cms::Exception(
"Configuration") <<
"ParticleReplacerClass does NOT support using primary vertex as the origin for taus" << std::endl;
178 HepMC::GenEvent * evt=0;
180 GenVertex * zvtx =
new GenVertex();
186 std::vector<reco::Particle> particles;
193 LogError(
"Replacer") <<
"the decay mode Z->mumu requires exactly two muons, aborting processing";
194 return std::auto_ptr<HepMC::GenEvent>(0);
204 particles.push_back(tau1);
205 particles.push_back(tau2);
212 LogError(
"Replacer") <<
"the decay mode Z->tautau requires exactly two muons, aborting processing";
213 return std::auto_ptr<HepMC::GenEvent>(0);
224 particles.push_back(tau1);
225 particles.push_back(tau2);
232 LogError(
"Replacer") <<
"the decay mode Z->tautau requires exactly two muons, aborting processing";
233 return std::auto_ptr<HepMC::GenEvent>(0);
244 particles.push_back(tau1);
245 particles.push_back(tau2);
252 LogError(
"Replacer") <<
"transformation mode mu-nu ->tau-nu - wrong input";
253 return std::auto_ptr<HepMC::GenEvent>(0);
258 int targetParticlePdgIDNu_ = 16;
263 muon1.
px()*muon1.
px()+
264 muon1.
py()*muon1.
py()+
271 particles.push_back(tau1);
276 particles.push_back(nutau);
283 if (particles.size()==0)
285 LogError(
"Replacer") <<
"the creation of the new particles failed somehow";
286 return std::auto_ptr<HepMC::GenEvent>(0);
290 LogInfo(
"Replacer") << particles.size() <<
" particles found, continue processing";
297 evt =
new HepMC::GenEvent(*genEvt);
301 GenVertex * vtx=(*p);
303 if ( vtx->particles_out_size()<=0 || vtx->particles_in_size()<=0)
306 bool temp_muon1=
false,temp_muon2=
false,temp_z=
false;
307 for (GenVertex::particles_out_const_iterator it = vtx->particles_out_const_begin(); it!=vtx->particles_out_const_end(); it++)
309 if ((*it)->pdg_id()==13) temp_muon1=
true;
310 if ((*it)->pdg_id()==-13) temp_muon2=
true;
311 if ((*it)->pdg_id()==23) temp_z=
true;
314 int mother_pdg=(*vtx->particles_in_const_begin())->pdg_id();
316 if ((vtx->particles_out_size()==2 && vtx->particles_in_size()>0
320 || (vtx->particles_out_size()>2 && vtx->particles_in_size()>0
322 && temp_muon1 && temp_muon2 && temp_z ))
332 for (GenEvent::particle_iterator it=evt->particles_begin();it!=evt->particles_end();it++)
333 (*it)->set_status(0);
335 for (std::vector<reco::Particle>::const_iterator it=particles.begin();it!=particles.end();it++)
337 zvtx->add_particle_out(
new HepMC::GenParticle((FourVector)it->p4(), it->pdgId(), 1, Flow(), Polarization(0,0)));
344 for (std::vector<reco::Particle>::const_iterator it=particles.begin();it!=particles.end();it++)
345 mother_particle_p4+=it->p4();
349 GenVertex* startVtx =
new GenVertex(FourVector(production_point.x()*10,production_point.y()*10,production_point.z()*10,0));
350 startVtx->add_particle_in(
new GenParticle( FourVector(0,0,7000,7000), 2212, 3 ) );
351 startVtx->add_particle_in(
new GenParticle( FourVector(0,0,-7000,7000), 2212, 3 ) );
353 GenVertex * decayvtx =
new GenVertex(FourVector(production_point.x()*10,production_point.y()*10,production_point.z()*10,0));
357 int muPDG = particles.begin()->pdgId();
359 mother_particle->set_pdg_id(
id);
362 startVtx->add_particle_out(mother_particle);
363 decayvtx->add_particle_in(mother_particle);
364 evt =
new HepMC::GenEvent();
365 for (std::vector<reco::Particle>::const_iterator it=particles.begin();it!=particles.end();it++)
368 decayvtx->add_particle_out(
new HepMC::GenParticle((FourVector)it->p4(), it->pdgId(), 1, Flow(), Polarization(0,0)));
371 evt->add_vertex(startVtx);
372 evt->add_vertex(decayvtx);
376 HepMC::GenEvent * retevt = 0;
377 HepMC::GenEvent * tempevt = 0;
382 unsigned int cntVisPt_all = 0;
383 unsigned int cntVisPt_pass = 0;
385 HepMC::IO_HEPEVT
conv;
391 LogError(
"Replacer") <<
"Pythia is currently not supported!";
392 return std::auto_ptr<HepMC::GenEvent>(evt);
397 conv.write_event(evt);
413 eventWeight = (double)cntVisPt_pass / (
double)cntVisPt_all;
414 std::cout <<
" " << cntVisPt_pass <<
"\t" << cntVisPt_all <<
"\n";
417 LogError(
"Replacer") <<
"failed to create an event which satisfies the minimum visible transverse momentum cuts ";
421 return std::auto_ptr<HepMC::GenEvent>(0);
429 for (GenEvent::particle_iterator it=retevt->particles_begin();it!=retevt->particles_end();it++)
431 if ((*it)->end_vertex())
432 (*it)->set_status(2);
434 (*it)->set_status(1);
438 std::auto_ptr<HepMC::GenEvent>
ret(retevt);
465 using namespace HepMC;
471 std::vector<double> mus;
472 std::vector<double> elecs;
473 std::vector<double> hads;
474 std::vector<double> taus;
476 for (GenEvent::particle_iterator it=evt->particles_begin();it!=evt->particles_end();it++)
478 if (
abs((*it)->pdg_id())==15 && (*it)->end_vertex())
480 FourVector vis_mom();
482 std::queue<const GenParticle *> decaying_particles;
483 decaying_particles.push(*it);
485 enum { ELEC, MU, HAD }
type = HAD;
486 while(!decaying_particles.empty() && (++t < 30))
489 decaying_particles.pop();
491 if (!front->end_vertex())
494 if (pdgId>10 && pdgId!=12 && pdgId!=14 && pdgId!=16)
497 if(pdgId == 11)
type = ELEC;
498 if(pdgId == 13)
type = MU;
502 GenVertex * temp_vert = front->end_vertex();
503 for (GenVertex::particles_out_const_iterator it2=temp_vert->particles_out_const_begin();it2!=temp_vert->particles_out_const_end();it2++)
504 decaying_particles.push((*it2));
508 double vis_pt = visible_momentum.pt();
509 taus.push_back(vis_pt);
510 if(
type == MU) mus.push_back(vis_pt);
511 if(
type == ELEC) elecs.push_back(vis_pt);
512 if(
type == HAD) hads.push_back(vis_pt);
516 std::sort(taus.begin(), taus.end(), std::greater<double>());
517 std::sort(elecs.begin(), elecs.end(), std::greater<double>());
518 std::sort(mus.begin(), mus.end(), std::greater<double>());
519 std::sort(hads.begin(), hads.end(), std::greater<double>());
523 std::vector<MinVisPtCut>::const_iterator iter2;
524 for(iter2 = iter->begin(); iter2 != iter->end(); ++iter2)
533 default: assert(
false);
break;
537 if(iter2->index_ >= collection->size() || (*collection)[iter2->index_] < iter2->pt_)
542 if(iter2 == iter->end())
546 LogInfo(
"Replacer") <<
"refusing the event as the sum of the visible transverse momenta is too small\n";
552 using namespace HepMC;
555 using namespace reco;
557 stack<HepMC::GenParticle *> deleteParticle;
559 stack<GenVertex *> deleteVertex;
560 stack<GenVertex *> queueVertex;
562 if (vtx->particles_out_size()>0)
564 for (GenVertex::particles_out_const_iterator it=vtx->particles_out_const_begin();it!=vtx->particles_out_const_end();it++)
566 deleteParticle.push(*it);
567 if ((*it)->end_vertex())
568 queueVertex.push((*it)->end_vertex());
572 while (!queueVertex.empty())
574 GenVertex * temp_vtx=queueVertex.top();
575 if (temp_vtx->particles_out_size()>0)
577 for (GenVertex::particles_out_const_iterator it=temp_vtx->particles_out_const_begin();it!=temp_vtx->particles_out_const_end();it++)
579 if ((*it)->end_vertex())
580 queueVertex.push((*it)->end_vertex());
584 deleteVertex.push(queueVertex.top());
588 while (!deleteVertex.empty())
590 evt->remove_vertex(deleteVertex.top());
594 while (!deleteParticle.empty())
596 delete vtx->remove_particle(deleteParticle.top());
597 deleteParticle.pop();
600 while (!deleteVertex.empty())
602 while (!queueVertex.empty())
610 using namespace HepMC;
615 while (!(*it)->suggest_barcode(-1*(++max_barc)))
619 for (GenEvent::particle_iterator it=evt->particles_begin();it!=evt->particles_end();it++)
620 while (!(*it)->suggest_barcode(++max_barc))
628 using namespace reco;
635 ROOT::Math::Boost booster(z_momentum.BoostToCM());
636 ROOT::Math::Boost invbooster(booster.Inverse());
645 double muonxb_mom2 = muon1b.x()*muon1b.x() + muon1b.y()*muon1b.y() + muon1b.z() * muon1b.z();
646 double tauxb_mom2 = 0.25 * Zb.t() * Zb.t() - tau_mass2;
648 float scaling1 =
sqrt(tauxb_mom2 / muonxb_mom2);
649 float scaling2 = scaling1;
651 float tauEnergy= Zb.t() / 2.;
653 if (tauEnergy*tauEnergy<tau_mass2)
668 assert(
std::abs((muon1_momentum+muon2_momentum).
mass()-(tau1_mom+tau2_mom).
mass())/(muon1_momentum+muon2_momentum).
mass()<0.001);
669 assert(
std::abs((muon1_momentum+muon2_momentum).pt()-(tau1_mom+tau2_mom).pt())/(muon1_momentum+muon2_momentum).pt()<0.001);
671 muon1->
setP4(tau1_mom);
672 muon2->
setP4(tau2_mom);
686 using namespace reco;
693 ROOT::Math::Boost booster(z_momentum.BoostToCM());
694 ROOT::Math::Boost invbooster(booster.Inverse());
698 const double breitWignerWidth_Z = 2.4952;
699 const double breitWignerWidth_W = 2.141;
700 const double knownMass_W = 80.398;
701 const double knownMass_Z = 91.1876;
703 double Wb_mass = ( Zb.mass() - knownMass_Z ) * ( breitWignerWidth_W / breitWignerWidth_Z ) + knownMass_W;
704 std::cout <<
"Wb_mass: " << Wb_mass <<
"\n";
711 double muonxb_mom2 = muon1b.x()*muon1b.x() + muon1b.y()*muon1b.y() + muon1b.z() * muon1b.z();
712 double tauxb_mom2 = 0.25 * Zb.t() * Zb.t() - tau_mass2;
714 float scaling1 =
sqrt(tauxb_mom2 / muonxb_mom2) * Wb_mass/Zb.mass();
715 float scaling2 = scaling1;
717 float tauEnergy= Zb.t() / 2.;
719 if (tauEnergy*tauEnergy<tau_mass2)
725 std::cout <<
"muon1b_momentum: " << muon1b <<
"\n";
726 std::cout <<
"muon2b_momentum: " << muon2b <<
"\n";
728 std::cout <<
"tau1b_momentum: " << tau1b_mom <<
"\n";
729 std::cout <<
"tau2b_momentum: " << tau2b_mom <<
"\n";
731 std::cout <<
"zb_momentum: " << Zb <<
"\n";
732 std::cout <<
"wb_momentum: " << (tau1b_mom+tau2b_mom) <<
"\n";
739 assert((muon1b.phi()-tau1b_mom.phi())/muon1b.phi()<0.001);
740 assert((muon2b.phi()-tau2b_mom.phi())/muon2b.phi()<0.001);
741 assert((muon1b.theta()-tau1b_mom.theta())/muon1b.theta()<0.001);
742 assert((muon2b.theta()-tau2b_mom.theta())/muon2b.theta()<0.001);
756 part1->
setP4(tau1_mom);
757 part2->
setP4(tau2_mom);
T getParameter(std::string const &) const
double targetParticleMass_
T getUntrackedParameter(std::string const &, T const &) const
virtual int pdgId() const
PDG identifier.
void repairBarcodes(HepMC::GenEvent *evt)
virtual void SetDecayRandomEngine(CLHEP::HepRandomEngine *decayRandomEngine)
void setP4(const LorentzVector &p4)
set 4-momentum
CLHEP::HepRandomEngine * decayRandomEngine
virtual void init(const edm::EventSetup &)
static HepMC::IO_HEPEVT conv
virtual const Point & vertex() const
vertex position
enum ParticleReplacerClass::MinVisPtCut::@449 type_
ParticleReplacerClass(const edm::ParameterSet &, bool)
bool testEvent(HepMC::GenEvent *evt)
unsigned int replacementMode_
replace mode:
const LorentzVector & p4() const
four-momentum Lorentz vector
virtual void statistics()
int pdgId() const
PDG identifier.
std::vector< Muon > MuonCollection
collection of Muon objects
PtEtaPhiMLorentzVectorD PtEtaPhiMLorentzVector
Lorentz vector with cartesian internal representation.
gen::TauolaInterfaceBase * tauola_
virtual int charge() const
electric charge
math::XYZPoint Point
point in the space
std::string generatorMode_
virtual double px() const
x coordinate of momentum vector
virtual HepMC::GenEvent * decay(HepMC::GenEvent *evt)
virtual double pz() const
z coordinate of momentum vector
unsigned int transformationMode_
math::XYZTLorentzVector LorentzVector
Lorentz vector.
void cleanEvent(HepMC::GenEvent *evt, HepMC::GenVertex *vtx)
void transformMuMu2TauNu(reco::Particle *muon1, reco::Particle *muon2)
transform a muon pair into tau nu (as coming from a W boson)
virtual void beginRun(edm::Run &iRun, const edm::EventSetup &iSetup)
VertexRefVector::iterator vertex_iterator
iterator over a vector of references to Vertex objects in the same collection
void transformMuMu2TauTau(reco::Particle *muon1, reco::Particle *muon2)
transform a muon pair into a tau pair
virtual std::auto_ptr< HepMC::GenEvent > produce(const reco::MuonCollection &muons, const reco::Vertex *pvtx=0, const HepMC::GenEvent *genEvt=0)
virtual const LorentzVector & p4() const
four-momentum Lorentz vector
math::XYZTLorentzVector LorentzVector
Lorentz vector.
void setStatus(int status)
set status word
virtual double py() const
y coordinate of momentum vector
std::vector< std::vector< MinVisPtCut > > minVisPtCuts_
T get(const Candidate &c)