CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Fourvec_Constrainer.cc
Go to the documentation of this file.
1 //
2 //
3 // File: src/Fourvec_Constrainer.cc
4 // Purpose: Do a kinematic fit for a set of 4-vectors, given a set
5 // of mass constraints.
6 // Created: Jul, 2000, sss, based on run 1 mass analysis code.
7 //
8 // CMSSW File : src/Fourvec_Constrainer.cc
9 // Original Author : Scott Stuart Snyder <snyder@bnl.gov> for D0
10 // Imported to CMSSW by Haryo Sumowidagdo <Suharyo.Sumowidagdo@cern.ch>
11 //
12 
13 
46 #include <cmath>
47 #include <iostream>
48 
49 
50 using std::sqrt;
51 using std::exp;
52 using std::cos;
53 using std::sin;
54 using std::ostream;
55 using std::vector;
56 
57 
58 namespace hitfit {
59 
60 
61 //*************************************************************************
62 // Argument handling.
63 //
64 
65 
67 //
68 // Purpose: Constructor.
69 //
70 // Inputs:
71 // defs - The Defaults instance from which to initialize.
72 //
73  : _use_e (defs.get_bool ("use_e")),
74  _e_com (defs.get_float ("e_com")),
75  _ignore_met (defs.get_bool ("ignore_met")),
76  _chisq_constrainer_args (defs)
77 {
78 }
79 
80 
82 //
83 // Purpose: Return the use_e parameter.
84 // See the header for documentation.
85 //
86 {
87  return _use_e;
88 }
89 
90 
92 //
93 // Purpose: Return the e_com parameter.
94 // See the header for documentation.
95 //
96 {
97  return _e_com;
98 }
99 
100 
102 //
103 // Purpose: Return the ignore_met parameter.
104 // See the header for documentation.
105 //
106 {
107  return _ignore_met;
108 }
109 
110 
113 //
114 // Purpose: Return the contained subobject parameters.
115 //
116 {
118 }
119 
120 
121 //*************************************************************************
122 // Variable layout.
123 //
124 // We need to map the quantities we fit onto the vectors of well- and
125 // poorly-measured quantities.
126 //
127 
128 //
129 // The well-measured variables consist of three variables for each
130 // object. If we are using transverse momentum constraints,
131 // these fill be followed by the two cartesian components of kt.
132 //
133 // Each object is represented by three variables: the momentum (or 1/p
134 // if the muon flag was set), and the two spherical angles, phi and eta.
135 // Here is how they're ordered.
136 //
140 typedef enum {
141  p_offs = 0,
142  phi_offs = 1,
144 } Offsets;
145 
150 typedef enum {
151  x_offs = 0,
152  y_offs = 1
153 } Kt_Offsets;
154 
155 //
156 // If there is a neutrino, then it is at index 1 of the poorly-measured
157 // set (otherwise, that set is empty).
158 //
163 typedef enum {
164  nu_z = 1
166 
167 
168 
169 namespace {
170 
171 
178 int obj_index (int i)
179 //
180 // Purpose: Return the starting variable index for object I.
181 //
182 // Inputs:
183 // i - The object index.
184 //
185 // Returns:
186 // The index in the well-measured set of the first variable
187 // for object I.
188 //
189 {
190  return i*3 + 1;
191 }
192 
193 
194 } // unnamed namespace
195 
196 
197 //*************************************************************************
198 // Object management.
199 //
200 
201 
203 //
204 // Purpose: Constructor.
205 //
206 // Inputs:
207 // args - The parameter settings for this instance.
208 //
209  : _args (args)
210 {
211 }
212 
213 
215 //
216 // Purpose: Specify an additional constraint S for the problem.
217 // The format for S is described in the header.
218 //
219 // Inputs:
220 // s - The constraint to add.
221 //
222 {
223  _constraints.push_back (Constraint (s));
224 }
225 
226 
228 //
229 // Purpose: Specify the combination of objects that will be returned by
230 // constrain() as the mass. The format of S is the same as for
231 // normal constraints. The LHS specifies the mass to calculate;
232 // the RHS should be zero.
233 // This should only be called once.
234 //
235 // Inputs:
236 // s - The constraint defining the mass.
237 //
238 {
239  assert (_mass_constraint.size() == 0);
240  _mass_constraint.push_back (Constraint (s));
241 }
242 
243 
253 std::ostream& operator<< (std::ostream& s, const Fourvec_Constrainer& c)
254 //
255 // Purpose: Print the object to S.
256 //
257 // Inputs:
258 // s - The stream to which to write.
259 // c - The object to write.
260 //
261 // Returns:
262 // The stream S.
263 //
264 {
265  s << "Constraints: (e_com = " << c._args.e_com() << ") ";
266  if (c._args.use_e())
267  s << "(E)";
268  s << "\n";
269 
270  for (std::vector<Constraint>::size_type i=0; i < c._constraints.size(); i++)
271  s << " " << c._constraints[i] << "\n";
272 
273  if (c._mass_constraint.size() > 0) {
274  s << "Mass constraint:\n";
275  s << c._mass_constraint[0] << "\n";
276  }
277  return s;
278 }
279 
280 
281 //*************************************************************************
282 // Event packing and unpacking.
283 //
284 
285 
286 namespace {
287 
288 
299 void adjust_fourvecs (Fourvec_Event& ev,
300  bool use_e_flag)
301 //
302 // Purpose: For all objects in EV, adjust their 4-momenta
303 // to have their requested masses.
304 //
305 // Inputs:
306 // ev - The event on which to operate.
307 // use_e_flag - If true, keep E and scale 3-momentum.
308 // If false, keep the 3-momentum and scale E.
309 //
310 {
311  int nobjs = ev.nobjs ();
312  for (int i=0; i < nobjs; i++) {
313  const FE_Obj& obj = ev.obj (i);
314  Fourvec p = obj.p;
315  if (use_e_flag)
316  adjust_p_for_mass (p, obj.mass);
317  else
318  adjust_e_for_mass (p, obj.mass);
319  ev.set_obj_p (i, p);
320  }
321 }
322 
323 
340 Fourvec get_p_eta_phi_vec (const Column_Vector& c,
341  int ndx,
342  const FE_Obj& obj,
343  bool use_e_flag)
344 //
345 // Purpose: Convert object NDX from its representation in the set
346 // of well-measured variables C to a 4-vector.
347 //
348 // Inputs:
349 // c - The vector of well-measured variables.
350 // ndx - The index of the object in which we're interested.
351 // obj - The object from the Fourvec_Event.
352 // use_e_flag - If true, we're using E as the fit variable, otherwise p.
353 //
354 // Returns:
355 // The object's 4-momentum.
356 //
357 {
358  // Get the energy and momentum of the object.
359  double e, p;
360 
361  if (use_e_flag) {
362  // We're using E as a fit variable. Get it directly.
363  e = c(ndx + p_offs);
364 
365  // Take into account the muon case.
366  if (obj.muon_p) e = 1/e;
367 
368  // Find the momentum given the energy.
369  if (obj.mass == 0)
370  p = e;
371  else {
372  double xx = e*e - obj.mass*obj.mass;
373  if (xx >= 0)
374  p = sqrt (xx);
375  else
376  p = 0;
377  }
378  }
379  else {
380  // We're using P as a fit variable. Fetch it.
381  p = c(ndx + p_offs);
382 
383  // Take into account the muon case.
384  if (obj.muon_p) p = 1/p;
385 
386  // Find the energy given the momentum.
387  e = (obj.mass == 0 ? p : sqrt (obj.mass*obj.mass + p*p));
388  }
389 
390  // Get angular variables.
391  double phi = c(ndx + phi_offs);
392  double eta = c(ndx + eta_offs);
393  if (fabs (eta) > 50) {
394  // Protect against ridiculously large etas
395  eta = eta > 0 ? 50 : -50;
396  }
397  double exp_eta = exp (eta);
398  double iexp_eta = 1/exp_eta;
399  double sin_theta = 2 / (exp_eta + iexp_eta);
400  double cos_theta = (exp_eta - iexp_eta) / (exp_eta + iexp_eta);
401 
402  // Form the 4-momentum.
403  return Fourvec (p * sin_theta * cos (phi),
404  p * sin_theta * sin (phi),
405  p * cos_theta,
406  e);
407 }
408 
409 
426 void set_p_eta_phi_vec (const FE_Obj& obj,
427  Column_Vector& c,
428  int ndx,
429  bool use_e_flag)
430 //
431 // Purpose: Initialize the variables in the well-measured set C describing
432 // object NDX from its Fourvec_Event representation OBJ.
433 //
434 // Inputs:
435 // obj - The object from the Fourvec_Event.
436 // c - The vector of well-measured variables.
437 // ndx - The index of the object in which we're interested.
438 // use_e_flag - If true, we're using E as the fit variable, otherwise p.
439 //
440 //
441 {
442  if (use_e_flag)
443  c(ndx + p_offs) = obj.p.e();
444  else
445  c(ndx + p_offs) = obj.p.vect().mag();
446 
447  if (obj.muon_p) c(ndx + p_offs) = 1/c(ndx + p_offs);
448  c(ndx + phi_offs) = obj.p.phi();
449  c(ndx + eta_offs) = obj.p.pseudoRapidity();
450 }
451 
452 
462 Matrix error_matrix (double p_sig,
463  double phi_sig,
464  double eta_sig)
465 //
466 // Purpose: Set up the 3x3 error matrix for an object.
467 //
468 // Inputs:
469 // p_sig - The momentum uncertainty.
470 // phi_sig - The phi uncertainty.
471 // eta_sig - The eta uncertainty.
472 //
473 // Returns:
474 // The object's error matrix.
475 //
476 {
477  Matrix err (3, 3, 0);
478  err(1+ p_offs, 1+ p_offs) = p_sig * p_sig;
479  err(1+phi_offs, 1+phi_offs) = phi_sig * phi_sig;
480  err(1+eta_offs, 1+eta_offs) = eta_sig * eta_sig;
481  return err;
482 }
483 
484 
510 void pack_event (const Fourvec_Event& ev,
511  bool use_e_flag,
512  bool use_kt_flag,
513  Column_Vector& xm,
514  Column_Vector& ym,
515  Matrix& G_i,
517 //
518 // Purpose: Take the information from a Fourvec_Event EV and pack
519 // it into vectors of well- and poorly-measured variables.
520 // Also set up the error matrices.
521 //
522 // Inputs:
523 // ev - The event to pack.
524 // use_e_flag - If true, we're using E as the fit variable, otherwise p.
525 // use_kt_flag - True if we're to pack kt variables.
526 //
527 // Outputs:
528 // xm - Vector of well-measured variables.
529 // ym - Vector of poorly-measured variables.
530 // G_i - Error matrix for well-measured variables.
531 // Y - Inverse error matrix for poorly-measured variables.
532 //
533 {
534  // Number of objects in the event.
535  int nobjs = ev.nobjs ();
536 
537  int n_measured_vars = nobjs * 3;
538  if (use_kt_flag)
539  n_measured_vars += 2;
540 
541  // Clear the error matrix.
542  G_i = Matrix (n_measured_vars, n_measured_vars, 0);
543 
544  // Loop over objects.
545  for (int i=0; i<nobjs; i++) {
546  const FE_Obj& obj = ev.obj (i);
547  int this_index = obj_index (i);
548  set_p_eta_phi_vec (obj, xm, this_index, use_e_flag);
549  G_i.sub (this_index, this_index, error_matrix (obj.p_error,
550  obj.phi_error,
551  obj.eta_error));
552 
553  }
554 
555  if (use_kt_flag) {
556  // Set up kt.
557  int kt_ndx = obj_index (nobjs);
558  xm (kt_ndx+x_offs) = ev.kt().x();
559  xm (kt_ndx+y_offs) = ev.kt().y();
560 
561  // And its error matrix.
562  G_i(kt_ndx+x_offs, kt_ndx+x_offs) = ev.kt_x_error() * ev.kt_x_error();
563  G_i(kt_ndx+y_offs, kt_ndx+y_offs) = ev.kt_y_error() * ev.kt_y_error();
564  G_i(kt_ndx+x_offs, kt_ndx+y_offs) = ev.kt_xy_covar();
565  G_i(kt_ndx+y_offs, kt_ndx+x_offs) = ev.kt_xy_covar();
566  }
567 
568  // Handle a neutrino.
569  if (ev.has_neutrino()) {
570  ym(nu_z) = ev.nu().z();
571  Y = Diagonal_Matrix (1, 0);
572  }
573 }
574 
575 
595 void unpack_event (Fourvec_Event& ev,
596  const Column_Vector& x,
597  const Column_Vector& y,
598  bool use_e_flag,
599  bool use_kt_flag)
600 //
601 // Purpose: Update the contents of EV from the variable sets X and Y.
602 //
603 // Inputs:
604 // ev - The event.
605 // x - Vector of well-measured variables.
606 // use_e_flag - If true, we're using E as the fit variable, otherwise p.
607 // use_kt_flag - True if we're too unpack kt variables.
608 //
609 // Outputs:
610 // ev - The event after updating.
611 //
612 {
613  // Do all the objects.
614  Fourvec sum = Fourvec(0);
615  for (int j=0; j<ev.nobjs(); j++) {
616  const FE_Obj& obj = ev.obj (j);
617  ev.set_obj_p (j, get_p_eta_phi_vec (x, obj_index (j), obj, use_e_flag));
618  sum += obj.p;
619  }
620 
621 
622  if (use_kt_flag) {
623  int kt_ndx = obj_index (ev.nobjs());
624  Fourvec kt = Fourvec (x(kt_ndx+x_offs), x(kt_ndx+y_offs), 0, 0);
625  Fourvec nu = kt - sum;
626  if (ev.has_neutrino()) {
627  nu.setPz (y(nu_z));
628  adjust_e_for_mass (nu, 0);
629  ev.set_nu_p (nu);
630  }
631  else {
632  adjust_e_for_mass (nu, 0);
633  ev.set_x_p (nu);
634  }
635  }
636 }
637 
638 
639 } // unnamed namespace
640 
641 
642 //*************************************************************************
643 // Constraint evaluation.
644 //
645 
646 
647 namespace {
648 
649 
678 double dot_and_gradient (const Fourvec& v1,
679  const Fourvec& v2,
680  bool use_e_flag,
681  double v1_x[3],
682  double v2_x[3],
683  bool& badflag)
684 //
685 // Purpose: Compute the dot product v1.v2 and its gradients wrt
686 // p, phi, and theta of each 4-vector.
687 //
688 // Inputs:
689 // v1 - The first 4-vector in the dot product.
690 // v2 - The second 4-vector in the dot product.
691 // use_e_flag - If true, we're using E as the fit variable, otherwise p.
692 //
693 // Outputs:
694 // v1_x - Gradients of the dot product wrt v1's p, phi, theta.
695 // v2_x - Gradients of the dot product wrt v2's p, phi, theta.
696 // badflag - Set to true for the singular case (vectors vanish).
697 //
698 // Returns:
699 // The dot product.
700 //
701 {
702  // Calculate the dot product.
703  double dot = v1 * v2;
704 
705  double p1 = v1.vect().mag();
706  double p2 = v2.vect().mag();
707  double e1 = v1.e();
708  double e2 = v2.e();
709  double pt1 = v1.vect().perp();
710  double pt2 = v2.vect().perp();
711 
712  // Protect against the singular case.
713  badflag = false;
714  if (p1 == 0 || p2 == 0 || e1 == 0 || e2 == 0 || pt1 == 0 || pt2 == 0) {
715  badflag = true;
716  v1_x[p_offs] = v1_x[phi_offs] = v1_x[eta_offs] = 0;
717  v2_x[p_offs] = v2_x[phi_offs] = v2_x[eta_offs] = 0;
718  return false;
719  }
720 
721  // Calculate the gradients.
722  v1_x[p_offs] = (dot - v1.m2() * e2 / e1) / p1;
723  v2_x[p_offs] = (dot - v2.m2() * e1 / e2) / p2;
724 
725  if (use_e_flag) {
726  v1_x[p_offs] *= e1 / p1;
727  v2_x[p_offs] *= e2 / p2;
728  }
729 
730  v1_x[phi_offs] = v1(1)*v2(0) - v1(0)*v2(1);
731  v2_x[phi_offs] = -v1_x[phi_offs];
732 
733  double fac = v1(0)*v2(0) + v1(1)*v2(1);
734  v1_x[eta_offs] = pt1*v2(2) - v1(2)/pt1 * fac;
735  v2_x[eta_offs] = pt2*v1(2) - v2(2)/pt2 * fac;
736 
737  return dot;
738 }
739 
740 
762 void addin_obj_gradient (int constraint_no,
763  int sign,
764  int base_index,
765  const double grad[],
766  Matrix& Bx)
767 //
768 // Purpose: Tally up the dot product gradients for an object
769 // into Bx.
770 //
771 // Inputs:
772 // constraint_no-The number of the constraint.
773 // base_index - The index in the well-measured variable list
774 // of the first variable for this object.
775 // sign - The sign with which these gradients should be
776 // added into Bx, either +1 or -1. (I.e., which
777 // side of the constraint equation.)
778 // grad - The gradients for this object, vs. p, phi, theta.
779 // Bx - The well-measured variable gradients.
780 //
781 // Outputs:
782 // Bx - The well-measured variable gradients, updated.
783 //
784 {
785  Bx(base_index + p_offs, constraint_no) += sign * grad[p_offs];
786  Bx(base_index + phi_offs, constraint_no) += sign * grad[phi_offs];
787  Bx(base_index + eta_offs, constraint_no) += sign * grad[eta_offs];
788 }
789 
790 
818 void addin_nu_gradient (int constraint_no,
819  int sign,
820  int kt_ndx,
821  const double grad[],
822  Matrix& Bx, Matrix& By)
823 //
824 // Purpose: Tally up the dot product gradients for a neutrino
825 // into Bx and By.
826 //
827 // Inputs:
828 // constraint_no-The number of the constraint.
829 // sign - The sign with which these gradients should be
830 // added into Bx, either +1 or -1. (I.e., which
831 // side of the constraint equation.)
832 // kt_ndx - The index of the kt variables in the variables array.
833 // grad - The gradients for this object, vs. p, phi, theta.
834 // Bx - The well-measured variable gradients.
835 // By - The poorly-measured variable gradients.
836 //
837 // Outputs:
838 // Bx - The well-measured variable gradients, updated.
839 // By - The poorly-measured variable gradients, updated.
840 //
841 {
842  Bx(kt_ndx+x_offs,constraint_no) += sign*grad[p_offs]; // Really p for now.
843  Bx(kt_ndx+y_offs,constraint_no) += sign*grad[phi_offs];// Really phi for now.
844  By(nu_z, constraint_no) += sign*grad[eta_offs]; // Really theta ...
845 }
846 
847 
874 void addin_gradient (const Fourvec_Event& ev,
875  int constraint_no,
876  int sign,
877  int obj_no,
878  const double grad[],
879  Matrix& Bx,
880  Matrix& By)
881 //
882 // Purpose: Tally up the dot product gradients for an object (which may
883 // or may not be a neutrino) into Bx and By.
884 //
885 // Inputs:
886 // ev - The event we're fitting.
887 // constraint_no-The number of the constraint.
888 // sign - The sign with which these gradients should be
889 // added into Bx, either +1 or -1. (I.e., which
890 // side of the constraint equation.)
891 // obj_no - The number of the object.
892 // grad - The gradients for this object, vs. p, phi, theta.
893 // Bx - The well-measured variable gradients.
894 // By - The poorly-measured variable gradients.
895 //
896 // Outputs:
897 // Bx - The well-measured variable gradients, updated.
898 // By - The poorly-measured variable gradients, updated.
899 //
900 {
901  if (obj_no >= ev.nobjs()) {
902  assert (obj_no == ev.nobjs());
903  addin_nu_gradient (constraint_no, sign, obj_index (obj_no), grad, Bx, By);
904  }
905  else
906  addin_obj_gradient (constraint_no, sign, obj_index (obj_no), grad, Bx);
907 }
908 
909 
940 void addin_gradients (const Fourvec_Event& ev,
941  int constraint_no,
942  int sign,
943  int i,
944  const double igrad[],
945  int j,
946  const double jgrad[],
947  Matrix& Bx,
948  Matrix& By)
949 //
950 // Purpose: Tally up the gradients from a single dot product into Bx and By.
951 //
952 // Inputs:
953 // ev - The event we're fitting.
954 // constraint_no-The number of the constraint.
955 // sign - The sign with which these gradients should be
956 // added into Bx, either +1 or -1. (I.e., which
957 // side of the constraint equation.)
958 // i - The number of the first object.
959 // igrad - The gradients for the first object, vs. p, phi, theta.
960 // j - The number of the second object.
961 // jgrad - The gradients for the second object, vs. p, phi, theta.
962 // Bx - The well-measured variable gradients.
963 // By - The poorly-measured variable gradients.
964 //
965 // Outputs:
966 // Bx - The well-measured variable gradients, updated.
967 // By - The poorly-measured variable gradients, updated.
968 //
969 {
970  addin_gradient (ev, constraint_no, sign, i, igrad, Bx, By);
971  addin_gradient (ev, constraint_no, sign, j, jgrad, Bx, By);
972 }
973 
974 
987 void add_mass_terms (const Fourvec_Event& ev,
988  const vector<Constraint>& constraints,
989  Row_Vector& F)
990 //
991 // Purpose: Add the m^2 terms into the constraint values.
992 //
993 // Inputs:
994 // ev - The event we're fitting.
995 // constraints - The list of constraints.
996 // F - Vector of constraint values.
997 //
998 // Outputs:
999 // F - Vector of constraint values, updated.
1000 //
1001 {
1002  for (std::vector<Constraint>::size_type i=0; i<constraints.size(); i++)
1003  F(i+1) += constraints[i].sum_mass_terms (ev);
1004 }
1005 
1006 
1040 bool calculate_mass_constraints (const Fourvec_Event& ev,
1041  const Pair_Table& pt,
1042  const vector<Constraint>& constraints,
1043  bool use_e_flag,
1044  Row_Vector& F,
1045  Matrix& Bx,
1046  Matrix& By)
1047 //
1048 // Purpose: Calculate the mass constraints and gradients.
1049 // Note: At this stage, the gradients are calculated not
1050 // quite with respect to the fit variables; instead, for
1051 // all objects (including the neutrino) we calculate
1052 // the gradients with respect to p, phi, theta. They'll
1053 // be converted via appropriate Jacobian transformations
1054 // later.
1055 //
1056 // Inputs:
1057 // ev - The event we're fitting.
1058 // pt - Table of cached pair assignments.
1059 // constraints - The list of constraints.
1060 // use_e_flag - If true, we're using E as the fit variable, otherwise p.
1061 //
1062 // Outputs:
1063 // (nb. these should be passed in correctly dimensioned.)
1064 // F - Vector of constraint values.
1065 // Bx - The well-measured variable gradients.
1066 // By - The poorly-measured variable gradients.
1067 //
1068 {
1069  int npairs = pt.npairs ();
1070  for (int p=0; p < npairs; p++) {
1071  const Objpair& objpair = pt.get_pair (p);
1072  int i = objpair.i();
1073  int j = objpair.j();
1074  double igrad[3], jgrad[3];
1075  bool badflag = false;
1076  double dot = dot_and_gradient (ev.obj (i).p,
1077  ev.obj (j).p,
1078  use_e_flag,
1079  igrad,
1080  jgrad,
1081  badflag);
1082  if (badflag)
1083  return false;
1084 
1085  for (std::vector<Constraint>::size_type k=0; k < constraints.size(); k++)
1086  if (objpair.for_constraint (k)) {
1087  F(k+1) += objpair.for_constraint (k) * dot;
1088  addin_gradients (ev, k+1, objpair.for_constraint (k),
1089  i, igrad, j, jgrad, Bx, By);
1090  }
1091 
1092  }
1093 
1094  add_mass_terms (ev, constraints, F);
1095  return true;
1096 }
1097 
1098 
1121 void add_nuterm (unsigned ndx,
1122  const Fourvec& v,
1123  Matrix& Bx,
1124  bool use_e_flag,
1125  int kt_ndx)
1126 //
1127 // Purpose: Carry out the Jacobian transformation for
1128 // (p_nu^x,p_nu_y) -> (kt^x, kt_y) for a given object.
1129 //
1130 // Inputs:
1131 // ndx - Index of the object for which to transform gradients.
1132 // v - The object's 4-momentum.
1133 // Bx - The well-measured variable gradients.
1134 // use_e_flag - If true, we're using E as the fit variable, otherwise p.
1135 // kt_ndx - The index of the kt variables in the variables array.
1136 //
1137 // Outputs:
1138 // Bx - The well-measured variable gradients, updated.
1139 //
1140 {
1141  double px = v.px();
1142  double py = v.py();
1143  double cot_theta = v.pz() / v.vect().perp();
1144 
1145  for (int j=1; j<=Bx.num_col(); j++) {
1146  double dxnu = Bx(kt_ndx+x_offs, j);
1147  double dynu = Bx(kt_ndx+y_offs, j);
1148 
1149  if (dxnu != 0 || dynu != 0) {
1150  double fac = 1 / v.vect().mag();
1151  if (use_e_flag)
1152  fac = v.e() * fac * fac;
1153  Bx(ndx + p_offs, j) -= (px*dxnu + py*dynu) * fac;
1154  Bx(ndx + phi_offs, j) += py*dxnu - px*dynu;
1155  Bx(ndx + eta_offs, j) -= (px*dxnu + py*dynu) * cot_theta;
1156  }
1157  }
1158 }
1159 
1160 
1184 void convert_neutrino (const Fourvec_Event& ev,
1185  bool use_e_flag,
1186  Matrix& Bx,
1187  Matrix& By)
1188 //
1189 // Purpose: Carry out the Jacobian transformations the neutrino.
1190 // First, convert from spherical (p, phi, theta) coordinates
1191 // to rectangular (x, y, z). Then convert from neutrino pt
1192 // components to kt components.
1193 //
1194 // Inputs:
1195 // ev - The event we're fitting.
1196 // use_e_flag - If true, we're using E as the fit variable, otherwise p.
1197 // Bx - The well-measured variable gradients.
1198 // By - The poorly-measured variable gradients.
1199 //
1200 // Outputs:
1201 // Bx - The well-measured variable gradients, updated.
1202 // By - The poorly-measured variable gradients, updated.
1203 //
1204 {
1205  int nconstraints = Bx.num_col ();
1206 
1207  const Fourvec& nu = ev.nu ();
1208 
1209  // convert neutrino from polar coordinates to rectangular.
1210  double pnu2 = nu.vect().mag2(); double pnu = sqrt (pnu2);
1211  double ptnu2 = nu.vect().perp2(); double ptnu = sqrt (ptnu2);
1212 
1213  // Doesn't matter whether we use E or P here, since nu is massless.
1214 
1215  double thfac = nu.z()/pnu2/ptnu;
1216  double fac[3][3];
1217  fac[0][0] = nu(0)/pnu; fac[0][1] = nu(1)/pnu; fac[0][2] = nu(2)/pnu;
1218  fac[1][0] = - nu(1)/ptnu2; fac[1][1] = nu(0)/ptnu2; fac[1][2] = 0;
1219  fac[2][0] = nu(0)*thfac; fac[2][1] = nu(1)*thfac; fac[2][2] = -ptnu/pnu2;
1220 
1221  int kt_ndx = obj_index (ev.nobjs());
1222  for (int j=1; j<=nconstraints; j++) {
1223  double tmp1 = fac[0][0]*Bx(kt_ndx+x_offs,j) +
1224  fac[1][0]*Bx(kt_ndx+y_offs,j) +
1225  fac[2][0]*By(nu_z,j);
1226  Bx(kt_ndx+y_offs,j) = fac[0][1]*Bx(kt_ndx+x_offs,j) +
1227  fac[1][1]*Bx(kt_ndx+y_offs,j) +
1228  fac[2][1]*By(nu_z,j);
1229  By(nu_z,j) = fac[0][2]*Bx(kt_ndx+x_offs,j) + fac[2][2]*By(nu_z,j);
1230 
1231  Bx(kt_ndx+x_offs,j) = tmp1;
1232  }
1233 
1234  // Add nu terms.
1235  for (int j=0; j<ev.nobjs(); j++) {
1236  add_nuterm (obj_index (j), ev.obj(j).p, Bx, use_e_flag, kt_ndx);
1237  }
1238 }
1239 
1240 
1261 void calculate_kt_constraints (const Fourvec_Event& ev,
1262  bool use_e_flag,
1263  Row_Vector& F,
1264  Matrix& Bx)
1265 //
1266 // Purpose: Calculate the overall kt constraints (kt = 0) for the case
1267 // where there is no neutrino.
1268 //
1269 // Inputs:
1270 // ev - The event we're fitting.
1271 // use_e_flag - If true, we're using E as the fit variable, otherwise p.
1272 // F - Vector of constraint values.
1273 // Bx - The well-measured variable gradients.
1274 //
1275 // Outputs:
1276 // F - Vector of constraint values, updated.
1277 // Bx - The well-measured variable gradients, updated.
1278 //
1279 {
1280  Fourvec tmp = Fourvec(0);
1281  int base = F.num_col() - 2;
1282  int nobjs = ev.nobjs();
1283  for (int j=0; j<nobjs; j++) {
1284  const Fourvec& obj = ev.obj (j).p;
1285  tmp += obj;
1286 
1287  int ndx = obj_index (j);
1288  double p = obj.vect().mag();
1289  double cot_theta = obj.z() / obj.vect().perp();
1290 
1291  Bx(ndx + p_offs, base+1) = obj(0) / p;
1292  Bx(ndx + phi_offs, base+1) = -obj(1);
1293  Bx(ndx + eta_offs, base+1) = obj(0) * cot_theta;
1294 
1295  Bx(ndx + p_offs, base+2) = obj(1) / p;
1296  Bx(ndx + phi_offs, base+2) = obj(0);
1297  Bx(ndx + eta_offs, base+2) = obj(1) * cot_theta;
1298 
1299  if (use_e_flag) {
1300  Bx(ndx + p_offs, base+1) *= obj.e() / p;
1301  Bx(ndx + p_offs, base+2) *= obj.e() / p;
1302  }
1303  }
1304 
1305  int kt_ndx = obj_index (nobjs);
1306  Bx(kt_ndx+x_offs, base+1) = -1;
1307  Bx(kt_ndx+y_offs, base+2) = -1;
1308 
1309  F(base+1) = tmp(0) - ev.kt().x ();
1310  F(base+2) = tmp(1) - ev.kt().y ();
1311 }
1312 
1313 
1327 void ddtheta_to_ddeta (int i, double cos_theta, Matrix& Bx)
1328 //
1329 // Purpose: Do the Jacobian transformation from theta -> eta
1330 // for a single object.
1331 //
1332 // Inputs:
1333 // i - The index of the object.
1334 // cos_theta - cos(theta) for the object.
1335 // Bx - The well-measured variable gradients.
1336 //
1337 // Outputs:
1338 // Bx - The well-measured variable gradients, updated.
1339 //
1340 {
1341  double sin_theta = sqrt (1 - cos_theta * cos_theta);
1342  for (int j=1; j<=Bx.num_col(); j++)
1343  Bx(i,j) *= - sin_theta; /* \sin\theta = 1 / \cosh\eta */
1344 }
1345 
1346 
1347 } // unnamed namespace
1348 
1349 
1376  : public Constraint_Calculator
1377 //
1378 // Purpose: Constraint evaluator.
1379 //
1380 {
1381 public:
1382  // Constructor, destructor.
1384  const vector<Constraint>& constraints,
1387 
1388  // Evaluate constraints at the point described by X and Y (well-measured
1389  // and poorly-measured variables, respectively). The results should
1390  // be stored in F. BX and BY should be set to the gradients of F with
1391  // respect to X and Y, respectively.
1392  //
1393  // Return true if the point X, Y is accepted.
1394  // Return false if it is rejected (i.e., in an unphysical region).
1395  // The constraints need not be evaluated in that case.
1396  virtual bool eval (const Column_Vector& x,
1397  const Column_Vector& y,
1398  Row_Vector& F,
1399  Matrix& Bx,
1400  Matrix& By) override;
1401 
1402 
1403  // Calculate the constraint functions and gradients.
1404  bool calculate_constraints (Row_Vector& F,
1405  Matrix& Bx,
1406  Matrix& By) const;
1407 
1408 private:
1409  // The event we're fitting.
1411 
1412  // Vector of constraints.
1413  const vector<Constraint>& _constraints;
1414 
1415  // Argument values.
1417 
1418  // The pair table.
1420 };
1421 
1422 
1434  const vector<Constraint>& constraints,
1436 //
1437 // Purpose: Constructor.
1438 //
1439 // Inputs:
1440 // ev - The event we're fitting.
1441 // constraints - The list of constraints.
1442 // args - The parameter settings.
1443 //
1444 //
1445  : Constraint_Calculator (constraints.size() +
1446  ((ev.has_neutrino() || args.ignore_met()) ? 0 : 2)),
1447  _ev (ev),
1448  _constraints (constraints),
1449  _args (args),
1450  _pt (constraints, ev)
1451 
1452 {
1453 }
1454 
1455 
1470 bool
1472  Matrix& Bx,
1473  Matrix& By) const
1474 //
1475 // Purpose: Calculate the constraint functions and gradients.
1476 //
1477 // Outputs:
1478 // F - Vector of constraint values.
1479 // Bx - The well-measured variable gradients.
1480 // By - The poorly-measured variable gradients.
1481 //
1482 {
1483  // Clear the matrices.
1484  Bx = Matrix (Bx.num_row(), Bx.num_col(), 0);
1485  By = Matrix (By.num_row(), By.num_col(), 0);
1486  F = Row_Vector (F.num_col(), 0);
1487 
1488  const double p_eps = 1e-10;
1489 
1490  if (_ev.has_neutrino() && _ev.nu().z() > _args.e_com()) {
1491  return false;
1492  }
1493 
1494  int nobjs = _ev.nobjs ();
1495 
1496  // Reject the point if any of the momenta get too small.
1497  for (int j=0; j<nobjs; j++) {
1498  if (_ev.obj(j).p.perp() <= p_eps || _ev.obj(j).p.e() <= p_eps) {
1499  return false;
1500  }
1501  }
1502 
1503  if (! calculate_mass_constraints (_ev, _pt, _constraints, _args.use_e(),
1504  F, Bx, By))
1505  return false;
1506 
1507  if (_ev.has_neutrino())
1508  convert_neutrino (_ev, _args.use_e(), Bx, By);
1509  else if (!_args.ignore_met())
1510  {
1511  /* kt constraints */
1512  calculate_kt_constraints (_ev, _args.use_e(), F, Bx);
1513  }
1514 
1515  /* convert d/dtheta to d/deta */
1516  for (int j=0; j<nobjs; j++) {
1517  ddtheta_to_ddeta (obj_index (j) + eta_offs,
1518  _ev.obj(j).p.cosTheta(),
1519  Bx);
1520  }
1521 
1522  /* handle muons */
1523  for (int j=0; j<nobjs; j++) {
1524  const FE_Obj& obj = _ev.obj (j);
1525  if (obj.muon_p) {
1526  // Again, E vs. P doesn't matter here since we assume muons to be massless.
1527  double pmu2 = obj.p.vect().mag2();
1528  int ndx = obj_index (j) + p_offs;
1529  for (int k=1; k<=Bx.num_col(); k++)
1530  Bx(ndx, k) = - Bx(ndx, k) * pmu2;
1531  }
1532  }
1533 
1534  return true;
1535 }
1536 
1537 
1568  const Column_Vector& y,
1569  Row_Vector& F,
1570  Matrix& Bx,
1571  Matrix& By)
1572 //
1573 // Purpose: Evaluate constraints at the point described by X and Y
1574 // (well-measured and poorly-measured variables, respectively).
1575 // The results should be stored in F. BX and BY should be set
1576 // to the gradients of F with respect to X and Y, respectively.
1577 //
1578 // Inputs:
1579 // x - Vector of well-measured variables.
1580 // y - Vector of poorly-measured variables.
1581 //
1582 // Outputs:
1583 // F - Vector of constraint values.
1584 // Bx - The well-measured variable gradients.
1585 // By - The poorly-measured variable gradients.
1586 //
1587 {
1588  int nobjs = _ev.nobjs();
1589 
1590  const double p_eps = 1e-10;
1591  const double eta_max = 10;
1592 
1593  // Give up if we've gone into an obviously unphysical region.
1594  for (int j=0; j<nobjs; j++)
1595  if (x(obj_index (j) + p_offs) < p_eps ||
1596  fabs(x(obj_index (j) + eta_offs)) > eta_max) {
1597  return false;
1598  }
1599 
1600  unpack_event (_ev, x, y, _args.use_e(),
1601  _ev.has_neutrino() || !_args.ignore_met());
1602 
1603  return calculate_constraints (F, Bx, By);
1604 }
1605 
1606 
1607 //*************************************************************************
1608 // Mass calculation.
1609 //
1610 
1611 
1612 namespace {
1613 
1614 
1632 double calculate_sigm (const Matrix& Q,
1633  const Matrix& R,
1634  const Matrix& S,
1635  const Column_Vector& Bx,
1636  const Column_Vector& By)
1637 //
1638 // Purpose: Do error propagation to find the uncertainty in the final mass.
1639 //
1640 // Inputs:
1641 // Q - The final error matrix for the well-measured variables.
1642 // R - The final error matrix for the poorly-measured variables.
1643 // S - The final cross error matrix for the two sets of variables.
1644 // Bx - The well-measured variable gradients.
1645 // By - The poorly-measured variable gradients.
1646 //
1647 {
1648  double sig2 = scalar (Bx.T() * Q * Bx);
1649 
1650  if (By.num_row() > 0) {
1651  sig2 += scalar (By.T() * R * By);
1652  sig2 += 2 * scalar (Bx.T() * S * By);
1653  }
1654 
1655  assert (sig2 >= 0);
1656  return sqrt (sig2);
1657 }
1658 
1659 
1682 void calculate_mass (Fourvec_Event& ev,
1683  const vector<Constraint>& mass_constraint,
1684  const Fourvec_Constrainer_Args& args,
1685  double chisq,
1686  const Matrix& Q,
1687  const Matrix& R,
1688  const Matrix& S,
1689  double& m,
1690  double& sigm)
1691 //
1692 // Purpose: Calculate the final requested mass and its uncertainty.
1693 //
1694 // Inputs:
1695 // ev - The event we're fitting.
1696 // mass_constraint- The description of the mass we're to calculate.
1697 // args - Parameter settings.
1698 // chisq - The chisq from the fit.
1699 // Q - The final error matrix for the well-measured variables.
1700 // R - The final error matrix for the poorly-measured variables.
1701 // S - The final cross error matrix for the two sets of variables.
1702 //
1703 // Outputs:
1704 // m - The mass.
1705 // sigm - Its uncertainty.
1706 //
1707 {
1708  // Don't do anything if the mass wasn't specified.
1709  if (mass_constraint.size () == 0) {
1710  m = 0;
1711  sigm = 0;
1712  return;
1713  }
1714 
1715  // Do the constraint calculation.
1716  int n_measured_vars = ev.nobjs()*3 + 2;
1717  int n_unmeasured_vars = 0;
1718  if (ev.has_neutrino ()) n_unmeasured_vars = 1;
1719 
1720  Row_Vector F(1);
1721  Matrix Bx (n_measured_vars, 1);
1722  Matrix By (n_unmeasured_vars, 1);
1723 
1724  Fourvec_Constraint_Calculator cc (ev, mass_constraint, args);
1725  cc.calculate_constraints (F, Bx, By);
1726 
1727  // Calculate the mass.
1728  //assert (F(1) >= 0);
1729  if (F(1) >= 0.)
1730  m = sqrt (F(1) * 2);
1731  else {
1732  m = 0.;
1733  chisq = -100.;
1734  }
1735 
1736  // And the uncertainty.
1737  // We can only do this if the fit converged.
1738  if (chisq < 0)
1739  sigm = 0;
1740  else {
1741  //assert (F(1) > 0);
1742  Bx = Bx / (m);
1743  By = By / (m);
1744 
1745  sigm = calculate_sigm (Q, R, S, Bx, By);
1746  }
1747 }
1748 
1749 
1750 } // unnamed namespace
1751 
1752 
1754  double& m,
1755  double& sigm,
1756  Column_Vector& pullx,
1757  Column_Vector& pully)
1758 //
1759 // Purpose: Do a constrained fit for EV. Returns the requested mass and
1760 // its error in M and SIGM, and the pull quantities in PULLX and
1761 // PULLY. Returns the chisq; this will be < 0 if the fit failed
1762 // to converge.
1763 //
1764 // Inputs:
1765 // ev - The event we're fitting.
1766 //
1767 // Outputs:
1768 // ev - The fitted event.
1769 // m - Requested invariant mass.
1770 // sigm - Uncertainty on m.
1771 // pullx - Pull quantities for well-measured variables.
1772 // pully - Pull quantities for poorly-measured variables.
1773 //
1774 // Returns:
1775 // The fit chisq, or < 0 if the fit didn't converge.
1776 //
1777 {
1778  adjust_fourvecs (ev, _args.use_e ());
1779 
1780  bool use_kt = ev.has_neutrino() || !_args.ignore_met();
1781  int nobjs = ev.nobjs ();
1782  int n_measured_vars = nobjs * 3;
1783  int n_unmeasured_vars = 0;
1784  int n_constraints = _constraints.size ();
1785 
1786  if (use_kt) {
1787  n_measured_vars += 2;
1788 
1789  if (ev.has_neutrino ())
1790  n_unmeasured_vars = 1;
1791  else
1792  n_constraints += 2;
1793  }
1794 
1795  Matrix G_i (n_measured_vars, n_measured_vars);
1796  Diagonal_Matrix Y (n_unmeasured_vars);
1797  Column_Vector x (n_measured_vars);
1798  Column_Vector y (n_unmeasured_vars);
1799  pack_event (ev, _args.use_e(), use_kt, x, y, G_i, Y);
1800 
1801  Column_Vector xm = x;
1802  Column_Vector ym = y;
1803 
1804  // ??? Should allow for using a different underlying fitter here.
1806 
1808 
1809  Matrix Q;
1810  Matrix R;
1811  Matrix S;
1812  double chisq = fitter.fit (cc,
1813  xm, x, ym, y, G_i, Y,
1814  pullx, pully,
1815  Q, R, S);
1816 
1817  unpack_event (ev, x, y, _args.use_e (), use_kt);
1818 
1819  calculate_mass (ev, _mass_constraint, _args, chisq, Q, R, S, m, sigm);
1820 
1821  return chisq;
1822 }
1823 
1824 
1825 } // namespace hitfit
const Objpair & get_pair(std::vector< Objpair >::size_type pairno) const
Get one of the pairs from the table, index starts from 0.
Definition: Pair_Table.cc:100
friend std::ostream & operator<<(std::ostream &s, const Fourvec_Constrainer &c)
Output stream operator, print the content of this Fourvec_Constrainer to an output stream...
int for_constraint(std::vector< signed char >::size_type k) const
Retrieve the value set for a constraint.
Define an abstract interface for getting parameter settings.
Represent a mass constraint equation. Mass constraints come in two varieties, either saying that the ...
Definition: Constraint.h:79
void adjust_e_for_mass(Fourvec &v, double mass)
Adjust the energy component of four-vector v (leaving the three-vector part unchanged) so that the fo...
Definition: fourvec.cc:123
Represent an event for kinematic fitting as a collection of four-momenta. Each object is represented ...
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Fourvec_Constrainer(const Fourvec_Constrainer_Args &args)
Constructor.
CLHEP::HepVector Column_Vector
Definition: matutil.h:66
int npairs() const
Return the number of pairs in the table.
Definition: Pair_Table.cc:88
bool ev
Concrete realization of the Constraint_Calculator class. Evaluate constraints at the point described ...
Define matrix types for the HitFit package, and supply a few additional operations.
const FE_Obj & obj(std::vector< FE_Obj >::size_type i) const
Access object at index i, with the convention that the index starts at 0.
uint16_t size_type
int nobjs() const
Return the number of objects in the event not including any neutrinos.
Abstract base class for evaluating constraints. Users derive from this class and implement the eval()...
CLHEP::HepMatrix Matrix
Definition: matutil.h:65
Fourvec_Constrainer_Args(const Defaults &defs)
Constructor, initialize from a Defaults object.
T x() const
Cartesian x coordinate.
bool calculate_constraints(Row_Vector &F, Matrix &Bx, Matrix &By) const
Calculate the constraint functions and gradients.
double kt_y_error() const
Return the y uncertainty in .
double kt_xy_covar() const
Return the xy covariance in .
int j() const
Return the index of the second object.
T sqrt(T t)
Definition: SSEVec.h:18
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Minimize a subject to a set of constraints. Based on the SQUAW algorithm.
Minimize a subject to a set of constraints, based on the SQUAW algorithm.
CLHEP::HepLorentzVector Fourvec
Typedef for a HepLorentzVector.
Definition: fourvec.h:57
const Chisq_Constrainer_Args & chisq_constrainer_args() const
base
Make Sure CMSSW is Setup ##.
Hold on to parameters for the Chisq_Constrainer class.
const Fourvec_Constrainer_Args _args
double p2[4]
Definition: TauolaWrapper.h:90
void adjust_p_for_mass(Fourvec &v, double mass)
Adjust the three-vector part of v, leaving the energy unchanged,.
Definition: fourvec.cc:98
std::vector< Constraint > _mass_constraint
Fourvec_Constraint_Calculator(Fourvec_Event &ev, const vector< Constraint > &constraints, const Fourvec_Constrainer_Args &args)
Constructor.
int k[5][pyjets_maxn]
void set_nu_p(const Fourvec &p)
Set the neutrino four-momentum to . This method adds a neutrino if there wasn&#39;t already one...
void set_obj_p(std::vector< FE_Obj >::size_type i, const Fourvec &p)
Set the four-momentum of object at index i to .
Float e1
Definition: deltaR.h:20
CLHEP::HepDiagMatrix Diagonal_Matrix
Definition: matutil.h:67
void set_x_p(const Fourvec &p)
Set the four-momentum of the object.
const Fourvec & kt() const
Access the four-momentum.
A lookup table to speed up constraint evaluation using Fourvec_Constrainer.
Definition: Pair_Table.h:99
virtual bool eval(const Column_Vector &x, const Column_Vector &y, Row_Vector &F, Matrix &Bx, Matrix &By) override
Evaluate constraints at the point described by x and y (well-measured and poorly-measured variables...
double S(const TLorentzVector &, const TLorentzVector &)
Definition: Particle.cc:99
Float e2
Definition: deltaR.h:21
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
Chisq_Constrainer_Args _chisq_constrainer_args
std::vector< Constraint > _constraints
double p1[4]
Definition: TauolaWrapper.h:89
T dot(const Basic3DVector &v) const
Scalar product, or "dot" product, with a vector of same type.
double constrain(Fourvec_Event &ev, double &m, double &sigm, Column_Vector &pullx, Column_Vector &pully)
Do a constrained fit for event ev. Returns the requested mass and its uncertainty in m and sigm...
if(dp >Float(M_PI)) dp-
Define an interface for getting parameter settings.
Definition: Defaults.h:61
const Fourvec_Constrainer_Args & _args
Do a kinematic fit for a set of four-vectors, given a set of mass constraints.
void mass_constraint(std::string s)
Specify a combination of objects that will be returned by the constrain() method as mass...
double kt_x_error() const
Return the x uncertainty in .
Represent a single object in a Fourvec_Event, this is just a dumb data container. Each object in a Fo...
Definition: Fourvec_Event.h:94
static uInt32 F(BLOWFISH_CTX *ctx, uInt32 x)
Definition: blowfish.cc:281
const vector< Constraint > & _constraints
void add_constraint(std::string s)
Specify an additional constraint s for the problem. The format for s is described in the class descri...
const Fourvec & nu() const
Access the neutrino four-momentum.
Do a kinematic fit for a set of four-momenta, given a set of mass constraints.
A lookup table to speed up constraint evaluation.
Represent a pair of objects in Pair_Table.
Definition: Objpair.h:62
Row-vector class. CLHEP doesn&#39;t have a row-vector class, so HitFit uses its own. This is only a simpl...
Definition: matutil.h:83
Hold on to parameters for the Fourvec_Constrainer class.
int i() const
Return the index of the first object.
double scalar(const CLHEP::HepGenMatrix &m)
Return the matrix as a scalar. Raise an assertion if the matris is not .
Definition: matutil.cc:183
Represent an event for kinematic fitting as a collection of four-momenta.
bool has_neutrino() const
Return TRUE is this event contains a neutrino, otherwise returns FALSE.