ParticleGuns/src/FlatRandomOneOverPtGunProducer.cc

0001 #include <ostream>
0002
0003 #include "IOMC/ParticleGuns/interface/FlatRandomOneOverPtGunProducer.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 #include <CLHEP/Units/SystemOfUnits.h>
0016 using namespace edm;
0017
0018 FlatRandomOneOverPtGunProducer::FlatRandomOneOverPtGunProducer(const edm::ParameterSet& pset)
0019     : BaseFlatGunProducer(pset) {
0020   edm::ParameterSet defpset;
0021   edm::ParameterSet pgun_params = pset.getParameter<ParameterSet>("PGunParameters");
0022
0023   fMinOneOverPt = pgun_params.getParameter<double>("MinOneOverPt");
0024   fMaxOneOverPt = pgun_params.getParameter<double>("MaxOneOverPt");
0025
0026   produces<HepMCProduct>("unsmeared");
0027   produces<GenEventInfoProduct>();
0028
0029   edm::LogInfo("ParticleGun") << "FlatRandomOneOverPtGunProducer: initialized with minimum and maximum 1/pt "
0030                               << fMinOneOverPt << ":" << fMaxOneOverPt;
0031 }
0032
0033 FlatRandomOneOverPtGunProducer::~FlatRandomOneOverPtGunProducer() {
0034   // no need to cleanup GenEvent memory - done in HepMCProduct
0035 }
0036
0037 void FlatRandomOneOverPtGunProducer::produce(Event& e, const EventSetup& es) {
0038   edm::Service<edm::RandomNumberGenerator> rng;
0039   CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
0040
0041   LogDebug("ParticleGun") << " FlatRandomOneOverPtGunProducer : Begin New Event Generation";
0042
0043   // event loop (well, another step in it...)
0044
0045   // no need to clean up GenEvent memory - done in HepMCProduct
0046   //
0047
0048   // here re-create fEvt (memory)
0049   //
0050   fEvt = new HepMC::GenEvent();
0051
0052   // now actualy, cook up the event from PDGTable and gun parameters
0053   //
0054   // 1st, primary vertex
0055   //
0056   HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(0., 0., 0.));
0057
0058   // loop over particles
0059   //
0060   int barcode = 1;
0061   for (unsigned int ip = 0; ip < fPartIDs.size(); ++ip) {
0062     double xx = CLHEP::RandFlat::shoot(engine, 0.0, 1.0);
0063     double pt = std::exp((1. - xx) * std::log(fMinOneOverPt) + xx * std::log(fMaxOneOverPt));
0064     double eta = CLHEP::RandFlat::shoot(engine, fMinEta, fMaxEta);
0065     double phi = CLHEP::RandFlat::shoot(engine, fMinPhi, fMaxPhi);
0066     if (pt != 0)
0067       pt = 1. / pt;
0068     int PartID = fPartIDs[ip];
0069     const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID)));
0070     double mass = PData->mass().value();
0071     double theta = 2. * atan(exp(-eta));
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     LogDebug("ParticleGun") << "FlatRandomOneOverPtGunProducer: Event generated with pt:eta:phi " << pt << " " << eta
0084                             << " " << phi << " (" << theta / CLHEP::deg << ":" << phi / CLHEP::deg << ")";
0085
0086     if (fAddAntiParticle) {
0087       HepMC::FourVector ap(-px, -py, -pz, energy);
0088       int APartID = -PartID;
0089       if (PartID == 22 || PartID == 23) {
0090         APartID = PartID;
0091       }
0092       HepMC::GenParticle* APart = new HepMC::GenParticle(ap, APartID, 1);
0093       APart->suggest_barcode(barcode);
0094       barcode++;
0095       Vtx->add_particle_out(APart);
0096     }
0097   }
0098
0099   fEvt->add_vertex(Vtx);
0100   fEvt->set_event_number(e.id().event());
0101   fEvt->set_signal_process_id(20);
0102
0103   if (fVerbosity > 0) {
0104     fEvt->print();
0105   }
0106
0107   std::unique_ptr<HepMCProduct> BProduct(new HepMCProduct());
0108   BProduct->addHepMCData(fEvt);
0109   e.put(std::move(BProduct), "unsmeared");
0110
0111   std::unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
0112   e.put(std::move(genEventInfo));
0113
0114   LogDebug("ParticleGun") << " FlatRandomOneOverPtGunProducer : Event Generation Done ";
0115 }