Pythia6Interface/plugins/Pythia6PtYDistGun.cc

0001
0002 #include <iostream>
0003
0004 #include "Pythia6PtYDistGun.h"
0005 #include "GeneratorInterface/Pythia6Interface/interface/PtYDistributor.h"
0006
0007 #include "FWCore/Utilities/interface/Exception.h"
0008
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010
0011 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0012
0013 #include "FWCore/Framework/interface/MakerMacros.h"
0014
0015 using namespace edm;
0016 using namespace gen;
0017
0018 Pythia6PtYDistGun::Pythia6PtYDistGun(const ParameterSet& pset) : Pythia6ParticleGun(pset), fPtYGenerator(nullptr) {
0019   ParameterSet pgun_params = pset.getParameter<ParameterSet>("PGunParameters");
0020
0021   fPtYGenerator = new PtYDistributor(pgun_params.getParameter<FileInPath>("kinematicsFile"),
0022                                      pgun_params.getParameter<double>("MaxPt"),
0023                                      pgun_params.getParameter<double>("MinPt"),
0024                                      pgun_params.getParameter<double>("MaxY"),
0025                                      pgun_params.getParameter<double>("MinY"),
0026                                      pgun_params.getParameter<int>("PtBinning"),
0027                                      pgun_params.getParameter<int>("YBinning"));
0028 }
0029
0030 Pythia6PtYDistGun::~Pythia6PtYDistGun() {
0031   if (fPtYGenerator)
0032     delete fPtYGenerator;
0033 }
0034
0035 void Pythia6PtYDistGun::generateEvent(CLHEP::HepRandomEngine* engine) {
0036   Pythia6Service::InstanceWrapper guard(fPy6Service);  // grab Py6 instance
0037
0038   // now actualy, start cooking up the event gun
0039   //
0040
0041   // 1st, primary vertex
0042   //
0043   HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(0., 0., 0.));
0044
0045   // here re-create fEvt (memory)
0046   //
0047   fEvt = new HepMC::GenEvent();
0048
0049   int ip = 1;
0050
0051   for (size_t i = 0; i < fPartIDs.size(); i++) {
0052     int particleID = fPartIDs[i];  // this is PDG - need to convert to Py6 !!!
0053     int py6PID = HepPID::translatePDTtoPythia(particleID);
0054
0055     int dum = 0;
0056     double pt = 0, y = 0, u = 0, ee = 0, the = 0;
0057
0058     pt = fPtYGenerator->firePt(engine);
0059     y = fPtYGenerator->fireY(engine);
0060     u = exp(y);
0061
0062     double mass = pymass_(py6PID);
0063
0064     // fill p(ip,5) (in PYJETS) with mass value right now,
0065     // because the (hardcoded) mstu(10)=1 will make py1ent
0066     // pick the mass from there
0067     pyjets.p[4][ip - 1] = mass;
0068
0069     ee = 0.5 * std::sqrt(mass * mass + pt * pt) * (u * u + 1) / u;
0070
0071     double pz = std::sqrt(ee * ee - pt * pt - mass * mass);
0072     if (y < 0.)
0073       pz *= -1;
0074
0075     double phi = (fMaxPhi - fMinPhi) * pyr_(&dum) + fMinPhi;
0076
0077     the = atan(pt / pz);
0078     if (pz < 0.)
0079       the += M_PI;
0080
0081     py1ent_(ip, py6PID, ee, the, phi);
0082
0083     /*
0084          double px     = pt*cos(phi) ;
0085          double py     = pt*sin(phi) ;
0086 */
0087     double px = pyjets.p[0][ip - 1];  // pt*cos(phi) ;
0088     double py = pyjets.p[1][ip - 1];  // pt*sin(phi) ;
0089     pz = pyjets.p[2][ip - 1];         // mom*cos(the) ;
0090
0091     HepMC::FourVector p(px, py, pz, ee);
0092     HepMC::GenParticle* Part = new HepMC::GenParticle(p, particleID, 1);
0093     Part->suggest_barcode(ip);
0094     Vtx->add_particle_out(Part);
0095
0096     ip++;
0097   }
0098
0099   fEvt->add_vertex(Vtx);
0100
0101   // run pythia
0102   pyexec_();
0103
0104   return;
0105 }
0106
0107 DEFINE_FWK_MODULE(Pythia6PtYDistGun);