File indexing completed on 2025-04-04 01:26:46
0001 #include <ostream>
0002
0003 #include "IOMC/ParticleGuns/interface/GaussRandomPThetaGunProducer.h"
0004
0005 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0006 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0007
0008 #include "FWCore/AbstractServices/interface/RandomNumberGenerator.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013
0014 #include "CLHEP/Random/RandGaussQ.h"
0015 #include "CLHEP/Random/RandFlat.h"
0016
0017
0018
0019 namespace CLHEP {
0020 class HepRandomEngine;
0021 }
0022
0023 using namespace edm;
0024
0025 GaussRandomPThetaGunProducer::GaussRandomPThetaGunProducer(const edm::ParameterSet& pset)
0026 : FlatBaseThetaGunProducer(pset) {
0027 edm::ParameterSet defpset;
0028 edm::ParameterSet pgun_params = pset.getParameter<edm::ParameterSet>("PGunParameters");
0029
0030
0031
0032 fMeanP = pgun_params.getParameter<double>("MeanP");
0033 fSigmaP = pgun_params.getParameter<double>("SigmaP");
0034 fMeanTheta = 0.5 * (fMinTheta + fMaxTheta);
0035 fSigmaTheta = 0.5 * (fMaxTheta - fMinTheta);
0036
0037 produces<HepMCProduct>("unsmeared");
0038 produces<GenEventInfoProduct>();
0039 }
0040
0041 GaussRandomPThetaGunProducer::~GaussRandomPThetaGunProducer() {}
0042
0043 void GaussRandomPThetaGunProducer::produce(edm::Event& e, const edm::EventSetup& es) {
0044 #ifdef DebugLog
0045 if (fVerbosity >= 0)
0046 std::cout << "GaussRandomPThetaGunProducer : Begin New Event Generation\n";
0047 #endif
0048 edm::Service<edm::RandomNumberGenerator> rng;
0049 CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
0050
0051
0052
0053
0054
0055
0056
0057 fEvt = new HepMC::GenEvent();
0058
0059
0060
0061
0062
0063
0064 HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(0., 0., 0.));
0065
0066
0067
0068 int barcode = 1;
0069 for (unsigned int ip = 0; ip < fPartIDs.size(); ip++) {
0070 double mom = CLHEP::RandGaussQ::shoot(engine, fMeanP, fSigmaP);
0071 double theta = CLHEP::RandGaussQ::shoot(engine, fMeanTheta, fSigmaTheta);
0072 double phi = CLHEP::RandFlat::shoot(engine, fMinPhi, fMaxPhi);
0073 int PartID = fPartIDs[ip];
0074 const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID)));
0075 double mass = PData->mass().value();
0076 double energy = sqrt(mom * mom + mass * mass);
0077 double px = mom * sin(theta) * cos(phi);
0078 double py = mom * sin(theta) * sin(phi);
0079 double pz = mom * cos(theta);
0080
0081 HepMC::FourVector p(px, py, pz, energy);
0082 HepMC::GenParticle* Part = new HepMC::GenParticle(p, PartID, 1);
0083 Part->suggest_barcode(barcode);
0084 barcode++;
0085 Vtx->add_particle_out(Part);
0086
0087 if (fAddAntiParticle) {
0088 HepMC::FourVector ap(-px, -py, -pz, energy);
0089 int APartID = -PartID;
0090 if (PartID == 22 || PartID == 23) {
0091 APartID = PartID;
0092 }
0093 HepMC::GenParticle* APart = new HepMC::GenParticle(ap, APartID, 1);
0094 APart->suggest_barcode(barcode);
0095 barcode++;
0096 Vtx->add_particle_out(APart);
0097 }
0098 }
0099 fEvt->add_vertex(Vtx);
0100 fEvt->set_event_number(e.id().event());
0101 fEvt->set_signal_process_id(20);
0102
0103 #ifdef DebugLog
0104 if (fVerbosity >= 0)
0105 fEvt->print();
0106 #endif
0107
0108 std::unique_ptr<HepMCProduct> BProduct(new HepMCProduct());
0109 BProduct->addHepMCData(fEvt);
0110 e.put(std::move(BProduct), "unsmeared");
0111
0112 std::unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
0113 e.put(std::move(genEventInfo));
0114
0115 #ifdef DebugLog
0116 if (fVerbosity >= 0)
0117 std::cout << "GaussRandomPThetaGunProducer : Event Generation Done\n";
0118 #endif
0119 }