ParticleGuns/src/ExpoRandomPGunProducer.cc

0001 /*
0002  *  \author Jean-Roch Vlimant
0003  *  modified by S.Abdullin 04/02/2011
0004  */
0005
0006 #include <ostream>
0007
0008 #include "IOMC/ParticleGuns/interface/ExpoRandomPGunProducer.h"
0009
0010 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0011 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0012
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/ServiceRegistry/interface/Service.h"
0016 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0017
0018 #include "CLHEP/Random/RandFlat.h"
0019
0020 using namespace edm;
0021 using namespace std;
0022
0023 ExpoRandomPGunProducer::ExpoRandomPGunProducer(const ParameterSet& pset) : BaseFlatGunProducer(pset) {
0024   ParameterSet defpset;
0025   ParameterSet pgun_params = pset.getParameter<ParameterSet>("PGunParameters");
0026
0027   fMinP = pgun_params.getParameter<double>("MinP");
0028   fMaxP = pgun_params.getParameter<double>("MaxP");
0029
0030   produces<HepMCProduct>("unsmeared");
0031   produces<GenEventInfoProduct>();
0032 }
0033
0034 ExpoRandomPGunProducer::~ExpoRandomPGunProducer() {
0035   // no need to cleanup GenEvent memory - done in HepMCProduct
0036 }
0037
0038 void ExpoRandomPGunProducer::produce(Event& e, const EventSetup& es) {
0039   edm::Service<edm::RandomNumberGenerator> rng;
0040   CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
0041
0042   if (fVerbosity > 0) {
0043     std::cout << " ExpoRandomPGunProducer : Begin New Event Generation" << std::endl;
0044   }
0045   // event loop (well, another step in it...)
0046
0047   // no need to clean up GenEvent memory - done in HepMCProduct
0048   //
0049
0050   // here re-create fEvt (memory)
0051   //
0052
0053   fEvt = new HepMC::GenEvent();
0054
0055   // now actualy, cook up the event from PDGTable and gun parameters
0056   //
0057   // 1st, primary vertex
0058   //
0059   //HepMC::GenVertex* Vtx = new HepMC::GenVertex(CLHEP::HepLorentzVector(0.,0.,0.));
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 pmom = CLHEP::RandFlat::shoot(engine, fMinP, fMaxP);
0067     double y = (1. / fMinP) * CLHEP::RandFlat::shoot(engine, 0.0, 1.0);
0068     double f = 1. / pmom;
0069     bool accpt = (y < f);
0070     //shoot until in the designated range
0071     while ((pmom < fMinP || pmom > fMaxP) || !accpt) {
0072       pmom = CLHEP::RandFlat::shoot(engine, fMinP, fMaxP);
0073       y = (1. / fMinP) * CLHEP::RandFlat::shoot(engine, 0.0, 1.0);
0074       f = 1. / pmom;
0075       accpt = (y < f);
0076     }
0077
0078     double eta = CLHEP::RandFlat::shoot(engine, fMinEta, fMaxEta);
0079     double phi = CLHEP::RandFlat::shoot(engine, fMinPhi, fMaxPhi);
0080     int PartID = fPartIDs[ip];
0081     const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID)));
0082     double mass = PData->mass().value();
0083     double theta = 2. * atan(exp(-eta));
0084     double mom = pmom;
0085     double pt = mom * sin(theta);
0086     double px = pt * cos(phi);
0087     double py = pt * sin(phi);
0088     double pz = mom * cos(theta);
0089     double energy2 = mom * mom + mass * mass;
0090     double energy = sqrt(energy2);
0091     //CLHEP::Hep3Vector p(px,py,pz) ;
0092     //HepMC::GenParticle* Part =
0093     //    new HepMC::GenParticle(CLHEP::HepLorentzVector(p,energy),PartID,1);
0094     HepMC::FourVector p(px, py, pz, energy);
0095     HepMC::GenParticle* Part = new HepMC::GenParticle(p, PartID, 1);
0096     Part->suggest_barcode(barcode);
0097     barcode++;
0098     Vtx->add_particle_out(Part);
0099
0100     if (fAddAntiParticle) {
0101       //CLHEP::Hep3Vector ap(-px,-py,-pz) ;
0102       HepMC::FourVector ap(-px, -py, -pz, energy);
0103       int APartID = -PartID;
0104       if (PartID == 22 || PartID == 23) {
0105         APartID = PartID;
0106       }
0107       //HepMC::GenParticle* APart =
0108       //   new HepMC::GenParticle(CLHEP::HepLorentzVector(ap,energy),APartID,1);
0109       HepMC::GenParticle* APart = new HepMC::GenParticle(ap, APartID, 1);
0110       APart->suggest_barcode(barcode);
0111       barcode++;
0112       Vtx->add_particle_out(APart);
0113     }
0114   }
0115
0116   fEvt->add_vertex(Vtx);
0117   fEvt->set_event_number(e.id().event());
0118   fEvt->set_signal_process_id(20);
0119
0120   if (fVerbosity > 0) {
0121     fEvt->print();
0122   }
0123
0124   unique_ptr<HepMCProduct> BProduct(new HepMCProduct());
0125   BProduct->addHepMCData(fEvt);
0126   e.put(std::move(BProduct), "unsmeared");
0127
0128   unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
0129   e.put(std::move(genEventInfo));
0130
0131   if (fVerbosity > 0) {
0132     // for testing purpose only
0133     // fEvt->print() ; // prints empty info after it's made into edm::Event
0134     std::cout << " FlatRandomPGunProducer : Event Generation Done " << std::endl;
0135   }
0136 }