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