CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
gentop.cc
Go to the documentation of this file.
1 //
2 //
3 // File: src/gentop.cc
4 // Purpose: Toy ttbar event generator for testing.
5 // Created: Jul, 2000, sss.
6 //
7 // CMSSW File : src/gentop.cc
8 // Original Author : Scott Stuart Snyder <snyder@bnl.gov> for D0
9 // Imported to CMSSW by Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>
10 //
11 
37 #include "CLHEP/Random/RandFlat.h"
38 #include "CLHEP/Random/RandExponential.h"
39 #include "CLHEP/Random/RandBreitWigner.h"
40 #include "CLHEP/Random/RandGauss.h"
41 #include "CLHEP/Units/PhysicalConstants.h"
45 #include <cmath>
46 #include <ostream>
47 
48 using std::ostream;
49 
50 // Event number counter.
54 namespace {
55  int next_evnum = 0;
56 }
57 
58 namespace hitfit {
59 
60  //**************************************************************************
61  // Argument handling.
62  //
63 
65  //
66  // Purpose: Constructor.
67  //
68  // Inputs:
69  // defs - The Defaults instance from which to initialize.
70  //
71  : _t_pt_mean(defs.get_float("t_pt_mean")),
72  _mt(defs.get_float("mt")),
73  _sigma_mt(defs.get_float("sigma_mt")),
74  _mh(defs.get_float("mh")),
75  _sigma_mh(defs.get_float("sigma_mh")),
76  _recoil_pt_mean(defs.get_float("recoil_pt_mean")),
77  _boost_sigma(defs.get_float("boost_sigma")),
78  _m_boost(defs.get_float("m_boost")),
79  _mb(defs.get_float("mb")),
80  _sigma_mb(defs.get_float("sigma_mb")),
81  _mw(defs.get_float("mw")),
82  _sigma_mw(defs.get_float("sigma_mw")),
83  _svx_tageff(defs.get_float("svx_tageff")),
84  _smear(defs.get_bool("smear")),
85  _smear_dir(defs.get_bool("smear_dir")),
86  _muon(defs.get_bool("muon")),
87  _ele_res_str(defs.get_string("ele_res_str")),
88  _muo_res_str(defs.get_string("muo_res_str")),
89  _jet_res_str(defs.get_string("jet_res_str")),
90  _kt_res_str(defs.get_string("kt_res_str")) {}
91 
93  //
94  // Purpose: Return the t_pt_mean parameter.
95  // See the header for documentation.
96  //
97  {
98  return _t_pt_mean;
99  }
100 
102  //
103  // Purpose: Return the mt parameter.
104  // See the header for documentation.
105  //
106  {
107  return _mt;
108  }
109 
111  //
112  // Purpose: Return the sigma_mt parameter.
113  // See the header for documentation.
114  //
115  {
116  return _sigma_mt;
117  }
118 
120  //
121  // Purpose: Return the svx_tageff parameter.
122  // See the header for documentation.
123  //
124  {
125  return _svx_tageff;
126  }
127 
129  //
130  // Purpose: Return the mh parameter.
131  // See the header for documentation.
132  //
133  {
134  return _mh;
135  }
136 
138  //
139  // Purpose: Return the sigma_mh parameter.
140  // See the header for documentation.
141  //
142  {
143  return _sigma_mh;
144  }
145 
147  //
148  // Purpose: Return the smear parameter.
149  // See the header for documentation.
150  //
151  {
152  return _smear;
153  }
154 
156  //
157  // Purpose: Return the smear_dir parameter.
158  // See the header for documentation.
159  //
160  {
161  return _smear_dir;
162  }
163 
165  //
166  // Purpose: Return the muon parameter.
167  // See the header for documentation.
168  //
169  {
170  return _muon;
171  }
172 
174  //
175  // Purpose: Return the recoil_pt_mean parameter.
176  // See the header for documentation.
177  //
178  {
179  return _recoil_pt_mean;
180  }
181 
183  //
184  // Purpose: Return the boost_sigma parameter.
185  // See the header for documentation.
186  //
187  {
188  return _boost_sigma;
189  }
190 
192  //
193  // Purpose: Return the m_boost parameter.
194  // See the header for documentation.
195  //
196  {
197  return _m_boost;
198  }
199 
201  //
202  // Purpose: Return the mb parameter.
203  // See the header for documentation.
204  //
205  {
206  return _mb;
207  }
208 
210  //
211  // Purpose: Return the sigma_mb parameter.
212  // See the header for documentation.
213  //
214  {
215  return _sigma_mb;
216  }
217 
219  //
220  // Purpose: Return the mw parameter.
221  // See the header for documentation.
222  //
223  {
224  return _mw;
225  }
226 
228  //
229  // Purpose: Return the sigma_mw parameter.
230  // See the header for documentation.
231  //
232  {
233  return _sigma_mw;
234  }
235 
237  //
238  // Purpose: Return the ele_res_str parameter.
239  // See the header for documentation.
240  //
241  {
242  return _ele_res_str;
243  }
244 
246  //
247  // Purpose: Return the muo_res_str parameter.
248  // See the header for documentation.
249  //
250  {
251  return _muo_res_str;
252  }
253 
255  //
256  // Purpose: Return the jet_res_str parameter.
257  // See the header for documentation.
258  //
259  {
260  return _jet_res_str;
261  }
262 
264  //
265  // Purpose: Return the kt_res_str parameter.
266  // See the header for documentation.
267  //
268  {
269  return _kt_res_str;
270  }
271 
272  //**************************************************************************
273  // Internal helper functions.
274  //
275 
276  namespace {
277 
284  Threevec rand_spher(CLHEP::HepRandomEngine& engine)
285  //
286  // Purpose: Return a unit vector with (theta, phi) uniformly distributed
287  // over a sphere.
288  //
289  // Inputs:
290  // engine - The underlying RNG.
291  //
292  // Returns:
293  // The generated vector.
294  //
295  {
296  CLHEP::RandFlat r(engine);
297 
298  Threevec v;
299 
300  double U = r.fire(0.0, 1.0);
301  double V = r.fire(0.0, 1.0);
302 
303  double theta = 2.0 * CLHEP::pi * U;
304  double phi = acos(2 * V - 1.0);
305 
306  double x = sin(theta) * cos(phi);
307  double y = sin(theta) * sin(phi);
308  double z = cos(theta);
309 
310  v = Threevec(x, y, z);
311 
312  return v.unit();
313  }
314 
329  Fourvec make_massive(const Threevec& p, double m_true, double sigma, CLHEP::HepRandomEngine& engine)
330  //
331  // Purpose: Given a direction, mass, and width, choose a mass from a
332  // Breit-Wigner and return a 4-vector with the chosen mass
333  // and the specified direction.
334  //
335  // Inputs:
336  // p - The direction.
337  // m_true - The mean for the Breit-Wigner.
338  // sigma - The width for the Breit-Wigner.
339  // engine - The underlying RNG.
340  //
341  // Returns:
342  // The generated 4-vector.
343  //
344  {
345  CLHEP::RandBreitWigner rbw(engine);
346  double m = rbw.fire(m_true, sigma);
347  return Fourvec(p, sqrt(m * m + p.mag2()));
348  }
349 
367  void decay(const Fourvec& v, double m1, double m2, CLHEP::HepRandomEngine& engine, Fourvec& vout1, Fourvec& vout2)
368  //
369  // Purpose: v decays into two particles w/masses m1, m2.
370  //
371  // Inputs:
372  // v - The incoming 4-vector.
373  // m1 - Mass of the first decay product.
374  // m2 - Mass of the second decay product.
375  // engine - The underlying RNG.
376  //
377  // Outputs:
378  // vout1 - Outgoing 4-vector of the first decay product.
379  // vout2 - Outgoing 4-vector of the second decay product.
380  //
381  {
382  // Construct a decay in the incoming particle's rest frame,
383  // uniformly distributed in direction.
384  Threevec p = rand_spher(engine);
385  double m0 = v.m();
386 
387  if (m1 + m2 > m0) {
388  // What ya gonna do?
389  double f = m0 / (m1 + m2);
390  m1 *= f;
391  m2 *= f;
392  }
393 
394  double m0_2 = m0 * m0;
395  double m1_2 = m1 * m1;
396  double m2_2 = m2 * m2;
397 
398  // Calculate the 3-momentum of each particle in the decay frame.
399  p *= 0.5 / m0 *
400  sqrt(m0_2 * m0_2 + m1_2 * m1_2 + m2_2 * m2_2 - 2 * m0_2 * m1_2 - 2 * m0_2 * m2_2 - 2 * m1_2 * m2_2);
401  double p2 = p.mag2();
402 
403  vout1 = Fourvec(p, sqrt(p2 + m1_2));
404  vout2 = Fourvec(-p, sqrt(p2 + m2_2));
405 
406  // Boost out of the rest frame.
407  vout1.boost(v.boostVector());
408  vout2.boost(v.boostVector());
409  }
410 
420  Threevec rand_pt(double pt_mean, CLHEP::HepRandomEngine& engine)
421  //
422  // Purpose: Generate a vector in a (uniformly) random direction,
423  // with pt chosen from an exponential distribution with mean pt_mean.
424  //
425  // Inputs:
426  // pt_mean - The mean of the distribution.
427  // engine - The underlying RNG.
428  //
429  // Returns:
430  // The generated vector.
431  //
432  {
433  CLHEP::RandExponential rexp(engine);
434 
435  // A random direction.
436  Threevec p = rand_spher(engine);
437 
438  // Scale by random pt.
439  p *= (rexp.fire(pt_mean) / p.perp());
440 
441  return p;
442  }
443 
451  Fourvec rand_boost(const Gentop_Args& args, CLHEP::HepRandomEngine& engine)
452  //
453  // Purpose: Generate a random boost for the event.
454  //
455  // Inputs:
456  // args - The parameter settings.
457  // engine - The underlying RNG.
458  //
459  // Returns:
460  // The generated boost.
461  //
462  {
463  CLHEP::RandExponential rexp(engine);
464  CLHEP::RandFlat rflat(engine);
465  CLHEP::RandGauss rgauss(engine);
466 
467  // Boost in pt and z.
468  Threevec p(1, 0, 0);
469  p.rotateZ(rflat.fire(0, 2 * M_PI));
470  p *= rexp.fire(args.recoil_pt_mean());
471  p.setZ(rgauss.fire(0, args.boost_sigma()));
472  return Fourvec(p, sqrt(p.mag2() + args.m_boost() * args.m_boost()));
473  }
474 
484  void tagsim(const Gentop_Args& args, Lepjets_Event& ev, CLHEP::HepRandomEngine& engine)
485  //
486  // Purpose: Simulate SVX tagging.
487  //
488  // Inputs:
489  // args - The parameter settings.
490  // ev - The event to tag.
491  // engine - The underlying RNG.
492  //
493  // Outputs:
494  // ev - The event with tags filled in.
495  //
496  {
497  CLHEP::RandFlat rflat(engine);
498  for (std::vector<Lepjets_Event_Jet>::size_type i = 0; i < ev.njets(); i++) {
499  int typ = ev.jet(i).type();
500  if (typ == hadb_label || typ == lepb_label || typ == higgs_label) {
501  if (rflat.fire() < args.svx_tageff())
502  ev.jet(i).svx_tag() = true;
503  }
504  }
505  }
506 
507  } // unnamed namespace
508 
509  //**************************************************************************
510  // External interface.
511  //
512 
513  Lepjets_Event gentop(const Gentop_Args& args, CLHEP::HepRandomEngine& engine)
514  //
515  // Purpose: Generate a ttbar -> ljets event.
516  //
517  // Inputs:
518  // args - The parameter settings.
519  // engine - The underlying RNG.
520  //
521  // Returns:
522  // The generated event.
523  //
524  {
525  CLHEP::RandBreitWigner rbw(engine);
526  CLHEP::RandGauss rgauss(engine);
527 
528  // Get the t decay momentum in the ttbar rest frame.
529  Threevec p = rand_pt(args.t_pt_mean(), engine);
530 
531  // Make the t/tbar vectors.
532  Fourvec lept = make_massive(p, args.mt(), args.sigma_mt(), engine);
533  Fourvec hadt = make_massive(-p, args.mt(), args.sigma_mt(), engine);
534 
535  // Boost the rest frame.
536  Fourvec boost = rand_boost(args, engine);
537  lept.boost(boost.boostVector());
538  hadt.boost(boost.boostVector());
539 
540  // Decay t -> b W, leptonic side.
541  Fourvec lepb, lepw;
542  double mlb = rgauss.fire(args.mb(), args.sigma_mb());
543  double mlw = rbw.fire(args.mw(), args.sigma_mw());
544  decay(lept, mlb, mlw, engine, lepb, lepw);
545 
546  // Decay t -> b W, hadronic side.
547  Fourvec hadb, hadw;
548  double mhb = rgauss.fire(args.mb(), args.sigma_mb());
549  double mhw = rbw.fire(args.mw(), args.sigma_mw());
550  decay(hadt, mhb, mhw, engine, hadb, hadw);
551 
552  // Decay W -> l nu.
553  Fourvec lep, nu;
554  decay(lepw, 0, 0, engine, lep, nu);
555 
556  // Decay W -> qqbar.
557  Fourvec q1, q2;
558  decay(hadw, 0, 0, engine, q1, q2);
559 
560  // Fill in the event.
561  Lepjets_Event ev(0, ++next_evnum);
562  Vector_Resolution lep_res(args.muon() ? args.muo_res_str() : args.ele_res_str());
563  Vector_Resolution jet_res(args.jet_res_str());
564  Resolution kt_res = (args.kt_res_str());
565 
566  ev.add_lep(Lepjets_Event_Lep(lep, args.muon() ? muon_label : electron_label, lep_res));
567 
568  ev.add_jet(Lepjets_Event_Jet(lepb, lepb_label, jet_res));
569  ev.add_jet(Lepjets_Event_Jet(hadb, hadb_label, jet_res));
570  ev.add_jet(Lepjets_Event_Jet(q1, hadw1_label, jet_res));
571  ev.add_jet(Lepjets_Event_Jet(q2, hadw2_label, jet_res));
572 
573  ev.met() = nu;
574  ev.kt_res() = kt_res;
575 
576  // Simulate SVX tagging.
577  tagsim(args, ev, engine);
578 
579  // Smear the event, if requested.
580  if (args.smear())
581  ev.smear(engine, args.smear_dir());
582 
583  // Done!
584  return ev;
585  }
586 
587  Lepjets_Event gentth(const Gentop_Args& args, CLHEP::HepRandomEngine& engine)
588  //
589  // Purpose: Generate a ttH -> ljets event.
590  //
591  // Inputs:
592  // args - The parameter settings.
593  // engine - The underlying RNG.
594  //
595  // Returns:
596  // The generated event.
597  //
598  {
599  CLHEP::RandBreitWigner rbw(engine);
600  CLHEP::RandGauss rgauss(engine);
601 
602  // Generate three-vectors for two tops.
603  Threevec p_t1 = rand_pt(args.t_pt_mean(), engine);
604  Threevec p_t2 = rand_pt(args.t_pt_mean(), engine);
605 
606  // Conserve momentum.
607  Threevec p_h = -(p_t1 + p_t2);
608 
609  // Construct the 4-vectors.
610  Fourvec lept = make_massive(p_t1, args.mt(), args.sigma_mt(), engine);
611  Fourvec hadt = make_massive(p_t2, args.mt(), args.sigma_mt(), engine);
612  Fourvec higgs = make_massive(p_h, args.mh(), args.sigma_mh(), engine);
613 
614  // Boost the rest frame.
615  Fourvec boost = rand_boost(args, engine);
616  lept.boost(boost.boostVector());
617  hadt.boost(boost.boostVector());
618  higgs.boost(boost.boostVector());
619 
620  // Decay t -> b W, leptonic side.
621  Fourvec lepb, lepw;
622  decay(lept, rgauss.fire(args.mb(), args.sigma_mb()), rbw.fire(args.mw(), args.sigma_mw()), engine, lepb, lepw);
623 
624  // Decay t -> b W, hadronic side.
625  Fourvec hadb, hadw;
626  decay(hadt, rgauss.fire(args.mb(), args.sigma_mb()), rbw.fire(args.mw(), args.sigma_mw()), engine, hadb, hadw);
627 
628  // Decay W -> l nu.
629  Fourvec lep, nu;
630  decay(lepw, 0, 0, engine, lep, nu);
631 
632  // Decay W -> qqbar.
633  Fourvec q1, q2;
634  decay(hadw, 0, 0, engine, q1, q2);
635 
636  // Decay H -> bbbar.
637  Fourvec hb1, hb2;
638  decay(higgs, rgauss.fire(args.mb(), args.sigma_mb()), rgauss.fire(args.mb(), args.sigma_mb()), engine, hb1, hb2);
639 
640  // Fill in the event.
641  Lepjets_Event ev(0, ++next_evnum);
642  Vector_Resolution lep_res(args.muon() ? args.muo_res_str() : args.ele_res_str());
643  Vector_Resolution jet_res(args.jet_res_str());
644  Resolution kt_res = (args.kt_res_str());
645 
646  ev.add_lep(Lepjets_Event_Lep(lep, args.muon() ? muon_label : electron_label, lep_res));
647 
648  ev.add_jet(Lepjets_Event_Jet(lepb, lepb_label, jet_res));
649  ev.add_jet(Lepjets_Event_Jet(hadb, hadb_label, jet_res));
650  ev.add_jet(Lepjets_Event_Jet(q1, hadw1_label, jet_res));
651  ev.add_jet(Lepjets_Event_Jet(q2, hadw2_label, jet_res));
652  ev.add_jet(Lepjets_Event_Jet(hb1, higgs_label, jet_res));
653  ev.add_jet(Lepjets_Event_Jet(hb2, higgs_label, jet_res));
654 
655  ev.met() = nu;
656  ev.kt_res() = kt_res;
657 
658  // Simulate SVX tagging.
659  tagsim(args, ev, engine);
660 
661  // Smear the event, if requested.
662  if (args.smear())
663  ev.smear(engine, args.smear_dir());
664 
665  // Done!
666  return ev;
667  }
668 
669 } // namespace hitfit
std::string _jet_res_str
Definition: gentop.h:337
Define an abstract interface for getting parameter settings.
double _sigma_mt
Definition: gentop.h:252
Definition: CLHEP.h:16
Represent a lepton in an instance of Lepjets_Event class. This class hold the following information: ...
Define three-vector and four-vector classes for the HitFit package, and supply a few additional opera...
Calculate and represent resolution for a physical quantity.
Definition: Resolution.h:98
double mw() const
Return the value of mw parameter.
Definition: gentop.cc:218
Gentop_Args(const Defaults &defs)
Constructor, initialize an instance of Gentop_Args from an instance of Defaults object.
Definition: gentop.cc:64
double _t_pt_mean
Definition: gentop.h:242
double boost_sigma() const
Return the value of boost_sigma parameter.
Definition: gentop.cc:182
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Lepjets_Event gentth(const Gentop_Args &args, CLHEP::HepRandomEngine &engine)
Generate a event.
Definition: gentop.cc:587
double _sigma_mh
Definition: gentop.h:262
double m_boost() const
Return the value of m_boost parameter.
Definition: gentop.cc:191
bool muon() const
Return the value of muon parameter.
Definition: gentop.cc:164
std::string ele_res_str() const
Return the value of ele_res_str parameter.
Definition: gentop.cc:236
double _svx_tageff
Definition: gentop.h:304
uint16_t size_type
std::string _kt_res_str
Definition: gentop.h:343
const Double_t pi
double _sigma_mb
Definition: gentop.h:288
double _recoil_pt_mean
Definition: gentop.h:268
uint32_t T const *__restrict__ uint32_t const *__restrict__ int32_t int Histo::index_type cudaStream_t V
Represent a simple event consisting of lepton(s) and jet(s).
double sigma_mt() const
Return the value of sigma_mt parameter.
Definition: gentop.cc:110
double mt() const
Return the value of mt parameter.
Definition: gentop.cc:101
T sqrt(T t)
Definition: SSEVec.h:19
Represent a simple event consisting of lepton(s) and jet(s). An instance of this class holds a list o...
Definition: Lepjets_Event.h:62
Lepjets_Event gentop(const Gentop_Args &args, CLHEP::HepRandomEngine &engine)
Generate a event.
Definition: gentop.cc:513
std::string _ele_res_str
Definition: gentop.h:327
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
CLHEP::HepLorentzVector Fourvec
Typedef for a HepLorentzVector.
Definition: fourvec.h:55
A class to represent a jet in an instance of Lepjets_Event class. The class is derived from the Lepje...
double f[11][100]
bool smear_dir() const
Return the value of smear_dir parameter.
Definition: gentop.cc:155
double _sigma_mw
Definition: gentop.h:298
Hold on to parameters for the toy event generator.
Definition: gentop.h:69
#define M_PI
bool smear() const
Return the value of smear parameter.
Definition: gentop.cc:146
CLHEP::Hep3Vector Threevec
Typedef for a Hep3Vector.
Definition: fourvec.h:60
std::string _muo_res_str
Definition: gentop.h:332
std::string jet_res_str() const
Return the value of jet_res_str parameter.
Definition: gentop.cc:254
double recoil_pt_mean() const
Return the value of recoil_pt_mean parameter.
Definition: gentop.cc:173
double t_pt_mean() const
Return the value of t_pt_mean parameter.
Definition: gentop.cc:92
double sigma_mw() const
Return the value of sigma_mw parameter.
Definition: gentop.cc:227
Define an interface for getting parameter settings.
Definition: Defaults.h:57
double _boost_sigma
Definition: gentop.h:273
std::string kt_res_str() const
Return the value of kt_res_str parameter.
Definition: gentop.cc:263
float x
double sigma_mb() const
Return the value of sigma_mb parameter.
Definition: gentop.cc:209
double sigma_mh() const
Return the value of sigma_mh parameter.
Definition: gentop.cc:137
Calculate and represent resolution for a vector of , pseudorapidity , and azimuthal angle ...
double mb() const
Return the value of mb parameter.
Definition: gentop.cc:200
double svx_tageff() const
Return the value of svx_tageff parameter.
Definition: gentop.cc:119
double mh() const
Return the value of mh parameter.
Definition: gentop.cc:128
Geom::Theta< T > theta() const
A toy event generator for events.
std::string muo_res_str() const
Return the value of muon_res_str parameter.
Definition: gentop.cc:245