File indexing completed on 2024-04-06 12:19:03
0001
0002
0003
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
0042 edm::Service<edm::RandomNumberGenerator> rng;
0043 engine_ = &rng->getEngine(e.streamID());
0044
0045 auto const &pdgTable = es.getData(tableToken_);
0046
0047
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
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
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 }