ATLAS Offline Software
Loading...
Searching...
No Matches
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
13template < class FourMom>
14FourMomentumError< 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
24template < class FourMom>
25FourMomentumError< 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
35template < class FourMom>
36FourMomentumError< 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
46template < class FourMom>
47FourMomentumError< 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
57template < class FourMom>
58FourMomentumError< 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
85template < class FourMom>
86FourMomentumError< 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
110template < class FourMom>
111FourMomentumError< 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
137template < class FourMom>
138const ErrorMatrixPxPyPzE* FourMomentumError< FourMom>::pxPyPzEMatrix() const
139{
140 if (!m_PxPyPzE) computePxPyPzEMatrix();
141 return m_PxPyPzE.get();
142}
143
144template < class FourMom>
145const ErrorMatrixEEtaPhiM* FourMomentumError< FourMom>::eEtaPhiMMatrix() const
146{
147 if (!m_EEtaPhiM) computeEEtaPhiMMatrix();
148 return m_EEtaPhiM.get();
149}
150
151template < class FourMom>
152const ErrorMatrixPtEtaPhiM* FourMomentumError< FourMom>::ptEtaPhiMMatrix() const
153{
154 if (!m_PtEtaPhiM) computePtEtaPhiMMatrix();
155 return m_PtEtaPhiM.get();
156}
157
158template < class FourMom>
159const ErrorMatrixPtCotThPhiM* FourMomentumError< FourMom>::ptCotThPhiMMatrix() const
160{
161 return m_PtCotThPhiM.get();
162 // FIXME not implemented yet
163}
164
165template < class FourMom>
166bool 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
184template < class FourMom>
185bool 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
203template < class FourMom>
204bool 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
220template < class FourMom>
221double 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
230template < class FourMom>
231double 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
240template < class FourMom>
241double 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
250template < class FourMom>
251double 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
260template < class FourMom>
261double FourMomentumError< FourMom>::m2Error() const
262{
263 return 2* m_p4.m() * mError();
264}
265
266template < class FourMom>
267double 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
301template < class FourMom>
302double FourMomentumError< FourMom>::p2Error() const
303{
304 return 2* m_p4.p() * pError();
305}
306
307template < class FourMom>
308double 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
326template < class FourMom>
327double 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
336template < class FourMom>
337double FourMomentumError< FourMom>::rapidityError() const
338{
339 return etaError();
340}
341
342template < class FourMom>
343double 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
360template < class FourMom>
361double FourMomentumError< FourMom>::etError() const
362{
363 return 0;// FIXME!
364}
365
366template < class FourMom>
367double 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
382template < class FourMom>
383double FourMomentumError< FourMom>::iPtError() const
384{
385 double p2tmp = m_p4.p2();
386 return p2tmp == 0? 0 : ptError() / p2tmp;
387}
388
389template < class FourMom>
390double 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
400template < class FourMom>
401double 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
411template < class FourMom>
412double 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
425template < class FourMom>
426double 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
439template < class FourMom>
440double 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
453template < class FourMom>
454double 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 for ATLAS Offline Software by doxygen 1.14.0