Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 //
0002 //
0003 // File: hitfit/Chisq_Constrainer.h
0004 // Purpose: Minimize a chisq subject to a set of constraints.
0005 //          Based on the SQUAW algorithm.
0006 // Created: Jul, 2000, sss, based on run 1 mass analysis code.
0007 //
0008 // For full details on the algorithm, see
0009 //
0010 //    @phdthesis{sssthesis,
0011 //      author =       "Scott Snyder",
0012 //      school =       "State University of New York at Stony Brook",
0013 //      month =        may,
0014 //      year =         "1995 (unpublished)"}
0015 //    @comment{  note =         "available from {\tt http://www-d0.fnal.gov/publications\_talks/thesis/ snyder/thesis-ps.html}"
0016 //    }
0017 //
0018 // CMSSW File      : interface/Chisq_Constrainer.h
0019 // Original Author : Scott Stuart Snyder <snyder@bnl.gov> for D0
0020 // Imported to CMSSW by Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>
0021 
0022 /**
0023 
0024     @file Chisq_Constrainer.h
0025 
0026     @brief Minimize a \f$\chi^{2}\f$ subject to a set of constraints,
0027     based on the SQUAW algorithm.
0028 
0029     For the full details of the algorithm see:
0030     - Scott Snyder, <i>Measurement of the top quark mass at D0</i>, PhD thesis,
0031     State University of New York at Stony Brook, (1995).  Available from
0032     <a href=
0033     http://lss.fnal.gov/archive/thesis/1900/fermilab-thesis-1995-27.shtml>
0034     http://lss.fnal.gov/archive/thesis/1900/fermilab-thesis-1995-27.shtml</a>.
0035 
0036     - B. Abbott <i>et al.</i> (D0 Collaboration), <i>Direct Measurement of
0037     the Top Quark Mass at D0</i>,
0038     <a href="http://prola.aps.org/abstract/PRD/v58/i5/e052001">
0039     Phys. Rev. <b>D58</b>, 052001  (1998)</a>,
0040     <a href="http://arxiv.org/abs/hep-ex/9801025">arXiv:hep-ex/9801025</a>.
0041 
0042     - O. I. Dahl, T. B. Day, F. T. Solnitz, and M. L. Gould, <i>SQUAW Kinematic
0043     Fitting Program</i>, Group A Programming  Note P-126, Lawrence Radiation
0044     Laboratory (1968).  Available from
0045     <a href="http://alvarezphysicsmemos.lbl.gov">
0046     Luis Alvarez Physics Memos</a>.
0047 
0048     @par Creation date.
0049     July 2000.
0050 
0051     @author
0052     Scott Stuart Snyder <snyder@bnl.gov>.
0053 
0054     @par Modification History:
0055     Apr 2009: Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>:
0056     Imported to CMSSW.<br>
0057     Oct 2009: Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>:
0058     Added Doxygen tags for automatic generation of documentation.
0059 
0060     @par Terms of Usage:
0061     With consent from the original author (Scott Snyder).
0062 
0063 */
0064 
0065 #ifndef HITFIT_CHISQ_CONSTRAINER_H
0066 #define HITFIT_CHISQ_CONSTRAINER_H
0067 
0068 #include "TopQuarkAnalysis/TopHitFit/interface/Base_Constrainer.h"
0069 #include "TopQuarkAnalysis/TopHitFit/interface/matutil.h"
0070 #include <iosfwd>
0071 
0072 namespace hitfit {
0073 
0074   class Defaults;
0075 
0076   /**
0077      @class Chisq_Constrainer_Args
0078      @brief Hold on to parameters for the Chisq_Constrainer class.
0079 */
0080   class Chisq_Constrainer_Args
0081   //
0082   // Purpose: Hold on to parameters for the Chisq_Constrainer class.
0083   //
0084   // Parameters controlling the operation of the fitter:
0085   //   bool printfit      - If true, print a trace of the fit to cout.
0086   //   bool use_G         - If true, check the chisq formula by computing
0087   //                        chisq directly from G.  This requires that G_i
0088   //                        be invertable.
0089   //
0090   // Parameters affecting the fit:
0091   //   float constraint_sum_eps - Convergence threshold for sum of constraints.
0092   //   float chisq_diff_eps - onvergence threshold for change in chisq.
0093   //   int maxit          - Maximum number of iterations permitted.
0094   //   int max_cut        - Maximum number of cut steps permitted.
0095   //   float cutsize      - Fraction by which to cut steps.
0096   //   float min_tot_cutsize - Smallest fractional cut step permitted.
0097   //
0098   // Parameters affecting testing modes:
0099   //   float chisq_test_eps - When use_G is true, the maximum relative
0100   //                          difference permitted between the two chisq
0101   //                          calculations.
0102   //
0103   {
0104   public:
0105     // Constructor.  Initialize from a Defaults object.
0106 
0107     /**
0108        @brief Constructor, creates an instance of Chisq_Constrainer_Args
0109        from a Defaults object.
0110 
0111        @param defs The Defaults instance from which to instantiate.
0112        The instance
0113        must contain variables of type and name:
0114        - bool <i>use_G</i>.
0115        - bool <i>printfit</i>.
0116        - double <i>constraint_sum_eps</i>.
0117        - double <i>chisq_diff_eps</i>.
0118        - int <i>maxit</i>.
0119        - int <i>max_cut</i>.
0120        - double <i>cutsize</i>.
0121        - double <i>min_tot_cutsize</i>.
0122        - double <i>chisq_test_eps</i>.
0123    */
0124     Chisq_Constrainer_Args(const Defaults& defs);
0125 
0126     // Retrieve parameter values.
0127 
0128     /**
0129      Return the <i>printfit</i> parameter.
0130    */
0131     bool printfit() const;
0132 
0133     /**
0134      Return the <i>use_G</i> parameter.
0135    */
0136     bool use_G() const;
0137 
0138     /**
0139      Return the <i>constraint_sum_eps</i> parameter.
0140    */
0141     double constraint_sum_eps() const;
0142 
0143     /**
0144      Return the <i>chisq_diff_eps</i> parameter.
0145    */
0146     double chisq_diff_eps() const;
0147 
0148     /**
0149      Return the <i>maxit</i> parameter.
0150    */
0151     unsigned maxit() const;
0152 
0153     /**
0154      Return the <i>max_cut</i> parameter.
0155    */
0156     unsigned max_cut() const;
0157 
0158     /**
0159      Return the <i>cutsize</i> parameter.
0160    */
0161     double cutsize() const;
0162 
0163     /**
0164      Return the <i>min_tot_cutsize</i> parameter.
0165    */
0166     double min_tot_cutsize() const;
0167 
0168     /**
0169      Return the <i>chisq_test_eps</i> parameter.
0170    */
0171     double chisq_test_eps() const;
0172 
0173     // Arguments for subobjects.
0174 
0175     /**
0176      Return the argument for the Base_Constrainer class.
0177    */
0178     const Base_Constrainer_Args& base_constrainer_args() const;
0179 
0180   private:
0181     // Hold on to parameter values.
0182 
0183     /**
0184      If true, print a trace of the fit to std::cout.
0185    */
0186     bool _printfit;
0187 
0188     /**
0189      If true, check the \f$\chi^{2}\f$ formula by computing the \f$\chi^{2}\f$
0190      directly from \f${\bf G}\f$. This requires that \f${\bf G}_{i}\f$ be
0191      invertible.
0192    */
0193     bool _use_G;
0194 
0195     /**
0196      Convergence threshold for sum of constraints.
0197    */
0198     double _constraint_sum_eps;
0199 
0200     /**
0201      Convergence threshold for change in \f$\chi^{2}\f$.
0202    */
0203     double _chisq_diff_eps;
0204 
0205     /**
0206      Maxium number of iterations permitted.
0207    */
0208     int _maxit;
0209 
0210     /**
0211      Maximum number of cut steps permitted.
0212    */
0213     int _max_cut;
0214 
0215     /**
0216      Fraction by which to cut steps.
0217    */
0218     double _cutsize;
0219 
0220     /**
0221      Smallest fractional cut step permitted.
0222    */
0223     double _min_tot_cutsize;
0224 
0225     /**
0226      When <i>use_G</i> is true, the maximum relative difference between
0227      the \f$\chi^{2}\f$ calculations.
0228    */
0229     double _chisq_test_eps;
0230 
0231     /**
0232      Parameters for the underlying base class Base_Constrainer.
0233    */
0234     const Base_Constrainer_Args _base_constrainer_args;
0235   };
0236 
0237   //*************************************************************************
0238 
0239   /**
0240     @class Chisq_Constrainer
0241     @brief Minimize a \f$\chi^{2}\f$ subject to a set of constraints.  Based
0242     on the SQUAW algorithm.
0243  */
0244   class Chisq_Constrainer
0245       //
0246       // Purpose: Minimize a chisq subject to a set of constraints.
0247       //          Based on the SQUAW algorithm.
0248       //
0249       : public Base_Constrainer {
0250   public:
0251     // Constructor, destructor.
0252     // ARGS holds the parameter settings for this instance.
0253 
0254     /**
0255      Constructor.  Create an instance of Chisq_Constrainer from a
0256      Chisq_Constrainer_Args object.
0257      @param args The parameter settings for this instance.
0258 
0259    */
0260     Chisq_Constrainer(const Chisq_Constrainer_Args& args);
0261 
0262     /**
0263      Destructor.
0264    */
0265     ~Chisq_Constrainer() override {}
0266 
0267     // Do the fit.
0268     // Call the number of well-measured variables Nw, the number of
0269     // poorly-measured variables Np, and the number of constraints Nc.
0270     // Inputs:
0271     //   CONSTRAINT_CALCULATOR is the object that will be used to evaluate
0272     //     the constraints.
0273     //   XM(Nw) and YM(Np) are the measured values of the well- and
0274     //     poorly-measured variables, respectively.
0275     //   X(Nw) and Y(Np) are the starting values for the fit.
0276     //   G_I(Nw,Nw) is the error matrix for the well-measured variables.
0277     //   Y(Np,Np) is the inverse error matrix for the poorly-measured variables.
0278     //
0279     // Outputs:
0280     //   X(Nw) and Y(Np) is the point at the minimum.
0281     //   PULLX(Nw) and PULLY(Np) are the pull quantities.
0282     //   Q(Nw,Nw), R(Np,Np), and S(Nw,Np) are the final error matrices
0283     //     between all the variables.
0284     //
0285     // The return value is the final chisq.  Returns a value < 0 if the
0286     // fit failed to converge.
0287 
0288     /**
0289      @ brief Do a constrained fit.  Call the number of well-measured variables
0290      <i>Nw</i>, the number of poorly-measured variables <i>Np</i>,
0291      and the number of constraints <i>Nc</i>.
0292 
0293      @param constraint_calculator The object that will be used to evaluate
0294      the constraints.
0295 
0296      @param xm Measured values of the well-measured variables, has dimension
0297      <i>Nw</i>.
0298 
0299      @param x  Before the fit: starting value of the well-measured variables
0300      for the fit, has dimension <i>Nw</i>.
0301      After the fit: final values of the well-measured variables.
0302 
0303      @param ym Measured values of the poorly-measured variables, has
0304      dimension <i>Np</i>.
0305 
0306      @param y  Before the fit: starting value of the poorly-measured variables
0307      for the fit, has dimension <i>Np</i>.
0308      After the fit: final values of the poorly-measured variables.
0309 
0310      @param G_i Error matrix for the well-measured variables, has dimension
0311      <i>Nw,Nw</i>.
0312 
0313      @param Y Inverse error matrix for the poorly-measured variables, has
0314      dimension <i>Np,Np</i>.
0315 
0316      @param pullx Pull quantities for the well-measured variables, has
0317      dimension <i>Nw</i>.
0318 
0319      @param pully Pull quantities for the poorly-measured variables, has
0320      dimension <i>Np</i>.
0321 
0322      @param Q Error matrix for the well-measured variables, has dimension
0323      <i>Nw,Nw</i>.
0324 
0325      @param R Error matrix for the poorly-measured variables, has dimension
0326      <i>Np,Np</i>.
0327 
0328      @param S Error matrix for the correlation between well-measured variables
0329      and poorly-measured variables, has dimension <i>Nw,Np</i>.
0330 
0331      @par Input:
0332      <i>constraint_calculator, xm, x, ym, y, G_i, y.</i>.
0333 
0334      @par Output:
0335      <i>x, y, pullx, pully, Q, R, S.</i>
0336      @par Return:
0337      \f$\chi^{2}\f$ of the fit.  Should returns a negative value if the fit
0338      does not converge.
0339    */
0340     double fit(Constraint_Calculator& constraint_calculator,
0341                const Column_Vector& xm,
0342                Column_Vector& x,
0343                const Column_Vector& ym,
0344                Column_Vector& y,
0345                const Matrix& G_i,
0346                const Diagonal_Matrix& Y,
0347                Column_Vector& pullx,
0348                Column_Vector& pully,
0349                Matrix& Q,
0350                Matrix& R,
0351                Matrix& S) override;
0352 
0353     // Print out any internal state to S.
0354     /**
0355      @brief Print the state of this instance of Chisq_Constrainer.
0356 
0357      @param s The output stream to which the output is sent.
0358  */
0359     std::ostream& print(std::ostream& s) const override;
0360 
0361   private:
0362     // Parameter settings.
0363     /**
0364          Parameter settings for this instance of Chisq_Constrainer.
0365    */
0366     const Chisq_Constrainer_Args _args;
0367   };
0368 
0369 }  // namespace hitfit
0370 
0371 #endif  // not HITFIT_CHISQ_CONSTRAINER_H