ATLAS Offline Software
Loading...
Searching...
No Matches
DistanceCalculatorSaggingOff.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
3*/
4
7
8//#define LWC_PARAM_ANGLE
9
10#ifdef LWC_PARAM_ANGLE
11#include <CxxUtils/sincos.h>
12#endif
13
14#include <signal.h>
15#include <cassert>
16#include <cmath>
17
19{
20
26
27 double DistanceCalculatorSaggingOff::DistanceToTheNeutralFibre_ref(const CLHEP::Hep3Vector& P, int /*fan_number*/) const
28 {
29 assert(P.y() > 0.);
30 double distance = 0.;
31 double z = P.z() - lwc()->m_StraightStartSection;
32 double x = P.x();
33
34#ifdef LWC_PARAM_ANGLE //old variant
35 const double alpha = lwc()->parameterized_slant_angle(P.y());
36 //double cos_a, sin_a;
37 //::sincos(alpha, &sin_a, &cos_a);
38 CxxUtils::sincos scalpha(alpha);
39 const double cos_a = scalpha.cs, sin_a = scalpha.sn;
40#else // parameterized sine
41 double cos_a, sin_a;
42 lwc()->m_vsincos_par.eval(P.y(), sin_a, cos_a);
43#endif
44 // determination of the nearest quarter-wave number
45 int nqwave = (z < 0.) ? 0 : int(z / lwc()->m_QuarterWaveLength);
46 //if(z < 0.) nqwave = 0;
47 //else nqwave = int(z / lwc()->m_QuarterWaveLength);
48
49 bool begin_qw = false;
50 if((nqwave % 2) != 0){
51 nqwave ++;
52 begin_qw = true;
53 }
54
55 nqwave /= 2;
56
57 // now nqwave is not the number of quarter-wave but the number of half-wave
58 // half-waves with last and zero numbers are not really half-waves but start
59 // and finish quarter-waves
60 // It means that half-wave number 1 begins after starting quarter-wave
61 if(nqwave != 0 && nqwave != lwc()->m_NumberOfHalfWaves){ // regular half-waves
62 z -= nqwave * lwc()->m_HalfWaveLength;
63 // there are some symmetries, so use them
64 if((nqwave % 2) == 0) x = -x;
65 if(begin_qw){
66 z = -z;
67 x = -x;
68 }
69 // certain situation: rising slope of wave, z is positive
70 // rotate to prime-coordinate system (see description)
71 const double z_prime = z * cos_a + x * sin_a;
72 const double x_prime = x * cos_a - z * sin_a;
73 const double straight_part = (lwc()->m_QuarterWaveLength - lwc()->m_FanFoldRadius * sin_a) / cos_a;
74 if(z_prime > straight_part){// fold region
75 const double dz = straight_part - z_prime;
76 const double dx = x_prime + lwc()->m_FanFoldRadius;
77 distance = sqrt(dz*dz + dx*dx) - lwc()->m_FanFoldRadius;
78 } else if(z_prime > -straight_part){
79 distance = x_prime; // straight part of the quarter-wave
80 } else {// fold region
81 const double dz = straight_part + z_prime;
82 const double dx = x_prime - lwc()->m_FanFoldRadius;
83 distance = lwc()->m_FanFoldRadius - sqrt(dz*dz + dx*dx);
84 }
85 // set correct sign for result
86 if(!begin_qw) distance = -distance;
87 if((nqwave % 2) == 0) distance = -distance;
88
89 } else { // start and finish quarter-waves
90 if(nqwave == 0) { // start quarter-wave
91 x = - x;
92 } else { // finish quarter-wave
93 z = lwc()->m_ActiveLength - z;
94 }
95
96 const double tan_beta = sin_a/(1.0 + cos_a); //tan(alpha * 0.5);
97 const double local_straight_section = lwc()->m_FanFoldRadius * tan_beta;
98 if( z < - local_straight_section &&
99 ( x < lwc()->m_FanFoldRadius ||
100 x < - lwc()->m_StraightStartSection * z / local_straight_section / tan_beta ) )
101 {
102 distance = - x;
103 }
104 else {
105 const double z_prime = z * cos_a + x * sin_a;
106 const double x_prime = x * cos_a - z * sin_a;
107 if (z_prime < local_straight_section) { // start fold region
108 const double dz = local_straight_section - z_prime;
109 const double dx = x_prime - lwc()->m_FanFoldRadius;
110 distance = sqrt(dz*dz + dx*dx) - lwc()->m_FanFoldRadius;
111 } else {
112 const double straight_part = (lwc()->m_QuarterWaveLength - lwc()->m_FanFoldRadius * sin_a) / cos_a;
113 if (z_prime <= straight_part) { // straight part of quarter-wave
114 distance = - x_prime;
115 } else { // regular fold region of the quarter-wave
116 const double dz = straight_part - z_prime;
117 const double dx = x_prime + lwc()->m_FanFoldRadius;
118 distance = lwc()->m_FanFoldRadius - sqrt(dz*dz + dx*dx);
119 }
120 }
121 }
122 // set correct sign
123 if (nqwave == 0) distance = -distance;
124 }
125#ifdef HARDDEBUG
126 double dd = DistanceToTheNeutralFibre_ref(P);
127 if(std::abs(dd - distance) > 0.000001){
128 //static int cnt = 0;
129 std::cout << "DTNF MISMATCH " << this << " " << P << " "
130 << dd << " vs " << distance << std::endl;
131 //cnt ++;
132 //if(cnt > 100) exit(0);
133 }
134#endif
135 return distance;
136 }
137
138 // IMPROVED PERFORMANCE
139 double DistanceCalculatorSaggingOff::DistanceToTheNeutralFibre(const CLHEP::Hep3Vector& P, int /*fan_number*/) const
140 {
141 double z = P.z() - lwc()->m_StraightStartSection;
142 double x = P.x();
143
144#ifdef LWC_PARAM_ANGLE //old variant
145 const double alpha = lwc()->parameterized_slant_angle(P.y());
146 CxxUtils::sincos scalpha(alpha);
147 double cos_a = scalpha.cs, sin_a = scalpha.sn;
148#else // parameterized sine
149 double cos_a, sin_a;
150 lwc()->m_vsincos_par.eval(P.y(), sin_a, cos_a);
151#endif
152
153 bool sqw = false;
154 if(z > lwc()->m_QuarterWaveLength){
155 if(z < m_EndQuarterWave){ // regular half-waves
156 unsigned int nhwave = (unsigned int)(z / lwc()->m_HalfWaveLength + 0.5);
157 z -= lwc()->m_HalfWaveLength * nhwave;
158 const double straight_part = (lwc()->m_QuarterWaveLength - lwc()->m_FanFoldRadius * sin_a) / cos_a;
159 nhwave &= 1U;
160 if(nhwave == 0) sin_a = - sin_a;
161 double z_prime = z * cos_a + x * sin_a;
162 const double x_prime = z * sin_a - x * cos_a;
163 if(z_prime > straight_part){ // up fold region
164 const double dz = z_prime - straight_part;
165 if(nhwave == 0){
166 const double dx = x_prime + lwc()->m_FanFoldRadius;
167 return sqrt(dz*dz + dx*dx) - lwc()->m_FanFoldRadius;
168 } else {
169 const double dx = x_prime - lwc()->m_FanFoldRadius;
170 return lwc()->m_FanFoldRadius - sqrt(dz*dz + dx*dx);
171 }
172 }
173 z_prime += straight_part;
174 if(z_prime > 0){
175 return x_prime; // straight part of the quarter-wave
176 } else { // low fold region
177 const double &dz = z_prime;
178 if(nhwave == 0){
179 const double dx = x_prime - lwc()->m_FanFoldRadius;
180 return lwc()->m_FanFoldRadius - sqrt(dz*dz + dx*dx);
181 } else {
182 const double dx = x_prime + lwc()->m_FanFoldRadius;
183 return sqrt(dz*dz + dx*dx) - lwc()->m_FanFoldRadius;
184 }
185 }
186 } else { // ending quarter-wave
187 z = lwc()->m_ActiveLength - z;
188 }
189 } else { // starting quarter-wave
190 x = - x;
191 sqw = true;
192 }
193
194 // start and finish quarter-waves
195 const double tan_beta = sin_a/(1.0 + cos_a); //tan(alpha * 0.5);
196 const double local_straight_section = lwc()->m_FanFoldRadius * tan_beta;
197 if(z < - local_straight_section &&
198 ( x < lwc()->m_FanFoldRadius ||
199 x < - lwc()->m_StraightStartSection * z / local_straight_section / tan_beta ))
200 {
201 return sqw? x: (-x);
202 }
203 else {
204 const double z_prime = z * cos_a + x * sin_a;
205 const double x_prime = x * cos_a - z * sin_a;
206 if (z_prime < local_straight_section) { // start fold region
207 const double dz = local_straight_section - z_prime;
208 const double dx = x_prime - lwc()->m_FanFoldRadius;
209 if(sqw) return lwc()->m_FanFoldRadius - sqrt(dz*dz + dx*dx);
210 else return sqrt(dz*dz + dx*dx) - lwc()->m_FanFoldRadius;
211 } else {
212 const double straight_part =
213 (lwc()->m_QuarterWaveLength - lwc()->m_FanFoldRadius * sin_a) / cos_a;
214 if (z_prime <= straight_part) { // straight part of quarter-wave
215 return sqw? x_prime: (-x_prime);
216 } else { // regular fold region of the quarter-wave
217 const double dz = straight_part - z_prime;
218 const double dx = x_prime + lwc()->m_FanFoldRadius;
219 if(sqw) return sqrt(dz*dz + dx*dx) - lwc()->m_FanFoldRadius;
220 else return lwc()->m_FanFoldRadius - sqrt(dz*dz + dx*dx);
221 }
222 }
223 }
224 }
225
226 CLHEP::Hep3Vector DistanceCalculatorSaggingOff::NearestPointOnNeutralFibre(const CLHEP::Hep3Vector &P, int /*fan_number*/) const
227 {
228 CLHEP::Hep3Vector result;
229 double z = P.z() - lwc()->m_StraightStartSection;
230 double x = P.x();
231 double y = P.y();
232
233#ifdef LWC_PARAM_ANGLE //old variant
234 const double alpha = lwc()->parameterized_slant_angle(P.y());
235 CxxUtils::sincos scalpha(alpha);
236 const double cos_a = scalpha.cs, sin_a = scalpha.sn;
237#else // parameterized sine
238 double cos_a, sin_a;
239 lwc()->m_vsincos_par.eval(P.y(), sin_a, cos_a);
240#endif
241
242 int nqwave;
243 if(z < 0.) nqwave = 0;
244 else nqwave = int(z / lwc()->m_QuarterWaveLength);
245 bool begin_qw = false;
246 if((nqwave % 2) != 0){
247 nqwave ++;
248 begin_qw = true;
249 }
250 nqwave /= 2;
251 if(nqwave != 0 && nqwave != lwc()->m_NumberOfHalfWaves){
252 z -= nqwave * lwc()->m_HalfWaveLength;
253 if((nqwave % 2) == 0) x = -x;
254 if(begin_qw){
255 z = -z;
256 x = -x;
257 }
258 const double z_prime = z * cos_a + x * sin_a;
259 const double x_prime = x * cos_a - z * sin_a;
260 const double straight_part = (lwc()->m_QuarterWaveLength - lwc()->m_FanFoldRadius * sin_a) / cos_a;
261 const double dz = straight_part - z_prime;
262 if (dz > 0) result.set(0., y, z_prime);
263 else {
264 double a = std::atan(std::abs(dz / (x_prime + lwc()->m_FanFoldRadius)));
265 result.set(lwc()->m_FanFoldRadius * (std::cos(a) - 1), y, straight_part + lwc()->m_FanFoldRadius * std::sin(a));
266 }
267 result.rotateY(std::asin(sin_a));
268 if(begin_qw){
269 result.setX(-result.x());
270 result.setZ(-result.z());
271 }
272 if((nqwave % 2) == 0) result.setX(-result.x());
273 result.setZ(result.z() + nqwave * lwc()->m_HalfWaveLength);
274 } else {
275 if(nqwave == 0) x = -x;
276 else z = lwc()->m_ActiveLength - z;
277 const double tan_beta = sin_a/(1.0+cos_a); //tan(alpha * 0.5);
278 const double local_straight_section = lwc()->m_FanFoldRadius * tan_beta;
279 if(z < - local_straight_section &&
280 ( x < lwc()->m_FanFoldRadius ||
281 x < - lwc()->m_StraightStartSection * z / local_straight_section / tan_beta))
282 {
283 result.set(0., y, z);
284 }
285 else {
286 const double z_prime = z * cos_a + x * sin_a;
287 const double x_prime = x * cos_a - z * sin_a;
288 if(z_prime < local_straight_section) {
289 double a = std::abs(std::atan((local_straight_section - z_prime) / (x_prime - lwc()->m_FanFoldRadius)));
290
291 result.set(lwc()->m_FanFoldRadius * (1 - std::cos(a)), y, local_straight_section - lwc()->m_FanFoldRadius * std::sin(a));
292 } else {
293 double straight_part = (lwc()->m_QuarterWaveLength - lwc()->m_FanFoldRadius * sin_a) / cos_a;
294 if(z_prime <= straight_part) {
295 result.set(0., y, z_prime);
296 } else {
297 double a = std::abs(std::atan((straight_part - z_prime) / (x_prime + lwc()->m_FanFoldRadius)) );
298 result.set(lwc()->m_FanFoldRadius * (std::cos(a) - 1), y, straight_part + lwc()->m_FanFoldRadius * std::sin(a));
299 }
300 }
301 result.rotateY(std::asin(sin_a));
302 }
303 if(nqwave != 0){
304 result.setZ(lwc()->m_ActiveLength - result.z());
305 } else {
306 result.setX(-result.x());
307 }
308 }
309 result.setZ(result.z() + lwc()->m_StraightStartSection);
310 return result;
311 }
312
313 // IMPROVED VERSION
314 CLHEP::Hep3Vector DistanceCalculatorSaggingOff::NearestPointOnNeutralFibre_ref(const CLHEP::Hep3Vector &P, int /*fan_number*/) const
315 {
316 CLHEP::Hep3Vector result;
317 double z = P.z() - lwc()->m_StraightStartSection;
318 double x = P.x();
319 double y = P.y();
320
321#ifdef LWC_PARAM_ANGLE //old variant
322 const double alpha = lwc()->parameterized_slant_angle(P.y());
323 CxxUtils::sincos scalpha(alpha);
324 double cos_a = scalpha.cs, sin_a = scalpha.sn;
325#else // parameterized sine
326 double cos_a, sin_a;
327 lwc()->m_vsincos_par.eval(P.y(), sin_a, cos_a);
328#endif
329
330 bool sqw = false;
331 if(z > lwc()->m_QuarterWaveLength){
332 if(z < m_EndQuarterWave){ // regular half-waves
333 unsigned int nhwave = (unsigned int)(z / lwc()->m_HalfWaveLength + 0.5);
334 const double zshift = lwc()->m_HalfWaveLength * nhwave;
335 z -= zshift;
336 const double straight_part =
337 (lwc()->m_QuarterWaveLength - lwc()->m_FanFoldRadius * sin_a) / cos_a;
338 nhwave &= 1U;
339 if(nhwave == 0) sin_a = - sin_a;
340 const double z_prime = z * cos_a + x * sin_a;
341 if(z_prime > straight_part){ // up fold
342 const double x_prime = x * cos_a - z * sin_a;
343 const double dz = straight_part - z_prime;
344 double a1 = std::atan(std::abs(dz / (x_prime + lwc()->m_FanFoldRadius)));
345 const double x1 = lwc()->m_FanFoldRadius * (std::cos(a1) - 1.);
346 const double z1 = straight_part + lwc()->m_FanFoldRadius * std::sin(a1);
347 result.set(x1*cos_a - z1*sin_a, y, z1*cos_a + z1*sin_a);
348 return result;
349 } else if(z_prime > -straight_part){ // straight part
350 result.set(-z_prime * sin_a, y, z_prime*cos_a + zshift);
351 return result;
352 } else { // low fold
353 const double x_prime = x * cos_a - z * sin_a;
354 const double dz = straight_part + z_prime;
355 double a1 = std::atan(std::abs(dz / (x_prime + lwc()->m_FanFoldRadius)));
356 const double x1 = lwc()->m_FanFoldRadius * (std::cos(a1) - 1.);
357 const double z1 = straight_part + lwc()->m_FanFoldRadius * std::sin(a1);
358 result.set(x1*cos_a - z1*sin_a, y, z1*cos_a + z1*sin_a);
359 return result;
360 }
361 } else { // ending quarter-wave
362 z = lwc()->m_ActiveLength - z;
363 }
364 } else { // starting quarter-wave
365 x = - x;
366 sqw = true;
367 }
368
369 // start and finish quarter-waves
370 const double tan_beta = sin_a / (1.0 + cos_a);
371 const double local_straight_section = lwc()->m_FanFoldRadius * tan_beta;
372 if(z < - local_straight_section &&
373 (x < lwc()->m_FanFoldRadius ||
374 x < - lwc()->m_StraightStartSection * z / local_straight_section / tan_beta))
375 {
376 result.set(0., y, z);
377 }
378 else {
379 const double z_prime = z * cos_a + x * sin_a;
380 const double x_prime = x * cos_a - z * sin_a;
381 if(z_prime < local_straight_section) {
382 double a = std::abs(std::atan((local_straight_section - z_prime) / (x_prime - lwc()->m_FanFoldRadius)));
383 result.set(lwc()->m_FanFoldRadius * (1 - std::cos(a)), y, local_straight_section - lwc()->m_FanFoldRadius * std::sin(a));
384 } else {
385 double straight_part = (lwc()->m_QuarterWaveLength - lwc()->m_FanFoldRadius * sin_a) / cos_a;
386 if(z_prime <= straight_part) {
387 result.set(0., y, z_prime);
388 } else {
389 double a = std::abs(std::atan((straight_part - z_prime) / (x_prime + lwc()->m_FanFoldRadius)) );
390 result.set(lwc()->m_FanFoldRadius * (std::cos(a) - 1), y, straight_part + lwc()->m_FanFoldRadius * std::sin(a));
391 }
392 }
393 result.rotateY(std::asin(sin_a));
394 }
395 if(sqw) result.setX(-result.x());
396 else result.setZ(lwc()->m_ActiveLength - result.z());
397 result.setZ(result.z() + lwc()->m_StraightStartSection);
398 return result;
399 }
400
401 /*
402 input is in local fan's coordinate system
403 side: < 0 - lower phi
404 > 0 - greater phi
405 = 0 - neutral fibre
406 */
407 double DistanceCalculatorSaggingOff::AmplitudeOfSurface(const CLHEP::Hep3Vector& P, int side, int /*fan_number*/) const
408 {
409 double result = 0.;
410 double rho = lwc()->m_FanFoldRadius;
411 double z = P.z() - lwc()->m_StraightStartSection;
412
413#ifdef LWC_PARAM_ANGLE //old variant
414 const double alpha = lwc()->parameterized_slant_angle(P.y());
415 //double cos_a, sin_a;
416 //::sincos(alpha, &sin_a, &cos_a);
417 CxxUtils::sincos scalpha(alpha);
418 const double cos_a = scalpha.cs, sin_a = scalpha.sn;
419 // parameterized sine
420#else
421 double cos_a, sin_a;
422 lwc()->m_vsincos_par.eval(P.y(), sin_a, cos_a);
423#endif
424
425 // determination of the nearest quarter-wave number
426 int nqwave;
427 if(z < 0.) nqwave = 0;
428 else nqwave = int(z / lwc()->m_QuarterWaveLength);
429 bool begin_qw = false;
430 if((nqwave % 2) != 0){
431 nqwave ++;
432 begin_qw = true;
433 }
434 nqwave /= 2;
435 // now nqwave is not the number of quarter-wave but the number of half-wave
436 // half-waves with last and zero numbers are not really half-waves but start
437 // and finish quarter-waves
438 // It means that half-wave number 1 begins after starting quarter-wave
439 if(nqwave != 0 && nqwave != lwc()->m_NumberOfHalfWaves){ // regular half-waves
440 z -= nqwave * lwc()->m_HalfWaveLength;
441 if(begin_qw) z = -z;
442 double dz = lwc()->m_QuarterWaveLength - z;
443
444 int local_side = 1;
445 if((nqwave % 2) == 0){
446 if(begin_qw) local_side = -1;
447 } else {
448 if(!begin_qw) local_side = -1;
449 }
450
451 rho += lwc()->m_FanHalfThickness * local_side * side;
452
453 if(dz >= rho * sin_a){
454 result = z * sin_a / cos_a; // straight part of the quarter-wave
455 } else { // fold region
456 result = (lwc()->m_QuarterWaveLength * sin_a - rho) / cos_a
457 + sqrt(rho * rho - dz * dz);
458 }
459 result *= -local_side;
460 if(side < 0) result += lwc()->m_FanHalfThickness / cos_a;
461 else if(side > 0) result -= lwc()->m_FanHalfThickness / cos_a;
462
463 } else { // start and finish quarter-waves
464 int local_side = 1;
465 if(nqwave == 0) { // start quarter-wave
466 local_side = -1;
467 } else { // finish quarter-wave
468 z = lwc()->m_ActiveLength - z;
469 }
470
471 const double rho1i = lwc()->m_FanFoldRadius;
472 const double tan_beta = sin_a/(1.0+cos_a); //tan(alpha * 0.5);
473 const double min_local_fold_region = rho1i * tan_beta;
474
475 if(z <= - min_local_fold_region){
476 result = - side * lwc()->m_FanHalfThickness;
477 } else {
478 const double rho1 = rho1i + lwc()->m_FanHalfThickness * side * local_side;
479
480 const double max_local_fold_region = rho1 * sin_a - min_local_fold_region;
481 //const double max_local_fold_region = rho1 * tan_beta * (1. + cos_a) - min_local_fold_region;
482 if(z < max_local_fold_region){ // start fold region
483 z += min_local_fold_region;
484 result = rho1 - sqrt(rho1 * rho1 - z * z);
485 if(nqwave == 0) result = -result;
486 if(side < 0) result += lwc()->m_FanHalfThickness;
487 else if(side > 0) result -= lwc()->m_FanHalfThickness;
488 } else {
489 rho -= lwc()->m_FanHalfThickness * local_side * side;
490 const double dz = lwc()->m_QuarterWaveLength - z;
491 if(dz >= rho * sin_a){
492 result = z * sin_a / cos_a; // straight part of the quarter-wave
493 } else { // regular fold region
494 result = (lwc()->m_QuarterWaveLength * sin_a - rho) / cos_a
495 + sqrt(rho * rho - dz * dz);
496 }
497 if(nqwave == 0) result = -result;
498 if(side < 0) result += lwc()->m_FanHalfThickness / cos_a;
499 else if(side > 0) result -= lwc()->m_FanHalfThickness / cos_a;
500 }
501 }
502 }
503 return result;
504 }
505
506}
static Double_t a
static Double_t P(Double_t *tt, Double_t *par)
#define y
#define x
#define z
virtual CLHEP::Hep3Vector NearestPointOnNeutralFibre_ref(const CLHEP::Hep3Vector &p, int fan_number) const
virtual double DistanceToTheNeutralFibre(const CLHEP::Hep3Vector &p, int fan_number) const
virtual CLHEP::Hep3Vector NearestPointOnNeutralFibre(const CLHEP::Hep3Vector &p, int fan_number) const
const LArWheelCalculator * lwc() const
Return the calculator:
virtual double AmplitudeOfSurface(const CLHEP::Hep3Vector &P, int side, int fan_number) const
virtual double DistanceToTheNeutralFibre_ref(const CLHEP::Hep3Vector &p, int fan_number) const
This class separates some of the geometry details of the LAr endcap.
double parameterized_slant_angle(double) const
Calculates wave slant angle using parametrization for current wheel for given distance from calorimet...
Helper to simultaneously calculate sin and cos of the same angle.
Helper to simultaneously calculate sin and cos of the same angle.
Definition sincos.h:39
void eval(const double r, double &ATH_RESTRICT sin_a, double &ATH_RESTRICT cos_a) const ATH_RESTRICT