ParticleGuns/src/FileRandomMultiParticlePGunProducer.cc

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   // event loop (well, another step in it...)
0084   // no need to clean up GenEvent memory - done in HepMCProduct
0085   // here re-create fEvt (memory)
0086   //
0087   fEvt = new HepMC::GenEvent();
0088
0089   // now actualy, cook up the event from PDGTable and gun parameters
0090   //
0091
0092   // 1st, primary vertex
0093   //
0094   HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(0., 0., 0.));
0095
0096   // Now p, eta, phi
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   // Now particle id
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 }