ParticleGuns/src/CloseByParticleGunProducer.cc

0001 #include <ostream>
0002 #include <cmath>
0003
0004 #include "IOMC/ParticleGuns/interface/CloseByParticleGunProducer.h"
0005
0006 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0007 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0008
0009 #include "DataFormats/Math/interface/Vector3D.h"
0010
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "FWCore/ServiceRegistry/interface/Service.h"
0014 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016
0017 #include <CLHEP/Random/RandFlat.h>
0018 #include <CLHEP/Units/SystemOfUnits.h>
0019 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0020 #include <CLHEP/Random/RandFlat.h>
0021
0022 using namespace edm;
0023 using namespace std;
0024
0025 CloseByParticleGunProducer::CloseByParticleGunProducer(const ParameterSet& pset)
0026     : BaseFlatGunProducer(pset), m_fieldToken(esConsumes()) {
0027   ParameterSet pgun_params = pset.getParameter<ParameterSet>("PGunParameters");
0028   fControlledByEta = pgun_params.getParameter<bool>("ControlledByEta");
0029   fVarMax = pgun_params.getParameter<double>("VarMax");
0030   fVarMin = pgun_params.getParameter<double>("VarMin");
0031   fMaxVarSpread = pgun_params.getParameter<bool>("MaxVarSpread");
0032   fFlatPtGeneration = pgun_params.getParameter<bool>("FlatPtGeneration");
0033   if (fVarMin < 1 && !fFlatPtGeneration)
0034     LogError("CloseByParticleGunProducer") << " Please choose a minimum energy greater than 1 GeV, otherwise time "
0035                                               "information may be invalid or not reliable";
0036   if (fControlledByEta) {
0037     fEtaMax = pgun_params.getParameter<double>("MaxEta");
0038     fEtaMin = pgun_params.getParameter<double>("MinEta");
0039     if (fEtaMax <= fEtaMin)
0040       LogError("CloseByParticleGunProducer") << " Please fix MinEta and MaxEta values in the configuration";
0041   } else {
0042     fRMax = pgun_params.getParameter<double>("RMax");
0043     fRMin = pgun_params.getParameter<double>("RMin");
0044     if (fRMax <= fRMin)
0045       LogError("CloseByParticleGunProducer") << " Please fix RMin and RMax values in the configuration";
0046   }
0047   fZMax = pgun_params.getParameter<double>("ZMax");
0048   fZMin = pgun_params.getParameter<double>("ZMin");
0049   fDelta = pgun_params.getParameter<double>("Delta");
0050   fPhiMin = pgun_params.getParameter<double>("MinPhi");
0051   fPhiMax = pgun_params.getParameter<double>("MaxPhi");
0052   fPointing = pgun_params.getParameter<bool>("Pointing");
0053   fOverlapping = pgun_params.getParameter<bool>("Overlapping");
0054   if (fFlatPtGeneration && !fPointing)
0055     LogError("CloseByParticleGunProducer")
0056         << " Can't generate non pointing FlatPt samples; please disable FlatPt generation or generate pointing sample";
0057   fRandomShoot = pgun_params.getParameter<bool>("RandomShoot");
0058   fNParticles = pgun_params.getParameter<int>("NParticles");
0059   fPartIDs = pgun_params.getParameter<vector<int>>("PartID");
0060
0061   // set dt between particles
0062   fUseDeltaT = pgun_params.getParameter<bool>("UseDeltaT");
0063   fTMax = pgun_params.getParameter<double>("TMax");
0064   fTMin = pgun_params.getParameter<double>("TMin");
0065   if (fTMax <= fTMin)
0066     LogError("CloseByParticleGunProducer") << " Please fix TMin and TMax values in the configuration";
0067   // set a fixed time offset for the particles
0068   fOffsetFirst = pgun_params.getParameter<double>("OffsetFirst");
0069
0070   produces<HepMCProduct>("unsmeared");
0071   produces<GenEventInfoProduct>();
0072 }
0073
0074 CloseByParticleGunProducer::~CloseByParticleGunProducer() {
0075   // no need to cleanup GenEvent memory - done in HepMCProduct
0076 }
0077
0078 void CloseByParticleGunProducer::fillDescriptions(ConfigurationDescriptions& descriptions) {
0079   edm::ParameterSetDescription desc;
0080   desc.add<bool>("AddAntiParticle", false);
0081   {
0082     edm::ParameterSetDescription psd0;
0083     psd0.add<bool>("ControlledByEta", false);
0084     psd0.add<double>("Delta", 10);
0085     psd0.add<double>("VarMax", 200.0);
0086     psd0.add<double>("VarMin", 25.0);
0087     psd0.add<bool>("MaxVarSpread", false);
0088     psd0.add<bool>("FlatPtGeneration", false);
0089     psd0.add<double>("MaxEta", 2.7);
0090     psd0.add<double>("MaxPhi", 3.14159265359);
0091     psd0.add<double>("MinEta", 1.7);
0092     psd0.add<double>("MinPhi", -3.14159265359);
0093     psd0.add<int>("NParticles", 2);
0094     psd0.add<bool>("Overlapping", false);
0095     psd0.add<std::vector<int>>("PartID",
0096                                {
0097                                    22,
0098                                });
0099     psd0.add<bool>("Pointing", true);
0100     psd0.add<double>("RMax", 120);
0101     psd0.add<double>("RMin", 60);
0102     psd0.add<bool>("RandomShoot", false);
0103     psd0.add<double>("ZMax", 321);
0104     psd0.add<double>("ZMin", 320);
0105     psd0.add<bool>("UseDeltaT", false);
0106     psd0.add<double>("TMin", 0.);
0107     psd0.add<double>("TMax", 0.05);
0108     psd0.add<double>("OffsetFirst", 0.);
0109     desc.add<edm::ParameterSetDescription>("PGunParameters", psd0);
0110   }
0111   desc.addUntracked<int>("Verbosity", 0);
0112   desc.addUntracked<unsigned int>("firstRun", 1);
0113   desc.add<std::string>("psethack", "random particles in phi and r windows");
0114   descriptions.add("CloseByParticleGunProducer", desc);
0115 }
0116
0117 void CloseByParticleGunProducer::produce(Event& e, const EventSetup& es) {
0118   edm::Service<edm::RandomNumberGenerator> rng;
0119   CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
0120
0121   if (fVerbosity > 0) {
0122     LogDebug("CloseByParticleGunProducer") << " CloseByParticleGunProducer : Begin New Event Generation" << endl;
0123   }
0124   fEvt = new HepMC::GenEvent();
0125
0126   auto const& field = es.getData(m_fieldToken);
0127
0128   int barcode = 1;
0129   unsigned int numParticles = fRandomShoot ? CLHEP::RandFlat::shoot(engine, 1, fNParticles) : fNParticles;
0130
0131   double phi = CLHEP::RandFlat::shoot(engine, fPhiMin, fPhiMax);
0132   double fZ = CLHEP::RandFlat::shoot(engine, fZMin, fZMax);
0133   double fR, fEta;
0134   double fT;
0135
0136   if (!fControlledByEta) {
0137     fR = CLHEP::RandFlat::shoot(engine, fRMin, fRMax);
0138     fEta = asinh(fZ / fR);
0139   } else {
0140     fEta = CLHEP::RandFlat::shoot(engine, fEtaMin, fEtaMax);
0141     fR = (fZ / sinh(fEta));
0142   }
0143
0144   if (fUseDeltaT) {
0145     fT = CLHEP::RandFlat::shoot(engine, fTMin, fTMax);
0146   } else {
0147     fT = 0.;
0148   }
0149
0150   double tmpPhi = phi;
0151   double tmpR = fR;
0152
0153   // Loop over particles
0154   for (unsigned int ip = 0; ip < numParticles; ++ip) {
0155     if (fOverlapping) {
0156       fR = CLHEP::RandFlat::shoot(engine, tmpR - fDelta, tmpR + fDelta);
0157       phi = CLHEP::RandFlat::shoot(engine, tmpPhi - fDelta / fR, tmpPhi + fDelta / fR);
0158     } else
0159       phi += fDelta / fR;
0160
0161     double fVar;
0162     if (numParticles > 1 && fMaxVarSpread)
0163       fVar = fVarMin + ip * (fVarMax - fVarMin) / (numParticles - 1);
0164     else
0165       fVar = CLHEP::RandFlat::shoot(engine, fVarMin, fVarMax);
0166
0167     int partIdx = CLHEP::RandFlat::shoot(engine, 0, fPartIDs.size());
0168     int PartID = fPartIDs[partIdx];
0169     const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID)));
0170     double mass = PData->mass().value();
0171
0172     double mom, px, py, pz;
0173     double energy;
0174
0175     if (!fFlatPtGeneration) {
0176       double mom2 = fVar * fVar - mass * mass;
0177       mom = 0.;
0178       if (mom2 > 0.) {
0179         mom = sqrt(mom2);
0180       }
0181       px = 0.;
0182       py = 0.;
0183       pz = mom;
0184       energy = fVar;
0185     } else {
0186       double theta = 2. * atan(exp(-fEta));
0187       mom = fVar / sin(theta);
0188       px = fVar * cos(phi);
0189       py = fVar * sin(phi);
0190       pz = mom * cos(theta);
0191       double energy2 = mom * mom + mass * mass;
0192       energy = sqrt(energy2);
0193     }
0194     // Compute Vertex Position
0195     double x = fR * cos(phi);
0196     double y = fR * sin(phi);
0197
0198     HepMC::FourVector p(px, py, pz, energy);
0199     // If we are requested to be pointing to (0,0,0), correct the momentum direction
0200     if (fPointing) {
0201       math::XYZVector direction(x, y, fZ);
0202       math::XYZVector momentum = direction.unit() * mom;
0203       p.setX(momentum.x());
0204       p.setY(momentum.y());
0205       p.setZ(momentum.z());
0206     }
0207
0208     // compute correct path assuming uniform magnetic field in CMS
0209     double pathLength = 0.;
0210     const double speed = p.pz() / p.e() * c_light / CLHEP::cm;
0211     if (PData->charge()) {
0212       // Radius [cm] = P[GeV/c] * 10^9 / (c[mm/ns] * 10^6 * q[C] * B[T]) * 100[cm/m]
0213       const double radius = std::sqrt(p.px() * p.px() + p.py() * p.py()) * std::pow(10, 5) /
0214                             (c_light * field.inTesla({0.f, 0.f, 0.f}).z());  // cm
0215       const double arc = 2 * asinf(std::sqrt(x * x + y * y) / (2 * radius)) * radius;
0216       pathLength = std::sqrt(arc * arc + fZ * fZ);
0217     } else {
0218       pathLength = std::sqrt(x * x + y * y + fZ * fZ);
0219     }
0220
0221     // if not pointing time doesn't mean a lot, keep the old way
0222     const double pathTime = fPointing ? (pathLength / speed) : (std::sqrt(x * x + y * y + fZ * fZ) / speed);
0223     double timeOffset = fOffsetFirst + (pathTime + ip * fT) * CLHEP::ns * c_light;
0224
0225     HepMC::GenVertex* Vtx =
0226         new HepMC::GenVertex(HepMC::FourVector(x * CLHEP::cm, y * CLHEP::cm, fZ * CLHEP::cm, timeOffset));
0227
0228     HepMC::GenParticle* Part = new HepMC::GenParticle(p, PartID, 1);
0229     Part->suggest_barcode(barcode);
0230     barcode++;
0231
0232     Vtx->add_particle_out(Part);
0233
0234     if (fVerbosity > 0) {
0235       Vtx->print();
0236       Part->print();
0237     }
0238     fEvt->add_vertex(Vtx);
0239   }
0240
0241   fEvt->set_event_number(e.id().event());
0242   fEvt->set_signal_process_id(20);
0243
0244   if (fVerbosity > 0) {
0245     fEvt->print();
0246   }
0247
0248   unique_ptr<HepMCProduct> BProduct(new HepMCProduct());
0249   BProduct->addHepMCData(fEvt);
0250   e.put(std::move(BProduct), "unsmeared");
0251
0252   unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
0253   e.put(std::move(genEventInfo));
0254
0255   if (fVerbosity > 0) {
0256     LogDebug("CloseByParticleGunProducer") << " CloseByParticleGunProducer : Event Generation Done " << endl;
0257   }
0258 }