Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 //
0002 //
0003 // File: src/Fourvec_Constrainer.cc
0004 // Purpose: Do a kinematic fit for a set of 4-vectors, given a set
0005 //          of mass constraints.
0006 // Created: Jul, 2000, sss, based on run 1 mass analysis code.
0007 //
0008 // CMSSW File      : src/Fourvec_Constrainer.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 
0015     @file Fourvec_Constrainer.cc
0016 
0017     @brief Do a kinematic fit for a set of four-vectors, given a set
0018     of mass constraints. Contains definitions of helper class
0019     Fourvec_Constraint_Calculator and some helper functions.
0020     See the documentation for the header file Fourvec_Constrainer.h for
0021     details.
0022 
0023     @author Scott Stuart Snyder <snyder@bnl.gov>
0024 
0025     @par Creation date:
0026     Jul 2000.
0027 
0028     @par Modification History:
0029     Apr 2009: Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>:
0030     Imported to CMSSW.<br>
0031     Nov 2009: Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>:
0032     Added doxygen tags for automatic generation of documentation.
0033 
0034     @par Terms of Usage:
0035     With consent for the original author (Scott Snyder).
0036 
0037 */
0038 
0039 #include "TopQuarkAnalysis/TopHitFit/interface/Fourvec_Constrainer.h"
0040 #include "TopQuarkAnalysis/TopHitFit/interface/Fourvec_Event.h"
0041 #include "TopQuarkAnalysis/TopHitFit/interface/Pair_Table.h"
0042 #include "TopQuarkAnalysis/TopHitFit/interface/Chisq_Constrainer.h"
0043 #include "TopQuarkAnalysis/TopHitFit/interface/matutil.h"
0044 #include "TopQuarkAnalysis/TopHitFit/interface/Defaults.h"
0045 #include <cmath>
0046 #include <iostream>
0047 
0048 using std::cos;
0049 using std::exp;
0050 using std::ostream;
0051 using std::sin;
0052 using std::sqrt;
0053 using std::vector;
0054 
0055 namespace hitfit {
0056 
0057   //*************************************************************************
0058   // Argument handling.
0059   //
0060 
0061   Fourvec_Constrainer_Args::Fourvec_Constrainer_Args(const Defaults& defs)
0062       //
0063       // Purpose: Constructor.
0064       //
0065       // Inputs:
0066       //   defs -        The Defaults instance from which to initialize.
0067       //
0068       : _use_e(defs.get_bool("use_e")),
0069         _e_com(defs.get_float("e_com")),
0070         _ignore_met(defs.get_bool("ignore_met")),
0071         _chisq_constrainer_args(defs) {}
0072 
0073   bool Fourvec_Constrainer_Args::use_e() const
0074   //
0075   // Purpose: Return the use_e parameter.
0076   //          See the header for documentation.
0077   //
0078   {
0079     return _use_e;
0080   }
0081 
0082   double Fourvec_Constrainer_Args::e_com() const
0083   //
0084   // Purpose: Return the e_com parameter.
0085   //          See the header for documentation.
0086   //
0087   {
0088     return _e_com;
0089   }
0090 
0091   bool Fourvec_Constrainer_Args::ignore_met() const
0092   //
0093   // Purpose: Return the ignore_met parameter.
0094   //          See the header for documentation.
0095   //
0096   {
0097     return _ignore_met;
0098   }
0099 
0100   const Chisq_Constrainer_Args& Fourvec_Constrainer_Args::chisq_constrainer_args() const
0101   //
0102   // Purpose: Return the contained subobject parameters.
0103   //
0104   {
0105     return _chisq_constrainer_args;
0106   }
0107 
0108   //*************************************************************************
0109   // Variable layout.
0110   //
0111   // We need to map the quantities we fit onto the vectors of well- and
0112   // poorly-measured quantities.
0113   //
0114 
0115   //
0116   // The well-measured variables consist of three variables for each
0117   // object.  If we are using transverse momentum constraints,
0118   // these fill be followed by the two cartesian components of kt.
0119   //
0120   // Each object is represented by three variables: the momentum (or 1/p
0121   // if the muon flag was set), and the two spherical angles, phi and eta.
0122   // Here is how they're ordered.
0123   //
0124   /**
0125     Offset indices for the component of four-momentum variables.
0126  */
0127   typedef enum { p_offs = 0, phi_offs = 1, eta_offs = 2 } Offsets;
0128 
0129   /**
0130     Offset indices for the components of missing transverse energy
0131     (or \f$k_{T}\f$) variables.
0132  */
0133   typedef enum { x_offs = 0, y_offs = 1 } Kt_Offsets;
0134 
0135   //
0136   // If there is a neutrino, then it is at index 1 of the poorly-measured
0137   // set (otherwise, that set is empty).
0138   //
0139   /**
0140     If there is a neutrino, then it is at index 1 of the vector of
0141     poorly-measured variables.
0142  */
0143   typedef enum { nu_z = 1 } Unmeasured_Variables;
0144 
0145   namespace {
0146 
0147     /**
0148     @brief Helper function: Return the starting variable index for object
0149     number <i>i</i>.
0150 
0151     @param i The object's index.
0152  */
0153     int obj_index(int i)
0154     //
0155     // Purpose: Return the starting variable index for object I.
0156     //
0157     // Inputs:
0158     //   i -           The object index.
0159     //
0160     // Returns:
0161     //   The index in the well-measured set of the first variable
0162     //   for object I.
0163     //
0164     {
0165       return i * 3 + 1;
0166     }
0167 
0168   }  // unnamed namespace
0169 
0170   //*************************************************************************
0171   // Object management.
0172   //
0173 
0174   Fourvec_Constrainer::Fourvec_Constrainer(const Fourvec_Constrainer_Args& args)
0175       //
0176       // Purpose: Constructor.
0177       //
0178       // Inputs:
0179       //   args -        The parameter settings for this instance.
0180       //
0181       : _args(args) {}
0182 
0183   void Fourvec_Constrainer::add_constraint(std::string s)
0184   //
0185   // Purpose: Specify an additional constraint S for the problem.
0186   //          The format for S is described in the header.
0187   //
0188   // Inputs:
0189   //   s -           The constraint to add.
0190   //
0191   {
0192     _constraints.push_back(Constraint(s));
0193   }
0194 
0195   void Fourvec_Constrainer::mass_constraint(std::string s)
0196   //
0197   // Purpose: Specify the combination of objects that will be returned by
0198   //          constrain() as the mass.  The format of S is the same as for
0199   //          normal constraints.  The LHS specifies the mass to calculate;
0200   //          the RHS should be zero.
0201   //          This should only be called once.
0202   //
0203   // Inputs:
0204   //   s -           The constraint defining the mass.
0205   //
0206   {
0207     assert(_mass_constraint.empty());
0208     _mass_constraint.push_back(Constraint(s));
0209   }
0210 
0211   /**
0212     @brief Output stream operator, print the content of this
0213     Fourvec_Constrainer to an output stream.
0214 
0215     @param s The output stream to which to write.
0216 
0217     @param c The instance of Fourvec_Constrainer to be printed.
0218 
0219 */
0220   std::ostream& operator<<(std::ostream& s, const Fourvec_Constrainer& c)
0221   //
0222   // Purpose: Print the object to S.
0223   //
0224   // Inputs:
0225   //   s -           The stream to which to write.
0226   //   c -           The object to write.
0227   //
0228   // Returns:
0229   //   The stream S.
0230   //
0231   {
0232     s << "Constraints: (e_com = " << c._args.e_com() << ") ";
0233     if (c._args.use_e())
0234       s << "(E)";
0235     s << "\n";
0236 
0237     for (std::vector<Constraint>::size_type i = 0; i < c._constraints.size(); i++)
0238       s << "  " << c._constraints[i] << "\n";
0239 
0240     if (!c._mass_constraint.empty()) {
0241       s << "Mass constraint:\n";
0242       s << c._mass_constraint[0] << "\n";
0243     }
0244     return s;
0245   }
0246 
0247   //*************************************************************************
0248   // Event packing and unpacking.
0249   //
0250 
0251   namespace {
0252 
0253     /**
0254     @brief Helper function: For all objects in the Fourvec_Event instance
0255     <i>ev</i>, adjust their four-momenta to have their requested masses.
0256 
0257     @param ev The event on which to operate.
0258 
0259     @param use_e_flag If TRUE, keep the energy and scale the
0260     three-momentum.<br>
0261     If FALSE, keep the three-momentum and scale the energy.
0262  */
0263     void adjust_fourvecs(Fourvec_Event& ev, bool use_e_flag)
0264     //
0265     // Purpose: For all objects in EV, adjust their 4-momenta
0266     //          to have their requested masses.
0267     //
0268     // Inputs:
0269     //   ev -          The event on which to operate.
0270     //   use_e_flag -  If true, keep E and scale 3-momentum.
0271     //                 If false, keep the 3-momentum and scale E.
0272     //
0273     {
0274       int nobjs = ev.nobjs();
0275       for (int i = 0; i < nobjs; i++) {
0276         const FE_Obj& obj = ev.obj(i);
0277         Fourvec p = obj.p;
0278         if (use_e_flag)
0279           adjust_p_for_mass(p, obj.mass);
0280         else
0281           adjust_e_for_mass(p, obj.mass);
0282         ev.set_obj_p(i, p);
0283       }
0284     }
0285 
0286     /**
0287     @brief Helper function: Convert object at index <i>ndx</i> from its
0288     representation in the vector of well-measured variables <i>c</i>
0289     to a four-momentum.
0290 
0291     @param c The vector of well-measured variables.
0292 
0293     @param ndx The index of the object in which we are interested.
0294 
0295     @param obj The object from the instance of Fourvec_Event.
0296 
0297     @param use_e_flag If TRUE, then we are using energy <i>E</i> as the
0298     fit variable.<br>
0299     If FALSE, we are using magnitude of three-momentum <i>p</i> as the
0300     fit variable.
0301  */
0302     Fourvec get_p_eta_phi_vec(const Column_Vector& c, int ndx, const FE_Obj& obj, bool use_e_flag)
0303     //
0304     // Purpose: Convert object NDX from its representation in the set
0305     //          of well-measured variables C to a 4-vector.
0306     //
0307     // Inputs:
0308     //   c -           The vector of well-measured variables.
0309     //   ndx -         The index of the object in which we're interested.
0310     //   obj -         The object from the Fourvec_Event.
0311     //   use_e_flag -  If true, we're using E as the fit variable, otherwise p.
0312     //
0313     // Returns:
0314     //   The object's 4-momentum.
0315     //
0316     {
0317       // Get the energy and momentum of the object.
0318       double e, p;
0319 
0320       if (use_e_flag) {
0321         // We're using E as a fit variable.  Get it directly.
0322         e = c(ndx + p_offs);
0323 
0324         // Take into account the muon case.
0325         if (obj.muon_p)
0326           e = 1 / e;
0327 
0328         // Find the momentum given the energy.
0329         if (obj.mass == 0)
0330           p = e;
0331         else {
0332           double xx = e * e - obj.mass * obj.mass;
0333           if (xx >= 0)
0334             p = sqrt(xx);
0335           else
0336             p = 0;
0337         }
0338       } else {
0339         // We're using P as a fit variable.  Fetch it.
0340         p = c(ndx + p_offs);
0341 
0342         // Take into account the muon case.
0343         if (obj.muon_p)
0344           p = 1 / p;
0345 
0346         // Find the energy given the momentum.
0347         e = (obj.mass == 0 ? p : sqrt(obj.mass * obj.mass + p * p));
0348       }
0349 
0350       // Get angular variables.
0351       double phi = c(ndx + phi_offs);
0352       double eta = c(ndx + eta_offs);
0353       if (fabs(eta) > 50) {
0354         // Protect against ridiculously large etas
0355         eta = eta > 0 ? 50 : -50;
0356       }
0357       double exp_eta = exp(eta);
0358       double iexp_eta = 1 / exp_eta;
0359       double sin_theta = 2 / (exp_eta + iexp_eta);
0360       double cos_theta = (exp_eta - iexp_eta) / (exp_eta + iexp_eta);
0361 
0362       // Form the 4-momentum.
0363       return Fourvec(p * sin_theta * cos(phi), p * sin_theta * sin(phi), p * cos_theta, e);
0364     }
0365 
0366     /**
0367     @brief Helper function: Initialize the variables in the vector of
0368     well-measured variables <i>c</i> describing object at index <i>ndx</i>
0369     from its Fourvec_Event representation <i>obj</i>.
0370 
0371     @param obj The object from the instance of Fourvec_Event.
0372 
0373     @param c The vector of well-measured variables.
0374 
0375     @param ndx The index of the object in which we are interested.
0376 
0377     @param use_e_flag If TRUE, then we are using energy <i>E</i> as the
0378     fit variable.<br>
0379     If FALSE, we are using magnitude of three-momentum <i>p</i> as the
0380     fit variable.
0381  */
0382     void set_p_eta_phi_vec(const FE_Obj& obj, Column_Vector& c, int ndx, bool use_e_flag)
0383     //
0384     // Purpose: Initialize the variables in the well-measured set C describing
0385     //          object NDX from its Fourvec_Event representation OBJ.
0386     //
0387     // Inputs:
0388     //   obj -         The object from the Fourvec_Event.
0389     //   c -           The vector of well-measured variables.
0390     //   ndx -         The index of the object in which we're interested.
0391     //   use_e_flag -  If true, we're using E as the fit variable, otherwise p.
0392     //
0393     //
0394     {
0395       if (use_e_flag)
0396         c(ndx + p_offs) = obj.p.e();
0397       else
0398         c(ndx + p_offs) = obj.p.vect().mag();
0399 
0400       if (obj.muon_p)
0401         c(ndx + p_offs) = 1 / c(ndx + p_offs);
0402       c(ndx + phi_offs) = obj.p.phi();
0403       c(ndx + eta_offs) = obj.p.pseudoRapidity();
0404     }
0405 
0406     /**
0407     @brief Helper function: Construct a 3 by 3 error matrix for an object.
0408 
0409     @param p_sig The uncertainty on the momentum \f$p\f$.
0410 
0411     @param phi_sig The uncertainty on the azimuthal angle \f$\phi\f$.
0412 
0413     @param eta_sig The uncertainty on the pseudorapidity \f$\eta\f$.
0414  */
0415     Matrix error_matrix(double p_sig, double phi_sig, double eta_sig)
0416     //
0417     // Purpose: Set up the 3x3 error matrix for an object.
0418     //
0419     // Inputs:
0420     //   p_sig -       The momentum uncertainty.
0421     //   phi_sig -     The phi uncertainty.
0422     //   eta_sig -     The eta uncertainty.
0423     //
0424     // Returns:
0425     //   The object's error matrix.
0426     //
0427     {
0428       Matrix err(3, 3, 0);
0429       err(1 + p_offs, 1 + p_offs) = p_sig * p_sig;
0430       err(1 + phi_offs, 1 + phi_offs) = phi_sig * phi_sig;
0431       err(1 + eta_offs, 1 + eta_offs) = eta_sig * eta_sig;
0432       return err;
0433     }
0434 
0435     /**
0436     @brief Helper function: Take the information from a
0437     Fourvec_Event instance <i>ev</i> and pack it into vectors
0438     of well- and poorly-measured variables.  Also set up
0439     the error matrices.
0440 
0441     @param ev The event to pack.
0442 
0443     @param use_e_flag If TRUE, then we are using energy <i>E</i> as the
0444     fit variable.<br>
0445     If FALSE, we are using magnitude of three-momentum <i>p</i> as the
0446     fit variable.
0447 
0448     @param use_kt_flag If TRUE, then we are also pack \f$k){T}\f$
0449     variables.<br>
0450     If FALSE, then we are not.
0451 
0452     @param xm The vector of well-measured variables.
0453 
0454     @param ym The vector of poorly-measured variables.
0455 
0456     @param G_i The error matrix for well-measured variables.
0457 
0458     @param Y Inverse error matrix for poorly-measured variables.
0459  */
0460     void pack_event(const Fourvec_Event& ev,
0461                     bool use_e_flag,
0462                     bool use_kt_flag,
0463                     Column_Vector& xm,
0464                     Column_Vector& ym,
0465                     Matrix& G_i,
0466                     Diagonal_Matrix& Y)
0467     //
0468     // Purpose: Take the information from a Fourvec_Event EV and pack
0469     //          it into vectors of well- and poorly-measured variables.
0470     //          Also set up the error matrices.
0471     //
0472     // Inputs:
0473     //   ev -          The event to pack.
0474     //   use_e_flag -  If true, we're using E as the fit variable, otherwise p.
0475     //   use_kt_flag - True if we're to pack kt variables.
0476     //
0477     // Outputs:
0478     //   xm -          Vector of well-measured variables.
0479     //   ym -          Vector of poorly-measured variables.
0480     //   G_i -         Error matrix for well-measured variables.
0481     //   Y -           Inverse error matrix for poorly-measured variables.
0482     //
0483     {
0484       // Number of objects in the event.
0485       int nobjs = ev.nobjs();
0486 
0487       int n_measured_vars = nobjs * 3;
0488       if (use_kt_flag)
0489         n_measured_vars += 2;
0490 
0491       // Clear the error matrix.
0492       G_i = Matrix(n_measured_vars, n_measured_vars, 0);
0493 
0494       // Loop over objects.
0495       for (int i = 0; i < nobjs; i++) {
0496         const FE_Obj& obj = ev.obj(i);
0497         int this_index = obj_index(i);
0498         set_p_eta_phi_vec(obj, xm, this_index, use_e_flag);
0499         G_i.sub(this_index, this_index, error_matrix(obj.p_error, obj.phi_error, obj.eta_error));
0500       }
0501 
0502       if (use_kt_flag) {
0503         // Set up kt.
0504         int kt_ndx = obj_index(nobjs);
0505         xm(kt_ndx + x_offs) = ev.kt().x();
0506         xm(kt_ndx + y_offs) = ev.kt().y();
0507 
0508         // And its error matrix.
0509         G_i(kt_ndx + x_offs, kt_ndx + x_offs) = ev.kt_x_error() * ev.kt_x_error();
0510         G_i(kt_ndx + y_offs, kt_ndx + y_offs) = ev.kt_y_error() * ev.kt_y_error();
0511         G_i(kt_ndx + x_offs, kt_ndx + y_offs) = ev.kt_xy_covar();
0512         G_i(kt_ndx + y_offs, kt_ndx + x_offs) = ev.kt_xy_covar();
0513       }
0514 
0515       // Handle a neutrino.
0516       if (ev.has_neutrino()) {
0517         ym(nu_z) = ev.nu().z();
0518         Y = Diagonal_Matrix(1, 0);
0519       }
0520     }
0521 
0522     /**
0523     @brief Helper function: Update the content of <i>ev</i> from the variable
0524     sets <i>x</i> and <i>y</i>.
0525 
0526     @param ev The event to update.
0527 
0528     @param x Vector of well-measured variables.
0529 
0530     @param y Vector of poorly-measured variables.
0531 
0532     @param use_e_flag If TRUE, then we are using energy <i>E</i> as the
0533     fit variable.<br>
0534     If FALSE, we are using magnitude of three-momentum <i>p</i> as the
0535     fit variable.
0536 
0537     @param use_kt_flag If TRUE, then we are also pack \f$k){T}\f$
0538     variables.<br>
0539     If FALSE, then \f$k){T}\f$ variables are not pack.
0540  */
0541     void unpack_event(
0542         Fourvec_Event& ev, const Column_Vector& x, const Column_Vector& y, bool use_e_flag, bool use_kt_flag)
0543     //
0544     // Purpose: Update the contents of EV from the variable sets X and Y.
0545     //
0546     // Inputs:
0547     //   ev -          The event.
0548     //   x -           Vector of well-measured variables.
0549     //   use_e_flag -  If true, we're using E as the fit variable, otherwise p.
0550     //   use_kt_flag - True if we're too unpack kt variables.
0551     //
0552     // Outputs:
0553     //   ev -          The event after updating.
0554     //
0555     {
0556       // Do all the objects.
0557       Fourvec sum = Fourvec(0);
0558       for (int j = 0; j < ev.nobjs(); j++) {
0559         const FE_Obj& obj = ev.obj(j);
0560         ev.set_obj_p(j, get_p_eta_phi_vec(x, obj_index(j), obj, use_e_flag));
0561         sum += obj.p;
0562       }
0563 
0564       if (use_kt_flag) {
0565         int kt_ndx = obj_index(ev.nobjs());
0566         Fourvec kt = Fourvec(x(kt_ndx + x_offs), x(kt_ndx + y_offs), 0, 0);
0567         Fourvec nu = kt - sum;
0568         if (ev.has_neutrino()) {
0569           nu.setPz(y(nu_z));
0570           adjust_e_for_mass(nu, 0);
0571           ev.set_nu_p(nu);
0572         } else {
0573           adjust_e_for_mass(nu, 0);
0574           ev.set_x_p(nu);
0575         }
0576       }
0577     }
0578 
0579   }  // unnamed namespace
0580 
0581   //*************************************************************************
0582   // Constraint evaluation.
0583   //
0584 
0585   namespace {
0586 
0587     /**
0588     @brief Helper function: Compute the dot product
0589     \f$\mathbf{v}_{1} \cdot \mathbf{v}_{2}\f$ and its gradient w.r.t.
0590     \f$p\f$, \f$\phi\f$, and \f$\theta\f$ of each four-vector.
0591 
0592     @param v1 The first four-vector in the dot product.
0593 
0594     @param v2 The second four-vector in the dot product.
0595 
0596     @param use_e_flag If TRUE, then we are using energy <i>E</i> as the
0597     fit variable.<br>
0598     If FALSE, we are using magnitude of three-momentum <i>p</i> as the
0599     fit variable.
0600 
0601     @param v1_x Gradients of the dot product w.r.t. to
0602     \f$p\f$, \f$\phi\f$, and \f$\theta\f$ of \f$\mathbf{v}_{1}\f$.
0603 
0604     @param v2_x Gradients of the dot product w.r.t. to
0605     \f$p\f$, \f$\phi\f$, and \f$\theta\f$ of \f$\mathbf{v}_{2}\f$.
0606 
0607     @param badflag Return status, set to TRUE for singular cases
0608     when calculating the dot product or gradient (zero values for
0609     any of these: energy \f$E\f$, magnitude of three-momentum \f$p\f$,
0610     magnitude of transverse component of three-momentum \f$p_{T}\f$.
0611 
0612     @par Returns:
0613     The dot product.
0614  */
0615     double dot_and_gradient(
0616         const Fourvec& v1, const Fourvec& v2, bool use_e_flag, double v1_x[3], double v2_x[3], bool& badflag)
0617     //
0618     // Purpose: Compute the dot product v1.v2 and its gradients wrt
0619     //          p, phi, and theta of each 4-vector.
0620     //
0621     // Inputs:
0622     //   v1 -          The first 4-vector in the dot product.
0623     //   v2 -          The second 4-vector in the dot product.
0624     //   use_e_flag -  If true, we're using E as the fit variable, otherwise p.
0625     //
0626     // Outputs:
0627     //   v1_x -        Gradients of the dot product wrt v1's p, phi, theta.
0628     //   v2_x -        Gradients of the dot product wrt v2's p, phi, theta.
0629     //   badflag -     Set to true for the singular case (vectors vanish).
0630     //
0631     // Returns:
0632     //   The dot product.
0633     //
0634     {
0635       // Calculate the dot product.
0636       double dot = v1 * v2;
0637 
0638       double p1 = v1.vect().mag();
0639       double p2 = v2.vect().mag();
0640       double e1 = v1.e();
0641       double e2 = v2.e();
0642       double pt1 = v1.vect().perp();
0643       double pt2 = v2.vect().perp();
0644 
0645       // Protect against the singular case.
0646       badflag = false;
0647       if (p1 == 0 || p2 == 0 || e1 == 0 || e2 == 0 || pt1 == 0 || pt2 == 0) {
0648         badflag = true;
0649         v1_x[p_offs] = v1_x[phi_offs] = v1_x[eta_offs] = 0;
0650         v2_x[p_offs] = v2_x[phi_offs] = v2_x[eta_offs] = 0;
0651         return false;
0652       }
0653 
0654       // Calculate the gradients.
0655       v1_x[p_offs] = (dot - v1.m2() * e2 / e1) / p1;
0656       v2_x[p_offs] = (dot - v2.m2() * e1 / e2) / p2;
0657 
0658       if (use_e_flag) {
0659         v1_x[p_offs] *= e1 / p1;
0660         v2_x[p_offs] *= e2 / p2;
0661       }
0662 
0663       v1_x[phi_offs] = v1(1) * v2(0) - v1(0) * v2(1);
0664       v2_x[phi_offs] = -v1_x[phi_offs];
0665 
0666       double fac = v1(0) * v2(0) + v1(1) * v2(1);
0667       v1_x[eta_offs] = pt1 * v2(2) - v1(2) / pt1 * fac;
0668       v2_x[eta_offs] = pt2 * v1(2) - v2(2) / pt2 * fac;
0669 
0670       return dot;
0671     }
0672 
0673     /**
0674     @brief Helper function: Tally up dot product gradients for an object into
0675     <i>Bx</i>, the gradients of the well-measured variables.
0676 
0677     @param constraint_no The number/index of the constraint.
0678 
0679     @param base_index The index in the vector of well-measured variables
0680     of the first variable for this object.
0681 
0682     @param sign The sign with which these gradients should be added into
0683     <i>Bx</i>, either \f$+1\f$ or \f$-1\f$ (that is, which side of
0684     the constraint equation).
0685 
0686     @param grad The gradient for this object w.r.t. to \f$p\f$,
0687     \f$\phi\f$, and \f$\theta\f$.
0688 
0689     @param Bx The gradients of well-measured variables.
0690 
0691     @par Output:
0692     - <i>Bx</i>.
0693  */
0694     void addin_obj_gradient(int constraint_no, int sign, int base_index, const double grad[], Matrix& Bx)
0695     //
0696     // Purpose: Tally up the dot product gradients for an object
0697     //          into Bx.
0698     //
0699     // Inputs:
0700     //   constraint_no-The number of the constraint.
0701     //   base_index -  The index in the well-measured variable list
0702     //                 of the first variable for this object.
0703     //   sign -        The sign with which these gradients should be
0704     //                 added into Bx, either +1 or -1.  (I.e., which
0705     //                 side of the constraint equation.)
0706     //   grad -        The gradients for this object, vs. p, phi, theta.
0707     //   Bx -          The well-measured variable gradients.
0708     //
0709     // Outputs:
0710     //   Bx -          The well-measured variable gradients, updated.
0711     //
0712     {
0713       Bx(base_index + p_offs, constraint_no) += sign * grad[p_offs];
0714       Bx(base_index + phi_offs, constraint_no) += sign * grad[phi_offs];
0715       Bx(base_index + eta_offs, constraint_no) += sign * grad[eta_offs];
0716     }
0717 
0718     /**
0719 
0720     @brief Helper function: Tally up the dot product gradients for a neutrino
0721     into <i>Bx</i> and <i>By</i>, the gradient vector of well-measured and
0722     poorly-measured variables, respectively.
0723 
0724     @param constraint_no The number/index of the constraint.
0725 
0726     @param sign The sign with which these gradients should be added into
0727     <i>Bx</i>, either \f$+1\f$ or \f$-1\f$ (that is, which side of
0728     the constraint equation).
0729 
0730     @param kt_ndx The index of the \f$k_{T}\f$ variables in the variables
0731     array.
0732 
0733     @param grad The gradient for this object w.r.t. to \f$p\f$,
0734     \f$\phi\f$, and \f$\theta\f$.
0735 
0736     @param Bx The gradients of well-measured variables.
0737 
0738     @param By The gradients of poorly-measured variables.
0739 
0740     @par Output:
0741     - <i>Bx</i> The updated gradients of the well-measured variables.
0742     - <i>By</i> The updated gradients of the poorly-measured variables.
0743 
0744  */
0745     void addin_nu_gradient(int constraint_no, int sign, int kt_ndx, const double grad[], Matrix& Bx, Matrix& By)
0746     //
0747     // Purpose: Tally up the dot product gradients for a neutrino
0748     //          into Bx and By.
0749     //
0750     // Inputs:
0751     //   constraint_no-The number of the constraint.
0752     //   sign -        The sign with which these gradients should be
0753     //                 added into Bx, either +1 or -1.  (I.e., which
0754     //                 side of the constraint equation.)
0755     //   kt_ndx -      The index of the kt variables in the variables array.
0756     //   grad -        The gradients for this object, vs. p, phi, theta.
0757     //   Bx -          The well-measured variable gradients.
0758     //   By -          The poorly-measured variable gradients.
0759     //
0760     // Outputs:
0761     //   Bx -          The well-measured variable gradients, updated.
0762     //   By -          The poorly-measured variable gradients, updated.
0763     //
0764     {
0765       Bx(kt_ndx + x_offs, constraint_no) += sign * grad[p_offs];    // Really p for now.
0766       Bx(kt_ndx + y_offs, constraint_no) += sign * grad[phi_offs];  // Really phi for now.
0767       By(nu_z, constraint_no) += sign * grad[eta_offs];             // Really theta ...
0768     }
0769 
0770     /**
0771     @brief Helper function: Tally up the dot product gradients for an object
0772     (which may or may not be a neutrino) into <i>Bx</i> and <i>By</i>.
0773 
0774     @param ev The event we are fitting.
0775 
0776     @param constraint_no The number/index of the constraint.
0777 
0778     @param sign The sign with which these gradients should be added into
0779     <i>Bx</i>, either \f$+1\f$ or \f$-1\f$ (that is, which side of
0780     the constraint equation).
0781 
0782     @param obj_no The number of the object.
0783 
0784     @param grad The gradient for this object w.r.t. to \f$p\f$,
0785     \f$\phi\f$, and \f$\theta\f$.
0786 
0787     @param Bx The gradients of the well-measured variables.
0788 
0789     @param By The gradients of the poorly-measured variables.
0790 
0791     @par Output:
0792     - <i>Bx</i> The updated gradients of the well-measured variables.<br>
0793     - <i>By</i> The updated gradients of the poorly-measured variables.
0794 
0795  */
0796     void addin_gradient(
0797         const Fourvec_Event& ev, int constraint_no, int sign, int obj_no, const double grad[], Matrix& Bx, Matrix& By)
0798     //
0799     // Purpose: Tally up the dot product gradients for an object (which may
0800     //          or may not be a neutrino) into Bx and By.
0801     //
0802     // Inputs:
0803     //   ev -          The event we're fitting.
0804     //   constraint_no-The number of the constraint.
0805     //   sign -        The sign with which these gradients should be
0806     //                 added into Bx, either +1 or -1.  (I.e., which
0807     //                 side of the constraint equation.)
0808     //   obj_no -      The number of the object.
0809     //   grad -        The gradients for this object, vs. p, phi, theta.
0810     //   Bx -          The well-measured variable gradients.
0811     //   By -          The poorly-measured variable gradients.
0812     //
0813     // Outputs:
0814     //   Bx -          The well-measured variable gradients, updated.
0815     //   By -          The poorly-measured variable gradients, updated.
0816     //
0817     {
0818       if (obj_no >= ev.nobjs()) {
0819         assert(obj_no == ev.nobjs());
0820         addin_nu_gradient(constraint_no, sign, obj_index(obj_no), grad, Bx, By);
0821       } else
0822         addin_obj_gradient(constraint_no, sign, obj_index(obj_no), grad, Bx);
0823     }
0824 
0825     /**
0826     @brief Helper function: Tally up the gradients from a single dot product
0827     into <i>Bx</i> and <i>By</i>.
0828 
0829     @param ev The event we are fitting.
0830 
0831     @param constraint_no The number/index of the constraint.
0832 
0833     @param sign The sign with which these gradients should be added into
0834     <i>Bx</i>, either \f$+1\f$ or \f$-1\f$ (that is, which side of
0835     the constraint equation).
0836 
0837     @param i The number/index of the first object.
0838 
0839     @param igrad The gradients of the first object w.r.t. to \f$p\f$,
0840     \f$\phi\f$, and \f$\theta\f$.
0841 
0842     @param j The number/index of the second object.
0843 
0844     @param jgrad The gradients of the second object w.r.t. to \f$p\f$,
0845     \f$\phi\f$, and \f$\theta\f$.
0846 
0847     @param Bx The gradients of the well-measured variables.
0848 
0849     @param By The gradients of the poorly-measured variables.
0850 
0851     @par Output:
0852     - <i>Bx</i> The updated gradients of the well-measured variables.<br>
0853     - <i>By</i> The updated gradients of the poorly-measured variables.
0854  */
0855     void addin_gradients(const Fourvec_Event& ev,
0856                          int constraint_no,
0857                          int sign,
0858                          int i,
0859                          const double igrad[],
0860                          int j,
0861                          const double jgrad[],
0862                          Matrix& Bx,
0863                          Matrix& By)
0864     //
0865     // Purpose: Tally up the gradients from a single dot product into Bx and By.
0866     //
0867     // Inputs:
0868     //   ev -          The event we're fitting.
0869     //   constraint_no-The number of the constraint.
0870     //   sign -        The sign with which these gradients should be
0871     //                 added into Bx, either +1 or -1.  (I.e., which
0872     //                 side of the constraint equation.)
0873     //   i -           The number of the first object.
0874     //   igrad -       The gradients for the first object, vs. p, phi, theta.
0875     //   j -           The number of the second object.
0876     //   jgrad -       The gradients for the second object, vs. p, phi, theta.
0877     //   Bx -          The well-measured variable gradients.
0878     //   By -          The poorly-measured variable gradients.
0879     //
0880     // Outputs:
0881     //   Bx -          The well-measured variable gradients, updated.
0882     //   By -          The poorly-measured variable gradients, updated.
0883     //
0884     {
0885       addin_gradient(ev, constraint_no, sign, i, igrad, Bx, By);
0886       addin_gradient(ev, constraint_no, sign, j, jgrad, Bx, By);
0887     }
0888 
0889     /**
0890     @brief Helper function: Add the \f$m^{2}\f$ into the constraint values.
0891 
0892     @param ev The event we are fitting.
0893 
0894     @param constraints The list of constraints.
0895 
0896     @param F Vector of constraint values.
0897 
0898     @par Output:
0899     - <i>F</i> The updated vector of constraint values.
0900  */
0901     void add_mass_terms(const Fourvec_Event& ev, const vector<Constraint>& constraints, Row_Vector& F)
0902     //
0903     // Purpose: Add the m^2 terms into the constraint values.
0904     //
0905     // Inputs:
0906     //   ev -          The event we're fitting.
0907     //   constraints - The list of constraints.
0908     //   F -           Vector of constraint values.
0909     //
0910     // Outputs:
0911     //   F -           Vector of constraint values, updated.
0912     //
0913     {
0914       for (std::vector<Constraint>::size_type i = 0; i < constraints.size(); i++)
0915         F(i + 1) += constraints[i].sum_mass_terms(ev);
0916     }
0917 
0918     /**
0919     @brief Helper function: Calculate the mass constraints and gradients.
0920     At this sage, the gradients are calculated not quite with respect to
0921     the fit variables, instead for all objects (including neutrino) we
0922     calculate the gradients with respect to \f$p\f$, \f$\phi\f$, and
0923     \f$\theta\f$.  They will be converted via the appropriate
0924     Jacobian transformation later
0925 
0926     @param ev The event we are fitting.
0927 
0928     @param pt Table of cached pair assigments.
0929 
0930     @param constraints The list of constraints.
0931 
0932     @param use_e_flag If TRUE, then we are using energy <i>E</i> as the
0933     fit variable.<br>
0934     If FALSE, we are using magnitude of three-momentum <i>p</i> as the
0935     fit variable.
0936 
0937     @param F The vector of constraint values, should be passed with the
0938     correct dimensionaly.
0939 
0940     @param Bx The gradients of the well-measured variables, should be passed
0941     with the correct dimensionaly.
0942 
0943     @param By The gradients of the poorly-measured variables, should be
0944     passed with the scorrect dimensionaly.
0945 
0946     @par Output:
0947     - <i>F</i> The vector of constraint variables.<br>
0948     - <i>Bx</i> The updated gradients of the well-measured variables.<br>
0949     - <i>By</i> The updated gradients of the poorly-measured variables.
0950  */
0951     bool calculate_mass_constraints(const Fourvec_Event& ev,
0952                                     const Pair_Table& pt,
0953                                     const vector<Constraint>& constraints,
0954                                     bool use_e_flag,
0955                                     Row_Vector& F,
0956                                     Matrix& Bx,
0957                                     Matrix& By)
0958     //
0959     // Purpose: Calculate the mass constraints and gradients.
0960     //          Note: At this stage, the gradients are calculated not
0961     //                quite with respect to the fit variables; instead, for
0962     //                all objects (including the neutrino) we calculate
0963     //                the gradients with respect to p, phi, theta.  They'll
0964     //                be converted via appropriate Jacobian transformations
0965     //                later.
0966     //
0967     // Inputs:
0968     //   ev -          The event we're fitting.
0969     //   pt -          Table of cached pair assignments.
0970     //   constraints - The list of constraints.
0971     //   use_e_flag -  If true, we're using E as the fit variable, otherwise p.
0972     //
0973     // Outputs:
0974     //   (nb. these should be passed in correctly dimensioned.)
0975     //   F -           Vector of constraint values.
0976     //   Bx -          The well-measured variable gradients.
0977     //   By -          The poorly-measured variable gradients.
0978     //
0979     {
0980       int npairs = pt.npairs();
0981       for (int p = 0; p < npairs; p++) {
0982         const Objpair& objpair = pt.get_pair(p);
0983         int i = objpair.i();
0984         int j = objpair.j();
0985         double igrad[3], jgrad[3];
0986         bool badflag = false;
0987         double dot = dot_and_gradient(ev.obj(i).p, ev.obj(j).p, use_e_flag, igrad, jgrad, badflag);
0988         if (badflag)
0989           return false;
0990 
0991         for (std::vector<Constraint>::size_type k = 0; k < constraints.size(); k++)
0992           if (objpair.for_constraint(k)) {
0993             F(k + 1) += objpair.for_constraint(k) * dot;
0994             addin_gradients(ev, k + 1, objpair.for_constraint(k), i, igrad, j, jgrad, Bx, By);
0995           }
0996       }
0997 
0998       add_mass_terms(ev, constraints, F);
0999       return true;
1000     }
1001 
1002     /**
1003     @brief Helper function: Perform the Jacobian tranformation for
1004     \f$(p_{\nu}^{x},p_{\nu}^{y}) \to (k_{T}^{x},k_{T}^{y})\f$ for
1005     a given object.
1006 
1007     @param ndx Index of the object for which to transform gradients.
1008 
1009     @param v The object's four-momentum.
1010 
1011     @param Bx The gradients of the well-measured variables.
1012 
1013     @param use_e_flag If TRUE, then we are using energy <i>E</i> as the
1014     fit variable.<br>
1015     If FALSE, we are using magnitude of three-momentum <i>p</i> as the
1016     fit variable.
1017 
1018     @param kt_ndx The index of the \f$k_{T}\f$ variables in the variables
1019     array.
1020 
1021     @par Output:
1022     - <i>Bx</i> The updated gradients of the well-measured variables.
1023  */
1024     void add_nuterm(unsigned ndx, const Fourvec& v, Matrix& Bx, bool use_e_flag, int kt_ndx)
1025     //
1026     // Purpose: Carry out the Jacobian transformation for
1027     //          (p_nu^x,p_nu_y) -> (kt^x, kt_y) for a given object.
1028     //
1029     // Inputs:
1030     //   ndx -         Index of the object for which to transform gradients.
1031     //   v -           The object's 4-momentum.
1032     //   Bx -          The well-measured variable gradients.
1033     //   use_e_flag -  If true, we're using E as the fit variable, otherwise p.
1034     //   kt_ndx -      The index of the kt variables in the variables array.
1035     //
1036     // Outputs:
1037     //   Bx -          The well-measured variable gradients, updated.
1038     //
1039     {
1040       double px = v.px();
1041       double py = v.py();
1042       double cot_theta = v.pz() / v.vect().perp();
1043 
1044       for (int j = 1; j <= Bx.num_col(); j++) {
1045         double dxnu = Bx(kt_ndx + x_offs, j);
1046         double dynu = Bx(kt_ndx + y_offs, j);
1047 
1048         if (dxnu != 0 || dynu != 0) {
1049           double fac = 1 / v.vect().mag();
1050           if (use_e_flag)
1051             fac = v.e() * fac * fac;
1052           Bx(ndx + p_offs, j) -= (px * dxnu + py * dynu) * fac;
1053           Bx(ndx + phi_offs, j) += py * dxnu - px * dynu;
1054           Bx(ndx + eta_offs, j) -= (px * dxnu + py * dynu) * cot_theta;
1055         }
1056       }
1057     }
1058 
1059     /**
1060     @brief Helper function: Carry out the Jacobian transformations for
1061     the neutrino.  First, convert from spherical coordinates
1062     \f$(p,\phi,\eta)\f$ to rectangular \f$(x,y,z)\f$.  Then convert from
1063     neutrino \f$p_{T}\f$ components to \f$k_{T}\f$ components.
1064 
1065     @param ev The event we are fitting.
1066 
1067     @param use_e_flag If TRUE, then we are using energy <i>E</i> as the
1068     fit variable.<br>
1069     If FALSE, we are using magnitude of three-momentum <i>p</i> as the
1070     fit variable.
1071 
1072     @param Bx The gradients of the well-measured variables, should be passed
1073     with the correct dimensionaly.
1074 
1075     @param By The gradients of the poorly-measured variables, should be
1076     passed with the scorrect dimensionaly.
1077 
1078     @par Output:
1079     - <i>Bx</i> The updated gradients of the well-measured variables.
1080     - <i>By</i> The updated gradients of the poorly-measured variables.
1081  */
1082     void convert_neutrino(const Fourvec_Event& ev, bool use_e_flag, Matrix& Bx, Matrix& By)
1083     //
1084     // Purpose: Carry out the Jacobian transformations the neutrino.
1085     //          First, convert from spherical (p, phi, theta) coordinates
1086     //          to rectangular (x, y, z).  Then convert from neutrino pt
1087     //          components to kt components.
1088     //
1089     // Inputs:
1090     //   ev -          The event we're fitting.
1091     //   use_e_flag -  If true, we're using E as the fit variable, otherwise p.
1092     //   Bx -          The well-measured variable gradients.
1093     //   By -          The poorly-measured variable gradients.
1094     //
1095     // Outputs:
1096     //   Bx -          The well-measured variable gradients, updated.
1097     //   By -          The poorly-measured variable gradients, updated.
1098     //
1099     {
1100       int nconstraints = Bx.num_col();
1101 
1102       const Fourvec& nu = ev.nu();
1103 
1104       // convert neutrino from polar coordinates to rectangular.
1105       double pnu2 = nu.vect().mag2();
1106       double pnu = sqrt(pnu2);
1107       double ptnu2 = nu.vect().perp2();
1108       double ptnu = sqrt(ptnu2);
1109 
1110       // Doesn't matter whether we use E or P here, since nu is massless.
1111 
1112       double thfac = nu.z() / pnu2 / ptnu;
1113       double fac[3][3];
1114       fac[0][0] = nu(0) / pnu;
1115       fac[0][1] = nu(1) / pnu;
1116       fac[0][2] = nu(2) / pnu;
1117       fac[1][0] = -nu(1) / ptnu2;
1118       fac[1][1] = nu(0) / ptnu2;
1119       fac[1][2] = 0;
1120       fac[2][0] = nu(0) * thfac;
1121       fac[2][1] = nu(1) * thfac;
1122       fac[2][2] = -ptnu / pnu2;
1123 
1124       int kt_ndx = obj_index(ev.nobjs());
1125       for (int j = 1; j <= nconstraints; j++) {
1126         double tmp1 = fac[0][0] * Bx(kt_ndx + x_offs, j) + fac[1][0] * Bx(kt_ndx + y_offs, j) + fac[2][0] * By(nu_z, j);
1127         Bx(kt_ndx + y_offs, j) =
1128             fac[0][1] * Bx(kt_ndx + x_offs, j) + fac[1][1] * Bx(kt_ndx + y_offs, j) + fac[2][1] * By(nu_z, j);
1129         By(nu_z, j) = fac[0][2] * Bx(kt_ndx + x_offs, j) + fac[2][2] * By(nu_z, j);
1130 
1131         Bx(kt_ndx + x_offs, j) = tmp1;
1132       }
1133 
1134       // Add nu terms.
1135       for (int j = 0; j < ev.nobjs(); j++) {
1136         add_nuterm(obj_index(j), ev.obj(j).p, Bx, use_e_flag, kt_ndx);
1137       }
1138     }
1139 
1140     /**
1141     @brief Helpere function: Calculate the overall \f$k_{T}\f$ constraints
1142     \f$(k_{T} = 0)\f$ for the case where there is no neutrino.
1143 
1144     @param ev The event we are fitting.
1145 
1146     @param use_e_flag If TRUE, then we are using energy <i>E</i> as the
1147     fit variable.<br>
1148     If FALSE, we are using magnitude of three-momentum <i>p</i> as the
1149     fit variable.
1150 
1151     @param F Vector of constraint values.
1152 
1153     @param Bx The gradients of the well-measured variables, should be passed
1154     with the correct dimensionaly.
1155 
1156     @par Output:
1157     - <i>F</i> The updated vectgor of constraints values.<br>
1158     - <i>Bx</i> The updated gradients of the well-measured variables.
1159  */
1160     void calculate_kt_constraints(const Fourvec_Event& ev, bool use_e_flag, Row_Vector& F, Matrix& Bx)
1161     //
1162     // Purpose: Calculate the overall kt constraints (kt = 0) for the case
1163     //          where there is no neutrino.
1164     //
1165     // Inputs:
1166     //   ev -          The event we're fitting.
1167     //   use_e_flag -  If true, we're using E as the fit variable, otherwise p.
1168     //   F -           Vector of constraint values.
1169     //   Bx -          The well-measured variable gradients.
1170     //
1171     // Outputs:
1172     //   F -           Vector of constraint values, updated.
1173     //   Bx -          The well-measured variable gradients, updated.
1174     //
1175     {
1176       Fourvec tmp = Fourvec(0);
1177       int base = F.num_col() - 2;
1178       int nobjs = ev.nobjs();
1179       for (int j = 0; j < nobjs; j++) {
1180         const Fourvec& obj = ev.obj(j).p;
1181         tmp += obj;
1182 
1183         int ndx = obj_index(j);
1184         double p = obj.vect().mag();
1185         double cot_theta = obj.z() / obj.vect().perp();
1186 
1187         Bx(ndx + p_offs, base + 1) = obj(0) / p;
1188         Bx(ndx + phi_offs, base + 1) = -obj(1);
1189         Bx(ndx + eta_offs, base + 1) = obj(0) * cot_theta;
1190 
1191         Bx(ndx + p_offs, base + 2) = obj(1) / p;
1192         Bx(ndx + phi_offs, base + 2) = obj(0);
1193         Bx(ndx + eta_offs, base + 2) = obj(1) * cot_theta;
1194 
1195         if (use_e_flag) {
1196           Bx(ndx + p_offs, base + 1) *= obj.e() / p;
1197           Bx(ndx + p_offs, base + 2) *= obj.e() / p;
1198         }
1199       }
1200 
1201       int kt_ndx = obj_index(nobjs);
1202       Bx(kt_ndx + x_offs, base + 1) = -1;
1203       Bx(kt_ndx + y_offs, base + 2) = -1;
1204 
1205       F(base + 1) = tmp(0) - ev.kt().x();
1206       F(base + 2) = tmp(1) - ev.kt().y();
1207     }
1208 
1209     /**
1210     @brief Do the Jacobian transformation from \f$\theta\f$ to \f$\eta\f$
1211     for a single object.
1212 
1213     @param i The index of the object.
1214 
1215     @param cos_theta The \f$\cos \theta\f$ for the object.
1216 
1217     @param Bx The gradients of the well-measured variables.
1218 
1219     @par Output:
1220     - <i>Bx</i> The updated gradients of the well-measured variables.
1221  */
1222     void ddtheta_to_ddeta(int i, double cos_theta, Matrix& Bx)
1223     //
1224     // Purpose: Do the Jacobian transformation from theta -> eta
1225     //          for a single object.
1226     //
1227     // Inputs:
1228     //   i -           The index of the object.
1229     //   cos_theta -   cos(theta) for the object.
1230     //   Bx -          The well-measured variable gradients.
1231     //
1232     // Outputs:
1233     //   Bx -          The well-measured variable gradients, updated.
1234     //
1235     {
1236       double sin_theta = sqrt(1 - cos_theta * cos_theta);
1237       for (int j = 1; j <= Bx.num_col(); j++)
1238         Bx(i, j) *= -sin_theta; /* \sin\theta = 1 / \cosh\eta */
1239     }
1240 
1241   }  // unnamed namespace
1242 
1243   /**
1244     @brief Concrete realization of the Constraint_Calculator class.
1245     Evaluate constraints at the point described by <i>x</i> and
1246     <i>y</i>
1247     (well-measured and poorly-measured variables, respectively).  The results
1248     should be stored in <i>F</i>.  <i>Bx</i> and <i>By</i> should be set to
1249     the gradients of <i>F</i> with respect to <i>x</i> and <i>y</i>,
1250     respectively.
1251 
1252     @param x  Column_Vector of well-measured variables.
1253 
1254     @param y  Column_Vector of poorly-measured variables.
1255 
1256     @param F  Row_Vector contains the results of the constraint evaluation.
1257 
1258     @param Bx Gradients of <i>F</i> with respect to <i>x</i>
1259 
1260     @param By Gradients of <i>F</i> with respect to <i>y</i>
1261 
1262     @par Return:
1263     <b>true</b> if the point <i>(x,y)</i> is accepted.<br>
1264     <b>false</b> if the point <i>(x,y)</i> is rejected
1265     (i.e., in an unphysical region).  The constraints need not be
1266     evaluated in that case.
1267  */
1268   class Fourvec_Constraint_Calculator : public Constraint_Calculator
1269   //
1270   // Purpose: Constraint evaluator.
1271   //
1272   {
1273   public:
1274     // Constructor, destructor.
1275     Fourvec_Constraint_Calculator(Fourvec_Event& ev,
1276                                   const vector<Constraint>& constraints,
1277                                   const Fourvec_Constrainer_Args& args);
1278     ~Fourvec_Constraint_Calculator() override {}
1279 
1280     // Evaluate constraints at the point described by X and Y (well-measured
1281     // and poorly-measured variables, respectively).  The results should
1282     // be stored in F.  BX and BY should be set to the gradients of F with
1283     // respect to X and Y, respectively.
1284     //
1285     // Return true if the point X, Y is accepted.
1286     // Return false if it is rejected (i.e., in an unphysical region).
1287     // The constraints need not be evaluated in that case.
1288     bool eval(const Column_Vector& x, const Column_Vector& y, Row_Vector& F, Matrix& Bx, Matrix& By) override;
1289 
1290     // Calculate the constraint functions and gradients.
1291     bool calculate_constraints(Row_Vector& F, Matrix& Bx, Matrix& By) const;
1292 
1293   private:
1294     // The event we're fitting.
1295     Fourvec_Event& _ev;
1296 
1297     // Vector of constraints.
1298     const vector<Constraint>& _constraints;
1299 
1300     // Argument values.
1301     const Fourvec_Constrainer_Args& _args;
1302 
1303     // The pair table.
1304     Pair_Table _pt;
1305   };
1306 
1307   /**
1308     @brief Constructor
1309 
1310     @param ev The event we are fitting.
1311 
1312     @param constraints The list of constraints.
1313 
1314     @param args The parameter settings for this instance.
1315  */
1316   Fourvec_Constraint_Calculator::Fourvec_Constraint_Calculator(Fourvec_Event& ev,
1317                                                                const vector<Constraint>& constraints,
1318                                                                const Fourvec_Constrainer_Args& args)
1319       //
1320       // Purpose: Constructor.
1321       //
1322       // Inputs:
1323       //   ev -          The event we're fitting.
1324       //   constraints - The list of constraints.
1325       //   args -        The parameter settings.
1326       //
1327       //
1328       : Constraint_Calculator(constraints.size() + ((ev.has_neutrino() || args.ignore_met()) ? 0 : 2)),
1329         _ev(ev),
1330         _constraints(constraints),
1331         _args(args),
1332         _pt(constraints, ev)
1333 
1334   {}
1335 
1336   /**
1337     @brief Calculate the constraint functions and gradients.
1338 
1339     @param F Vector of constraint values.
1340 
1341     @param Bx The gradient of well-measured variables.
1342 
1343     @param By The gradient of poorly-measured variables.
1344 
1345     @par Output:
1346     - <i>F</i>.
1347     - <i>Bx</i>.
1348     - <i>By</i>.
1349  */
1350   bool Fourvec_Constraint_Calculator::calculate_constraints(Row_Vector& F, Matrix& Bx, Matrix& By) const
1351   //
1352   // Purpose: Calculate the constraint functions and gradients.
1353   //
1354   // Outputs:
1355   //   F -           Vector of constraint values.
1356   //   Bx -          The well-measured variable gradients.
1357   //   By -          The poorly-measured variable gradients.
1358   //
1359   {
1360     // Clear the matrices.
1361     Bx = Matrix(Bx.num_row(), Bx.num_col(), 0);
1362     By = Matrix(By.num_row(), By.num_col(), 0);
1363     F = Row_Vector(F.num_col(), 0);
1364 
1365     const double p_eps = 1e-10;
1366 
1367     if (_ev.has_neutrino() && _ev.nu().z() > _args.e_com()) {
1368       return false;
1369     }
1370 
1371     int nobjs = _ev.nobjs();
1372 
1373     // Reject the point if any of the momenta get too small.
1374     for (int j = 0; j < nobjs; j++) {
1375       if (_ev.obj(j).p.perp() <= p_eps || _ev.obj(j).p.e() <= p_eps) {
1376         return false;
1377       }
1378     }
1379 
1380     if (!calculate_mass_constraints(_ev, _pt, _constraints, _args.use_e(), F, Bx, By))
1381       return false;
1382 
1383     if (_ev.has_neutrino())
1384       convert_neutrino(_ev, _args.use_e(), Bx, By);
1385     else if (!_args.ignore_met()) {
1386       /* kt constraints */
1387       calculate_kt_constraints(_ev, _args.use_e(), F, Bx);
1388     }
1389 
1390     /* convert d/dtheta to d/deta */
1391     for (int j = 0; j < nobjs; j++) {
1392       ddtheta_to_ddeta(obj_index(j) + eta_offs, _ev.obj(j).p.cosTheta(), Bx);
1393     }
1394 
1395     /* handle muons */
1396     for (int j = 0; j < nobjs; j++) {
1397       const FE_Obj& obj = _ev.obj(j);
1398       if (obj.muon_p) {
1399         // Again, E vs. P doesn't matter here since we assume muons to be massless.
1400         double pmu2 = obj.p.vect().mag2();
1401         int ndx = obj_index(j) + p_offs;
1402         for (int k = 1; k <= Bx.num_col(); k++)
1403           Bx(ndx, k) = -Bx(ndx, k) * pmu2;
1404       }
1405     }
1406 
1407     return true;
1408   }
1409 
1410   /**
1411     @brief Evaluate constraints at the point described by <i>x</i> and
1412     <i>y</i>
1413     (well-measured and poorly-measured variables, respectively).  The results
1414     should be stored in <i>F</i>.  <i>Bx</i> and <i>By</i> should be set to
1415     the gradients of <i>F</i> with respect to <i>x</i> and <i>y</i>,
1416     respectively.
1417 
1418     @param x  Column_Vector of well-measured variables.
1419 
1420     @param y  Column_Vector of poorly-measured variables.
1421 
1422     @param F  Row_Vector contains the results of the constraint evaluation.
1423 
1424     @param Bx Gradients of <i>F</i> with respect to <i>x</i>
1425 
1426     @param By Gradients of <i>F</i> with respect to <i>y</i>
1427 
1428     @par Output:
1429     - <i>F</i>.
1430     - <i>Bx</i>.
1431     - <i>By</i>.
1432 
1433     @par Return:
1434     <b>true</b> if the point <i>(x,y)</i> is accepted.<br>
1435     <b>false</b> if the point <i>(x,y)</i> is rejected
1436     (i.e., in an unphysical region).  The constraints need not be
1437     evaluated in that case.
1438  */
1439   bool Fourvec_Constraint_Calculator::eval(
1440       const Column_Vector& x, const Column_Vector& y, Row_Vector& F, Matrix& Bx, Matrix& By)
1441   //
1442   // Purpose: Evaluate constraints at the point described by X and Y
1443   //          (well-measured and poorly-measured variables, respectively).
1444   //          The results should be stored in F.  BX and BY should be set
1445   //          to the gradients of F with respect to X and Y, respectively.
1446   //
1447   // Inputs:
1448   //   x -           Vector of well-measured variables.
1449   //   y -           Vector of poorly-measured variables.
1450   //
1451   // Outputs:
1452   //   F -           Vector of constraint values.
1453   //   Bx -          The well-measured variable gradients.
1454   //   By -          The poorly-measured variable gradients.
1455   //
1456   {
1457     int nobjs = _ev.nobjs();
1458 
1459     const double p_eps = 1e-10;
1460     const double eta_max = 10;
1461 
1462     // Give up if we've gone into an obviously unphysical region.
1463     for (int j = 0; j < nobjs; j++)
1464       if (x(obj_index(j) + p_offs) < p_eps || fabs(x(obj_index(j) + eta_offs)) > eta_max) {
1465         return false;
1466       }
1467 
1468     unpack_event(_ev, x, y, _args.use_e(), _ev.has_neutrino() || !_args.ignore_met());
1469 
1470     return calculate_constraints(F, Bx, By);
1471   }
1472 
1473   //*************************************************************************
1474   // Mass calculation.
1475   //
1476 
1477   namespace {
1478 
1479     /**
1480     @brief Helper function: Calculate the error propagation to find the
1481     uncertainty in the final mass.
1482 
1483     @param Q The final error matrix for the well-measured variables.
1484 
1485     @param R The final error matrix for the poorly-measured variables.
1486 
1487     @param S The final cross error matrix for the two sets of variables.
1488 
1489     @param Bx The gradient of the well-measured variables.
1490 
1491     @param By The gradient of the poorly-measured variables.
1492 
1493     @par Return:
1494     The uncertainty in the final mass.
1495  */
1496     double calculate_sigm(
1497         const Matrix& Q, const Matrix& R, const Matrix& S, const Column_Vector& Bx, const Column_Vector& By)
1498     //
1499     // Purpose: Do error propagation to find the uncertainty in the final mass.
1500     //
1501     // Inputs:
1502     //   Q -           The final error matrix for the well-measured variables.
1503     //   R -           The final error matrix for the poorly-measured variables.
1504     //   S -           The final cross error matrix for the two sets of variables.
1505     //   Bx -          The well-measured variable gradients.
1506     //   By -          The poorly-measured variable gradients.
1507     //
1508     {
1509       double sig2 = scalar(Bx.T() * Q * Bx);
1510 
1511       if (By.num_row() > 0) {
1512         sig2 += scalar(By.T() * R * By);
1513         sig2 += 2 * scalar(Bx.T() * S * By);
1514       }
1515 
1516       assert(sig2 >= 0);
1517       return sqrt(sig2);
1518     }
1519 
1520     /**
1521     @brief Helper function: Calculate the final requested mass and
1522     its uncertainty.
1523 
1524     @param ev The event we are fitting.
1525 
1526     @param mass_constraint The description of the mass to be calculated.
1527 
1528     @param args Parameter settings.
1529 
1530     @param chisq The \f$\chi^{2}\f$ from the fit.
1531 
1532     @param Q The final error matrix for the well-measured variables.
1533 
1534     @param R The final error matrix for the poorly-measured variables.
1535 
1536     @param S The final cross error matrix for the two sets of variables.
1537 
1538     @param m The mass (output).
1539 
1540     @param sigm The uncertainty on the mass (output).
1541  */
1542     void calculate_mass(Fourvec_Event& ev,
1543                         const vector<Constraint>& mass_constraint,
1544                         const Fourvec_Constrainer_Args& args,
1545                         double chisq,
1546                         const Matrix& Q,
1547                         const Matrix& R,
1548                         const Matrix& S,
1549                         double& m,
1550                         double& sigm)
1551     //
1552     // Purpose: Calculate the final requested mass and its uncertainty.
1553     //
1554     // Inputs:
1555     //   ev -          The event we're fitting.
1556     //   mass_constraint- The description of the mass we're to calculate.
1557     //   args -        Parameter settings.
1558     //   chisq -       The chisq from the fit.
1559     //   Q -           The final error matrix for the well-measured variables.
1560     //   R -           The final error matrix for the poorly-measured variables.
1561     //   S -           The final cross error matrix for the two sets of variables.
1562     //
1563     // Outputs:
1564     //   m -           The mass.
1565     //   sigm -        Its uncertainty.
1566     //
1567     {
1568       // Don't do anything if the mass wasn't specified.
1569       if (mass_constraint.empty()) {
1570         m = 0;
1571         sigm = 0;
1572         return;
1573       }
1574 
1575       // Do the constraint calculation.
1576       int n_measured_vars = ev.nobjs() * 3 + 2;
1577       int n_unmeasured_vars = 0;
1578       if (ev.has_neutrino())
1579         n_unmeasured_vars = 1;
1580 
1581       Row_Vector F(1);
1582       Matrix Bx(n_measured_vars, 1);
1583       Matrix By(n_unmeasured_vars, 1);
1584 
1585       Fourvec_Constraint_Calculator cc(ev, mass_constraint, args);
1586       cc.calculate_constraints(F, Bx, By);
1587 
1588       // Calculate the mass.
1589       //assert (F(1) >= 0);
1590       if (F(1) >= 0.)
1591         m = sqrt(F(1) * 2);
1592       else {
1593         m = 0.;
1594         chisq = -100.;
1595       }
1596 
1597       // And the uncertainty.
1598       // We can only do this if the fit converged.
1599       if (chisq < 0)
1600         sigm = 0;
1601       else {
1602         //assert (F(1) > 0);
1603         Bx = Bx / (m);
1604         By = By / (m);
1605 
1606         sigm = calculate_sigm(Q, R, S, Bx, By);
1607       }
1608     }
1609 
1610   }  // unnamed namespace
1611 
1612   double Fourvec_Constrainer::constrain(
1613       Fourvec_Event& ev, double& m, double& sigm, Column_Vector& pullx, Column_Vector& pully)
1614   //
1615   // Purpose: Do a constrained fit for EV.  Returns the requested mass and
1616   //          its error in M and SIGM, and the pull quantities in PULLX and
1617   //          PULLY.  Returns the chisq; this will be < 0 if the fit failed
1618   //          to converge.
1619   //
1620   // Inputs:
1621   //   ev -          The event we're fitting.
1622   //
1623   // Outputs:
1624   //   ev -          The fitted event.
1625   //   m -           Requested invariant mass.
1626   //   sigm -        Uncertainty on m.
1627   //   pullx -       Pull quantities for well-measured variables.
1628   //   pully -       Pull quantities for poorly-measured variables.
1629   //
1630   // Returns:
1631   //   The fit chisq, or < 0 if the fit didn't converge.
1632   //
1633   {
1634     adjust_fourvecs(ev, _args.use_e());
1635 
1636     bool use_kt = ev.has_neutrino() || !_args.ignore_met();
1637     int nobjs = ev.nobjs();
1638     int n_measured_vars = nobjs * 3;
1639     int n_unmeasured_vars = 0;
1640     int n_constraints = _constraints.size();
1641 
1642     if (use_kt) {
1643       n_measured_vars += 2;
1644 
1645       if (ev.has_neutrino())
1646         n_unmeasured_vars = 1;
1647       else
1648         n_constraints += 2;
1649     }
1650 
1651     Matrix G_i(n_measured_vars, n_measured_vars);
1652     Diagonal_Matrix Y(n_unmeasured_vars);
1653     Column_Vector x(n_measured_vars);
1654     Column_Vector y(n_unmeasured_vars);
1655     pack_event(ev, _args.use_e(), use_kt, x, y, G_i, Y);
1656 
1657     Column_Vector xm = x;
1658     Column_Vector ym = y;
1659 
1660     // ??? Should allow for using a different underlying fitter here.
1661     Chisq_Constrainer fitter(_args.chisq_constrainer_args());
1662 
1663     Fourvec_Constraint_Calculator cc(ev, _constraints, _args);
1664 
1665     Matrix Q;
1666     Matrix R;
1667     Matrix S;
1668     double chisq = fitter.fit(cc, xm, x, ym, y, G_i, Y, pullx, pully, Q, R, S);
1669 
1670     unpack_event(ev, x, y, _args.use_e(), use_kt);
1671 
1672     calculate_mass(ev, _mass_constraint, _args, chisq, Q, R, S, m, sigm);
1673 
1674     return chisq;
1675   }
1676 
1677 }  // namespace hitfit