CMS 3D CMS Logo

CloseByParticleGunProducer.cc
Go to the documentation of this file.
1 #include <ostream>
2 
4 
7 
9 
15 
16 #include "CLHEP/Random/RandFlat.h"
17 #include "CLHEP/Units/GlobalSystemOfUnits.h"
18 #include "CLHEP/Units/GlobalPhysicalConstants.h"
19 #include "CLHEP/Random/RandFlat.h"
20 
21 using namespace edm;
22 using namespace std;
23 
26 {
27 
28  ParameterSet defpset ;
29  ParameterSet pgun_params =
30  pset.getParameter<ParameterSet>("PGunParameters") ;
31 
32  fEnMax = pgun_params.getParameter<double>("EnMax");
33  fEnMin = pgun_params.getParameter<double>("EnMin");
34  fRMax = pgun_params.getParameter<double>("RMax");
35  fRMin = pgun_params.getParameter<double>("RMin");
36  fZMax = pgun_params.getParameter<double>("ZMax");
37  fZMin = pgun_params.getParameter<double>("ZMin");
38  fDelta = pgun_params.getParameter<double>("Delta");
39  fPhiMin = pgun_params.getParameter<double>("MinPhi");
40  fPhiMax = pgun_params.getParameter<double>("MaxPhi");
41  fPointing = pgun_params.getParameter<bool>("Pointing");
42  fOverlapping = pgun_params.getParameter<bool>("Overlapping");
43  fRandomShoot = pgun_params.getParameter<bool>("RandomShoot");
44  fNParticles = pgun_params.getParameter<int>("NParticles");
45  fPartIDs = pgun_params.getParameter< vector<int> >("PartID");
46 
47  produces<HepMCProduct>("unsmeared");
48  produces<GenEventInfoProduct>();
49 
50 }
51 
53 {
54  // no need to cleanup GenEvent memory - done in HepMCProduct
55 }
56 
58 {
60  CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
61 
62  if ( fVerbosity > 0 )
63  {
64  LogDebug("CloseByParticleGunProducer") << " CloseByParticleGunProducer : Begin New Event Generation" << endl ;
65  }
66  fEvt = new HepMC::GenEvent() ;
67 
68  // loop over particles
69  //
70  int barcode = 1 ;
71  int numParticles = fRandomShoot ? CLHEP::RandFlat::shoot(engine, 1, fNParticles) : fNParticles;
72  std::vector<int> particles;
73 
74  for(int i=0; i<numParticles; i++){
75  int partIdx = CLHEP::RandFlat::shoot(engine, 0, fPartIDs.size());
76  particles.push_back(fPartIDs[partIdx]);
77  }
78 
79  double phi = CLHEP::RandFlat::shoot(engine, fPhiMin, fPhiMax);
80  double fR = CLHEP::RandFlat::shoot(engine,fRMin,fRMax);
81  double fZ = CLHEP::RandFlat::shoot(engine,fZMin,fZMax);
82  double tmpPhi = phi;
83  double tmpR = fR;
84 
85  for (unsigned int ip=0; ip<particles.size(); ++ip)
86  {
87  if(fOverlapping)
88  {
89  fR = CLHEP::RandFlat::shoot(engine,tmpR-fDelta,tmpR+fDelta);
90  phi = CLHEP::RandFlat::shoot(engine, tmpPhi-fDelta/fR, tmpPhi+fDelta/fR);
91  }
92  else
93  phi += fDelta/fR;
94 
95  double fEn = CLHEP::RandFlat::shoot(engine,fEnMin,fEnMax);
96  int PartID = particles[ip] ;
97  const HepPDT::ParticleData *PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID))) ;
98  double mass = PData->mass().value() ;
99  double mom = sqrt(fEn*fEn-mass*mass);
100  double px = 0.;
101  double py = 0.;
102  double pz = mom;
103  double energy = fEn;
104 
105  // Compute Vertex Position
106  double x=fR*cos(phi);
107  double y=fR*sin(phi);
108  constexpr double c= 2.99792458e+1; // cm/ns
109  double timeOffset = sqrt(x*x + y*y + fZ*fZ)/c*ns*c_light;
110  HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(x*cm,y*cm,fZ*cm,timeOffset));
111 
112  HepMC::FourVector p(px,py,pz,energy) ;
113  // If we are requested to be pointing to (0,0,0), correct the momentum direction
114  if (fPointing) {
115  math::XYZVector direction(x,y,fZ);
116  math::XYZVector momentum = direction.unit() * mom;
117  p.setX(momentum.x());
118  p.setY(momentum.y());
119  p.setZ(momentum.z());
120  }
121  HepMC::GenParticle* Part = new HepMC::GenParticle(p,PartID,1);
122  Part->suggest_barcode( barcode );
123  barcode++;
124 
125  Vtx->add_particle_out(Part);
126 
127  if (fVerbosity > 0) {
128  Vtx->print();
129  Part->print();
130  }
131  fEvt->add_vertex(Vtx);
132  }
133 
134 
135  fEvt->set_event_number(e.id().event());
136  fEvt->set_signal_process_id(20);
137 
138  if ( fVerbosity > 0 )
139  {
140  fEvt->print();
141  }
142 
143  unique_ptr<HepMCProduct> BProduct(new HepMCProduct());
144  BProduct->addHepMCData( fEvt );
145  e.put(std::move(BProduct), "unsmeared");
146 
147  unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
148  e.put(std::move(genEventInfo));
149 
150  if ( fVerbosity > 0 )
151  {
152  LogDebug("CloseByParticleGunProducer") << " CloseByParticleGunProducer : Event Generation Done " << endl;
153  }
154 
155  particles.clear();
156 }
#define LogDebug(id)
T getParameter(std::string const &) const
EventNumber_t event() const
Definition: EventID.h:41
static unsigned int partIdx(const InputGenJetsParticleSelector::ParticleVector &p, const reco::Candidate *particle)
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:125
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
virtual CLHEP::HepRandomEngine & getEngine(StreamID const &)=0
Use this engine in event methods.
ESHandle< HepPDT::ParticleDataTable > fPDGTable
T sqrt(T t)
Definition: SSEVec.h:18
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
void produce(Event &e, const EventSetup &es) override
HepPDT::ParticleData ParticleData
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:30
edm::EventID id() const
Definition: EventBase.h:59
HLT enums.
StreamID streamID() const
Definition: Event.h:95
def move(src, dest)
Definition: eostools.py:511
#define constexpr