Add new displaced muon gun

Configuration/Generator/python/FlatRandomPtAndDxyGunProducer_cfi.py

@@ -0,0 +1,24 @@
+import FWCore.ParameterSet.Config as cms
+
+generator = cms.EDProducer("FlatRandomPtAndDxyGunProducer",
+    PGunParameters = cms.PSet(
+        PartID = cms.vint32(-13),
+        MinPt  = cms.double(2.00),
+        MaxPt  = cms.double(50.00),
+        MinEta = cms.double(-3.0),
+        MaxEta = cms.double(3.0),
+        MinPhi = cms.double(-3.14159265359),
+        MaxPhi = cms.double(3.14159265359),
+        LxyMax = cms.double(3000.0),#make sure most muons generated before Muon system, Gauss distribution
+        LzMax = cms.double(5000.0),#make sure most muons generated before Muon system, Gauss distribution
+	ConeRadius = cms.double(1000.0),
+	ConeH = cms.double(3000.0),
+	DistanceToAPEX = cms.double(2000.0),
+        dxyMin = cms.double(0.0),
+        dxyMax = cms.double(500.0)
+    ),
+    Verbosity = cms.untracked.int32(0), ## set to 1 (or greater)  for printouts
+    psethack = cms.string('displaced muon'),
+    AddAntiParticle = cms.bool(True),
+    firstRun = cms.untracked.uint32(1)
+)

IOMC/ParticleGuns/interface/FlatRandomPtAndDxyGunProducer.h

@@ -0,0 +1,50 @@
+#ifndef FlatRandomPtAndDxyGunProducer_H
+#define FlatRandomPtAndDxyGunProducer_H
+
+/** \class FlatRandomPtAndDxyGunProducer
+ *
+ * Generates single particle gun in HepMC format
+ * The d0 is taken by convention equal to dxy
+ * v1: first version of displaced Muon Gun
+ * flat in dxy(0, 50), flat in lz(-50, 50)
+ * v2: bugfixed, unit in HepMC is mm, rather than cm
+ * v3: check pz*vz>0, make sure muon is not flying 
+ * back to beamspot, Lz is generated by Gauss distribution
+ * v4: dxy is produced with sign, (-500, 500) [mm]
+ * v5: check vx*px+vy*py>0, and muons inside beamline-symmetric 
+ * cone is vetoed by lxy^2/c^2<=(vz-d) ^2 && vz-d>0
+
+ * Contact Sven Dildick, Tao Huang
+ ***************************************/
+
+#include "IOMC/ParticleGuns/interface/BaseFlatGunProducer.h"
+
+namespace edm
+{
+
+    class FlatRandomPtAndDxyGunProducer : public BaseFlatGunProducer
+    {
+
+    public:
+        FlatRandomPtAndDxyGunProducer(const ParameterSet & pset);
+        virtual ~FlatRandomPtAndDxyGunProducer();
+
+        virtual void produce(Event & e, const EventSetup& es) override;
+
+    private:
+
+        // data members
+
+        double            fMinPt   ;
+        double            fMaxPt   ;
+        double            dxyMin_   ;
+        double            dxyMax_   ;
+        double            lxyMax_   ;
+	double 		  lzMax_   ;
+        double            ConeRadius_ ;
+	double 		  ConeH_    ;
+	double	 	  DistanceToAPEX_; 
+    };
+}
+
+#endif

IOMC/ParticleGuns/src/FlatRandomPtAndDxyGunProducer.cc

@@ -0,0 +1,168 @@
+#include <ostream>
+
+#include "IOMC/ParticleGuns/interface/FlatRandomPtAndDxyGunProducer.h"
+
+#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
+#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.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 "CLHEP/Random/RandFlat.h"
+
+using namespace edm;
+using namespace std;
+
+FlatRandomPtAndDxyGunProducer::FlatRandomPtAndDxyGunProducer(const ParameterSet& pset) :
+BaseFlatGunProducer(pset)
+{
+  ParameterSet defpset ;
+  ParameterSet pgun_params =
+    pset.getParameter<ParameterSet>("PGunParameters") ;
+
+  fMinPt = pgun_params.getParameter<double>("MinPt");
+  fMaxPt = pgun_params.getParameter<double>("MaxPt");
+  dxyMin_ = pgun_params.getParameter<double>("dxyMin");
+  dxyMax_ = pgun_params.getParameter<double>("dxyMax");
+  lxyMax_ = pgun_params.getParameter<double>("LxyMax");
+  lzMax_ = pgun_params.getParameter<double>("LzMax");
+  ConeRadius_ = pgun_params.getParameter<double>("ConeRadius");
+  ConeH_ = pgun_params.getParameter<double>("ConeH");
+  DistanceToAPEX_ = pgun_params.getParameter<double>("DistanceToAPEX");
+
+  produces<HepMCProduct>("unsmeared");
+  produces<GenEventInfoProduct>();
+}
+
+FlatRandomPtAndDxyGunProducer::~FlatRandomPtAndDxyGunProducer()
+{
+    // no need to cleanup GenEvent memory - done in HepMCProduct
+}
+
+void FlatRandomPtAndDxyGunProducer::produce(Event &e, const EventSetup& es)
+{
+   edm::Service<edm::RandomNumberGenerator> rng;
+   CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
+
+  if ( fVerbosity > 0 ){
+    cout << " FlatRandomPtAndDxyGunProducer : Begin New Event Generation" << endl ;
+  }
+  // event loop (well, another step in it...)
+
+  // no need to clean up GenEvent memory - done in HepMCProduct
+  //
+
+  // here re-create fEvt (memory)
+  //
+  fEvt = new HepMC::GenEvent() ;
+
+  // now actualy, cook up the event from PDGTable and gun parameters
+  int barcode = 1 ;
+  for (unsigned int ip=0; ip<fPartIDs.size(); ++ip){
+
+    double phi_vtx = 0;
+    double dxy = 0;
+    double pt = 0;
+    double eta = 0;
+    double px = 0;
+    double py = 0;
+    double pz = 0;
+    double vx = 0;
+    double vy = 0;
+    double vz = 0;
+    double lxy = 0;
+
+    bool passLoop = false;
+    while (not passLoop) {
+
+      bool passLxy = false;
+      bool passLz = false;
+      phi_vtx = CLHEP::RandFlat::shoot(engine, fMinPhi, fMaxPhi);
+      dxy = CLHEP::RandFlat::shoot(engine, dxyMin_, dxyMax_);
+      float dxysign = CLHEP::RandFlat::shoot(engine, -1, 1);
+      if (dxysign < 0)
+	  dxy = -dxy;
+      pt = CLHEP::RandFlat::shoot(engine, fMinPt, fMaxPt);
+      px = pt*cos(phi_vtx);
+      py = pt*sin(phi_vtx);
+      for (int i=0; i<10000; i++){
+        vx = CLHEP::RandFlat::shoot(engine, -lxyMax_, lxyMax_);
+        vy = (pt*dxy + vx * py)/px;
+        lxy = sqrt(vx*vx + vy*vy);
+        if (lxy<abs(lxyMax_) and (vx*px + vy*py)>0){
+          passLxy = true;
+          break;
+        }
+      }
+      eta = CLHEP::RandFlat::shoot(engine, fMinEta, fMaxEta);
+      pz = pt * sinh(eta);
+      //vz = fabs(fRandomGaussGenerator->fire(0.0, LzWidth_/2.0));
+      float ConeTheta = ConeRadius_/ConeH_;
+      for (int j=0; j<100; j++){
+	  vz = CLHEP::RandFlat::shoot(engine, 0.0, lzMax_);// this is abs(vz)
+	  float v0 = vz - DistanceToAPEX_;
+	  if (v0>0 and lxy*lxy/(ConeTheta*ConeTheta)<=(v0*v0))
+	     continue;
+	  else{
+	     passLz = true;
+	     break;
+	  }
+      }
+      if (pz<0)
+	  vz = -vz;
+      passLoop = (passLxy and passLz);
+
+      if (passLoop)
+        break;
+    }
+
+    HepMC::GenVertex* Vtx1  = new HepMC::GenVertex(HepMC::FourVector(vx,vy,vz));
+
+    int PartID = fPartIDs[ip] ;
+    const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID))) ;
+    double mass   = PData->mass().value() ;
+    double energy2= px*px + py*py + pz*pz + mass*mass ;
+    double energy = sqrt(energy2) ;
+    HepMC::FourVector p(px,py,pz,energy) ;
+    HepMC::GenParticle* Part = new HepMC::GenParticle(p,PartID,1);
+    Part->suggest_barcode( barcode ) ;
+    barcode++;
+    Vtx1->add_particle_out(Part);
+    fEvt->add_vertex(Vtx1) ;
+
+    if ( fAddAntiParticle ) {
+      HepMC::GenVertex* Vtx2  = new HepMC::GenVertex(HepMC::FourVector(-vx,-vy,-vz));
+      HepMC::FourVector ap(-px,-py,-pz,energy) ;
+      int APartID = -PartID ;
+      if ( PartID == 22 || PartID == 23 )
+        {
+          APartID = PartID ;
+        }
+      HepMC::GenParticle* APart = new HepMC::GenParticle(ap,APartID,1);
+      APart->suggest_barcode( barcode ) ;
+      barcode++ ;
+      Vtx2->add_particle_out(APart) ;
+      fEvt->add_vertex(Vtx2) ;
+    }
+  }
+  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 ){
+    cout << " FlatRandomPtAndDxyGunProducer : End New Event Generation" << endl ;
+    fEvt->print();
+  }
+}