ParticleGuns/src/FlatRandomEGunProducer.cc

0001 /*
0002  *  \author Julia Yarba
0003  */
0004
0005 #include <ostream>
0006
0007 #include "IOMC/ParticleGuns/interface/FlatRandomEGunProducer.h"
0008
0009 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0010 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0011
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/ServiceRegistry/interface/Service.h"
0015 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0016
0017 #include "CLHEP/Random/RandFlat.h"
0018
0019 using namespace edm;
0020 using namespace std;
0021
0022 FlatRandomEGunProducer::FlatRandomEGunProducer(const ParameterSet& pset) : BaseFlatGunProducer(pset) {
0023   ParameterSet defpset;
0024   // ParameterSet pgun_params = pset.getParameter<ParameterSet>("PGunParameters") ;
0025   ParameterSet pgun_params = pset.getParameter<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   cout << "Internal FlatRandomEGun is initialzed" << endl;
0036   // cout << "It is going to generate " << remainingEvents() << "events" << endl ;
0037 }
0038
0039 FlatRandomEGunProducer::~FlatRandomEGunProducer() {
0040   // no need to cleanup fEvt since it's done in HepMCProduct
0041 }
0042
0043 void FlatRandomEGunProducer::produce(Event& e, const EventSetup& es) {
0044   edm::Service<edm::RandomNumberGenerator> rng;
0045   CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
0046
0047   if (fVerbosity > 0) {
0048     cout << " FlatRandomEGunProducer : Begin New Event Generation" << endl;
0049   }
0050
0051   // event loop (well, another step in it...)
0052
0053   // no need to clean up GenEvent memory - done in HepMCProduct
0054
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 = 1;
0069   for (unsigned int ip = 0; ip < fPartIDs.size(); ip++) {
0070     double energy = CLHEP::RandFlat::shoot(engine, fMinE, fMaxE);
0071     double eta = CLHEP::RandFlat::shoot(engine, fMinEta, fMaxEta);
0072     double phi = CLHEP::RandFlat::shoot(engine, fMinPhi, fMaxPhi);
0073     int PartID = fPartIDs[ip];
0074     const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID)));
0075     double mass = PData->mass().value();
0076     double mom2 = energy * energy - mass * mass;
0077     double mom = 0.;
0078     if (mom2 > 0.) {
0079       mom = sqrt(mom2);
0080     } else {
0081       mom = 0.;
0082     }
0083     double theta = 2. * atan(exp(-eta));
0084     double px = mom * sin(theta) * cos(phi);
0085     double py = mom * sin(theta) * sin(phi);
0086     double pz = mom * cos(theta);
0087
0088     HepMC::FourVector p(px, py, pz, energy);
0089     HepMC::GenParticle* Part = new HepMC::GenParticle(p, PartID, 1);
0090     Part->suggest_barcode(barcode);
0091     barcode++;
0092     Vtx->add_particle_out(Part);
0093
0094     if (fAddAntiParticle) {
0095       HepMC::FourVector ap(-px, -py, -pz, energy);
0096       int APartID = -PartID;
0097       if (PartID == 22 || PartID == 23) {
0098         APartID = PartID;
0099       }
0100       HepMC::GenParticle* APart = new HepMC::GenParticle(ap, APartID, 1);
0101       APart->suggest_barcode(barcode);
0102       barcode++;
0103       Vtx->add_particle_out(APart);
0104     }
0105   }
0106   fEvt->add_vertex(Vtx);
0107   fEvt->set_event_number(e.id().event());
0108   fEvt->set_signal_process_id(20);
0109
0110   if (fVerbosity > 0) {
0111     fEvt->print();
0112   }
0113
0114   unique_ptr<HepMCProduct> BProduct(new HepMCProduct());
0115   BProduct->addHepMCData(fEvt);
0116   e.put(std::move(BProduct), "unsmeared");
0117
0118   unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
0119   e.put(std::move(genEventInfo));
0120
0121   if (fVerbosity > 0) {
0122     // for testing purpose only
0123     //fEvt->print() ;  // for some strange reason, it prints NO event info
0124     // after it's been put into edm::Event...
0125     cout << " FlatRandomEGunProducer : Event Generation Done " << endl;
0126   }
0127 }