ParticleGuns/src/ExpoRandomPtGunProducer.cc

0001 /*
0002  *  \author Jean-Roch Vlimant
0003  */
0004
0005 #include <ostream>
0006
0007 #include "IOMC/ParticleGuns/interface/ExpoRandomPtGunProducer.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/RandExponential.h"
0018 #include "CLHEP/Random/RandFlat.h"
0019
0020 using namespace edm;
0021 using namespace std;
0022
0023 ExpoRandomPtGunProducer::ExpoRandomPtGunProducer(const ParameterSet& pset) : BaseFlatGunProducer(pset) {
0024   ParameterSet defpset;
0025   ParameterSet pgun_params = pset.getParameter<ParameterSet>("PGunParameters");
0026
0027   fMinPt = pgun_params.getParameter<double>("MinPt");
0028   fMaxPt = pgun_params.getParameter<double>("MaxPt");
0029
0030   fMeanPt = pgun_params.getParameter<double>("MeanPt");
0031
0032   produces<HepMCProduct>("unsmeared");
0033   produces<GenEventInfoProduct>();
0034 }
0035
0036 ExpoRandomPtGunProducer::~ExpoRandomPtGunProducer() {
0037   // no need to cleanup GenEvent memory - done in HepMCProduct
0038 }
0039
0040 void ExpoRandomPtGunProducer::produce(Event& e, const EventSetup& es) {
0041   edm::Service<edm::RandomNumberGenerator> rng;
0042   CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
0043
0044   if (fVerbosity > 0) {
0045     cout << " ExpoRandomPtGunProducer : Begin New Event Generation" << endl;
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(CLHEP::HepLorentzVector(0.,0.,0.));
0061   HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(0., 0., 0.));
0062
0063   // loop over particles
0064   //
0065   int barcode = 1;
0066   for (unsigned int ip = 0; ip < fPartIDs.size(); ++ip) {
0067     //the max is to ensure you don't generate at 0
0068     //the 90% is to get rid of edge effect
0069
0070     double pt = std::max(0.00001, 0.90 * fMinPt) + CLHEP::RandExponential::shoot(engine, fMeanPt);
0071     //shoot until in the designated range
0072     while (pt < fMinPt || pt > fMaxPt) {
0073       pt = std::max(0.00001, 0.90 * fMinPt) + CLHEP::RandExponential::shoot(engine, fMeanPt);
0074     }
0075
0076     double eta = CLHEP::RandFlat::shoot(engine, fMinEta, fMaxEta);
0077     double phi = CLHEP::RandFlat::shoot(engine, fMinPhi, fMaxPhi);
0078     int PartID = fPartIDs[ip];
0079     const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID)));
0080     double mass = PData->mass().value();
0081     double theta = 2. * atan(exp(-eta));
0082     double mom = pt / sin(theta);
0083     double px = pt * cos(phi);
0084     double py = pt * sin(phi);
0085     double pz = mom * cos(theta);
0086     double energy2 = mom * mom + mass * mass;
0087     double energy = sqrt(energy2);
0088     //CLHEP::Hep3Vector p(px,py,pz) ;
0089     //HepMC::GenParticle* Part =
0090     //    new HepMC::GenParticle(CLHEP::HepLorentzVector(p,energy),PartID,1);
0091     HepMC::FourVector p(px, py, pz, energy);
0092     HepMC::GenParticle* Part = new HepMC::GenParticle(p, PartID, 1);
0093     Part->suggest_barcode(barcode);
0094     barcode++;
0095     Vtx->add_particle_out(Part);
0096
0097     if (fAddAntiParticle) {
0098       //CLHEP::Hep3Vector ap(-px,-py,-pz) ;
0099       HepMC::FourVector ap(-px, -py, -pz, energy);
0100       int APartID = -PartID;
0101       if (PartID == 22 || PartID == 23) {
0102         APartID = PartID;
0103       }
0104       //HepMC::GenParticle* APart =
0105       //   new HepMC::GenParticle(CLHEP::HepLorentzVector(ap,energy),APartID,1);
0106       HepMC::GenParticle* APart = new HepMC::GenParticle(ap, APartID, 1);
0107       APart->suggest_barcode(barcode);
0108       barcode++;
0109       Vtx->add_particle_out(APart);
0110     }
0111   }
0112
0113   fEvt->add_vertex(Vtx);
0114   fEvt->set_event_number(e.id().event());
0115   fEvt->set_signal_process_id(20);
0116
0117   if (fVerbosity > 0) {
0118     fEvt->print();
0119   }
0120
0121   unique_ptr<HepMCProduct> BProduct(new HepMCProduct());
0122   BProduct->addHepMCData(fEvt);
0123   e.put(std::move(BProduct), "unsmeared");
0124
0125   unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
0126   e.put(std::move(genEventInfo));
0127
0128   if (fVerbosity > 0) {
0129     // for testing purpose only
0130     // fEvt->print() ; // prints empty info after it's made into edm::Event
0131     cout << " FlatRandomPtGunProducer : Event Generation Done " << endl;
0132   }
0133 }