ParticleGuns/src/RandomMultiParticlePGunProducer.cc

0001 #include <ostream>
0002
0003 #include "IOMC/ParticleGuns/interface/RandomMultiParticlePGunProducer.h"
0004
0005 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0006 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0007
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/ServiceRegistry/interface/Service.h"
0012 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0013 #include "FWCore/Utilities/interface/EDMException.h"
0014 #include "CLHEP/Random/RandFlat.h"
0015
0016 //#define DebugLog
0017 using namespace edm;
0018
0019 RandomMultiParticlePGunProducer::RandomMultiParticlePGunProducer(const ParameterSet& pset) : BaseFlatGunProducer(pset) {
0020   ParameterSet pgunParams = pset.getParameter<ParameterSet>("PGunParameters");
0021   fProbParticle_ = pgunParams.getParameter<std::vector<double> >("ProbParts");
0022   fProbP_ = pgunParams.getParameter<std::vector<double> >("ProbP");
0023   fMinP_ = pgunParams.getParameter<double>("MinP");
0024   fMaxP_ = pgunParams.getParameter<double>("MaxP");
0025   if (fProbParticle_.size() != fPartIDs.size())
0026     throw cms::Exception("Configuration") << "Not all probabilities given for all particle types "
0027                                           << fProbParticle_.size() << ":" << fPartIDs.size() << " need them to match\n";
0028
0029   produces<HepMCProduct>("unsmeared");
0030   produces<GenEventInfoProduct>();
0031   fBinsP_ = (int)(fProbP_.size());
0032   fDeltaP_ = (fMaxP_ - fMinP_) / fBinsP_;
0033 #ifdef DebugLog
0034   edm::LogVerbatim("IOMC") << "Internal FlatRandomPGun is initialzed for " << fPartIDs.size() << " particles in "
0035                            << fBinsP_ << " bins within momentum range " << fMinP_ << ":" << fMaxP_;
0036   for (unsigned int k = 0; k < fPartIDs.size(); ++k)
0037     edm::LogVerbatim("IOMC") << " [" << k << "] " << fPartIDs[k] << ":" << fProbParticle_[k];
0038   edm::LogVerbatim("IOMC") << "Momentum distribution is given by";
0039   for (int k = 0; k < fBinsP_; ++k) {
0040     double p = fMinP_ + k * fDeltaP_;
0041     edm::LogVerbatim("IOMC") << " Bin[" << k << "] " << p << ":" << p + fDeltaP_ << " --> " << fProbP_[k];
0042   }
0043 #endif
0044   for (unsigned int k = 1; k < fProbParticle_.size(); ++k)
0045     fProbParticle_[k] += fProbParticle_[k - 1];
0046   for (unsigned int k = 0; k < fProbParticle_.size(); ++k)
0047     fProbParticle_[k] /= fProbParticle_[fProbParticle_.size() - 1];
0048 #ifdef DebugLog
0049   edm::LogVerbatim("IOMC") << "Corrected probabilities for particle type:";
0050   for (unsigned int k = 0; k < fProbParticle_.size(); ++k)
0051     edm::LogVerbatim("IOMC") << "  [" << k << "]: " << fProbParticle_[k];
0052 #endif
0053   for (int k = 1; k < fBinsP_; ++k)
0054     fProbP_[k] += fProbP_[k - 1];
0055   for (int k = 0; k < fBinsP_; ++k)
0056     fProbP_[k] /= fProbP_[fBinsP_ - 1];
0057 #ifdef DebugLog
0058   edm::LogVerbatim("IOMC") << "Corrected probabilities for momentum:";
0059   for (int k = 0; k < fBinsP_; ++k) {
0060     double p = fMinP_ + k * fDeltaP_;
0061     edm::LogVerbatim("IOMC") << " Bin[" << k << "] " << p << ":" << p + fDeltaP_ << " --> " << fProbP_[k];
0062   }
0063 #endif
0064 }
0065
0066 void RandomMultiParticlePGunProducer::produce(edm::Event& e, const edm::EventSetup& es) {
0067   edm::Service<edm::RandomNumberGenerator> rng;
0068   CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
0069
0070 #ifdef DebugLog
0071   if (fVerbosity > 0)
0072     edm::LogVerbatim("IOMC") << "RandomMultiParticlePGunProducer: "
0073                              << "Begin New Event Generation";
0074 #endif
0075
0076   // event loop (well, another step in it...)
0077   // no need to clean up GenEvent memory - done in HepMCProduct
0078   // here re-create fEvt (memory)
0079   //
0080   fEvt = new HepMC::GenEvent();
0081
0082   // now actualy, cook up the event from PDGTable and gun parameters
0083   //
0084
0085   // 1st, primary vertex
0086   //
0087   HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(0., 0., 0.));
0088
0089   // loop over particles
0090   //
0091   int barcode(0), PartID(fPartIDs[0]);
0092   double r1 = CLHEP::RandFlat::shoot(engine, 0., 1.);
0093   for (unsigned int ip = 0; ip < fPartIDs.size(); ip++) {
0094     if (r1 <= fProbParticle_[ip]) {
0095       PartID = fPartIDs[ip];
0096       break;
0097     }
0098   }
0099   double minP(fMinP_);
0100   double r2 = CLHEP::RandFlat::shoot(engine, 0., 1.);
0101   for (int ip = 0; ip < fBinsP_; ip++) {
0102     if (r2 <= fProbP_[ip]) {
0103       minP = fMinP_ + ip * fDeltaP_;
0104       break;
0105     }
0106   }
0107   double maxP = minP + fDeltaP_;
0108 #ifdef DebugLog
0109   if (fVerbosity > 0)
0110     edm::LogVerbatim("IOMC") << "Random " << r1 << " PartID " << PartID << " and for p " << r2 << " range " << minP
0111                              << ":" << maxP;
0112 #endif
0113   double mom = CLHEP::RandFlat::shoot(engine, minP, maxP);
0114   double eta = CLHEP::RandFlat::shoot(engine, fMinEta, fMaxEta);
0115   double phi = CLHEP::RandFlat::shoot(engine, fMinPhi, fMaxPhi);
0116   const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID)));
0117   double mass = PData->mass().value();
0118   double energy = sqrt(mom * mom + mass * mass);
0119   double theta = 2. * atan(exp(-eta));
0120   double px = mom * sin(theta) * cos(phi);
0121   double py = mom * sin(theta) * sin(phi);
0122   double pz = mom * cos(theta);
0123
0124   HepMC::FourVector p(px, py, pz, energy);
0125   HepMC::GenParticle* Part = new HepMC::GenParticle(p, PartID, 1);
0126   barcode++;
0127   Part->suggest_barcode(barcode);
0128   Vtx->add_particle_out(Part);
0129
0130   if (fAddAntiParticle) {
0131     HepMC::FourVector ap(-px, -py, -pz, energy);
0132     int APartID = (PartID == 22 || PartID == 23) ? PartID : -PartID;
0133     HepMC::GenParticle* APart = new HepMC::GenParticle(ap, APartID, 1);
0134     barcode++;
0135     APart->suggest_barcode(barcode);
0136     Vtx->add_particle_out(APart);
0137   }
0138
0139   fEvt->add_vertex(Vtx);
0140   fEvt->set_event_number(e.id().event());
0141   fEvt->set_signal_process_id(20);
0142
0143 #ifdef DebugLog
0144   if (fVerbosity > 0)
0145     fEvt->print();
0146 #endif
0147
0148   std::unique_ptr<HepMCProduct> BProduct(new HepMCProduct());
0149   BProduct->addHepMCData(fEvt);
0150   e.put(std::move(BProduct), "unsmeared");
0151
0152   std::unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
0153   e.put(std::move(genEventInfo));
0154 #ifdef DebugLog
0155   if (fVerbosity > 0)
0156     edm::LogVerbatim("IOMC") << " RandomMultiParticlePGunProducer : Event Generation Done ";
0157 #endif
0158 }