CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 
16  using namespace std;
17 
18  typedef pat::CandKinResolution::Parametrization Parametrization;
19 
20  static map<string,Parametrization> const parMaps = {
22  {"ECart" , pat::CandKinResolution::ECart},
23  {"MCCart" , pat::CandKinResolution::MCCart},
24  {"Spher" , pat::CandKinResolution::Spher},
25  {"ESpher" , pat::CandKinResolution::ESpher},
26  {"MCSpher" , pat::CandKinResolution::MCSpher},
27  {"MCPInvSpher" , pat::CandKinResolution::MCPInvSpher},
28  {"EtEtaPhi" , pat::CandKinResolution::EtEtaPhi},
29  {"EtThetaPhi" , pat::CandKinResolution::EtThetaPhi},
30  {"MomDev" , pat::CandKinResolution::MomDev},
31  {"EMomDev" , pat::CandKinResolution::EMomDev},
32  {"MCMomDev" , pat::CandKinResolution::MCMomDev},
33  {"EScaledMomDev" , pat::CandKinResolution::EScaledMomDev}
34  };
35  map<string,Parametrization>::const_iterator itP = parMaps.find(name);
36  if (itP == parMaps.end()) {
37  throw cms::Exception("StringResolutionProvider") << "Bad parametrization '" << name.c_str() << "'";
38  }
39  return itP->second;
40 }
41 
42 const char *
44  switch (param) {
45  case pat::CandKinResolution::Cart: return "Cart";
46  case pat::CandKinResolution::ECart: return "ECart";
47  case pat::CandKinResolution::MCCart: return "MCCart";
48  case pat::CandKinResolution::Spher: return "Spher";
49  case pat::CandKinResolution::ESpher: return "ESpher";
50  case pat::CandKinResolution::MCSpher: return "MCSpher";
51  case pat::CandKinResolution::MCPInvSpher: return "MCPInvSpher";
52  case pat::CandKinResolution::EtEtaPhi: return "EtEtaPhi";
53  case pat::CandKinResolution::EtThetaPhi: return "EtThetaPhi";
54  case pat::CandKinResolution::MomDev: return "MomDev";
55  case pat::CandKinResolution::EMomDev: return "EMomDev";
56  case pat::CandKinResolution::MCMomDev: return "MCMomDev";
57  case pat::CandKinResolution::EScaledMomDev: return "EScaledMomDev";
58  case pat::CandKinResolution::Invalid: return "Invalid";
59  default: return "UNKNOWN";
60  }
61 }
62 
66  const math::PtEtaPhiMLorentzVector &initialP4) {
68  ROOT::Math::CylindricalEta3D<double> converter;
69  double m2;
70  switch (parametrization) {
71  // ======= CARTESIAN ==========
73  ret = math::XYZTLorentzVector(parameters[0], parameters[1], parameters[2],
74  sqrt(parameters[0]*parameters[0] +
75  parameters[1]*parameters[1] +
76  parameters[2]*parameters[2] +
77  parameters[3]*parameters[3]) );
78  ret.SetM(parameters[3]); // to be sure about roundoffs
79  break;
81  ret = math::XYZTLorentzVector(parameters[0], parameters[1], parameters[2],
82  sqrt(parameters[0]*parameters[0] +
83  parameters[1]*parameters[1] +
84  parameters[2]*parameters[2] +
85  initialP4.mass()*initialP4.mass()) );
86  ret.SetM(initialP4.mass()); // to be sure about roundoffs
87  break;
89  ret = math::XYZTLorentzVector(parameters[0], parameters[1], parameters[2], parameters[3]);
90  break;
91  // ======= SPHERICAL ==========
93  converter = ROOT::Math::Polar3D<double>(parameters[0], parameters[1], 0);
94  ret.SetCoordinates(converter.Rho(),converter.Eta(),parameters[2],parameters[3]);
95  break;
96  case pat::CandKinResolution::MCSpher: // same as above
97  converter = ROOT::Math::Polar3D<double>(parameters[0], parameters[1], 0);
98  ret.SetCoordinates(converter.Rho(),converter.Eta(),parameters[2],initialP4.mass());
99  break;
101  converter = ROOT::Math::Polar3D<double>(parameters[0], parameters[1], 0);
102  m2 = - parameters[0]*parameters[0] + parameters[3]*parameters[3];
103  ret.SetCoordinates(converter.Rho(),converter.Eta(),parameters[2],(m2 > 0 ? sqrt(m2) : 0.0));
104  break;
106  converter = ROOT::Math::Polar3D<double>(1.0/parameters[0], parameters[1], 0);
107  ret.SetCoordinates(converter.Rho(),converter.Eta(),parameters[2],initialP4.mass());
108  break;
109  // ======= HEP CYLINDRICAL ==========
111  converter = ROOT::Math::Polar3D<double>(1.0, parameters[1], 0);
112  ret.SetCoordinates(parameters[0],converter.Eta(),parameters[2],0);
113  break;
114  case pat::CandKinResolution::EtEtaPhi: // as simple as that
115  ret.SetCoordinates(parameters[0], parameters[1], parameters[2], 0);
116  break;
117  // ======= MomentumDeviates ==========
122  {
123  ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D<double> > p = initialP4.Vect(), uz(0,0,1), uph, uth;
124  uph = uz.Cross(p).Unit();
125  uth = p.Cross(uph).Unit();
126  p *= parameters[0];
127  p += uth * parameters[1] + uph * parameters[2];
128  if (parametrization == pat::CandKinResolution::MomDev) {
129  ret.SetCoordinates(p.Rho(), p.Eta(), p.Phi(), initialP4.mass() * parameters[3]);
130  } else if (parametrization == pat::CandKinResolution::EMomDev) {
131  double m2 = ROOT::Math::Square(parameters[3] * initialP4.energy()) - p.Mag2();
132  ret.SetCoordinates(p.Rho(), p.Eta(), p.Phi(), (m2 > 0 ? sqrt(m2) : 0.0));
133  } else if (parametrization == pat::CandKinResolution::EScaledMomDev) {
134  double m2 = ROOT::Math::Square(p.R()*initialP4.E()/initialP4.P()) - p.Mag2();
135  ret.SetCoordinates(p.Rho(), p.Eta(), p.Phi(), (m2 > 0 ? sqrt(m2) : 0.0));
136  } else if (parametrization == pat::CandKinResolution::MCMomDev) {
137  ret.SetCoordinates(p.Rho(), p.Eta(), p.Phi(), initialP4.mass());
138  }
139  break;
140  }
141  // ======= OTHER ==========
143  throw cms::Exception("Invalid parametrization") << parametrization;
144  default:
145  throw cms::Exception("Not Implemented") << "getResolEta not yet implemented for parametrization " << parametrization ;
146  }
147  return ret;
148 }
149 
153  const math::XYZTLorentzVector &initialP4) {
155  switch (parametrization) {
156  // ======= CARTESIAN ==========
158  ret.SetCoordinates(parameters[0], parameters[1], parameters[2],
159  sqrt(parameters[0]*parameters[0] +
160  parameters[1]*parameters[1] +
161  parameters[2]*parameters[2] +
162  parameters[3]*parameters[3]) );
163  break;
164  case pat::CandKinResolution::MCCart: // same as above
165  ret.SetCoordinates(parameters[0], parameters[1], parameters[2],
166  sqrt(parameters[0]*parameters[0] +
167  parameters[1]*parameters[1] +
168  parameters[2]*parameters[2] +
169  initialP4.mass()*initialP4.mass()) );
170  break;
172  ret.SetCoordinates(parameters[0], parameters[1], parameters[2], parameters[3]);
173  break;
174  // ======= MomentumDeviates ==========
179  {
180  ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D<double> > p = initialP4.Vect(), uz(0,0,1), uph, uth;
181  uph = uz.Cross(p).Unit();
182  uth = p.Cross(uph).Unit();
183  p *= parameters[0];
184  p += uth * parameters[1] + uph * parameters[2];
185  if (parametrization == pat::CandKinResolution::MomDev) {
186  ret.SetCoordinates(p.X(), p.Y(), p.Z(), sqrt(p.Mag2() + ROOT::Math::Square(initialP4.mass() * parameters[3])) );
187  } else if (parametrization == pat::CandKinResolution::EMomDev) {
188  ret.SetCoordinates(p.X(), p.Y(), p.Z(), parameters[3] * initialP4.energy());
189  } else if (parametrization == pat::CandKinResolution::EMomDev) {
190  ret.SetCoordinates(p.X(), p.Y(), p.Z(), p.R() * initialP4.E()/initialP4.P());
191  } else {
192  ret.SetCoordinates(p.X(), p.Y(), p.Z(), sqrt(p.Mag2() + initialP4.mass()*initialP4.mass()));
193  }
194  break;
195  }
196  // ======= ALL OTHERS ==========
197  default:
198  ret = polarP4fromParameters(parametrization, parameters, math::PtEtaPhiMLorentzVector(initialP4));
199  }
200  return ret;
201 }
202 
203 template <typename T>
204 void
207  const T &p4) {
208  switch (parametrization) {
209  // ======= CARTESIAN ==========
211  ret[0] = p4.px(); ret[1] = p4.py(); ret[2] = p4.pz(); ret[3] = p4.mass();
212  break;
214  ret[0] = p4.px(); ret[1] = p4.py(); ret[2] = p4.pz(); ret[3] = p4.mass();
215  break;
217  ret[0] = p4.px(); ret[1] = p4.py(); ret[2] = p4.pz(); ret[3] = p4.energy();
218  break;
219  // ======= SPHERICAL ==========
221  ret[0] = p4.P(); ret[1] = p4.theta(); ret[2] = p4.phi(); ret[3] = p4.mass();
222  break;
224  ret[0] = p4.P(); ret[1] = p4.theta(); ret[2] = p4.phi(); ret[3] = p4.mass();
225  break;
227  ret[0] = p4.P(); ret[1] = p4.theta(); ret[2] = p4.phi(); ret[3] = p4.energy();
228  break;
230  ret[0] = 1.0/p4.P(); ret[1] = p4.theta(); ret[2] = p4.phi(); ret[3] = p4.mass();
231  break;
232  // ======= HEP CYLINDRICAL ==========
234  ret[0] = p4.pt(); ret[1] = p4.theta(); ret[2] = p4.phi(); ret[3] = 0;
235  break;
237  ret[0] = p4.pt(); ret[1] = p4.eta(); ret[2] = p4.phi(); ret[3] = 0;
238  break;
239  // ======= DEVIATES ==========
242  ret[0] = 1.0; ret[1] = 0.0; ret[2] = 0.0; ret[3] = 1.0;
243  break;
245  ret[0] = 1.0; ret[1] = 0.0; ret[2] = 0.0; ret[3] = p4.mass();
246  break;
248  ret[0] = 1.0; ret[1] = 0.0; ret[2] = 0.0; ret[3] = p4.E()/p4.P();
249  break;
250  // ======= OTHER ==========
252  throw cms::Exception("Invalid parametrization") << parametrization;
253  default:
254  throw cms::Exception("Not Implemented") << "getResolEta not yet implemented for parametrization " << parametrization ;
255  }
256 }
257 
258 void
261  setParametersFromAnyVector(parametrization, v, p4);
262 }
263 
264 void
267  setParametersFromAnyVector(parametrization, v, p4);
268 }
269 
273  setParametersFromP4(parametrization, ret, p4);
274  return ret;
275 }
276 
280  setParametersFromP4(parametrization, ret, p4);
281  return ret;
282 }
283 
287 {
289  switch (parametrization) {
293  ret = parametersFromP4(parametrization,p4fin) - parametersFromP4(parametrization,p4ini);
294  break;
301  ret = parametersFromP4(parametrization,p4fin) - parametersFromP4(parametrization,p4ini);
302  while(ret[2] > +M_PI) ret[2] -= (2*M_PI);
303  while(ret[2] < -M_PI) ret[2] += (2*M_PI);
304  break;
309  return diffToParameters(parametrization,
311  math::XYZTLorentzVector(p4fin));
313  throw cms::Exception("Invalid parametrization") << parametrization;
314  default:
315  throw cms::Exception("Not Implemented") << "diffToParameters not yet implemented for parametrization " << parametrization ;
316  }
317  return ret;
318 }
319 
320 
323  const math::XYZTLorentzVector &p4ini, const math::XYZTLorentzVector &p4fin)
324 {
326  switch (parametrization) {
331  ret = parametersFromP4(parametrization,p4fin) - parametersFromP4(parametrization,p4ini);
332  break;
338  ret = parametersFromP4(parametrization,p4fin) - parametersFromP4(parametrization,p4ini);
339  while(ret[2] > +M_PI) ret[2] -= (2*M_PI);
340  while(ret[2] < -M_PI) ret[2] += (2*M_PI);
341  break;
346  {
347  typedef ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D<double> > V3Cart;
348  V3Cart p1 = p4ini.Vect(), p2 = p4fin.Vect();
349  V3Cart ur = p1.Unit();
350  V3Cart uz(0,0,1);
351  V3Cart uph = uz.Cross(ur).Unit();
352  V3Cart uth = ur.Cross(uph).Unit();
353  ret[0] = p2.Dot(ur)/p1.R() - 1.0;
354  ret[1] = (p2 - p1).Dot(uth);
355  ret[2] = (p2 - p1).Dot(uph);
356  if (parametrization == pat::CandKinResolution::MomDev) {
357  ret[3] = p4fin.mass()/p4ini.mass() - 1.0;
358  } else if (parametrization == pat::CandKinResolution::EMomDev) {
359  ret[3] = p4fin.energy()/p4ini.energy() - 1.0;
360  }
361  }
362  break;
364  throw cms::Exception("Invalid parametrization") << parametrization;
365  default:
366  throw cms::Exception("Not Implemented") << "diffToParameters not yet implemented for parametrization " << parametrization ;
367  }
368  return ret;
369 }
370 
371 bool
373 {
374  switch (parametrization) {
386  return false;
389  return true;
391  throw cms::Exception("Invalid parametrization") << parametrization;
392  default:
393  throw cms::Exception("Not Implemented") << "isAlwaysMassless not yet implemented for parametrization " << parametrization ;
394  }
395 }
396 
397 bool
399 {
400  switch (parametrization) {
413  return false;
415  return true;
417  throw cms::Exception("Invalid parametrization") << parametrization;
418  default:
419  throw cms::Exception("Not Implemented") << "isAlwaysMassless not yet implemented for parametrization " << parametrization ;
420  }
421 }
422 
423 bool
425 {
426  switch (parametrization) {
436  return false;
441  return true;
443  throw cms::Exception("Invalid parametrization") << parametrization;
444  default:
445  throw cms::Exception("Not Implemented") << "isAlwaysMassless not yet implemented for parametrization " << parametrization ;
446  }
447 }
448 
449 bool
452  const math::PtEtaPhiMLorentzVector &initialP4) {
453  switch (parametrization) {
454  // ======= CARTESIAN ==========
456  return parameters[3] >= 0; // M >= 0
458  return true;
460  return (parameters[0]*parameters[0] + parameters[1]*parameters[1] + parameters[2]*parameters[2] <= parameters[3]*parameters[3]); // E >= P
462  return (parameters[0] >= 0 ) && // P >= 0
463  (parameters[3] >= 0 ) && // M >= 0
464  (parameters[1] >= 0 ) && // theta >= 0
465  (parameters[1] <= M_PI); // theta <= pi
467  return (parameters[0] >= 0 ) && // P >= 0
468  (parameters[1] >= 0 ) && // theta >= 0
469  (parameters[1] <= M_PI); // theta <= pi
471  return (parameters[0] >= 0 ) && // P >= 0
472  (parameters[3] >= parameters[0]) && // E >= P
473  (parameters[1] >= 0 ) && // theta >= 0
474  (parameters[1] <= M_PI ) ; // theta <= PI
476  return (parameters[0] > 0 ) && // 1/P > 0
477  (parameters[1] >= 0 ) && // theta >= 0
478  (parameters[1] <= M_PI); // theta <= pi
479  // ======= HEP CYLINDRICAL ==========
481  return (parameters[0] > 0 ) && // Et >= 0
482  (parameters[1] >= 0 ) && // theta >= 0
483  (parameters[1] <= M_PI); // theta <= pi
484  case pat::CandKinResolution::EtEtaPhi: // as simple as that
485  return (parameters[0] > 0); // Et >= 0
486  // ======= MomentumDeviates ==========
488  return (parameters[0] >= 0) && // P >= 0
489  (parameters[3] >= 0); // m/M0 >= 0
491  return (parameters[0] >= 0); // P >= 0
494  {
495  if (parameters[0] <= 0) return false;
496  ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D<double> > p = initialP4.Vect(), uz(0,0,1), uph, uth;
497  uph = uz.Cross(p).Unit();
498  uth = p.Cross(uph).Unit();
499  p *= parameters[0];
500  p += uth * parameters[1] + uph * parameters[2];
501  if (parametrization == pat::CandKinResolution::EMomDev) {
502  if (parameters[3] < 0) return false;
503  double m2 = ROOT::Math::Square(parameters[3] * initialP4.energy()) - p.Mag2();
504  if (m2 < 0) return false;
505  } else if (parametrization == pat::CandKinResolution::EScaledMomDev) {
506  if (parameters[3] < 0) return false;
507  double m2 = ROOT::Math::Square(p.R()*initialP4.E()/initialP4.P()) - p.Mag2();
508  if (m2 < 0) return false;
509  }
510  return true;
511  }
512  // ======= OTHER ==========
514  throw cms::Exception("Invalid parametrization") << parametrization;
515  default:
516  throw cms::Exception("Not Implemented") << "getResolEta not yet implemented for parametrization " << parametrization ;
517  }
518 }
519 
520 
tuple ret
prodAgent to be discontinued
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.
PtEtaPhiMLorentzVectorD PtEtaPhiMLorentzVector
Lorentz vector with cartesian internal representation.
Definition: LorentzVector.h:25
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
T sqrt(T t)
Definition: SSEVec.h:18
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
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.
ROOT::Math::SVector< double, 4 > AlgebraicVector4
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) ...