ParticleGuns/src/FlatRandomPtGunProducer.cc

0001 /*
0002  *  \author Julia Yarba
0003  */
0004
0005 #include <ostream>
0006
0007 #include "IOMC/ParticleGuns/interface/FlatRandomPtGunProducer.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 FlatRandomPtGunProducer::FlatRandomPtGunProducer(const ParameterSet& pset) : BaseFlatGunProducer(pset) {
0023   ParameterSet defpset;
0024   ParameterSet pgun_params = pset.getParameter<ParameterSet>("PGunParameters");
0025
0026   fMinPt = pgun_params.getParameter<double>("MinPt");
0027   fMaxPt = pgun_params.getParameter<double>("MaxPt");
0028
0029   produces<HepMCProduct>("unsmeared");
0030   produces<GenEventInfoProduct>();
0031 }
0032
0033 FlatRandomPtGunProducer::~FlatRandomPtGunProducer() {
0034   // no need to cleanup GenEvent memory - done in HepMCProduct
0035 }
0036
0037 void FlatRandomPtGunProducer::produce(Event& e, const EventSetup& es) {
0038   edm::Service<edm::RandomNumberGenerator> rng;
0039   CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
0040
0041   if (fVerbosity > 0) {
0042     cout << " FlatRandomPtGunProducer : Begin New Event Generation" << endl;
0043   }
0044   // event loop (well, another step in it...)
0045
0046   // no need to clean up GenEvent memory - done in HepMCProduct
0047   //
0048
0049   // here re-create fEvt (memory)
0050   //
0051   fEvt = new HepMC::GenEvent();
0052
0053   // now actualy, cook up the event from PDGTable and gun parameters
0054   //
0055   // 1st, primary vertex
0056   //
0057   HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(0., 0., 0.));
0058
0059   // loop over particles
0060   //
0061   int barcode = 1;
0062   for (unsigned int ip = 0; ip < fPartIDs.size(); ++ip) {
0063     double pt = CLHEP::RandFlat::shoot(engine, fMinPt, fMaxPt);
0064     double eta = CLHEP::RandFlat::shoot(engine, fMinEta, fMaxEta);
0065     double phi = CLHEP::RandFlat::shoot(engine, fMinPhi, fMaxPhi);
0066     int PartID = fPartIDs[ip];
0067     const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID)));
0068     double mass = PData->mass().value();
0069     double theta = 2. * atan(exp(-eta));
0070     double mom = pt / sin(theta);
0071     double px = pt * cos(phi);
0072     double py = pt * sin(phi);
0073     double pz = mom * cos(theta);
0074     double energy2 = mom * mom + mass * mass;
0075     double energy = sqrt(energy2);
0076     HepMC::FourVector p(px, py, pz, energy);
0077     HepMC::GenParticle* Part = new HepMC::GenParticle(p, PartID, 1);
0078     Part->suggest_barcode(barcode);
0079     barcode++;
0080     Vtx->add_particle_out(Part);
0081
0082     if (fAddAntiParticle) {
0083       HepMC::FourVector ap(-px, -py, -pz, energy);
0084       int APartID = -PartID;
0085       if (PartID == 22 || PartID == 23) {
0086         APartID = PartID;
0087       }
0088       HepMC::GenParticle* APart = new HepMC::GenParticle(ap, APartID, 1);
0089       APart->suggest_barcode(barcode);
0090       barcode++;
0091       Vtx->add_particle_out(APart);
0092     }
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   unique_ptr<HepMCProduct> BProduct(new HepMCProduct());
0104   BProduct->addHepMCData(fEvt);
0105   e.put(std::move(BProduct), "unsmeared");
0106
0107   unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
0108   e.put(std::move(genEventInfo));
0109
0110   if (fVerbosity > 0) {
0111     // for testing purpose only
0112     // fEvt->print() ; // prints empty info after it's made into edm::Event
0113     cout << " FlatRandomPtGunProducer : Event Generation Done " << endl;
0114   }
0115 }
0116 //#include "FWCore/Framework/interface/MakerMacros.h"
0117 //DEFINE_FWK_MODULE(FlatRandomPtGunProducer);