Pythia6Interface/plugins/Pythia6EGun.cc

0001
0002 #include <iostream>
0003
0004 #include "Pythia6EGun.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 Pythia6EGun::Pythia6EGun(const ParameterSet& pset) : Pythia6ParticleGun(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");                              // ,  0.);
0023   fMaxE = pgun_params.getParameter<double>("MaxE");                              // ,  0.);
0024   fAddAntiParticle = pgun_params.getParameter<bool>("AddAntiParticle");          //, false) ;
0025 }
0026
0027 Pythia6EGun::~Pythia6EGun() {}
0028
0029 void Pythia6EGun::generateEvent(CLHEP::HepRandomEngine*) {
0030   Pythia6Service::InstanceWrapper guard(fPy6Service);  // grab Py6 instance
0031
0032   // now actualy, start cooking up the event gun
0033   //
0034
0035   // 1st, primary vertex
0036   //
0037   HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(0., 0., 0.));
0038
0039   // here re-create fEvt (memory)
0040   //
0041   fEvt = new HepMC::GenEvent();
0042
0043   int ip = 1;
0044   for (size_t i = 0; i < fPartIDs.size(); i++) {
0045     int particleID = fPartIDs[i];  // this is PDG - need to convert to Py6 !!!
0046     int py6PID = HepPID::translatePDTtoPythia(particleID);
0047     int dum = 0;
0048     double ee = 0, the = 0, eta = 0;
0049     double mass = pymass_(particleID);
0050
0051     // fill p(ip,5) (in PYJETS) with mass value right now,
0052     // because the (hardcoded) mstu(10)=1 will make py1ent
0053     // pick the mass from there
0054     pyjets.p[4][ip - 1] = mass;
0055
0056     double phi = (fMaxPhi - fMinPhi) * pyr_(&dum) + fMinPhi;
0057     ee = (fMaxE - fMinE) * pyr_(&dum) + fMinE;
0058     eta = (fMaxEta - fMinEta) * pyr_(&dum) + fMinEta;
0059     the = 2. * atan(exp(-eta));
0060
0061     py1ent_(ip, py6PID, ee, the, phi);
0062
0063     double px = pyjets.p[0][ip - 1];  // pt*cos(phi) ;
0064     double py = pyjets.p[1][ip - 1];  // pt*sin(phi) ;
0065     double pz = pyjets.p[2][ip - 1];  // mom*cos(the) ;
0066
0067     HepMC::FourVector p(px, py, pz, ee);
0068     HepMC::GenParticle* Part = new HepMC::GenParticle(p, particleID, 1);
0069     Part->suggest_barcode(ip);
0070     Vtx->add_particle_out(Part);
0071
0072     if (fAddAntiParticle) {
0073       ip = ip + 1;
0074       HepMC::GenParticle* APart = addAntiParticle(ip, particleID, ee, eta, phi);
0075       if (APart)
0076         Vtx->add_particle_out(APart);
0077     }
0078     ip++;
0079   }
0080
0081   fEvt->add_vertex(Vtx);
0082
0083   // run pythia
0084   pyexec_();
0085
0086   return;
0087 }
0088
0089 DEFINE_FWK_MODULE(Pythia6EGun);