Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:31:16

0001 //
0002 //
0003 // File: src/fourvec.cc
0004 // Purpose: Define 3- and 4-vector types for the hitfit package, and
0005 //          supply a few additional operations.
0006 // Created: Jul, 2000, sss, based on run 1 mass analysis code.
0007 //
0008 // CMSSW File      : src/fourvec.cc
0009 // Original Author : Scott Stuart Snyder <snyder@bnl.gov> for D0
0010 // Imported to CMSSW by Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>
0011 //
0012 
0013 /**
0014     @file fourvec.cc
0015 
0016     @brief Define three-vector and four-vector classes for the HitFit
0017     package, and supply a few additional operations. See the documentation
0018     for the header file fourvec.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 from the original author (Scott Snyder).
0033  */
0034 
0035 #include "TopQuarkAnalysis/TopHitFit/interface/fourvec.h"
0036 #include <cmath>
0037 
0038 using std::atan;
0039 using std::atan2;
0040 using std::cos;
0041 using std::exp;
0042 using std::fabs;
0043 using std::log;
0044 using std::sin;
0045 using std::sqrt;
0046 using std::tan;
0047 
0048 namespace {  // unnamed namespace
0049 
0050   double cal_th(double theta, double z)
0051   //
0052   // Purpose: Helper for deteta.  Find `calorimeter theta' given
0053   //          physics theta and z-vertex.  Copied from run 1 code.
0054   //
0055   // Inputs:
0056   //   theta -       Physics theta.
0057   //   z -           Z-vertex.
0058   //
0059   // Returns:
0060   //   Calorimeter theta.
0061   //
0062   {
0063     const double R_CC = 91.6;
0064     const double Z_EC = 178.9;
0065     const double BIGG = 1e8;
0066 
0067     double tanth;
0068     if (fabs(cos(theta)) < 1 / BIGG)
0069       tanth = BIGG * sin(theta);
0070     else
0071       tanth = tan(theta);
0072 
0073     double z_cc = R_CC / tanth + z;
0074 
0075     if (fabs(z_cc) < Z_EC)
0076       theta = atan2(R_CC, z_cc);
0077 
0078     else {
0079       double zz = Z_EC;
0080       if (tanth < 0)
0081         zz = -zz;
0082       double r_ec = (zz - z) * tanth;
0083       theta = atan2(r_ec, zz);
0084     }
0085 
0086     if (theta < 0)
0087       theta += 2 * M_PI;
0088     return theta;
0089   }
0090 
0091 }  // unnamed namespace
0092 
0093 namespace hitfit {
0094 
0095   void adjust_p_for_mass(Fourvec& v, double mass)
0096   //
0097   // Purpose: Adjust the 3-vector part of V (leaving the energy unchanged)
0098   //          so that it has mass MASS.  (Provided that is possible.)
0099   //
0100   // Inputs:
0101   //   v -           The 4-vector to scale.
0102   //   mass -        The desired mass of the 4-vector.
0103   //
0104   // Outputs:
0105   //   v -           The scaled 4-vector.
0106   //
0107   {
0108     CLHEP::Hep3Vector vect = v.vect();
0109     double old_p2 = vect.mag2();
0110     if (old_p2 == 0)
0111       return;
0112     double new_p2 = v.e() * v.e() - mass * mass;
0113     if (new_p2 < 0)
0114       new_p2 = 0;
0115     vect *= sqrt(new_p2 / old_p2);
0116     v.setVect(vect);
0117   }
0118 
0119   void adjust_e_for_mass(Fourvec& v, double mass)
0120   //
0121   // Purpose: Adjust the energy component of V (leaving the 3-vector part
0122   //          unchanged) so that it has mass MASS.
0123   //
0124   // Inputs:
0125   //   v -           The 4-vector to scale.
0126   //   mass -        The desired mass of the 4-vector.
0127   //
0128   // Outputs:
0129   //   v -           The scaled 4-vector.
0130   //
0131   {
0132     v.setE(sqrt(v.vect().mag2() + mass * mass));
0133   }
0134 
0135   void rottheta(Fourvec& v, double theta)
0136   //
0137   // Purpose: Rotate V through polar angle THETA.
0138   //
0139   // Inputs:
0140   //   v -           The 4-vector to rotate.
0141   //   theta -       The rotation angle.
0142   //
0143   // Outputs:
0144   //   v -           The rotated 4-vector.
0145   //
0146   {
0147     double s = sin(theta), c = cos(theta);
0148     double old_pt = v.perp();
0149     double new_pt = old_pt * c - v.z() * s;
0150     v.setZ(old_pt * s + v.z() * c);
0151 
0152     v.setX(v.x() * new_pt / old_pt);
0153     v.setY(v.y() * new_pt / old_pt);
0154   }
0155 
0156   void roteta(Fourvec& v, double eta)
0157   //
0158   // Purpose: Rotate a Fourvec through a polar angle such that
0159   //          its pseudorapidity changes by ETA.
0160   //
0161   // Inputs:
0162   //   v -           The 4-vector to rotate.
0163   //   eta -         The rotation angle.
0164   //
0165   // Outputs:
0166   //   v -           The rotated 4-vector.
0167   //
0168   {
0169     double theta1 = v.theta();
0170     double eta1 = theta_to_eta(theta1);
0171     double eta2 = eta1 + eta;
0172     double theta2 = eta_to_theta(eta2);
0173 
0174     rottheta(v, theta1 - theta2);
0175   }
0176 
0177   double eta_to_theta(double eta)
0178   //
0179   // Purpose: Convert psuedorapidity to polar angle.
0180   //
0181   // Inputs:
0182   //   eta -         Pseudorapidity.
0183   //
0184   // Returns:
0185   //   Polar angle.
0186   //
0187   {
0188     return 2 * atan(exp(-eta));
0189   }
0190 
0191   double theta_to_eta(double theta)
0192   //
0193   // Purpose: Convert polar angle to psuedorapidity.
0194   //
0195   // Inputs:
0196   //   theta -       Polar angle.
0197   //
0198   // Returns:
0199   //   Pseudorapidity.
0200   //
0201   {
0202     return -log(tan(theta / 2));
0203   }
0204 
0205   double deteta(const Fourvec& v, double zvert)
0206   //
0207   // Purpose: Get the detector eta (D0-specific).
0208   //
0209   // Inputs:
0210   //   v -           Vector on which to operate.
0211   //   zvert -       Z-vertex.
0212   //
0213   // Returns:
0214   //   Detector eta of V.
0215   //
0216   {
0217     return theta_to_eta(cal_th(v.theta(), zvert));
0218   }
0219 
0220   double phidiff(double phi)
0221   //
0222   // Purpose: Handle wraparound for a difference in azimuthal angles.
0223   //
0224   // Inputs:
0225   //   phi -         Azimuthal angle.
0226   //
0227   // Returns:
0228   //   PHI normalized to the range -pi .. pi.
0229   //
0230   {
0231     while (phi < -M_PI)
0232       phi += 2 * M_PI;
0233     while (phi > M_PI)
0234       phi -= 2 * M_PI;
0235     return phi;
0236   }
0237 
0238   double delta_r(const Fourvec& a, const Fourvec& b)
0239   //
0240   // Purpose: Find the distance in R between two four-vectors.
0241   //
0242   // Inputs:
0243   //   a -           First four-vector.
0244   //   b -           Second four-vector.
0245   //
0246   // Returns:
0247   //   the distance in R between A and B.
0248   //
0249   {
0250     double deta = a.pseudoRapidity() - b.pseudoRapidity();
0251     double dphi = phidiff(a.phi() - b.phi());
0252     return sqrt(deta * deta + dphi * dphi);
0253   }
0254 
0255 }  // namespace hitfit