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