ParticleGuns/src/FileRandomKEThetaGunProducer.cc

0001 #include <iostream>
0002
0003 #include "IOMC/ParticleGuns/interface/FileRandomKEThetaGunProducer.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/Utilities/interface/Exception.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0014
0015 #include "CLHEP/Random/RandFlat.h"
0016 #include "CLHEP/Random/RandomEngine.h"
0017
0018 using namespace edm;
0019
0020 FileRandomKEThetaGunProducer::FileRandomKEThetaGunProducer(const edm::ParameterSet& pset)
0021     : FlatBaseThetaGunProducer(pset) {
0022   edm::ParameterSet defpset;
0023   edm::ParameterSet pgun_params = pset.getParameter<edm::ParameterSet>("PGunParameters");
0024
0025   edm::FileInPath fp = pgun_params.getParameter<edm::FileInPath>("File");
0026   std::string file = fp.fullPath();
0027   particleN = pgun_params.getParameter<int>("Particles");
0028   if (particleN <= 0)
0029     particleN = 1;
0030   edm::LogInfo("FlatThetaGun") << "Internal FileRandomKEThetaGun is initialzed"
0031                                << " with data read from " << file << " and " << particleN << " particles created/event";
0032   std::ifstream is(file.c_str(), std::ios::in);
0033   if (is) {
0034     double energy, elem, sum = 0;
0035     while (!is.eof()) {
0036       is >> energy >> elem;
0037       kineticE.push_back(0.001 * energy);
0038       fdistn.push_back(elem);
0039       sum += elem;
0040       if (fVerbosity > 0)
0041         LogDebug("FlatThetaGun") << "KE " << energy << " GeV Count rate " << elem;
0042     }
0043     is.close();
0044     double last = 0;
0045     for (unsigned int i = 0; i < fdistn.size(); i++) {
0046       fdistn[i] /= sum;
0047       fdistn[i] += last;
0048       last = fdistn[i];
0049       if (fVerbosity > 0)
0050         LogDebug("FlatThetaGun") << "Bin [" << i << "]: KE " << kineticE[i] << " Distn " << fdistn[i];
0051     }
0052   }
0053   if (kineticE.size() < 2)
0054     throw cms::Exception("FileNotFound", "Not enough data point found in file ") << file << "\n";
0055
0056   produces<HepMCProduct>("unsmeared");
0057   produces<GenEventInfoProduct>();
0058 }
0059
0060 FileRandomKEThetaGunProducer::~FileRandomKEThetaGunProducer() {}
0061
0062 void FileRandomKEThetaGunProducer::produce(edm::Event& e, const edm::EventSetup& es) {
0063   if (fVerbosity > 0)
0064     LogDebug("FlatThetaGun") << "FileRandomKEThetaGunProducer : Begin New Event Generation";
0065
0066   edm::Service<edm::RandomNumberGenerator> rng;
0067   CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
0068
0069   // event loop (well, another step in it...)
0070
0071   // no need to clean up GenEvent memory - done in HepMCProduct
0072
0073   // here re-create fEvt (memory)
0074   //
0075   fEvt = new HepMC::GenEvent();
0076
0077   // now actualy, cook up the event from PDGTable and gun parameters
0078   //
0079
0080   // 1st, primary vertex
0081   //
0082   HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(0., 0., 0.));
0083
0084   // loop over particles
0085   //
0086   int barcode = 1;
0087   for (int ip = 0; ip < particleN; ip++) {
0088     double keMin = kineticE[0], keMax = kineticE[1];
0089     double rMin = fdistn[0], rMax = fdistn[1];
0090     double r1 = engine->flat();
0091     for (unsigned int ii = kineticE.size() - 2; ii > 0; --ii) {
0092       if (r1 > fdistn[ii]) {
0093         keMin = kineticE[ii];
0094         keMax = kineticE[ii + 1];
0095         rMin = fdistn[ii];
0096         rMax = fdistn[ii + 1];
0097         break;
0098       }
0099     }
0100     double ke = (keMin * (rMax - r1) + keMax * (r1 - rMin)) / (rMax - rMin);
0101     if (fVerbosity > 1)
0102       LogDebug("FlatThetaGun") << "FileRandomKEThetaGunProducer: KE " << ke << " in range " << keMin << ":" << keMax
0103                                << " with " << r1 << " in " << rMin << ":" << rMax;
0104     double theta = CLHEP::RandFlat::shoot(engine, fMinTheta, fMaxTheta);
0105     double phi = CLHEP::RandFlat::shoot(engine, fMinPhi, fMaxPhi);
0106     int PartID = fPartIDs[0];
0107     const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID)));
0108     double mass = PData->mass().value();
0109     double energy = ke + mass;
0110     double mom2 = ke * ke + 2. * ke * mass;
0111     double mom = std::sqrt(mom2);
0112     double px = mom * sin(theta) * cos(phi);
0113     double py = mom * sin(theta) * sin(phi);
0114     double pz = mom * cos(theta);
0115
0116     HepMC::FourVector p(px, py, pz, energy);
0117     HepMC::GenParticle* Part = new HepMC::GenParticle(p, PartID, 1);
0118     Part->suggest_barcode(barcode);
0119     barcode++;
0120     Vtx->add_particle_out(Part);
0121   }
0122   fEvt->add_vertex(Vtx);
0123   fEvt->set_event_number(e.id().event());
0124   fEvt->set_signal_process_id(20);
0125
0126   if (fVerbosity > 0) {
0127     fEvt->print();
0128   }
0129
0130   std::unique_ptr<HepMCProduct> BProduct(new HepMCProduct());
0131   BProduct->addHepMCData(fEvt);
0132   e.put(std::move(BProduct), "unsmeared");
0133
0134   std::unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
0135   e.put(std::move(genEventInfo));
0136
0137   if (fVerbosity > 0)
0138     LogDebug("FlatThetaGun") << "FileRandomKEThetaGunProducer : Event Generation Done";
0139 }