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