ATLAS Offline Software
sincos_poly.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "GeoSpecialShapes/LArWheelCalculator.h"
6 
7 #include "CxxUtils/sincos.h"
8 
9 #ifdef PORTABLE_LAR_SHAPE
10 #define SINCOSPOLY_USE_EIGEN 1
11 #endif
12 
13 #ifndef SINCOSPOLY_USE_EIGEN
14 #include "TMath.h"
15 #include "TMatrixD.h"
16 #include "TVectorD.h"
17 #include "TMatrixDLazy.h"
18 #include "TDecompLU.h"
19 #include "TDecompSVD.h"
20 typedef TVectorD VECTOR;
21 typedef TMatrixD MATRIX;
22 typedef TMatrixDSym SYMMATRIX;
23 #else
24 #include "Eigen/Dense"
25 typedef Eigen::VectorXd VECTOR;
26 typedef Eigen::MatrixXd MATRIX;
27 typedef Eigen::MatrixXd SYMMATRIX;
28 #define Int_t int
29 #define Double_t double
30 #endif
31 
32 #include <sys/time.h>
33 #include <iostream>
34 #include <iomanip>
35 #include <math.h>
36 #include <mutex>
37 
38 #define DEBUGPRINT 0
39 #include "GeoModelKernel/Units.h"
40 using namespace GeoModelKernelUnits;
41 
42 #ifndef SINCOSPOLY_USE_EIGEN
43 template<typename T, typename Q>
44 std::ostream & operator << (std::ostream & ostr, const TVectorT<T> & v)
45 {
46  std::ios_base::fmtflags save_flags(ostr.flags());
47  ostr << '[';
48  ostr << std::scientific;
49  for(Int_t idx=v.GetLwb();idx<v.GetUpb();idx++) {
50  ostr << v[idx] << ", ";
51  }
52  ostr << v[v.GetUpb()];
53  ostr << ']';
54  ostr.flags(save_flags);
55  return ostr;
56 }
57 #endif
58 
59 // find best approximation of y values using linear combination of basis functions in bf
60 static VECTOR findLinearApproximation(
61  const Int_t dataLen, const Int_t nBasisFuntions,
62  const VECTOR &y, const MATRIX & bf)
63 {
64 #ifndef SINCOSPOLY_USE_EIGEN
65  SYMMATRIX A(nBasisFuntions);
66  VECTOR vY(nBasisFuntions);
67 #else
68  SYMMATRIX A=SYMMATRIX::Zero(nBasisFuntions,nBasisFuntions);
69  VECTOR vY=VECTOR::Zero(nBasisFuntions);
70 #endif
71  for(Int_t j = 0; j < nBasisFuntions; ++ j){
72  for(Int_t k = 0; k < nBasisFuntions; ++ k){
73  Double_t Ajk = 0.0;
74  for(Int_t i = 0; i < dataLen; ++ i){
75  Ajk += bf(j, i) * bf(k, i);
76  }
77  A(j, k) = Ajk;
78  }
79  }
80 
81  for(Int_t k = 0; k < nBasisFuntions; ++ k){
82  Double_t vYk = 0.0;
83  for(Int_t i = 0; i < dataLen; ++ i){
84  vYk += y[i]*bf(k,i);
85  }
86  vY[k] = vYk;
87  }
88 #ifndef SINCOSPOLY_USE_EIGEN
89  SYMMATRIX Ainv(A);
90  Ainv.Invert();
91  return Ainv*vY;
92 #else
93  return A.inverse()*vY;
94 #endif
95 }
96 
97 using namespace CLHEP;
98 
100 namespace {
101  const Int_t nBasisFunctions = LARWC_SINCOS_POLY + 1;
102 
104  struct sincos_params {
105  std::array<double, nBasisFunctions> sin{};
106  std::array<double, nBasisFunctions> cos{};
107  };
108 
110  sincos_params calc_params(LArWheelCalculator& lwc, double Rmin, double Rmax) {
111 
112  sincos_params p;
113 
114  const Double_t Rstep = 1.*mm;
115  const Int_t nrPoints = (Rmax - Rmin) * (1./Rstep);
116  const Int_t dataLen = nrPoints + 1;
117 #ifndef SINCOSPOLY_USE_EIGEN
118  VECTOR x(dataLen); // angle points
119  VECTOR ysin(dataLen); // to be approximated function values at angle points - sin
120  VECTOR ycos(dataLen); // to be approximated function values at angle points - cos
121  MATRIX bf(nBasisFunctions, dataLen); // Matrix of values of basis functions at angle points
122 #else
123  VECTOR x=VECTOR::Zero(dataLen); // angle points
124  VECTOR ysin=VECTOR::Zero(dataLen); // to be approximated function values at angle points - sin
125  VECTOR ycos=VECTOR::Zero(dataLen); // to be approximated function values at angle points - cos
126  MATRIX bf=MATRIX::Zero(nBasisFunctions, dataLen); // Matrix of values of basis functions at angle points
127 #endif
128  for(Int_t i = 0; i < dataLen; ++ i){
129  const Double_t a = Rmin + i * Rstep;
130  x[i] = a;
131  CxxUtils::sincos scalpha(lwc.parameterized_slant_angle(a));
132  ysin[i] = scalpha.sn;
133  ycos[i] = scalpha.cs;
134  for(Int_t n = 0; n < nBasisFunctions; ++ n) {
135  bf(n, i) = pow(a, n);
136  }
137  }
138 
139  VECTOR params_sin =
140  findLinearApproximation(dataLen, nBasisFunctions, ysin, bf);
141  VECTOR params_cos =
142  findLinearApproximation(dataLen, nBasisFunctions, ycos, bf);
143 
144  for(Int_t i = 0; i < nBasisFunctions; ++ i){
145  p.sin[i] = params_sin[i];
146  p.cos[i] = params_cos[i];
147  }
148 
149  //
150  // Nothing below is needed unless for debugging
151  //
152 #if DEBUGPRINT
153  std::cout << "LARWC_SINCOS_POLY: " << LARWC_SINCOS_POLY << std::endl;
154  std::cout << "sin params:" << params_sin << std::endl;
155  std::cout << "cos params:" << params_cos << std::endl;
156 
157  double dsinr = 0., dcosr = 0.;
158  double dtrigr = 0;
159 
160  double dsin = 0., dcos = 0.;
161  double dtrig = 0.;
162  for(double r = Rmin + 40.; r < Rmax - 40.; r += Rstep / 10.){
163  CxxUtils::sincos scalpha(lwc.parameterized_slant_angle(r));
164  double sin_a, cos_a;
165  double sin_a_v, cos_a_v;
166  lwc.parameterized_sincos(r, sin_a, cos_a);
167  lwc.m_vsincos_par.eval(r, sin_a_v, cos_a_v);
168  std::streamsize ss = std::cout.precision();
169  std::cout.precision(16);
170  std::cout << "def: " << r << " " << sin_a << " " << cos_a << std::endl;
171  std::cout << "vec: " << r << " " << sin_a_v << " " << cos_a_v << std::endl;
172  std::cout << "dif: " << r << " " << (sin_a - sin_a_v) / sin_a << " " << (cos_a - cos_a_v) / cos_a << std::endl;
173  std::cout.precision(ss);
174  double ds = fabs(scalpha.sn - sin_a);
175  if(ds > dsin){
176  dsin = ds;
177  dsinr = r;
178  }
179  double dc = fabs(scalpha.cs - cos_a);
180  if(dc > dcos){
181  dcos = dc;
182  dcosr = r;
183  }
184  double dt = fabs(sin_a*sin_a + cos_a*cos_a - 1.);
185  if(dt > dtrig){
186  dtrig = dt;
187  dtrigr = r;
188  }
189  }
190 
191  std::cout << "Max. difference: " << std::endl
192  << "\tsin: " << dsin << " at " << dsinr << std::endl
193  << "\tcos: " << dcos << " at " << dcosr << std::endl
194  << "\tsin^2+cos^2: " << dtrig << " at " << dtrigr << std::endl;
195 #endif
196 
197 #ifdef HARDDEBUG
198  TVectorD y_test(dataLen);
199  const Int_t nIter=10000;
200  std::cout << "Perfomance test started, " << nIter << " iterations" << std::endl;
201 
202  double y_testsin[dataLen] ={};
203  double y_testcos[dataLen] ={};
204  struct timeval tvsincos_start, tvsincos_stop;
205  gettimeofday(&tvsincos_start, 0);
206  for(Int_t iIter=0;iIter<nIter;iIter++) {
207  for(Int_t i=0;i<dataLen;i++) {
208  sincos(parameterized_slant_angle(x[i]), &y_testsin[i], &y_testcos[i]);
209  }
210  }
211  gettimeofday(&tvsincos_stop, 0);
212  double timeSinCos=(tvsincos_stop.tv_sec-tvsincos_start.tv_sec + 1E-6*(tvsincos_stop.tv_usec-tvsincos_start.tv_usec))/nIter;
213 
214  std::cout.unsetf ( std::ios::fixed | std::ios::scientific );
215  std::cout << "Time to fill 2x" << dataLen << " elements using sincos function: " << timeSinCos << std::endl;
216 
217  struct timeval tvpoly_start, tvpoly_stop;
218  gettimeofday(&tvpoly_start, 0);
219  for(Int_t iIter=0;iIter<nIter;iIter++) {
220  for(Int_t i=0;i<dataLen;i++) {
221  lwc.parameterized_sincos(x[i], y_testsin[i], y_testcos[i]);
222  }
223  }
224  gettimeofday(&tvpoly_stop, 0);
225  double timePoly=(tvpoly_stop.tv_sec - tvpoly_start.tv_sec + 1E-6*(tvpoly_stop.tv_usec - tvpoly_start.tv_usec))/nIter;
226  std::cout << "Time to fill 2x" << dataLen << " elements using approximation sin&cos: " << timePoly << std::endl;
227  std::cout.unsetf ( std::ios::fixed | std::ios::scientific );
228  std::cout << "Approximation is " << timeSinCos/timePoly << " faster " << std::endl;
229 #endif
230 
231  return p;
232  }
233 
235  const sincos_params& inner_params(LArWheelCalculator& lwc) {
236  static const sincos_params p = calc_params(lwc, 250.*mm, 750.*mm);
237  return p;
238  }
239 
241  const sincos_params& outer_params(LArWheelCalculator& lwc) {
242  static const sincos_params p = calc_params(lwc, 560.*mm, 2090.*mm);
243  return p;
244  }
245 }
246 
247 void LArWheelCalculator::fill_sincos_parameterization()
248 {
249  // The polynomial approximations are calculated once per side, and stored
250  // as statics for reuse in successive calculator instances.
251  const sincos_params& p = m_isInner ? inner_params(*this) : outer_params(*this);
252 
253  // Fill the parameters into our instances variables
254  m_sin_parametrization = p.sin;
255  m_cos_parametrization = p.cos;
256 
257  // Parameterization for the vectorized sincos calculation
258  // see ATLASSIM-4753 for details
259  // s4, s5, c4, c5
260  // s2, s3, c2, c3
261  // s0, s1, c0, c1
262  m_vsincos_par.param_0[0] = m_sin_parametrization[4];
263  m_vsincos_par.param_0[1] = m_sin_parametrization[5];
264  m_vsincos_par.param_0[2] = m_cos_parametrization[4];
265  m_vsincos_par.param_0[3] = m_cos_parametrization[5];
266  m_vsincos_par.param_1[0] = m_sin_parametrization[2];
267  m_vsincos_par.param_1[1] = m_sin_parametrization[3];
268  m_vsincos_par.param_1[2] = m_cos_parametrization[2];
269  m_vsincos_par.param_1[3] = m_cos_parametrization[3];
270  m_vsincos_par.param_2[0] = m_sin_parametrization[0];
271  m_vsincos_par.param_2[1] = m_sin_parametrization[1];
272  m_vsincos_par.param_2[2] = m_cos_parametrization[0];
273  m_vsincos_par.param_2[3] = m_cos_parametrization[1];
274 }
beamspotman.r
def r
Definition: beamspotman.py:676
checkxAOD.ds
ds
Definition: Tools/PyUtils/bin/checkxAOD.py:260
PowhegControl_ttHplus_NLO.ss
ss
Definition: PowhegControl_ttHplus_NLO.py:83
CaloSwCorrections.lwc
def lwc(flags, cells_name, *args, **kw)
Definition: CaloSwCorrections.py:215
sincos.h
Helper to simultaneously calculate sin and cos of the same angle.
MATRIX
TMatrixD MATRIX
Definition: sincos_poly.cxx:21
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
x
#define x
dqt_zlumi_alleff_HIST.A
A
Definition: dqt_zlumi_alleff_HIST.py:110
A
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
lumiFormat.i
int i
Definition: lumiFormat.py:85
beamspotman.n
n
Definition: beamspotman.py:731
CaloNoise_fillDB.dt
dt
Definition: CaloNoise_fillDB.py:58
CLHEP
STD'S.
Definition: IAtRndmGenSvc.h:19
LARWC_SINCOS_POLY
#define LARWC_SINCOS_POLY
Definition: LArWheelCalculator.h:35
LArWheelCalculator
Definition: LArWheelCalculator.h:60
SYMMATRIX
TMatrixDSym SYMMATRIX
Definition: sincos_poly.cxx:22
LArWheelCalculator::fill_sincos_parameterization
void fill_sincos_parameterization()
Definition: sincos_poly.cxx:247
VP1PartSpect::E
@ E
Definition: VP1PartSpectFlags.h:21
python.SystemOfUnits.mm
int mm
Definition: SystemOfUnits.py:83
python.PyAthena.v
v
Definition: PyAthena.py:154
a
TList * a
Definition: liststreamerinfos.cxx:10
y
#define y
VECTOR
TVectorD VECTOR
Definition: sincos_poly.cxx:20
LArNewCalib_DelayDump_OFC_Cali.idx
idx
Definition: LArNewCalib_DelayDump_OFC_Cali.py:69
CxxUtils::sincos
Helper to simultaneously calculate sin and cos of the same angle.
Definition: sincos.h:76
Rmin
double Rmin
Definition: LArDetectorConstructionTBEC.cxx:56
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
operator<<
std::ostream & operator<<(std::ostream &ostr, const TVectorT< T > &v)
Definition: sincos_poly.cxx:44
fitman.k
k
Definition: fitman.py:528
generate::Zero
void Zero(TH1D *hin)
Definition: generate.cxx:32
LArWheelCalculator.h
Generated on Sun Dec 22 2024 21:18:14 for ATLAS Offline Software by   doxygen 1.8.18