| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4ConvergenceTester class implementation
27 //
28 // Author: Koi, Tatsumi (SLAC/SCCS)
29 // --------------------------------------------------------------------
30
31 #include "G4ConvergenceTester.hh"
32 #include "G4AutoLock.hh"
33 #include <iomanip>
34
35 namespace
36 {
37 G4Mutex aMutex = G4MUTEX_INITIALIZER;
38 }
39
40 G4ConvergenceTester::G4ConvergenceTester(const G4String& theName)
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
65 G4ConvergenceTester::~G4ConvergenceTester()
66 {
67 delete timer;
68 }
69
70 void G4ConvergenceTester::AddScore(G4double x)
71 {
72 G4AutoLock l(&aMutex);
73
74 timer->Stop();
75 cpu_time.push_back(timer->GetSystemElapsed() + timer->GetUserElapsed());
76
77 if(x < 0.0)
78 {
79 std::ostringstream message;
80 message << "Expecting zero or positive number as inputs,\n"
81 << "but received a negative number.";
82 G4Exception("G4ConvergenceTester::AddScore()", "Warning",
83 JustWarning, message);
84 }
85
86 if(x == 0.0)
87 {
88 }
89 else
90 {
91 nonzero_histories.insert(std::pair<G4int, G4double>(n, x));
92 if(x > largest_scores.back())
93 {
94 // Following search should become faster if begin from bottom.
95 for(auto it = largest_scores.begin(); it != largest_scores.end(); ++it)
96 {
97 if(x > *it)
98 {
99 largest_scores.insert(it, x);
100 break;
101 }
102 }
103
104 if(largest_scores.size() > 201)
105 {
106 largest_scores.pop_back();
107 }
108 }
109 sum += x;
110 }
111
112 // Data has been added so statistics have now been updated to new values
113 statsAreUpdated = false;
114 ++n;
115 l.unlock();
116 return;
117 }
118
119 void G4ConvergenceTester::calStat()
120 {
121 if (n == 0) return;
122
123 efficiency = G4double(nonzero_histories.size()) / n;
124
125 mean = sum / n;
126
127 G4double sum_x2 = 0.0;
128 var = 0.0;
129 shift = 0.0;
130 vov = 0.0;
131
132 G4double xi;
133 for(const auto& nonzero_historie : nonzero_histories)
134 {
135 xi = nonzero_historie.second;
136 sum_x2 += xi * xi;
137 var += (xi - mean) * (xi - mean);
138 shift += (xi - mean) * (xi - mean) * (xi - mean);
139 vov += (xi - mean) * (xi - mean) * (xi - mean) * (xi - mean);
140 }
141
142 var += (n - nonzero_histories.size()) * mean * mean;
143 shift += (n - nonzero_histories.size()) * mean * mean * mean * (-1);
144 vov += (n - nonzero_histories.size()) * mean * mean * mean * 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 / (efficiency * n);
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_histories)
170 {
171 if(std::abs(nonzero_historie.second) > largest)
172 {
173 largest = nonzero_historie.second;
174 largest_score_happened = nonzero_historie.first;
175 spend_time_of_largest =
176 cpu_time[nonzero_historie.first + 1] - cpu_time[nonzero_historie.first];
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_histories)
191 {
192 xi = nonzero_historie.second;
193 var_1 += (xi - mean_1) * (xi - mean_1);
194 shift_1 += (xi - mean_1) * (xi - mean_1) * (xi - mean_1);
195 vov_1 += (xi - mean_1) * (xi - mean_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) * (xi - mean_1);
200 vov_1 += (xi - mean_1) * (xi - mean_1) * (xi - mean_1) * (xi - mean_1);
201
202 var_1 += (n - nonzero_histories.size()) * mean_1 * mean_1;
203
204 if(var_1 != 0.0)
205 {
206 shift_1 += (n - nonzero_histories.size()) * mean_1 * mean_1 * mean_1 * (-1);
207 vov_1 += (n - nonzero_histories.size()) * mean_1 * mean_1 * mean_1 * mean_1;
208
209 vov_1 = vov_1 / (var_1 * var_1) - 1.0 / (n + 1);
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 + 1));
216
217 shift_1 = shift_1 / (2 * var_1 * (n + 1));
218
219 fom_1 = 1 / (r * r) / (cpu_time.back() + spend_time_of_largest);
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 function does not need
243 // to be called again until data has been added
244 statsAreUpdated = true;
245 }
246
247 void G4ConvergenceTester::calc_grid_point_of_history()
248 {
249 // histroy_grid [ 0,,,15 ]
250 // history_grid [0] 1/16 ,,, history_grid [15] 16/16
251 // if number of event is x then history_grid [15] become x-1.
252 // 16 -> noBinOfHisotry
253
254 for(G4int i = 1; i <= noBinOfHistory; ++i)
255 {
256 history_grid[i - 1] = G4int(n / (G4double(noBinOfHistory)) * i - 0.1);
257 }
258 }
259
260 void G4ConvergenceTester::calc_stat_history()
261 {
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
276 for(const auto& itr : nonzero_histories)
277 {
278 if(itr.first <= ith)
279 {
280 xi = itr.second;
281 mean_till_ith += xi;
282 ++nonzero_till_ith;
283 }
284 }
285
286 if(nonzero_till_ith == 0)
287 {
288 continue;
289 }
290
291 mean_till_ith = mean_till_ith / (ith + 1);
292 mean_history[i] = mean_till_ith;
293
294 G4double sum_x2_till_ith = 0.0;
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 {
303 xi = itr.second;
304 sum_x2_till_ith += std::pow(xi, 2.0);
305 var_till_ith += std::pow(xi - mean_till_ith, 2.0);
306 shift_till_ith += std::pow(xi - mean_till_ith, 3.0);
307 vov_till_ith += std::pow(xi - mean_till_ith, 4.0);
308 }
309 }
310
311 var_till_ith +=
312 ((ith + 1) - nonzero_till_ith) * std::pow(mean_till_ith, 2.0);
313 vov_till_ith +=
314 ((ith + 1) - nonzero_till_ith) * std::pow(mean_till_ith, 4.0);
315
316 G4double sum_till_ith = mean_till_ith * (ith + 1);
317
318 if(!(std::fabs(var_till_ith) > 0.0))
319 {
320 continue;
321 }
322 if(!(std::fabs(mean_till_ith) > 0.0))
323 {
324 continue;
325 }
326 if(!(std::fabs(sum_till_ith) > 0.0))
327 {
328 continue;
329 }
330
331 vov_till_ith = vov_till_ith / std::pow(var_till_ith, 2.0) - 1.0 / (ith + 1);
332 vov_history[i] = vov_till_ith;
333
334 var_till_ith = var_till_ith / (ith + 1 - 1);
335 var_history[i] = var_till_ith;
336 sd_history[i] = std::sqrt(var_till_ith);
337 r_history[i] =
338 std::sqrt(var_till_ith) / mean_till_ith / std::sqrt(1.0 * (ith + 1));
339
340 if(std::fabs(cpu_time[ith]) > 0.0 && std::fabs(r_history[i]) > 0.0)
341 {
342 fom_history[i] = 1.0 / std::pow(r_history[i], 2.0) / cpu_time[ith];
343 }
344 else
345 {
346 fom_history[i] = 0.0;
347 }
348
349 shift_till_ith +=
350 ((ith + 1) - nonzero_till_ith) * std::pow(mean_till_ith, 3.0) * (-1.0);
351 shift_till_ith = shift_till_ith / (2 * var_till_ith * (ith + 1));
352 shift_history[i] = shift_till_ith;
353
354 e_history[i] = 1.0 * nonzero_till_ith / (ith + 1);
355 if(std::fabs(e_history[i]) > 0.0)
356 {
357 r2eff_history[i] = (1 - e_history[i]) / (e_history[i] * (ith + 1));
358
359 r2int_history[i] = (sum_x2_till_ith) / std::pow(sum_till_ith, 2.0) -
360 1 / (e_history[i] * (ith + 1));
361 }
362 }
363 }
364
365 void G4ConvergenceTester::ShowResult(std::ostream& out)
366 {
367 // if data has been added since the last computation of the statistical values
368 // (not statsAreUpdated) call calStat to recompute the statistical values
369 if(!statsAreUpdated)
370 {
371 calStat();
372 }
373
374 out << std::setprecision(6);
375
376 out << G4endl;
377 out << "G4ConvergenceTester Output Result of " << name << G4endl;
378 out << std::setw(20) << "EFFICIENCY = " << std::setw(13) << efficiency
379 << G4endl;
380 out << std::setw(20) << "MEAN = " << std::setw(13) << mean << G4endl;
381 out << std::setw(20) << "VAR = " << std::setw(13) << var << G4endl;
382 out << std::setw(20) << "SD = " << std::setw(13) << sd << G4endl;
383 out << std::setw(20) << "R = " << std::setw(13) << r << G4endl;
384 out << std::setw(20) << "SHIFT = " << std::setw(13) << shift << G4endl;
385 out << std::setw(20) << "VOV = " << std::setw(13) << vov << G4endl;
386 out << std::setw(20) << "FOM = " << std::setw(13) << fom << G4endl;
387
388 out << std::setw(20) << "THE LARGEST SCORE = " << std::setw(13) << largest
389 << " and it happened at " << largest_score_happened << "th event"
390 << G4endl;
391 if(mean != 0)
392 {
393 out << std::setw(20) << "Affected Mean = " << std::setw(13) << mean_1
394 << " and its ratio to original is " << mean_1 / mean << G4endl;
395 }
396 else
397 {
398 out << std::setw(20) << "Affected Mean = " << std::setw(13) << mean_1
399 << G4endl;
400 }
401 if(var != 0)
402 {
403 out << std::setw(20) << "Affected VAR = " << std::setw(13) << var_1
404 << " and its ratio to original is " << var_1 / var << G4endl;
405 }
406 else
407 {
408 out << std::setw(20) << "Affected VAR = " << std::setw(13) << var_1
409 << G4endl;
410 }
411 if(r != 0)
412 {
413 out << std::setw(20) << "Affected R = " << std::setw(13) << r_1
414 << " and its ratio to original is " << r_1 / r << G4endl;
415 }
416 else
417 {
418 out << std::setw(20) << "Affected R = " << std::setw(13) << r_1 << G4endl;
419 }
420 if(shift != 0)
421 {
422 out << std::setw(20) << "Affected SHIFT = " << std::setw(13) << shift_1
423 << " and its ratio to original is " << shift_1 / shift << G4endl;
424 }
425 else
426 {
427 out << std::setw(20) << "Affected SHIFT = " << std::setw(13) << shift_1
428 << G4endl;
429 }
430 if(fom != 0)
431 {
432 out << std::setw(20) << "Affected FOM = " << std::setw(13) << fom_1
433 << " and its ratio to original is " << fom_1 / fom << G4endl;
434 }
435 else
436 {
437 out << std::setw(20) << "Affected FOM = " << std::setw(13) << fom_1
438 << G4endl;
439 }
440
441 if(!showHistory)
442 {
443 out << "Number of events of this run is too small to do convergence tests."
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" << G4endl;
461 }
462 }
463 else
464 {
465 out << "Number of non zero history too small to calculate SLOPE" << G4endl;
466 }
467
468 out << "This result passes " << noPass << " / " << noTotal
469 << " Convergence Test." << G4endl;
470 out << G4endl;
471 }
472
473 void G4ConvergenceTester::ShowHistory(std::ostream& out)
474 {
475 if(!showHistory)
476 {
477 out << "Number of events of this run is too small to show history."
478 << G4endl;
479 return;
480 }
481
482 out << std::setprecision(6);
483
484 out << G4endl;
485 out << "G4ConvergenceTester Output History of " << name << G4endl;
486 out << "i/" << noBinOfHistory << " till_ith mean" << std::setw(13)
487 << "var" << std::setw(13) << "sd" << std::setw(13) << "r" << std::setw(13)
488 << "vov" << std::setw(13) << "fom" << std::setw(13) << "shift"
489 << std::setw(13) << "e" << std::setw(13) << "r2eff" << std::setw(13)
490 << "r2int" << G4endl;
491 for(G4int i = 1; i <= noBinOfHistory; i++)
492 {
493 out << std::setw(4) << i << " " << std::setw(5) << history_grid[i - 1]
494 << std::setw(13) << mean_history[i - 1] << std::setw(13)
495 << var_history[i - 1] << std::setw(13) << sd_history[i - 1]
496 << std::setw(13) << r_history[i - 1] << std::setw(13)
497 << vov_history[i - 1] << std::setw(13) << fom_history[i - 1]
498 << std::setw(13) << shift_history[i - 1] << std::setw(13)
499 << e_history[i - 1] << std::setw(13) << r2eff_history[i - 1]
500 << std::setw(13) << r2int_history[i - 1] << G4endl;
501 }
502 }
503
504 void G4ConvergenceTester::check_stat_history(std::ostream& out)
505 {
506 // 1 sigma rejection for null hypothesis
507
508 std::vector<G4double> first_ally;
509 std::vector<G4double> second_ally;
510
511 // use 2nd half of hisories
512 std::size_t N = mean_history.size() / 2;
513 std::size_t i;
514
515 G4double pearson_r;
516 G4double t;
517
518 first_ally.resize(N);
519 second_ally.resize(N);
520
521 G4double sum_of_var =
522 std::accumulate(var_history.begin(), var_history.end(), 0.0);
523 if(sum_of_var == 0.0)
524 {
525 out << "Variances in all historical grids are zero." << G4endl;
526 out << "Terminating checking behavior of statistics numbers." << G4endl;
527 return;
528 }
529
530 // Mean
531
532 for(i = 0; i < N; ++i)
533 {
534 first_ally[i] = history_grid[N + i];
535 second_ally[i] = mean_history[N + i];
536 }
537
538 pearson_r = calc_Pearson_r((G4int)N, first_ally, second_ally);
539 t = pearson_r * std::sqrt((N - 2) / (1 - pearson_r * pearson_r));
540
541 if(t < 0.429318) // Student t of (Degree of freedom = N-2 )
542 {
543 out << "MEAN distribution is RANDOM" << G4endl;
544 noPass++;
545 }
546 else
547 {
548 out << "MEAN distribution is not RANDOM" << G4endl;
549 }
550
551 // R
552
553 for(i = 0; i < N; ++i)
554 {
555 first_ally[i] = 1.0 / std::sqrt(G4double(history_grid[N + i]));
556 second_ally[i] = r_history[N + i];
557 }
558
559 pearson_r = calc_Pearson_r(G4int(N), first_ally, second_ally);
560 t = pearson_r * std::sqrt((N - 2) / (1 - pearson_r * pearson_r));
561
562 if(t > 1.090546)
563 {
564 out << "r follows 1/std::sqrt(N)" << G4endl;
565 noPass++;
566 }
567 else
568 {
569 out << "r does not follow 1/std::sqrt(N)" << G4endl;
570 }
571
572 if(is_monotonically_decrease(second_ally))
573 {
574 out << "r is monotonically decrease " << G4endl;
575 }
576 else
577 {
578 out << "r is NOT monotonically decrease " << G4endl;
579 }
580
581 if(r_history.back() < 0.1)
582 {
583 out << "r is less than 0.1. r = " << r_history.back() << G4endl;
584 noPass++;
585 }
586 else
587 {
588 out << "r is NOT less than 0.1. r = " << r_history.back() << G4endl;
589 }
590
591 // VOV
592 for(i = 0; i < N; ++i)
593 {
594 first_ally[i] = 1.0 / history_grid[N + i];
595 second_ally[i] = vov_history[N + i];
596 }
597
598 pearson_r = calc_Pearson_r(G4int(N), first_ally, second_ally);
599 t = pearson_r * std::sqrt((N - 2) / (1 - pearson_r * pearson_r));
600
601 if(t > 1.090546)
602 {
603 out << "VOV follows 1/std::sqrt(N)" << G4endl;
604 noPass++;
605 }
606 else
607 {
608 out << "VOV does not follow 1/std::sqrt(N)" << G4endl;
609 }
610
611 if(is_monotonically_decrease(second_ally))
612 {
613 out << "VOV is monotonically decrease " << G4endl;
614 }
615 else
616 {
617 out << "VOV is NOT monotonically decrease " << G4endl;
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::move(first_ally), std::move(second_ally));
629 t = pearson_r * std::sqrt((N - 2) / (1 - pearson_r * pearson_r));
630
631 if(t < 0.429318)
632 {
633 out << "FOM distribution is RANDOM" << G4endl;
634 noPass++;
635 }
636 else
637 {
638 out << "FOM distribution is not RANDOM" << G4endl;
639 }
640 }
641
642 G4double G4ConvergenceTester::calc_Pearson_r(G4int N,
643 std::vector<G4double> first_ally,
644 std::vector<G4double> second_ally)
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) * (second_ally[i] - second_mean);
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) * (first_ally[i] - first_mean);
669 b2 += (second_ally[i] - second_mean) * (second_ally[i] - second_mean);
670 }
671
672 G4double rds = a / std::sqrt(b1 * b2);
673
674 return rds;
675 }
676
677 G4bool G4ConvergenceTester::is_monotonically_decrease(
678 const std::vector<G4double>& ally)
679 {
680 for(auto it = ally.cbegin(); it != ally.cend() - 1; ++it)
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 std::vector<G4double>&)
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 reached
711 slope = 10.0;
712 return;
713 }
714
715 std::vector<G4double> pdf_grid;
716
717 pdf_grid.resize(noBinOfPDF + 1); // no grid = no bins + 1
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 - log10_delta / 10.0 * (i));
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 {
737 pdf[i] += 1.0 / (pdf_grid[i] - pdf_grid[i + 1]) / n;
738 break;
739 }
740 }
741 }
742
743 f_xi.resize(noBinOfPDF);
744 f_yi.resize(noBinOfPDF);
745 for(G4int i = 0; i < noBinOfPDF; ++i)
746 {
747 f_xi[i] = (pdf_grid[i] + pdf_grid[i + 1]) / 2;
748 f_yi[i] = pdf[i];
749 }
750
751 // number of variables ( a and k )
752 minimizer = new G4SimplexDownhill<G4ConvergenceTester>(this, 2);
753 // G4double minimum = minimizer->GetMinimum();
754 std::vector<G4double> mp = minimizer->GetMinimumPoint();
755 G4double k = mp[1];
756
757 // G4cout << "SLOPE " << 1/mp[1]+1 << G4endl;
758 // G4cout << "SLOPE a " << mp[0] << G4endl;
759 // G4cout << "SLOPE k " << mp[1] << G4endl;
760 // G4cout << "SLOPE minimum " << minimizer->GetMinimum() << G4endl;
761
762 slope = 1 / mp[1] + 1;
763 if(k < 1.0 / 9) // Please look Pareto distribution with "sigma=a" and "k"
764 {
765 slope = 10;
766 }
767 if(slope > 10)
768 {
769 slope = 10;
770 }
771 }
772
773 G4double G4ConvergenceTester::slope_fitting_function(std::vector<G4double> x)
774 {
775 G4double a = x[0];
776 G4double k = x[1];
777
778 if(a <= 0)
779 {
780 return 3.402823466e+38; // FLOAT_MAX
781 }
782 if(k == 0)
783 {
784 return 3.402823466e+38; // FLOAT_MAX
785 }
786
787 // f_xi and f_yi is filled at "calc_slope_fit"
788
789 G4double y = 0.0;
790 for(G4int i = 0; i < G4int(f_yi.size()); ++i)
791 {
792 // if ( 1/a * ( 1 + k * f_xi [ i ] / a ) < 0 )
793 if((1 + k * f_xi[i] / a) < 0)
794 {
795 y += 3.402823466e+38; // FLOAT_MAX
796 }
797 else
798 {
799 y += (f_yi[i] - 1 / a * std::pow(1 + k * f_xi[i] / a, -1 / k - 1)) *
800 (f_yi[i] - 1 / a * std::pow(1 + k * f_xi[i] / a, -1 / k - 1));
801 }
802 }
803
804 return y;
805 }
806