File indexing completed on 2022-01-29 01:30:47
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
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
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
0115 double x = fR * cos(phi);
0116 double y = fR * sin(phi);
0117 constexpr double c = 2.99792458e+1;
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
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 }