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