ATLAS Offline Software
Loading...
Searching...
No Matches
Typedefs | Functions
FourMomUtils Namespace Reference

Typedefs

typedef INavigable4MomentumCollection::const_iterator I4MomIter_t

Functions

double forwardTerm (const I4MomIter_t &iBeg, const I4MomIter_t &iEnd, double m_central, bool useThreeD=false)
double forwardTerm (const INavigable4MomentumCollection *particles, double central, bool useThreeD=false)
bool foxWolfram (const I4MomIter_t &iBeg, const I4MomIter_t &iEnd, std::vector< double > &H, unsigned int order=5)
bool foxWolfram (const INavigable4MomentumCollection *theParticles, std::vector< double > &H, unsigned int order=5)
bool jetBroadening (const I4MomIter_t &iBeg, const I4MomIter_t &iEnd, double &wideJetBroadening, double &totalJetBroadening, CLHEP::Hep3Vector thrust, bool useThreeD=false)
bool jetBroadening (const INavigable4MomentumCollection *theParticles, double &wideJetBroadening, double &totalJetBroadening, CLHEP::Hep3Vector thrust, bool useThreeD=false)
bool jetMasses (const I4MomIter_t &iBeg, const I4MomIter_t &iEnd, double &heavyJetMass, double &lightJetMass, const CLHEP::Hep3Vector &thrust)
bool jetMasses (const INavigable4MomentumCollection *theParticles, double &heavyJetMass, double &lightJetMass, CLHEP::Hep3Vector thrust)
template<class I4MomIter>
std::ostream & dump (std::ostream &out, const I4MomIter iBeg, const I4MomIter iEnd)
 Helper to stream out a range of I4Momentum objects.
template<class Container>
std::ostream & dump (std::ostream &out, const Container &c)
 Helper to stream out a collection of I4Momentum objects.
CLHEP::Hep3Vector thrust (const I4MomIter_t &iBeg, const I4MomIter_t &iEnd, double &thrust_major, double &thrust_minor, bool useThreeD=false)
CLHEP::Hep3Vector thrust (const INavigable4MomentumCollection *theParticles, double &thrust_major, double &thrust_minor, bool useThreeD=false)

Typedef Documentation

◆ I4MomIter_t

typedef INavigable4MomentumCollection::const_iterator FourMomUtils::I4MomIter_t

Definition at line 16 of file ForwardTerm.h.

Function Documentation

◆ dump() [1/2]

template<class Container>
std::ostream & FourMomUtils::dump ( std::ostream & out,
const Container & c )
inline

Helper to stream out a collection of I4Momentum objects.

Definition at line 35 of file P4Dumper.h.

35 {
36 return dump( out, c.begin(), c.end() );
37 }
dump
-event-from-file

◆ dump() [2/2]

template<class I4MomIter>
std::ostream & FourMomUtils::dump ( std::ostream & out,
const I4MomIter iBeg,
const I4MomIter iEnd )
inline

Helper to stream out a range of I4Momentum objects.

Definition at line 24 of file P4Dumper.h.

25 {
26 for ( I4MomIter i = iBeg; i != iEnd; ++i ) {
27 out << (**i) << "\n";
28 }
29 return out;
30 }

◆ forwardTerm() [1/2]

double FourMomUtils::forwardTerm ( const I4MomIter_t & iBeg,
const I4MomIter_t & iEnd,
double m_central,
bool useThreeD = false )

Definition at line 18 of file ForwardTerm.cxx.

21{
22 double Q = 0;
23 double suppressed = 0;
24 double eta = 0;
25
26 // determine average eta for event
27 for ( I4MomIter_t itr = iBeg; itr != iEnd; ++itr )
28 {
29 if(abs((*itr)->eta())<central)
30 {
31 // default : standard definition of z component zero
32 double z=0;
33 if(useThreeD)
34 z=(*itr)->pz();
35
36 CLHEP::Hep3Vector c( (*itr)->px(), (*itr)->py(), z );
37 Q += c.mag();
38 eta += (*itr)->eta() * c.mag();
39 }
40 }
41
42 if (Q <= 0)
43 return 0;
44 const double inv_Q = 1. / Q;
45
46 eta *= inv_Q;
47
48 for ( I4MomIter_t itr = iBeg; itr != iEnd; ++itr )
49 {
50 if(abs((*itr)->eta())>central)
51 {
52 // default : standard definition of z component zero
53 double z=0;
54 if(useThreeD)
55 z=(*itr)->pz();
56
57 CLHEP::Hep3Vector c( (*itr)->px(), (*itr)->py(), z );
58 suppressed += c.mag() * exp( -abs( eta - (*itr)->eta() ));
59 }
60 }
61
62 return suppressed * inv_Q;
63}
eta
Scalar eta() const
pseudorapidity method
Definition AmgMatrixBasePlugin.h:83
z
#define z
FourMomUtils::I4MomIter_t
INavigable4MomentumCollection::const_iterator I4MomIter_t
Definition ForwardTerm.h:16

◆ forwardTerm() [2/2]

double FourMomUtils::forwardTerm ( const INavigable4MomentumCollection * particles,
double central,
bool useThreeD = false )
inline

Definition at line 24 of file ForwardTerm.h.

25 {
26 return forwardTerm( particles->begin(),
27 particles->end(),
28 central, useThreeD );
29 }
FourMomUtils::forwardTerm
double forwardTerm(const I4MomIter_t &iBeg, const I4MomIter_t &iEnd, double m_central, bool useThreeD=false)
Definition ForwardTerm.cxx:18

◆ foxWolfram() [1/2]

bool FourMomUtils::foxWolfram ( const I4MomIter_t & iBeg,
const I4MomIter_t & iEnd,
std::vector< double > & H,
unsigned int order = 5 )

Definition at line 20 of file Event/FourMomUtils/src/FoxWolfram.cxx.

22{
23 // this is the vector of results
24 H.resize(order, 0);
25 if(order==0) return true;
26
27 double N0=0;
28
29 // store the first four values of the legendre polynomials in vector,
30 // then use recursive formula to calculate others (just need two previous elements)
31 std::vector<double> P(order, 1);
32
33 for ( I4MomIter_t itr_i = iBeg; itr_i != iEnd; ++itr_i )
34 {
35 CLHEP::Hep3Vector ci( (*itr_i)->px(), (*itr_i)->py(), (*itr_i)->pz() );
36
37 for ( I4MomIter_t itr_j = iBeg; itr_j != iEnd; ++itr_j )
38 {
39 CLHEP::Hep3Vector cj( (*itr_j)->px(), (*itr_j)->py(), (*itr_j)->pz() );
40 double x=cos(ci.angle(cj));
41
42 double P0=1;
43 double P1=x;
44 double P2=0.5*(3.0*x*x-1);
45 double P3=0.5*(5.0*x*x*x-3.0*x);
46 double P4=0.125*(35.0*x*x*x*x-30*x*x+3);
47
48 H[0]+=abs(ci.mag())*abs(cj.mag())*P0;
49 if(order>=1)
50 {
51 H[1]+=abs(ci.mag())*abs(cj.mag())*P1;
52 P[1]=P1; // never used
53 }
54 if(order>=2)
55 {
56 H[2]+=abs(ci.mag())*abs(cj.mag())*P2;
57 P[2]=P2; // never used
58 }
59 if(order>=3)
60 {
61 H[3]+=abs(ci.mag())*abs(cj.mag())*P3;
62 P[3]=P3;
63 }
64 if(order>=4)
65 {
66 H[4]+=abs(ci.mag())*abs(cj.mag())*P4;
67 P[4]=P4;
68 }
69
70 for ( unsigned int loop=5; loop<order; ++loop )
71 {
72 P[loop] = 1.0/loop*((2.0*loop-1)*x*P[loop-1]-
73 (loop-1)*P[loop-2]);
74 H[loop]=abs(ci.mag())*abs(ci.mag())*P[loop];
75 }
76 }
77 N0+=(*itr_i)->e();
78 }
79
80 N0=N0*N0;
81
82 // and normalize
83 const double inv_N0 = N0!=0 ? (1. / N0) : 1;
84 for ( unsigned int loop=5; loop<order; ++loop ) {
85 H[loop] *= inv_N0;
86 }
87
88 return true;
89}
P
static Double_t P(Double_t *tt, Double_t *par)
Definition LArPhysWaveHECTool.cxx:280
H
#define H(x, y, z)
Definition MD5.cxx:114
x
#define x
mc.order
order
Configure Herwig7.
Definition mc.Herwig7_Dijet.py:12

◆ foxWolfram() [2/2]

bool FourMomUtils::foxWolfram ( const INavigable4MomentumCollection * theParticles,
std::vector< double > & H,
unsigned int order = 5 )
inline

Definition at line 25 of file Event/FourMomUtils/FourMomUtils/FoxWolfram.h.

26 {
27 return foxWolfram( theParticles->begin(), theParticles->end(),
28 H, order );
29 }
DataVector::end
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
DataVector::begin
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
FourMomUtils::foxWolfram
bool foxWolfram(const I4MomIter_t &iBeg, const I4MomIter_t &iEnd, std::vector< double > &H, unsigned int order=5)
Definition Event/FourMomUtils/src/FoxWolfram.cxx:20

◆ jetBroadening() [1/2]

bool FourMomUtils::jetBroadening ( const I4MomIter_t & iBeg,
const I4MomIter_t & iEnd,
double & wideJetBroadening,
double & totalJetBroadening,
CLHEP::Hep3Vector thrust,
bool useThreeD = false )

Definition at line 17 of file JetBroadening.cxx.

20{
21 if(abs(thrust.mag()-1)>0.01)
22 return false;
23
24 double Qu=0;
25 double Qd=0;
26 double etau=0;
27 double etad=0;
28 double phiu=0;
29 double phid=0;
30
31 // ensure z component is zero
32 if(!useThreeD)
33 thrust.setZ(0);
34
35 // determine average eta/phi for each hemisphere
36 for ( I4MomIter_t itr = iBeg; itr != iEnd; ++itr )
37 {
38 // default : standard definition of z component zero
39 double z=0;
40 if(useThreeD)
41 z=(*itr)->pz();
42
43 CLHEP::Hep3Vector c( (*itr)->px(), (*itr)->py(), z );
44 if(c.dot(thrust)>0)
45 {
46 Qu += c.mag();
47 etau += (*itr)->eta() * c.mag();
48 phiu += (*itr)->phi() * c.mag();
49 }else{
50 Qd += c.mag();
51 etad += (*itr)->eta() * c.mag();
52 phid += (*itr)->phi() * c.mag();
53 }
54 }
55
56 const double inv_Qu = Qu!=0 ? (1. / Qu) : 1;
57 const double inv_Qd = Qd!=0 ? (1. / Qd) : 1;
58
59 etau *= inv_Qu;
60 etad *= inv_Qd;
61
62 phiu *= inv_Qu;
63 phid *= inv_Qd;
64
65 double Bu=0;
66 double Bd=0;
67
68 // now detemine JetBroadenings in each hemisphere
69 for ( I4MomIter_t itr = iBeg; itr != iEnd; ++itr )
70 {
71 // default : standard definition of z component zero
72 double z=0;
73 if(useThreeD)
74 z=(*itr)->pz();
75
76 CLHEP::Hep3Vector c( (*itr)->px(), (*itr)->py(), z );
77
78 if(c.dot(thrust)>0)
79 {
80 Bu += c.mag() *
81 sqrt((etau-(*itr)->eta())*(etau-(*itr)->eta())+
82 (phiu-(*itr)->phi())*(phiu-(*itr)->phi()));
83 }else{
84 Bd += c.mag() *
85 sqrt((etad-(*itr)->eta())*(etad-(*itr)->eta())+
86 (phid-(*itr)->phi())*(phid-(*itr)->phi()));
87 }
88 }
89
90 if (Qu+Qd != 0) {
91 Bu /= 2*(Qu+Qd);
92 Bd /= 2*(Qu+Qd);
93 }
94
95 totalJetBroadening=Bu+Bd;
96 wideJetBroadening=Bu;
97 if(Bd>Bu)
98 wideJetBroadening=Bd;
99
100 return true;
101}
FourMomUtils::thrust
CLHEP::Hep3Vector thrust(const I4MomIter_t &iBeg, const I4MomIter_t &iEnd, double &thrust_major, double &thrust_minor, bool useThreeD=false)
Definition Event/FourMomUtils/src/Thrust.cxx:19

◆ jetBroadening() [2/2]

bool FourMomUtils::jetBroadening ( const INavigable4MomentumCollection * theParticles,
double & wideJetBroadening,
double & totalJetBroadening,
CLHEP::Hep3Vector thrust,
bool useThreeD = false )
inline

Definition at line 23 of file JetBroadening.h.

25 {
26 return jetBroadening( theParticles->begin(), theParticles->end(),
27 wideJetBroadening, totalJetBroadening,
28 thrust, useThreeD );
29 }
FourMomUtils::jetBroadening
bool jetBroadening(const I4MomIter_t &iBeg, const I4MomIter_t &iEnd, double &wideJetBroadening, double &totalJetBroadening, CLHEP::Hep3Vector thrust, bool useThreeD=false)
Definition JetBroadening.cxx:17

◆ jetMasses() [1/2]

bool FourMomUtils::jetMasses ( const I4MomIter_t & iBeg,
const I4MomIter_t & iEnd,
double & heavyJetMass,
double & lightJetMass,
const CLHEP::Hep3Vector & thrust )

Definition at line 18 of file JetMasses.cxx.

21{
22 if(abs(thrust.mag()-1)>0.01)
23 return false;
24
25 // split event into upper and lower hemisphere
26 CLHEP::Hep3Vector up(0,0,0);
27 CLHEP::Hep3Vector down(0,0,0);
28
29 double Q=0;
30
31 for ( I4MomIter_t itr = iBeg; itr != iEnd; ++itr )
32 {
33 CLHEP::Hep3Vector c( (*itr)->px(), (*itr)->py(), (*itr)->pz() );
34 Q += sqrt( (*itr)->px()*(*itr)->px()+(*itr)->py()*(*itr)->py() );
35
36 if(c.dot(thrust)>0)
37 {
38 up += c;
39 }else{
40 down += c;
41 }
42 }
43
44 const double inv_Q2 = Q!=0 ? (1. / (Q*Q)) : 0;
45 lightJetMass=down.mag2()*inv_Q2;
46 heavyJetMass=up.mag2()*inv_Q2;
47 if ( lightJetMass > heavyJetMass )
48 {
49 lightJetMass=up.mag2()*inv_Q2;
50 heavyJetMass=down.mag2()*inv_Q2;
51 }
52
53 return true;
54}

◆ jetMasses() [2/2]

bool FourMomUtils::jetMasses ( const INavigable4MomentumCollection * theParticles,
double & heavyJetMass,
double & lightJetMass,
CLHEP::Hep3Vector thrust )
inline

Definition at line 23 of file JetMasses.h.

25 {
26 return jetMasses( theParticles->begin(), theParticles->end(),
27 heavyJetMass, lightJetMass, thrust );
28 }
FourMomUtils::jetMasses
bool jetMasses(const I4MomIter_t &iBeg, const I4MomIter_t &iEnd, double &heavyJetMass, double &lightJetMass, const CLHEP::Hep3Vector &thrust)
Definition JetMasses.cxx:18

◆ thrust() [1/2]

CLHEP::Hep3Vector FourMomUtils::thrust ( const I4MomIter_t & iBeg,
const I4MomIter_t & iEnd,
double & thrust_major,
double & thrust_minor,
bool useThreeD = false )

Definition at line 19 of file Event/FourMomUtils/src/Thrust.cxx.

21{
22 /*
23 Finding the thrust axis in an event is not trivial.
24
25 Here, we follow the procedure described in the PYTHIA manual JHEP 05 (2006) 026,
26 also hep-ph/0603175. The approach is to use an iterative method, which usually
27 converges quickly. As the minimization can find just a local minimum, different
28 starting points for the thrust axis are tried. By default, first the direction
29 of the four most energetic particles are tried, if their result disagrees, all
30 16 permutations of the sum of all 4 particles are tried (with coefficients +- 1)
31
32 Note, that the thrust is calculated for _ALL_ particles. If you want only a subset
33 of particles, you have to apply a cut beforehand.
34 See e.g. Reconstruction/EventShapes/EventShapeTools for examples.
35 */
36
37 CLHEP::Hep3Vector thrust(0,0,0);
38
39 int agree=0;
40 int disagree=0;
41
42 CLHEP::Hep3Vector n_0[20];
43 short add0[20] = { 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1,-1,-1 };
44 short add1[20] = { 0, 1, 0, 0, 1, 1, 1, 1,-1,-1,-1,-1, 1, 1, 1, 1,-1,-1,-1,-1 };
45 short add2[20] = { 0, 0, 1, 0, 1, 1,-1,-1, 1, 1,-1,-1, 1, 1,-1,-1, 1, 1,-1,-1 };
46 short add3[20] = { 0, 0, 0, 1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1 };
47
48 std::vector<const INavigable4Momentum*> v_copy(4);
49
50 // partial_sort_copy sorts a copy of the collection according to energy and
51 // returns only the first four elements (minimum of input collection size and
52 // pre-allocated output vector v_copy
53 partial_sort_copy( iBeg, iEnd,
54 v_copy.begin(), v_copy.end(),
55 P4Sorters::Descending::Ene() );
56
57 int n_tests=0;
58 int max_tests=std::min<int>(20, std::distance(iBeg, iEnd));
59 do{
60 n_0[n_tests]=CLHEP::Hep3Vector(0,0,0);
61
62 // assign direction of four most energetic particles
63 n_0[n_tests] +=
64 add0[n_tests] * CLHEP::Hep3Vector(v_copy.at(0)->px(), v_copy.at(0)->py(), v_copy.at(0)->pz()) +
65 add1[n_tests] * CLHEP::Hep3Vector(v_copy.at(1)->px(), v_copy.at(1)->py(), v_copy.at(1)->pz()) +
66 add2[n_tests] * CLHEP::Hep3Vector(v_copy.at(2)->px(), v_copy.at(2)->py(), v_copy.at(2)->pz()) +
67 add3[n_tests] * CLHEP::Hep3Vector(v_copy.at(3)->px(), v_copy.at(3)->py(), v_copy.at(3)->pz());
68
69 if(!useThreeD)
70 n_0[n_tests].setZ(0);
71
72 /* // my convention : x is always positive (thrust axis has two fold ambiguity)
73 if(n_0[n_tests].x()<0)
74 n_0[n_tests] = - n_0[n_tests]; */
75
76 // protect against small number of input particles (smaller than 4!)
77 if(n_0[n_tests].mag()>0)
78 n_0[n_tests] /= n_0[n_tests].mag();
79
80 int loop=0;
81 bool run=false;
82 do{
83 CLHEP::Hep3Vector n_1(0,0,0);
84 for ( I4MomIter_t itr = iBeg; itr != iEnd; ++itr )
85 {
86 // if(((*itr)->hlv()).vect().dot(n_0[n_tests])>0)
87 if((*itr)->px() * n_0[n_tests].x() +
88 (*itr)->py() * n_0[n_tests].y() +
89 (*itr)->pz() * n_0[n_tests].z() > 0 )
90 n_1 += CLHEP::Hep3Vector((*itr)->px(), (*itr)->py(), (*itr)->pz() );
91 else
92 n_1 -= CLHEP::Hep3Vector((*itr)->px(), (*itr)->py(), (*itr)->pz() );
93 }
94
95 if(!useThreeD)
96 n_1.setZ(0);
97
98 // protect against few number of input particles (smaller than 4!)
99 if(n_1.mag()>0)
100 n_1 /= n_1.mag();
101
102 // has axis changed ? if so, try at most ten times (thrust axis has two fold ambiguity)
103 run = ( n_0[n_tests]!=n_1 ) && ( -n_0[n_tests]!=n_1 ) && loop++<10;
104 n_0[n_tests] = n_1;
105 }while( run );
106
107 // agrees or disagrees with first result ?
108 // thrust has a sign ambiguity
109 if( n_tests>0 && ( n_0[0] == n_0[n_tests] || n_0[0] == -n_0[n_tests] ) ) agree++;
110 if( n_tests>0 && n_0[0] != n_0[n_tests] && n_0[0] != -n_0[n_tests] ) disagree++;
111
112 // stop if four first tries give same result (no test for first try! )
113 // if not, try at most max_tests combinations
114 }while( ( disagree>0 || agree<4 ) && ++n_tests < max_tests);
115
116 // now that we have the thrust axis, we determine the thrust value
117 // if the various calculations of the thrust axes disagree, try all
118 // and take the maximum, calculate minor and mayor axis
119 n_tests=0;
120 do{
121 double denominator = 0;
122 double numerator_t = 0;
123 double numerator_m = 0;
124 for ( I4MomIter_t itr = iBeg; itr != iEnd; ++itr )
125 {
126 CLHEP::Hep3Vector c((*itr)->hlv().vect());
127 c.setZ(0);
128 numerator_t += abs(c.dot(n_0[n_tests]));
129 numerator_m += (c.cross(n_0[n_tests])).mag();
130 denominator += c.mag();
131 }
132 const double inv_denominator = denominator!= 0 ? (1. / denominator) : 1;
133 if( numerator_t * inv_denominator > thrust_major )
134 {
135 thrust_major = numerator_t * inv_denominator;
136 thrust_minor = numerator_m * inv_denominator;
137 thrust=n_0[n_tests];
138 }
139 }while(disagree>0 && ++n_tests < max_tests);
140
141 /* std::cout << "Calculation of Thrust gave: ( "
142 << thrust.x() << " | "
143 << thrust.y() << " | "
144 << thrust.z() << " )\n"; */
145
146 // return StatusCode::SUCCESS;
147 return thrust;
148}
mag
Scalar mag() const
mag method
Definition AmgMatrixBasePlugin.h:26

◆ thrust() [2/2]

CLHEP::Hep3Vector FourMomUtils::thrust ( const INavigable4MomentumCollection * theParticles,
double & thrust_major,
double & thrust_minor,
bool useThreeD = false )
inline

Definition at line 24 of file Event/FourMomUtils/FourMomUtils/Thrust.h.

26 {
27 return thrust( theParticles->begin(), theParticles->end(),
28 thrust_major, thrust_minor,
29 useThreeD );
30 }
Generated on for ATLAS Offline Software by doxygen 1.14.0