| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4QSS2 26 // G4QSS2
27 // 27 //
28 // G4QSS2 simulator 28 // G4QSS2 simulator
29 29
30 // Authors: Lucio Santi, Rodrigo Castro (Univ. << 30 // Authors: Lucio Santi, Rodrigo Castro - 2018-2021
31 // ------------------------------------------- 31 // --------------------------------------------------------------------
32 #ifndef _G4QSS2_H_ 32 #ifndef _G4QSS2_H_
33 #define _G4QSS2_H_ 1 << 33 #define _G4QSS2_H_
34 34
35 #include "G4Types.hh" << 35 #include "G4Types.hh" // For G4int, G4double
36 #include "G4qss_misc.hh" 36 #include "G4qss_misc.hh"
37 37
38 #include <cmath> 38 #include <cmath>
39 #include <cassert> 39 #include <cassert>
40 40
41 #define REPORT_CRITICAL_PROBLEM 1 41 #define REPORT_CRITICAL_PROBLEM 1
42 42
43 #ifdef REPORT_CRITICAL_PROBLEM 43 #ifdef REPORT_CRITICAL_PROBLEM
44 #include <cassert> 44 #include <cassert>
45 #include "G4Log.hh" 45 #include "G4Log.hh"
46 #endif 46 #endif
47 47
48 class G4QSS2 48 class G4QSS2
49 { 49 {
50 public: 50 public:
51 51
52 G4QSS2(QSS_simulator sim) : simulator(sim) 52 G4QSS2(QSS_simulator sim) : simulator(sim) {}
53 53
54 inline QSS_simulator getSimulator() const 54 inline QSS_simulator getSimulator() const { return this->simulator; }
55 55
56 inline G4int order() const { return 2; } 56 inline G4int order() const { return 2; }
57 57
58 inline void full_definition(G4double coeff 58 inline void full_definition(G4double coeff)
59 { 59 {
60 G4double* const x = simulator->q; 60 G4double* const x = simulator->q;
61 G4double* const dx = simulator->x; 61 G4double* const dx = simulator->x;
62 G4double* const alg = simulator->alg; 62 G4double* const alg = simulator->alg;
63 63
64 dx[1] = x[9]; 64 dx[1] = x[9];
65 dx[2] = 0; 65 dx[2] = 0;
66 66
67 dx[4] = x[12]; 67 dx[4] = x[12];
68 dx[5] = 0; 68 dx[5] = 0;
69 69
70 dx[7] = x[15]; 70 dx[7] = x[15];
71 dx[8] = 0; 71 dx[8] = 0;
72 72
73 dx[10] = coeff * (alg[2] * x[12] - alg[1 73 dx[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
74 dx[11] = 0; 74 dx[11] = 0;
75 75
76 dx[13] = coeff * (alg[0] * x[15] - alg[2 76 dx[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
77 dx[14] = 0; 77 dx[14] = 0;
78 78
79 dx[16] = coeff * (alg[1] * x[9] - alg[0] 79 dx[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
80 dx[17] = 0; 80 dx[17] = 0;
81 } 81 }
82 82
83 inline void dependencies(G4int i, G4double 83 inline void dependencies(G4int i, G4double coeff)
84 { 84 {
85 G4double* const x = simulator->q; 85 G4double* const x = simulator->q;
86 G4double* const der = simulator->x; 86 G4double* const der = simulator->x;
87 G4double* const alg = simulator->alg; 87 G4double* const alg = simulator->alg;
88 88
89 switch (i) 89 switch (i)
90 { 90 {
91 case 0: 91 case 0:
92 der[10] = coeff * (alg[2] * x[12] - 92 der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
93 der[11] = ((alg[2] * x[13] - x[16] * 93 der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
94 94
95 der[13] = coeff * (alg[0] * x[15] - 95 der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
96 der[14] = ((alg[0] * x[16] - alg[2] 96 der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
97 97
98 der[16] = coeff * (alg[1] * x[9] - a 98 der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
99 der[17] = (-coeff * (alg[0] * x[13] 99 der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
100 return; 100 return;
101 case 1: 101 case 1:
102 der[10] = coeff * (alg[2] * x[12] - 102 der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
103 der[11] = ((alg[2] * x[13] - x[16] * 103 der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
104 104
105 der[13] = coeff * (alg[0] * x[15] - 105 der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
106 der[14] = ((alg[0] * x[16] - alg[2] 106 der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
107 107
108 der[16] = coeff * (alg[1] * x[9] - a 108 der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
109 der[17] = (-coeff * (alg[0] * x[13] 109 der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
110 return; 110 return;
111 case 2: 111 case 2:
112 der[10] = coeff * (alg[2] * x[12] - 112 der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
113 der[11] = ((alg[2] * x[13] - x[16] * 113 der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
114 114
115 der[13] = coeff * (alg[0] * x[15] - 115 der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
116 der[14] = ((alg[0] * x[16] - alg[2] 116 der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
117 117
118 der[16] = coeff * (alg[1] * x[9] - a 118 der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
119 der[17] = (-coeff * (alg[0] * x[13] 119 der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
120 return; 120 return;
121 case 3: 121 case 3:
122 der[1] = x[9]; 122 der[1] = x[9];
123 der[2] = (x[10]) / 2; 123 der[2] = (x[10]) / 2;
124 124
125 der[13] = coeff * (alg[0] * x[15] - 125 der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
126 der[14] = ((alg[0] * x[16] - alg[2] 126 der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
127 127
128 der[16] = coeff * (alg[1] * x[9] - a 128 der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
129 der[17] = (-coeff * (alg[0] * x[13] 129 der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
130 return; 130 return;
131 case 4: 131 case 4:
132 der[4] = x[12]; 132 der[4] = x[12];
133 der[5] = (x[13]) / 2; 133 der[5] = (x[13]) / 2;
134 134
135 der[10] = coeff * (alg[2] * x[12] - 135 der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
136 der[11] = ((alg[2] * x[13] - x[16] * 136 der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
137 137
138 der[16] = coeff * (alg[1] * x[9] - a 138 der[16] = coeff * (alg[1] * x[9] - alg[0] * x[12]);
139 der[17] = (-coeff * (alg[0] * x[13] 139 der[17] = (-coeff * (alg[0] * x[13] - x[10] * alg[1])) / 2;
140 return; 140 return;
141 case 5: 141 case 5:
142 der[7] = x[15]; 142 der[7] = x[15];
143 der[8] = (x[16]) / 2; 143 der[8] = (x[16]) / 2;
144 144
145 der[10] = coeff * (alg[2] * x[12] - 145 der[10] = coeff * (alg[2] * x[12] - alg[1] * x[15]);
146 der[11] = ((alg[2] * x[13] - x[16] * 146 der[11] = ((alg[2] * x[13] - x[16] * alg[1]) * coeff) / 2;
147 147
148 der[13] = coeff * (alg[0] * x[15] - 148 der[13] = coeff * (alg[0] * x[15] - alg[2] * x[9]);
149 der[14] = ((alg[0] * x[16] - alg[2] 149 der[14] = ((alg[0] * x[16] - alg[2] * x[10]) * coeff) / 2;
150 return; 150 return;
151 } 151 }
152 } 152 }
153 153
154 inline void recompute_next_times(G4int* in 154 inline void recompute_next_times(G4int* inf, G4double t)
155 { 155 {
156 G4int i; 156 G4int i;
157 G4double* x = simulator->x; 157 G4double* x = simulator->x;
158 G4double* q = simulator->q; 158 G4double* q = simulator->q;
159 G4double* lqu = simulator->lqu; 159 G4double* lqu = simulator->lqu;
160 G4double* time = simulator->nextStateTim 160 G4double* time = simulator->nextStateTime;
161 161
162 for (i = 0; i < 3; i++) 162 for (i = 0; i < 3; i++)
163 { 163 {
164 const G4int var = inf[i]; 164 const G4int var = inf[i];
165 const G4int icf0 = 3 * var; 165 const G4int icf0 = 3 * var;
166 const G4int icf1 = icf0 + 1; 166 const G4int icf1 = icf0 + 1;
167 const G4int icf2 = icf1 + 1; 167 const G4int icf2 = icf1 + 1;
168 168
169 time[var] = t; 169 time[var] = t;
170 170
171 if (std::fabs(q[icf0] - x[icf0]) < lqu 171 if (std::fabs(q[icf0] - x[icf0]) < lqu[var])
172 { 172 {
173 G4double mpr = -1, mpr2; 173 G4double mpr = -1, mpr2;
174 G4double cf0 = q[icf0] + lqu[var] - 174 G4double cf0 = q[icf0] + lqu[var] - x[icf0];
175 G4double cf1 = q[icf1] - x[icf1]; 175 G4double cf1 = q[icf1] - x[icf1];
176 G4double cf2 = -x[icf2]; 176 G4double cf2 = -x[icf2];
177 G4double cf0Alt = q[icf0] - lqu[var] 177 G4double cf0Alt = q[icf0] - lqu[var] - x[icf0];
178 178
179 if (unlikely(cf2 == 0 || (1000 * std 179 if (unlikely(cf2 == 0 || (1000 * std::fabs(cf2)) < std::fabs(cf1)))
180 { 180 {
181 if (cf1 == 0) { 181 if (cf1 == 0) {
182 mpr = Qss_misc::INF; 182 mpr = Qss_misc::INF;
183 } else 183 } else
184 { 184 {
185 mpr = -cf0 / cf1; 185 mpr = -cf0 / cf1;
186 mpr2 = -cf0Alt / cf1; 186 mpr2 = -cf0Alt / cf1;
187 if (mpr < 0 || (mpr2 > 0 && mpr2 187 if (mpr < 0 || (mpr2 > 0 && mpr2 < mpr)) { mpr = mpr2; }
188 } 188 }
189 189
190 if (mpr < 0) { mpr = Qss_misc::INF 190 if (mpr < 0) { mpr = Qss_misc::INF; }
191 } 191 }
192 else 192 else
193 { 193 {
194 static G4ThreadLocal unsigned long 194 static G4ThreadLocal unsigned long long okCalls=0LL, badCalls= 0LL;
195 constexpr G4double dangerZone = 1. 195 constexpr G4double dangerZone = 1.0e+30;
196 static G4ThreadLocal G4double bigC 196 static G4ThreadLocal G4double bigCf1_pr = dangerZone,
197 bigC 197 bigCf2_pr = dangerZone;
198 static G4ThreadLocal G4double bigC 198 static G4ThreadLocal G4double bigCf1 = 0.0, bigCf2 = 0.0;
199 if( std::abs(cf1) > dangerZone || 199 if( std::abs(cf1) > dangerZone || std::fabs(cf2) > dangerZone )
200 { 200 {
201 badCalls++; 201 badCalls++;
202 if( badCalls == 1 202 if( badCalls == 1
203 || ( badCalls < 1000 && badCa 203 || ( badCalls < 1000 && badCalls % 20 == 0 )
204 || ( 1000 < badCalls && bad 204 || ( 1000 < badCalls && badCalls < 10000 && badCalls % 100 == 0 )
205 || ( 10000 < badCalls && bad 205 || ( 10000 < badCalls && badCalls < 100000 && badCalls % 1000 == 0 )
206 || ( 100000 < badCalls && 206 || ( 100000 < badCalls && badCalls % 10000 == 0 )
207 || ( std::fabs(cf1) > 1.5 * b 207 || ( std::fabs(cf1) > 1.5 * bigCf1_pr || std::fabs(cf2) > 1.5 * bigCf2_pr )
208 ) 208 )
209 { 209 {
210 std::cout << " cf1 = " << std: 210 std::cout << " cf1 = " << std::setw(15) << cf1 << " cf2= " << std::setw(15) << cf2
211 << " badCall # " << 211 << " badCall # " << badCalls << " of " << badCalls + okCalls
212 << " fraction = " < 212 << " fraction = " << double(badCalls) / double(badCalls+okCalls);
213 213
214 if( std::fabs(cf1) > 1.5 * big 214 if( std::fabs(cf1) > 1.5 * bigCf1_pr ) { bigCf1_pr = std::fabs(cf1); std::cout << " Bigger cf1 "; }
215 if( std::fabs(cf2) > 1.5 * big 215 if( std::fabs(cf2) > 1.5 * bigCf2_pr ) { bigCf2_pr = std::fabs(cf2); std::cout << " Bigger cf2 "; }
216 std::cout << std::endl; 216 std::cout << std::endl;
217 } 217 }
218 if( std::fabs(cf1) > 1.5 * bigCf 218 if( std::fabs(cf1) > 1.5 * bigCf1 ) { bigCf1 = std::fabs(cf1); }
219 if( std::fabs(cf2) > 1.5 * bigCf 219 if( std::fabs(cf2) > 1.5 * bigCf2 ) { bigCf2 = std::fabs(cf2); }
220 } 220 }
221 else 221 else
222 { 222 {
223 okCalls++; 223 okCalls++;
224 } 224 }
225 225
226 #ifdef REPORT_CRITICAL_PROBLEM 226 #ifdef REPORT_CRITICAL_PROBLEM
227 constexpr unsigned int exp_limit= 227 constexpr unsigned int exp_limit= 140;
228 constexpr G4double limit= 1.0e+140 228 constexpr G4double limit= 1.0e+140; // std::pow(10,exp_limit));
229 assert( std::fabs( std::pow(10, ex 229 assert( std::fabs( std::pow(10, exp_limit) - limit ) < 1.0e-14*limit );
230 G4bool bad_cf2fac= G4Log(std::fabs 230 G4bool bad_cf2fac= G4Log(std::fabs(cf2))
231 + G4Log(std::max( 231 + G4Log(std::max( std::fabs(cf0), std::fabs(cf0Alt))) > 2*limit;
232 if( std::fabs(cf1) > limit 232 if( std::fabs(cf1) > limit
233 || G4Log(std::fabs(cf2)) 233 || G4Log(std::fabs(cf2))
234 + G4Log(std::max( std::fabs(cf 234 + G4Log(std::max( std::fabs(cf0), std::fabs(cf0Alt))) > 2*exp_limit )
235 { 235 {
236 G4ExceptionDescription ermsg; 236 G4ExceptionDescription ermsg;
237 ermsg << "QSS2: Coefficients exc 237 ermsg << "QSS2: Coefficients exceed tolerable values -- beyond " << limit << G4endl;
238 if( std::fabs(cf1) > limit ) 238 if( std::fabs(cf1) > limit )
239 { 239 {
240 ermsg << " |cf1| = " << cf1 << 240 ermsg << " |cf1| = " << cf1 << " is > " << limit << " (limit)";
241 } 241 }
242 if( bad_cf2fac) 242 if( bad_cf2fac)
243 { 243 {
244 ermsg << " bad cf2-factor: cf 244 ermsg << " bad cf2-factor: cf2 = " << cf2
245 << " product is > " << 2 245 << " product is > " << 2*limit << " (limit)";
246 } 246 }
247 G4Exception("QSS2::recompute_nex 247 G4Exception("QSS2::recompute_next_times",
248 "Field/Qss2-", Fatal 248 "Field/Qss2-", FatalException, ermsg );
249 } 249 }
250 #endif 250 #endif
251 G4double cf1_2 = cf1 * cf1; 251 G4double cf1_2 = cf1 * cf1;
252 G4double cf2_4 = 4 * cf2; 252 G4double cf2_4 = 4 * cf2;
253 G4double disc1 = cf1_2 - cf2_4 * c 253 G4double disc1 = cf1_2 - cf2_4 * cf0;
254 G4double disc2 = cf1_2 - cf2_4 * c 254 G4double disc2 = cf1_2 - cf2_4 * cf0Alt;
255 G4double cf2_d2 = 2 * cf2; 255 G4double cf2_d2 = 2 * cf2;
256 256
257 if (unlikely(disc1 < 0 && disc2 < 257 if (unlikely(disc1 < 0 && disc2 < 0)) // no real roots
258 { 258 {
259 mpr = Qss_misc::INF; 259 mpr = Qss_misc::INF;
260 } 260 }
261 else if (disc2 < 0) 261 else if (disc2 < 0)
262 { 262 {
263 G4double sd, r1; 263 G4double sd, r1;
264 sd = std::sqrt(disc1); 264 sd = std::sqrt(disc1);
265 r1 = (-cf1 + sd) / cf2_d2; 265 r1 = (-cf1 + sd) / cf2_d2;
266 if (r1 > 0) { 266 if (r1 > 0) {
267 mpr = r1; 267 mpr = r1;
268 } else { 268 } else {
269 mpr = Qss_misc::INF; 269 mpr = Qss_misc::INF;
270 } 270 }
271 r1 = (-cf1 - sd) / cf2_d2; 271 r1 = (-cf1 - sd) / cf2_d2;
272 if ((r1 > 0) && (r1 < mpr)) { mp 272 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
273 } 273 }
274 else if (disc1 < 0) 274 else if (disc1 < 0)
275 { 275 {
276 G4double sd, r1; 276 G4double sd, r1;
277 sd = std::sqrt(disc2); 277 sd = std::sqrt(disc2);
278 r1 = (-cf1 + sd) / cf2_d2; 278 r1 = (-cf1 + sd) / cf2_d2;
279 if (r1 > 0) { 279 if (r1 > 0) {
280 mpr = r1; 280 mpr = r1;
281 } else { 281 } else {
282 mpr = Qss_misc::INF; 282 mpr = Qss_misc::INF;
283 } 283 }
284 r1 = (-cf1 - sd) / cf2_d2; 284 r1 = (-cf1 - sd) / cf2_d2;
285 if ((r1 > 0) && (r1 < mpr)) { mp 285 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
286 } 286 }
287 else 287 else
288 { 288 {
289 G4double sd1, r1, sd2, r2; 289 G4double sd1, r1, sd2, r2;
290 sd1 = std::sqrt(disc1); 290 sd1 = std::sqrt(disc1);
291 sd2 = std::sqrt(disc2); 291 sd2 = std::sqrt(disc2);
292 r1 = (-cf1 + sd1) / cf2_d2; 292 r1 = (-cf1 + sd1) / cf2_d2;
293 r2 = (-cf1 + sd2) / cf2_d2; 293 r2 = (-cf1 + sd2) / cf2_d2;
294 if (r1 > 0) { mpr = r1; } 294 if (r1 > 0) { mpr = r1; }
295 else { mpr = Qss_misc::INF; } 295 else { mpr = Qss_misc::INF; }
296 r1 = (-cf1 - sd1) / cf2_d2; 296 r1 = (-cf1 - sd1) / cf2_d2;
297 if ((r1 > 0) && (r1 < mpr)) { mp 297 if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
298 if (r2 > 0 && r2 < mpr) { mpr = 298 if (r2 > 0 && r2 < mpr) { mpr = r2; }
299 r2 = (-cf1 - sd2) / cf2_d2; 299 r2 = (-cf1 - sd2) / cf2_d2;
300 if ((r2 > 0) && (r2 < mpr)) { mp 300 if ((r2 > 0) && (r2 < mpr)) { mpr = r2; }
301 } 301 }
302 } 302 }
303 time[var] += mpr; 303 time[var] += mpr;
304 } 304 }
305 } 305 }
306 } 306 }
307 307
308 inline void recompute_all_state_times(G4do 308 inline void recompute_all_state_times(G4double t)
309 { 309 {
310 G4double mpr; 310 G4double mpr;
311 G4double* const x = simulator->x; 311 G4double* const x = simulator->x;
312 G4double* const lqu = simulator->lqu; 312 G4double* const lqu = simulator->lqu;
313 G4double* const time = simulator->nextSt 313 G4double* const time = simulator->nextStateTime;
314 314
315 for (G4int var = 0, icf0 = 0; var < 6; v 315 for (G4int var = 0, icf0 = 0; var < 6; var++, icf0 += 3)
316 { 316 {
317 const G4int icf1 = icf0 + 1; 317 const G4int icf1 = icf0 + 1;
318 318
319 if (x[icf1] == 0) 319 if (x[icf1] == 0)
320 { 320 {
321 time[var] = Qss_misc::INF; 321 time[var] = Qss_misc::INF;
322 } 322 }
323 else 323 else
324 { 324 {
325 mpr = lqu[var] / x[icf1]; 325 mpr = lqu[var] / x[icf1];
326 if (mpr < 0) { mpr *= -1; } 326 if (mpr < 0) { mpr *= -1; }
327 time[var] = t + mpr; 327 time[var] = t + mpr;
328 } 328 }
329 } 329 }
330 } 330 }
331 331
332 inline void next_time(G4int var, G4double 332 inline void next_time(G4int var, G4double t)
333 { 333 {
334 const G4int cf2 = var * 3 + 2; 334 const G4int cf2 = var * 3 + 2;
335 G4double* const x = simulator->x; 335 G4double* const x = simulator->x;
336 G4double* const lqu = simulator->lqu; 336 G4double* const lqu = simulator->lqu;
337 G4double* const time = simulator->nextSt 337 G4double* const time = simulator->nextStateTime;
338 338
339 if (x[cf2] != 0.0) { 339 if (x[cf2] != 0.0) {
340 time[var] = t + std::sqrt(lqu[var] / s 340 time[var] = t + std::sqrt(lqu[var] / std::fabs(x[cf2]));
341 } else { 341 } else {
342 time[var] = Qss_misc::INF; 342 time[var] = Qss_misc::INF;
343 } 343 }
344 } 344 }
345 345
346 inline void update_quantized_state(G4int i 346 inline void update_quantized_state(G4int i)
347 { 347 {
348 const G4int cf0 = i * 3, cf1 = cf0 + 1; 348 const G4int cf0 = i * 3, cf1 = cf0 + 1;
349 G4double* const q = simulator->q; 349 G4double* const q = simulator->q;
350 G4double* const x = simulator->x; 350 G4double* const x = simulator->x;
351 351
352 q[cf0] = x[cf0]; 352 q[cf0] = x[cf0];
353 q[cf1] = x[cf1]; 353 q[cf1] = x[cf1];
354 } 354 }
355 355
356 inline void reset_state(G4int i, G4double 356 inline void reset_state(G4int i, G4double value)
357 { 357 {
358 G4double* const x = simulator->x; 358 G4double* const x = simulator->x;
359 G4double* const q = simulator->q; 359 G4double* const q = simulator->q;
360 G4double* const tq = simulator->tq; 360 G4double* const tq = simulator->tq;
361 G4double* const tx = simulator->tx; 361 G4double* const tx = simulator->tx;
362 const G4int idx = 3 * i; 362 const G4int idx = 3 * i;
363 363
364 x[idx] = value; 364 x[idx] = value;
365 365
366 simulator->lqu[i] = simulator->dQRel[i] 366 simulator->lqu[i] = simulator->dQRel[i] * std::fabs(value);
367 if (simulator->lqu[i] < simulator->dQMin 367 if (simulator->lqu[i] < simulator->dQMin[i])
368 { 368 {
369 simulator->lqu[i] = simulator->dQMin[i 369 simulator->lqu[i] = simulator->dQMin[i];
370 } 370 }
371 371
372 q[idx] = value; 372 q[idx] = value;
373 q[idx + 1] = tq[i] = tx[i] = 0; 373 q[idx + 1] = tq[i] = tx[i] = 0;
374 } 374 }
375 375
376 inline G4double evaluate_x_poly(G4int i, G 376 inline G4double evaluate_x_poly(G4int i, G4double dt, G4double* p)
377 { 377 {
378 return (p[i + 2] * dt + p[i + 1]) * dt + 378 return (p[i + 2] * dt + p[i + 1]) * dt + p[i];
379 } 379 }
380 380
381 inline void advance_time_q(G4int i, G4doub 381 inline void advance_time_q(G4int i, G4double dt) // __attribute__((hot))
382 { 382 {
383 G4double* const p = simulator->q; 383 G4double* const p = simulator->q;
384 p[i] = p[i] + dt * p[i + 1]; 384 p[i] = p[i] + dt * p[i + 1];
385 } 385 }
386 386
387 inline void advance_time_x(G4int i, G4doub 387 inline void advance_time_x(G4int i, G4double dt) // __attribute__((hot))
388 { 388 {
389 G4double* const p = simulator->x; 389 G4double* const p = simulator->x;
390 const G4int i0 = i, i1 = i0 + 1, i2 = i1 390 const G4int i0 = i, i1 = i0 + 1, i2 = i1 + 1;
391 p[i0] = (p[i2] * dt + p[i1]) * dt + p[i0 391 p[i0] = (p[i2] * dt + p[i1]) * dt + p[i0];
392 p[i1] = p[i1] + 2 * dt * p[i2]; 392 p[i1] = p[i1] + 2 * dt * p[i2];
393 } 393 }
394 394
395 private: 395 private:
396 396
397 QSS_simulator simulator; 397 QSS_simulator simulator;
398 }; 398 };
399 399
400 #endif 400 #endif
401 401