ParticleGuns/src/FlatRandomEThetaGunProducer.cc

0001 #include <ostream>
0002
0003 #include "IOMC/ParticleGuns/interface/FlatRandomEThetaGunProducer.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/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013
0014 #include "CLHEP/Random/RandFlat.h"
0015
0016 namespace CLHEP {
0017   class HepRandomEngine;
0018 }
0019
0020 using namespace edm;
0021
0022 FlatRandomEThetaGunProducer::FlatRandomEThetaGunProducer(const edm::ParameterSet& pset)
0023     : FlatBaseThetaGunProducer(pset) {
0024   edm::ParameterSet defpset;
0025   edm::ParameterSet pgun_params = pset.getParameter<edm::ParameterSet>("PGunParameters");
0026
0027   // doesn't seem necessary to check if pset is empty - if this
0028   // is the case, default values will be taken for params
0029   fMinE = pgun_params.getParameter<double>("MinE");
0030   fMaxE = pgun_params.getParameter<double>("MaxE");
0031
0032   produces<HepMCProduct>("unsmeared");
0033   produces<GenEventInfoProduct>();
0034 }
0035
0036 FlatRandomEThetaGunProducer::~FlatRandomEThetaGunProducer() {}
0037
0038 void FlatRandomEThetaGunProducer::produce(edm::Event& e, const edm::EventSetup& es) {
0039   if (fVerbosity > 0) {
0040     edm::LogVerbatim("FlatThetaGun") << "FlatRandomEThetaGunProducer : Begin New Event Generation";
0041   }
0042
0043   edm::Service<edm::RandomNumberGenerator> rng;
0044   CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
0045
0046   // event loop (well, another step in it...)
0047
0048   // no need to clean up GenEvent memory - done in HepMCProduct
0049
0050   // here re-create fEvt (memory)
0051   //
0052   fEvt = new HepMC::GenEvent();
0053
0054   // now actualy, cook up the event from PDGTable and gun parameters
0055   //
0056
0057   // 1st, primary vertex
0058   //
0059   HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(0., 0., 0.));
0060
0061   // loop over particles
0062   //
0063   int barcode = 1;
0064   for (unsigned int ip = 0; ip < fPartIDs.size(); ip++) {
0065     double energy = CLHEP::RandFlat::shoot(engine, fMinE, fMaxE);
0066     double theta = CLHEP::RandFlat::shoot(engine, fMinTheta, fMaxTheta);
0067     double phi = CLHEP::RandFlat::shoot(engine, fMinPhi, fMaxPhi);
0068     int PartID = fPartIDs[ip];
0069     const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID)));
0070     double mass = PData->mass().value();
0071     double mom2 = energy * energy - mass * mass;
0072     double mom = (mom2 > 0. ? std::sqrt(mom2) : 0.);
0073     double px = mom * sin(theta) * cos(phi);
0074     double py = mom * sin(theta) * sin(phi);
0075     double pz = mom * cos(theta);
0076
0077     HepMC::FourVector p(px, py, pz, energy);
0078     HepMC::GenParticle* Part = new HepMC::GenParticle(p, PartID, 1);
0079     Part->suggest_barcode(barcode);
0080     barcode++;
0081     Vtx->add_particle_out(Part);
0082
0083     if (fAddAntiParticle) {
0084       HepMC::FourVector ap(-px, -py, -pz, energy);
0085       int APartID = -PartID;
0086       if (PartID == 22 || PartID == 23) {
0087         APartID = PartID;
0088       }
0089       HepMC::GenParticle* APart = new HepMC::GenParticle(ap, APartID, 1);
0090       APart->suggest_barcode(barcode);
0091       barcode++;
0092       Vtx->add_particle_out(APart);
0093     }
0094   }
0095   fEvt->add_vertex(Vtx);
0096   fEvt->set_event_number(e.id().event());
0097   fEvt->set_signal_process_id(20);
0098
0099   if (fVerbosity > 0) {
0100     fEvt->print();
0101   }
0102
0103   std::unique_ptr<HepMCProduct> BProduct(new HepMCProduct());
0104   BProduct->addHepMCData(fEvt);
0105   e.put(std::move(BProduct), "unsmeared");
0106
0107   std::unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
0108   e.put(std::move(genEventInfo));
0109
0110   if (fVerbosity > 0) {
0111     edm::LogVerbatim("FlatThetaGun") << "FlatRandomEThetaGunProducer : Event Generation Done";
0112   }
0113 }