Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:31:19

0001 //
0002 //
0003 // File: src/Top_Decaykin.cc
0004 // Purpose: Calculate some kinematic quantities for ttbar events.
0005 // Created: Jul, 2000, sss, based on run 1 mass analysis code.
0006 //
0007 // CMSSW File      : src/Top_Decaykin.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 Top_Decaykin.cc
0014 
0015     @brief A class to hold functions to calculate kinematic quantities
0016     of interest in \f$t\bar{t} \to \ell + 4 \mathrm{jets}\f$ events.
0017     See the documentation for the header file Top_Decaykin.h for details.
0018 
0019     @author Scott Stuart Snyder <snyder@bnl.gov>
0020 
0021     @par Creation date:
0022     Jul 2000.
0023 
0024     @par Modification History:
0025     Apr 2009: Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>:
0026     Imported to CMSSW.<br>
0027     Nov 2009: Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>:
0028     Added doxygen tags for automatic generation of documentation.
0029 
0030     @par Terms of Usage:
0031     With consent for the original author (Scott Snyder).
0032 
0033  */
0034 
0035 #include "TopQuarkAnalysis/TopHitFit/interface/Top_Decaykin.h"
0036 #include "TopQuarkAnalysis/TopHitFit/interface/Lepjets_Event.h"
0037 #include "TopQuarkAnalysis/TopHitFit/interface/fourvec.h"
0038 #include <cmath>
0039 #include <algorithm>
0040 #include <ostream>
0041 
0042 using std::abs;
0043 using std::ostream;
0044 using std::sqrt;
0045 using std::swap;
0046 
0047 namespace hitfit {
0048 
0049   namespace {
0050 
0051     /**
0052     @brief Sum the four-momenta of all leptons in an event.
0053 
0054     @param ev The event.
0055  */
0056     Fourvec leptons(const Lepjets_Event& ev)
0057     //
0058     // Purpose: Sum all leptons in EV.
0059     //
0060     // Inputs:
0061     //   ev -          The event.
0062     //
0063     // Returns:
0064     //   The sum of all leptons in EV.
0065     //
0066     {
0067       return (ev.sum(lepton_label) + ev.sum(electron_label) + ev.sum(muon_label));
0068     }
0069 
0070   }  // unnamed namespace
0071 
0072   bool Top_Decaykin::solve_nu_tmass(const Lepjets_Event& ev, double tmass, double& nuz1, double& nuz2)
0073   //
0074   // Purpose: Solve for the neutrino longitudinal z-momentum that makes
0075   //          the leptonic top have mass TMASS.
0076   //
0077   // Inputs:
0078   //   ev -          The event to solve.
0079   //   tmass -       The desired top mass.
0080   //
0081   // Outputs:
0082   //   nuz1 -        First solution (smaller absolute value).
0083   //   nuz2 -        Second solution.
0084   //
0085   // Returns:
0086   //   True if there was a real solution.  False if there were only
0087   //   imaginary solutions.  (In that case, we just set the imaginary
0088   //   part to zero.)
0089   //
0090   {
0091     bool discrim_flag = true;
0092 
0093     const Fourvec& vnu = ev.met();
0094     Fourvec cprime = leptons(ev) + ev.sum(lepb_label);
0095     double alpha1 = tmass * tmass - cprime.m2();
0096     double a = 2 * 4 * (cprime.z() * cprime.z() - cprime.e() * cprime.e());
0097     double alpha = alpha1 + 2 * (cprime.x() * vnu.x() + cprime.y() * vnu.y());
0098     double b = 4 * alpha * cprime.z();
0099     double c = alpha * alpha - 4 * cprime.e() * cprime.e() * vnu.vect().perp2();
0100     double d = b * b - 2 * a * c;
0101     if (d < 0) {
0102       discrim_flag = false;
0103       d = 0;
0104     }
0105 
0106     double dd = sqrt(d);
0107     nuz1 = (-b + dd) / a;
0108     nuz2 = (-b - dd) / a;
0109     if (abs(nuz1) > abs(nuz2))
0110       swap(nuz1, nuz2);
0111 
0112     return discrim_flag;
0113   }
0114 
0115   bool Top_Decaykin::solve_nu_tmass(
0116       const Lepjets_Event& ev, double tmass, double& re_nuz1, double& im_nuz1, double& re_nuz2, double& im_nuz2)
0117   //
0118   // Purpose: Solve for the neutrino longitudinal z-momentum that makes
0119   //          the leptonic top have mass TMASS, including the imaginary
0120   //          component of the z-component of the neutrino'somentum.
0121   //
0122   // Inputs:
0123   //   ev -          The event to solve.
0124   //   tmass -       The desired top mass.
0125   //
0126   // Outputs:
0127   //   re_nuz1 -     Real component of the first solution.
0128   //   im_nuz1 -     Imaginary component of the first solution (in the lower half of
0129   //                 the complex plane).
0130   //   re_nuz2 -     Real component of the second solution.
0131   //   im_nuz2 -     Imaginary component of the second solution (in the upper half of
0132   //                 the complex plane).
0133   //
0134   // Returns:
0135   //   True if there was a real solution.  False if there were only
0136   //   complex solutions.
0137   //
0138   {
0139     bool discrim_flag = true;
0140 
0141     const Fourvec& vnu = ev.met();
0142     Fourvec cprime = leptons(ev) + ev.sum(lepb_label);
0143     double alpha1 = tmass * tmass - cprime.m2();
0144     double a = 2 * 4 * (cprime.z() * cprime.z() - cprime.e() * cprime.e());
0145     // Haryo's note: Here a is equivalent to '2a' in the quadratic
0146     // equation ax^2 + bx + c = 0
0147     double alpha = alpha1 + 2 * (cprime.x() * vnu.x() + cprime.y() * vnu.y());
0148     double b = 4 * alpha * cprime.z();
0149     double c = alpha * alpha - 4 * cprime.e() * cprime.e() * vnu.vect().perp2();
0150     double d = b * b - 2 * a * c;
0151     if (d < 0) {
0152       discrim_flag = false;
0153     }
0154 
0155     if (discrim_flag) {
0156       re_nuz1 = (-b + sqrt(d)) / a;
0157       im_nuz1 = 0;
0158       re_nuz2 = (-b - sqrt(d)) / a;
0159       im_nuz2 = 0;
0160       if (abs(re_nuz1) > abs(re_nuz2)) {
0161         swap(re_nuz1, re_nuz2);
0162       }
0163 
0164     } else {
0165       // Take absolute value of the imaginary component of nuz, in case
0166       // a is negative, before multiplying by +1 or -1 to get the upper-half
0167       // or lower-half imaginary value.
0168       re_nuz1 = -b / a;
0169       im_nuz1 = -fabs(sqrt(-d) / a);
0170       re_nuz2 = -b / a;
0171       im_nuz2 = fabs(sqrt(-d) / a);
0172     }
0173 
0174     return discrim_flag;
0175   }
0176 
0177   bool Top_Decaykin::solve_nu(const Lepjets_Event& ev, double wmass, double& nuz1, double& nuz2)
0178   //
0179   // Purpose: Solve for the neutrino longitudinal z-momentum that makes
0180   //          the leptonic W have mass WMASS.
0181   //
0182   // Inputs:
0183   //   ev -          The event to solve.
0184   //   wmass -       The desired W mass.
0185   //
0186   // Outputs:
0187   //   nuz1 -        First solution (smaller absolute value).
0188   //   nuz2 -        Second solution.
0189   //
0190   // Returns:
0191   //   True if there was a real solution.  False if there were only
0192   //   imaginary solutions.  (In that case, we just set the imaginary
0193   //   part to zero.)
0194   //
0195   {
0196     bool discrim_flag = true;
0197 
0198     const Fourvec& vnu = ev.met();
0199     Fourvec vlep = leptons(ev);
0200 
0201     double x = vlep.x() * vnu.x() + vlep.y() * vnu.y() + wmass * wmass / 2;
0202     double a = vlep.z() * vlep.z() - vlep.e() * vlep.e();
0203     double b = 2 * x * vlep.z();
0204     double c = x * x - vnu.perp2() * vlep.e() * vlep.e();
0205 
0206     double d = b * b - 4 * a * c;
0207     if (d < 0) {
0208       d = 0;
0209       discrim_flag = false;
0210     }
0211 
0212     nuz1 = (-b + sqrt(d)) / 2 / a;
0213     nuz2 = (-b - sqrt(d)) / 2 / a;
0214     if (abs(nuz1) > abs(nuz2))
0215       swap(nuz1, nuz2);
0216 
0217     return discrim_flag;
0218   }
0219 
0220   bool Top_Decaykin::solve_nu(
0221       const Lepjets_Event& ev, double wmass, double& re_nuz1, double& im_nuz1, double& re_nuz2, double& im_nuz2)
0222   //
0223   // Purpose: Solve for the neutrino longitudinal z-momentum that makes
0224   //          the leptonic W have mass WMASS, including the imaginary
0225   //          component of the z-component of the neutrino'somentum.
0226   //
0227   // Inputs:
0228   //   ev -          The event to solve.
0229   //   wmass -       The desired W mass.
0230   //
0231   // Outputs:
0232   //   re_nuz1 -     Real component of the first solution.
0233   //   im_nuz1 -     Imaginary component of the first solution  (in the lower half of
0234   //                 the complex plane).
0235   //   re_nuz2 -     Real component of the second solution.
0236   //   im_nuz2 -     Imaginary component of the second solution  (in the upper half of
0237   //                 the complex plane).
0238   //
0239   // Returns:
0240   //   True if there was a real solution.  False if there were only
0241   //   complex solutions.
0242   //x
0243   {
0244     bool discrim_flag = true;
0245 
0246     const Fourvec& vnu = ev.met();
0247     Fourvec vlep = leptons(ev);
0248 
0249     double x = vlep.x() * vnu.x() + vlep.y() * vnu.y() + wmass * wmass / 2;
0250     double a = vlep.z() * vlep.z() - vlep.e() * vlep.e();
0251     double b = 2 * x * vlep.z();
0252     double c = x * x - vnu.perp2() * vlep.e() * vlep.e();
0253 
0254     double d = b * b - 4 * a * c;
0255     if (d < 0) {
0256       discrim_flag = false;
0257     }
0258 
0259     if (discrim_flag) {
0260       re_nuz1 = (-b + sqrt(d)) / 2 / a;
0261       im_nuz1 = 0.0;
0262       re_nuz2 = (-b - sqrt(d)) / 2 / a;
0263       im_nuz2 = 0.0;
0264       if (fabs(re_nuz1) > fabs(re_nuz2)) {
0265         swap(re_nuz1, re_nuz2);
0266       }
0267 
0268     } else {
0269       // Take absolute value of the imaginary component of nuz, in case
0270       // a is negative, before multiplying by +1 or -1 to get the upper-half
0271       // or lower-half imaginary value.
0272 
0273       re_nuz1 = -b / 2 / a;
0274       im_nuz1 = -fabs(sqrt(-d) / a);
0275       re_nuz2 = -b / 2 / a;
0276       im_nuz2 = fabs(sqrt(-d) / a);
0277     }
0278 
0279     return discrim_flag;
0280   }
0281 
0282   Fourvec Top_Decaykin::hadw(const Lepjets_Event& ev)
0283   //
0284   // Purpose: Sum up the appropriate 4-vectors to find the hadronic W.
0285   //
0286   // Inputs:
0287   //   ev -          The event.
0288   //
0289   // Returns:
0290   //   The hadronic W.
0291   //
0292   {
0293     return (ev.sum(hadw1_label) + ev.sum(hadw2_label));
0294   }
0295 
0296   Fourvec Top_Decaykin::hadw1(const Lepjets_Event& ev)
0297   //
0298   // Purpose:
0299   //
0300   // Inputs:
0301   //   ev -          The event.
0302   //
0303   // Returns:
0304   //   The higher-pT hadronic jet from W
0305   //
0306   {
0307     return ev.sum(hadw1_label);
0308   }
0309 
0310   Fourvec Top_Decaykin::hadw2(const Lepjets_Event& ev)
0311   //
0312   // Purpose:
0313   //
0314   // Inputs:
0315   //   ev -          The event.
0316   //
0317   // Returns:
0318   //   The lower-pT hadronic jet from W
0319   //
0320   //
0321   {
0322     return ev.sum(hadw2_label);
0323   }
0324 
0325   Fourvec Top_Decaykin::lepw(const Lepjets_Event& ev)
0326   //
0327   // Purpose: Sum up the appropriate 4-vectors to find the leptonic W.
0328   //
0329   // Inputs:
0330   //   ev -          The event.
0331   //
0332   // Returns:
0333   //   The leptonic W.
0334   //
0335   {
0336     return (leptons(ev) + ev.met());
0337   }
0338 
0339   Fourvec Top_Decaykin::hadt(const Lepjets_Event& ev)
0340   //
0341   // Purpose: Sum up the appropriate 4-vectors to find the hadronic t.
0342   //
0343   // Inputs:
0344   //   ev -          The event.
0345   //
0346   // Returns:
0347   //   The hadronic t.
0348   //
0349   {
0350     return (ev.sum(hadb_label) + hadw(ev));
0351   }
0352 
0353   Fourvec Top_Decaykin::lept(const Lepjets_Event& ev)
0354   //
0355   // Purpose: Sum up the appropriate 4-vectors to find the leptonic t.
0356   //
0357   // Inputs:
0358   //   ev -          The event.
0359   //
0360   // Returns:
0361   //   The leptonic t.
0362   //
0363   {
0364     return (ev.sum(lepb_label) + lepw(ev));
0365   }
0366 
0367   ostream& Top_Decaykin::dump_ev(std::ostream& s, const Lepjets_Event& ev)
0368   //
0369   // Purpose: Print kinematic information for EV.
0370   //
0371   // Inputs:
0372   //   s -           The stream to which to write.
0373   //   ev -          The event to dump.
0374   //
0375   // Returns:
0376   //   The stream S.
0377   //
0378   {
0379     s << ev;
0380     Fourvec p;
0381 
0382     p = lepw(ev);
0383     s << "lepw " << p << " " << p.m() << "\n";
0384     p = lept(ev);
0385     s << "lept " << p << " " << p.m() << "\n";
0386     p = hadw(ev);
0387     s << "hadw " << p << " " << p.m() << "\n";
0388     p = hadt(ev);
0389     s << "hadt " << p << " " << p.m() << "\n";
0390 
0391     return s;
0392   }
0393 
0394   double Top_Decaykin::cos_theta_star(const Fourvec& fermion, const Fourvec& W, const Fourvec& top)
0395   //
0396   // Purpose: Calculate cos theta star in top decay
0397   //
0398   // Inputs:
0399   //   fermion -     The four momentum of fermion from W
0400   //   W -           The four momentum of W from top
0401   //   top -         The four momentum of top
0402   // Returns:
0403   //   cos theta star
0404   //
0405   {
0406     if (W.isLightlike() || W.isSpacelike()) {
0407       return 100.0;
0408     }
0409 
0410     CLHEP::HepBoost BoostWCM(W.findBoostToCM());
0411 
0412     CLHEP::Hep3Vector boost_v3fermion = BoostWCM(fermion).vect();
0413     CLHEP::Hep3Vector boost_v3top = BoostWCM(top).vect();
0414 
0415     double costhetastar = boost_v3fermion.cosTheta(-boost_v3top);
0416 
0417     return costhetastar;
0418   }
0419 
0420   double Top_Decaykin::cos_theta_star(const Lepjets_Event& ev)
0421   //
0422   // Purpose: Calculate lepton cos theta star in top decay
0423   //
0424   // Inputs:
0425   //   ev -          A lepton+jets event
0426   // Returns:
0427   //   cos theta star of lepton
0428   //
0429   {
0430     return cos_theta_star(leptons(ev), lepw(ev), lept(ev));
0431   }
0432 
0433   double Top_Decaykin::cos_theta_star_lept(const Lepjets_Event& ev)
0434   //
0435   // Purpose: Calculate lepton cos theta star in top decay
0436   //
0437   // Inputs:
0438   //   ev -          A lepton+jets event
0439   // Returns:
0440   //   cos theta star of lepton
0441   //
0442   {
0443     return cos_theta_star(ev);
0444   }
0445 
0446   double Top_Decaykin::cos_theta_star_hadt(const Lepjets_Event& ev)
0447   //
0448   // Purpose: Calculate absolute value of cos theta star of
0449   //          one of the hadronic W jet from hadronic top.
0450   //
0451   // Inputs:
0452   //   ev -          A lepton+jets event
0453   // Returns:
0454   //   absolute value of cos theta star
0455   //
0456   {
0457     return fabs(cos_theta_star(hadw1(ev), hadw(ev), hadt(ev)));
0458   }
0459 
0460 }  // namespace hitfit