Pythia6PtYDistGun.cc

Go to the documentation of this file.

 
#include <iostream>
 
#include "Pythia6PtYDistGun.h"
#include "GeneratorInterface/Pythia6Interface/interface/PtYDistributor.h"
 
#include "FWCore/Utilities/interface/Exception.h"
 
#include "FWCore/Framework/interface/EDProducer.h"
#include "FWCore/Framework/interface/EventSetup.h"
 
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
 
#include "FWCore/Framework/interface/MakerMacros.h"
 
using namespace edm;
using namespace gen;
 
Pythia6PtYDistGun::Pythia6PtYDistGun(const ParameterSet& pset) : Pythia6ParticleGun(pset), fPtYGenerator(nullptr) {
  ParameterSet pgun_params = pset.getParameter<ParameterSet>("PGunParameters");
 
  fPtYGenerator = new PtYDistributor(pgun_params.getParameter<FileInPath>("kinematicsFile"),
                                     pgun_params.getParameter<double>("MaxPt"),
                                     pgun_params.getParameter<double>("MinPt"),
                                     pgun_params.getParameter<double>("MaxY"),
                                     pgun_params.getParameter<double>("MinY"),
                                     pgun_params.getParameter<int>("PtBinning"),
                                     pgun_params.getParameter<int>("YBinning"));
}
 
Pythia6PtYDistGun::~Pythia6PtYDistGun() {
  if (fPtYGenerator)
    delete fPtYGenerator;
}
 
void Pythia6PtYDistGun::generateEvent(CLHEP::HepRandomEngine* engine) {
  Pythia6Service::InstanceWrapper guard(fPy6Service);  // grab Py6 instance
 
  // now actualy, start cooking up the event gun
  //
 
  // 1st, primary vertex
  //
  HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(0., 0., 0.));
 
  // here re-create fEvt (memory)
  //
  fEvt = new HepMC::GenEvent();
 
  int ip = 1;
 
  for (size_t i = 0; i < fPartIDs.size(); i++) {
    int particleID = fPartIDs[i];  // this is PDG - need to convert to Py6 !!!
    int py6PID = HepPID::translatePDTtoPythia(particleID);
 
    int dum = 0;
    double pt = 0, y = 0, u = 0, ee = 0, the = 0;
 
    pt = fPtYGenerator->firePt(engine);
    y = fPtYGenerator->fireY(engine);
    u = exp(y);
 
    double mass = pymass_(py6PID);
 
    // fill p(ip,5) (in PYJETS) with mass value right now,
    // because the (hardcoded) mstu(10)=1 will make py1ent
    // pick the mass from there
    pyjets.p[4][ip - 1] = mass;
 
    ee = 0.5 * std::sqrt(mass * mass + pt * pt) * (u * u + 1) / u;
 
    double pz = std::sqrt(ee * ee - pt * pt - mass * mass);
    if (y < 0.)
      pz *= -1;
 
    double phi = (fMaxPhi - fMinPhi) * pyr_(&dum) + fMinPhi;
 
    the = atan(pt / pz);
    if (pz < 0.)
      the += M_PI;
 
    py1ent_(ip, py6PID, ee, the, phi);
 
    /*
         double px     = pt*cos(phi) ;
         double py     = pt*sin(phi) ;
*/
    double px = pyjets.p[0][ip - 1];  // pt*cos(phi) ;
    double py = pyjets.p[1][ip - 1];  // pt*sin(phi) ;
    pz = pyjets.p[2][ip - 1];         // mom*cos(the) ;
 
    HepMC::FourVector p(px, py, pz, ee);
    HepMC::GenParticle* Part = new HepMC::GenParticle(p, particleID, 1);
    Part->suggest_barcode(ip);
    Vtx->add_particle_out(Part);
 
    ip++;
  }
 
  fEvt->add_vertex(Vtx);
 
  // run pythia
  pyexec_();
 
  return;
}
 
DEFINE_FWK_MODULE(Pythia6PtYDistGun);