TopHitFit/src/Resolution.cc

0001 //
0002 //
0003 // File: src/Resolution.cc
0004 // Purpose: Calculate resolutions for a quantity.
0005 // Created: Jul, 2000, sss, based on run 1 mass analysis code.
0006 //
0007 // CMSSW File      : src/Resolution.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 Resolution.cc
0014
0015     @brief Calculate and represent resolution for a physical quantity.
0016
0017     @author Scott Stuart Snyder <snyder@bnl.gov>
0018
0019     @par Creation date:
0020     Jul 2000.
0021
0022     @par Modification History:
0023     Apr 2009: Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>:
0024     Imported to CMSSW.<br>
0025     Nov 2009: Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>:
0026     Added doxygen tags for automatic generation of documentation.
0027
0028     @par Terms of Usage:
0029     With consent for the original author (Scott Snyder).
0030
0031  */
0032
0033 #include "TopQuarkAnalysis/TopHitFit/interface/Resolution.h"
0034 #include "CLHEP/Random/RandGauss.h"
0035 #include <cmath>
0036 #include <iostream>
0037 #include <cctype>
0038 #include <cstdlib>
0039
0040 using std::isdigit;
0041 using std::isspace;
0042 using std::ostream;
0043 using std::sqrt;
0044 using std::string;
0045 #ifndef __GNUC__
0046 using std::atof;
0047 #endif
0048
0049 namespace {
0050
0051   bool get_field(string s, string::size_type i, double& x)
0052   //
0053   // Purpose: Scan string S starting at position I.
0054   //          Find the value of the first floating-point number we
0055   //          find there.
0056   //
0057   // Inputs:
0058   //   s -           The string to scan.
0059   //   i -           Starting character position in the string.
0060   //
0061   // Outputs:
0062   //   x -           The value of the number we found.
0063   //
0064   // Returns:
0065   //   True if we found something that looks like a number, false otherwise.
0066   //
0067   {
0068     string::size_type j = i;
0069     while (j < s.size() && s[j] != ',' && !isdigit(s[j]) && s[j] != '.')
0070       ++j;
0071     if (j < s.size() && (isdigit(s[j]) || s[j] == '.')) {
0072       x = atof(s.c_str() + j);
0073       return true;
0074     }
0075     return false;
0076   }
0077
0078 }  // unnamed namespace
0079
0080 namespace hitfit {
0081
0082   Resolution::Resolution(std::string s /*= ""*/)
0083       //
0084       // Purpose: Constructor.
0085       //
0086       // Inputs:
0087       //   s -           String encoding the resolution parameters, as described
0088       //                 in the comments in the header.
0089       //
0090       : _resolution_exponent(0) {
0091     _inverse = false;
0092     _constant_sigma = 0;
0093     _resolution_sigma = 0;
0094     _noise_sigma = 0;
0095
0096     // Skip spaces.
0097     double x;
0098     string::size_type i = 0;
0099     while (i < s.size() && isspace(s[i]))
0100       ++i;
0101
0102     // Check for the inverse flag.
0103     if (s[i] == '-') {
0104       _inverse = true;
0105       ++i;
0106     } else if (s[i] == '+') {
0107       ++i;
0108     }
0109
0110     // Get the constant term.
0111     if (get_field(s, i, x))
0112       _constant_sigma = x;
0113     i = s.find(',', i);
0114
0115     // Look for a resolution term.
0116     if (i != string::npos) {
0117       ++i;
0118       if (get_field(s, i, x))
0119         _resolution_sigma = x;
0120
0121       // Look for a noise term.
0122       i = s.find(',', i);
0123       if (i != string::npos) {
0124         if (get_field(s, i + 1, x))
0125           _noise_sigma = x;
0126       }
0127     }
0128   }
0129
0130   Resolution::Resolution(double C, double R, double m, double N, bool inverse /*= false*/)
0131       : _constant_sigma(C), _resolution_sigma(R), _resolution_exponent(m), _noise_sigma(N), _inverse(inverse) {}
0132
0133   Resolution::Resolution(double res, bool inverse /*= false*/)
0134       //
0135       // Purpose: Constructor, to give a constant resolution.
0136       //          I.e., sigma() will always return RES.
0137       //
0138       // Inputs:
0139       //   res -         The resolution value.
0140       //   inverse -     The inverse flag.
0141       //
0142       : _constant_sigma(0), _resolution_sigma(0), _resolution_exponent(0), _noise_sigma(res), _inverse(inverse) {}
0143
0144   bool Resolution::inverse() const
0145   //
0146   // Purpose: Return the inverse flag.
0147   //
0148   // Returns:
0149   //   The inverse flag.
0150   //
0151   {
0152     return _inverse;
0153   }
0154
0155   double Resolution::C() const { return _constant_sigma; }
0156
0157   double Resolution::R() const { return _resolution_sigma; }
0158
0159   double Resolution::m() const { return _resolution_exponent; }
0160
0161   double Resolution::N() const { return _noise_sigma; }
0162
0163   double Resolution::sigma(double p) const
0164   //
0165   // Purpose: Return the uncertainty for a momentum P.
0166   //
0167   // Inputs:
0168   //    p -          The momentum
0169   //
0170   // Returns:
0171   //   The uncertainty for a momentum P.
0172   //
0173   {
0174     p = std::fabs(p);
0175     if (_inverse)
0176       p = 1 / p;
0177
0178     return sqrt((_constant_sigma * _constant_sigma * p + _resolution_sigma * _resolution_sigma) * p +
0179                 _noise_sigma * _noise_sigma);
0180   }
0181
0182   double Resolution::pick(double x, double p, CLHEP::HepRandomEngine& engine) const
0183   //
0184   // Purpose: Given a value X, measured for an object with momentum P,
0185   //          pick a new value from a Gaussian distribution
0186   //          described by this resolution --- with mean X and width sigma(P).
0187   //
0188   // Inputs:
0189   //    x -          The quantity value (distribution mean).
0190   //    p -          The momentum, for calculating the width.
0191   //    engine -     The underlying RNG.
0192   //
0193   // Returns:
0194   //   A smeared value of X.
0195   //
0196   {
0197     CLHEP::RandGauss gen(engine);
0198     if (_inverse)
0199       return 1 / gen.fire(1 / x, sigma(p));
0200     else
0201       return gen.fire(x, sigma(p));
0202   }
0203
0204   /**
0205     @brief Output stream operator, print the content of this Resolution
0206     to an output stream.
0207
0208     @param s The stream to which to write.
0209
0210     @param r The instance of Resolution to be printed.
0211  */
0212   std::ostream& operator<<(std::ostream& s, const Resolution& r)
0213   //
0214   // Purpose: Dump this object to S.
0215   //
0216   // Inputs:
0217   //    s -          The stream to which to write.
0218   //    r -          The object to dump.
0219   //
0220   // Returns:
0221   //   The stream S.
0222   //
0223   {
0224     if (r._inverse)
0225       s << "-";
0226     s << r._constant_sigma << "," << r._resolution_sigma << "," << r._noise_sigma;
0227     return s;
0228   }
0229
0230 }  // namespace hitfit