Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4QSS3.hh
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// G4QSS3
27//
28// G4QSS3 simulator
29
30// Authors: Lucio Santi, Rodrigo Castro (Univ. Buenos Aires), 2018-2021
31// --------------------------------------------------------------------
32#ifndef G4QSS3_HH
33#define G4QSS3_HH
34
35#include "G4Types.hh"
36#include "G4qss_misc.hh"
37
38#include <cmath>
39
40/**
41 * @brief G4QSS3 defines the QSS3 simulator engine used in QSS field stepper.
42 */
43
44class G4QSS3
45{
46 public:
47
48 G4QSS3(QSS_simulator);
49
50 inline QSS_simulator getSimulator() const { return this->simulator; }
51
52 inline G4int order() const { return 3; }
53
54 inline void full_definition(G4double coeff)
55 {
56 G4double* const x = simulator->q;
57 G4double* const dx = simulator->x;
58 G4double* const alg = simulator->alg;
59
60 dx[1] = x[12];
61 dx[2] = 0;
62 dx[3] = 0;
63
64 dx[5] = x[16];
65 dx[6] = 0;
66 dx[7] = 0;
67
68 dx[9] = x[20];
69 dx[10] = 0;
70 dx[11] = 0;
71
72 dx[13] = coeff * (alg[2] * x[16] - alg[1] * x[20]);
73 dx[14] = 0;
74 dx[15] = 0;
75
76 dx[17] = coeff * (alg[0] * x[20] - alg[2] * x[12]);
77 dx[18] = 0;
78 dx[19] = 0;
79
80 dx[21] = coeff * (alg[1] * x[12] - alg[0] * x[16]);
81 dx[22] = 0;
82 dx[23] = 0;
83 }
84
85 inline void dependencies(G4int i, G4double coeff)
86 {
87 G4double* const x = simulator->q;
88 G4double* const der = simulator->x;
89 G4double* const alg = simulator->alg;
90
91 switch (i)
92 {
93 case 0:
94 der[13] = coeff * (alg[2] * x[16] - alg[1] * x[20]);
95 der[14] = ((alg[2] * x[17] - x[21] * alg[1]) * coeff) / 2;
96 der[15] = (coeff * (alg[2] * x[18] - x[22] * alg[1])) / 3;
97
98 der[17] = coeff * (alg[0] * x[20] - alg[2] * x[12]);
99 der[18] = ((alg[0] * x[21] - alg[2] * x[13]) * coeff) / 2;
100 der[19] = (coeff * (alg[0] * x[22] - alg[2] * x[14])) / 3;
101
102 der[21] = coeff * (alg[1] * x[12] - alg[0] * x[16]);
103 der[22] = (coeff * (x[13] * alg[1] - alg[0] * x[17])) / 2;
104 der[23] = (coeff * (alg[1] * x[14] - x[18] * alg[0])) / 3;
105 return;
106 case 1:
107 der[13] = coeff * (alg[2] * x[16] - alg[1] * x[20]);
108 der[14] = ((alg[2] * x[17] - x[21] * alg[1]) * coeff) / 2;
109 der[15] = (coeff * (alg[2] * x[18] - x[22] * alg[1])) / 3;
110
111 der[17] = coeff * (alg[0] * x[20] - alg[2] * x[12]);
112 der[18] = ((alg[0] * x[21] - alg[2] * x[13]) * coeff) / 2;
113 der[19] = (coeff * (alg[0] * x[22] - alg[2] * x[14])) / 3;
114
115 der[21] = coeff * (alg[1] * x[12] - alg[0] * x[16]);
116 der[22] = (coeff * (x[13] * alg[1] - alg[0] * x[17])) / 2;
117 der[23] = (coeff * (alg[1] * x[14] - x[18] * alg[0])) / 3;
118 return;
119 case 2:
120 der[13] = coeff * (alg[2] * x[16] - alg[1] * x[20]);
121 der[14] = ((alg[2] * x[17] - x[21] * alg[1]) * coeff) / 2;
122 der[15] = (coeff * (alg[2] * x[18] - x[22] * alg[1])) / 3;
123
124 der[17] = coeff * (alg[0] * x[20] - alg[2] * x[12]);
125 der[18] = ((alg[0] * x[21] - alg[2] * x[13]) * coeff) / 2;
126 der[19] = (coeff * (alg[0] * x[22] - alg[2] * x[14])) / 3;
127
128 der[21] = coeff * (alg[1] * x[12] - alg[0] * x[16]);
129 der[22] = (coeff * (x[13] * alg[1] - alg[0] * x[17])) / 2;
130 der[23] = (coeff * (alg[1] * x[14] - x[18] * alg[0])) / 3;
131 return;
132 case 3:
133 der[1] = x[12];
134 der[2] = x[13] / 2;
135 der[3] = x[14] / 3;
136
137 der[17] = coeff * (alg[0] * x[20] - alg[2] * x[12]);
138 der[18] = ((alg[0] * x[21] - alg[2] * x[13]) * coeff) / 2;
139 der[19] = (coeff * (alg[0] * x[22] - alg[2] * x[14])) / 3;
140
141 der[21] = coeff * (alg[1] * x[12] - alg[0] * x[16]);
142 der[22] = (coeff * (x[13] * alg[1] - alg[0] * x[17])) / 2;
143 der[23] = (coeff * (alg[1] * x[14] - x[18] * alg[0])) / 3;
144 return;
145 case 4:
146 der[5] = x[16];
147 der[6] = x[17] / 2;
148 der[7] = x[18] / 3;
149
150 der[13] = coeff * (alg[2] * x[16] - alg[1] * x[20]);
151 der[14] = ((alg[2] * x[17] - x[21] * alg[1]) * coeff) / 2;
152 der[15] = (coeff * (alg[2] * x[18] - x[22] * alg[1])) / 3;
153
154 der[21] = coeff * (alg[1] * x[12] - alg[0] * x[16]);
155 der[22] = (coeff * (x[13] * alg[1] - alg[0] * x[17])) / 2;
156 der[23] = (coeff * (alg[1] * x[14] - x[18] * alg[0])) / 3;
157 return;
158 case 5:
159 der[9] = x[20];
160 der[10] = x[21] / 2;
161 der[11] = x[22] / 3;
162
163 der[13] = coeff * (alg[2] * x[16] - alg[1] * x[20]);
164 der[14] = ((alg[2] * x[17] - x[21] * alg[1]) * coeff) / 2;
165 der[15] = (coeff * (alg[2] * x[18] - x[22] * alg[1])) / 3;
166
167 der[17] = coeff * (alg[0] * x[20] - alg[2] * x[12]);
168 der[18] = ((alg[0] * x[21] - alg[2] * x[13]) * coeff) / 2;
169 der[19] = (coeff * (alg[0] * x[22] - alg[2] * x[14])) / 3;
170 return;
171 }
172 }
173
174 void recompute_next_times(G4int* inf, G4double t); // __attribute__((hot));
175
176 inline void recompute_all_state_times(G4double t)
177 {
178 G4double mpr;
179 G4double* const x = simulator->x;
180 G4double* const lqu = simulator->lqu;
181 G4double* const time = simulator->nextStateTime;
182
183 for (G4int var = 0, icf0 = 0; var < 6; var++, icf0 += 4)
184 {
185 const G4int icf1 = icf0 + 1;
186
187 if (x[icf1] == 0)
188 {
189 time[var] = Qss_misc::INF;
190 }
191 else
192 {
193 mpr = lqu[var] / x[icf1];
194 if (mpr < 0) { mpr *= -1; }
195 time[var] = t + mpr;
196 }
197 }
198 }
199
200 inline void next_time(G4int i, G4double t)
201 {
202 const G4int cf3 = 4 * i + 3;
203 G4double* const x = simulator->x;
204 G4double* const lqu = simulator->lqu;
205 G4double* const time = simulator->nextStateTime;
206
207 if (likely(x[cf3])) {
208 time[i] = t + std::cbrt(lqu[i] / std::fabs(x[cf3]));
209 } else {
210 time[i] = Qss_misc::INF;
211 }
212 }
213
214 inline void update_quantized_state(G4int i)
215 {
216 const G4int cf0 = i * 4, cf1 = cf0 + 1, cf2 = cf1 + 1;
217 G4double* const q = simulator->q;
218 G4double* const x = simulator->x;
219
220 q[cf0] = x[cf0];
221 q[cf1] = x[cf1];
222 q[cf2] = x[cf2];
223 }
224
225 inline void reset_state(G4int i, G4double value)
226 {
227 G4double* const x = simulator->x;
228 G4double* const q = simulator->q;
229 G4double* const tq = simulator->tq;
230 G4double* const tx = simulator->tx;
231 const G4int idx = 4 * i;
232
233 x[idx] = value;
234
235 simulator->lqu[i] = simulator->dQRel[i] * std::fabs(value);
236 if (simulator->lqu[i] < simulator->dQMin[i])
237 {
238 simulator->lqu[i] = simulator->dQMin[i];
239 }
240 q[idx] = value;
241 q[idx + 1] = q[idx + 2] = tq[i] = tx[i] = 0;
242 }
243
244 inline G4double evaluate_x_poly(G4int i, G4double dt, G4double* p)
245 {
246 return ((p[i + 3] * dt + p[i + 2]) * dt + p[i + 1]) * dt + p[i];
247 }
248
249 inline void advance_time_q(G4int i, G4double dt) // __attribute__((hot))
250 {
251 G4double* const p = simulator->q;
252 const G4int i0 = i, i1 = i0 + 1, i2 = i1 + 1;
253 p[i0] = (p[i2] * dt + p[i1]) * dt + p[i0];
254 p[i1] = p[i1] + 2 * dt * p[i2];
255 }
256
257 inline void advance_time_x(G4int i, G4double dt) // __attribute__((hot))
258 {
259 G4double* const p = simulator->x;
260 const G4int i0 = i, i1 = i0 + 1, i2 = i1 + 1, i3 = i2 + 1;
261 p[i0] = ((p[i3] * dt + p[i2]) * dt + p[i1]) * dt + p[i0];
262 p[i1] = (3 * p[i3] * dt + 2 * p[i2]) * dt + p[i1];
263 p[i2] = p[i2] + 3 * dt * p[i3];
264 }
265
266 G4double min_pos_root(G4double* coeff, G4int order);
267
268 inline G4double min_pos_root_2(G4double* coeff)
269 {
270 G4double mpr = Qss_misc::INF;
271
272 if (coeff[2] == 0 || (1000 * std::fabs(coeff[2])) < std::fabs(coeff[1]))
273 {
274 if (coeff[1] == 0) {
275 mpr = Qss_misc::INF;
276 } else {
277 mpr = -coeff[0] / coeff[1];
278 }
279
280 if (mpr < 0) { mpr = Qss_misc::INF; }
281 }
282 else
283 {
284 G4double disc;
285 disc = coeff[1] * coeff[1] - 4 * coeff[2] * coeff[0];
286 if (disc < 0) // no real roots
287 {
288 mpr = Qss_misc::INF;
289 }
290 else
291 {
292 G4double sd, r1;
293 G4double cf2_d2 = 2 * coeff[2];
294
295 sd = std::sqrt(disc);
296 r1 = (-coeff[1] + sd) / cf2_d2;
297 if (r1 > 0) {
298 mpr = r1;
299 } else {
300 mpr = Qss_misc::INF;
301 }
302 r1 = (-coeff[1] - sd) / cf2_d2;
303 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
304 }
305 }
306
307 return mpr;
308 } // __attribute__((hot))
309
310 inline G4double min_pos_root_3(G4double* coeff)
311 {
312 G4double mpr = Qss_misc::INF;
313 static const G4double sqrt3 = std::sqrt(3);
314
315 if ((coeff[3] == 0) || (1000 * std::fabs(coeff[3]) < std::fabs(coeff[2])))
316 {
317 mpr = min_pos_root_2(coeff);
318 }
319 else if (coeff[0] == 0)
320 {
321 if (coeff[1] == 0)
322 {
323 mpr = -coeff[2] / coeff[3];
324 }
325 else
326 {
327 coeff[0] = coeff[1];
328 coeff[1] = coeff[2];
329 coeff[2] = coeff[3];
330 mpr = min_pos_root_2(coeff);
331 }
332 }
333 else
334 {
335 G4double q, r, disc, q3;
336 G4double val = coeff[2] / 3 / coeff[3];
337 G4double cf32 = coeff[3] * coeff[3];
338 G4double cf22 = coeff[2] * coeff[2];
339 G4double denq = 9 * cf32;
340 G4double denr = 6 * coeff[3] * denq;
341 G4double rcomm = 9 * coeff[3] * coeff[2] * coeff[1] - 2 * cf22 * coeff[2];
342
343 q = (3 * coeff[3] * coeff[1] - cf22) / denq;
344 q3 = q * q * q;
345
346 r = (rcomm - 27 * cf32 * coeff[0]) / denr;
347 disc = q3 + r * r;
348 mpr = Qss_misc::INF;
349
350 if (disc >= 0)
351 {
352 G4double sd, sx, t, r1, rsd;
353 sd = std::sqrt(disc);
354 rsd = r + sd;
355 if (rsd > 0) {
356 sx = std::cbrt(rsd);
357 } else {
358 sx = -std::cbrt(std::fabs(rsd));
359 }
360
361 rsd = r - sd;
362 if (rsd > 0) {
363 t = std::cbrt(rsd);
364 } else {
365 t = -std::cbrt(std::fabs(rsd));
366 }
367
368 r1 = sx + t - val;
369
370 if (r1 > 0) { mpr = r1; }
371 }
372 else
373 {
374 // three real roots
375 G4double rho, th, rho13, costh3, sinth3, spt, smti32, r1;
376 rho = std::sqrt(-q3);
377 th = std::acos(r / rho);
378 rho13 = std::cbrt(rho);
379 costh3 = std::cos(th / 3);
380 sinth3 = std::sin(th / 3);
381 spt = rho13 * 2 * costh3;
382 smti32 = -rho13 * sinth3 * sqrt3;
383 r1 = spt - val;
384 if (r1 > 0) { mpr = r1; }
385 r1 = -spt / 2 - val + smti32;
386 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
387 r1 = r1 - 2 * smti32;
388 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
389 }
390 }
391
392 return mpr;
393 } // __attribute__((hot))
394
395 inline G4double min_pos_root_2_alt(G4double* coeff, G4double cf0Alt)
396 {
397 G4double mpr = Qss_misc::INF;
398 G4double mpr2;
399
400 if (coeff[2] == 0 || (1000 * std::fabs(coeff[2])) < std::fabs(coeff[1]))
401 {
402 if (coeff[1] == 0)
403 {
404 mpr = Qss_misc::INF;
405 }
406 else
407 {
408 mpr = -coeff[0] / coeff[1];
409 mpr2 = -cf0Alt / coeff[1];
410 if (mpr < 0 || (mpr2 > 0 && mpr2 < mpr)) { mpr = mpr2; }
411 }
412
413 if (mpr < 0) { mpr = Qss_misc::INF; }
414 }
415 else
416 {
417 G4double cf1_2 = coeff[1] * coeff[1];
418 G4double cf2_4 = 4 * coeff[2];
419 G4double disc1 = cf1_2 - cf2_4 * coeff[0];
420 G4double disc2 = cf1_2 - cf2_4 * cf0Alt;
421 G4double cf2_d2 = 2 * coeff[2];
422
423 if (unlikely(disc1 < 0 && disc2 < 0))
424 {
425 mpr = Qss_misc::INF;
426 }
427 else if (disc2 < 0)
428 {
429 G4double sd, r1;
430 sd = std::sqrt(disc1);
431 r1 = (-coeff[1] + sd) / cf2_d2;
432 if (r1 > 0) {
433 mpr = r1;
434 } else {
435 mpr = Qss_misc::INF;
436 }
437 r1 = (-coeff[1] - sd) / cf2_d2;
438 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
439 }
440 else if (disc1 < 0)
441 {
442 G4double sd, r1;
443 sd = std::sqrt(disc2);
444 r1 = (-coeff[1] + sd) / cf2_d2;
445 if (r1 > 0) {
446 mpr = r1;
447 } else {
448 mpr = Qss_misc::INF;
449 }
450 r1 = (-coeff[1] - sd) / cf2_d2;
451 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
452 }
453 else
454 {
455 G4double sd1, r1, sd2, r2;
456 sd1 = std::sqrt(disc1);
457 sd2 = std::sqrt(disc2);
458 r1 = (-coeff[1] + sd1) / cf2_d2;
459 r2 = (-coeff[1] + sd2) / cf2_d2;
460
461 if (r1 > 0) {
462 mpr = r1;
463 } else {
464 mpr = Qss_misc::INF;
465 }
466 r1 = (-coeff[1] - sd1) / cf2_d2;
467 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
468
469 if (r2 > 0 && r2 < mpr) { mpr = r2; }
470 r2 = (-coeff[1] - sd2) / cf2_d2;
471 if ((r2 > 0) && (r2 < mpr)) { mpr = r2; }
472 }
473 }
474
475 return mpr;
476 } // __attribute__((hot))
477
478 inline G4double min_pos_root_3_alt(G4double* coeff, G4double cf0Alt)
479 {
480 G4double mpr = Qss_misc::INF;
481 static const G4double sqrt3 = std::sqrt(3);
482
483 if ((coeff[3] == 0) || (1000 * std::fabs(coeff[3]) < std::fabs(coeff[2])))
484 {
485 mpr = min_pos_root_2_alt(coeff, cf0Alt);
486 }
487 else if (coeff[0] == 0)
488 {
489 G4double mpr2;
490 coeff[0] = cf0Alt;
491 mpr = min_pos_root_3(coeff);
492
493 if (coeff[1] == 0)
494 {
495 mpr2 = -coeff[2] / coeff[3];
496 }
497 else
498 {
499 coeff[0] = coeff[1];
500 coeff[1] = coeff[2];
501 coeff[2] = coeff[3];
502 mpr2 = min_pos_root_2(coeff);
503 }
504
505 if (mpr2 > 0 && mpr2 < mpr) { mpr = mpr2; }
506 }
507 else if (cf0Alt == 0)
508 {
509 G4double mpr2;
510 mpr = min_pos_root_3(coeff);
511
512 if (coeff[1] == 0)
513 {
514 mpr2 = -coeff[2] / coeff[3];
515 }
516 else
517 {
518 coeff[0] = coeff[1];
519 coeff[1] = coeff[2];
520 coeff[2] = coeff[3];
521 mpr2 = min_pos_root_2(coeff);
522 }
523
524 if (mpr2 > 0 && mpr2 < mpr) { mpr = mpr2; }
525 }
526 else
527 {
528 G4double q, r, rAlt, disc, discAlt, q3;
529 G4double val = coeff[2] / 3 / coeff[3];
530 G4double cf32 = coeff[3] * coeff[3];
531 G4double cf22 = coeff[2] * coeff[2];
532 G4double denq = 9 * cf32;
533 G4double denr = 6 * coeff[3] * denq;
534 G4double rcomm = 9 * coeff[3] * coeff[2] * coeff[1] - 2 * cf22 * coeff[2];
535
536 q = (3 * coeff[3] * coeff[1] - cf22) / denq;
537 q3 = q * q * q;
538
539 r = (rcomm - 27 * cf32 * coeff[0]) / denr;
540 rAlt = (rcomm - 27 * cf32 * cf0Alt) / denr;
541
542 disc = q3 + r * r;
543 discAlt = q3 + rAlt * rAlt;
544 mpr = Qss_misc::INF;
545
546 if (disc >= 0)
547 {
548 G4double sd, sx, t, r1, rsd;
549 sd = std::sqrt(disc);
550 rsd = r + sd;
551 if (rsd > 0) {
552 sx = std::cbrt(rsd);
553 } else {
554 sx = -std::cbrt(std::fabs(rsd));
555 }
556
557 rsd = r - sd;
558 if (rsd > 0) {
559 t = std::cbrt(rsd);
560 } else {
561 t = -std::cbrt(std::fabs(rsd));
562 }
563
564 r1 = sx + t - val;
565
566 if (r1 > 0) { mpr = r1; }
567
568 if (discAlt >= 0)
569 {
570 G4double sdAlt, sAlt, tAlt, r1Alt, rsdAlt;
571 sdAlt = std::sqrt(discAlt);
572 rsdAlt = rAlt + sdAlt;
573 if (rsdAlt > 0) {
574 sAlt = std::cbrt(rsdAlt);
575 } else {
576 sAlt = -std::cbrt(std::fabs(rsdAlt));
577 }
578
579 rsdAlt = rAlt - sdAlt;
580 if (rsdAlt > 0) {
581 tAlt = std::cbrt(rsdAlt);
582 } else {
583 tAlt = -std::cbrt(std::fabs(rsdAlt));
584 }
585
586 r1Alt = sAlt + tAlt - val;
587
588 if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
589 }
590 else
591 {
592 G4double rho, th, rho13, costh3, sinth3, spt, smti32, r1Alt;
593
594 rho = std::sqrt(-q3);
595 th = std::acos(rAlt / rho);
596 rho13 = std::cbrt(rho);
597 costh3 = std::cos(th / 3);
598 sinth3 = std::sin(th / 3);
599 spt = rho13 * 2 * costh3;
600 smti32 = -rho13 * sinth3 * sqrt3;
601 r1Alt = spt - val;
602 if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
603 r1Alt = -spt / 2 - val + smti32;
604 if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
605 r1Alt = r1Alt - 2 * smti32;
606 if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
607 }
608 }
609 else
610 {
611 G4double rho, th, rho13, costh3, sinth3, spt, smti32, r1;
612
613 rho = std::sqrt(-q3);
614 th = std::acos(r / rho);
615 rho13 = std::cbrt(rho);
616 costh3 = std::cos(th / 3);
617 sinth3 = std::sin(th / 3);
618 spt = rho13 * 2 * costh3;
619 smti32 = -rho13 * sinth3 * sqrt3;
620 r1 = spt - val;
621 if (r1 > 0) { mpr = r1; }
622 r1 = -spt / 2 - val + smti32;
623 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
624 r1 = r1 - 2 * smti32;
625 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
626
627 if (discAlt >= 0)
628 {
629 G4double sdAlt, sAlt, tAlt, r1Alt, rsdAlt;
630 sdAlt = std::sqrt(discAlt);
631 rsdAlt = rAlt + sdAlt;
632 if (rsdAlt > 0) {
633 sAlt = std::cbrt(rsdAlt);
634 } else {
635 sAlt = -std::cbrt(std::fabs(rsdAlt));
636 }
637
638 rsdAlt = rAlt - sdAlt;
639 if (rsdAlt > 0) {
640 tAlt = std::cbrt(rsdAlt);
641 } else {
642 tAlt = -std::cbrt(std::fabs(rsdAlt));
643 }
644
645 r1Alt = sAlt + tAlt - val;
646
647 if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
648 }
649 else
650 {
651 G4double thAlt, costh3Alt, sinth3Alt, sptAlt, smti32Alt, r1Alt;
652 thAlt = std::acos(rAlt / rho);
653 costh3Alt = std::cos(thAlt / 3);
654 sinth3Alt = std::sin(thAlt / 3);
655 sptAlt = rho13 * 2 * costh3Alt;
656 smti32Alt = -rho13 * sinth3Alt * sqrt3;
657 r1Alt = sptAlt - val;
658 if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
659 r1Alt = -sptAlt / 2 - val + smti32Alt;
660 if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
661 r1Alt = r1Alt - 2 * smti32Alt;
662 if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
663 }
664 }
665 }
666
667 return mpr;
668 }
669
670 private:
671
672 QSS_simulator simulator;
673};
674
675#endif
G4double
double G4double
Definition G4Types.hh:83
G4int
int G4int
Definition G4Types.hh:85
G4qss_misc.hh
QSS_simulator
struct QSS_simulator_ * QSS_simulator
Definition G4qss_misc.hh:37
likely
#define likely(x)
Definition G4qss_misc.hh:68
unlikely
#define unlikely(x)
Definition G4qss_misc.hh:67
G4QSS3::recompute_all_state_times
void recompute_all_state_times(G4double t)
Definition G4QSS3.hh:176
G4QSS3::min_pos_root_3
G4double min_pos_root_3(G4double *coeff)
Definition G4QSS3.hh:310
G4QSS3::full_definition
void full_definition(G4double coeff)
Definition G4QSS3.hh:54
G4QSS3::min_pos_root_3_alt
G4double min_pos_root_3_alt(G4double *coeff, G4double cf0Alt)
Definition G4QSS3.hh:478
G4QSS3::min_pos_root_2
G4double min_pos_root_2(G4double *coeff)
Definition G4QSS3.hh:268
G4QSS3::next_time
void next_time(G4int i, G4double t)
Definition G4QSS3.hh:200
G4QSS3::recompute_next_times
void recompute_next_times(G4int *inf, G4double t)
Definition G4QSS3.cc:38
G4QSS3::evaluate_x_poly
G4double evaluate_x_poly(G4int i, G4double dt, G4double *p)
Definition G4QSS3.hh:244
G4QSS3::getSimulator
QSS_simulator getSimulator() const
Definition G4QSS3.hh:50
G4QSS3::update_quantized_state
void update_quantized_state(G4int i)
Definition G4QSS3.hh:214
G4QSS3::min_pos_root_2_alt
G4double min_pos_root_2_alt(G4double *coeff, G4double cf0Alt)
Definition G4QSS3.hh:395
G4QSS3::G4QSS3
G4QSS3(QSS_simulator)
Definition G4QSS3.cc:33
G4QSS3::advance_time_q
void advance_time_q(G4int i, G4double dt)
Definition G4QSS3.hh:249
G4QSS3::advance_time_x
void advance_time_x(G4int i, G4double dt)
Definition G4QSS3.hh:257
G4QSS3::dependencies
void dependencies(G4int i, G4double coeff)
Definition G4QSS3.hh:85
G4QSS3::reset_state
void reset_state(G4int i, G4double value)
Definition G4QSS3.hh:225
G4QSS3::min_pos_root
G4double min_pos_root(G4double *coeff, G4int order)
G4QSS3::order
G4int order() const
Definition G4QSS3.hh:52
Qss_misc::INF
constexpr G4double INF
Definition G4qss_misc.hh:60