Pythia6Interface/plugins/Pythia6JetGun.cc

0001
0002 #include <iostream>
0003
0004 #include "Pythia6JetGun.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 Pythia6JetGun::Pythia6JetGun(const ParameterSet& pset)
0018     : Pythia6ParticleGun(pset), fMinEta(0.), fMaxEta(0.), fMinE(0.), fMaxE(0.), fMinP(0.), fMaxP(0.) {
0019   ParameterSet pgun_params = pset.getParameter<ParameterSet>("PGunParameters");
0020   fMinEta = pgun_params.getParameter<double>("MinEta");
0021   fMaxEta = pgun_params.getParameter<double>("MaxEta");
0022   fMinE = pgun_params.getParameter<double>("MinE");
0023   fMaxE = pgun_params.getParameter<double>("MaxE");
0024   fMinP = pgun_params.getParameter<double>("MinP");
0025   fMaxP = pgun_params.getParameter<double>("MaxP");
0026 }
0027
0028 Pythia6JetGun::~Pythia6JetGun() {}
0029
0030 void Pythia6JetGun::generateEvent(CLHEP::HepRandomEngine*) {
0031   Pythia6Service::InstanceWrapper guard(fPy6Service);  // grab Py6 instance
0032
0033   // now actualy, start cooking up the event gun
0034   //
0035
0036   // 1st, primary vertex
0037   //
0038   HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(0., 0., 0.));
0039
0040   // here re-create fEvt (memory)
0041   //
0042   fEvt = new HepMC::GenEvent();
0043
0044   int ip = 1;
0045   double totPx = 0.;
0046   double totPy = 0.;
0047   double totPz = 0.;
0048   double totE = 0.;
0049   double totM = 0.;
0050   double phi, eta, the, ee, pp;
0051   int dum = 0;
0052   for (size_t i = 0; i < fPartIDs.size(); i++) {
0053     int particleID = fPartIDs[i];  // this is PDG - need to convert to Py6 !!!
0054     int py6PID = HepPID::translatePDTtoPythia(particleID);
0055
0056     // internal numbers
0057     //
0058     phi = 2. * M_PI * pyr_(&dum);
0059     the = std::acos(-1. + 2. * pyr_(&dum));
0060
0061     // from input
0062     //
0063     ee = (fMaxE - fMinE) * pyr_(&dum) + fMinE;
0064
0065     // fill p(ip,5) (in PYJETS) with mass value right now,
0066     // because the (hardcoded) mstu(10)=1 will make py1ent
0067     // pick the mass from there
0068     double mass = pymass_(py6PID);
0069     pyjets.p[4][ip - 1] = mass;
0070
0071     // add entry to py6
0072     //
0073     py1ent_(ip, py6PID, ee, the, phi);
0074
0075     // values for computing total mass
0076     //
0077     totPx += pyjets.p[0][ip - 1];
0078     totPy += pyjets.p[1][ip - 1];
0079     totPz += pyjets.p[2][ip - 1];
0080     totE += pyjets.p[3][ip - 1];
0081
0082     ip++;
0083
0084   }  // end forming up py6 record of the jet
0085
0086   // compute total mass
0087   //
0088   totM = std::sqrt(totE * totE - (totPx * totPx + totPy * totPy + totPz * totPz));
0089
0090   //now the boost (from input params)
0091   //
0092   pp = (fMaxP - fMinP) * pyr_(&dum) + fMinP;
0093   ee = std::sqrt(totM * totM + pp * pp);
0094
0095   //the boost direction (from input params)
0096   //
0097   phi = (fMaxPhi - fMinPhi) * pyr_(&dum) + fMinPhi;
0098   eta = (fMaxEta - fMinEta) * pyr_(&dum) + fMinEta;
0099   the = 2. * atan(exp(-eta));
0100
0101   double betaX = pp / ee * std::sin(the) * std::cos(phi);
0102   double betaY = pp / ee * std::sin(the) * std::sin(phi);
0103   double betaZ = pp / ee * std::cos(the);
0104
0105   // boost all particles
0106   // the first 2 params (-1) tell to boost all (fisrt-to-last),
0107   // and the next 2 params (0.) tell no rotation
0108   //
0109   int first = -1, last = -1;
0110   double rothe = 0, rophi = 0.;
0111
0112   pyrobo_(first, last, rothe, rophi, betaX, betaY, betaZ);
0113
0114   // event should be formed from boosted record !!!
0115   // that's why additional loop
0116   //
0117   for (int i = 0; i < pyjets.n; i++) {
0118     HepMC::FourVector p(pyjets.p[0][i], pyjets.p[1][i], pyjets.p[2][i], pyjets.p[3][i]);
0119     HepMC::GenParticle* Part = new HepMC::GenParticle(p, HepPID::translatePythiatoPDT(pyjets.k[1][i]), 1);
0120     Part->suggest_barcode(i + 1);
0121     Vtx->add_particle_out(Part);
0122   }
0123   fEvt->add_vertex(Vtx);
0124
0125   // run pythia
0126   pyexec_();
0127
0128   return;
0129 }
0130
0131 DEFINE_FWK_MODULE(Pythia6JetGun);