diff --git a/Configuration/Generator/python/DisplacedParticleGun_cfi.py b/Configuration/Generator/python/DisplacedParticleGun_cfi.py new file mode 100644 index 0000000000000..8ada6dd95972e --- /dev/null +++ b/Configuration/Generator/python/DisplacedParticleGun_cfi.py @@ -0,0 +1,41 @@ +import FWCore.ParameterSet.Config as cms + +generator = cms.EDProducer( + "DisplacedParticleGunProducer", + PGunParameters = cms.PSet( + # particle direction + MinPt = cms.double(5.0), + MaxPt = cms.double(100.0), + MinPhi = cms.double(-3.141592653589793), + MaxPhi = cms.double(+3.141592653589793), + + # displaced vertex in transverse plane (cm) + RMin = cms.double(0.0), + RMax = cms.double(10.0), + MinVtxPhi = cms.double(0.0), + MaxVtxPhi = cms.double(2.0 * 3.141592653589793), + ZVtx = cms.double(0.0), + + NParticles = cms.int32(1), + PartID = cms.vint32(22), # photon + + # how to sample the vertex radius + UniformDensityInR = cms.bool(False), + + # if True: derive theta to hit a region of the HGCAL front surface + # if False: use MinTheta/MaxTheta + PointingToHGCAL = cms.bool(True), + + # only used if PointingToHGCAL == True (cm), + # corresponds to the central third of HGCAL's front surface + RminFrontSurfaceHGCAL = cms.double(58.79), + RmaxFrontSurfaceHGCAL = cms.double(91.58), + + # only used if PointingToHGCAL == False (radians) + MinTheta = cms.double(0.), + MaxTheta = cms.double(3.141592653589793), + ), + Verbosity = cms.untracked.int32(0), + firstRun = cms.untracked.uint32(1), + psethack = cms.string("displaced gun with theta, optionally pointing to hard-coded HGCAL front surface"), +) diff --git a/Configuration/PyReleaseValidation/python/relval_steps.py b/Configuration/PyReleaseValidation/python/relval_steps.py index 6b8d833fafc84..77aca43759288 100644 --- a/Configuration/PyReleaseValidation/python/relval_steps.py +++ b/Configuration/PyReleaseValidation/python/relval_steps.py @@ -4950,6 +4950,9 @@ def gen2024HiMix(fragment,howMuch): upgradeStepDict['GenSimCloseBy'][k] = deepcopy(upgradeStepDict['GenSim'][k]) upgradeStepDict['GenSimCloseBy'][k]['--beamspot'] = 'CloseBy' + + upgradeStepDict['GenSimDisplaced'][k] = deepcopy(upgradeStepDict['GenSim'][k]) + upgradeStepDict['GenSimDisplaced'][k]['--beamspot'] = 'CloseBy' upgradeStepDict['GenSimHLBeamSpot'][k] = {'-s' : 'GEN,SIM', '-n' : 10, @@ -4963,7 +4966,7 @@ def gen2024HiMix(fragment,howMuch): upgradeStepDict['GenSimHLBeamSpot14'][k] = deepcopy(upgradeStepDict['GenSimHLBeamSpot'][k]) upgradeStepDict['GenSimHLBeamSpot14'][k]['--conditions'] = gt - upgradeStepDict['GenSimHLBeamSpotCloseBy'][k] = upgradeStepDict['GenSimCloseBy'][k] + upgradeStepDict['GenSimHLBeamSpotCloseBy'][k] = deepcopy(upgradeStepDict['GenSimCloseBy'][k]) upgradeStepDict['Sim'][k] = {'-s' : 'SIM', '-n' : 10, @@ -5201,7 +5204,7 @@ def gen2024HiMix(fragment,howMuch): # in case special WF has PU-specific changes: apply *after* basic PU step is created specialWF.setupPU(upgradeStepDict, k, upgradeProperties[year][k]) - + for step in upgradeStepDict.keys(): # we need to do this for each fragment if ('Sim' in step and ('Fast' not in step and step != 'Sim')) or ('Premix' in step) or ('Sim' not in step and 'Gen' in step): diff --git a/Configuration/PyReleaseValidation/python/relval_upgrade.py b/Configuration/PyReleaseValidation/python/relval_upgrade.py index 57907606f90b4..330c7ae4f0729 100644 --- a/Configuration/PyReleaseValidation/python/relval_upgrade.py +++ b/Configuration/PyReleaseValidation/python/relval_upgrade.py @@ -45,6 +45,8 @@ def notForGenOnly(key,specialType): step = 'GenSimHLBeamSpotCloseBy' elif 'CloseBy' in frag or 'CE_E' in frag or 'CE_H' in frag: step = 'GenSimCloseBy' + elif 'Displaced' in frag: + step = 'GenSimDisplaced' stepMaker = makeStepNameSim elif 'Gen' in step: if 'HLBeamSpot' in step: diff --git a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py index a8dae80702034..74d05ee3fb422 100644 --- a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py +++ b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py @@ -224,6 +224,7 @@ def condition(self, fragment, stepList, key, hasHarvest): 'Sim', 'GenSim', 'GenSimCloseBy', + 'GenSimDisplaced', 'GenSimHLBeamSpot', 'GenSimHLBeamSpot14', 'GenSimHLBeamSpotCloseBy', @@ -4081,4 +4082,5 @@ def __init__(self, howMuch, dataset): ('Hydjet_Quenched_MinBias_5519GeV_cfi', UpgradeFragment(U2000by1,'HydjetQMinBias_5519GeV')), ('SingleMuPt15Eta0_0p4_cfi', UpgradeFragment(Kby(9,100),'SingleMuPt15Eta0p_0p4')), ('CloseByPGun_Barrel_Front_cfi', UpgradeFragment(Kby(9,100),'CloseByPGun_Barrel_Front')), + ('DisplacedParticleGun_cfi', UpgradeFragment(Kby(9,100),'DisplacedPGun')), ]) diff --git a/IOMC/ParticleGuns/src/DisplacedParticleGunProducer.cc b/IOMC/ParticleGuns/src/DisplacedParticleGunProducer.cc new file mode 100644 index 0000000000000..4ed916b7c6c91 --- /dev/null +++ b/IOMC/ParticleGuns/src/DisplacedParticleGunProducer.cc @@ -0,0 +1,420 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include "HepMC/GenEvent.h" + +#include "FWCore/AbstractServices/interface/RandomNumberGenerator.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/Run.h" +#include "FWCore/Framework/interface/one/EDProducer.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h" +#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h" +#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h" + +namespace edm { + + namespace { + + // --- Hard-coded HGCAL front surface --- + constexpr double kHGCalZ = 296.50; + constexpr double kHGCalRMin = 26.00; + constexpr double kHGCalRMax = 124.37; + + // Sample r uniformly in area between [rmin, rmax] (uniform point density) + double shootUniformDensity(CLHEP::HepRandomEngine* eng, double rmin, double rmax) { + const double r2 = CLHEP::RandFlat::shoot(eng, rmin * rmin, rmax * rmax); + return std::sqrt(r2); + } + + // Sample r uniformly between [rmin, rmax] + double shootUniformR(CLHEP::HepRandomEngine* eng, double rmin, double rmax) { + return CLHEP::RandFlat::shoot(eng, rmin, rmax); + } + + // Ensure the particle hits HGCAL's front surface within [rMin; rMax] + bool hitsZPlaneWithinR(double x0, + double y0, + double z0, // vertex coordinates + double px, + double py, + double pz, // particle momentum + double zPlane, + double rMin, + double rMax) { // hard-coded HGCAL geometry + // Must move towards the plane: require (zPlane - z0) / pz > 0 + const double dz = zPlane - z0; + if (pz < 1e-6) { + return false; + } + const double t = dz / pz; + if (t <= 0.0) { + return false; + } + + const double xHit = x0 + t * px; + const double yHit = y0 + t * py; + const double rHit = std::hypot(xHit, yHit); + + return (rHit >= rMin && rHit <= rMax); + } + + // Compute the theta range (in radians) that intersects a plane at z=zFront, + // within radial bounds [rMin, rMax], from a vertex at transverse radius R0 and z0. + std::pair thetaRangeToPointToHGCALSurface( + double R, double z, double zHGCAL, double rminHGCAL, double rmaxHGCAL) { + const double dz = zHGCAL - z; + if (dz <= 0. or rminHGCAL > rmaxHGCAL) { + return {0., 0.}; + } + + double thetaMin = std::atan((rminHGCAL - R) / dz); + double thetaMax = std::atan((rmaxHGCAL - R) / dz); + + // For +z endcap, theta should be in (0, pi/2) + // clamping removes instabilities for pz + thetaMin = std::clamp(thetaMin, 1e-6, 0.5 * std::numbers::pi - 1e-6); + thetaMax = std::clamp(thetaMax, 1e-6, 0.5 * std::numbers::pi - 1e-6); + + if (thetaMax <= thetaMin) { + return {0., 0.}; + } + return {thetaMin, thetaMax}; + } + + std::tuple computeMomentum(double pt, double theta, double phi) { + double px = pt * std::cos(phi); + double py = pt * std::sin(phi); + double pz = 0.; + if (theta < 1e-6) { + throw cms::Exception("DisplacedParticleGunProducer") << "Theta is too small. Unstable pz."; + } else { + pz = pt / std::tan(theta); + } + return {px, py, pz}; + } + + } // namespace + + class DisplacedParticleGunProducer : public one::EDProducer { + public: + explicit DisplacedParticleGunProducer(const ParameterSet&); + ~DisplacedParticleGunProducer() override = default; + + void beginRun(const edm::Run& r, const edm::EventSetup&) override; + void endRun(edm::Run const& r, const edm::EventSetup&) override; + void endRunProduce(edm::Run& r, const edm::EventSetup&) override; + + static void fillDescriptions(ConfigurationDescriptions& descriptions); + + private: + void produce(Event& e, const EventSetup& es) override; + + double fPtMin = 0.; + double fPtMax = 0.; + double fPhiMin = 0.; + double fPhiMax = 0.; + double fRMin = 0.; + double fRMax = 0.; + double fPhiVtxMin = 0.; + double fPhiVtxMax = 0.; + double fZVtx = 0.; + int fNParticles = 1; + std::vector fPartIDs; + bool fUniformDensityInR = false; + unsigned int fMaxTries = 1000; + + // If true: derive theta range from hard-coded HGCAL front surface envelope + // (with R in [RminFrontSurfaceHGCAL, RmaxFrontSurfaceHGCAL]) and vertex rho + // If false: sample theta uniformly in [MinTheta, MaxTheta] + bool fPointingToHGCAL = true; + std::optional fRminFrontSurfaceHGCAL = kHGCalRMin; + std::optional fRmaxFrontSurfaceHGCAL = kHGCalRMax; + std::optional fThetaMin = 0.; + std::optional fThetaMax = 0.; + + ESHandle fPDGTable; + const ESGetToken fPDGTableToken; + HepMC::GenEvent* fEvt; + int fVerbosity; + }; + + DisplacedParticleGunProducer::DisplacedParticleGunProducer(const ParameterSet& pset) + : fPDGTableToken(esConsumes()), fEvt(nullptr) { + Service rng; + if (!rng.isAvailable()) { + throw cms::Exception("Configuration") + << "The RandomNumberProducer module requires the RandomNumberGeneratorService\n" + "which appears to be absent. Please add that service to your configuration\n" + "or remove the modules that require it."; + } + + const auto pgun = pset.getParameter("PGunParameters"); + + fPtMin = pgun.getParameter("MinPt"); + fPtMax = pgun.getParameter("MaxPt"); + fPhiMin = pgun.getParameter("MinPhi"); + fPhiMax = pgun.getParameter("MaxPhi"); + fPhiVtxMin = pgun.getParameter("MinVtxPhi"); + fPhiVtxMax = pgun.getParameter("MaxVtxPhi"); + fRMin = pgun.getParameter("RMin"); + fRMax = pgun.getParameter("RMax"); + fZVtx = pgun.getParameter("ZVtx"); + fNParticles = pgun.getParameter("NParticles"); + fPartIDs = pgun.getParameter>("PartID"); + fUniformDensityInR = pgun.getParameter("UniformDensityInR"); + fMaxTries = pgun.getParameter("MaxTries"); + fPointingToHGCAL = pgun.getParameter("PointingToHGCAL"); + + if (!fPointingToHGCAL) { + fThetaMin = pgun.getParameter("MinTheta"); + fThetaMax = pgun.getParameter("MaxTheta"); + if (*fThetaMax <= *fThetaMin) { + throw cms::Exception("DisplacedParticleGunProducer") << "Please ensure MinTheta <= MaxTheta."; + } + fRminFrontSurfaceHGCAL = std::nullopt; + fRmaxFrontSurfaceHGCAL = std::nullopt; + } else { + fRminFrontSurfaceHGCAL = pgun.getParameter("RminFrontSurfaceHGCAL"); + fRmaxFrontSurfaceHGCAL = pgun.getParameter("RmaxFrontSurfaceHGCAL"); + if (*fRmaxFrontSurfaceHGCAL <= *fRminFrontSurfaceHGCAL) { + throw cms::Exception("DisplacedParticleGunProducer") + << "Please ensure RmaxFrontSurfaceHGCAL > RminFrontSurfaceHGCAL."; + } + if (*fRmaxFrontSurfaceHGCAL > kHGCalRMax || *fRminFrontSurfaceHGCAL < kHGCalRMin) { + throw cms::Exception("DisplacedParticleGunProducer") + << "Please ensure RmaxFrontSurfaceHGCAL <= kHGCalRMax and RminFrontSurfaceHGCAL >= kHGCalRMin."; + } + fThetaMin = std::nullopt; + fThetaMax = std::nullopt; + } + + if (fPtMax <= fPtMin) { + throw cms::Exception("DisplacedParticleGunProducer") << "Please fix MinPt/MaxPt"; + } + if (fPhiMax <= fPhiMin) { + throw cms::Exception("DisplacedParticleGunProducer") << "Please fix MinPhi/MaxPhi"; + } + if (fPhiVtxMax <= fPhiVtxMin) { + throw cms::Exception("DisplacedParticleGunProducer") << "Please fix MinVtxPhi/MaxVtxPhi"; + } + if (fRMax <= fRMin) { + throw cms::Exception("DisplacedParticleGunProducer") << "Please fix RMin/RMax"; + } + if (fPartIDs.empty()) { + throw cms::Exception("DisplacedParticleGunProducer") << "PartID must be non-empty"; + } + if (fMaxTries == 0) { + throw cms::Exception("DisplacedParticleGunProducer") << "MaxTries must be > 0"; + } + + produces("unsmeared"); + produces(); + } + + void DisplacedParticleGunProducer::beginRun(const edm::Run& r, const EventSetup& es) { + fPDGTable = es.getHandle(fPDGTableToken); + return; + } + + void DisplacedParticleGunProducer::endRun(const Run& run, const EventSetup& es) {} + + void DisplacedParticleGunProducer::endRunProduce(Run& run, const EventSetup& es) {} + + void DisplacedParticleGunProducer::fillDescriptions(ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("AddAntiParticle", false); + + edm::ParameterSetDescription pgun; + + // particle direction + pgun.add("MinPt", 5.); + pgun.add("MaxPt", 100.); + pgun.add("MinPhi", -std::numbers::pi); + pgun.add("MaxPhi", +std::numbers::pi); + + // vertex displacement (cm) + pgun.add("RMin", 0.); + pgun.add("RMax", 10.); + pgun.add("MinVtxPhi", 0.); + pgun.add("MaxVtxPhi", 2 * std::numbers::pi); + pgun.add("ZVtx", 0.); + + pgun.add("NParticles", 1); + pgun.add>("PartID", {22}); + + pgun.add("UniformDensityInR", false); + + pgun.add("MaxTries", 1000u); + + // A particle shot at the extremities of the HGCAL surface will not traverse a substantial fraction of the detector. + // We use RminFrontSurfaceHGCAL and RmaxFrontSurfaceHGCAL to expose only a given region of the HGCAL surface + // THe default arguments correspond to the inner third of HGCAL's surface (R in ~[58.79, 91.58]cm). + // At (R=200,z=0)cm, an uncharged particle pointing to R=58.79cm at the HGCAL surface (the most extreme case) exits the calorimeter at the + // back face of the CE-E, crossing all its layers. + // For R>200cm there is no guarantee all CE-E layers will be crossed, so a tighter R range might be needed. + // For z>0cm the angles will become more extreme, so a tighter R range might be needed. + // The above reasoning breaks for charged particles, since the bending under the magnetic filed can enormously extend the particle's reach. + pgun.add("PointingToHGCAL", true); + pgun.addOptionalNode(edm::ParameterDescription("RminFrontSurfaceHGCAL", 58.79, true), true); + pgun.addOptionalNode(edm::ParameterDescription("RmaxFrontSurfaceHGCAL", 91.58, true), true); + pgun.addOptionalNode(edm::ParameterDescription("MinTheta", 0.2, true), true); + pgun.addOptionalNode(edm::ParameterDescription("MaxTheta", 1.2, true), true); + + desc.add("PGunParameters", pgun); + + desc.addUntracked("Verbosity", 0); + desc.addUntracked("firstRun", 1); + desc.add("psethack", + "displaced gun with theta, optionally pointing to hard-coded HGCAL front surface"); + + descriptions.add("DisplacedParticleGunProducer", desc); + } + + void DisplacedParticleGunProducer::produce(Event& e, const EventSetup& /*es*/) { + edm::Service rng; + CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID()); + + if (fVerbosity > 0) { + LogDebug("DisplacedParticleGunProducer") + << " DisplacedParticleGunProducer : Begin New Event Generation" << std::endl; + } + + fEvt = new HepMC::GenEvent(); + + const double zFront = kHGCalZ; + const double rMinFace = fRminFrontSurfaceHGCAL.value(); + const double rMaxFace = fRmaxFrontSurfaceHGCAL.value(); + + if (fPointingToHGCAL) { + if (!(rMaxFace > rMinFace) || zFront <= fZVtx) { + throw cms::Exception("DisplacedParticleGunProducer") + << "Invalid hard-coded HGCAL surface envelope: " + << "zFront =" << zFront << "cm, rMin =" << rMinFace << "cm, rMax=" << rMaxFace << "cm" + << " (check ZVtx)."; + } + } + + int barcode = 1; + + for (int ip = 0; ip < fNParticles; ++ip) { + // --- Sample displaced vertex in transverse annulus (z fixed) --- + const double RVtx = + fUniformDensityInR ? shootUniformDensity(engine, fRMin, fRMax) : shootUniformR(engine, fRMin, fRMax); + const double phiVtx = CLHEP::RandFlat::shoot(engine, fPhiVtxMin, fPhiVtxMax); + const double xVtx = RVtx * std::cos(phiVtx); + const double yVtx = RVtx * std::sin(phiVtx); + const double zVtx = fZVtx; + + const auto idx = static_cast(CLHEP::RandFlat::shoot(engine, 0, fPartIDs.size())); + const int partID = fPartIDs[idx]; + + const HepPDT::ParticleData* pData = fPDGTable->particle(HepPDT::ParticleID(std::abs(partID))); + if (!pData) { + throw cms::Exception("DisplacedParticleGunProducer") << "Particle ID " << partID << " not found in PDG table"; + } + const double mass = pData->mass().value(); + + if (pData->charge() != 0 && fPointingToHGCAL) { + throw cms::Exception("DisplacedParticleGunProducer") + << "The logic that points particles to HGCAL's front face assumes that particles move in straight lines."; + } + + double theta = 0., phi = phiVtx, px = 0., py = 0., pz = 0.; + const double pt = CLHEP::RandFlat::shoot(engine, fPtMin, fPtMax); + if (fPointingToHGCAL) { + const double RVtx0 = std::hypot(xVtx, yVtx); + auto [thetaMin, thetaMax] = thetaRangeToPointToHGCALSurface(RVtx0, zVtx, zFront, rMinFace, rMaxFace); + if (!(thetaMax >= thetaMin)) { + throw cms::Exception("DisplacedParticleGunProducer") + << "No valid theta window for this vertex at RVtx=" << RVtx0 << "cm within HGCAL's front surface."; + } + + bool accepted = false; + for (unsigned int itry = 0; itry < fMaxTries; ++itry) { + theta = CLHEP::RandFlat::shoot(engine, thetaMin, thetaMax); + phi = CLHEP::RandFlat::shoot(engine, fPhiMin, fPhiMax); + std::tie(px, py, pz) = computeMomentum(pt, theta, phi); + if (std::isnan(pz) || std::isinf(pz) || pz <= 0.0) { + continue; // must go towards +z plane + } + + if (hitsZPlaneWithinR( + xVtx, yVtx, zVtx, px, py, pz, zFront, *fRminFrontSurfaceHGCAL, *fRmaxFrontSurfaceHGCAL)) { + accepted = true; + break; + } + } + if (!accepted) { + throw cms::Exception("DisplacedParticleGunProducer") + << "Failed to generate a particle intersecting HGCAL front surface after MaxTries=" << fMaxTries + << ". Vertex: (R=" << RVtx << "cm, phiVtx=" << phiVtx << ", z=" << zVtx << "cm). HGCAL band: [" + << *fRminFrontSurfaceHGCAL << "," << *fRmaxFrontSurfaceHGCAL << "]cm at z=" << zFront << "cm."; + } + } else { // if (fPointingToHGCAL) + theta = CLHEP::RandFlat::shoot(engine, fThetaMin.value(), fThetaMax.value()); + phi = CLHEP::RandFlat::shoot(engine, fPhiMin, fPhiMax); + std::tie(px, py, pz) = computeMomentum(pt, theta, phi); + } + + const double p2 = px * px + py * py + pz * pz; + const double energy = std::sqrt(p2 + mass * mass); + + HepMC::FourVector p(px, py, pz, energy); + + HepMC::GenVertex* vtx = + new HepMC::GenVertex(HepMC::FourVector(xVtx * CLHEP::cm, yVtx * CLHEP::cm, zVtx * CLHEP::cm, 0.0)); + + HepMC::GenParticle* part = new HepMC::GenParticle(p, partID, 1); + part->suggest_barcode(barcode++); + + vtx->add_particle_out(part); + fEvt->add_vertex(vtx); + + if (fVerbosity > 0) { + vtx->print(); + part->print(); + } + } + + fEvt->set_event_number(e.id().event()); + fEvt->set_signal_process_id(20); + + if (fVerbosity > 0) { + fEvt->print(); + } + + auto bProduct = std::make_unique(); + bProduct->addHepMCData(fEvt); + e.put(std::move(bProduct), "unsmeared"); + + auto genEventInfo = std::make_unique(fEvt); + e.put(std::move(genEventInfo)); + + if (fVerbosity > 0) { + LogDebug("DisplacedParticleGunProducer") << " DisplacedParticleGunProducer : Event Generation Done " << std::endl; + } + } + +} // namespace edm + +#include "FWCore/Framework/interface/MakerMacros.h" +using edm::DisplacedParticleGunProducer; +DEFINE_FWK_MODULE(DisplacedParticleGunProducer);