File indexing completed on 2025-04-04 01:26:45
0001 #include <ostream>
0002
0003 #include "IOMC/ParticleGuns/interface/FileRandomMultiParticlePGunProducer.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 const unsigned int np = 6;
0019 const unsigned int kfactor = 100;
0020
0021 FileRandomMultiParticlePGunProducer::FileRandomMultiParticlePGunProducer(const ParameterSet& pset)
0022 : BaseFlatGunProducer(pset) {
0023 ParameterSet pgunParams = pset.getParameter<ParameterSet>("PGunParameters");
0024 fMinP_ = pgunParams.getParameter<double>("MinP");
0025 fMaxP_ = pgunParams.getParameter<double>("MaxP");
0026 edm::FileInPath fp = pgunParams.getParameter<edm::FileInPath>("FileName");
0027 const std::string& file = fp.fullPath();
0028
0029 produces<HepMCProduct>("unsmeared");
0030 produces<GenEventInfoProduct>();
0031 edm::LogVerbatim("ParticleGun") << "FileRandomMultiParticlePGun is initialzed with i/p file " << file
0032 << " and use momentum range " << fMinP_ << ":" << fMaxP_;
0033
0034 if (fPartIDs.size() != np)
0035 throw cms::Exception("ParticleGun") << "Invalid list of partices: " << fPartIDs.size() << " should be " << np
0036 << "\n";
0037
0038 std::ifstream is(file.c_str(), std::ios::in);
0039 if (!is) {
0040 throw cms::Exception("Configuration") << "Cannot find the file " << file << "\n";
0041 } else {
0042 double xl, xh;
0043 is >> fPBin_ >> xl >> xh >> fEtaBin_ >> fEtaMin_ >> fEtaBinWidth_;
0044 fP_.emplace_back(xl);
0045 edm::LogVerbatim("ParticleGun") << "FileRandomMultiParticlePGun: p: " << fPBin_ << ":" << xl << ":" << xh
0046 << " Eta: " << fEtaBin_ << ":" << fEtaMin_ << ":" << fEtaBinWidth_;
0047 for (int ip = 0; ip < fPBin_; ++ip) {
0048 for (int ie = 0; ie < fEtaBin_; ++ie) {
0049 double totprob(0);
0050 std::vector<double> prob(np, 0);
0051 int je;
0052 is >> xl >> xh >> je >> prob[0] >> prob[1] >> prob[2] >> prob[3] >> prob[4] >> prob[5];
0053 if (ie == 0)
0054 fP_.emplace_back(xh);
0055 for (unsigned int k = 0; k < np; ++k) {
0056 totprob += prob[k];
0057 if (k > 0)
0058 prob[k] += prob[k - 1];
0059 }
0060 for (unsigned int k = 0; k < np; ++k)
0061 prob[k] /= totprob;
0062 int indx = (ip + 1) * kfactor + ie;
0063 fProbParticle_[indx] = prob;
0064 if (fVerbosity > 0)
0065 edm::LogVerbatim("ParticleGun")
0066 << "FileRandomMultiParticlePGun [" << ip << "," << ie << ", " << indx << "] Probability " << prob[0]
0067 << ", " << prob[1] << ", " << prob[2] << ", " << prob[3] << ", " << prob[4] << ", " << prob[5];
0068 }
0069 }
0070 is.close();
0071 }
0072 }
0073
0074 FileRandomMultiParticlePGunProducer::~FileRandomMultiParticlePGunProducer() {}
0075
0076 void FileRandomMultiParticlePGunProducer::produce(edm::Event& e, const edm::EventSetup& es) {
0077 edm::Service<edm::RandomNumberGenerator> rng;
0078 CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
0079
0080 if (fVerbosity > 0)
0081 edm::LogVerbatim("ParticleGun") << "FileRandomMultiParticlePGunProducer: Begin New Event Generation";
0082
0083
0084
0085
0086
0087 fEvt = new HepMC::GenEvent();
0088
0089
0090
0091
0092
0093
0094 HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(0., 0., 0.));
0095
0096
0097 double mom = CLHEP::RandFlat::shoot(engine, fMinP_, fMaxP_);
0098 double eta = CLHEP::RandFlat::shoot(engine, fMinEta, fMaxEta);
0099 double phi = CLHEP::RandFlat::shoot(engine, fMinPhi, fMaxPhi);
0100 int ieta = static_cast<int>((eta - fEtaMin_) / fEtaBinWidth_);
0101 auto ipp = std::lower_bound(fP_.begin(), fP_.end(), mom);
0102 if (ipp == fP_.end())
0103 --ipp;
0104 int ip = static_cast<int>(ipp - fP_.begin());
0105 int indx = ip * kfactor + ieta;
0106 if (fVerbosity > 0)
0107 edm::LogVerbatim("ParticleGun") << "FileRandomMultiParticlePGunProducer: p " << mom << " Eta " << eta << " Phi "
0108 << phi << " Index " << indx;
0109
0110
0111
0112 int barcode(0), partID(fPartIDs[0]);
0113 double r1 = CLHEP::RandFlat::shoot(engine, 0., 1.);
0114 for (unsigned int ip = 0; ip < fPartIDs.size(); ip++) {
0115 if (r1 <= fProbParticle_[indx][ip])
0116 break;
0117 partID = fPartIDs[ip];
0118 }
0119 if (fVerbosity > 0)
0120 edm::LogVerbatim("ParticleGun") << "Random " << r1 << " PartID " << partID;
0121 const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(partID));
0122 double mass = PData->mass().value();
0123 double energy = sqrt(mom * mom + mass * mass);
0124 double theta = 2. * atan(exp(-eta));
0125 double px = mom * sin(theta) * cos(phi);
0126 double py = mom * sin(theta) * sin(phi);
0127 double pz = mom * cos(theta);
0128
0129 HepMC::FourVector p(px, py, pz, energy);
0130 HepMC::GenParticle* Part = new HepMC::GenParticle(p, partID, 1);
0131 barcode++;
0132 Part->suggest_barcode(barcode);
0133 Vtx->add_particle_out(Part);
0134
0135 fEvt->add_vertex(Vtx);
0136 fEvt->set_event_number(e.id().event());
0137 fEvt->set_signal_process_id(20);
0138
0139 if (fVerbosity > 1)
0140 fEvt->print();
0141
0142 std::unique_ptr<HepMCProduct> BProduct(new HepMCProduct());
0143 BProduct->addHepMCData(fEvt);
0144 e.put(std::move(BProduct), "unsmeared");
0145
0146 std::unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
0147 e.put(std::move(genEventInfo));
0148 if (fVerbosity > 0)
0149 edm::LogVerbatim("ParticleGun") << "FileRandomMultiParticlePGunProducer : Event Generation Done";
0150 }