ParticleGuns/src/FlatRandomPtThetaGunProducer.cc

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