CMS 3D CMS Logo

ParametrizationHelper.cc
Go to the documentation of this file.
2 #include <cmath>
3 #include <Math/CylindricalEta3D.h>
4 #include <Math/Polar3D.h>
5 #include <Math/Cartesian3D.h>
6 #include <Math/DisplacementVector3D.h>
7 #include <Math/Functions.h>
8 
9 using namespace std;
10 
11 //pat::helper::ParametrizationHelper::POLAR_ZERO();
12 //pat::helper::ParametrizationHelper::CARTESIAN_ZERO();
13 
15  using namespace std;
16 
17  typedef pat::CandKinResolution::Parametrization Parametrization;
18 
19  static map<string, Parametrization> const parMaps = {{"Cart", pat::CandKinResolution::Cart},
25  {"MCPInvSpher", pat::CandKinResolution::MCPInvSpher},
27  {"EtThetaPhi", pat::CandKinResolution::EtThetaPhi},
31  {"EScaledMomDev", pat::CandKinResolution::EScaledMomDev}};
32  map<string, Parametrization>::const_iterator itP = parMaps.find(name);
33  if (itP == parMaps.end()) {
34  throw cms::Exception("StringResolutionProvider") << "Bad parametrization '" << name.c_str() << "'";
35  }
36  return itP->second;
37 }
38 
40  switch (param) {
42  return "Cart";
44  return "ECart";
46  return "MCCart";
48  return "Spher";
50  return "ESpher";
52  return "MCSpher";
54  return "MCPInvSpher";
56  return "EtEtaPhi";
58  return "EtThetaPhi";
60  return "MomDev";
62  return "EMomDev";
64  return "MCMomDev";
66  return "EScaledMomDev";
68  return "Invalid";
69  default:
70  return "UNKNOWN";
71  }
72 }
73 
77  const math::PtEtaPhiMLorentzVector &initialP4) {
79  ROOT::Math::CylindricalEta3D<double> converter;
80  double m2;
81  switch (parametrization) {
82  // ======= CARTESIAN ==========
84  ret = math::XYZTLorentzVector(parameters[0],
85  parameters[1],
86  parameters[2],
87  sqrt(parameters[0] * parameters[0] + parameters[1] * parameters[1] +
88  parameters[2] * parameters[2] + parameters[3] * parameters[3]));
89  ret.SetM(parameters[3]); // to be sure about roundoffs
90  break;
92  ret = math::XYZTLorentzVector(parameters[0],
93  parameters[1],
94  parameters[2],
95  sqrt(parameters[0] * parameters[0] + parameters[1] * parameters[1] +
96  parameters[2] * parameters[2] + initialP4.mass() * initialP4.mass()));
97  ret.SetM(initialP4.mass()); // to be sure about roundoffs
98  break;
100  ret = math::XYZTLorentzVector(parameters[0], parameters[1], parameters[2], parameters[3]);
101  break;
102  // ======= SPHERICAL ==========
104  converter = ROOT::Math::Polar3D<double>(parameters[0], parameters[1], 0);
105  ret.SetCoordinates(converter.Rho(), converter.Eta(), parameters[2], parameters[3]);
106  break;
107  case pat::CandKinResolution::MCSpher: // same as above
108  converter = ROOT::Math::Polar3D<double>(parameters[0], parameters[1], 0);
109  ret.SetCoordinates(converter.Rho(), converter.Eta(), parameters[2], initialP4.mass());
110  break;
112  converter = ROOT::Math::Polar3D<double>(parameters[0], parameters[1], 0);
113  m2 = -parameters[0] * parameters[0] + parameters[3] * parameters[3];
114  ret.SetCoordinates(converter.Rho(), converter.Eta(), parameters[2], (m2 > 0 ? sqrt(m2) : 0.0));
115  break;
117  converter = ROOT::Math::Polar3D<double>(1.0 / parameters[0], parameters[1], 0);
118  ret.SetCoordinates(converter.Rho(), converter.Eta(), parameters[2], initialP4.mass());
119  break;
120  // ======= HEP CYLINDRICAL ==========
122  converter = ROOT::Math::Polar3D<double>(1.0, parameters[1], 0);
123  ret.SetCoordinates(parameters[0], converter.Eta(), parameters[2], 0);
124  break;
125  case pat::CandKinResolution::EtEtaPhi: // as simple as that
126  ret.SetCoordinates(parameters[0], parameters[1], parameters[2], 0);
127  break;
128  // ======= MomentumDeviates ==========
133  ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<double> > p = initialP4.Vect(), uz(0, 0, 1), uph, uth;
134  uph = uz.Cross(p).Unit();
135  uth = p.Cross(uph).Unit();
136  p *= parameters[0];
137  p += uth * parameters[1] + uph * parameters[2];
138  if (parametrization == pat::CandKinResolution::MomDev) {
139  ret.SetCoordinates(p.Rho(), p.Eta(), p.Phi(), initialP4.mass() * parameters[3]);
140  } else if (parametrization == pat::CandKinResolution::EMomDev) {
141  double m2 = ROOT::Math::Square(parameters[3] * initialP4.energy()) - p.Mag2();
142  ret.SetCoordinates(p.Rho(), p.Eta(), p.Phi(), (m2 > 0 ? sqrt(m2) : 0.0));
143  } else if (parametrization == pat::CandKinResolution::EScaledMomDev) {
144  double m2 = ROOT::Math::Square(p.R() * initialP4.E() / initialP4.P()) - p.Mag2();
145  ret.SetCoordinates(p.Rho(), p.Eta(), p.Phi(), (m2 > 0 ? sqrt(m2) : 0.0));
146  } else if (parametrization == pat::CandKinResolution::MCMomDev) {
147  ret.SetCoordinates(p.Rho(), p.Eta(), p.Phi(), initialP4.mass());
148  }
149  break;
150  }
151  // ======= OTHER ==========
153  throw cms::Exception("Invalid parametrization") << parametrization;
154  default:
155  throw cms::Exception("Not Implemented")
156  << "getResolEta not yet implemented for parametrization " << parametrization;
157  }
158  return ret;
159 }
160 
164  const math::XYZTLorentzVector &initialP4) {
166  switch (parametrization) {
167  // ======= CARTESIAN ==========
169  ret.SetCoordinates(parameters[0],
170  parameters[1],
171  parameters[2],
172  sqrt(parameters[0] * parameters[0] + parameters[1] * parameters[1] +
173  parameters[2] * parameters[2] + parameters[3] * parameters[3]));
174  break;
175  case pat::CandKinResolution::MCCart: // same as above
176  ret.SetCoordinates(parameters[0],
177  parameters[1],
178  parameters[2],
179  sqrt(parameters[0] * parameters[0] + parameters[1] * parameters[1] +
180  parameters[2] * parameters[2] + initialP4.mass() * initialP4.mass()));
181  break;
183  ret.SetCoordinates(parameters[0], parameters[1], parameters[2], parameters[3]);
184  break;
185  // ======= MomentumDeviates ==========
190  ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<double> > p = initialP4.Vect(), uz(0, 0, 1), uph, uth;
191  uph = uz.Cross(p).Unit();
192  uth = p.Cross(uph).Unit();
193  p *= parameters[0];
194  p += uth * parameters[1] + uph * parameters[2];
195  if (parametrization == pat::CandKinResolution::MomDev) {
196  ret.SetCoordinates(p.X(), p.Y(), p.Z(), sqrt(p.Mag2() + ROOT::Math::Square(initialP4.mass() * parameters[3])));
197  } else if (parametrization == pat::CandKinResolution::EMomDev) {
198  ret.SetCoordinates(p.X(), p.Y(), p.Z(), parameters[3] * initialP4.energy());
199  } else if (parametrization == pat::CandKinResolution::EMomDev) {
200  ret.SetCoordinates(p.X(), p.Y(), p.Z(), p.R() * initialP4.E() / initialP4.P());
201  } else {
202  ret.SetCoordinates(p.X(), p.Y(), p.Z(), sqrt(p.Mag2() + initialP4.mass() * initialP4.mass()));
203  }
204  break;
205  }
206  // ======= ALL OTHERS ==========
207  default:
208  ret = polarP4fromParameters(parametrization, parameters, math::PtEtaPhiMLorentzVector(initialP4));
209  }
210  return ret;
211 }
212 
213 template <typename T>
216  switch (parametrization) {
217  // ======= CARTESIAN ==========
219  ret[0] = p4.px();
220  ret[1] = p4.py();
221  ret[2] = p4.pz();
222  ret[3] = p4.mass();
223  break;
225  ret[0] = p4.px();
226  ret[1] = p4.py();
227  ret[2] = p4.pz();
228  ret[3] = p4.mass();
229  break;
231  ret[0] = p4.px();
232  ret[1] = p4.py();
233  ret[2] = p4.pz();
234  ret[3] = p4.energy();
235  break;
236  // ======= SPHERICAL ==========
238  ret[0] = p4.P();
239  ret[1] = p4.theta();
240  ret[2] = p4.phi();
241  ret[3] = p4.mass();
242  break;
244  ret[0] = p4.P();
245  ret[1] = p4.theta();
246  ret[2] = p4.phi();
247  ret[3] = p4.mass();
248  break;
250  ret[0] = p4.P();
251  ret[1] = p4.theta();
252  ret[2] = p4.phi();
253  ret[3] = p4.energy();
254  break;
256  ret[0] = 1.0 / p4.P();
257  ret[1] = p4.theta();
258  ret[2] = p4.phi();
259  ret[3] = p4.mass();
260  break;
261  // ======= HEP CYLINDRICAL ==========
263  ret[0] = p4.pt();
264  ret[1] = p4.theta();
265  ret[2] = p4.phi();
266  ret[3] = 0;
267  break;
269  ret[0] = p4.pt();
270  ret[1] = p4.eta();
271  ret[2] = p4.phi();
272  ret[3] = 0;
273  break;
274  // ======= DEVIATES ==========
277  ret[0] = 1.0;
278  ret[1] = 0.0;
279  ret[2] = 0.0;
280  ret[3] = 1.0;
281  break;
283  ret[0] = 1.0;
284  ret[1] = 0.0;
285  ret[2] = 0.0;
286  ret[3] = p4.mass();
287  break;
289  ret[0] = 1.0;
290  ret[1] = 0.0;
291  ret[2] = 0.0;
292  ret[3] = p4.E() / p4.P();
293  break;
294  // ======= OTHER ==========
296  throw cms::Exception("Invalid parametrization") << parametrization;
297  default:
298  throw cms::Exception("Not Implemented")
299  << "getResolEta not yet implemented for parametrization " << parametrization;
300  }
301 }
302 
306  setParametersFromAnyVector(parametrization, v, p4);
307 }
308 
311  const math::XYZTLorentzVector &p4) {
312  setParametersFromAnyVector(parametrization, v, p4);
313 }
314 
318  setParametersFromP4(parametrization, ret, p4);
319  return ret;
320 }
321 
325  setParametersFromP4(parametrization, ret, p4);
326  return ret;
327 }
328 
331  const math::PtEtaPhiMLorentzVector &p4ini,
332  const math::PtEtaPhiMLorentzVector &p4fin) {
334  switch (parametrization) {
338  ret = parametersFromP4(parametrization, p4fin) - parametersFromP4(parametrization, p4ini);
339  break;
346  ret = parametersFromP4(parametrization, p4fin) - parametersFromP4(parametrization, p4ini);
347  while (ret[2] > +M_PI)
348  ret[2] -= (2 * M_PI);
349  while (ret[2] < -M_PI)
350  ret[2] += (2 * M_PI);
351  break;
356  return diffToParameters(parametrization, math::XYZTLorentzVector(p4ini), math::XYZTLorentzVector(p4fin));
358  throw cms::Exception("Invalid parametrization") << parametrization;
359  default:
360  throw cms::Exception("Not Implemented")
361  << "diffToParameters not yet implemented for parametrization " << parametrization;
362  }
363  return ret;
364 }
365 
368  const math::XYZTLorentzVector &p4ini,
369  const math::XYZTLorentzVector &p4fin) {
371  switch (parametrization) {
376  ret = parametersFromP4(parametrization, p4fin) - parametersFromP4(parametrization, p4ini);
377  break;
383  ret = parametersFromP4(parametrization, p4fin) - parametersFromP4(parametrization, p4ini);
384  while (ret[2] > +M_PI)
385  ret[2] -= (2 * M_PI);
386  while (ret[2] < -M_PI)
387  ret[2] += (2 * M_PI);
388  break;
393  typedef ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<double> > V3Cart;
394  V3Cart p1 = p4ini.Vect(), p2 = p4fin.Vect();
395  V3Cart ur = p1.Unit();
396  V3Cart uz(0, 0, 1);
397  V3Cart uph = uz.Cross(ur).Unit();
398  V3Cart uth = ur.Cross(uph).Unit();
399  ret[0] = p2.Dot(ur) / p1.R() - 1.0;
400  ret[1] = (p2 - p1).Dot(uth);
401  ret[2] = (p2 - p1).Dot(uph);
402  if (parametrization == pat::CandKinResolution::MomDev) {
403  ret[3] = p4fin.mass() / p4ini.mass() - 1.0;
404  } else if (parametrization == pat::CandKinResolution::EMomDev) {
405  ret[3] = p4fin.energy() / p4ini.energy() - 1.0;
406  }
407  } break;
409  throw cms::Exception("Invalid parametrization") << parametrization;
410  default:
411  throw cms::Exception("Not Implemented")
412  << "diffToParameters not yet implemented for parametrization " << parametrization;
413  }
414  return ret;
415 }
416 
418  switch (parametrization) {
430  return false;
433  return true;
435  throw cms::Exception("Invalid parametrization") << parametrization;
436  default:
437  throw cms::Exception("Not Implemented")
438  << "isAlwaysMassless not yet implemented for parametrization " << parametrization;
439  }
440 }
441 
443  switch (parametrization) {
456  return false;
458  return true;
460  throw cms::Exception("Invalid parametrization") << parametrization;
461  default:
462  throw cms::Exception("Not Implemented")
463  << "isAlwaysMassless not yet implemented for parametrization " << parametrization;
464  }
465 }
466 
468  switch (parametrization) {
478  return false;
483  return true;
485  throw cms::Exception("Invalid parametrization") << parametrization;
486  default:
487  throw cms::Exception("Not Implemented")
488  << "isAlwaysMassless not yet implemented for parametrization " << parametrization;
489  }
490 }
491 
494  const math::PtEtaPhiMLorentzVector &initialP4) {
495  switch (parametrization) {
496  // ======= CARTESIAN ==========
498  return parameters[3] >= 0; // M >= 0
500  return true;
502  return (parameters[0] * parameters[0] + parameters[1] * parameters[1] + parameters[2] * parameters[2] <=
503  parameters[3] * parameters[3]); // E >= P
505  return (parameters[0] >= 0) && // P >= 0
506  (parameters[3] >= 0) && // M >= 0
507  (parameters[1] >= 0) && // theta >= 0
508  (parameters[1] <= M_PI); // theta <= pi
510  return (parameters[0] >= 0) && // P >= 0
511  (parameters[1] >= 0) && // theta >= 0
512  (parameters[1] <= M_PI); // theta <= pi
514  return (parameters[0] >= 0) && // P >= 0
515  (parameters[3] >= parameters[0]) && // E >= P
516  (parameters[1] >= 0) && // theta >= 0
517  (parameters[1] <= M_PI); // theta <= PI
519  return (parameters[0] > 0) && // 1/P > 0
520  (parameters[1] >= 0) && // theta >= 0
521  (parameters[1] <= M_PI); // theta <= pi
522  // ======= HEP CYLINDRICAL ==========
524  return (parameters[0] > 0) && // Et >= 0
525  (parameters[1] >= 0) && // theta >= 0
526  (parameters[1] <= M_PI); // theta <= pi
527  case pat::CandKinResolution::EtEtaPhi: // as simple as that
528  return (parameters[0] > 0); // Et >= 0
529  // ======= MomentumDeviates ==========
531  return (parameters[0] >= 0) && // P >= 0
532  (parameters[3] >= 0); // m/M0 >= 0
534  return (parameters[0] >= 0); // P >= 0
537  if (parameters[0] <= 0)
538  return false;
539  ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<double> > p = initialP4.Vect(), uz(0, 0, 1), uph, uth;
540  uph = uz.Cross(p).Unit();
541  uth = p.Cross(uph).Unit();
542  p *= parameters[0];
543  p += uth * parameters[1] + uph * parameters[2];
544  if (parametrization == pat::CandKinResolution::EMomDev) {
545  if (parameters[3] < 0)
546  return false;
547  double m2 = ROOT::Math::Square(parameters[3] * initialP4.energy()) - p.Mag2();
548  if (m2 < 0)
549  return false;
550  } else if (parametrization == pat::CandKinResolution::EScaledMomDev) {
551  if (parameters[3] < 0)
552  return false;
553  double m2 = ROOT::Math::Square(p.R() * initialP4.E() / initialP4.P()) - p.Mag2();
554  if (m2 < 0)
555  return false;
556  }
557  return true;
558  }
559  // ======= OTHER ==========
561  throw cms::Exception("Invalid parametrization") << parametrization;
562  default:
563  throw cms::Exception("Not Implemented")
564  << "getResolEta not yet implemented for parametrization " << parametrization;
565  }
566 }
AlgebraicVector4 parametersFromP4(pat::CandKinResolution::Parametrization parametrization, const math::XYZTLorentzVector &p4)
Returns a vector of coordinates values given a coordinate frame and a 4-vector.
const char * name(pat::CandKinResolution::Parametrization param)
Convert a number into a string.
ret
prodAgent to be discontinued
PtEtaPhiMLorentzVectorD PtEtaPhiMLorentzVector
Lorentz vector with cartesian internal representation.
Definition: LorentzVector.h:25
parametrization
specify parametrization (see SWGuidePATKinematicResolutions for more details)
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
T sqrt(T t)
Definition: SSEVec.h:19
double p4[4]
Definition: TauolaWrapper.h:92
void setParametersFromP4(pat::CandKinResolution::Parametrization parametrization, AlgebraicVector4 &pars, const math::XYZTLorentzVector &p4)
Set the values of the parameters for a given 4-momentum.
double p2[4]
Definition: TauolaWrapper.h:90
ROOT::Math::SVector< double, 4 > AlgebraicVector4
bool isAlwaysMassive(pat::CandKinResolution::Parametrization parametrization)
Is this parametrization usable only with massive objects?
#define M_PI
math::PtEtaPhiMLorentzVector polarP4fromParameters(pat::CandKinResolution::Parametrization parametrization, const AlgebraicVector4 &parameters, const math::PtEtaPhiMLorentzVector &initialP4)
AlgebraicVector4 diffToParameters(pat::CandKinResolution::Parametrization parametrization, const math::XYZTLorentzVector &p4ini, const math::XYZTLorentzVector &p4fin)
Expresses the difference between two 4-momentum vectors as a shift in coordinates in a given frame...
bool isPhysical(pat::CandKinResolution::Parametrization parametrization, const AlgebraicVector4 &v4, const math::PtEtaPhiMLorentzVector &initialP4)
math::XYZTLorentzVector p4fromParameters(pat::CandKinResolution::Parametrization parametrization, const AlgebraicVector4 &parameters, const math::XYZTLorentzVector &initialP4)
double p1[4]
Definition: TauolaWrapper.h:89
void setParametersFromAnyVector(pat::CandKinResolution::Parametrization parametrization, AlgebraicVector4 &pars, const T &p4)
For internal use only, so we provide only the interface. Use the &#39;setParametersFromP4&#39;.
pat::CandKinResolution::Parametrization fromString(const std::string &name)
Convert a name into a parametrization code.
long double T
bool isAlwaysMassless(pat::CandKinResolution::Parametrization parametrization)
Is this parametrization usable only with massless objects?
bool isMassConstrained(pat::CandKinResolution::Parametrization parametrization)
If this parametrization has a mass constraint (including the &#39;isAlwaysMassless&#39; case) ...