ParticleGuns/src/FlatRandomPtAndDxyGunProducer.cc

0001 #include <ostream>
0002
0003 #include "IOMC/ParticleGuns/interface/FlatRandomPtAndDxyGunProducer.h"
0004
0005 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0006 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0007
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "FWCore/ServiceRegistry/interface/Service.h"
0011 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0012
0013 #include "CLHEP/Random/RandFlat.h"
0014
0015 using namespace edm;
0016 using namespace std;
0017
0018 FlatRandomPtAndDxyGunProducer::FlatRandomPtAndDxyGunProducer(const ParameterSet& pset) : BaseFlatGunProducer(pset) {
0019   ParameterSet defpset;
0020   ParameterSet pgun_params = pset.getParameter<ParameterSet>("PGunParameters");
0021
0022   fMinPt = pgun_params.getParameter<double>("MinPt");
0023   fMaxPt = pgun_params.getParameter<double>("MaxPt");
0024   dxyMin_ = pgun_params.getParameter<double>("dxyMin");
0025   dxyMax_ = pgun_params.getParameter<double>("dxyMax");
0026   lxyMax_ = pgun_params.getParameter<double>("LxyMax");
0027   lzMax_ = pgun_params.getParameter<double>("LzMax");
0028   ConeRadius_ = pgun_params.getParameter<double>("ConeRadius");
0029   ConeH_ = pgun_params.getParameter<double>("ConeH");
0030   DistanceToAPEX_ = pgun_params.getParameter<double>("DistanceToAPEX");
0031
0032   produces<HepMCProduct>("unsmeared");
0033   produces<GenEventInfoProduct>();
0034 }
0035
0036 FlatRandomPtAndDxyGunProducer::~FlatRandomPtAndDxyGunProducer() {
0037   // no need to cleanup GenEvent memory - done in HepMCProduct
0038 }
0039
0040 void FlatRandomPtAndDxyGunProducer::produce(Event& e, const EventSetup& es) {
0041   edm::Service<edm::RandomNumberGenerator> rng;
0042   CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
0043
0044   if (fVerbosity > 0) {
0045     cout << " FlatRandomPtAndDxyGunProducer : Begin New Event Generation" << endl;
0046   }
0047   // event loop (well, another step in it...)
0048
0049   // no need to clean up GenEvent memory - done in HepMCProduct
0050   //
0051
0052   // here re-create fEvt (memory)
0053   //
0054   fEvt = new HepMC::GenEvent();
0055
0056   // now actualy, cook up the event from PDGTable and gun parameters
0057   int barcode = 1;
0058   for (unsigned int ip = 0; ip < fPartIDs.size(); ++ip) {
0059     double phi_vtx = 0;
0060     double dxy = 0;
0061     double pt = 0;
0062     double eta = 0;
0063     double px = 0;
0064     double py = 0;
0065     double pz = 0;
0066     double vx = 0;
0067     double vy = 0;
0068     double vz = 0;
0069     double lxy = 0;
0070
0071     bool passLoop = false;
0072     while (not passLoop) {
0073       bool passLxy = false;
0074       bool passLz = false;
0075       phi_vtx = CLHEP::RandFlat::shoot(engine, fMinPhi, fMaxPhi);
0076       dxy = CLHEP::RandFlat::shoot(engine, dxyMin_, dxyMax_);
0077       float dxysign = CLHEP::RandFlat::shoot(engine, -1, 1);
0078       if (dxysign < 0)
0079         dxy = -dxy;
0080       pt = CLHEP::RandFlat::shoot(engine, fMinPt, fMaxPt);
0081       px = pt * cos(phi_vtx);
0082       py = pt * sin(phi_vtx);
0083       for (int i = 0; i < 10000; i++) {
0084         vx = CLHEP::RandFlat::shoot(engine, -lxyMax_, lxyMax_);
0085         vy = (pt * dxy + vx * py) / px;
0086         lxy = sqrt(vx * vx + vy * vy);
0087         if (lxy < abs(lxyMax_) and (vx * px + vy * py) > 0) {
0088           passLxy = true;
0089           break;
0090         }
0091       }
0092       eta = CLHEP::RandFlat::shoot(engine, fMinEta, fMaxEta);
0093       pz = pt * sinh(eta);
0094       //vz = fabs(fRandomGaussGenerator->fire(0.0, LzWidth_/2.0));
0095       float ConeTheta = ConeRadius_ / ConeH_;
0096       for (int j = 0; j < 100; j++) {
0097         vz = CLHEP::RandFlat::shoot(engine, 0.0, lzMax_);  // this is abs(vz)
0098         float v0 = vz - DistanceToAPEX_;
0099         if (v0 <= 0 or lxy * lxy / (ConeTheta * ConeTheta) > v0 * v0) {
0100           passLz = true;
0101           break;
0102         }
0103       }
0104       if (pz < 0)
0105         vz = -vz;
0106       passLoop = (passLxy and passLz);
0107
0108       if (passLoop)
0109         break;
0110     }
0111
0112     float time = sqrt(vx * vx + vy * vy + vz * vz);
0113
0114     HepMC::GenVertex* Vtx1 = new HepMC::GenVertex(HepMC::FourVector(vx, vy, vz, time));
0115
0116     int PartID = fPartIDs[ip];
0117     const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID)));
0118     double mass = PData->mass().value();
0119     double energy2 = px * px + py * py + pz * pz + mass * mass;
0120     double energy = sqrt(energy2);
0121     HepMC::FourVector p(px, py, pz, energy);
0122     HepMC::GenParticle* Part = new HepMC::GenParticle(p, PartID, 1);
0123     Part->suggest_barcode(barcode);
0124     barcode++;
0125     Vtx1->add_particle_out(Part);
0126     fEvt->add_vertex(Vtx1);
0127
0128     if (fAddAntiParticle) {
0129       HepMC::GenVertex* Vtx2 = new HepMC::GenVertex(HepMC::FourVector(-vx, -vy, -vz, time));
0130       HepMC::FourVector ap(-px, -py, -pz, energy);
0131       int APartID = -PartID;
0132       if (PartID == 22 || PartID == 23) {
0133         APartID = PartID;
0134       }
0135       HepMC::GenParticle* APart = new HepMC::GenParticle(ap, APartID, 1);
0136       APart->suggest_barcode(barcode);
0137       barcode++;
0138       Vtx2->add_particle_out(APart);
0139       fEvt->add_vertex(Vtx2);
0140     }
0141   }
0142   fEvt->set_event_number(e.id().event());
0143   fEvt->set_signal_process_id(20);
0144
0145   if (fVerbosity > 0) {
0146     fEvt->print();
0147   }
0148
0149   unique_ptr<HepMCProduct> BProduct(new HepMCProduct());
0150   BProduct->addHepMCData(fEvt);
0151   e.put(std::move(BProduct), "unsmeared");
0152
0153   unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
0154   e.put(std::move(genEventInfo));
0155
0156   if (fVerbosity > 0) {
0157     cout << " FlatRandomPtAndDxyGunProducer : End New Event Generation" << endl;
0158     fEvt->print();
0159   }
0160 }