TopHitFit/src/gentop.cc

0001 //
0002 //
0003 // File: src/gentop.cc
0004 // Purpose: Toy ttbar event generator for testing.
0005 // Created: Jul, 2000, sss.
0006 //
0007 // CMSSW File      : src/gentop.cc
0008 // Original Author : Scott Stuart Snyder <snyder@bnl.gov> for D0
0009 // Imported to CMSSW by Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>
0010 //
0011
0012 /**
0013     @file gentop.cc
0014
0015     @brief A toy event generator for
0016     \f$t\bar{t} \to \ell + 4~\mathrm{jets}\f$ events. This file also contains
0017     helper function for generation of random toy events.
0018     See the documentation for the header file gentop.h for details.
0019
0020     @author Scott Stuart Snyder <snyder@bnl.gov>
0021
0022     @par Creation date:
0023     Jul 2000.
0024
0025     @par Modification History:
0026     Apr 2009: Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>:
0027     Imported to CMSSW.<br>
0028     Nov 2009: Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>:
0029     Added doxygen tags for automatic generation of documentation.
0030
0031     @par Terms of Usage:
0032     With consent for the original author (Scott Snyder).
0033
0034  */
0035
0036 #include "TopQuarkAnalysis/TopHitFit/interface/gentop.h"
0037 #include "CLHEP/Random/RandFlat.h"
0038 #include "CLHEP/Random/RandExponential.h"
0039 #include "CLHEP/Random/RandBreitWigner.h"
0040 #include "CLHEP/Random/RandGauss.h"
0041 #include "CLHEP/Units/PhysicalConstants.h"
0042 #include "TopQuarkAnalysis/TopHitFit/interface/fourvec.h"
0043 #include "TopQuarkAnalysis/TopHitFit/interface/Lepjets_Event.h"
0044 #include "TopQuarkAnalysis/TopHitFit/interface/Defaults.h"
0045 #include <cmath>
0046 #include <ostream>
0047
0048 using std::ostream;
0049
0050 // Event number counter.
0051 /**
0052     @brief Event number counter.
0053  */
0054 namespace {
0055   int next_evnum = 0;
0056 }
0057
0058 namespace hitfit {
0059
0060   //**************************************************************************
0061   // Argument handling.
0062   //
0063
0064   Gentop_Args::Gentop_Args(const Defaults& defs)
0065       //
0066       // Purpose: Constructor.
0067       //
0068       // Inputs:
0069       //   defs -        The Defaults instance from which to initialize.
0070       //
0071       : _t_pt_mean(defs.get_float("t_pt_mean")),
0072         _mt(defs.get_float("mt")),
0073         _sigma_mt(defs.get_float("sigma_mt")),
0074         _mh(defs.get_float("mh")),
0075         _sigma_mh(defs.get_float("sigma_mh")),
0076         _recoil_pt_mean(defs.get_float("recoil_pt_mean")),
0077         _boost_sigma(defs.get_float("boost_sigma")),
0078         _m_boost(defs.get_float("m_boost")),
0079         _mb(defs.get_float("mb")),
0080         _sigma_mb(defs.get_float("sigma_mb")),
0081         _mw(defs.get_float("mw")),
0082         _sigma_mw(defs.get_float("sigma_mw")),
0083         _svx_tageff(defs.get_float("svx_tageff")),
0084         _smear(defs.get_bool("smear")),
0085         _smear_dir(defs.get_bool("smear_dir")),
0086         _muon(defs.get_bool("muon")),
0087         _ele_res_str(defs.get_string("ele_res_str")),
0088         _muo_res_str(defs.get_string("muo_res_str")),
0089         _jet_res_str(defs.get_string("jet_res_str")),
0090         _kt_res_str(defs.get_string("kt_res_str")) {}
0091
0092   double Gentop_Args::t_pt_mean() const
0093   //
0094   // Purpose: Return the t_pt_mean parameter.
0095   //          See the header for documentation.
0096   //
0097   {
0098     return _t_pt_mean;
0099   }
0100
0101   double Gentop_Args::mt() const
0102   //
0103   // Purpose: Return the mt parameter.
0104   //          See the header for documentation.
0105   //
0106   {
0107     return _mt;
0108   }
0109
0110   double Gentop_Args::sigma_mt() const
0111   //
0112   // Purpose: Return the sigma_mt parameter.
0113   //          See the header for documentation.
0114   //
0115   {
0116     return _sigma_mt;
0117   }
0118
0119   double Gentop_Args::svx_tageff() const
0120   //
0121   // Purpose: Return the svx_tageff parameter.
0122   //          See the header for documentation.
0123   //
0124   {
0125     return _svx_tageff;
0126   }
0127
0128   double Gentop_Args::mh() const
0129   //
0130   // Purpose: Return the mh parameter.
0131   //          See the header for documentation.
0132   //
0133   {
0134     return _mh;
0135   }
0136
0137   double Gentop_Args::sigma_mh() const
0138   //
0139   // Purpose: Return the sigma_mh parameter.
0140   //          See the header for documentation.
0141   //
0142   {
0143     return _sigma_mh;
0144   }
0145
0146   bool Gentop_Args::smear() const
0147   //
0148   // Purpose: Return the smear parameter.
0149   //          See the header for documentation.
0150   //
0151   {
0152     return _smear;
0153   }
0154
0155   bool Gentop_Args::smear_dir() const
0156   //
0157   // Purpose: Return the smear_dir parameter.
0158   //          See the header for documentation.
0159   //
0160   {
0161     return _smear_dir;
0162   }
0163
0164   bool Gentop_Args::muon() const
0165   //
0166   // Purpose: Return the muon parameter.
0167   //          See the header for documentation.
0168   //
0169   {
0170     return _muon;
0171   }
0172
0173   double Gentop_Args::recoil_pt_mean() const
0174   //
0175   // Purpose: Return the recoil_pt_mean parameter.
0176   //          See the header for documentation.
0177   //
0178   {
0179     return _recoil_pt_mean;
0180   }
0181
0182   double Gentop_Args::boost_sigma() const
0183   //
0184   // Purpose: Return the boost_sigma parameter.
0185   //          See the header for documentation.
0186   //
0187   {
0188     return _boost_sigma;
0189   }
0190
0191   double Gentop_Args::m_boost() const
0192   //
0193   // Purpose: Return the m_boost parameter.
0194   //          See the header for documentation.
0195   //
0196   {
0197     return _m_boost;
0198   }
0199
0200   double Gentop_Args::mb() const
0201   //
0202   // Purpose: Return the mb parameter.
0203   //          See the header for documentation.
0204   //
0205   {
0206     return _mb;
0207   }
0208
0209   double Gentop_Args::sigma_mb() const
0210   //
0211   // Purpose: Return the sigma_mb parameter.
0212   //          See the header for documentation.
0213   //
0214   {
0215     return _sigma_mb;
0216   }
0217
0218   double Gentop_Args::mw() const
0219   //
0220   // Purpose: Return the mw parameter.
0221   //          See the header for documentation.
0222   //
0223   {
0224     return _mw;
0225   }
0226
0227   double Gentop_Args::sigma_mw() const
0228   //
0229   // Purpose: Return the sigma_mw parameter.
0230   //          See the header for documentation.
0231   //
0232   {
0233     return _sigma_mw;
0234   }
0235
0236   std::string Gentop_Args::ele_res_str() const
0237   //
0238   // Purpose: Return the ele_res_str parameter.
0239   //          See the header for documentation.
0240   //
0241   {
0242     return _ele_res_str;
0243   }
0244
0245   std::string Gentop_Args::muo_res_str() const
0246   //
0247   // Purpose: Return the muo_res_str parameter.
0248   //          See the header for documentation.
0249   //
0250   {
0251     return _muo_res_str;
0252   }
0253
0254   std::string Gentop_Args::jet_res_str() const
0255   //
0256   // Purpose: Return the jet_res_str parameter.
0257   //          See the header for documentation.
0258   //
0259   {
0260     return _jet_res_str;
0261   }
0262
0263   std::string Gentop_Args::kt_res_str() const
0264   //
0265   // Purpose: Return the kt_res_str parameter.
0266   //          See the header for documentation.
0267   //
0268   {
0269     return _kt_res_str;
0270   }
0271
0272   //**************************************************************************
0273   // Internal helper functions.
0274   //
0275
0276   namespace {
0277
0278     /**
0279     @brief Generate a unit vector with \f$(\theta,\phi)\f$ uniformly
0280     distributed over the sphere.
0281
0282     @param engine The underlying random number generator.
0283  */
0284     Threevec rand_spher(CLHEP::HepRandomEngine& engine)
0285     //
0286     // Purpose: Return a unit vector with (theta, phi) uniformly distributed
0287     //          over a sphere.
0288     //
0289     // Inputs:
0290     //   engine -      The underlying RNG.
0291     //
0292     // Returns:
0293     //   The generated vector.
0294     //
0295     {
0296       CLHEP::RandFlat r(engine);
0297
0298       Threevec v;
0299
0300       double U = r.fire(0.0, 1.0);
0301       double V = r.fire(0.0, 1.0);
0302
0303       double theta = 2.0 * CLHEP::pi * U;
0304       double phi = acos(2 * V - 1.0);
0305
0306       double x = sin(theta) * cos(phi);
0307       double y = sin(theta) * sin(phi);
0308       double z = cos(theta);
0309
0310       v = Threevec(x, y, z);
0311
0312       return v.unit();
0313     }
0314
0315     /**
0316     @brief Given a direction, mass, and width, chose a mass from a
0317     Breit-Wigner distribution and return a four-momentum with the chosen
0318     mass and the specified direction.
0319
0320     @param p The direction three-momenta
0321     (not necessary to have unit magnitude).
0322
0323     @param m_true The mean for the Breit-Wigner distribution.
0324
0325     @param sigma The width for the Breit-Wigner distribution.
0326
0327     @param engine The underlying random number generator.
0328  */
0329     Fourvec make_massive(const Threevec& p, double m_true, double sigma, CLHEP::HepRandomEngine& engine)
0330     //
0331     // Purpose: Given a direction, mass, and width, choose a mass from a
0332     //          Breit-Wigner and return a 4-vector with the chosen mass
0333     //          and the specified direction.
0334     //
0335     // Inputs:
0336     //   p -           The direction.
0337     //   m_true -      The mean for the Breit-Wigner.
0338     //   sigma -       The width for the Breit-Wigner.
0339     //   engine -      The underlying RNG.
0340     //
0341     // Returns:
0342     //   The generated 4-vector.
0343     //
0344     {
0345       CLHEP::RandBreitWigner rbw(engine);
0346       double m = rbw.fire(m_true, sigma);
0347       return Fourvec(p, sqrt(m * m + p.mag2()));
0348     }
0349
0350     /**
0351     @brief Decay a particle with initial four-momentum \f$v\f$ into
0352     two particles with mass \f$m_{1}\f$ and \f$m_{2}\f$.
0353
0354     @param v The initial four-momentum.
0355
0356     @param m1 Mass of the first decay product.
0357
0358     @param m2 Mass of the second decay product.
0359
0360     @param engine The underlying random number generator.
0361
0362     @param vout1 Output, the outgoing four-momentum of the first decay product.
0363
0364     @param vout2 Output, the outgoing four-momentum of the second decay
0365     product.
0366  */
0367     void decay(const Fourvec& v, double m1, double m2, CLHEP::HepRandomEngine& engine, Fourvec& vout1, Fourvec& vout2)
0368     //
0369     // Purpose: v decays into two particles w/masses m1, m2.
0370     //
0371     // Inputs:
0372     //   v -           The incoming 4-vector.
0373     //   m1 -          Mass of the first decay product.
0374     //   m2 -          Mass of the second decay product.
0375     //   engine -      The underlying RNG.
0376     //
0377     // Outputs:
0378     //   vout1 -       Outgoing 4-vector of the first decay product.
0379     //   vout2 -       Outgoing 4-vector of the second decay product.
0380     //
0381     {
0382       // Construct a decay in the incoming particle's rest frame,
0383       // uniformly distributed in direction.
0384       Threevec p = rand_spher(engine);
0385       double m0 = v.m();
0386
0387       if (m1 + m2 > m0) {
0388         // What ya gonna do?
0389         double f = m0 / (m1 + m2);
0390         m1 *= f;
0391         m2 *= f;
0392       }
0393
0394       double m0_2 = m0 * m0;
0395       double m1_2 = m1 * m1;
0396       double m2_2 = m2 * m2;
0397
0398       // Calculate the 3-momentum of each particle in the decay frame.
0399       p *= 0.5 / m0 *
0400            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);
0401       double p2 = p.mag2();
0402
0403       vout1 = Fourvec(p, sqrt(p2 + m1_2));
0404       vout2 = Fourvec(-p, sqrt(p2 + m2_2));
0405
0406       // Boost out of the rest frame.
0407       vout1.boost(v.boostVector());
0408       vout2.boost(v.boostVector());
0409     }
0410
0411     /**
0412     @brief Generate a vector in a spherically uniform random direction,
0413     with transverse momentum \f$p_{T}\f$ drawn from an exponential distribution
0414     with mean <i>pt_mean</i>.
0415
0416     @param pt_mean The mean of the distribution.
0417
0418     @param engine The underlying random number generator.
0419  */
0420     Threevec rand_pt(double pt_mean, CLHEP::HepRandomEngine& engine)
0421     //
0422     // Purpose: Generate a vector in a (uniformly) random direction,
0423     //          with pt chosen from an exponential distribution with mean pt_mean.
0424     //
0425     // Inputs:
0426     //   pt_mean -     The mean of the distribution.
0427     //   engine -      The underlying RNG.
0428     //
0429     // Returns:
0430     //   The generated vector.
0431     //
0432     {
0433       CLHEP::RandExponential rexp(engine);
0434
0435       // A random direction.
0436       Threevec p = rand_spher(engine);
0437
0438       // Scale by random pt.
0439       p *= (rexp.fire(pt_mean) / p.perp());
0440
0441       return p;
0442     }
0443
0444     /**
0445     @brief Generate a random boost for the event.
0446
0447     @param args The parameter settings.
0448
0449     @param engine The underlying random number generator.
0450  */
0451     Fourvec rand_boost(const Gentop_Args& args, CLHEP::HepRandomEngine& engine)
0452     //
0453     // Purpose: Generate a random boost for the event.
0454     //
0455     // Inputs:
0456     //   args -        The parameter settings.
0457     //   engine -      The underlying RNG.
0458     //
0459     // Returns:
0460     //   The generated boost.
0461     //
0462     {
0463       CLHEP::RandExponential rexp(engine);
0464       CLHEP::RandFlat rflat(engine);
0465       CLHEP::RandGauss rgauss(engine);
0466
0467       // Boost in pt and z.
0468       Threevec p(1, 0, 0);
0469       p.rotateZ(rflat.fire(0, 2 * M_PI));
0470       p *= rexp.fire(args.recoil_pt_mean());
0471       p.setZ(rgauss.fire(0, args.boost_sigma()));
0472       return Fourvec(p, sqrt(p.mag2() + args.m_boost() * args.m_boost()));
0473     }
0474
0475     /**
0476     @brief Simulate SVX (Secondary Vertex) tagging.
0477
0478     @param args The parameter settings.
0479
0480     @param ev The event to tag.
0481
0482     @param engine The underlying random number generator.
0483  */
0484     void tagsim(const Gentop_Args& args, Lepjets_Event& ev, CLHEP::HepRandomEngine& engine)
0485     //
0486     // Purpose: Simulate SVX tagging.
0487     //
0488     // Inputs:
0489     //   args -        The parameter settings.
0490     //   ev -          The event to tag.
0491     //   engine -      The underlying RNG.
0492     //
0493     // Outputs:
0494     //   ev -          The event with tags filled in.
0495     //
0496     {
0497       CLHEP::RandFlat rflat(engine);
0498       for (std::vector<Lepjets_Event_Jet>::size_type i = 0; i < ev.njets(); i++) {
0499         int typ = ev.jet(i).type();
0500         if (typ == hadb_label || typ == lepb_label || typ == higgs_label) {
0501           if (rflat.fire() < args.svx_tageff())
0502             ev.jet(i).svx_tag() = true;
0503         }
0504       }
0505     }
0506
0507   }  // unnamed namespace
0508
0509   //**************************************************************************
0510   // External interface.
0511   //
0512
0513   Lepjets_Event gentop(const Gentop_Args& args, CLHEP::HepRandomEngine& engine)
0514   //
0515   // Purpose: Generate a ttbar -> ljets event.
0516   //
0517   // Inputs:
0518   //   args -        The parameter settings.
0519   //   engine -      The underlying RNG.
0520   //
0521   // Returns:
0522   //   The generated event.
0523   //
0524   {
0525     CLHEP::RandBreitWigner rbw(engine);
0526     CLHEP::RandGauss rgauss(engine);
0527
0528     // Get the t decay momentum in the ttbar rest frame.
0529     Threevec p = rand_pt(args.t_pt_mean(), engine);
0530
0531     // Make the t/tbar vectors.
0532     Fourvec lept = make_massive(p, args.mt(), args.sigma_mt(), engine);
0533     Fourvec hadt = make_massive(-p, args.mt(), args.sigma_mt(), engine);
0534
0535     // Boost the rest frame.
0536     Fourvec boost = rand_boost(args, engine);
0537     lept.boost(boost.boostVector());
0538     hadt.boost(boost.boostVector());
0539
0540     // Decay t -> b W, leptonic side.
0541     Fourvec lepb, lepw;
0542     double mlb = rgauss.fire(args.mb(), args.sigma_mb());
0543     double mlw = rbw.fire(args.mw(), args.sigma_mw());
0544     decay(lept, mlb, mlw, engine, lepb, lepw);
0545
0546     // Decay t -> b W, hadronic side.
0547     Fourvec hadb, hadw;
0548     double mhb = rgauss.fire(args.mb(), args.sigma_mb());
0549     double mhw = rbw.fire(args.mw(), args.sigma_mw());
0550     decay(hadt, mhb, mhw, engine, hadb, hadw);
0551
0552     // Decay W -> l nu.
0553     Fourvec lep, nu;
0554     decay(lepw, 0, 0, engine, lep, nu);
0555
0556     // Decay W -> qqbar.
0557     Fourvec q1, q2;
0558     decay(hadw, 0, 0, engine, q1, q2);
0559
0560     // Fill in the event.
0561     Lepjets_Event ev(0, ++next_evnum);
0562     Vector_Resolution lep_res(args.muon() ? args.muo_res_str() : args.ele_res_str());
0563     Vector_Resolution jet_res(args.jet_res_str());
0564     Resolution kt_res = (args.kt_res_str());
0565
0566     ev.add_lep(Lepjets_Event_Lep(lep, args.muon() ? muon_label : electron_label, lep_res));
0567
0568     ev.add_jet(Lepjets_Event_Jet(lepb, lepb_label, jet_res));
0569     ev.add_jet(Lepjets_Event_Jet(hadb, hadb_label, jet_res));
0570     ev.add_jet(Lepjets_Event_Jet(q1, hadw1_label, jet_res));
0571     ev.add_jet(Lepjets_Event_Jet(q2, hadw2_label, jet_res));
0572
0573     ev.met() = nu;
0574     ev.kt_res() = kt_res;
0575
0576     // Simulate SVX tagging.
0577     tagsim(args, ev, engine);
0578
0579     // Smear the event, if requested.
0580     if (args.smear())
0581       ev.smear(engine, args.smear_dir());
0582
0583     // Done!
0584     return ev;
0585   }
0586
0587   Lepjets_Event gentth(const Gentop_Args& args, CLHEP::HepRandomEngine& engine)
0588   //
0589   // Purpose: Generate a ttH -> ljets event.
0590   //
0591   // Inputs:
0592   //   args -        The parameter settings.
0593   //   engine -      The underlying RNG.
0594   //
0595   // Returns:
0596   //   The generated event.
0597   //
0598   {
0599     CLHEP::RandBreitWigner rbw(engine);
0600     CLHEP::RandGauss rgauss(engine);
0601
0602     // Generate three-vectors for two tops.
0603     Threevec p_t1 = rand_pt(args.t_pt_mean(), engine);
0604     Threevec p_t2 = rand_pt(args.t_pt_mean(), engine);
0605
0606     // Conserve momentum.
0607     Threevec p_h = -(p_t1 + p_t2);
0608
0609     // Construct the 4-vectors.
0610     Fourvec lept = make_massive(p_t1, args.mt(), args.sigma_mt(), engine);
0611     Fourvec hadt = make_massive(p_t2, args.mt(), args.sigma_mt(), engine);
0612     Fourvec higgs = make_massive(p_h, args.mh(), args.sigma_mh(), engine);
0613
0614     // Boost the rest frame.
0615     Fourvec boost = rand_boost(args, engine);
0616     lept.boost(boost.boostVector());
0617     hadt.boost(boost.boostVector());
0618     higgs.boost(boost.boostVector());
0619
0620     // Decay t -> b W, leptonic side.
0621     Fourvec lepb, lepw;
0622     decay(lept, rgauss.fire(args.mb(), args.sigma_mb()), rbw.fire(args.mw(), args.sigma_mw()), engine, lepb, lepw);
0623
0624     // Decay t -> b W, hadronic side.
0625     Fourvec hadb, hadw;
0626     decay(hadt, rgauss.fire(args.mb(), args.sigma_mb()), rbw.fire(args.mw(), args.sigma_mw()), engine, hadb, hadw);
0627
0628     // Decay W -> l nu.
0629     Fourvec lep, nu;
0630     decay(lepw, 0, 0, engine, lep, nu);
0631
0632     // Decay W -> qqbar.
0633     Fourvec q1, q2;
0634     decay(hadw, 0, 0, engine, q1, q2);
0635
0636     // Decay H -> bbbar.
0637     Fourvec hb1, hb2;
0638     decay(higgs, rgauss.fire(args.mb(), args.sigma_mb()), rgauss.fire(args.mb(), args.sigma_mb()), engine, hb1, hb2);
0639
0640     // Fill in the event.
0641     Lepjets_Event ev(0, ++next_evnum);
0642     Vector_Resolution lep_res(args.muon() ? args.muo_res_str() : args.ele_res_str());
0643     Vector_Resolution jet_res(args.jet_res_str());
0644     Resolution kt_res = (args.kt_res_str());
0645
0646     ev.add_lep(Lepjets_Event_Lep(lep, args.muon() ? muon_label : electron_label, lep_res));
0647
0648     ev.add_jet(Lepjets_Event_Jet(lepb, lepb_label, jet_res));
0649     ev.add_jet(Lepjets_Event_Jet(hadb, hadb_label, jet_res));
0650     ev.add_jet(Lepjets_Event_Jet(q1, hadw1_label, jet_res));
0651     ev.add_jet(Lepjets_Event_Jet(q2, hadw2_label, jet_res));
0652     ev.add_jet(Lepjets_Event_Jet(hb1, higgs_label, jet_res));
0653     ev.add_jet(Lepjets_Event_Jet(hb2, higgs_label, jet_res));
0654
0655     ev.met() = nu;
0656     ev.kt_res() = kt_res;
0657
0658     // Simulate SVX tagging.
0659     tagsim(args, ev, engine);
0660
0661     // Smear the event, if requested.
0662     if (args.smear())
0663       ev.smear(engine, args.smear_dir());
0664
0665     // Done!
0666     return ev;
0667   }
0668
0669 }  // namespace hitfit