Pythia6Interface/plugins/Pythia6PartonEGun.cc

0001
0002 #include <iostream>
0003
0004 #include "Pythia6PartonEGun.h"
0005
0006 #include "FWCore/Utilities/interface/Exception.h"
0007
0008 #include "FWCore/Framework/interface/EventSetup.h"
0009
0010 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0011
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013
0014 using namespace edm;
0015 using namespace gen;
0016
0017 Pythia6PartonEGun::Pythia6PartonEGun(const ParameterSet& pset) : Pythia6PartonGun(pset) {
0018   // ParameterSet defpset ;
0019   ParameterSet pgun_params = pset.getParameter<ParameterSet>("PGunParameters");  //, defpset ) ;
0020   fMinEta = pgun_params.getParameter<double>("MinEta");                          // ,-2.2);
0021   fMaxEta = pgun_params.getParameter<double>("MaxEta");                          // , 2.2);
0022   fMinE = pgun_params.getParameter<double>("MinE");                              // ,  20.);
0023   fMaxE = pgun_params.getParameter<double>("MaxE");                              // , 420.);
0024 }
0025
0026 Pythia6PartonEGun::~Pythia6PartonEGun() {}
0027
0028 void Pythia6PartonEGun::generateEvent(CLHEP::HepRandomEngine*) {
0029   Pythia6Service::InstanceWrapper guard(fPy6Service);  // grab Py6 instance
0030
0031   // now actualy, start cooking up the event gun
0032   //
0033
0034   // 1st, primary vertex
0035   //
0036   HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(0., 0., 0.));
0037
0038   // here re-create fEvt (memory)
0039   //
0040   fEvt = new HepMC::GenEvent();
0041
0042   int ip = 1;
0043
0044   int py6PID = HepPID::translatePDTtoPythia(fPartonID);
0045   int dum = 0;
0046   double ee = 0, the = 0, eta = 0;
0047   double mass = pymass_(py6PID);
0048
0049   // fill p(ip,5) (in PYJETS) with mass value right now,
0050   // because the (hardcoded) mstu(10)=1 will make py1ent
0051   // pick the mass from there
0052   pyjets.p[4][ip - 1] = mass;
0053
0054   double phi = (fMaxPhi - fMinPhi) * pyr_(&dum) + fMinPhi;
0055   ee = (fMaxE - fMinE) * pyr_(&dum) + fMinE;
0056   eta = (fMaxEta - fMinEta) * pyr_(&dum) + fMinEta;
0057   the = 2. * atan(exp(-eta));
0058
0059   py1ent_(ip, py6PID, ee, the, phi);
0060
0061   double px = pyjets.p[0][ip - 1];  // pt*cos(phi) ;
0062   double py = pyjets.p[1][ip - 1];  // pt*sin(phi) ;
0063   double pz = pyjets.p[2][ip - 1];  // mom*cos(the) ;
0064
0065   HepMC::FourVector p(px, py, pz, ee);
0066   HepMC::GenParticle* Part = new HepMC::GenParticle(p, fPartonID, 1);
0067   Part->suggest_barcode(ip);
0068   Vtx->add_particle_out(Part);
0069
0070   // now add anti-quark
0071   ip = ip + 1;
0072   HepMC::GenParticle* APart = addAntiParticle(ip, fPartonID, ee, eta, phi);
0073   if (APart) {
0074     Vtx->add_particle_out(APart);
0075   } else {
0076     // otherwise it should throw !
0077   }
0078
0079   // this should probably be configurable...
0080   //
0081   double qmax = 2. * ee;
0082
0083   joinPartons(qmax);
0084
0085   fEvt->add_vertex(Vtx);
0086
0087   // run pythia
0088   pyexec_();
0089
0090   return;
0091 }
0092
0093 DEFINE_FWK_MODULE(Pythia6PartonEGun);