ATLAS Offline Software
P4Helpers.h
Go to the documentation of this file.
1 
3 /*
4  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
5 */
6 
7 // P4Helpers.h
8 // Header file for P4 helper functions
9 // Author: S.Binet<binet@cern.ch>
11 #ifndef FOURMOMUTILS_P4HELPERS_H
12 #define FOURMOMUTILS_P4HELPERS_H
13 
23 #include <cmath>
24 
25 namespace P4Helpers
26 {
28  inline
29  double deltaPhi( double phiA, double phiB )
30  {
31  return -remainder( -phiA + phiB, 2*M_PI );
32  }
33 }
34 
35 // AthAnalysisBase/ManaCore doesn't currently include the Trigger Service
36 #ifndef XAOD_ANALYSIS
37 
38 
39 // STL includes
40 #include <algorithm> // for std::sort
41 #include <limits> // for std::numeric_limits
42 
43 // CxxUtils includes
44 #include "CxxUtils/fpcompare.h" // for fpcompare::less
45 
46 // EventKernel includes
47 #include "EventKernel/I4Momentum.h"
48 
49 namespace P4Helpers
50 {
52  inline
53  double deltaEta( const I4Momentum& p1, const I4Momentum& p2 )
54  {
55  double dEta(999);
56  if ( p1.kind() == I4Momentum::P4PXPYPZE && p2.kind() == I4Momentum::P4PXPYPZE )
57  {
58  const double mom1 = p1.p(); const double pz1 = p1.pz();
59  const double mom2 = p2.p(); const double pz2 = p2.pz();
60  if ( ( mom1 + pz1 ) * ( mom2 - pz2 ) > 0 )
61  dEta = -0.5 * std::log( ( mom1 - pz1 ) / ( mom1 + pz1 )
62  * ( mom2 + pz2 ) / ( mom2 - pz2 ) );
63  }
64  else if ( p1.kind() == I4Momentum::P4PXPYPZE )
65  {
66  const double mom1 = p1.p(); const double pz1 = p1.pz();
67  if ( mom1 + pz1 > 0 )
68  dEta = -0.5 * std::log( ( mom1 - pz1 ) / ( mom1 + pz1 ) ) - p2.eta();
69  }
70  else if ( p2.kind() == I4Momentum::P4PXPYPZE )
71  {
72  const double mom2 = p2.p(); const double pz2 = p2.pz();
73  if ( mom2 - pz2 > 0 )
74  dEta = p1.eta() - 0.5 * std::log( ( mom2 + pz2 ) / ( mom2 - pz2 ) );
75  }
76  else
77  {
78  dEta = p1.eta() - p2.eta();
79  }
80  return dEta;
81  }
82 
84  inline
85  double deltaEtaSq( const I4Momentum& p1, const I4Momentum& p2 )
86  {
87  double dEtaSq(999);
88  if ( p1.kind() == I4Momentum::P4PXPYPZE && p2.kind() == I4Momentum::P4PXPYPZE )
89  {
90  const double mom1 = p1.p(); const double pz1 = p1.pz();
91  const double mom2 = p2.p(); const double pz2 = p2.pz();
92  if ( ( mom1 + pz1 ) * ( mom2 - pz2 ) > 0 )
93  {
94  dEtaSq = std::log( ( mom1 - pz1 ) / ( mom1 + pz1 )
95  * ( mom2 + pz2 ) / ( mom2 - pz2 ) );
96  dEtaSq = 0.25 * dEtaSq * dEtaSq;
97  }
98  }
99  else if ( p1.kind() == I4Momentum::P4PXPYPZE )
100  {
101  const double mom1 = p1.p(); const double pz1 = p1.pz();
102  if ( mom1 + pz1 > 0 )
103  {
104  dEtaSq = - 0.5 * std::log( ( mom1 - pz1 ) / ( mom1 + pz1 ) ) - p2.eta();
105  dEtaSq = dEtaSq * dEtaSq;
106  }
107  }
108  else if ( p2.kind() == I4Momentum::P4PXPYPZE )
109  {
110  const double mom2 = p2.p(); const double pz2 = p2.pz();
111  if ( mom2 - pz2 > 0 )
112  {
113  dEtaSq = p1.eta() - 0.5 * std::log( ( mom2 + pz2 ) / ( mom2 - pz2 ) );
114  dEtaSq = dEtaSq * dEtaSq;
115  }
116  }
117  else
118  {
119  dEtaSq = p1.eta() - p2.eta();
120  dEtaSq = dEtaSq * dEtaSq;
121  }
122  return dEtaSq;
123  }
124 
126  inline
127  double deltaEta( const I4Momentum * const p1, const I4Momentum * const p2 )
128  {
129  return deltaEta(*p1, *p2);
130  }
131 
133  inline
134  double deltaPhi( const I4Momentum& p4, const double phi )
135  {
136  return P4Helpers::deltaPhi( p4.phi(), phi );
137  }
138 
140  inline
141  double deltaPhi( const I4Momentum & pA, const I4Momentum & pB )
142  {
143  double dPhi(999);
144  if ( pA.kind() == I4Momentum::P4PXPYPZE && pB.kind() == I4Momentum::P4PXPYPZE )
145  {
146  double xp = pA.px() * pB.px() ;
147  double xyp1 = 0;
148  if ( xp != 0 )
149  {
150  xyp1 = 1.0 + pA.py() * pB.py() / xp ;
151  dPhi = 0;
152  if ( xyp1 != 0 )
153  dPhi = atan( ( pA.py() / pA.px() - pB.py() / pB.px() ) / xyp1 );
154  if ( xyp1 < 0 )
155  {
156  if ( pA.py() / pA.px() > 0 )
157  dPhi += M_PI;
158  else
159  dPhi -= M_PI;
160  }
161  }
162  }
163  else
164  {
165  dPhi = remainder( -pA.phi() + pB.phi(), 2*M_PI );
166  }
167  return dPhi;
168  }
169 
171  inline
172  double deltaPhiSq( const I4Momentum & pA, const I4Momentum & pB )
173  {
174  double dPhiSq(999);
175  if ( pA.kind() == I4Momentum::P4PXPYPZE && pB.kind() == I4Momentum::P4PXPYPZE )
176  {
177  double cosPhi = ( pA.px() * pB.px() + pA.py() * pB.py() ) / pA.pt() / pB.pt();
178  double phi = acos(cosPhi);
179  dPhiSq = phi*phi;
180  }
181  else
182  {
183  dPhiSq = remainder( -pA.phi() + pB.phi(), 2*M_PI );
184  dPhiSq = dPhiSq * dPhiSq;
185  }
186  return dPhiSq;
187  }
188 
190  inline
191  double deltaPhi( const I4Momentum * const pA, const I4Momentum * const pB )
192  { return deltaPhi( pA->phi(), pB->phi() ); }
193 
195  inline
196  double deltaR( const I4Momentum& p4, double eta, double phi )
197  {
198  using std::sqrt;
199  const double dEta = p4.eta() - eta;
200  const double dPhi = P4Helpers::deltaPhi( p4, phi );
201  return sqrt( dEta*dEta + dPhi*dPhi );
202  }
203 
205  inline
206  double deltaR( const I4Momentum& pA, const I4Momentum& pB )
207  {
208  using std::sqrt;
209  const double dEtaSq = P4Helpers::deltaEtaSq( pA, pB );
210  const double dPhiSq = P4Helpers::deltaPhiSq( pA, pB );
211  return sqrt( dEtaSq + dPhiSq );
212  }
213 
215  inline
216  double deltaR( const I4Momentum * const pA, const I4Momentum * const pB )
217  { return P4Helpers::deltaR( *pA, *pB ); }
218 
222  inline
223  bool isInDeltaR( const I4Momentum& p1, const I4Momentum& p2,
224  double dR )
225  {
226  using std::abs;
227  using std::sqrt;
228  const double dPhi = abs( P4Helpers::deltaPhi(p1,p2) ); // in [-pi,pi)
229  if ( dPhi > dR ) return false; // <==
230  const double dEta = abs( P4Helpers::deltaEta(p1,p2) );
231  if ( dEta > dR ) return false; // <==
232  const double deltaR2 = dEta*dEta + dPhi*dPhi;
233  if ( deltaR2 > dR*dR ) return false; // <==
234  return true;
235  }
236 
238  inline
239  double invMass( const I4Momentum & pA, const I4Momentum & pB )
240  { return (pA.hlv()+pB.hlv()).m(); }
241 
243  inline
244  double invMass( const I4Momentum * const pA, const I4Momentum * const pB )
245  { return invMass( *pA, *pB ); }
246 
248  inline
249  double invMass( const I4Momentum & pA, const I4Momentum & pB,
250  const I4Momentum & pC, const I4Momentum & pD )
251  { return (pA.hlv()+pB.hlv()+pC.hlv()+pD.hlv()).m(); }
252 
254  inline
255  double invMass( const I4Momentum * const pA, const I4Momentum * const pB,
256  const I4Momentum * const pC, const I4Momentum * const pD )
257  { return invMass( *pA, *pB, *pC, *pD ); }
258 
259  // -------------------------------------------------------------
260  // --------------- Sorting functions ---------------------------
261  // -------------------------------------------------------------
262 
264  template<class Iterator_t, class Predicate_t>
265  inline
266  void sort( Iterator_t itr, Iterator_t itrEnd, Predicate_t p )
267  {
268  // Koenig's look-up at our rescue: handle correctly DataVector's
269  // special sort method. => We inject the std::sort into our namespace
270  using std::sort;
271  return sort( itr, itrEnd, p );
272  }
273 
275  template<class Container_t, class Predicate_t>
276  inline
277  void sort( Container_t& container, Predicate_t p )
278  {
279  return P4Helpers::sort( container.begin(), container.end(), p );
280  }
281 
282  // -------------------------------------------------------------
283  // --------------- Matching functions --------------------------
284  // -------------------------------------------------------------
285 
290  template <class Container_t>
291  inline
292  bool closestDeltaR( const double eta, const double phi,
293  Container_t& coll, std::size_t& index, double& deltaR )
294  {
295  deltaR = std::numeric_limits<double>::max(); // big value
296  bool l_return = false;
297  std::size_t l_idx = 0;
298  typename Container_t::const_iterator it = coll.begin();
299  typename Container_t::const_iterator itE = coll.end();
300  for (; it != itE; ++it,++l_idx) {
301  double rtu = P4Helpers::deltaR(**it, eta, phi);
302  if ( CxxUtils::fpcompare::less(rtu, deltaR) ) {
303  index = l_idx;
304  deltaR = rtu;
305  l_return = true;
306  }
307  }
308  return l_return;
309  }
310 
315  template <class Container_t>
316  inline
317  bool closestDeltaR( const I4Momentum& p4,
318  Container_t& coll, std::size_t& index, double& deltaR )
319  {
320  return P4Helpers::closestDeltaR( p4.eta(), p4.phi(),
321  coll, index, deltaR );
322  }
323 
329  template <class Container_t>
330  inline
331  bool closestDeltaR( const double eta, const double phi, const double ene,
332  Container_t& coll, std::size_t& index,
333  double &deltaR, double &deltaE )
334  {
335  using std::abs;
336  deltaR = deltaE = std::numeric_limits<double>::max(); // big value
337  bool l_return = false;
338  std::size_t l_idx = 0;
339  typename Container_t::const_iterator it = coll.begin();
340  typename Container_t::const_iterator itE = coll.end();
341  for (; it != itE; ++it,++l_idx) {
342  const double dE = abs( ene - (*it)->e() );
343  if ( CxxUtils::fpcompare::less(dE, deltaE) ) {
344  const double rtu = P4Helpers::deltaR(**it,eta,phi);
345  if ( CxxUtils::fpcompare::less(rtu, deltaR) ) {
346  index = l_idx;
347  deltaR = rtu;
348  deltaE = dE;
349  l_return = true;
350  }
351  }
352  }
353  return l_return;
354  }
355 
361  template <class Container_t>
362  inline
363  bool closestDeltaR( const I4Momentum& p4,
364  Container_t& coll, std::size_t& index,
365  double &deltaR, double &deltaE )
366  {
367  return P4Helpers::closestDeltaR( p4.eta(), p4.phi(), p4.e(),
368  coll, index, deltaR, deltaE );
369  }
370 
371 } //> namespace P4Helpers
372 
373 
374 #endif
375 
376 #endif // FOURMOMUTILS_P4HELPERS_H
I4Momentum::py
virtual double py() const =0
y component of momentum
I4Momentum
Definition: I4Momentum.h:31
P4Helpers::deltaEtaSq
double deltaEtaSq(const I4Momentum &p1, const I4Momentum &p2)
Computes efficiently .
Definition: P4Helpers.h:85
P4Helpers::closestDeltaR
bool closestDeltaR(const double eta, const double phi, Container_t &coll, std::size_t &index, double &deltaR)
Find the closest element in a collection to an I4Momentum.
Definition: P4Helpers.h:292
python.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
max
#define max(a, b)
Definition: cfImp.cxx:41
P4Helpers::invMass
double invMass(const I4Momentum &pA, const I4Momentum &pB)
invariant mass from two I4momentum references
Definition: P4Helpers.h:239
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:64
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:79
index
Definition: index.py:1
I4Momentum::p
virtual double p() const =0
momentum magnitude
skel.it
it
Definition: skel.GENtoEVGEN.py:423
M_PI
#define M_PI
Definition: ActiveFraction.h:11
P4Helpers
P4Helpers provides static helper functions for kinematic calculation on objects deriving from I4Momen...
Definition: P4Helpers.h:26
xAOD::P4Helpers::deltaR2
double deltaR2(double rapidity1, double phi1, double rapidity2, double phi2)
from bare rapidity,phi
Definition: xAODP4Helpers.h:111
P4Helpers::isInDeltaR
bool isInDeltaR(const I4Momentum &p1, const I4Momentum &p2, double dR)
Check if 2 I4Momentum are in a cone.
Definition: P4Helpers.h:223
I4Momentum::pt
virtual double pt() const =0
transverse momentum
I4Momentum::hlv
virtual CLHEP::HepLorentzVector hlv() const =0
CLHEP HepLorentzVector.
I4Momentum.h
TruthTest.itE
itE
Definition: TruthTest.py:25
drawFromPickle.atan
atan
Definition: drawFromPickle.py:36
P4Helpers::deltaPhi
double deltaPhi(double phiA, double phiB)
delta Phi in range [-pi,pi[
Definition: P4Helpers.h:29
I4Momentum::pz
virtual double pz() const =0
z component of momentum
P4Helpers::deltaEta
double deltaEta(const I4Momentum &p1, const I4Momentum &p2)
Computes efficiently .
Definition: P4Helpers.h:53
I4Momentum::e
virtual double e() const =0
energy
TauGNNUtils::Variables::Track::dPhi
bool dPhi(const xAOD::TauJet &tau, const xAOD::TauTrack &track, double &out)
Definition: TauGNNUtils.cxx:530
fpcompare.h
Workaround x86 precision issues for FP inequality comparisons.
I4Momentum::eta
virtual double eta() const =0
pseudo rapidity
P4Helpers::deltaR
double deltaR(const I4Momentum &p4, double eta, double phi)
from 1 I4Momentum
Definition: P4Helpers.h:196
I4Momentum::phi
virtual double phi() const =0
phi in [-pi,pi[
I4Momentum::kind
virtual Kind kind() const =0
tells what kind of P4XYZT this is
P4Helpers::deltaPhiSq
double deltaPhiSq(const I4Momentum &pA, const I4Momentum &pB)
delta Phi squared in range ([-pi,pi[)^2 from two I4momentum references
Definition: P4Helpers.h:172
remainder
std::vector< std::string > remainder(const std::vector< std::string > &v1, const std::vector< std::string > &v2)
list of entries in a vector that are not in another
Definition: compareFlatTrees.cxx:44
CxxUtils::fpcompare::less
bool less(double a, double b)
Compare two FP numbers, working around x87 precision issues.
Definition: fpcompare.h:166
P4Helpers::sort
void sort(Iterator_t itr, Iterator_t itrEnd, Predicate_t p)
sort a container according to the given predicate
Definition: P4Helpers.h:266
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
P4Helpers::sort
void sort(Container_t &container, Predicate_t p)
sort a container according to the given predicate
Definition: P4Helpers.h:277
TauGNNUtils::Variables::Track::dEta
bool dEta(const xAOD::TauJet &tau, const xAOD::TauTrack &track, double &out)
Definition: TauGNNUtils.cxx:525
I4Momentum::px
virtual double px() const =0
x component of momentum
I4Momentum::P4PXPYPZE
@ P4PXPYPZE
Definition: I4Momentum.h:33
Generated on Tue Jun 25 2024 21:23:45 for ATLAS Offline Software by   doxygen 1.8.18