File indexing completed on 2024-04-06 12:19:02
0001
0002
0003
0004
0005
0006 #include <ostream>
0007
0008 #include "IOMC/ParticleGuns/interface/ExpoRandomPGunProducer.h"
0009
0010 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0011 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0012
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/ServiceRegistry/interface/Service.h"
0016 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0017
0018 #include "CLHEP/Random/RandFlat.h"
0019
0020 using namespace edm;
0021 using namespace std;
0022
0023 ExpoRandomPGunProducer::ExpoRandomPGunProducer(const ParameterSet& pset) : BaseFlatGunProducer(pset) {
0024 ParameterSet defpset;
0025 ParameterSet pgun_params = pset.getParameter<ParameterSet>("PGunParameters");
0026
0027 fMinP = pgun_params.getParameter<double>("MinP");
0028 fMaxP = pgun_params.getParameter<double>("MaxP");
0029
0030 produces<HepMCProduct>("unsmeared");
0031 produces<GenEventInfoProduct>();
0032 }
0033
0034 ExpoRandomPGunProducer::~ExpoRandomPGunProducer() {
0035
0036 }
0037
0038 void ExpoRandomPGunProducer::produce(Event& e, const EventSetup& es) {
0039 edm::Service<edm::RandomNumberGenerator> rng;
0040 CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
0041
0042 if (fVerbosity > 0) {
0043 std::cout << " ExpoRandomPGunProducer : Begin New Event Generation" << std::endl;
0044 }
0045
0046
0047
0048
0049
0050
0051
0052
0053 fEvt = new HepMC::GenEvent();
0054
0055
0056
0057
0058
0059
0060 HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(0., 0., 0.));
0061
0062
0063
0064 int barcode = 1;
0065 for (unsigned int ip = 0; ip < fPartIDs.size(); ++ip) {
0066 double pmom = CLHEP::RandFlat::shoot(engine, fMinP, fMaxP);
0067 double y = (1. / fMinP) * CLHEP::RandFlat::shoot(engine, 0.0, 1.0);
0068 double f = 1. / pmom;
0069 bool accpt = (y < f);
0070
0071 while ((pmom < fMinP || pmom > fMaxP) || !accpt) {
0072 pmom = CLHEP::RandFlat::shoot(engine, fMinP, fMaxP);
0073 y = (1. / fMinP) * CLHEP::RandFlat::shoot(engine, 0.0, 1.0);
0074 f = 1. / pmom;
0075 accpt = (y < f);
0076 }
0077
0078 double eta = CLHEP::RandFlat::shoot(engine, fMinEta, fMaxEta);
0079 double phi = CLHEP::RandFlat::shoot(engine, fMinPhi, fMaxPhi);
0080 int PartID = fPartIDs[ip];
0081 const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID)));
0082 double mass = PData->mass().value();
0083 double theta = 2. * atan(exp(-eta));
0084 double mom = pmom;
0085 double pt = mom * sin(theta);
0086 double px = pt * cos(phi);
0087 double py = pt * sin(phi);
0088 double pz = mom * cos(theta);
0089 double energy2 = mom * mom + mass * mass;
0090 double energy = sqrt(energy2);
0091
0092
0093
0094 HepMC::FourVector p(px, py, pz, energy);
0095 HepMC::GenParticle* Part = new HepMC::GenParticle(p, PartID, 1);
0096 Part->suggest_barcode(barcode);
0097 barcode++;
0098 Vtx->add_particle_out(Part);
0099
0100 if (fAddAntiParticle) {
0101
0102 HepMC::FourVector ap(-px, -py, -pz, energy);
0103 int APartID = -PartID;
0104 if (PartID == 22 || PartID == 23) {
0105 APartID = PartID;
0106 }
0107
0108
0109 HepMC::GenParticle* APart = new HepMC::GenParticle(ap, APartID, 1);
0110 APart->suggest_barcode(barcode);
0111 barcode++;
0112 Vtx->add_particle_out(APart);
0113 }
0114 }
0115
0116 fEvt->add_vertex(Vtx);
0117 fEvt->set_event_number(e.id().event());
0118 fEvt->set_signal_process_id(20);
0119
0120 if (fVerbosity > 0) {
0121 fEvt->print();
0122 }
0123
0124 unique_ptr<HepMCProduct> BProduct(new HepMCProduct());
0125 BProduct->addHepMCData(fEvt);
0126 e.put(std::move(BProduct), "unsmeared");
0127
0128 unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
0129 e.put(std::move(genEventInfo));
0130
0131 if (fVerbosity > 0) {
0132
0133
0134 std::cout << " FlatRandomPGunProducer : Event Generation Done " << std::endl;
0135 }
0136 }