ParticleGuns/src/RandomXiThetaGunProducer.cc

0001 /****************************************************************************
0002  * Authors:
0003  *   Jan Kašpar
0004  ****************************************************************************/
0005
0006 #include "IOMC/ParticleGuns/interface/RandomXiThetaGunProducer.h"
0007
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include "FWCore/ServiceRegistry/interface/Service.h"
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0015
0016 using namespace edm;
0017 using namespace std;
0018
0019 //----------------------------------------------------------------------------------------------------
0020
0021 RandomXiThetaGunProducer::RandomXiThetaGunProducer(const edm::ParameterSet &iConfig)
0022     : tableToken_(esConsumes()),
0023       verbosity_(iConfig.getUntrackedParameter<unsigned int>("verbosity", 0)),
0024       particleId_(iConfig.getParameter<unsigned int>("particleId")),
0025       energy_(iConfig.getParameter<double>("energy")),
0026       xi_min_(iConfig.getParameter<double>("xi_min")),
0027       xi_max_(iConfig.getParameter<double>("xi_max")),
0028       theta_x_mean_(iConfig.getParameter<double>("theta_x_mean")),
0029       theta_x_sigma_(iConfig.getParameter<double>("theta_x_sigma")),
0030       theta_y_mean_(iConfig.getParameter<double>("theta_y_mean")),
0031       theta_y_sigma_(iConfig.getParameter<double>("theta_y_sigma")),
0032       nParticlesSector45_(iConfig.getParameter<unsigned int>("nParticlesSector45")),
0033       nParticlesSector56_(iConfig.getParameter<unsigned int>("nParticlesSector56")),
0034       engine_(nullptr) {
0035   produces<HepMCProduct>("unsmeared");
0036 }
0037
0038 //----------------------------------------------------------------------------------------------------
0039
0040 void RandomXiThetaGunProducer::produce(edm::Event &e, const edm::EventSetup &es) {
0041   // get conditions
0042   edm::Service<edm::RandomNumberGenerator> rng;
0043   engine_ = &rng->getEngine(e.streamID());
0044
0045   auto const &pdgTable = es.getData(tableToken_);
0046
0047   // prepare HepMC event
0048   HepMC::GenEvent *fEvt = new HepMC::GenEvent();
0049   fEvt->set_event_number(e.id().event());
0050
0051   HepMC::GenVertex *vtx = new HepMC::GenVertex(HepMC::FourVector(0., 0., 0., 0.));
0052   fEvt->add_vertex(vtx);
0053
0054   const HepPDT::ParticleData *pData = pdgTable.particle(HepPDT::ParticleID(particleId_));
0055   double mass = pData->mass().value();
0056
0057   // generate particles
0058   unsigned int barcode = 0;
0059
0060   for (unsigned int i = 0; i < nParticlesSector45_; ++i)
0061     generateParticle(+1., mass, ++barcode, vtx);
0062
0063   for (unsigned int i = 0; i < nParticlesSector56_; ++i)
0064     generateParticle(-1., mass, ++barcode, vtx);
0065
0066   // save output
0067   std::unique_ptr<HepMCProduct> output(new HepMCProduct());
0068   output->addHepMCData(fEvt);
0069   e.put(std::move(output), "unsmeared");
0070 }
0071
0072 //----------------------------------------------------------------------------------------------------
0073
0074 void RandomXiThetaGunProducer::generateParticle(double z_sign,
0075                                                 double mass,
0076                                                 unsigned int barcode,
0077                                                 HepMC::GenVertex *vtx) const {
0078   const double xi = CLHEP::RandFlat::shoot(engine_, xi_min_, xi_max_);
0079   const double theta_x = CLHEP::RandGauss::shoot(engine_, theta_x_mean_, theta_x_sigma_);
0080   const double theta_y = CLHEP::RandGauss::shoot(engine_, theta_y_mean_, theta_y_sigma_);
0081
0082   if (verbosity_)
0083     LogInfo("RandomXiThetaGunProducer") << "xi = " << xi << ", theta_x = " << theta_x << ", theta_y" << theta_y
0084                                         << ", z_sign = " << z_sign;
0085
0086   const double cos_theta = sqrt(1. - theta_x * theta_x - theta_y * theta_y);
0087
0088   const double p_nom = sqrt(energy_ * energy_ - mass * mass);
0089   const double p = p_nom * (1. - xi);
0090   const double e = sqrt(p * p + mass * mass);
0091
0092   HepMC::FourVector momentum(p * theta_x, p * theta_y, z_sign * p * cos_theta, e);
0093
0094   HepMC::GenParticle *particle = new HepMC::GenParticle(momentum, particleId_, 1);
0095   particle->suggest_barcode(barcode);
0096   vtx->add_particle_out(particle);
0097 }