ATLAS Offline Software
T2TrackBSLLPoly.cxx
Go to the documentation of this file.
1 /*
2 Copyright (C) 2019-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include <cmath>
6 #include <cstdint>
7 
8 #include "T2TrackBSLLPoly.h"
9 #include <stdexcept>
10 
11 using namespace PESA;
12 
13 namespace {
14 
15 // Ordering of the monomials in Bx, By, tx, ty powers,
16 // array indices correspond to power of those variables
17 const std::int8_t g_order[3][3][3][3] = {
18 
19 {
20  {{ 0, 1, 2 }, // Bx**0, By**0, tx**0, ty**0..ty**2
21  { 3, 4, -1 }, // Bx**0, By**0, tx**1, ty**0..ty**2
22  { 5, -1, -1 }}, // Bx**0, By**0, tx**2, ty**0..ty**2
23 
24  {{ 6, 7, -1 }, // Bx**0, By**1, tx**0, ty**0..ty**2
25  { 8, -1, -1 }, // Bx**0, By**1, tx**1, ty**0..ty**2
26  { -1, -1, -1 }}, // Bx**0, By**1, tx**2, ty**0..ty**2
27 
28  {{ 9, -1, -1 }, // Bx**0, By**2, tx**0, ty**0..ty**2
29  { -1, -1, -1 }, // Bx**0, By**2, tx**1, ty**0..ty**2
30  { -1, -1, -1 }}, // Bx**0, By**2, tx**2, ty**0..ty**2
31 },
32 
33 {
34  {{ 10, 11, -1 }, // Bx**1, By**0, tx**0, ty**0..ty**2
35  { 12, -1, -1 }, // Bx**1, By**0, tx**1, ty**0..ty**2
36  { -1, -1, -1 }}, // Bx**1, By**0, tx**2, ty**0..ty**2
37 
38  {{ 13, -1, -1 }, // Bx**1, By**1, tx**0, ty**0..ty**2
39  { -1, -1, -1 }, // Bx**1, By**1, tx**1, ty**0..ty**2
40  { -1, -1, -1 }}, // Bx**1, By**1, tx**2, ty**0..ty**2
41 
42  {{ -1, -1, -1 }, // Bx**1, By**2, tx**0, ty**0..ty**2
43  { -1, -1, -1 }, // Bx**1, By**2, tx**1, ty**0..ty**2
44  { -1, -1, -1 }}, // Bx**1, By**2, tx**2, ty**0..ty**2
45 },
46 
47 {
48  {{ 14, -1, -1 }, // Bx**2, By**0, tx**0, ty**0..ty**2
49  { -1, -1, -1 }, // Bx**2, By**0, tx**1, ty**0..ty**2
50  { -1, -1, -1 }}, // Bx**2, By**0, tx**2, ty**0..ty**2
51 
52  {{ -1, -1, -1 }, // Bx**2, By**1, tx**0, ty**0..ty**2
53  { -1, -1, -1 }, // Bx**2, By**1, tx**1, ty**0..ty**2
54  { -1, -1, -1 }}, // Bx**2, By**1, tx**2, ty**0..ty**2
55 
56  {{ -1, -1, -1 }, // Bx**2, By**2, tx**0, ty**0..ty**2
57  { -1, -1, -1 }, // Bx**2, By**2, tx**1, ty**0..ty**2
58  { -1, -1, -1 }}, // Bx**2, By**2, tx**2, ty**0..ty**2
59 }
60 };
61 
62 // number of possible combinations of all powers,
63 // 15 is for squares, 1 for log
64 const unsigned g_size = 15 + 1;
65 
66 // Ordering of the monomials in omegax, omegay powers,
67 // array indices correspond to power of those variables
68 const std::int8_t g_order2[3][3] = {
69  { 0, 1, 2 }, // Bx**0, By**0, tx**0, ty**0..ty**2
70  { 3, 4, -1 }, // Bx**0, By**0, tx**1, ty**0..ty**2
71  { 5, -1, -1 }, // Bx**0, By**0, tx**2, ty**0..ty**2
72 };
73 
74 // number of possible combinations of omega powers
75 const unsigned g_size2 = 6;
76 
77 } // namespace
78 
79 unsigned
80 T2TrackBSLLPoly::nbins()
81 {
82  // two extra bins to count number of tracks and beam_size*n_tracks
83  return
84  g_size*g_size2 // all polynomial coefficients
85  + 1 // Sum(z0)
86  + 1 // Sum(z0**2)
87  + 1 // Sum(1)
88  + 1 // Sum(beam_size)
89  ;
90 }
91 
92 int
93 T2TrackBSLLPoly::idx(unsigned power_Bx, unsigned power_By,
94  unsigned power_tx, unsigned power_ty,
95  unsigned power_omegax, unsigned power_omegay)
96 {
97  // Previous error condition returned -1, which would be used as an array index in
98  // T2TrackBSLLPoly::update and return nonsense
99  if (power_Bx > 2 or power_By > 2
100  or power_tx > 2 or power_ty > 2
101  or power_omegax > 2 or power_omegay > 2) {
102  throw std::out_of_range("Power>2: Calculated array index is out of range in T2TrackBSLLPoly::idx");
103  }
104 
105  int idx = g_order[power_Bx][power_By][power_tx][power_ty];
106  if (idx < 0) throw std::out_of_range("idx<0: Calculated array index is out of range in T2TrackBSLLPoly::idx");
107 
108  int idx2 = g_order2[power_omegax][power_omegay];
109  if (idx2 < 0) throw std::out_of_range("idx2<0: Calculated array index is out of range in T2TrackBSLLPoly::idx");
110 
111  return idx*g_size2 + idx2;
112 }
113 
114 void
115 T2TrackBSLLPoly::update(double z_0, double d_0, double phi, double var_d0, std::vector<double>& coeff)
116 {
117  if (coeff.empty()) {
118  coeff.resize(nbins(), 0.);
119  }
120 
121  double cos_phi = std::cos(phi);
122  double sin_phi = std::sin(phi);
123  double var_b = m_beam_size*m_beam_size;
124  double var_bd = var_d0 + var_b;
125 
126  // z-0 and its square
127  coeff[g_size*g_size2] += z_0;
128  coeff[g_size*g_size2+1] += z_0*z_0;
129 
130  // store number of tracks and beam_size*n_tracks in last two bins
131  coeff[g_size*g_size2+2] += 1;
132  coeff[g_size*g_size2+3] += m_beam_size;
133 
134  // this code is generated by notebook, no point in changing it here
135  // unless performance is a consideration. Example: pow(cos_phi, 2) is calculated
136  // 23 times
137  using std::pow;
138  coeff[idx(0, 0, 0, 0, 0, 0)] += -pow(d_0, 2)/(2*var_bd);
139  coeff[idx(0, 0, 0, 0, 0, 1)] += pow(cos_phi, 2)*pow(d_0, 2)/(2*pow(var_bd, 2));
140  coeff[idx(0, 0, 0, 0, 0, 2)] += -pow(cos_phi, 4)*pow(d_0, 2)/(2*pow(var_bd, 3));
141  coeff[idx(0, 0, 0, 0, 1, 0)] += pow(d_0, 2)*pow(sin_phi, 2)/(2*pow(var_bd, 2));
142  coeff[idx(0, 0, 0, 0, 1, 1)] += -pow(cos_phi, 2)*pow(d_0, 2)*pow(sin_phi, 2)/pow(var_bd, 3);
143  coeff[idx(0, 0, 0, 0, 2, 0)] += -pow(d_0, 2)*pow(sin_phi, 4)/(2*pow(var_bd, 3));
144  coeff[idx(0, 0, 0, 1, 0, 0)] += cos_phi*d_0*z_0/var_bd;
145  coeff[idx(0, 0, 0, 1, 0, 1)] += -pow(cos_phi, 3)*d_0*z_0/pow(var_bd, 2);
146  coeff[idx(0, 0, 0, 1, 0, 2)] += pow(cos_phi, 5)*d_0*z_0/pow(var_bd, 3);
147  coeff[idx(0, 0, 0, 1, 1, 0)] += -cos_phi*d_0*pow(sin_phi, 2)*z_0/pow(var_bd, 2);
148  coeff[idx(0, 0, 0, 1, 1, 1)] += 2*pow(cos_phi, 3)*d_0*pow(sin_phi, 2)*z_0/pow(var_bd, 3);
149  coeff[idx(0, 0, 0, 1, 2, 0)] += cos_phi*d_0*pow(sin_phi, 4)*z_0/pow(var_bd, 3);
150  coeff[idx(0, 0, 0, 2, 0, 0)] += -pow(cos_phi, 2)*pow(z_0, 2)/(2*var_bd);
151  coeff[idx(0, 0, 0, 2, 0, 1)] += pow(cos_phi, 4)*pow(z_0, 2)/(2*pow(var_bd, 2));
152  coeff[idx(0, 0, 0, 2, 0, 2)] += -pow(cos_phi, 6)*pow(z_0, 2)/(2*pow(var_bd, 3));
153  coeff[idx(0, 0, 0, 2, 1, 0)] += pow(cos_phi, 2)*pow(sin_phi, 2)*pow(z_0, 2)/(2*pow(var_bd, 2));
154  coeff[idx(0, 0, 0, 2, 1, 1)] += -pow(cos_phi, 4)*pow(sin_phi, 2)*pow(z_0, 2)/pow(var_bd, 3);
155  coeff[idx(0, 0, 0, 2, 2, 0)] += -pow(cos_phi, 2)*pow(sin_phi, 4)*pow(z_0, 2)/(2*pow(var_bd, 3));
156  coeff[idx(0, 0, 1, 0, 0, 0)] += -d_0*sin_phi*z_0/var_bd;
157  coeff[idx(0, 0, 1, 0, 0, 1)] += pow(cos_phi, 2)*d_0*sin_phi*z_0/pow(var_bd, 2);
158  coeff[idx(0, 0, 1, 0, 0, 2)] += -pow(cos_phi, 4)*d_0*sin_phi*z_0/pow(var_bd, 3);
159  coeff[idx(0, 0, 1, 0, 1, 0)] += d_0*pow(sin_phi, 3)*z_0/pow(var_bd, 2);
160  coeff[idx(0, 0, 1, 0, 1, 1)] += -2*pow(cos_phi, 2)*d_0*pow(sin_phi, 3)*z_0/pow(var_bd, 3);
161  coeff[idx(0, 0, 1, 0, 2, 0)] += -d_0*pow(sin_phi, 5)*z_0/pow(var_bd, 3);
162  coeff[idx(0, 0, 1, 1, 0, 0)] += cos_phi*sin_phi*pow(z_0, 2)/var_bd;
163  coeff[idx(0, 0, 1, 1, 0, 1)] += -pow(cos_phi, 3)*sin_phi*pow(z_0, 2)/pow(var_bd, 2);
164  coeff[idx(0, 0, 1, 1, 0, 2)] += pow(cos_phi, 5)*sin_phi*pow(z_0, 2)/pow(var_bd, 3);
165  coeff[idx(0, 0, 1, 1, 1, 0)] += -cos_phi*pow(sin_phi, 3)*pow(z_0, 2)/pow(var_bd, 2);
166  coeff[idx(0, 0, 1, 1, 1, 1)] += 2*pow(cos_phi, 3)*pow(sin_phi, 3)*pow(z_0, 2)/pow(var_bd, 3);
167  coeff[idx(0, 0, 1, 1, 2, 0)] += cos_phi*pow(sin_phi, 5)*pow(z_0, 2)/pow(var_bd, 3);
168  coeff[idx(0, 0, 2, 0, 0, 0)] += -pow(sin_phi, 2)*pow(z_0, 2)/(2*var_bd);
169  coeff[idx(0, 0, 2, 0, 0, 1)] += pow(cos_phi, 2)*pow(sin_phi, 2)*pow(z_0, 2)/(2*pow(var_bd, 2));
170  coeff[idx(0, 0, 2, 0, 0, 2)] += -pow(cos_phi, 4)*pow(sin_phi, 2)*pow(z_0, 2)/(2*pow(var_bd, 3));
171  coeff[idx(0, 0, 2, 0, 1, 0)] += pow(sin_phi, 4)*pow(z_0, 2)/(2*pow(var_bd, 2));
172  coeff[idx(0, 0, 2, 0, 1, 1)] += -pow(cos_phi, 2)*pow(sin_phi, 4)*pow(z_0, 2)/pow(var_bd, 3);
173  coeff[idx(0, 0, 2, 0, 2, 0)] += -pow(sin_phi, 6)*pow(z_0, 2)/(2*pow(var_bd, 3));
174  coeff[idx(0, 1, 0, 0, 0, 0)] += cos_phi*d_0/var_bd;
175  coeff[idx(0, 1, 0, 0, 0, 1)] += -pow(cos_phi, 3)*d_0/pow(var_bd, 2);
176  coeff[idx(0, 1, 0, 0, 0, 2)] += pow(cos_phi, 5)*d_0/pow(var_bd, 3);
177  coeff[idx(0, 1, 0, 0, 1, 0)] += -cos_phi*d_0*pow(sin_phi, 2)/pow(var_bd, 2);
178  coeff[idx(0, 1, 0, 0, 1, 1)] += 2*pow(cos_phi, 3)*d_0*pow(sin_phi, 2)/pow(var_bd, 3);
179  coeff[idx(0, 1, 0, 0, 2, 0)] += cos_phi*d_0*pow(sin_phi, 4)/pow(var_bd, 3);
180  coeff[idx(0, 1, 0, 1, 0, 0)] += -pow(cos_phi, 2)*z_0/var_bd;
181  coeff[idx(0, 1, 0, 1, 0, 1)] += pow(cos_phi, 4)*z_0/pow(var_bd, 2);
182  coeff[idx(0, 1, 0, 1, 0, 2)] += -pow(cos_phi, 6)*z_0/pow(var_bd, 3);
183  coeff[idx(0, 1, 0, 1, 1, 0)] += pow(cos_phi, 2)*pow(sin_phi, 2)*z_0/pow(var_bd, 2);
184  coeff[idx(0, 1, 0, 1, 1, 1)] += -2*pow(cos_phi, 4)*pow(sin_phi, 2)*z_0/pow(var_bd, 3);
185  coeff[idx(0, 1, 0, 1, 2, 0)] += -pow(cos_phi, 2)*pow(sin_phi, 4)*z_0/pow(var_bd, 3);
186  coeff[idx(0, 1, 1, 0, 0, 0)] += cos_phi*sin_phi*z_0/var_bd;
187  coeff[idx(0, 1, 1, 0, 0, 1)] += -pow(cos_phi, 3)*sin_phi*z_0/pow(var_bd, 2);
188  coeff[idx(0, 1, 1, 0, 0, 2)] += pow(cos_phi, 5)*sin_phi*z_0/pow(var_bd, 3);
189  coeff[idx(0, 1, 1, 0, 1, 0)] += -cos_phi*pow(sin_phi, 3)*z_0/pow(var_bd, 2);
190  coeff[idx(0, 1, 1, 0, 1, 1)] += 2*pow(cos_phi, 3)*pow(sin_phi, 3)*z_0/pow(var_bd, 3);
191  coeff[idx(0, 1, 1, 0, 2, 0)] += cos_phi*pow(sin_phi, 5)*z_0/pow(var_bd, 3);
192  coeff[idx(0, 2, 0, 0, 0, 0)] += -pow(cos_phi, 2)/(2*var_bd);
193  coeff[idx(0, 2, 0, 0, 0, 1)] += pow(cos_phi, 4)/(2*pow(var_bd, 2));
194  coeff[idx(0, 2, 0, 0, 0, 2)] += -pow(cos_phi, 6)/(2*pow(var_bd, 3));
195  coeff[idx(0, 2, 0, 0, 1, 0)] += pow(cos_phi, 2)*pow(sin_phi, 2)/(2*pow(var_bd, 2));
196  coeff[idx(0, 2, 0, 0, 1, 1)] += -pow(cos_phi, 4)*pow(sin_phi, 2)/pow(var_bd, 3);
197  coeff[idx(0, 2, 0, 0, 2, 0)] += -pow(cos_phi, 2)*pow(sin_phi, 4)/(2*pow(var_bd, 3));
198  coeff[idx(1, 0, 0, 0, 0, 0)] += -d_0*sin_phi/var_bd;
199  coeff[idx(1, 0, 0, 0, 0, 1)] += pow(cos_phi, 2)*d_0*sin_phi/pow(var_bd, 2);
200  coeff[idx(1, 0, 0, 0, 0, 2)] += -pow(cos_phi, 4)*d_0*sin_phi/pow(var_bd, 3);
201  coeff[idx(1, 0, 0, 0, 1, 0)] += d_0*pow(sin_phi, 3)/pow(var_bd, 2);
202  coeff[idx(1, 0, 0, 0, 1, 1)] += -2*pow(cos_phi, 2)*d_0*pow(sin_phi, 3)/pow(var_bd, 3);
203  coeff[idx(1, 0, 0, 0, 2, 0)] += -d_0*pow(sin_phi, 5)/pow(var_bd, 3);
204  coeff[idx(1, 0, 0, 1, 0, 0)] += cos_phi*sin_phi*z_0/var_bd;
205  coeff[idx(1, 0, 0, 1, 0, 1)] += -pow(cos_phi, 3)*sin_phi*z_0/pow(var_bd, 2);
206  coeff[idx(1, 0, 0, 1, 0, 2)] += pow(cos_phi, 5)*sin_phi*z_0/pow(var_bd, 3);
207  coeff[idx(1, 0, 0, 1, 1, 0)] += -cos_phi*pow(sin_phi, 3)*z_0/pow(var_bd, 2);
208  coeff[idx(1, 0, 0, 1, 1, 1)] += 2*pow(cos_phi, 3)*pow(sin_phi, 3)*z_0/pow(var_bd, 3);
209  coeff[idx(1, 0, 0, 1, 2, 0)] += cos_phi*pow(sin_phi, 5)*z_0/pow(var_bd, 3);
210  coeff[idx(1, 0, 1, 0, 0, 0)] += -pow(sin_phi, 2)*z_0/var_bd;
211  coeff[idx(1, 0, 1, 0, 0, 1)] += pow(cos_phi, 2)*pow(sin_phi, 2)*z_0/pow(var_bd, 2);
212  coeff[idx(1, 0, 1, 0, 0, 2)] += -pow(cos_phi, 4)*pow(sin_phi, 2)*z_0/pow(var_bd, 3);
213  coeff[idx(1, 0, 1, 0, 1, 0)] += pow(sin_phi, 4)*z_0/pow(var_bd, 2);
214  coeff[idx(1, 0, 1, 0, 1, 1)] += -2*pow(cos_phi, 2)*pow(sin_phi, 4)*z_0/pow(var_bd, 3);
215  coeff[idx(1, 0, 1, 0, 2, 0)] += -pow(sin_phi, 6)*z_0/pow(var_bd, 3);
216  coeff[idx(1, 1, 0, 0, 0, 0)] += cos_phi*sin_phi/var_bd;
217  coeff[idx(1, 1, 0, 0, 0, 1)] += -pow(cos_phi, 3)*sin_phi/pow(var_bd, 2);
218  coeff[idx(1, 1, 0, 0, 0, 2)] += pow(cos_phi, 5)*sin_phi/pow(var_bd, 3);
219  coeff[idx(1, 1, 0, 0, 1, 0)] += -cos_phi*pow(sin_phi, 3)/pow(var_bd, 2);
220  coeff[idx(1, 1, 0, 0, 1, 1)] += 2*pow(cos_phi, 3)*pow(sin_phi, 3)/pow(var_bd, 3);
221  coeff[idx(1, 1, 0, 0, 2, 0)] += cos_phi*pow(sin_phi, 5)/pow(var_bd, 3);
222  coeff[idx(2, 0, 0, 0, 0, 0)] += -pow(sin_phi, 2)/(2*var_bd);
223  coeff[idx(2, 0, 0, 0, 0, 1)] += pow(cos_phi, 2)*pow(sin_phi, 2)/(2*pow(var_bd, 2));
224  coeff[idx(2, 0, 0, 0, 0, 2)] += -pow(cos_phi, 4)*pow(sin_phi, 2)/(2*pow(var_bd, 3));
225  coeff[idx(2, 0, 0, 0, 1, 0)] += pow(sin_phi, 4)/(2*pow(var_bd, 2));
226  coeff[idx(2, 0, 0, 0, 1, 1)] += -pow(cos_phi, 2)*pow(sin_phi, 4)/pow(var_bd, 3);
227  coeff[idx(2, 0, 0, 0, 2, 0)] += -pow(sin_phi, 6)/(2*pow(var_bd, 3));
228  coeff[idx(0, 0, 0, 0, 0, 0)+90] += -log(var_bd)/2;
229  coeff[idx(0, 0, 0, 0, 0, 1)+90] += -pow(cos_phi, 2)/(2*var_bd);
230  coeff[idx(0, 0, 0, 0, 0, 2)+90] += pow(cos_phi, 4)/(4*pow(var_bd, 2));
231  coeff[idx(0, 0, 0, 0, 1, 0)+90] += -pow(sin_phi, 2)/(2*var_bd);
232  coeff[idx(0, 0, 0, 0, 1, 1)+90] += pow(cos_phi, 2)*pow(sin_phi, 2)/(2*pow(var_bd, 2));
233  coeff[idx(0, 0, 0, 0, 2, 0)+90] += pow(sin_phi, 4)/(4*pow(var_bd, 2));
234 }
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
T2TrackBSLLPoly.h
PESA::T2TrackBSLLPoly::update
void update(double z0, double d0, double phi0, double d0_var, std::vector< double > &coeff)
Update polynomial coefficients with track data.
Definition: T2TrackBSLLPoly.cxx:115
PESA
Local tools.
Definition: T2BeamSpot.cxx:13
PESA::T2TrackBSLLPoly::idx
static int idx(unsigned power_Bx, unsigned power_By, unsigned power_tx, unsigned power_ty, unsigned power_omegax, unsigned power_omegay)
Return bin number (0-based) for a monomial coefficient given the powers of the variables in a monomia...
Definition: T2TrackBSLLPoly.cxx:93
PESA::T2TrackBSLLPoly::nbins
static unsigned nbins()
Return number of bins in the histogram for all plynomial coefficients (and few other numbers).
Definition: T2TrackBSLLPoly.cxx:80
MuonCalib::Legendre::coeff
constexpr double coeff(unsigned int l, unsigned int k)
Calculates the n-th coefficient of the legendre polynomial series.
Definition: LegendrePoly.h:73
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
PESA::T2TrackBSLLPoly::m_beam_size
double m_beam_size
Definition: T2TrackBSLLPoly.h:88
Generated on Thu Nov 7 2024 21:26:51 for ATLAS Offline Software by   doxygen 1.8.18