CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/TopQuarkAnalysis/TopHitFit/src/Vector_Resolution.cc

Go to the documentation of this file.
00001 //
00002 // $Id: Vector_Resolution.cc,v 1.1 2011/05/26 09:47:00 mseidel Exp $
00003 //
00004 // File: src/Vector_Resolution.cc
00005 // Purpose: Calculate resolutions in p, phi, eta.
00006 // Created: Jul, 2000, sss, based on run 1 mass analysis code.
00007 //
00008 // CMSSW File      : src/Vector_Resolution.cc
00009 // Original Author : Scott Stuart Snyder <snyder@bnl.gov> for D0
00010 // Imported to CMSSW by Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>
00011 //
00012 
00013 
00038 #include "TopQuarkAnalysis/TopHitFit/interface/Vector_Resolution.h"
00039 #include "TopQuarkAnalysis/TopHitFit/interface/fourvec.h"
00040 #include <cmath>
00041 #include <ostream>
00042 #include <cctype>
00043 
00044 
00045 using std::ostream;
00046 using std::string;
00047 using std::isspace;
00048 
00049 
00050 namespace {
00051 
00060 string field (string s, int i)
00061 //
00062 // Purpose: Extract the Ith slash-separated field from S.
00063 //
00064 // Inputs:
00065 //   s -           The string to analyze.
00066 //   i -           Field number (starting with 0).
00067 //
00068 // Returns:
00069 //   Field I of S.
00070 //
00071 {
00072   string::size_type pos = 0;
00073   while (i > 0) {
00074     pos = s.find ('/', pos);
00075     if (pos == string::npos)
00076       return "";
00077     ++pos;
00078     --i;
00079   }
00080 
00081   string::size_type pos2 = s.find ('/', pos+1);
00082   if (pos2 != string::npos)
00083     pos2 = pos2 - pos;
00084 
00085   while (pos < pos2 && isspace (s[pos]))
00086     ++pos;
00087 
00088   return s.substr (pos, pos2);
00089 }
00090 
00091 
00092 } // unnamed namespace
00093 
00094 
00095 namespace hitfit {
00096 
00097 
00098 Vector_Resolution::Vector_Resolution ()
00099 //
00100 // Purpose: Default constructor.  Create a vector resolution object with
00101 //          infinite precision.v
00102 //
00103   : _p_res   ("0,0,0"),
00104     _eta_res ("0,0,0"),
00105     _phi_res ("0,0,0"),
00106     _use_et  (false)
00107 
00108 {
00109 }
00110 
00111 
00112 Vector_Resolution::Vector_Resolution (std::string s)
00113 //
00114 // Purpose: Constructor.
00115 //
00116 // Inputs:
00117 //   s -           String encoding the resolution parameters, as described
00118 //                 in the comments in the header.
00119 //
00120   : _p_res   (field (s, 0)),
00121     _eta_res (field (s, 1)),
00122     _phi_res (field (s, 2)),
00123     _use_et  (field (s, 3) == "et")
00124 {
00125 }
00126 
00127 
00128 Vector_Resolution::Vector_Resolution (const Resolution& p_res,
00129                                       const Resolution& eta_res,
00130                                       const Resolution& phi_res,
00131                                       bool use_et /*= false*/)
00132 //
00133 // Purpose: Constructor from individual resolution objects.
00134 //
00135 // Inputs:
00136 //   p_res -       Momentum resolution object.
00137 //   eta_res -     Eta resolution object.
00138 //   phi_res -     Phi resolution object.
00139 //   use_et -      If true, momentum resolution is based on pt, not p.
00140 //
00141   : _p_res   (p_res),
00142     _eta_res (eta_res),
00143     _phi_res (phi_res),
00144     _use_et  (use_et)
00145 {
00146 }
00147 
00148 
00149 const Resolution& Vector_Resolution::p_res () const
00150 //
00151 // Purpose: Return the momentum resolution object.
00152 //
00153 // Returns:
00154 //   The momentum resolution object.
00155 //
00156 {
00157   return _p_res;
00158 }
00159 
00160 
00161 const Resolution& Vector_Resolution::eta_res () const
00162 //
00163 // Purpose: Return the eta resolution object.
00164 //
00165 // Returns:
00166 //   The eta resolution object.
00167 //
00168 {
00169   return _eta_res;
00170 }
00171 
00172 
00173 const Resolution& Vector_Resolution::phi_res () const
00174 //
00175 // Purpose: Return the phi resolution object.
00176 //
00177 // Returns:
00178 //   The phi resolution object.
00179 //
00180 {
00181   return _phi_res;
00182 }
00183 
00184 
00185 bool Vector_Resolution::use_et () const
00186 //
00187 // Purpose: Return the use_et flag.
00188 //
00189 // Returns:
00190 //   The use_et flag.
00191 //
00192 {
00193   return _use_et;
00194 }
00195 
00196 
00197 namespace {
00198 
00199 
00211 double find_sigma (const Fourvec& v, const Resolution& res, bool use_et)
00212 //
00213 // Purpose: Helper for *_sigma functions below.
00214 //
00215 // Inputs:
00216 //   v -           The 4-vector.
00217 //   res -         The resolution object.
00218 //   use_et -      Use_et flag.
00219 //
00220 // Returns:
00221 //   The result of res.sigma() (not corrected for e/et difference).
00222 //
00223 {
00224   double ee = use_et ? v.perp() : v.e(); // ??? is perp() correct here?
00225   return res.sigma (ee);
00226 }
00227 
00228 
00229 } // unnamed namespace
00230 
00231 
00232 double Vector_Resolution::p_sigma (const Fourvec& v) const
00233 //
00234 // Purpose: Calculate momentum resolution for 4-vector V.
00235 //
00236 // Inputs:
00237 //   v -           The 4-vector.
00238 //
00239 // Returns:
00240 //   The momentum resolution for 4-vector V.
00241 //
00242 {
00243   double sig = find_sigma (v, _p_res, _use_et);
00244   if (_use_et) {
00245     if(_p_res.inverse()){
00246       sig *= v.perp() / v.e();
00247     }else{
00248       sig *= v.e() / v.perp();
00249     }
00250   }
00251   return sig;
00252 }
00253 
00254 
00255 double Vector_Resolution::eta_sigma (const Fourvec& v) const
00256 //
00257 // Purpose: Calculate eta resolution for 4-vector V.
00258 //
00259 // Inputs:
00260 //   v -           The 4-vector.
00261 //
00262 // Returns:
00263 //   The eta resolution for 4-vector V.
00264 //
00265 {
00266   return find_sigma (v, _eta_res, _use_et);
00267 }
00268 
00269 
00270 double Vector_Resolution::phi_sigma (const Fourvec& v) const
00271 //
00272 // Purpose: Calculate phi resolution for 4-vector V.
00273 //
00274 // Inputs:
00275 //   v -           The 4-vector.
00276 //
00277 // Returns:
00278 //   The phi resolution for 4-vector V.
00279 //
00280 {
00281   return find_sigma (v, _phi_res, _use_et);
00282 }
00283 
00284 
00285 
00286 namespace {
00287 
00288 
00300 void smear_eta (Fourvec& v, double ee,
00301                 const Resolution& res, CLHEP::HepRandomEngine& engine)
00302 //
00303 // Purpose: Smear the eta direction of V.
00304 //
00305 // Inputs:
00306 //   v -           The 4-vector to smear.
00307 //   ee -          Energy for sigma calculation.
00308 //   res -         The resolution object, giving the amount of smearing.
00309 //   engine -      The underlying RNG.
00310 //
00311 // Outputs:
00312 //   v -           The smeared 4-vector.
00313 //
00314 {
00315   double rot = res.pick (0, ee, engine);
00316   roteta (v, rot);
00317 }
00318 
00319 
00331 void smear_phi (Fourvec& v, double ee,
00332                 const Resolution& res, CLHEP::HepRandomEngine& engine)
00333 //
00334 // Purpose: Smear the phi direction of V.
00335 //
00336 // Inputs:
00337 //   v -           The 4-vector to smear.
00338 //   ee -          Energy for sigma calculation.
00339 //   res -         The resolution object, giving the amount of smearing.
00340 //   engine -      The underlying RNG.
00341 //
00342 // Outputs:
00343 //   v -           The smeared 4-vector.
00344 // 
00345 {
00346   double rot = res.pick (0, ee, engine);
00347   v.rotateZ (rot);
00348 }
00349 
00350 
00351 } // unnamed namespace
00352 
00353 
00354 void Vector_Resolution::smear (Fourvec& v,
00355                                CLHEP::HepRandomEngine& engine,
00356                                bool do_smear_dir /*= false*/) const
00357 //
00358 // Purpose: Smear a 4-vector according to the resolutions.
00359 //
00360 // Inputs:
00361 //   v -           The 4-vector to smear.
00362 //   engine -      The underlying RNG.
00363 //   do_smear_dir- If false, only smear the energy.
00364 //
00365 // Outputs:
00366 //   v -           The smeared 4-vector.
00367 //
00368 {
00369   double ee = _use_et ? v.perp() : v.e(); // ??? is perp() correct?
00370   v *= _p_res.pick (ee, ee, engine) / ee;
00371 
00372   if (do_smear_dir) {
00373     smear_eta (v, ee, _eta_res, engine);
00374     smear_phi (v, ee, _phi_res, engine);
00375   }
00376 }
00377 
00378 
00387 std::ostream& operator<< (std::ostream& s, const Vector_Resolution& r)
00388 //
00389 // Purpose: Dump this object to S.
00390 //
00391 // Inputs:
00392 //    s -          The stream to which to write.
00393 //    r -          The object to dump.
00394 //
00395 // Returns:
00396 //   The stream S.
00397 //   
00398 {
00399   s << r._p_res << "/ " << r._eta_res << "/ " << r._phi_res;
00400   if (r._use_et)
00401     s << "/et";
00402   s << "\n";
00403   return s;
00404 }
00405 
00406 
00407 } // namespace hitfit