Pythia6Interface/plugins/Pythia6PtGun.cc

0001
0002 #include <iostream>
0003
0004 #include "Pythia6PtGun.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 Pythia6PtGun::Pythia6PtGun(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   fMinPt = pgun_params.getParameter<double>("MinPt");                            // ,  20.);
0023   fMaxPt = pgun_params.getParameter<double>("MaxPt");                            // , 420.);
0024   fAddAntiParticle = pgun_params.getParameter<bool>("AddAntiParticle");          //, false) ;
0025 }
0026
0027 Pythia6PtGun::~Pythia6PtGun() {}
0028
0029 void Pythia6PtGun::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 pt = 0, mom = 0, ee = 0, the = 0, eta = 0;
0049     double mass = pymass_(py6PID);
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
0058     eta = (fMaxEta - fMinEta) * pyr_(&dum) + fMinEta;
0059
0060     the = 2. * atan(exp(-eta));
0061
0062     pt = (fMaxPt - fMinPt) * pyr_(&dum) + fMinPt;
0063
0064     mom = pt / sin(the);
0065     ee = sqrt(mom * mom + mass * mass);
0066
0067     py1ent_(ip, py6PID, ee, the, phi);
0068
0069     double px = pyjets.p[0][ip - 1];  // pt*cos(phi) ;
0070     double py = pyjets.p[1][ip - 1];  // pt*sin(phi) ;
0071     double pz = pyjets.p[2][ip - 1];  // mom*cos(the) ;
0072
0073     HepMC::FourVector p(px, py, pz, ee);
0074     HepMC::GenParticle* Part = new HepMC::GenParticle(p, particleID, 1);
0075     Part->suggest_barcode(ip);
0076     Vtx->add_particle_out(Part);
0077
0078     if (fAddAntiParticle) {
0079       ip = ip + 1;
0080       HepMC::GenParticle* APart = addAntiParticle(ip, particleID, ee, eta, phi);
0081       if (APart)
0082         Vtx->add_particle_out(APart);
0083     }
0084     ip++;
0085   }
0086
0087   fEvt->add_vertex(Vtx);
0088
0089   // run pythia
0090   pyexec_();
0091
0092   return;
0093 }
0094
0095 DEFINE_FWK_MODULE(Pythia6PtGun);