File indexing completed on 2025-04-04 01:26:46
0001 #include <ostream>
0002
0003 #include "IOMC/ParticleGuns/interface/FlatRandomMultiParticlePGunProducer.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/MessageLogger/interface/MessageLogger.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "FWCore/Utilities/interface/EDMException.h"
0014 #include "CLHEP/Random/RandFlat.h"
0015
0016 using namespace edm;
0017
0018 FlatRandomMultiParticlePGunProducer::FlatRandomMultiParticlePGunProducer(const ParameterSet& pset)
0019 : BaseFlatGunProducer(pset) {
0020 ParameterSet pgunParams = pset.getParameter<ParameterSet>("PGunParameters");
0021 fProbParticle_ = pgunParams.getParameter<std::vector<double> >("ProbParts");
0022 fMinP_ = pgunParams.getParameter<double>("MinP");
0023 fMaxP_ = pgunParams.getParameter<double>("MaxP");
0024 if (fProbParticle_.size() != fPartIDs.size())
0025 throw cms::Exception("Configuration") << "Not all probabilities given for all particle types "
0026 << fProbParticle_.size() << ":" << fPartIDs.size() << " need them to match\n";
0027
0028 produces<HepMCProduct>("unsmeared");
0029 produces<GenEventInfoProduct>();
0030
0031 edm::LogVerbatim("ParticleGun") << "FlatRandomMultiParticlePGun is initialzed for " << fPartIDs.size()
0032 << " particles in momentum range " << fMinP_ << ":" << fMaxP_;
0033 for (unsigned int k = 0; k < fPartIDs.size(); ++k)
0034 edm::LogVerbatim("ParticleGun") << " [" << k << "] " << fPartIDs[k] << ":" << fProbParticle_[k];
0035
0036 for (unsigned int k = 1; k < fProbParticle_.size(); ++k)
0037 fProbParticle_[k] += fProbParticle_[k - 1];
0038 for (unsigned int k = 0; k < fProbParticle_.size(); ++k) {
0039 fProbParticle_[k] /= fProbParticle_[fProbParticle_.size() - 1];
0040 edm::LogVerbatim("ParticleGun") << "Corrected probaility for " << fPartIDs[k] << ":" << fProbParticle_[k];
0041 }
0042 }
0043
0044 FlatRandomMultiParticlePGunProducer::~FlatRandomMultiParticlePGunProducer() {}
0045
0046 void FlatRandomMultiParticlePGunProducer::produce(edm::Event& e, const edm::EventSetup& es) {
0047 edm::Service<edm::RandomNumberGenerator> rng;
0048 CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
0049
0050 if (fVerbosity > 0)
0051 edm::LogVerbatim("ParticleGun") << "FlatRandomMultiParticlePGunProducer: Begin New Event Generation";
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(0), PartID(fPartIDs[0]);
0069 double r1 = CLHEP::RandFlat::shoot(engine, 0., 1.);
0070 for (unsigned int ip = 0; ip < fPartIDs.size(); ip++) {
0071 if (r1 <= fProbParticle_[ip])
0072 break;
0073 PartID = fPartIDs[ip];
0074 }
0075 if (fVerbosity > 0)
0076 edm::LogVerbatim("ParticleGun") << "Random " << r1 << " PartID " << PartID;
0077
0078 double mom = CLHEP::RandFlat::shoot(engine, fMinP_, fMaxP_);
0079 double eta = CLHEP::RandFlat::shoot(engine, fMinEta, fMaxEta);
0080 double phi = CLHEP::RandFlat::shoot(engine, fMinPhi, fMaxPhi);
0081 const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID)));
0082 double mass = PData->mass().value();
0083 double energy = sqrt(mom * mom + mass * mass);
0084 double theta = 2. * atan(exp(-eta));
0085 double px = mom * sin(theta) * cos(phi);
0086 double py = mom * sin(theta) * sin(phi);
0087 double pz = mom * cos(theta);
0088
0089 HepMC::FourVector p(px, py, pz, energy);
0090 HepMC::GenParticle* Part = new HepMC::GenParticle(p, PartID, 1);
0091 barcode++;
0092 Part->suggest_barcode(barcode);
0093 Vtx->add_particle_out(Part);
0094
0095 if (fAddAntiParticle) {
0096 HepMC::FourVector ap(-px, -py, -pz, energy);
0097 int APartID = (PartID == 22 || PartID == 23) ? PartID : -PartID;
0098 HepMC::GenParticle* APart = new HepMC::GenParticle(ap, APartID, 1);
0099 barcode++;
0100 APart->suggest_barcode(barcode);
0101 Vtx->add_particle_out(APart);
0102 }
0103
0104 fEvt->add_vertex(Vtx);
0105 fEvt->set_event_number(e.id().event());
0106 fEvt->set_signal_process_id(20);
0107
0108 if (fVerbosity > 1)
0109 fEvt->print();
0110
0111 std::unique_ptr<HepMCProduct> BProduct(new HepMCProduct());
0112 BProduct->addHepMCData(fEvt);
0113 e.put(std::move(BProduct), "unsmeared");
0114
0115 std::unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
0116 e.put(std::move(genEventInfo));
0117
0118 if (fVerbosity > 0)
0119 edm::LogVerbatim("ParticleGun") << " FlatRandomMultiParticlePGunProducer : Event Generation Done";
0120 }