CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
Pythia6JetGun.cc
Go to the documentation of this file.
1 
2 #include <iostream>
3 
4 #include "Pythia6JetGun.h"
5 
7 
9 
11 
13 
14 using namespace edm;
15 using namespace gen;
16 
17 Pythia6JetGun::Pythia6JetGun(const ParameterSet& pset)
18  : Pythia6ParticleGun(pset), fMinEta(0.), fMaxEta(0.), fMinE(0.), fMaxE(0.), fMinP(0.), fMaxP(0.) {
19  ParameterSet pgun_params = pset.getParameter<ParameterSet>("PGunParameters");
20  fMinEta = pgun_params.getParameter<double>("MinEta");
21  fMaxEta = pgun_params.getParameter<double>("MaxEta");
22  fMinE = pgun_params.getParameter<double>("MinE");
23  fMaxE = pgun_params.getParameter<double>("MaxE");
24  fMinP = pgun_params.getParameter<double>("MinP");
25  fMaxP = pgun_params.getParameter<double>("MaxP");
26 }
27 
29 
30 void Pythia6JetGun::generateEvent(CLHEP::HepRandomEngine*) {
31  Pythia6Service::InstanceWrapper guard(fPy6Service); // grab Py6 instance
32 
33  // now actualy, start cooking up the event gun
34  //
35 
36  // 1st, primary vertex
37  //
38  HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(0., 0., 0.));
39 
40  // here re-create fEvt (memory)
41  //
42  fEvt = new HepMC::GenEvent();
43 
44  int ip = 1;
45  double totPx = 0.;
46  double totPy = 0.;
47  double totPz = 0.;
48  double totE = 0.;
49  double totM = 0.;
50  double phi, eta, the, ee, pp;
51  int dum = 0;
52  for (size_t i = 0; i < fPartIDs.size(); i++) {
53  int particleID = fPartIDs[i]; // this is PDG - need to convert to Py6 !!!
54  int py6PID = HepPID::translatePDTtoPythia(particleID);
55 
56  // internal numbers
57  //
58  phi = 2. * M_PI * pyr_(&dum);
59  the = std::acos(-1. + 2. * pyr_(&dum));
60 
61  // from input
62  //
63  ee = (fMaxE - fMinE) * pyr_(&dum) + fMinE;
64 
65  // fill p(ip,5) (in PYJETS) with mass value right now,
66  // because the (hardcoded) mstu(10)=1 will make py1ent
67  // pick the mass from there
68  double mass = pymass_(py6PID);
69  pyjets.p[4][ip - 1] = mass;
70 
71  // add entry to py6
72  //
73  py1ent_(ip, py6PID, ee, the, phi);
74 
75  // values for computing total mass
76  //
77  totPx += pyjets.p[0][ip - 1];
78  totPy += pyjets.p[1][ip - 1];
79  totPz += pyjets.p[2][ip - 1];
80  totE += pyjets.p[3][ip - 1];
81 
82  ip++;
83 
84  } // end forming up py6 record of the jet
85 
86  // compute total mass
87  //
88  totM = std::sqrt(totE * totE - (totPx * totPx + totPy * totPy + totPz * totPz));
89 
90  //now the boost (from input params)
91  //
92  pp = (fMaxP - fMinP) * pyr_(&dum) + fMinP;
93  ee = std::sqrt(totM * totM + pp * pp);
94 
95  //the boost direction (from input params)
96  //
97  phi = (fMaxPhi - fMinPhi) * pyr_(&dum) + fMinPhi;
98  eta = (fMaxEta - fMinEta) * pyr_(&dum) + fMinEta;
99  the = 2. * atan(exp(-eta));
100 
101  double betaX = pp / ee * std::sin(the) * std::cos(phi);
102  double betaY = pp / ee * std::sin(the) * std::sin(phi);
103  double betaZ = pp / ee * std::cos(the);
104 
105  // boost all particles
106  // the first 2 params (-1) tell to boost all (fisrt-to-last),
107  // and the next 2 params (0.) tell no rotation
108  //
109  int first = -1, last = -1;
110  double rothe = 0, rophi = 0.;
111 
112  pyrobo_(first, last, rothe, rophi, betaX, betaY, betaZ);
113 
114  // event should be formed from boosted record !!!
115  // that's why additional loop
116  //
117  for (int i = 0; i < pyjets.n; i++) {
118  HepMC::FourVector p(pyjets.p[0][i], pyjets.p[1][i], pyjets.p[2][i], pyjets.p[3][i]);
119  HepMC::GenParticle* Part = new HepMC::GenParticle(p, HepPID::translatePythiatoPDT(pyjets.k[1][i]), 1);
120  Part->suggest_barcode(i + 1);
121  Vtx->add_particle_out(Part);
122  }
123  fEvt->add_vertex(Vtx);
124 
125  // run pythia
126  pyexec_();
127 
128  return;
129 }
130 
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
double pyr_(int *idummy)
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
HepMC::GenEvent * fEvt
Definition: Pythia6Gun.h:68
void generateEvent(CLHEP::HepRandomEngine *) override
double p[5][pyjets_maxn]
double fMinPhi
Definition: Pythia6Gun.h:63
T sqrt(T t)
Definition: SSEVec.h:19
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
double pymass_(int &)
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
std::vector< int > fPartIDs
#define M_PI
double fMaxPhi
Definition: Pythia6Gun.h:64
void py1ent_(int &ip, int &kf, double &pe, double &the, double &phi)
HLT enums.
Pythia6Service * fPy6Service
Definition: Pythia6Gun.h:58
~Pythia6JetGun() override
void pyrobo_(int &, int &, double &, double &, double &, double &, double &)
void pyexec_()