| 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 // G4ConvergenceTester class implementation <<
27 // 26 //
28 // Author: Koi, Tatsumi (SLAC/SCCS) << 27 //
29 // ------------------------------------------- << 28 // Convergence Tests for Monte Carlo results.
>> 29 //
>> 30 // Reference
>> 31 // MCNP(TM) -A General Monte Carlo N-Particle Transport Code
>> 32 // Version 4B
>> 33 // Judith F. Briesmeister, Editor
>> 34 // LA-12625-M, Issued: March 1997, UC 705 and UC 700
>> 35 // CHAPTER 2. GEOMETRY, DATA, PHYSICS, AND MATHEMATICS
>> 36 // VI. ESTIMATION OF THE MONTE CARLO PRECISION
>> 37 //
>> 38 // Positives numbers are assumed for inputs
>> 39 //
>> 40 // Koi, Tatsumi (SLAC/SCCS)
>> 41 //
30 42
31 #include "G4ConvergenceTester.hh" 43 #include "G4ConvergenceTester.hh"
32 #include "G4AutoLock.hh" <<
33 #include <iomanip> <<
34 44
35 namespace << 45 G4ConvergenceTester::G4ConvergenceTester()
36 { << 46 : n(0), sum(0.), mean(0.), var(0.), sd(0.), r(0.), efficiency(0.),
37 G4Mutex aMutex = G4MUTEX_INITIALIZER; << 47 r2eff(0.), r2int(0.), shift(0.), vov(0.), fom(0.), largest(0.),
>> 48 largest_score_happened(0), mean_1(0.), var_1(0.), sd_1(0.), r_1(0.),
>> 49 shift_1(0.), vov_1(0.), fom_1(0.), noBinOfHistory(16), slope(0.),
>> 50 noBinOfPDF(10), minimizer(0), noPass(0), noTotal(8)
>> 51 {
>> 52 nonzero_histories.clear();
>> 53 largest_scores.clear();
>> 54 largest_scores.push_back( 0.0 );
>> 55
>> 56 history_grid.resize( noBinOfHistory , 0 );
>> 57 mean_history.resize( noBinOfHistory , 0.0 );
>> 58 var_history.resize( noBinOfHistory , 0.0 );
>> 59 sd_history.resize( noBinOfHistory , 0.0 );
>> 60 r_history.resize( noBinOfHistory , 0.0 );
>> 61 vov_history.resize( noBinOfHistory , 0.0 );
>> 62 fom_history.resize( noBinOfHistory , 0.0 );
>> 63 shift_history.resize( noBinOfHistory , 0.0 );
>> 64 e_history.resize( noBinOfHistory , 0.0 );
>> 65 r2eff_history.resize( noBinOfHistory , 0.0 );
>> 66 r2int_history.resize( noBinOfHistory , 0.0 );
>> 67
>> 68 timer = new G4Timer();
>> 69 timer->Start();
>> 70 cpu_time.clear();
>> 71 cpu_time.push_back( 0.0 );
38 } 72 }
39 73
40 G4ConvergenceTester::G4ConvergenceTester(const << 74
41 : name(theName) <<
42 { <<
43 nonzero_histories.clear(); <<
44 largest_scores.clear(); <<
45 largest_scores.push_back(0.0); <<
46 <<
47 history_grid.resize(noBinOfHistory, 0); <<
48 mean_history.resize(noBinOfHistory, 0.0); <<
49 var_history.resize(noBinOfHistory, 0.0); <<
50 sd_history.resize(noBinOfHistory, 0.0); <<
51 r_history.resize(noBinOfHistory, 0.0); <<
52 vov_history.resize(noBinOfHistory, 0.0); <<
53 fom_history.resize(noBinOfHistory, 0.0); <<
54 shift_history.resize(noBinOfHistory, 0.0); <<
55 e_history.resize(noBinOfHistory, 0.0); <<
56 r2eff_history.resize(noBinOfHistory, 0.0); <<
57 r2int_history.resize(noBinOfHistory, 0.0); <<
58 <<
59 timer = new G4Timer(); <<
60 timer->Start(); <<
61 cpu_time.clear(); <<
62 cpu_time.push_back(0.0); <<
63 } <<
64 75
65 G4ConvergenceTester::~G4ConvergenceTester() 76 G4ConvergenceTester::~G4ConvergenceTester()
66 { 77 {
67 delete timer; << 78 delete timer;
68 } 79 }
69 80
70 void G4ConvergenceTester::AddScore(G4double x) <<
71 { <<
72 G4AutoLock l(&aMutex); <<
73 81
74 timer->Stop(); <<
75 cpu_time.push_back(timer->GetSystemElapsed() <<
76 <<
77 if(x < 0.0) <<
78 { <<
79 std::ostringstream message; <<
80 message << "Expecting zero or positive num <<
81 << "but received a negative number <<
82 G4Exception("G4ConvergenceTester::AddScore <<
83 JustWarning, message); <<
84 } <<
85 <<
86 if(x == 0.0) <<
87 { <<
88 } <<
89 else <<
90 { <<
91 nonzero_histories.insert(std::pair<G4int, <<
92 if(x > largest_scores.back()) <<
93 { <<
94 // Following search should become faster <<
95 for(auto it = largest_scores.begin(); it <<
96 { <<
97 if(x > *it) <<
98 { <<
99 largest_scores.insert(it, x); <<
100 break; <<
101 } <<
102 } <<
103 82
104 if(largest_scores.size() > 201) << 83 void G4ConvergenceTester::AddScore( G4double x )
105 { << 84 {
106 largest_scores.pop_back(); << 85
107 } << 86 //G4cout << x << G4endl;
108 } << 87
109 sum += x; << 88 timer->Stop();
110 } << 89 cpu_time.push_back( timer->GetSystemElapsed() + timer->GetUserElapsed() );
111 << 90
112 // Data has been added so statistics have no << 91 if ( x == 0.0 )
113 statsAreUpdated = false; << 92 {
114 ++n; << 93 }
115 l.unlock(); << 94 else
116 return; << 95 {
>> 96 nonzero_histories.insert( std::pair< G4int , G4double > ( n , x ) );
>> 97 if ( x > largest_scores.back() )
>> 98 {
>> 99 // Following serch should become faster if begin from bottom.
>> 100 std::vector< G4double >::iterator it;
>> 101 for ( it = largest_scores.begin() ; it != largest_scores.end() ; it++ )
>> 102 {
>> 103 if ( x > *it )
>> 104 {
>> 105 largest_scores.insert( it , x );
>> 106 break;
>> 107 }
>> 108 }
>> 109
>> 110 if ( largest_scores.size() > 201 )
>> 111 {
>> 112 largest_scores.pop_back();
>> 113 }
>> 114 //G4cout << largest_scores.size() << " " << largest_scores.front() << " " << largest_scores.back() << G4endl;
>> 115 }
>> 116 sum += x;
>> 117 }
>> 118
>> 119 n++;
>> 120 return;
117 } 121 }
118 122
>> 123
>> 124
119 void G4ConvergenceTester::calStat() 125 void G4ConvergenceTester::calStat()
120 { 126 {
121 if (n == 0) return; << 127
122 128
123 efficiency = G4double(nonzero_histories.size << 129 efficiency = double( nonzero_histories.size() ) / n;
124 130
125 mean = sum / n; << 131 mean = sum / n;
>> 132
>> 133 G4double sum_x2 = 0.0;
>> 134 var = 0.0;
>> 135 shift = 0.0;
>> 136 vov = 0.0;
>> 137
>> 138 G4double xi;
>> 139 std::map< G4int , G4double >::iterator it;
>> 140 for ( it = nonzero_histories.begin() ; it != nonzero_histories.end() ; it++ )
>> 141 {
>> 142 xi = it->second;
>> 143 sum_x2 += xi * xi;
>> 144 var += ( xi - mean ) * ( xi - mean );
>> 145 shift += ( xi - mean ) * ( xi - mean ) * ( xi - mean );
>> 146 vov += ( xi - mean ) * ( xi - mean ) * ( xi - mean ) * ( xi - mean );
>> 147 }
>> 148
>> 149 var += ( n - nonzero_histories.size() ) * mean * mean;
>> 150 shift += ( n - nonzero_histories.size() ) * mean * mean * mean * ( -1 );
>> 151 vov += ( n - nonzero_histories.size() ) * mean * mean * mean * mean;
>> 152
>> 153 vov = vov / ( var * var ) - 1.0 / n;
>> 154
>> 155 var = var/(n-1);
>> 156
>> 157 sd = std::sqrt ( var );
>> 158
>> 159 r = sd / mean / std::sqrt ( G4double(n) );
>> 160
>> 161 r2eff = ( 1 - efficiency ) / ( efficiency * n );
>> 162 r2int = sum_x2 / ( sum * sum ) - 1 / ( efficiency * n );
>> 163
>> 164 shift = shift / ( 2 * var * n );
>> 165
>> 166 fom = 1 / (r*r) / cpu_time.back();
>> 167
>> 168 // Find Largest History
>> 169 //G4double largest = 0.0;
>> 170 largest = 0.0;
>> 171 largest_score_happened = 0;
>> 172 G4double spend_time_of_largest = 0.0;
>> 173 for ( it = nonzero_histories.begin() ; it != nonzero_histories.end() ; it++ )
>> 174 {
>> 175 if ( std::abs ( it->second ) > largest )
>> 176 {
>> 177 largest = it->second;
>> 178 largest_score_happened = it->first;
>> 179 spend_time_of_largest = cpu_time [ it->first+1 ] - cpu_time [ it->first ];
>> 180 }
>> 181 }
>> 182
>> 183 mean_1 = 0.0;
>> 184 var_1 = 0.0;
>> 185 shift_1 = 0.0;
>> 186 vov_1 = 0.0;
>> 187 sd_1 = 0.0;
>> 188 r_1 = 0.0;
>> 189 vov_1 = 0.0;
>> 190
>> 191 // G4cout << "The largest history = " << largest << G4endl;
>> 192
>> 193 mean_1 = ( sum + largest ) / ( n + 1 );
>> 194
>> 195 for ( it = nonzero_histories.begin() ; it != nonzero_histories.end() ; it++ )
>> 196 {
>> 197 xi = it->second;
>> 198 var_1 += ( xi - mean_1 ) * ( xi - mean_1 );
>> 199 shift_1 += ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 );
>> 200 vov_1 += ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 );
>> 201 }
>> 202 xi = largest;
>> 203 var_1 += ( xi - mean_1 ) * ( xi - mean_1 );
>> 204 shift_1 += ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 );
>> 205 vov_1 += ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 );
>> 206
>> 207 var_1 += ( n - nonzero_histories.size() ) * mean_1 * mean_1;
>> 208 shift_1 += ( n - nonzero_histories.size() ) * mean_1 * mean_1 * mean_1 * ( -1 );
>> 209 vov_1 += ( n - nonzero_histories.size() ) * mean_1 * mean_1 * mean_1 * mean_1;
>> 210
>> 211 vov_1 = vov_1 / ( var_1 * var_1 ) - 1.0 / ( n + 1 );
>> 212
>> 213 var_1 = var_1 / n ;
>> 214
>> 215 sd_1 = std::sqrt ( var_1 );
>> 216
>> 217 r_1 = sd_1 / mean_1 / std::sqrt ( G4double(n + 1) );
>> 218
>> 219 shift_1 = shift_1 / ( 2 * var_1 * ( n + 1 ) );
>> 220
>> 221 fom_1 = 1 / ( r * r ) / ( cpu_time.back() + spend_time_of_largest );
>> 222
>> 223 if ( nonzero_histories.size() < 500 )
>> 224 {
>> 225 G4cout << "Number of non zero history too small to calcuarte SLOPE" << G4endl;
>> 226 }
>> 227 else
>> 228 {
>> 229 G4int i = int ( nonzero_histories.size() );
>> 230
>> 231 // 5% criterion
>> 232 G4int j = int ( i * 0.05 );
>> 233 while ( int( largest_scores.size() ) > j )
>> 234 {
>> 235 largest_scores.pop_back();
>> 236 }
>> 237 calc_slope_fit( largest_scores );
>> 238 }
126 239
127 G4double sum_x2 = 0.0; << 240 calc_grid_point_of_history();
128 var = 0.0; << 241 calc_stat_history();
129 shift = 0.0; <<
130 vov = 0.0; <<
131 <<
132 G4double xi; <<
133 for(const auto& nonzero_historie : nonzero_h <<
134 { <<
135 xi = nonzero_historie.second; <<
136 sum_x2 += xi * xi; <<
137 var += (xi - mean) * (xi - mean); <<
138 shift += (xi - mean) * (xi - mean) * (xi - <<
139 vov += (xi - mean) * (xi - mean) * (xi - m <<
140 } <<
141 <<
142 var += (n - nonzero_histories.size()) * mean <<
143 shift += (n - nonzero_histories.size()) * me <<
144 vov += (n - nonzero_histories.size()) * mean <<
145 <<
146 if(var != 0.0) <<
147 { <<
148 vov = vov / (var * var) - 1.0 / n; <<
149 <<
150 var = var / (n - 1); <<
151 <<
152 sd = std::sqrt(var); <<
153 <<
154 r = sd / mean / std::sqrt(G4double(n)); <<
155 <<
156 r2eff = (1 - efficiency) / (efficiency * n <<
157 r2int = sum_x2 / (sum * sum) - 1 / (effici <<
158 <<
159 shift = shift / (2 * var * n); <<
160 <<
161 fom = 1 / (r * r) / cpu_time.back(); <<
162 } <<
163 <<
164 // Find Largest History <<
165 // G4double largest = 0.0; <<
166 largest = 0.0; <<
167 largest_score_happened = 0; <<
168 G4double spend_time_of_largest = 0.0; <<
169 for(const auto& nonzero_historie : nonzero_h <<
170 { <<
171 if(std::abs(nonzero_historie.second) > lar <<
172 { <<
173 largest = nonzero_histori <<
174 largest_score_happened = nonzero_histori <<
175 spend_time_of_largest = <<
176 cpu_time[nonzero_historie.first + 1] - <<
177 } <<
178 } <<
179 <<
180 mean_1 = 0.0; <<
181 var_1 = 0.0; <<
182 shift_1 = 0.0; <<
183 vov_1 = 0.0; <<
184 sd_1 = 0.0; <<
185 r_1 = 0.0; <<
186 vov_1 = 0.0; <<
187 <<
188 mean_1 = (sum + largest) / (n + 1); <<
189 <<
190 for(const auto& nonzero_historie : nonzero_h <<
191 { <<
192 xi = nonzero_historie.second; <<
193 var_1 += (xi - mean_1) * (xi - mean_1); <<
194 shift_1 += (xi - mean_1) * (xi - mean_1) * <<
195 vov_1 += (xi - mean_1) * (xi - mean_1) * ( <<
196 } <<
197 xi = largest; <<
198 var_1 += (xi - mean_1) * (xi - mean_1); <<
199 shift_1 += (xi - mean_1) * (xi - mean_1) * ( <<
200 vov_1 += (xi - mean_1) * (xi - mean_1) * (xi <<
201 <<
202 var_1 += (n - nonzero_histories.size()) * me <<
203 <<
204 if(var_1 != 0.0) <<
205 { <<
206 shift_1 += (n - nonzero_histories.size()) <<
207 vov_1 += (n - nonzero_histories.size()) * <<
208 <<
209 vov_1 = vov_1 / (var_1 * var_1) - 1.0 / (n <<
210 <<
211 var_1 = var_1 / n; <<
212 <<
213 sd_1 = std::sqrt(var_1); <<
214 <<
215 r_1 = sd_1 / mean_1 / std::sqrt(G4double(n <<
216 <<
217 shift_1 = shift_1 / (2 * var_1 * (n + 1)); <<
218 <<
219 fom_1 = 1 / (r * r) / (cpu_time.back() + s <<
220 } <<
221 <<
222 if(nonzero_histories.size() < 500) <<
223 { <<
224 calcSLOPE = false; <<
225 } <<
226 else <<
227 { <<
228 auto i = G4int(nonzero_histories.size()); <<
229 <<
230 // 5% criterion <<
231 auto j = G4int(i * 0.05); <<
232 while(G4int(largest_scores.size()) > j) <<
233 { <<
234 largest_scores.pop_back(); <<
235 } <<
236 calc_slope_fit(largest_scores); <<
237 } <<
238 <<
239 calc_grid_point_of_history(); <<
240 calc_stat_history(); <<
241 <<
242 // statistics have been calculated and this <<
243 // to be called again until data has been ad <<
244 statsAreUpdated = true; <<
245 } 242 }
246 243
>> 244
>> 245
247 void G4ConvergenceTester::calc_grid_point_of_h 246 void G4ConvergenceTester::calc_grid_point_of_history()
248 { 247 {
249 // histroy_grid [ 0,,,15 ] << 248
250 // history_grid [0] 1/16 ,,, history_grid [1 << 249 // histroy_grid [ 0,,,15 ]
251 // if number of event is x then history_grid << 250 // history_grid [0] 1/16 ,,, history_grid [15] 16/16
252 // 16 -> noBinOfHisotry << 251 // if number of event is x then history_grid [15] become x-1.
253 << 252 // 16 -> noBinOfHisotry
254 for(G4int i = 1; i <= noBinOfHistory; ++i) << 253
255 { << 254 G4int i;
256 history_grid[i - 1] = G4int(n / (G4double( << 255 for ( i = 1 ; i <= noBinOfHistory ; i++ )
257 } << 256 {
>> 257 history_grid [ i-1 ] = int ( n / ( double( noBinOfHistory ) ) * i - 0.1 );
>> 258 //G4cout << "history_grid " << i-1 << " " << history_grid [ i-1 ] << G4endl;
>> 259 }
>> 260
258 } 261 }
259 262
>> 263
>> 264
260 void G4ConvergenceTester::calc_stat_history() 265 void G4ConvergenceTester::calc_stat_history()
261 { 266 {
262 if(history_grid[0] == 0) << 267 // G4cout << "i/16 till_ith mean var sd r vov fom shift e r2eff r2int" << G4endl;
263 { << 268
264 showHistory = false; << 269 G4int i;
265 return; << 270 for ( i = 1 ; i <= noBinOfHistory ; i++ )
266 } << 271 {
267 << 272
268 for(G4int i = 0; i < noBinOfHistory; ++i) << 273 G4int ith = history_grid [ i-1 ];
269 { << 274
270 G4int ith = history_grid[i]; << 275 G4int nonzero_till_ith = 0;
271 << 276 G4double xi;
272 G4int nonzero_till_ith = 0; << 277 G4double mean_till_ith = 0.0;
273 G4double xi; << 278 std::map< G4int , G4double >::iterator it;
274 G4double mean_till_ith = 0.0; << 279
275 << 280 for ( it = nonzero_histories.begin() ; it !=nonzero_histories.end() ; it++ )
276 for(const auto& itr : nonzero_histories) <<
277 { <<
278 if(itr.first <= ith) <<
279 { 281 {
280 xi = itr.second; << 282 if ( it->first <= ith )
281 mean_till_ith += xi; << 283 {
282 ++nonzero_till_ith; << 284 xi = it->second;
>> 285 mean_till_ith += xi;
>> 286 nonzero_till_ith++;
>> 287 }
283 } 288 }
284 } <<
285 289
286 if(nonzero_till_ith == 0) << 290 mean_till_ith = mean_till_ith / ( ith+1 );
287 { << 291 mean_history [ i-1 ] = mean_till_ith;
288 continue; << 292
289 } << 293 G4double sum_x2_till_ith = 0.0;
290 << 294 G4double var_till_ith = 0.0;
291 mean_till_ith = mean_till_ith / (ith + 1 << 295 G4double vov_till_ith = 0.0;
292 mean_history[i] = mean_till_ith; << 296 G4double shift_till_ith = 0.0;
293 << 297
294 G4double sum_x2_till_ith = 0.0; << 298 for ( it = nonzero_histories.begin() ; it !=nonzero_histories.end() ; it++ )
295 G4double var_till_ith = 0.0; <<
296 G4double vov_till_ith = 0.0; <<
297 G4double shift_till_ith = 0.0; <<
298 <<
299 for(const auto& itr : nonzero_histories) <<
300 { <<
301 if(itr.first <= ith) <<
302 { 299 {
303 xi = itr.second; << 300 if ( it->first <= ith )
304 sum_x2_till_ith += std::pow(xi, 2.0); << 301 {
305 var_till_ith += std::pow(xi - mean_til << 302 xi = it->second;
306 shift_till_ith += std::pow(xi - mean_t << 303 sum_x2_till_ith += xi * xi;
307 vov_till_ith += std::pow(xi - mean_til << 304 var_till_ith += ( xi - mean_till_ith ) * ( xi - mean_till_ith );
>> 305 shift_till_ith += ( xi - mean_till_ith ) * ( xi - mean_till_ith ) * ( xi - mean_till_ith );
>> 306 vov_till_ith += ( xi - mean_till_ith ) * ( xi - mean_till_ith ) * ( xi - mean_till_ith ) * ( xi - mean_till_ith );
>> 307 }
308 } 308 }
309 } <<
310 309
311 var_till_ith += << 310 var_till_ith += ( (ith+1) - nonzero_till_ith ) * mean_till_ith * mean_till_ith;
312 ((ith + 1) - nonzero_till_ith) * std::po << 311
313 vov_till_ith += << 312 vov_till_ith += ( (ith+1) - nonzero_till_ith ) * mean_till_ith * mean_till_ith * mean_till_ith * mean_till_ith ;
314 ((ith + 1) - nonzero_till_ith) * std::po << 313 vov_till_ith = vov_till_ith / ( var_till_ith * var_till_ith ) - 1.0 / (ith+1);
315 << 314 vov_history [ i-1 ] = vov_till_ith;
316 G4double sum_till_ith = mean_till_ith * (i << 315
317 << 316 var_till_ith = var_till_ith / ( ith+1 - 1 );
318 if(!(std::fabs(var_till_ith) > 0.0)) << 317 var_history [ i-1 ] = var_till_ith;
319 { << 318 sd_history [ i-1 ] = std::sqrt( var_till_ith );
320 continue; << 319 r_history [ i-1 ] = std::sqrt( var_till_ith ) / mean_till_ith / std::sqrt ( 1.0*(ith+1) );
321 } << 320
322 if(!(std::fabs(mean_till_ith) > 0.0)) << 321 fom_history [ i-1 ] = 1 / ( r_history [ i-1 ] * r_history [ i-1 ] ) / cpu_time [ ith ];
323 { << 322
324 continue; << 323 shift_till_ith += ( (ith+1) - nonzero_till_ith ) * mean_till_ith * mean_till_ith * mean_till_ith * ( -1 );
325 } << 324 shift_till_ith = shift_till_ith / ( 2 * var_till_ith * (ith+1) );
326 if(!(std::fabs(sum_till_ith) > 0.0)) << 325 shift_history [ i-1 ] = shift_till_ith;
327 { << 326
328 continue; << 327 e_history [ i-1 ] = 1.0*nonzero_till_ith / (ith+1);
329 } << 328 r2eff_history [ i-1 ] = ( 1 - e_history [ i-1 ] ) / ( e_history [ i-1 ] * (ith+1) );
330 << 329
331 vov_till_ith = vov_till_ith / std::pow(var << 330 G4double sum_till_ith = mean_till_ith * (ith+1);
332 vov_history[i] = vov_till_ith; << 331 r2int_history [ i-1 ] = ( sum_x2_till_ith ) / ( sum_till_ith * sum_till_ith ) - 1 / ( e_history [ i-1 ] * (ith+1) );
333 << 332
334 var_till_ith = var_till_ith / (ith + 1 - << 333 }
335 var_history[i] = var_till_ith; << 334
336 sd_history[i] = std::sqrt(var_till_ith); << 335 }
337 r_history[i] = << 336
338 std::sqrt(var_till_ith) / mean_till_ith << 337
339 << 338
340 if(std::fabs(cpu_time[ith]) > 0.0 && std:: << 339 void G4ConvergenceTester::ShowResult()
341 { << 340 {
342 fom_history[i] = 1.0 / std::pow(r_histor << 341 calStat();
343 } << 342
344 else << 343 G4cout << "EFFICIENCY = " << efficiency << G4endl;
345 { << 344 G4cout << "MEAN = " << mean << G4endl;
346 fom_history[i] = 0.0; << 345 G4cout << "VAR = " << var << G4endl;
347 } << 346 G4cout << "SD = " << sd << G4endl;
348 << 347 G4cout << "R = "<< r << G4endl;
349 shift_till_ith += << 348 G4cout << "SHIFT = "<< shift << G4endl;
350 ((ith + 1) - nonzero_till_ith) * std::po << 349 G4cout << "VOV = "<< vov << G4endl;
351 shift_till_ith = shift_till_ith / (2 * v << 350 G4cout << "FOM = "<< fom << G4endl;
352 shift_history[i] = shift_till_ith; << 351
353 << 352 G4cout << "THE LARGEST SCORE = " << largest << " and it happend at " << largest_score_happened << "th event" << G4endl;
354 e_history[i] = 1.0 * nonzero_till_ith / (i << 353 G4cout << "Affected Mean = " << mean_1 << " and its ratio to orignal is " << mean_1/mean << G4endl;
355 if(std::fabs(e_history[i]) > 0.0) << 354 G4cout << "Affected VAR = " << var_1 << " and its ratio to orignal is " << var_1/var << G4endl;
356 { << 355 G4cout << "Affected R = " << r_1 << " and its ratio to orignal is " << r_1/r << G4endl;
357 r2eff_history[i] = (1 - e_history[i]) / << 356 G4cout << "Affected SHIFT = " << shift_1 << " and its ratio to orignal is " << shift_1/shift << G4endl;
358 << 357 G4cout << "Affected FOM = " << fom_1 << " and its ratio to orignal is " << fom_1/fom << G4endl;
359 r2int_history[i] = (sum_x2_till_ith) / s << 358
360 1 / (e_history[i] * ( << 359 check_stat_history();
361 } << 360
362 } << 361 // check SLOPE and output result
363 } << 362 if ( slope >= 3 )
364 << 363 {
365 void G4ConvergenceTester::ShowResult(std::ostr << 364 noPass++;
366 { << 365 G4cout << "SLOPE is large enough" << G4endl;
367 // if data has been added since the last com << 366 }
368 // (not statsAreUpdated) call calStat to rec << 367 else
369 if(!statsAreUpdated) << 368 {
370 { << 369 G4cout << "SLOPE is not large enough" << G4endl;
371 calStat(); << 370 }
372 } << 371
373 << 372 G4cout << "This result passes " << noPass << " / "<< noTotal << " Convergence Test." << G4endl;
374 out << std::setprecision(6); << 373 G4cout << G4endl;
375 << 374
376 out << G4endl; << 375 }
377 out << "G4ConvergenceTester Output Result of << 376
378 out << std::setw(20) << "EFFICIENCY = " << s << 377 void G4ConvergenceTester::ShowHistory()
379 << G4endl; << 378 {
380 out << std::setw(20) << "MEAN = " << std::se << 379 G4cout << "i/" << noBinOfHistory << " till_ith mean var sd r vov fom shift e r2eff r2int" << G4endl;
381 out << std::setw(20) << "VAR = " << std::set << 380 for ( G4int i = 1 ; i <= noBinOfHistory ; i++ )
382 out << std::setw(20) << "SD = " << std::setw << 381 {
383 out << std::setw(20) << "R = " << std::setw( << 382 G4cout << i << " "
384 out << std::setw(20) << "SHIFT = " << std::s << 383 << history_grid [ i-1 ] << " "
385 out << std::setw(20) << "VOV = " << std::set << 384 << mean_history [ i-1 ] << " "
386 out << std::setw(20) << "FOM = " << std::set << 385 << var_history [ i-1 ] << " "
387 << 386 << sd_history [ i-1 ] << " "
388 out << std::setw(20) << "THE LARGEST SCORE = << 387 << r_history [ i-1 ] << " "
389 << " and it happened at " << largest_sco << 388 << vov_history [ i-1 ] << " "
390 << G4endl; << 389 << fom_history [ i-1 ] << " "
391 if(mean != 0) << 390 << shift_history [ i-1 ] << " "
392 { << 391 << e_history [ i-1 ] << " "
393 out << std::setw(20) << "Affected Mean = " << 392 << r2eff_history [ i-1 ] << " "
394 << " and its ratio to original is " << << 393 << r2int_history [ i-1 ] << " "
395 } << 394 << G4endl;
396 else << 395 }
397 { << 396 }
398 out << std::setw(20) << "Affected Mean = " << 397
399 << G4endl; << 398 void G4ConvergenceTester::check_stat_history()
400 } << 399 {
401 if(var != 0) << 400
402 { << 401 // 1 sigma rejection for null hypothesis
403 out << std::setw(20) << "Affected VAR = " << 402
404 << " and its ratio to original is " << << 403 std::vector<G4double> first_ally;
405 } << 404 std::vector<G4double> second_ally;
406 else << 405
407 { << 406 // use 2nd half of hisories
408 out << std::setw(20) << "Affected VAR = " << 407 G4int N = mean_history.size() / 2;
409 << G4endl; << 408 G4int i;
410 } << 409
411 if(r != 0) << 410 G4double pearson_r;
412 { << 411 G4double t;
413 out << std::setw(20) << "Affected R = " << << 412
414 << " and its ratio to original is " << << 413 first_ally.resize( N );
415 } << 414 second_ally.resize( N );
416 else << 415
417 { << 416 // Mean
418 out << std::setw(20) << "Affected R = " << << 417
419 } << 418 for ( i = 0 ; i < N ; i++ )
420 if(shift != 0) << 419 {
421 { << 420 first_ally [ i ] = history_grid [ N + i ];
422 out << std::setw(20) << "Affected SHIFT = << 421 second_ally [ i ] = mean_history [ N + i ];
423 << " and its ratio to original is " << << 422 }
424 } << 423
425 else << 424 pearson_r = calc_Pearson_r ( N , first_ally , second_ally );
426 { << 425 t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) );
427 out << std::setw(20) << "Affected SHIFT = << 426
428 << G4endl; << 427 if ( t < 0.429318 ) // Student t of (Degree of freedom = N-2 )
429 } << 428 {
430 if(fom != 0) << 429 G4cout << "MEAN distribution is RANDOM" << G4endl;
431 { << 430 noPass++;
432 out << std::setw(20) << "Affected FOM = " << 431 }
433 << " and its ratio to original is " << << 432 else
434 } << 433 {
435 else << 434 G4cout << "MEAN distribution is not RANDOM" << G4endl;
436 { << 435 }
437 out << std::setw(20) << "Affected FOM = " << 436
438 << G4endl; << 437
439 } << 438 // R
440 << 439
441 if(!showHistory) << 440 for ( i = 0 ; i < N ; i++ )
442 { << 441 {
443 out << "Number of events of this run is to << 442 first_ally [ i ] = 1.0 / std::sqrt ( G4double(history_grid [ N + i ]) );
444 << G4endl; << 443 second_ally [ i ] = r_history [ N + i ];
445 return; << 444 }
446 } << 445
447 << 446 pearson_r = calc_Pearson_r ( N , first_ally , second_ally );
448 check_stat_history(out); << 447 t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) );
449 << 448
450 // check SLOPE and output result << 449 if ( t > 1.090546 )
451 if(calcSLOPE) << 450 {
452 { << 451 G4cout << "r follows 1/std::sqrt(N)" << G4endl;
453 if(slope >= 3) <<
454 { <<
455 noPass++; 452 noPass++;
456 out << "SLOPE is large enough" << G4endl << 453 }
457 } << 454 else
458 else << 455 {
459 { << 456 G4cout << "r does not follow 1/std::sqrt(N)" << G4endl;
460 out << "SLOPE is not large enough" << G4 << 457 }
461 } << 458
462 } << 459 if ( is_monotonically_decrease( second_ally ) == true )
463 else << 460 {
464 { << 461 G4cout << "r is monotonically decrease " << G4endl;
465 out << "Number of non zero history too sma << 462 }
466 } << 463 else
467 << 464 {
468 out << "This result passes " << noPass << " << 465 G4cout << "r is NOT monotonically decrease " << G4endl;
469 << " Convergence Test." << G4endl; << 466 }
470 out << G4endl; << 467
471 } << 468 if ( r_history.back() < 0.1 )
472 << 469 {
473 void G4ConvergenceTester::ShowHistory(std::ost << 470 G4cout << "r is less than 0.1. r = " << r_history.back() << G4endl;
474 { << 471 noPass++;
475 if(!showHistory) << 472 }
476 { << 473 else
477 out << "Number of events of this run is to << 474 {
478 << G4endl; << 475 G4cout << "r is NOT less than 0.1. r = " << r_history.back() << G4endl;
479 return; << 476 }
480 } << 477
481 << 478
482 out << std::setprecision(6); << 479 // VOV
483 << 480 for ( i = 0 ; i < N ; i++ )
484 out << G4endl; << 481 {
485 out << "G4ConvergenceTester Output History o << 482 first_ally [ i ] = 1.0 / history_grid [ N + i ];
486 out << "i/" << noBinOfHistory << " till_ith << 483 second_ally [ i ] = vov_history [ N + i ];
487 << "var" << std::setw(13) << "sd" << std << 484 }
488 << "vov" << std::setw(13) << "fom" << st << 485
489 << std::setw(13) << "e" << std::setw(13) << 486 pearson_r = calc_Pearson_r ( N , first_ally , second_ally );
490 << "r2int" << G4endl; << 487 t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) );
491 for(G4int i = 1; i <= noBinOfHistory; i++) << 488
492 { << 489 if ( t > 1.090546 )
493 out << std::setw(4) << i << " " << std::se << 490 {
494 << std::setw(13) << mean_history[i - 1 << 491 G4cout << "VOV follows 1/std::sqrt(N)" << G4endl;
495 << var_history[i - 1] << std::setw(13) << 492 noPass++;
496 << std::setw(13) << r_history[i - 1] < << 493 }
497 << vov_history[i - 1] << std::setw(13) << 494 else
498 << std::setw(13) << shift_history[i - << 495 {
499 << e_history[i - 1] << std::setw(13) < << 496 G4cout << "VOV does not follow 1/std::sqrt(N)" << G4endl;
500 << std::setw(13) << r2int_history[i - << 497 }
501 } << 498
502 } << 499 if ( is_monotonically_decrease( second_ally ) == true )
503 << 500 {
504 void G4ConvergenceTester::check_stat_history(s << 501 G4cout << "VOV is monotonically decrease " << G4endl;
505 { << 502 }
506 // 1 sigma rejection for null hypothesis << 503 else
507 << 504 {
508 std::vector<G4double> first_ally; << 505 G4cout << "VOV is NOT monotonically decrease " << G4endl;
509 std::vector<G4double> second_ally; << 506 }
510 << 507
511 // use 2nd half of hisories << 508 // FOM
512 std::size_t N = mean_history.size() / 2; << 509
513 std::size_t i; << 510 for ( i = 0 ; i < N ; i++ )
514 << 511 {
515 G4double pearson_r; << 512 first_ally [ i ] = history_grid [ N + i ];
516 G4double t; << 513 second_ally [ i ] = fom_history [ N + i ];
517 << 514 }
518 first_ally.resize(N); << 515
519 second_ally.resize(N); << 516 pearson_r = calc_Pearson_r ( N , first_ally , second_ally );
520 << 517 t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) );
521 G4double sum_of_var = << 518
522 std::accumulate(var_history.begin(), var_h << 519 if ( t < 0.429318 )
523 if(sum_of_var == 0.0) << 520 {
524 { << 521 G4cout << "FOM distribution is RANDOM" << G4endl;
525 out << "Variances in all historical grids << 522 noPass++;
526 out << "Terminating checking behavior of s << 523 }
527 return; << 524 else
528 } << 525 {
529 << 526 G4cout << "FOM distribution is not RANDOM" << G4endl;
530 // Mean << 527 }
531 << 528
532 for(i = 0; i < N; ++i) << 529 }
533 { << 530
534 first_ally[i] = history_grid[N + i]; << 531
535 second_ally[i] = mean_history[N + i]; << 532
536 } << 533 G4double G4ConvergenceTester::calc_Pearson_r ( G4int N , std::vector<G4double> first_ally , std::vector<G4double> second_ally )
537 << 534 {
538 pearson_r = calc_Pearson_r((G4int)N, first_a << 535 G4double first_mean = 0.0;
539 t = pearson_r * std::sqrt((N - 2) / << 536 G4double second_mean = 0.0;
540 << 537
541 if(t < 0.429318) // Student t of (Degree of << 538 G4int i;
542 { << 539 for ( i = 0 ; i < N ; i++ )
543 out << "MEAN distribution is RANDOM" << G << 540 {
544 noPass++; << 541 first_mean += first_ally [ i ];
545 } << 542 second_mean += second_ally [ i ];
546 else << 543 }
547 { << 544 first_mean = first_mean / N;
548 out << "MEAN distribution is not RANDOM" < << 545 second_mean = second_mean / N;
549 } << 546
550 << 547 G4double a = 0.0;
551 // R << 548 for ( i = 0 ; i < N ; i++ )
552 << 549 {
553 for(i = 0; i < N; ++i) << 550 a += ( first_ally [ i ] - first_mean ) * ( second_ally [ i ] - second_mean );
554 { << 551 }
555 first_ally[i] = 1.0 / std::sqrt(G4double( << 552
556 second_ally[i] = r_history[N + i]; << 553 G4double b1 = 0.0;
557 } << 554 G4double b2 = 0.0;
558 << 555 for ( i = 0 ; i < N ; i++ )
559 pearson_r = calc_Pearson_r(G4int(N), first_a << 556 {
560 t = pearson_r * std::sqrt((N - 2) / << 557 b1 += ( first_ally [ i ] - first_mean ) * ( first_ally [ i ] - first_mean );
561 << 558 b2 += ( second_ally [ i ] - second_mean ) * ( second_ally [ i ] - second_mean );
562 if(t > 1.090546) << 559 }
563 { << 560
564 out << "r follows 1/std::sqrt(N)" << G4end << 561 G4double r = a / std::sqrt ( b1 * b2 );
565 noPass++; << 562
566 } << 563 return r;
567 else << 564 }
568 { << 565
569 out << "r does not follow 1/std::sqrt(N)" << 566
570 } << 567
571 << 568 G4bool G4ConvergenceTester::is_monotonically_decrease ( std::vector<G4double> ally )
572 if(is_monotonically_decrease(second_ally)) << 569 {
573 { << 570
574 out << "r is monotonically decrease " << G << 571 std::vector<G4double>::iterator it;
575 } << 572 for ( it = ally.begin() ; it != ally.end() - 1 ; it++ )
576 else << 573 {
577 { << 574 if ( *it < *(it+1) ) return FALSE;
578 out << "r is NOT monotonically decrease " << 575 }
579 } << 576
580 << 577 noPass++;
581 if(r_history.back() < 0.1) << 578 return TRUE;
582 { << 579 }
583 out << "r is less than 0.1. r = " << r_his << 580
584 noPass++; << 581
585 } << 582
586 else << 583 //void G4ConvergenceTester::calc_slope_fit ( std::vector<G4double> largest_socres )
587 { << 584 void G4ConvergenceTester::calc_slope_fit ( std::vector<G4double> )
588 out << "r is NOT less than 0.1. r = " << r << 585 {
589 } << 586
590 << 587 // create PDF bins
591 // VOV << 588 G4double max = largest_scores.front();
592 for(i = 0; i < N; ++i) << 589 G4int last = int ( largest_scores.size() );
593 { << 590 G4double min = 0.0;
594 first_ally[i] = 1.0 / history_grid[N + i] << 591 if ( largest_scores.back() != 0 )
595 second_ally[i] = vov_history[N + i]; << 592 {
596 } << 593 min = largest_scores.back();
597 << 594 }
598 pearson_r = calc_Pearson_r(G4int(N), first_a << 595 else
599 t = pearson_r * std::sqrt((N - 2) / << 596 {
600 << 597 min = largest_scores[ last-1 ];
601 if(t > 1.090546) << 598 last = last - 1;
602 { << 599 }
603 out << "VOV follows 1/std::sqrt(N)" << G4e << 600
604 noPass++; << 601 //G4cout << "largest " << max << G4endl;
605 } << 602 //G4cout << "last " << min << G4endl;
606 else << 603
607 { << 604 if ( max*0.99 < min )
608 out << "VOV does not follow 1/std::sqrt(N) << 605 {
609 } << 606 // upper limit is assumed to have been reached
610 << 607 slope = 10.0;
611 if(is_monotonically_decrease(second_ally)) << 608 return;
612 { << 609 }
613 out << "VOV is monotonically decrease " << << 610
614 } << 611 std::vector < G4double > pdf_grid;
615 else << 612
616 { << 613 pdf_grid.resize( noBinOfPDF+1 ); // no grid = no bins + 1
617 out << "VOV is NOT monotonically decrease << 614 pdf_grid[ 0 ] = max;
618 } << 615 pdf_grid[ noBinOfPDF ] = min;
619 << 616 G4double log10_max = std::log10( max );
620 // FOM << 617 G4double log10_min = std::log10( min );
621 << 618 G4double log10_delta = log10_max - log10_min;
622 for(i = 0; i < N; ++i) << 619 for ( G4int i = 1 ; i < noBinOfPDF ; i++ )
623 { << 620 {
624 first_ally[i] = history_grid[N + i]; << 621 pdf_grid[i] = std::pow ( 10.0 , log10_max - log10_delta/10.0*(i) );
625 second_ally[i] = fom_history[N + i]; << 622 //G4cout << "pdf i " << i << " " << pdf_grid[i] << G4endl;
626 } << 623 }
627 << 624
628 pearson_r = calc_Pearson_r(G4int(N), std::mo << 625 std::vector < G4double > pdf;
629 t = pearson_r * std::sqrt((N - 2) / << 626 pdf.resize( noBinOfPDF );
630 << 627
631 if(t < 0.429318) << 628 for ( G4int j=0 ; j < last ; j ++ )
632 { << 629 {
633 out << "FOM distribution is RANDOM" << G4e << 630 for ( G4int i = 0 ; i < 11 ; i++ )
634 noPass++; <<
635 } <<
636 else <<
637 { <<
638 out << "FOM distribution is not RANDOM" << <<
639 } <<
640 } <<
641 <<
642 G4double G4ConvergenceTester::calc_Pearson_r(G <<
643 s <<
644 s <<
645 { <<
646 G4double first_mean = 0.0; <<
647 G4double second_mean = 0.0; <<
648 <<
649 G4int i; <<
650 for(i = 0; i < N; i++) <<
651 { <<
652 first_mean += first_ally[i]; <<
653 second_mean += second_ally[i]; <<
654 } <<
655 first_mean = first_mean / N; <<
656 second_mean = second_mean / N; <<
657 <<
658 G4double a = 0.0; <<
659 for(i = 0; i < N; ++i) <<
660 { <<
661 a += (first_ally[i] - first_mean) * (secon <<
662 } <<
663 <<
664 G4double b1 = 0.0; <<
665 G4double b2 = 0.0; <<
666 for(i = 0; i < N; ++i) <<
667 { <<
668 b1 += (first_ally[i] - first_mean) * (firs <<
669 b2 += (second_ally[i] - second_mean) * (se <<
670 } <<
671 <<
672 G4double rds = a / std::sqrt(b1 * b2); <<
673 <<
674 return rds; <<
675 } <<
676 <<
677 G4bool G4ConvergenceTester::is_monotonically_d <<
678 const std::vector<G4double>& ally) <<
679 { <<
680 for(auto it = ally.cbegin(); it != ally.cend <<
681 { <<
682 if(*it < *(it + 1)) <<
683 { <<
684 return FALSE; <<
685 } <<
686 } <<
687 <<
688 ++noPass; <<
689 return TRUE; <<
690 } <<
691 <<
692 void G4ConvergenceTester::calc_slope_fit(const <<
693 { <<
694 // create PDF bins <<
695 G4double max = largest_scores.front(); <<
696 auto last = G4int(largest_scores.size()); <<
697 G4double min = 0.0; <<
698 if(largest_scores.back() != 0) <<
699 { <<
700 min = largest_scores.back(); <<
701 } <<
702 else <<
703 { <<
704 min = largest_scores[last - 1]; <<
705 last = last - 1; <<
706 } <<
707 <<
708 if(max * 0.99 < min) <<
709 { <<
710 // upper limit is assumed to have been rea <<
711 slope = 10.0; <<
712 return; <<
713 } <<
714 <<
715 std::vector<G4double> pdf_grid; <<
716 <<
717 pdf_grid.resize(noBinOfPDF + 1); // no grid <<
718 pdf_grid[0] = max; <<
719 pdf_grid[noBinOfPDF] = min; <<
720 G4double log10_max = std::log10(max); <<
721 G4double log10_min = std::log10(min); <<
722 G4double log10_delta = log10_max - log10_min <<
723 for(G4int i = 1; i < noBinOfPDF; ++i) <<
724 { <<
725 pdf_grid[i] = std::pow(10.0, log10_max - l <<
726 } <<
727 <<
728 std::vector<G4double> pdf; <<
729 pdf.resize(noBinOfPDF); <<
730 <<
731 for(G4int j = 0; j < last; ++j) <<
732 { <<
733 for(G4int i = 0; i < 11; ++i) <<
734 { <<
735 if(largest_scores[j] >= pdf_grid[i + 1]) <<
736 { 631 {
737 pdf[i] += 1.0 / (pdf_grid[i] - pdf_gri << 632 if ( largest_scores[j] >= pdf_grid[i+1] )
738 break; << 633 {
>> 634 pdf[i] += 1.0 / ( pdf_grid[i] - pdf_grid[i+1] ) / n;
>> 635 //G4cout << "pdf " << j << " " << i << " " << largest_scores[j] << " " << G4endl;
>> 636 break;
>> 637 }
739 } 638 }
740 } << 639 }
741 } <<
742 640
743 f_xi.resize(noBinOfPDF); << 641 f_xi.resize( noBinOfPDF );
744 f_yi.resize(noBinOfPDF); << 642 f_yi.resize( noBinOfPDF );
745 for(G4int i = 0; i < noBinOfPDF; ++i) << 643 for ( G4int i = 0 ; i < noBinOfPDF ; i++ )
746 { << 644 {
747 f_xi[i] = (pdf_grid[i] + pdf_grid[i + 1]) << 645 //G4cout << "pdf i " << i << " " << (pdf_grid[i]+pdf_grid[i+1])/2 << " " << pdf[i] << G4endl;
748 f_yi[i] = pdf[i]; << 646 f_xi[i] = (pdf_grid[i]+pdf_grid[i+1])/2;
749 } << 647 f_yi[i] = pdf[i];
750 << 648 }
751 // number of variables ( a and k ) << 649
752 minimizer = new G4SimplexDownhill<G4Converge << 650 // number of variables ( a and k )
753 // G4double minimum = minimizer->GetMinimum << 651 minimizer = new G4SimplexDownhill<G4ConvergenceTester> ( this , 2 );
754 std::vector<G4double> mp = minimizer->GetMin << 652 //G4double minimum = minimizer->GetMinimum();
755 G4double k = mp[1]; << 653 std::vector<G4double> mp = minimizer->GetMinimumPoint();
756 << 654 G4double k = mp[1];
757 // G4cout << "SLOPE " << 1/mp[1]+1 << G4endl << 655
758 // G4cout << "SLOPE a " << mp[0] << G4endl; << 656 //G4cout << "SLOPE " << 1/mp[1]+1 << G4endl;
759 // G4cout << "SLOPE k " << mp[1] << G4endl; << 657 //G4cout << "SLOPE a " << mp[0] << G4endl;
760 // G4cout << "SLOPE minimum " << minimizer- << 658 //G4cout << "SLOPE k " << mp[1] << G4endl;
761 << 659 //G4cout << "SLOPE minimum " << minimizer->GetMinimum() << G4endl;
762 slope = 1 / mp[1] + 1; << 660
763 if(k < 1.0 / 9) // Please look Pareto distr << 661 slope = 1/mp[1]+1;
764 { << 662 if ( k < 1.0/9 ) // Please look Pareto distribution with "sigma=a" and "k"
765 slope = 10; << 663 {
766 } << 664 slope = 10;
767 if(slope > 10) << 665 }
768 { << 666 if ( slope > 10 )
769 slope = 10; << 667 {
770 } << 668 slope = 10;
771 } << 669 }
772 << 670 }
773 G4double G4ConvergenceTester::slope_fitting_fu << 671
774 { << 672
775 G4double a = x[0]; << 673
776 G4double k = x[1]; << 674 G4double G4ConvergenceTester::slope_fitting_function ( std::vector< G4double > x )
777 << 675 {
778 if(a <= 0) << 676
779 { << 677 G4double a = x[0];
780 return 3.402823466e+38; // FLOAT_MAX << 678 G4double k = x[1];
781 } << 679
782 if(k == 0) << 680 if ( a <= 0 )
783 { << 681 {
784 return 3.402823466e+38; // FLOAT_MAX << 682 return 3.402823466e+38; // FLOAT_MAX
785 } << 683 }
786 << 684 if ( k == 0 )
787 // f_xi and f_yi is filled at "calc_slope_fi << 685 {
788 << 686 return 3.402823466e+38; // FLOAT_MAX
789 G4double y = 0.0; << 687 }
790 for(G4int i = 0; i < G4int(f_yi.size()); ++i << 688
791 { << 689 // f_xi and f_yi is filled at "calc_slope_fit"
792 // if ( 1/a * ( 1 + k * f_xi [ i ] / a ) < << 690
793 if((1 + k * f_xi[i] / a) < 0) << 691 G4double y = 0.0;
794 { << 692 G4int i;
795 y += 3.402823466e+38; // FLOAT_MAX << 693 for ( i = 0 ; i < int ( f_yi.size() ) ; i++ )
796 } << 694 {
797 else << 695 //if ( 1/a * ( 1 + k * f_xi [ i ] / a ) < 0 )
798 { << 696 if ( ( 1 + k * f_xi [ i ] / a ) < 0 )
799 y += (f_yi[i] - 1 / a * std::pow(1 + k * << 697 {
800 (f_yi[i] - 1 / a * std::pow(1 + k * << 698 y +=3.402823466e+38; // FLOAT_MAX
801 } << 699 }
802 } << 700 else
>> 701 {
>> 702 y += ( f_yi [ i ] - 1/a*std::pow ( 1 + k * f_xi [ i ] / a , - 1/k - 1 ) ) * ( f_yi [ i ] - 1/a*std::pow ( 1 + k * f_xi [ i ] / a , - 1/k - 1 ) );
>> 703 }
>> 704 }
>> 705 // G4cout << "y = " << y << G4endl;
803 706
804 return y; << 707 return y;
805 } 708 }
806 709