ParticleGuns/src/FlatRandomMultiParticlePGunProducer.cc

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   // event loop (well, another step in it...)
0054   // no need to clean up GenEvent memory - done in HepMCProduct
0055   // here re-create fEvt (memory)
0056   //
0057   fEvt = new HepMC::GenEvent();
0058
0059   // now actualy, cook up the event from PDGTable and gun parameters
0060   //
0061
0062   // 1st, primary vertex
0063   //
0064   HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(0., 0., 0.));
0065
0066   // loop over particles
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 }