ATLAS Offline Software
FourMomentumError.icc
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "FourMom/ErrorMatrixBase.h"
6 #include "FourMom/ErrorMatrixPxPyPzE.h"
7 #include "FourMom/ErrorMatrixEEtaPhiM.h"
8 #include "FourMom/ErrorMatrixPtEtaPhiM.h"
9 #include "FourMom/ErrorMatrixPtCotThPhiM.h"
10 #include "FourMom/P4ErrorTransforms.h"
11 
12 
13 template < class FourMom>
14 FourMomentumError< FourMom>::FourMomentumError(const ErrorMatrixPxPyPzE& m,
15  const FourMom& p4) :
16  m_p4(p4),
17  m_PxPyPzE(std::make_unique<ErrorMatrixPxPyPzE>(m)),
18  m_EEtaPhiM(nullptr),
19  m_PtEtaPhiM(nullptr),
20  m_PtCotThPhiM(nullptr),
21  m_originalType(PxPyPzE)
22 {}
23 
24 template < class FourMom>
25 FourMomentumError< FourMom>::FourMomentumError(const ErrorMatrixEEtaPhiM& m,
26  const FourMom& p4) :
27  m_p4(p4),
28  m_PxPyPzE(nullptr),
29  m_EEtaPhiM(std::make_unique<ErrorMatrixEEtaPhiM>(m)),
30  m_PtEtaPhiM(nullptr),
31  m_PtCotThPhiM(nullptr),
32  m_originalType(EEtaPhiM)
33 {}
34 
35 template < class FourMom>
36 FourMomentumError< FourMom>::FourMomentumError(const ErrorMatrixPtEtaPhiM& m,
37  const FourMom& p4) :
38  m_p4(p4),
39  m_PxPyPzE(nullptr),
40  m_EEtaPhiM(nullptr),
41  m_PtEtaPhiM(std::make_unique<ErrorMatrixPtEtaPhiM>(m)),
42  m_PtCotThPhiM(nullptr),
43  m_originalType(PtEtaPhiM)
44 {}
45 
46 template < class FourMom>
47 FourMomentumError< FourMom>::FourMomentumError(const ErrorMatrixPtCotThPhiM& m,
48  const FourMom& p4) :
49  m_p4(p4),
50  m_PxPyPzE(nullptr),
51  m_EEtaPhiM(nullptr),
52  m_PtEtaPhiM(nullptr),
53  m_PtCotThPhiM(std::make_unique<ErrorMatrixPtCotThPhiM>(m)),
54  m_originalType(PtCotThPhiM)
55 {}
56 
57 template < class FourMom>
58 FourMomentumError< FourMom>::FourMomentumError( const CLHEP::HepSymMatrix& m,
59  ParamType t,
60  const FourMom& p4) :
61  m_p4(p4),
62  m_PxPyPzE(nullptr),
63  m_EEtaPhiM(nullptr),
64  m_PtEtaPhiM(nullptr),
65  m_PtCotThPhiM(nullptr)
66 {
67  if ( t == PxPyPzE) {
68  m_PxPyPzE.store (std::make_unique<ErrorMatrixPxPyPzE>( m));
69  m_originalType = PxPyPzE;
70  }
71  if ( t == EEtaPhiM) {
72  m_EEtaPhiM.store (std::make_unique<ErrorMatrixEEtaPhiM>( m));
73  m_originalType = EEtaPhiM;
74  }
75  if ( t == PtEtaPhiM) {
76  m_PtEtaPhiM.store (std::make_unique<ErrorMatrixPtEtaPhiM>( m));
77  m_originalType = PtEtaPhiM;
78  }
79  if ( t == PtCotThPhiM) {
80  m_PtCotThPhiM.store (std::make_unique<ErrorMatrixPtCotThPhiM>( m));
81  m_originalType = PtCotThPhiM;
82  }
83 }
84 
85 template < class FourMom>
86 FourMomentumError< FourMom>::FourMomentumError( const FourMomentumError& other) :
87  I4MomentumError(),
88  m_p4( other.m_p4),
89  m_PxPyPzE(nullptr),
90  m_EEtaPhiM(nullptr),
91  m_PtEtaPhiM(nullptr),
92  m_PtCotThPhiM(nullptr),
93  m_originalType( other.originalType())
94 {
95  if (other.originalType() == PxPyPzE) {
96  m_PxPyPzE.store (std::make_unique<ErrorMatrixPxPyPzE>( *other.pxPyPzEMatrix()));
97  }
98  else if (other.originalType() == EEtaPhiM) {
99  m_EEtaPhiM.store (std::make_unique<ErrorMatrixEEtaPhiM>(*other.eEtaPhiMMatrix()));
100  }
101  else if (other.originalType() == PtEtaPhiM) {
102  m_PtEtaPhiM.store (std::make_unique<ErrorMatrixPtEtaPhiM>( *other.ptEtaPhiMMatrix()));
103  }
104  else if (other.originalType() == PtCotThPhiM) {
105  m_PtCotThPhiM.store (std::make_unique<ErrorMatrixPtCotThPhiM>( *other.ptCotThPhiMMatrix()));
106  }
107 }
108 
109 
110 template < class FourMom>
111 FourMomentumError< FourMom> &FourMomentumError<FourMom>::operator=( const FourMomentumError& other)
112 {
113  if (this!=&other){
114  m_p4( other.m_p4);
115  m_PxPyPzE.store (std::make_unique<ErrorMatrixPxPyPzE>());
116  m_EEtaPhiM.store (std::make_unique<ErrorMatrixEEtaPhiM>());
117  m_PtEtaPhiM.store (std::make_unique<ErrorMatrixPtEtaPhiM>());
118  m_PtCotThPhiM.store (std::make_unique<ErrorMatrixPtCotThPhiM>());
119  m_originalType = other.originalType();
120 
121  if (other.originalType() == PxPyPzE) {
122  m_PxPyPzE.store (std::make_unique<ErrorMatrixPxPyPzE>( *other.pxPyPzEMatrix()));
123  }
124  else if (other.originalType() == EEtaPhiM) {
125  m_EEtaPhiM.store (std::make_unique<ErrorMatrixEEtaPhiM>(*other.eEtaPhiMMatrix()));
126  }
127  else if (other.originalType() == PtEtaPhiM) {
128  m_PtEtaPhiM.store (std::make_unique<ErrorMatrixPtEtaPhiM>( *other.ptEtaPhiMMatrix()));
129  }
130  else if (other.originalType() == PtCotThPhiM) {
131  m_PtCotThPhiM.store (std::make_unique<ErrorMatrixPtCotThPhiM>( *other.ptCotThPhiMMatrix()));
132  }
133  }
134  return *this;
135 }
136 
137 template < class FourMom>
138 const ErrorMatrixPxPyPzE* FourMomentumError< FourMom>::pxPyPzEMatrix() const
139 {
140  if (!m_PxPyPzE) computePxPyPzEMatrix();
141  return m_PxPyPzE.get();
142 }
143 
144 template < class FourMom>
145 const ErrorMatrixEEtaPhiM* FourMomentumError< FourMom>::eEtaPhiMMatrix() const
146 {
147  if (!m_EEtaPhiM) computeEEtaPhiMMatrix();
148  return m_EEtaPhiM.get();
149 }
150 
151 template < class FourMom>
152 const ErrorMatrixPtEtaPhiM* FourMomentumError< FourMom>::ptEtaPhiMMatrix() const
153 {
154  if (!m_PtEtaPhiM) computePtEtaPhiMMatrix();
155  return m_PtEtaPhiM.get();
156 }
157 
158 template < class FourMom>
159 const ErrorMatrixPtCotThPhiM* FourMomentumError< FourMom>::ptCotThPhiMMatrix() const
160 {
161  return m_PtCotThPhiM.get();
162  // FIXME not implemented yet
163 }
164 
165 template < class FourMom>
166 bool FourMomentumError< FourMom>::computePxPyPzEMatrix() const
167 {
168  using namespace P4ErrorTransforms;
169 
170  if (m_originalType == EEtaPhiM) {
171  m_PxPyPzE.set(toPxPyPzE( *m_EEtaPhiM, m_p4.e(), m_p4.eta(), m_p4.phi(), m_p4.m()));
172  return true;
173  }
174  else if (m_originalType == PtEtaPhiM) {
175  //FIXME m_PxPyPzE = m_helper->pxPyPzE( *m_PtEtaPhiM);
176  return false;
177  }
178  else if (m_originalType == PtCotThPhiM) {
179  return false; // FIXME
180  }
181  else return false; // should never happen
182 }
183 
184 template < class FourMom>
185 bool FourMomentumError< FourMom>::computeEEtaPhiMMatrix() const
186 {
187  using namespace P4ErrorTransforms;
188 
189  if (m_originalType == PxPyPzE) {
190  m_EEtaPhiM.set(toEEtaPhiM( *m_PxPyPzE, m_p4.px(), m_p4.py(), m_p4.pz(), m_p4.e()));
191  return true;
192  }
193  else if (m_originalType == PtEtaPhiM) {
194  //FIXME m_EEtaPhiM = m_helper->eEtaPhiM( *m_PtEtaPhiM);
195  return false;
196  }
197  else if (m_originalType == PtCotThPhiM) {
198  return false; // FIXME
199  }
200  else return false; // should never happen
201 }
202 
203 template < class FourMom>
204 bool FourMomentumError< FourMom>::computePtEtaPhiMMatrix() const
205 {
206  if (m_originalType == PxPyPzE) {
207  //FIXME m_PtEtaPhiM = m_helper->ptEtaPhiM( *m_PxPyPzE);
208  return false;
209  }
210  else if (m_originalType == EEtaPhiM) {
211  //FIXME m_PtEtaPhiM = m_helper->ptEtaPhiM( *m_EEtaPhiM);
212  return false;
213  }
214  else if (m_originalType == PtCotThPhiM) {
215  return false; // FIXME
216  }
217  else return false; // should never happen
218 }
219 
220 template < class FourMom>
221 double FourMomentumError< FourMom>::pxError() const
222 {
223  const ErrorMatrixPxPyPzE* em = pxPyPzEMatrix();
224  if (em != 0) {
225  return std::sqrt( (*em)(ErrorMatrixPxPyPzE::ipx));
226  }
227  else return 0;
228 }
229 
230 template < class FourMom>
231 double FourMomentumError< FourMom>::pyError() const
232 {
233  const ErrorMatrixPxPyPzE* em = pxPyPzEMatrix();
234  if (em != 0) {
235  return std::sqrt( (*em)(ErrorMatrixPxPyPzE::ipy));
236  }
237  else return 0;
238 }
239 
240 template < class FourMom>
241 double FourMomentumError< FourMom>::pzError() const
242 {
243  const ErrorMatrixPxPyPzE* em = pxPyPzEMatrix();
244  if (em != 0) {
245  return std::sqrt( (*em)(ErrorMatrixPxPyPzE::ipz));
246  }
247  else return 0;
248 }
249 
250 template < class FourMom>
251 double FourMomentumError< FourMom>::mError() const
252 {
253  const ErrorMatrixEEtaPhiM* em = eEtaPhiMMatrix();
254  if (em != 0) {
255  return std::sqrt( (*em)(ErrorMatrixEEtaPhiM::im));
256  }
257  else return 0;
258 }
259 
260 template < class FourMom>
261 double FourMomentumError< FourMom>::m2Error() const
262 {
263  return 2* m_p4.m() * mError();
264 }
265 
266 template < class FourMom>
267 double FourMomentumError< FourMom>::pError() const
268 {
269  if (m_originalType == PxPyPzE) {
270  const ErrorMatrixPxPyPzE* em = pxPyPzEMatrix();
271  if (em != 0) {
272  double pxtmp = m_p4.px(); // avoid calling virtual methods multiple times
273  double pytmp = m_p4.py();
274  double pztmp = m_p4.pz();
275  double sigma2 =
276  sqr(pxtmp)*(*em)(ErrorMatrixPxPyPzE::ipx) +
277  sqr(pytmp)*(*em)(ErrorMatrixPxPyPzE::ipy) +
278  sqr(pztmp)*(*em)(ErrorMatrixPxPyPzE::ipz) +
279  2*pxtmp*pytmp*(*em)(ErrorMatrixPxPyPzE::ipx,ErrorMatrixPxPyPzE::ipy) +
280  2*pxtmp*pztmp*(*em)(ErrorMatrixPxPyPzE::ipx,ErrorMatrixPxPyPzE::ipz) +
281  2*pytmp*pztmp*(*em)(ErrorMatrixPxPyPzE::ipy,ErrorMatrixPxPyPzE::ipz);
282  return std::sqrt(sigma2/m_p4.p2()); // divide by p2 globally
283  }
284  }
285  else {
286  const ErrorMatrixEEtaPhiM* em = eEtaPhiMMatrix();
287  if (em != 0) {
288  double E = m_p4.e(); // avoid calling virtual methods multiple times
289  double M = m_p4.m();
290  double sigma2 =
291  E*E*(*em)(ErrorMatrixEEtaPhiM::ie) +
292  M*M*(*em)(ErrorMatrixEEtaPhiM::im) -
293  E*M*(*em)(ErrorMatrixEEtaPhiM::ie,ErrorMatrixEEtaPhiM::im);
294  return std::sqrt(sigma2/m_p4.p2()); // divide by p2 globally
295  }
296  }
297  // FIXME: no optimized version for PtEtaPhiM
298  return 0; // in case no error matrix, should never happen
299 }
300 
301 template < class FourMom>
302 double FourMomentumError< FourMom>::p2Error() const
303 {
304  return 2* m_p4.p() * pError();
305 }
306 
307 template < class FourMom>
308 double FourMomentumError< FourMom>::eError() const
309 {
310  if (m_originalType == PxPyPzE) {
311  const ErrorMatrixPxPyPzE* em = pxPyPzEMatrix();
312  if (em != 0) {
313  return std::sqrt( (*em)(ErrorMatrixPxPyPzE::ie));
314  }
315  else return 0;
316  }
317  else {
318  const ErrorMatrixEEtaPhiM* em = eEtaPhiMMatrix();
319  if (em != 0) {
320  return std::sqrt( (*em)(ErrorMatrixEEtaPhiM::ie));
321  }
322  else return 0;
323  }
324 }
325 
326 template < class FourMom>
327 double FourMomentumError< FourMom>::etaError() const
328 {
329  const ErrorMatrixEEtaPhiM* em = eEtaPhiMMatrix();
330  if (em != 0) {
331  return std::sqrt( (*em)(ErrorMatrixEEtaPhiM::ieta));
332  }
333  else return 0;
334 }
335 
336 template < class FourMom>
337 double FourMomentumError< FourMom>::rapidityError() const
338 {
339  return etaError();
340 }
341 
342 template < class FourMom>
343 double FourMomentumError< FourMom>::phiError() const
344 {
345  if (m_originalType == PtEtaPhiM) {
346  const ErrorMatrixPtEtaPhiM* em = ptEtaPhiMMatrix();
347  if (em != 0) {
348  return std::sqrt( (*em)(ErrorMatrixPtEtaPhiM::iphi));
349  }
350  }
351  else {
352  const ErrorMatrixEEtaPhiM* em = eEtaPhiMMatrix();
353  if (em != 0) {
354  return std::sqrt( (*em)(ErrorMatrixEEtaPhiM::iphi));
355  }
356  }
357  return 0; // in case no error matrix, should never happen
358 }
359 
360 template < class FourMom>
361 double FourMomentumError< FourMom>::etError() const
362 {
363  return 0;// FIXME!
364 }
365 
366 template < class FourMom>
367 double FourMomentumError< FourMom>::ptError() const
368 {
369  const ErrorMatrixPxPyPzE* em = pxPyPzEMatrix();
370  if (em == 0) return 0;
371  else {
372  double pxtmp = m_p4.px(); // avoid calling virtual methods multiple times
373  double pytmp = m_p4.py();
374  double sigma2 =
375  sqr(pxtmp)*(*em)(ErrorMatrixPxPyPzE::ipx) +
376  sqr(pytmp)*(*em)(ErrorMatrixPxPyPzE::ipy) +
377  2*pxtmp*pytmp*(*em)(ErrorMatrixPxPyPzE::ipx,ErrorMatrixPxPyPzE::ipy);
378  return std::sqrt(sigma2)/m_p4.pt(); // divide by pt globally
379  }
380 }
381 
382 template < class FourMom>
383 double FourMomentumError< FourMom>::iPtError() const
384 {
385  double p2tmp = m_p4.p2();
386  return p2tmp == 0? 0 : ptError() / p2tmp;
387 }
388 
389 template < class FourMom>
390 double FourMomentumError< FourMom>::cosPhiError() const
391 {
392  const ErrorMatrixEEtaPhiM* em = eEtaPhiMMatrix();
393  if (em == 0) return 0;
394  else {
395  return -m_p4.sinPhi() * std::sqrt( (*em)( ErrorMatrixEEtaPhiM::iphi));
396  }
397  // FIXME: missing optimization for ptEtaPhiMMatrix
398 }
399 
400 template < class FourMom>
401 double FourMomentumError< FourMom>::sinPhiError() const
402 {
403  const ErrorMatrixEEtaPhiM* em = eEtaPhiMMatrix();
404  if (em == 0) return 0;
405  else {
406  return m_p4.cosPhi() * std::sqrt( (*em)( ErrorMatrixEEtaPhiM::iphi));
407  }
408  // FIXME: missing optimization for ptEtaPhiMMatrix
409 }
410 
411 template < class FourMom>
412 double FourMomentumError< FourMom>::cosThError() const
413 {
414  const ErrorMatrixEEtaPhiM* em = eEtaPhiMMatrix();
415  if (em == 0) return 0;
416  else {
417  double t = exp(-m_p4.eta()); // intermediate variable
418  // derivative d(cosTheta)/d(theta/2)
419  double deriv_t = -4*t/sqr(1+t*t);
420  return std::abs( deriv_t * t) * std::sqrt( (*em)( ErrorMatrixEEtaPhiM::ieta));
421  }
422  // FIXME: missing optimization for ptEtaPhiMMatrix
423 }
424 
425 template < class FourMom>
426 double FourMomentumError< FourMom>::sinThError() const
427 {
428  const ErrorMatrixEEtaPhiM* em = eEtaPhiMMatrix();
429  if (em == 0) return 0;
430  else {
431  double t = exp(-m_p4.eta()); // intermediate variable
432  // derivative d(cosTheta)/d(theta/2)
433  double deriv_t = 2*(1-t*t)/sqr(1+t*t);
434  return std::abs( deriv_t * t) * std::sqrt( (*em)( ErrorMatrixEEtaPhiM::ieta));
435  }
436  // FIXME: missing optimization for ptEtaPhiMMatrix
437 }
438 
439 template < class FourMom>
440 double FourMomentumError< FourMom>::tanThError() const
441 {
442  const ErrorMatrixEEtaPhiM* em = eEtaPhiMMatrix();
443  if (em == 0) return 0;
444  else {
445  double t = exp(-m_p4.eta()); // intermediate variable
446  // derivative d(cosTheta)/d(theta/2)
447  double deriv_t = 2*(1+t*t)/sqr(1-t*t);
448  return std::abs( deriv_t * t) * std::sqrt( (*em)( ErrorMatrixEEtaPhiM::ieta));
449  }
450  // FIXME: missing optimization for ptEtaPhiMMatrix
451 }
452 
453 template < class FourMom>
454 double FourMomentumError< FourMom>::cotThError() const
455 {
456  const ErrorMatrixEEtaPhiM* em = eEtaPhiMMatrix();
457  if (em == 0) return 0;
458  else {
459  double t = exp(-m_p4.eta()); // intermediate variable
460  // derivative d(cosTheta)/d(theta/2)
461  double deriv_t = -(1+t*t)/(2*t*t);
462  return std::abs( deriv_t * t) * std::sqrt( (*em)( ErrorMatrixEEtaPhiM::ieta));
463  }
464  // FIXME: missing optimization for ptEtaPhiMMatrix
465 }
Generated on Sun Dec 22 2024 21:10:34 for ATLAS Offline Software by   doxygen 1.8.18