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