Pythia6Interface/plugins/Pythia6PartonPtGun.cc

0001
0002 #include <iostream>
0003
0004 #include "Pythia6PartonPtGun.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 Pythia6PartonPtGun::Pythia6PartonPtGun(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   fMinPt = pgun_params.getParameter<double>("MinPt");                            // ,  20.);
0023   fMaxPt = pgun_params.getParameter<double>("MaxPt");                            // , 420.);
0024 }
0025
0026 Pythia6PartonPtGun::~Pythia6PartonPtGun() {}
0027
0028 void Pythia6PartonPtGun::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 pt = 0, mom = 0, 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
0056   eta = (fMaxEta - fMinEta) * pyr_(&dum) + fMinEta;
0057
0058   the = 2. * atan(exp(-eta));
0059
0060   pt = (fMaxPt - fMinPt) * pyr_(&dum) + fMinPt;
0061
0062   mom = pt / sin(the);
0063   ee = sqrt(mom * mom + mass * mass);
0064
0065   py1ent_(ip, py6PID, ee, the, phi);
0066
0067   double px = pyjets.p[0][ip - 1];  // pt*cos(phi) ;
0068   double py = pyjets.p[1][ip - 1];  // pt*sin(phi) ;
0069   double pz = pyjets.p[2][ip - 1];  // mom*cos(the) ;
0070
0071   HepMC::FourVector p(px, py, pz, ee);
0072   HepMC::GenParticle* Part = new HepMC::GenParticle(p, fPartonID, 1);
0073   Part->suggest_barcode(ip);
0074   Vtx->add_particle_out(Part);
0075
0076   // now add anti-quark
0077   ip = ip + 1;
0078   HepMC::GenParticle* APart = addAntiParticle(ip, fPartonID, ee, eta, phi);
0079   if (APart) {
0080     Vtx->add_particle_out(APart);
0081   } else {
0082     // otherwise it should throw !
0083   }
0084
0085   // this should probably be configurable...
0086   //
0087   double qmax = 2. * ee;
0088
0089   joinPartons(qmax);
0090
0091   fEvt->add_vertex(Vtx);
0092
0093   // run pythia
0094   pyexec_();
0095
0096   return;
0097 }
0098
0099 DEFINE_FWK_MODULE(Pythia6PartonPtGun);