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