| 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(); << 87
107 } << 88 //G4cout << x << G4endl;
108 } << 89
109 sum += x; << 90 timer->Stop();
110 } << 91 cpu_time.push_back( timer->GetSystemElapsed() + timer->GetUserElapsed() );
111 << 92
112 // Data has been added so statistics have no << 93 if ( x < 0.0 ) {
113 statsAreUpdated = false; << 94 G4cout << "Warning: G4convergenceTester expects zero or positive number as inputs, but received a negative number." << G4endl;
114 ++n; << 95 }
115 l.unlock(); << 96
116 return; << 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
>> 136 efficiency = double( nonzero_histories.size() ) / n;
122 137
123 efficiency = G4double(nonzero_histories.size << 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 }
124 192
125 mean = sum / n; << 193 mean_1 = 0.0;
>> 194 var_1 = 0.0;
>> 195 shift_1 = 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 }
126 252
127 G4double sum_x2 = 0.0; << 253 calc_grid_point_of_history();
128 var = 0.0; << 254 calc_stat_history();
129 shift = 0.0; << 255
130 vov = 0.0; << 256 // statistics have been calculated and this function does not need
131 << 257 // to be called again until data has been added
132 G4double xi; << 258 statsAreUpdated = true;
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 } 259 }
246 260
>> 261
>> 262
247 void G4ConvergenceTester::calc_grid_point_of_h 263 void G4ConvergenceTester::calc_grid_point_of_history()
248 { 264 {
249 // histroy_grid [ 0,,,15 ] << 265
250 // history_grid [0] 1/16 ,,, history_grid [1 << 266 // histroy_grid [ 0,,,15 ]
251 // if number of event is x then history_grid << 267 // history_grid [0] 1/16 ,,, history_grid [15] 16/16
252 // 16 -> noBinOfHisotry << 268 // if number of event is x then history_grid [15] become x-1.
253 << 269 // 16 -> noBinOfHisotry
254 for(G4int i = 1; i <= noBinOfHistory; ++i) << 270
255 { << 271 G4int i;
256 history_grid[i - 1] = G4int(n / (G4double( << 272 for ( i = 1 ; i <= noBinOfHistory ; i++ )
257 } << 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 }
>> 277
258 } 278 }
259 279
>> 280
>> 281
260 void G4ConvergenceTester::calc_stat_history() 282 void G4ConvergenceTester::calc_stat_history()
261 { 283 {
262 if(history_grid[0] == 0) << 284 // G4cout << "i/16 till_ith mean var sd r vov fom shift e r2eff r2int" << G4endl;
263 { << 285
264 showHistory = false; << 286 if ( history_grid [ 0 ] == 0 ) {
265 return; << 287 showHistory=false;
266 } << 288 return;
267 << 289 }
268 for(G4int i = 0; i < noBinOfHistory; ++i) << 290
269 { << 291 G4int i;
270 G4int ith = history_grid[i]; << 292 for ( i = 1 ; i <= noBinOfHistory ; i++ )
271 << 293 {
272 G4int nonzero_till_ith = 0; << 294
273 G4double xi; << 295 G4int ith = history_grid [ i-1 ];
274 G4double mean_till_ith = 0.0; << 296
275 << 297 G4int nonzero_till_ith = 0;
276 for(const auto& itr : nonzero_histories) << 298 G4double xi;
277 { << 299 G4double mean_till_ith = 0.0;
278 if(itr.first <= ith) << 300 std::map< G4int , G4double >::iterator it;
>> 301
>> 302 for ( it = nonzero_histories.begin() ; it !=nonzero_histories.end() ; it++ )
279 { 303 {
280 xi = itr.second; << 304 if ( it->first <= ith )
281 mean_till_ith += xi; << 305 {
282 ++nonzero_till_ith; << 306 xi = it->second;
>> 307 mean_till_ith += xi;
>> 308 nonzero_till_ith++;
>> 309 }
283 } 310 }
284 } <<
285 311
286 if(nonzero_till_ith == 0) << 312 if ( nonzero_till_ith == 0 ) continue;
287 { << 313
288 continue; << 314 mean_till_ith = mean_till_ith / ( ith+1 );
289 } << 315 mean_history [ i-1 ] = mean_till_ith;
290 << 316
291 mean_till_ith = mean_till_ith / (ith + 1 << 317 G4double sum_x2_till_ith = 0.0;
292 mean_history[i] = mean_till_ith; << 318 G4double var_till_ith = 0.0;
293 << 319 G4double vov_till_ith = 0.0;
294 G4double sum_x2_till_ith = 0.0; << 320 G4double shift_till_ith = 0.0;
295 G4double var_till_ith = 0.0; << 321
296 G4double vov_till_ith = 0.0; << 322 for ( it = nonzero_histories.begin() ; it !=nonzero_histories.end() ; it++ )
297 G4double shift_till_ith = 0.0; <<
298 <<
299 for(const auto& itr : nonzero_histories) <<
300 { <<
301 if(itr.first <= ith) <<
302 { 323 {
303 xi = itr.second; << 324 if ( it->first <= ith )
304 sum_x2_till_ith += std::pow(xi, 2.0); << 325 {
305 var_till_ith += std::pow(xi - mean_til << 326 xi = it->second;
306 shift_till_ith += std::pow(xi - mean_t << 327 sum_x2_till_ith += xi * xi;
307 vov_till_ith += std::pow(xi - mean_til << 328 var_till_ith += ( xi - mean_till_ith ) * ( xi - mean_till_ith );
>> 329 shift_till_ith += ( xi - mean_till_ith ) * ( xi - mean_till_ith ) * ( xi - mean_till_ith );
>> 330 vov_till_ith += ( xi - mean_till_ith ) * ( xi - mean_till_ith ) * ( xi - mean_till_ith ) * ( xi - mean_till_ith );
>> 331 }
308 } 332 }
309 } <<
310 333
311 var_till_ith += << 334 var_till_ith += ( (ith+1) - nonzero_till_ith ) * mean_till_ith * mean_till_ith;
312 ((ith + 1) - nonzero_till_ith) * std::po << 335 vov_till_ith += ( (ith+1) - nonzero_till_ith ) * mean_till_ith * mean_till_ith * mean_till_ith * mean_till_ith ;
313 vov_till_ith += << 336
314 ((ith + 1) - nonzero_till_ith) * std::po << 337
315 << 338 if ( var_till_ith == 0 ) continue;
316 G4double sum_till_ith = mean_till_ith * (i << 339 vov_till_ith = vov_till_ith / ( var_till_ith * var_till_ith ) - 1.0 / (ith+1);
317 << 340 vov_history [ i-1 ] = vov_till_ith;
318 if(!(std::fabs(var_till_ith) > 0.0)) << 341
319 { << 342 var_till_ith = var_till_ith / ( ith+1 - 1 );
320 continue; << 343 var_history [ i-1 ] = var_till_ith;
321 } << 344
322 if(!(std::fabs(mean_till_ith) > 0.0)) << 345 sd_history [ i-1 ] = std::sqrt( var_till_ith );
323 { << 346 r_history [ i-1 ] = std::sqrt( var_till_ith ) / mean_till_ith / std::sqrt ( 1.0*(ith+1) );
324 continue; << 347
325 } << 348 fom_history [ i-1 ] = 1 / ( r_history [ i-1 ] * r_history [ i-1 ] ) / cpu_time [ ith ];
326 if(!(std::fabs(sum_till_ith) > 0.0)) << 349
327 { << 350 shift_till_ith += ( (ith+1) - nonzero_till_ith ) * mean_till_ith * mean_till_ith * mean_till_ith * ( -1 );
328 continue; << 351 shift_till_ith = shift_till_ith / ( 2 * var_till_ith * (ith+1) );
329 } << 352 shift_history [ i-1 ] = shift_till_ith;
330 << 353
331 vov_till_ith = vov_till_ith / std::pow(var << 354 e_history [ i-1 ] = 1.0*nonzero_till_ith / (ith+1);
332 vov_history[i] = vov_till_ith; << 355 r2eff_history [ i-1 ] = ( 1 - e_history [ i-1 ] ) / ( e_history [ i-1 ] * (ith+1) );
333 << 356
334 var_till_ith = var_till_ith / (ith + 1 - << 357 G4double sum_till_ith = mean_till_ith * (ith+1);
335 var_history[i] = var_till_ith; << 358 r2int_history [ i-1 ] = ( sum_x2_till_ith ) / ( sum_till_ith * sum_till_ith ) - 1 / ( e_history [ i-1 ] * (ith+1) );
336 sd_history[i] = std::sqrt(var_till_ith); << 359
337 r_history[i] = << 360 }
338 std::sqrt(var_till_ith) / mean_till_ith << 361
339 <<
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 } <<
348 <<
349 shift_till_ith += <<
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 <<
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 <<
359 r2int_history[i] = (sum_x2_till_ith) / s <<
360 1 / (e_history[i] * ( <<
361 } <<
362 } <<
363 } 362 }
364 363
>> 364
>> 365
365 void G4ConvergenceTester::ShowResult(std::ostr 366 void G4ConvergenceTester::ShowResult(std::ostream& out)
366 { 367 {
367 // if data has been added since the last com << 368 // if data has been added since the last computation of the statistical values (not statsAreUpdated)
368 // (not statsAreUpdated) call calStat to rec << 369 // call calStat to recompute the statistical values
369 if(!statsAreUpdated) << 370 if(!statsAreUpdated) { calStat(); }
370 { << 371
371 calStat(); << 372 out << std::setprecision( 6 );
372 } << 373
373 << 374 out << G4endl;
374 out << std::setprecision(6); << 375 out << "G4ConvergenceTester Output Result of " << name << G4endl;
375 << 376 out << std::setw(20) << "EFFICIENCY = " << std::setw(13) << efficiency << G4endl;
376 out << G4endl; << 377 out << std::setw(20) << "MEAN = " << std::setw(13) << mean << G4endl;
377 out << "G4ConvergenceTester Output Result of << 378 out << std::setw(20) << "VAR = " << std::setw(13) << var << G4endl;
378 out << std::setw(20) << "EFFICIENCY = " << s << 379 out << std::setw(20) << "SD = " << std::setw(13) << sd << G4endl;
379 << G4endl; << 380 out << std::setw(20) << "R = " << std::setw(13) << r << G4endl;
380 out << std::setw(20) << "MEAN = " << std::se << 381 out << std::setw(20) << "SHIFT = "<< std::setw(13) << shift << G4endl;
381 out << std::setw(20) << "VAR = " << std::set << 382 out << std::setw(20) << "VOV = "<< std::setw(13) << vov << G4endl;
382 out << std::setw(20) << "SD = " << std::setw << 383 out << std::setw(20) << "FOM = "<< std::setw(13) << fom << G4endl;
383 out << std::setw(20) << "R = " << std::setw( << 384
384 out << std::setw(20) << "SHIFT = " << std::s << 385 out << std::setw(20) << "THE LARGEST SCORE = " << std::setw(13) << largest << " and it happend at " << largest_score_happened << "th event" << G4endl;
385 out << std::setw(20) << "VOV = " << std::set << 386 if ( mean!=0 ) {
386 out << std::setw(20) << "FOM = " << std::set << 387 out << std::setw(20) << "Affected Mean = " << std::setw(13) << mean_1 << " and its ratio to orignal is " << mean_1/mean << G4endl;
387 << 388 } else {
388 out << std::setw(20) << "THE LARGEST SCORE = << 389 out << std::setw(20) << "Affected Mean = " << std::setw(13) << mean_1 << G4endl;
389 << " and it happened at " << largest_sco << 390 }
390 << G4endl; << 391 if ( var!=0 ) {
391 if(mean != 0) << 392 out << std::setw(20) << "Affected VAR = " << std::setw(13) << var_1 << " and its ratio to orignal is " << var_1/var << G4endl;
392 { << 393 } else {
393 out << std::setw(20) << "Affected Mean = " << 394 out << std::setw(20) << "Affected VAR = " << std::setw(13) << var_1 << G4endl;
394 << " and its ratio to original is " << << 395 }
395 } << 396 if ( r!=0 ) {
396 else << 397 out << std::setw(20) << "Affected R = " << std::setw(13) << r_1 << " and its ratio to orignal is " << r_1/r << G4endl;
397 { << 398 } else {
398 out << std::setw(20) << "Affected Mean = " << 399 out << std::setw(20) << "Affected R = " << std::setw(13) << r_1 << G4endl;
399 << G4endl; << 400 }
400 } << 401 if ( shift!=0 ) {
401 if(var != 0) << 402 out << std::setw(20) << "Affected SHIFT = " << std::setw(13) << shift_1 << " and its ratio to orignal is " << shift_1/shift << G4endl;
402 { << 403 } else {
403 out << std::setw(20) << "Affected VAR = " << 404 out << std::setw(20) << "Affected SHIFT = " << std::setw(13) << shift_1 << G4endl;
404 << " and its ratio to original is " << << 405 }
405 } << 406 if ( fom!=0 ) {
406 else << 407 out << std::setw(20) << "Affected FOM = " << std::setw(13) << fom_1 << " and its ratio to orignal is " << fom_1/fom << G4endl;
407 { << 408 } else {
408 out << std::setw(20) << "Affected VAR = " << 409 out << std::setw(20) << "Affected FOM = " << std::setw(13) << fom_1 << G4endl;
409 << G4endl; << 410 }
410 } << 411
411 if(r != 0) << 412 if ( !showHistory ) {
412 { << 413 out << "Number of events of this run is too small to do convergence tests." << G4endl;
413 out << std::setw(20) << "Affected R = " << << 414 return;
414 << " and its ratio to original is " << << 415 }
415 } << 416
416 else << 417 check_stat_history(out);
417 { << 418
418 out << std::setw(20) << "Affected R = " << << 419 // check SLOPE and output result
419 } << 420 if ( calcSLOPE ) {
420 if(shift != 0) << 421 if ( slope >= 3 )
421 { << 422 {
422 out << std::setw(20) << "Affected SHIFT = << 423 noPass++;
423 << " and its ratio to original is " << << 424 out << "SLOPE is large enough" << G4endl;
424 } << 425 }
425 else << 426 else
426 { << 427 {
427 out << std::setw(20) << "Affected SHIFT = << 428 out << "SLOPE is not large enough" << G4endl;
428 << G4endl; << 429 }
429 } << 430 } else {
430 if(fom != 0) << 431 out << "Number of non zero history too small to calculate SLOPE" << G4endl;
431 { << 432 }
432 out << std::setw(20) << "Affected FOM = " << 433
433 << " and its ratio to original is " << << 434 out << "This result passes " << noPass << " / "<< noTotal << " Convergence Test." << G4endl;
434 } << 435 out << G4endl;
435 else << 436
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 } 437 }
472 438
473 void G4ConvergenceTester::ShowHistory(std::ost 439 void G4ConvergenceTester::ShowHistory(std::ostream& out)
474 { 440 {
475 if(!showHistory) << 441
476 { << 442 if ( !showHistory ) {
477 out << "Number of events of this run is to << 443 out << "Number of events of this run is too small to show history." << G4endl;
478 << G4endl; << 444 return;
479 return; << 445 }
480 } << 446
481 << 447 out << std::setprecision( 6 );
482 out << std::setprecision(6); << 448
483 << 449 out << G4endl;
484 out << G4endl; << 450 out << "G4ConvergenceTester Output History of " << name << G4endl;
485 out << "G4ConvergenceTester Output History o << 451 out << "i/" << noBinOfHistory << " till_ith mean"
486 out << "i/" << noBinOfHistory << " till_ith << 452 << std::setw(13) << "var"
487 << "var" << std::setw(13) << "sd" << std << 453 << std::setw(13) << "sd"
488 << "vov" << std::setw(13) << "fom" << st << 454 << std::setw(13) << "r"
489 << std::setw(13) << "e" << std::setw(13) << 455 << std::setw(13) << "vov"
490 << "r2int" << G4endl; << 456 << std::setw(13) << "fom"
491 for(G4int i = 1; i <= noBinOfHistory; i++) << 457 << std::setw(13) << "shift"
492 { << 458 << std::setw(13) << "e"
493 out << std::setw(4) << i << " " << std::se << 459 << std::setw(13) << "r2eff"
494 << std::setw(13) << mean_history[i - 1 << 460 << std::setw(13) << "r2int"
495 << var_history[i - 1] << std::setw(13) << 461 << G4endl;
496 << std::setw(13) << r_history[i - 1] < << 462 for ( G4int i = 1 ; i <= noBinOfHistory ; i++ )
497 << vov_history[i - 1] << std::setw(13) << 463 {
498 << std::setw(13) << shift_history[i - << 464 out << std::setw( 4) << i << " "
499 << e_history[i - 1] << std::setw(13) < << 465 << std::setw( 5) << history_grid [ i-1 ]
500 << std::setw(13) << r2int_history[i - << 466 << std::setw(13) << mean_history [ i-1 ]
501 } << 467 << std::setw(13) << var_history [ i-1 ]
>> 468 << std::setw(13) << sd_history [ i-1 ]
>> 469 << std::setw(13) << r_history [ i-1 ]
>> 470 << std::setw(13) << vov_history [ i-1 ]
>> 471 << std::setw(13) << fom_history [ i-1 ]
>> 472 << std::setw(13) << shift_history [ i-1 ]
>> 473 << std::setw(13) << e_history [ i-1 ]
>> 474 << std::setw(13) << r2eff_history [ i-1 ]
>> 475 << std::setw(13) << r2int_history [ i-1 ]
>> 476 << G4endl;
>> 477 }
502 } 478 }
503 479
504 void G4ConvergenceTester::check_stat_history(s 480 void G4ConvergenceTester::check_stat_history(std::ostream& out)
505 { 481 {
506 // 1 sigma rejection for null hypothesis <<
507 482
508 std::vector<G4double> first_ally; << 483 // 1 sigma rejection for null hypothesis
509 std::vector<G4double> second_ally; <<
510 484
511 // use 2nd half of hisories << 485 std::vector<G4double> first_ally;
512 std::size_t N = mean_history.size() / 2; << 486 std::vector<G4double> second_ally;
513 std::size_t i; << 487
514 << 488 // use 2nd half of hisories
515 G4double pearson_r; << 489 G4int N = mean_history.size() / 2;
516 G4double t; << 490 G4int i;
517 << 491
518 first_ally.resize(N); << 492 G4double pearson_r;
519 second_ally.resize(N); << 493 G4double t;
520 << 494
521 G4double sum_of_var = << 495 first_ally.resize( N );
522 std::accumulate(var_history.begin(), var_h << 496 second_ally.resize( N );
523 if(sum_of_var == 0.0) << 497
524 { << 498 //
525 out << "Variances in all historical grids << 499 G4double sum_of_var = std::accumulate ( var_history.begin() , var_history.end() , 0.0 );
526 out << "Terminating checking behavior of s << 500 if ( sum_of_var == 0.0 ) {
527 return; << 501 out << "Variances in all historical grids are zero." << G4endl;
528 } << 502 out << "Terminating checking behavior of statistics numbers." << G4endl;
529 << 503 return;
530 // Mean << 504 }
531 << 505
532 for(i = 0; i < N; ++i) << 506 // Mean
533 { << 507
534 first_ally[i] = history_grid[N + i]; << 508 for ( i = 0 ; i < N ; i++ )
535 second_ally[i] = mean_history[N + i]; << 509 {
536 } << 510 first_ally [ i ] = history_grid [ N + i ];
537 << 511 second_ally [ i ] = mean_history [ N + i ];
538 pearson_r = calc_Pearson_r((G4int)N, first_a << 512 }
539 t = pearson_r * std::sqrt((N - 2) / << 513
540 << 514 pearson_r = calc_Pearson_r ( N , first_ally , second_ally );
541 if(t < 0.429318) // Student t of (Degree of << 515 t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) );
542 { << 516
543 out << "MEAN distribution is RANDOM" << G << 517 if ( t < 0.429318 ) // Student t of (Degree of freedom = N-2 )
544 noPass++; << 518 {
545 } << 519 out << "MEAN distribution is RANDOM" << G4endl;
546 else << 520 noPass++;
547 { << 521 }
548 out << "MEAN distribution is not RANDOM" < << 522 else
549 } << 523 {
550 << 524 out << "MEAN distribution is not RANDOM" << G4endl;
551 // R << 525 }
552 << 526
553 for(i = 0; i < N; ++i) << 527
554 { << 528 // R
555 first_ally[i] = 1.0 / std::sqrt(G4double( << 529
556 second_ally[i] = r_history[N + i]; << 530 for ( i = 0 ; i < N ; i++ )
557 } << 531 {
558 << 532 first_ally [ i ] = 1.0 / std::sqrt ( G4double(history_grid [ N + i ]) );
559 pearson_r = calc_Pearson_r(G4int(N), first_a << 533 second_ally [ i ] = r_history [ N + i ];
560 t = pearson_r * std::sqrt((N - 2) / << 534 }
561 << 535
562 if(t > 1.090546) << 536 pearson_r = calc_Pearson_r ( N , first_ally , second_ally );
563 { << 537 t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) );
564 out << "r follows 1/std::sqrt(N)" << G4end << 538
565 noPass++; << 539 if ( t > 1.090546 )
566 } << 540 {
567 else << 541 out << "r follows 1/std::sqrt(N)" << G4endl;
568 { << 542 noPass++;
569 out << "r does not follow 1/std::sqrt(N)" << 543 }
570 } << 544 else
571 << 545 {
572 if(is_monotonically_decrease(second_ally)) << 546 out << "r does not follow 1/std::sqrt(N)" << G4endl;
573 { << 547 }
574 out << "r is monotonically decrease " << G << 548
575 } << 549 if ( is_monotonically_decrease( second_ally ) == true )
576 else << 550 {
577 { << 551 out << "r is monotonically decrease " << G4endl;
578 out << "r is NOT monotonically decrease " << 552 }
579 } << 553 else
580 << 554 {
581 if(r_history.back() < 0.1) << 555 out << "r is NOT monotonically decrease " << G4endl;
582 { << 556 }
583 out << "r is less than 0.1. r = " << r_his << 557
584 noPass++; << 558 if ( r_history.back() < 0.1 )
585 } << 559 {
586 else << 560 out << "r is less than 0.1. r = " << r_history.back() << G4endl;
587 { << 561 noPass++;
588 out << "r is NOT less than 0.1. r = " << r << 562 }
589 } << 563 else
590 << 564 {
591 // VOV << 565 out << "r is NOT less than 0.1. r = " << r_history.back() << G4endl;
592 for(i = 0; i < N; ++i) << 566 }
593 { << 567
594 first_ally[i] = 1.0 / history_grid[N + i] << 568
595 second_ally[i] = vov_history[N + i]; << 569 // VOV
596 } << 570 for ( i = 0 ; i < N ; i++ )
597 << 571 {
598 pearson_r = calc_Pearson_r(G4int(N), first_a << 572 first_ally [ i ] = 1.0 / history_grid [ N + i ];
599 t = pearson_r * std::sqrt((N - 2) / << 573 second_ally [ i ] = vov_history [ N + i ];
600 << 574 }
601 if(t > 1.090546) << 575
602 { << 576 pearson_r = calc_Pearson_r ( N , first_ally , second_ally );
603 out << "VOV follows 1/std::sqrt(N)" << G4e << 577 t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) );
604 noPass++; << 578
605 } << 579 if ( t > 1.090546 )
606 else << 580 {
607 { << 581 out << "VOV follows 1/std::sqrt(N)" << G4endl;
608 out << "VOV does not follow 1/std::sqrt(N) << 582 noPass++;
609 } << 583 }
610 << 584 else
611 if(is_monotonically_decrease(second_ally)) << 585 {
612 { << 586 out << "VOV does not follow 1/std::sqrt(N)" << G4endl;
613 out << "VOV is monotonically decrease " << << 587 }
614 } << 588
615 else << 589 if ( is_monotonically_decrease( second_ally ) == true )
616 { << 590 {
617 out << "VOV is NOT monotonically decrease << 591 out << "VOV is monotonically decrease " << G4endl;
618 } << 592 }
619 << 593 else
620 // FOM << 594 {
621 << 595 out << "VOV is NOT monotonically decrease " << G4endl;
622 for(i = 0; i < N; ++i) << 596 }
623 { << 597
624 first_ally[i] = history_grid[N + i]; << 598 // FOM
625 second_ally[i] = fom_history[N + i]; << 599
626 } << 600 for ( i = 0 ; i < N ; i++ )
627 << 601 {
628 pearson_r = calc_Pearson_r(G4int(N), std::mo << 602 first_ally [ i ] = history_grid [ N + i ];
629 t = pearson_r * std::sqrt((N - 2) / << 603 second_ally [ i ] = fom_history [ N + i ];
630 << 604 }
631 if(t < 0.429318) << 605
632 { << 606 pearson_r = calc_Pearson_r ( N , first_ally , second_ally );
633 out << "FOM distribution is RANDOM" << G4e << 607 t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) );
634 noPass++; << 608
635 } << 609 if ( t < 0.429318 )
636 else << 610 {
637 { << 611 out << "FOM distribution is RANDOM" << G4endl;
638 out << "FOM distribution is not RANDOM" << << 612 noPass++;
639 } << 613 }
640 } << 614 else
641 << 615 {
642 G4double G4ConvergenceTester::calc_Pearson_r(G << 616 out << "FOM distribution is not RANDOM" << G4endl;
643 s << 617 }
644 s << 618
645 { << 619 }
646 G4double first_mean = 0.0; << 620
647 G4double second_mean = 0.0; << 621
648 << 622
649 G4int i; << 623 G4double G4ConvergenceTester::calc_Pearson_r ( G4int N , std::vector<G4double> first_ally , std::vector<G4double> second_ally )
650 for(i = 0; i < N; i++) << 624 {
651 { << 625 G4double first_mean = 0.0;
652 first_mean += first_ally[i]; << 626 G4double second_mean = 0.0;
653 second_mean += second_ally[i]; << 627
654 } << 628 G4int i;
655 first_mean = first_mean / N; << 629 for ( i = 0 ; i < N ; i++ )
656 second_mean = second_mean / N; << 630 {
657 << 631 first_mean += first_ally [ i ];
658 G4double a = 0.0; << 632 second_mean += second_ally [ i ];
659 for(i = 0; i < N; ++i) << 633 }
660 { << 634 first_mean = first_mean / N;
661 a += (first_ally[i] - first_mean) * (secon << 635 second_mean = second_mean / N;
662 } << 636
663 << 637 G4double a = 0.0;
664 G4double b1 = 0.0; << 638 for ( i = 0 ; i < N ; i++ )
665 G4double b2 = 0.0; << 639 {
666 for(i = 0; i < N; ++i) << 640 a += ( first_ally [ i ] - first_mean ) * ( second_ally [ i ] - second_mean );
667 { << 641 }
668 b1 += (first_ally[i] - first_mean) * (firs << 642
669 b2 += (second_ally[i] - second_mean) * (se << 643 G4double b1 = 0.0;
670 } << 644 G4double b2 = 0.0;
671 << 645 for ( i = 0 ; i < N ; i++ )
672 G4double rds = a / std::sqrt(b1 * b2); << 646 {
673 << 647 b1 += ( first_ally [ i ] - first_mean ) * ( first_ally [ i ] - first_mean );
674 return rds; << 648 b2 += ( second_ally [ i ] - second_mean ) * ( second_ally [ i ] - second_mean );
675 } << 649 }
676 << 650
677 G4bool G4ConvergenceTester::is_monotonically_d << 651 G4double rds = a / std::sqrt ( b1 * b2 );
678 const std::vector<G4double>& ally) << 652
679 { << 653 return rds;
680 for(auto it = ally.cbegin(); it != ally.cend << 654 }
681 { << 655
682 if(*it < *(it + 1)) << 656
683 { << 657
684 return FALSE; << 658 G4bool G4ConvergenceTester::is_monotonically_decrease ( std::vector<G4double> ally )
685 } << 659 {
686 } << 660
687 << 661 std::vector<G4double>::iterator it;
688 ++noPass; << 662 for ( it = ally.begin() ; it != ally.end() - 1 ; it++ )
689 return TRUE; << 663 {
690 } << 664 if ( *it < *(it+1) ) return FALSE;
691 << 665 }
692 void G4ConvergenceTester::calc_slope_fit(const << 666
693 { << 667 noPass++;
694 // create PDF bins << 668 return TRUE;
695 G4double max = largest_scores.front(); << 669 }
696 auto last = G4int(largest_scores.size()); << 670
697 G4double min = 0.0; << 671
698 if(largest_scores.back() != 0) << 672
699 { << 673 //void G4ConvergenceTester::calc_slope_fit ( std::vector<G4double> largest_socres )
700 min = largest_scores.back(); << 674 void G4ConvergenceTester::calc_slope_fit ( std::vector<G4double> )
701 } << 675 {
702 else << 676
703 { << 677 // create PDF bins
704 min = largest_scores[last - 1]; << 678 G4double max = largest_scores.front();
705 last = last - 1; << 679 G4int last = int ( largest_scores.size() );
706 } << 680 G4double min = 0.0;
707 << 681 if ( largest_scores.back() != 0 )
708 if(max * 0.99 < min) << 682 {
709 { << 683 min = largest_scores.back();
710 // upper limit is assumed to have been rea << 684 }
711 slope = 10.0; << 685 else
712 return; << 686 {
713 } << 687 min = largest_scores[ last-1 ];
714 << 688 last = last - 1;
715 std::vector<G4double> pdf_grid; << 689 }
716 << 690
717 pdf_grid.resize(noBinOfPDF + 1); // no grid << 691 //G4cout << "largest " << max << G4endl;
718 pdf_grid[0] = max; << 692 //G4cout << "last " << min << G4endl;
719 pdf_grid[noBinOfPDF] = min; << 693
720 G4double log10_max = std::log10(max); << 694 if ( max*0.99 < min )
721 G4double log10_min = std::log10(min); << 695 {
722 G4double log10_delta = log10_max - log10_min << 696 // upper limit is assumed to have been reached
723 for(G4int i = 1; i < noBinOfPDF; ++i) << 697 slope = 10.0;
724 { << 698 return;
725 pdf_grid[i] = std::pow(10.0, log10_max - l << 699 }
726 } << 700
727 << 701 std::vector < G4double > pdf_grid;
728 std::vector<G4double> pdf; << 702
729 pdf.resize(noBinOfPDF); << 703 pdf_grid.resize( noBinOfPDF+1 ); // no grid = no bins + 1
730 << 704 pdf_grid[ 0 ] = max;
731 for(G4int j = 0; j < last; ++j) << 705 pdf_grid[ noBinOfPDF ] = min;
732 { << 706 G4double log10_max = std::log10( max );
733 for(G4int i = 0; i < 11; ++i) << 707 G4double log10_min = std::log10( min );
734 { << 708 G4double log10_delta = log10_max - log10_min;
735 if(largest_scores[j] >= pdf_grid[i + 1]) << 709 for ( G4int i = 1 ; i < noBinOfPDF ; i++ )
>> 710 {
>> 711 pdf_grid[i] = std::pow ( 10.0 , log10_max - log10_delta/10.0*(i) );
>> 712 //G4cout << "pdf i " << i << " " << pdf_grid[i] << G4endl;
>> 713 }
>> 714
>> 715 std::vector < G4double > pdf;
>> 716 pdf.resize( noBinOfPDF );
>> 717
>> 718 for ( G4int j=0 ; j < last ; j ++ )
>> 719 {
>> 720 for ( G4int i = 0 ; i < 11 ; i++ )
736 { 721 {
737 pdf[i] += 1.0 / (pdf_grid[i] - pdf_gri << 722 if ( largest_scores[j] >= pdf_grid[i+1] )
738 break; << 723 {
>> 724 pdf[i] += 1.0 / ( pdf_grid[i] - pdf_grid[i+1] ) / n;
>> 725 //G4cout << "pdf " << j << " " << i << " " << largest_scores[j] << " " << G4endl;
>> 726 break;
>> 727 }
739 } 728 }
740 } << 729 }
741 } <<
742 730
743 f_xi.resize(noBinOfPDF); << 731 f_xi.resize( noBinOfPDF );
744 f_yi.resize(noBinOfPDF); << 732 f_yi.resize( noBinOfPDF );
745 for(G4int i = 0; i < noBinOfPDF; ++i) << 733 for ( G4int i = 0 ; i < noBinOfPDF ; i++ )
746 { << 734 {
747 f_xi[i] = (pdf_grid[i] + pdf_grid[i + 1]) << 735 //G4cout << "pdf i " << i << " " << (pdf_grid[i]+pdf_grid[i+1])/2 << " " << pdf[i] << G4endl;
748 f_yi[i] = pdf[i]; << 736 f_xi[i] = (pdf_grid[i]+pdf_grid[i+1])/2;
749 } << 737 f_yi[i] = pdf[i];
750 << 738 }
751 // number of variables ( a and k ) << 739
752 minimizer = new G4SimplexDownhill<G4Converge << 740 // number of variables ( a and k )
753 // G4double minimum = minimizer->GetMinimum << 741 minimizer = new G4SimplexDownhill<G4ConvergenceTester> ( this , 2 );
754 std::vector<G4double> mp = minimizer->GetMin << 742 //G4double minimum = minimizer->GetMinimum();
755 G4double k = mp[1]; << 743 std::vector<G4double> mp = minimizer->GetMinimumPoint();
756 << 744 G4double k = mp[1];
757 // G4cout << "SLOPE " << 1/mp[1]+1 << G4endl << 745
758 // G4cout << "SLOPE a " << mp[0] << G4endl; << 746 //G4cout << "SLOPE " << 1/mp[1]+1 << G4endl;
759 // G4cout << "SLOPE k " << mp[1] << G4endl; << 747 //G4cout << "SLOPE a " << mp[0] << G4endl;
760 // G4cout << "SLOPE minimum " << minimizer- << 748 //G4cout << "SLOPE k " << mp[1] << G4endl;
761 << 749 //G4cout << "SLOPE minimum " << minimizer->GetMinimum() << G4endl;
762 slope = 1 / mp[1] + 1; << 750
763 if(k < 1.0 / 9) // Please look Pareto distr << 751 slope = 1/mp[1]+1;
764 { << 752 if ( k < 1.0/9 ) // Please look Pareto distribution with "sigma=a" and "k"
765 slope = 10; << 753 {
766 } << 754 slope = 10;
767 if(slope > 10) << 755 }
768 { << 756 if ( slope > 10 )
769 slope = 10; << 757 {
770 } << 758 slope = 10;
771 } << 759 }
772 << 760 }
773 G4double G4ConvergenceTester::slope_fitting_fu << 761
774 { << 762
775 G4double a = x[0]; << 763
776 G4double k = x[1]; << 764 G4double G4ConvergenceTester::slope_fitting_function ( std::vector< G4double > x )
777 << 765 {
778 if(a <= 0) << 766
779 { << 767 G4double a = x[0];
780 return 3.402823466e+38; // FLOAT_MAX << 768 G4double k = x[1];
781 } << 769
782 if(k == 0) << 770 if ( a <= 0 )
783 { << 771 {
784 return 3.402823466e+38; // FLOAT_MAX << 772 return 3.402823466e+38; // FLOAT_MAX
785 } << 773 }
786 << 774 if ( k == 0 )
787 // f_xi and f_yi is filled at "calc_slope_fi << 775 {
788 << 776 return 3.402823466e+38; // FLOAT_MAX
789 G4double y = 0.0; << 777 }
790 for(G4int i = 0; i < G4int(f_yi.size()); ++i << 778
791 { << 779 // f_xi and f_yi is filled at "calc_slope_fit"
792 // if ( 1/a * ( 1 + k * f_xi [ i ] / a ) < << 780
793 if((1 + k * f_xi[i] / a) < 0) << 781 G4double y = 0.0;
794 { << 782 G4int i;
795 y += 3.402823466e+38; // FLOAT_MAX << 783 for ( i = 0 ; i < int ( f_yi.size() ) ; i++ )
796 } << 784 {
797 else << 785 //if ( 1/a * ( 1 + k * f_xi [ i ] / a ) < 0 )
798 { << 786 if ( ( 1 + k * f_xi [ i ] / a ) < 0 )
799 y += (f_yi[i] - 1 / a * std::pow(1 + k * << 787 {
800 (f_yi[i] - 1 / a * std::pow(1 + k * << 788 y +=3.402823466e+38; // FLOAT_MAX
801 } << 789 }
802 } << 790 else
>> 791 {
>> 792 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 ) );
>> 793 }
>> 794 }
>> 795 // G4cout << "y = " << y << G4endl;
803 796
804 return y; << 797 return y;
805 } 798 }
806 799