ParticleGuns/src/CloseByParticleGunProducer.cc

0001 #include <ostream>
0002
0003 #include "IOMC/ParticleGuns/interface/CloseByParticleGunProducer.h"
0004
0005 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0006 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0007
0008 #include "DataFormats/Math/interface/Vector3D.h"
0009
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016 #include "CLHEP/Random/RandFlat.h"
0017 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0018 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0019 #include "CLHEP/Random/RandFlat.h"
0020
0021 using namespace edm;
0022 using namespace std;
0023
0024 CloseByParticleGunProducer::CloseByParticleGunProducer(const ParameterSet& pset) : BaseFlatGunProducer(pset) {
0025   ParameterSet pgun_params = pset.getParameter<ParameterSet>("PGunParameters");
0026   fControlledByEta = pgun_params.getParameter<bool>("ControlledByEta");
0027   fEnMax = pgun_params.getParameter<double>("EnMax");
0028   fEnMin = pgun_params.getParameter<double>("EnMin");
0029   fMaxEnSpread = pgun_params.getParameter<bool>("MaxEnSpread");
0030   if (fControlledByEta) {
0031     fEtaMax = pgun_params.getParameter<double>("MaxEta");
0032     fEtaMin = pgun_params.getParameter<double>("MinEta");
0033     if (fEtaMax <= fEtaMin)
0034       LogError("CloseByParticleGunProducer") << " Please fix MinEta and MaxEta values in the configuration";
0035   } else {
0036     fRMax = pgun_params.getParameter<double>("RMax");
0037     fRMin = pgun_params.getParameter<double>("RMin");
0038     if (fRMax <= fRMin)
0039       LogError("CloseByParticleGunProducer") << " Please fix RMin and RMax values in the configuration";
0040   }
0041   fZMax = pgun_params.getParameter<double>("ZMax");
0042   fZMin = pgun_params.getParameter<double>("ZMin");
0043   fDelta = pgun_params.getParameter<double>("Delta");
0044   fPhiMin = pgun_params.getParameter<double>("MinPhi");
0045   fPhiMax = pgun_params.getParameter<double>("MaxPhi");
0046   fPointing = pgun_params.getParameter<bool>("Pointing");
0047   fOverlapping = pgun_params.getParameter<bool>("Overlapping");
0048   fRandomShoot = pgun_params.getParameter<bool>("RandomShoot");
0049   fNParticles = pgun_params.getParameter<int>("NParticles");
0050   fPartIDs = pgun_params.getParameter<vector<int> >("PartID");
0051
0052   produces<HepMCProduct>("unsmeared");
0053   produces<GenEventInfoProduct>();
0054 }
0055
0056 CloseByParticleGunProducer::~CloseByParticleGunProducer() {
0057   // no need to cleanup GenEvent memory - done in HepMCProduct
0058 }
0059
0060 void CloseByParticleGunProducer::produce(Event& e, const EventSetup& es) {
0061   edm::Service<edm::RandomNumberGenerator> rng;
0062   CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
0063
0064   if (fVerbosity > 0) {
0065     LogDebug("CloseByParticleGunProducer") << " CloseByParticleGunProducer : Begin New Event Generation" << endl;
0066   }
0067   fEvt = new HepMC::GenEvent();
0068
0069   int barcode = 1;
0070   unsigned int numParticles = fRandomShoot ? CLHEP::RandFlat::shoot(engine, 1, fNParticles) : fNParticles;
0071
0072   double phi = CLHEP::RandFlat::shoot(engine, fPhiMin, fPhiMax);
0073   double fZ = CLHEP::RandFlat::shoot(engine, fZMin, fZMax);
0074   double fR;
0075
0076   if (!fControlledByEta) {
0077     fR = CLHEP::RandFlat::shoot(engine, fRMin, fRMax);
0078   } else {
0079     double fEta = CLHEP::RandFlat::shoot(engine, fEtaMin, fEtaMax);
0080     fR = (fZ / sinh(fEta));
0081   }
0082
0083   double tmpPhi = phi;
0084   double tmpR = fR;
0085
0086   // Loop over particles
0087   for (unsigned int ip = 0; ip < numParticles; ++ip) {
0088     if (fOverlapping) {
0089       fR = CLHEP::RandFlat::shoot(engine, tmpR - fDelta, tmpR + fDelta);
0090       phi = CLHEP::RandFlat::shoot(engine, tmpPhi - fDelta / fR, tmpPhi + fDelta / fR);
0091     } else
0092       phi += fDelta / fR;
0093
0094     double fEn;
0095     if (numParticles > 1 && fMaxEnSpread)
0096       fEn = fEnMin + ip * (fEnMax - fEnMin) / (numParticles - 1);
0097     else
0098       fEn = CLHEP::RandFlat::shoot(engine, fEnMin, fEnMax);
0099
0100     int partIdx = CLHEP::RandFlat::shoot(engine, 0, fPartIDs.size());
0101     int PartID = fPartIDs[partIdx];
0102     const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID)));
0103     double mass = PData->mass().value();
0104     double mom2 = fEn * fEn - mass * mass;
0105     double mom = 0.;
0106     if (mom2 > 0.) {
0107       mom = sqrt(mom2);
0108     }
0109     double px = 0.;
0110     double py = 0.;
0111     double pz = mom;
0112     double energy = fEn;
0113
0114     // Compute Vertex Position
0115     double x = fR * cos(phi);
0116     double y = fR * sin(phi);
0117     constexpr double c = 2.99792458e+1;  // cm/ns
0118     double timeOffset = sqrt(x * x + y * y + fZ * fZ) / c * ns * c_light;
0119     HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(x * cm, y * cm, fZ * cm, timeOffset));
0120
0121     HepMC::FourVector p(px, py, pz, energy);
0122     // If we are requested to be pointing to (0,0,0), correct the momentum direction
0123     if (fPointing) {
0124       math::XYZVector direction(x, y, fZ);
0125       math::XYZVector momentum = direction.unit() * mom;
0126       p.setX(momentum.x());
0127       p.setY(momentum.y());
0128       p.setZ(momentum.z());
0129     }
0130     HepMC::GenParticle* Part = new HepMC::GenParticle(p, PartID, 1);
0131     Part->suggest_barcode(barcode);
0132     barcode++;
0133
0134     Vtx->add_particle_out(Part);
0135
0136     if (fVerbosity > 0) {
0137       Vtx->print();
0138       Part->print();
0139     }
0140     fEvt->add_vertex(Vtx);
0141   }
0142
0143   fEvt->set_event_number(e.id().event());
0144   fEvt->set_signal_process_id(20);
0145
0146   if (fVerbosity > 0) {
0147     fEvt->print();
0148   }
0149
0150   unique_ptr<HepMCProduct> BProduct(new HepMCProduct());
0151   BProduct->addHepMCData(fEvt);
0152   e.put(std::move(BProduct), "unsmeared");
0153
0154   unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
0155   e.put(std::move(genEventInfo));
0156
0157   if (fVerbosity > 0) {
0158     LogDebug("CloseByParticleGunProducer") << " CloseByParticleGunProducer : Event Generation Done " << endl;
0159   }
0160 }