ATLAS Offline Software
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. More...
 
template<class Container >
std::ostream & dump (std::ostream &out, const Container &c)
 Helper to stream out a collection of I4Momentum objects. More...
 
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() [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 }

◆ 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  }

◆ 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 = 1. / N0;
84  for ( unsigned int loop=5; loop<order; ++loop ) {
85  H[loop] *= inv_N0;
86  }
87 
88  return true;
89 }

◆ 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  }

◆ 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 = 1. / Qu;
57  const double inv_Qd = 1. / Qd;
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  Bu /= 2*(Qu+Qd);
91  Bd /= 2*(Qu+Qd);
92 
93  totalJetBroadening=Bu+Bd;
94  wideJetBroadening=Bu;
95  if(Bd>Bu)
96  wideJetBroadening=Bd;
97 
98  return true;
99 }

◆ 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  }

◆ 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 = 1. / (Q*Q);
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  }

◆ 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 = 1. / denominator;
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 }

◆ 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  }
FourMomUtils::forwardTerm
double forwardTerm(const INavigable4MomentumCollection *particles, double central, bool useThreeD=false)
Definition: ForwardTerm.h:24
FourMomUtils::dump
std::ostream & dump(std::ostream &out, const Container &c)
Helper to stream out a collection of I4Momentum objects.
Definition: P4Dumper.h:35
FourMomUtils::jetBroadening
bool jetBroadening(const INavigable4MomentumCollection *theParticles, double &wideJetBroadening, double &totalJetBroadening, CLHEP::Hep3Vector thrust, bool useThreeD=false)
Definition: JetBroadening.h:23
checkPlugins.suppressed
suppressed
Definition: checkPlugins.py:206
P4Sorters::Descending::Ene
Definition: P4DescendingSorters.h:105
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:83
DMTest::P
P_v1 P
Definition: P.h:23
python.AthDsoLogger.out
out
Definition: AthDsoLogger.py:71
FourMomUtils::I4MomIter_t
INavigable4MomentumCollection::const_iterator I4MomIter_t
Definition: ForwardTerm.h:16
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
drawFromPickle.exp
exp
Definition: drawFromPickle.py:36
x
#define x
H
#define H(x, y, z)
Definition: MD5.cxx:114
FourMomUtils::jetMasses
bool jetMasses(const INavigable4MomentumCollection *theParticles, double &heavyJetMass, double &lightJetMass, CLHEP::Hep3Vector thrust)
Definition: JetMasses.h:23
lumiFormat.i
int i
Definition: lumiFormat.py:85
z
#define z
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
mc.order
order
Configure Herwig7.
Definition: mc.Herwig7_Dijet.py:12
CalibCoolCompareRT.up
up
Definition: CalibCoolCompareRT.py:109
run
Definition: run.py:1
FourMomUtils::thrust
CLHEP::Hep3Vector thrust(const INavigable4MomentumCollection *theParticles, double &thrust_major, double &thrust_minor, bool useThreeD=false)
Definition: Event/FourMomUtils/FourMomUtils/Thrust.h:24
ReadTripsProbsFromCool.denominator
denominator
Definition: ReadTripsProbsFromCool.py:96
python.pydraw.loop
def loop(arg)
Definition: pydraw.py:1242
FourMomUtils::foxWolfram
bool foxWolfram(const INavigable4MomentumCollection *theParticles, std::vector< double > &H, unsigned int order=5)
Definition: Event/FourMomUtils/FourMomUtils/FoxWolfram.h:25
DataVector::end
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
LArG4FSStartPointFilter.particles
list particles
Definition: LArG4FSStartPointFilter.py:84
Amg::distance
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
Definition: GeoPrimitivesHelpers.h:54
python.compressB64.c
def c
Definition: compressB64.py:93
mag
Scalar mag() const
mag method
Definition: AmgMatrixBasePlugin.h:26
DataVector::begin
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
Generated on Thu Nov 7 2024 22:14:34 for ATLAS Offline Software by   doxygen 1.8.18