Merge pull request #26065 from rovere/CloseByGunGenerator

Close by gun generator

IOMC/EventVertexGenerators/interface/PassThroughEvtVtxGenerator.h

@@ -0,0 +1,44 @@
+#ifndef IOMC_EventVertexGenerators_PassThroughEvtVtxGenerator_H
+#define IOMC_EventVertexGenerators_PassThroughEvtVtxGenerator_H
+/*
+*/
+
+#include "IOMC/EventVertexGenerators/interface/BaseEvtVtxGenerator.h"
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/Utilities/interface/EDGetToken.h"
+
+#include "TMatrixD.h"
+
+namespace HepMC {
+   class FourVector ;
+}
+
+namespace CLHEP {
+   class HepRandomEngine;
+}
+
+namespace edm {
+   class HepMCProduct;
+}
+
+class PassThroughEvtVtxGenerator : public BaseEvtVtxGenerator
+{
+   public:
+
+   // ctor & dtor
+   explicit PassThroughEvtVtxGenerator( const edm::ParameterSet& );
+   ~PassThroughEvtVtxGenerator() override;
+
+   void produce( edm::Event&, const edm::EventSetup&) override;
+
+   HepMC::FourVector newVertex(CLHEP::HepRandomEngine*) const override;
+
+   TMatrixD const* GetInvLorentzBoost() const override { return nullptr;};
+
+   private :
+
+   edm::EDGetTokenT<edm::HepMCProduct> sourceToken;
+
+};
+
+#endif

IOMC/EventVertexGenerators/src/PassThroughEvtVtxGenerator.cc

@@ -0,0 +1,67 @@
+
+/*
+*/
+
+#include "IOMC/EventVertexGenerators/interface/PassThroughEvtVtxGenerator.h"
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
+#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
+
+#include "FWCore/ServiceRegistry/interface/Service.h"
+#include "FWCore/Utilities/interface/RandomNumberGenerator.h"
+
+#include "FWCore/Utilities/interface/Exception.h"
+
+#include "DataFormats/Provenance/interface/Provenance.h"
+#include "FWCore/Utilities/interface/EDMException.h"
+
+//#include "HepMC/GenEvent.h"
+// #include "CLHEP/Vector/ThreeVector.h"
+// #include "HepMC/SimpleVector.h"
+
+using namespace edm;
+using namespace CLHEP;
+//using namespace HepMC;
+
+
+PassThroughEvtVtxGenerator::PassThroughEvtVtxGenerator( const ParameterSet& pset )
+  : BaseEvtVtxGenerator(pset)
+{
+   Service<RandomNumberGenerator> rng;
+   if ( ! rng.isAvailable()) {
+     throw cms::Exception("Configuration")
+       << "The PassThroughEvtVtxGenerator requires the RandomNumberGeneratorService\n"
+          "which is not present in the configuration file. \n"
+          "You must add the service\n"
+          "in the configuration file or remove the modules that require it.";
+   }
+   sourceToken=consumes<edm::HepMCProduct>(pset.getParameter<edm::InputTag>("src"));
+}
+
+PassThroughEvtVtxGenerator::~PassThroughEvtVtxGenerator()
+{
+}
+
+HepMC::FourVector PassThroughEvtVtxGenerator::newVertex(CLHEP::HepRandomEngine*) const  {
+  return HepMC::FourVector(0.,0.,0.);
+}
+
+void PassThroughEvtVtxGenerator::produce( Event& evt, const EventSetup& )
+{
+   edm::Service<edm::RandomNumberGenerator> rng;
+
+   Handle<HepMCProduct> HepUnsmearedMCEvt ;
+
+   evt.getByToken( sourceToken, HepUnsmearedMCEvt ) ;
+
+   // Copy the HepMC::GenEvent
+   HepMC::GenEvent* genevt = new HepMC::GenEvent(*HepUnsmearedMCEvt->GetEvent());
+   std::unique_ptr<edm::HepMCProduct> HepMCEvt(new edm::HepMCProduct(genevt));
+
+   evt.put(std::move(HepMCEvt)) ;
+
+   return ;
+}

IOMC/EventVertexGenerators/src/module.cc

@@ -3,6 +3,7 @@
 //#include "IOMC/EventVertexGenerators/interface/VertexGenerator.h"
 
 #include "IOMC/EventVertexGenerators/interface/BaseEvtVtxGenerator.h"
+#include "IOMC/EventVertexGenerators/interface/PassThroughEvtVtxGenerator.h"
 #include "IOMC/EventVertexGenerators/interface/GaussEvtVtxGenerator.h"
 #include "IOMC/EventVertexGenerators/interface/FlatEvtVtxGenerator.h"
 #include "IOMC/EventVertexGenerators/interface/BeamProfileVtxGenerator.h"
@@ -18,6 +19,7 @@
 
 //using edm::VertexGenerator;
 //DEFINE_FWK_MODULE(VertexGenerator) ;
+DEFINE_FWK_MODULE(PassThroughEvtVtxGenerator) ;
 DEFINE_FWK_MODULE(GaussEvtVtxGenerator) ;
 DEFINE_FWK_MODULE(FlatEvtVtxGenerator) ;
 DEFINE_FWK_MODULE(BeamProfileVtxGenerator) ;

IOMC/ParticleGuns/interface/CloseByParticleGunProducer.h

@@ -0,0 +1,29 @@
+#ifndef IOMC_ParticleGun_CloseByParticleGunProducer_H
+#define IOMC_ParticleGun_CloseByParticleGunProducer_H
+
+#include "IOMC/ParticleGuns/interface/BaseFlatGunProducer.h"
+
+namespace edm
+{
+
+  class CloseByParticleGunProducer : public BaseFlatGunProducer
+  {
+
+  public:
+    CloseByParticleGunProducer(const ParameterSet &);
+    ~CloseByParticleGunProducer() override;
+
+  private:
+
+    void produce(Event & e, const EventSetup& es) override;
+
+  protected :
+
+    // data members
+    double fEn,fR,fZ,fDelta;
+    bool fPointing = false;
+    std::vector<int> fPartIDs;
+  };
+}
+
+#endif

IOMC/ParticleGuns/src/CloseByParticleGunProducer.cc

@@ -0,0 +1,125 @@
+#include <ostream>
+
+#include "IOMC/ParticleGuns/interface/CloseByParticleGunProducer.h"
+
+#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
+#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
+
+#include "DataFormats/Math/interface/Vector3D.h"
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ServiceRegistry/interface/Service.h"
+#include "FWCore/Utilities/interface/RandomNumberGenerator.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+
+#include "CLHEP/Random/RandFlat.h"
+#include "CLHEP/Units/GlobalSystemOfUnits.h"
+#include "CLHEP/Units/GlobalPhysicalConstants.h"
+#include "CLHEP/Random/RandFlat.h"
+
+using namespace edm;
+using namespace std;
+
+CloseByParticleGunProducer::CloseByParticleGunProducer(const ParameterSet& pset) :
+   BaseFlatGunProducer(pset)
+{
+
+  ParameterSet defpset ;
+  ParameterSet pgun_params =
+    pset.getParameter<ParameterSet>("PGunParameters") ;
+
+  fEn = pgun_params.getParameter<double>("En");
+  fR = pgun_params.getParameter<double>("R");
+  fZ = pgun_params.getParameter<double>("Z");
+  fDelta = pgun_params.getParameter<double>("Delta");
+  fPartIDs = pgun_params.getParameter< vector<int> >("PartID");
+  fPointing = pgun_params.getParameter<bool>("Pointing");
+
+  produces<HepMCProduct>("unsmeared");
+  produces<GenEventInfoProduct>();
+}
+
+CloseByParticleGunProducer::~CloseByParticleGunProducer()
+{
+   // no need to cleanup GenEvent memory - done in HepMCProduct
+}
+
+void CloseByParticleGunProducer::produce(Event &e, const EventSetup& es)
+{
+   edm::Service<edm::RandomNumberGenerator> rng;
+   CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
+
+   if ( fVerbosity > 0 )
+     {
+       LogDebug("CloseByParticleGunProducer") << " CloseByParticleGunProducer : Begin New Event Generation" << endl ;
+     }
+   fEvt = new HepMC::GenEvent() ;
+
+   // loop over particles
+   //
+   int barcode = 1 ;
+   double phi = CLHEP::RandFlat::shoot(engine, -3.14159265358979323846, 3.14159265358979323846);
+   for (unsigned int ip=0; ip<fPartIDs.size(); ++ip, phi += fDelta/fR)
+   {
+
+     int PartID = fPartIDs[ip] ;
+     const HepPDT::ParticleData *PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID))) ;
+     double mass   = PData->mass().value() ;
+     double mom    = sqrt(fEn*fEn-mass*mass);
+     double px     = 0.;
+     double py     = 0.;
+     double pz     = mom;
+     double energy = fEn;
+
+     // Compute Vertex Position
+     double x=fR*cos(phi);
+     double y=fR*sin(phi);
+     constexpr double c= 2.99792458e+1; // cm/ns
+     double timeOffset = sqrt(x*x + y*y + fZ*fZ)/c*ns*c_light;
+     HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(x*cm,y*cm,fZ*cm,timeOffset));
+
+     HepMC::FourVector p(px,py,pz,energy) ;
+     // If we are requested to be pointing to (0,0,0), correct the momentum direction
+     if (fPointing) {
+       math::XYZVector direction(x,y,fZ);
+       math::XYZVector momentum = direction.unit() * mom;
+       p.setX(momentum.x());
+       p.setY(momentum.y());
+       p.setZ(momentum.z());
+     }
+     HepMC::GenParticle* Part = new HepMC::GenParticle(p,PartID,1);
+     Part->suggest_barcode( barcode );
+     barcode++;
+
+     Vtx->add_particle_out(Part);
+
+     if (fVerbosity > 0) {
+       Vtx->print();
+       Part->print();
+     }
+     fEvt->add_vertex(Vtx);
+   }
+
+
+   fEvt->set_event_number(e.id().event());
+   fEvt->set_signal_process_id(20);
+
+   if ( fVerbosity > 0 )
+   {
+      fEvt->print();
+   }
+
+   unique_ptr<HepMCProduct> BProduct(new HepMCProduct());
+   BProduct->addHepMCData( fEvt );
+   e.put(std::move(BProduct), "unsmeared");
+
+   unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
+   e.put(std::move(genEventInfo));
+
+   if ( fVerbosity > 0 )
+     {
+       LogDebug("CloseByParticleGunProducer") << " CloseByParticleGunProducer : Event Generation Done " << endl;
+     }
+}
+

IOMC/ParticleGuns/src/SealModule.cc

@@ -19,12 +19,13 @@
 #include "IOMC/ParticleGuns/interface/RandomtXiGunProducer.h"
 #include "IOMC/ParticleGuns/interface/FlatRandomPtAndDxyGunProducer.h"
 #include "IOMC/ParticleGuns/interface/RandomMultiParticlePGunProducer.h"
+#include "IOMC/ParticleGuns/interface/CloseByParticleGunProducer.h"
 
 
 // particle gun prototypes
 //
-  
-  
+
+
 /*
 using edm::FlatEGunASCIIWriter;
 DEFINE_FWK_MODULE(FlatEGunASCIIWriter);
@@ -58,3 +59,5 @@ using edm::FlatRandomPtAndDxyGunProducer;
 DEFINE_FWK_MODULE(FlatRandomPtAndDxyGunProducer);
 using edm::RandomMultiParticlePGunProducer;
 DEFINE_FWK_MODULE(RandomMultiParticlePGunProducer);
+using edm::CloseByParticleGunProducer;
+DEFINE_FWK_MODULE(CloseByParticleGunProducer);