![]() |
|
|||
File indexing completed on 2024-04-06 12:31:18
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 1641 if (use_kt) { 1642 n_measured_vars += 2; 1643 1644 if (ev.has_neutrino()) 1645 n_unmeasured_vars = 1; 1646 } 1647 1648 Matrix G_i(n_measured_vars, n_measured_vars); 1649 Diagonal_Matrix Y(n_unmeasured_vars); 1650 Column_Vector x(n_measured_vars); 1651 Column_Vector y(n_unmeasured_vars); 1652 pack_event(ev, _args.use_e(), use_kt, x, y, G_i, Y); 1653 1654 Column_Vector xm = x; 1655 Column_Vector ym = y; 1656 1657 // ??? Should allow for using a different underlying fitter here. 1658 Chisq_Constrainer fitter(_args.chisq_constrainer_args()); 1659 1660 Fourvec_Constraint_Calculator cc(ev, _constraints, _args); 1661 1662 Matrix Q; 1663 Matrix R; 1664 Matrix S; 1665 double chisq = fitter.fit(cc, xm, x, ym, y, G_i, Y, pullx, pully, Q, R, S); 1666 1667 unpack_event(ev, x, y, _args.use_e(), use_kt); 1668 1669 calculate_mass(ev, _mass_constraint, _args, chisq, Q, R, S, m, sigm); 1670 1671 return chisq; 1672 } 1673 1674 } // namespace hitfit
[ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
This page was automatically generated by the 2.2.1 LXR engine. The LXR team |
![]() ![]() |