source: trunk/source/global/HEPNumerics/src/G4ConvergenceTester.cc

Last change on this file was 1350, checked in by garnier, 15 years ago

update to last version 4.9.4
File size: 20.0 KB
Line 
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//
27//
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//
42
43#include "G4ConvergenceTester.hh"
44
45G4ConvergenceTester::G4ConvergenceTester()
46 : n(0), sum(0.), mean(0.), var(0.), sd(0.), r(0.), efficiency(0.),
47 r2eff(0.), r2int(0.), shift(0.), vov(0.), fom(0.), largest(0.),
48 largest_score_happened(0), mean_1(0.), var_1(0.), sd_1(0.), r_1(0.),
49 shift_1(0.), vov_1(0.), fom_1(0.), noBinOfHistory(16), slope(0.),
50 noBinOfPDF(10), minimizer(0), noPass(0), noTotal(8)
51{
52 nonzero_histories.clear();
53 largest_scores.clear();
54 largest_scores.push_back( 0.0 );
55
56 history_grid.resize( noBinOfHistory , 0 );
57 mean_history.resize( noBinOfHistory , 0.0 );
58 var_history.resize( noBinOfHistory , 0.0 );
59 sd_history.resize( noBinOfHistory , 0.0 );
60 r_history.resize( noBinOfHistory , 0.0 );
61 vov_history.resize( noBinOfHistory , 0.0 );
62 fom_history.resize( noBinOfHistory , 0.0 );
63 shift_history.resize( noBinOfHistory , 0.0 );
64 e_history.resize( noBinOfHistory , 0.0 );
65 r2eff_history.resize( noBinOfHistory , 0.0 );
66 r2int_history.resize( noBinOfHistory , 0.0 );
67
68 timer = new G4Timer();
69 timer->Start();
70 cpu_time.clear();
71 cpu_time.push_back( 0.0 );
72}
73
74
75
76G4ConvergenceTester::~G4ConvergenceTester()
77{
78 delete timer;
79}
80
81
82
83void G4ConvergenceTester::AddScore( G4double x )
84{
85
86 //G4cout << x << G4endl;
87
88 timer->Stop();
89 cpu_time.push_back( timer->GetSystemElapsed() + timer->GetUserElapsed() );
90
91 if ( x == 0.0 )
92 {
93 }
94 else
95 {
96 nonzero_histories.insert( std::pair< G4int , G4double > ( n , x ) );
97 if ( x > largest_scores.back() )
98 {
99// Following serch should become faster if begin from bottom.
100 std::vector< G4double >::iterator it;
101 for ( it = largest_scores.begin() ; it != largest_scores.end() ; it++ )
102 {
103 if ( x > *it )
104 {
105 largest_scores.insert( it , x );
106 break;
107 }
108 }
109
110 if ( largest_scores.size() > 201 )
111 {
112 largest_scores.pop_back();
113 }
114 //G4cout << largest_scores.size() << " " << largest_scores.front() << " " << largest_scores.back() << G4endl;
115 }
116 sum += x;
117 }
118
119 n++;
120 return;
121}
122
123
124
125void G4ConvergenceTester::calStat()
126{
127
128
129 efficiency = double( nonzero_histories.size() ) / n;
130
131 mean = sum / n;
132
133 G4double sum_x2 = 0.0;
134 var = 0.0;
135 shift = 0.0;
136 vov = 0.0;
137
138 G4double xi;
139 std::map< G4int , G4double >::iterator it;
140 for ( it = nonzero_histories.begin() ; it != nonzero_histories.end() ; it++ )
141 {
142 xi = it->second;
143 sum_x2 += xi * xi;
144 var += ( xi - mean ) * ( xi - mean );
145 shift += ( xi - mean ) * ( xi - mean ) * ( xi - mean );
146 vov += ( xi - mean ) * ( xi - mean ) * ( xi - mean ) * ( xi - mean );
147 }
148
149 var += ( n - nonzero_histories.size() ) * mean * mean;
150 shift += ( n - nonzero_histories.size() ) * mean * mean * mean * ( -1 );
151 vov += ( n - nonzero_histories.size() ) * mean * mean * mean * mean;
152
153 vov = vov / ( var * var ) - 1.0 / n;
154
155 var = var/(n-1);
156
157 sd = std::sqrt ( var );
158
159 r = sd / mean / std::sqrt ( G4double(n) );
160
161 r2eff = ( 1 - efficiency ) / ( efficiency * n );
162 r2int = sum_x2 / ( sum * sum ) - 1 / ( efficiency * n );
163
164 shift = shift / ( 2 * var * n );
165
166 fom = 1 / (r*r) / cpu_time.back();
167
168// Find Largest History
169 //G4double largest = 0.0;
170 largest = 0.0;
171 largest_score_happened = 0;
172 G4double spend_time_of_largest = 0.0;
173 for ( it = nonzero_histories.begin() ; it != nonzero_histories.end() ; it++ )
174 {
175 if ( std::abs ( it->second ) > largest )
176 {
177 largest = it->second;
178 largest_score_happened = it->first;
179 spend_time_of_largest = cpu_time [ it->first+1 ] - cpu_time [ it->first ];
180 }
181 }
182
183 mean_1 = 0.0;
184 var_1 = 0.0;
185 shift_1 = 0.0;
186 vov_1 = 0.0;
187 sd_1 = 0.0;
188 r_1 = 0.0;
189 vov_1 = 0.0;
190
191// G4cout << "The largest history = " << largest << G4endl;
192
193 mean_1 = ( sum + largest ) / ( n + 1 );
194
195 sum_x2 = 0.0;
196 for ( it = nonzero_histories.begin() ; it != nonzero_histories.end() ; it++ )
197 {
198 xi = it->second;
199 sum_x2 += xi * xi;
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 sum_x2 += xi * xi;
206 var_1 += ( xi - mean_1 ) * ( xi - mean_1 );
207 shift_1 += ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 );
208 vov_1 += ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 );
209
210 var_1 += ( n - nonzero_histories.size() ) * mean_1 * mean_1;
211 shift_1 += ( n - nonzero_histories.size() ) * mean_1 * mean_1 * mean_1 * ( -1 );
212 vov_1 += ( n - nonzero_histories.size() ) * mean_1 * mean_1 * mean_1 * mean_1;
213
214 vov_1 = vov_1 / ( var_1 * var_1 ) - 1.0 / ( n + 1 );
215
216 var_1 = var_1 / n ;
217
218 sd_1 = std::sqrt ( var_1 );
219
220 r_1 = sd_1 / mean_1 / std::sqrt ( G4double(n + 1) );
221
222 shift_1 = shift_1 / ( 2 * var_1 * ( n + 1 ) );
223
224 fom_1 = 1 / ( r * r ) / ( cpu_time.back() + spend_time_of_largest );
225
226 if ( nonzero_histories.size() < 500 )
227 {
228 G4cout << "Number of non zero history too small to calcuarte SLOPE" << G4endl;
229 }
230 else
231 {
232 G4int i = int ( nonzero_histories.size() );
233
234 // 5% criterion
235 G4int j = int ( i * 0.05 );
236 while ( int( largest_scores.size() ) > j )
237 {
238 largest_scores.pop_back();
239 }
240 calc_slope_fit( largest_scores );
241 }
242
243 calc_grid_point_of_history();
244 calc_stat_history();
245}
246
247
248
249void G4ConvergenceTester::calc_grid_point_of_history()
250{
251
252// histroy_grid [ 0,,,15 ]
253// history_grid [0] 1/16 ,,, history_grid [15] 16/16
254// if number of event is x then history_grid [15] become x-1.
255// 16 -> noBinOfHisotry
256
257 G4int i;
258 for ( i = 1 ; i <= noBinOfHistory ; i++ )
259 {
260 history_grid [ i-1 ] = int ( n / ( double( noBinOfHistory ) ) * i - 0.1 );
261 //G4cout << "history_grid " << i-1 << " " << history_grid [ i-1 ] << G4endl;
262 }
263
264}
265
266
267
268void G4ConvergenceTester::calc_stat_history()
269{
270// G4cout << "i/16 till_ith mean var sd r vov fom shift e r2eff r2int" << G4endl;
271
272 G4int i;
273 for ( i = 1 ; i <= noBinOfHistory ; i++ )
274 {
275
276 G4int ith = history_grid [ i-1 ];
277
278 G4int nonzero_till_ith = 0;
279 G4double xi;
280 G4double mean_till_ith = 0.0;
281 std::map< G4int , G4double >::iterator it;
282
283 for ( it = nonzero_histories.begin() ; it !=nonzero_histories.end() ; it++ )
284 {
285 if ( it->first <= ith )
286 {
287 xi = it->second;
288 mean_till_ith += xi;
289 nonzero_till_ith++;
290 }
291 }
292
293 mean_till_ith = mean_till_ith / ( ith+1 );
294 mean_history [ i-1 ] = mean_till_ith;
295
296 G4double sum_x2_till_ith = 0.0;
297 G4double var_till_ith = 0.0;
298 G4double vov_till_ith = 0.0;
299 G4double shift_till_ith = 0.0;
300
301 for ( it = nonzero_histories.begin() ; it !=nonzero_histories.end() ; it++ )
302 {
303 if ( it->first <= ith )
304 {
305 xi = it->second;
306 sum_x2_till_ith += xi * xi;
307 var_till_ith += ( xi - mean_till_ith ) * ( xi - mean_till_ith );
308 shift_till_ith += ( xi - mean_till_ith ) * ( xi - mean_till_ith ) * ( xi - mean_till_ith );
309 vov_till_ith += ( xi - mean_till_ith ) * ( xi - mean_till_ith ) * ( xi - mean_till_ith ) * ( xi - mean_till_ith );
310 }
311 }
312
313 var_till_ith += ( (ith+1) - nonzero_till_ith ) * mean_till_ith * mean_till_ith;
314
315 vov_till_ith += ( (ith+1) - nonzero_till_ith ) * mean_till_ith * mean_till_ith * mean_till_ith * mean_till_ith ;
316 vov_till_ith = vov_till_ith / ( var_till_ith * var_till_ith ) - 1.0 / (ith+1);
317 vov_history [ i-1 ] = vov_till_ith;
318
319 var_till_ith = var_till_ith / ( ith+1 - 1 );
320 var_history [ i-1 ] = var_till_ith;
321 sd_history [ i-1 ] = std::sqrt( var_till_ith );
322 r_history [ i-1 ] = std::sqrt( var_till_ith ) / mean_till_ith / std::sqrt ( 1.0*(ith+1) );
323
324 fom_history [ i-1 ] = 1 / ( r_history [ i-1 ] * r_history [ i-1 ] ) / cpu_time [ ith ];
325
326 shift_till_ith += ( (ith+1) - nonzero_till_ith ) * mean_till_ith * mean_till_ith * mean_till_ith * ( -1 );
327 shift_till_ith = shift_till_ith / ( 2 * var_till_ith * (ith+1) );
328 shift_history [ i-1 ] = shift_till_ith;
329
330 e_history [ i-1 ] = 1.0*nonzero_till_ith / (ith+1);
331 r2eff_history [ i-1 ] = ( 1 - e_history [ i-1 ] ) / ( e_history [ i-1 ] * (ith+1) );
332
333 G4double sum_till_ith = mean_till_ith * (ith+1);
334 r2int_history [ i-1 ] = ( sum_x2_till_ith ) / ( sum_till_ith * sum_till_ith ) - 1 / ( e_history [ i-1 ] * (ith+1) );
335
336 }
337
338}
339
340
341
342void G4ConvergenceTester::ShowResult()
343{
344 calStat();
345
346 G4cout << "EFFICIENCY = " << efficiency << G4endl;
347 G4cout << "MEAN = " << mean << G4endl;
348 G4cout << "VAR = " << var << G4endl;
349 G4cout << "SD = " << sd << G4endl;
350 G4cout << "R = "<< r << G4endl;
351 G4cout << "SHIFT = "<< shift << G4endl;
352 G4cout << "VOV = "<< vov << G4endl;
353 G4cout << "FOM = "<< fom << G4endl;
354
355 G4cout << "THE LARGEST SCORE = " << largest << " and it happend at " << largest_score_happened << "th event" << G4endl;
356 G4cout << "Affected Mean = " << mean_1 << " and its ratio to orignal is " << mean_1/mean << G4endl;
357 G4cout << "Affected VAR = " << var_1 << " and its ratio to orignal is " << var_1/var << G4endl;
358 G4cout << "Affected R = " << r_1 << " and its ratio to orignal is " << r_1/r << G4endl;
359 G4cout << "Affected SHIFT = " << shift_1 << " and its ratio to orignal is " << shift_1/shift << G4endl;
360 G4cout << "Affected FOM = " << fom_1 << " and its ratio to orignal is " << fom_1/fom << G4endl;
361
362 check_stat_history();
363
364// check SLOPE and output result
365 if ( slope >= 3 )
366 {
367 noPass++;
368 G4cout << "SLOPE is large enough" << G4endl;
369 }
370 else
371 {
372 G4cout << "SLOPE is not large enough" << G4endl;
373 }
374
375 G4cout << "This result passes " << noPass << " / "<< noTotal << " Convergence Test." << G4endl;
376 G4cout << G4endl;
377
378}
379
380void G4ConvergenceTester::ShowHistory()
381{
382 G4cout << "i/" << noBinOfHistory << " till_ith mean var sd r vov fom shift e r2eff r2int" << G4endl;
383 for ( G4int i = 1 ; i <= noBinOfHistory ; i++ )
384 {
385 G4cout << i << " "
386 << history_grid [ i-1 ] << " "
387 << mean_history [ i-1 ] << " "
388 << var_history [ i-1 ] << " "
389 << sd_history [ i-1 ] << " "
390 << r_history [ i-1 ] << " "
391 << vov_history [ i-1 ] << " "
392 << fom_history [ i-1 ] << " "
393 << shift_history [ i-1 ] << " "
394 << e_history [ i-1 ] << " "
395 << r2eff_history [ i-1 ] << " "
396 << r2int_history [ i-1 ] << " "
397 << G4endl;
398 }
399}
400
401void G4ConvergenceTester::check_stat_history()
402{
403
404// 1 sigma rejection for null hypothesis
405
406 std::vector<G4double> first_ally;
407 std::vector<G4double> second_ally;
408
409// use 2nd half of hisories
410 G4int N = mean_history.size() / 2;
411 G4int i;
412
413 G4double pearson_r;
414 G4double t;
415
416 first_ally.resize( N );
417 second_ally.resize( N );
418
419// Mean
420
421 for ( i = 0 ; i < N ; i++ )
422 {
423 first_ally [ i ] = history_grid [ N + i ];
424 second_ally [ i ] = mean_history [ N + i ];
425 }
426
427 pearson_r = calc_Pearson_r ( N , first_ally , second_ally );
428 t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) );
429
430 if ( t < 0.429318 ) // Student t of (Degree of freedom = N-2 )
431 {
432 G4cout << "MEAN distribution is RANDOM" << G4endl;
433 noPass++;
434 }
435 else
436 {
437 G4cout << "MEAN distribution is not RANDOM" << G4endl;
438 }
439
440
441// R
442
443 for ( i = 0 ; i < N ; i++ )
444 {
445 first_ally [ i ] = 1.0 / std::sqrt ( G4double(history_grid [ N + i ]) );
446 second_ally [ i ] = r_history [ N + i ];
447 }
448
449 pearson_r = calc_Pearson_r ( N , first_ally , second_ally );
450 t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) );
451
452 if ( t > 1.090546 )
453 {
454 G4cout << "r follows 1/std::sqrt(N)" << G4endl;
455 noPass++;
456 }
457 else
458 {
459 G4cout << "r does not follow 1/std::sqrt(N)" << G4endl;
460 }
461
462 if ( is_monotonically_decrease( second_ally ) == true )
463 {
464 G4cout << "r is monotonically decrease " << G4endl;
465 }
466 else
467 {
468 G4cout << "r is NOT monotonically decrease " << G4endl;
469 }
470
471 if ( r_history.back() < 0.1 )
472 {
473 G4cout << "r is less than 0.1. r = " << r_history.back() << G4endl;
474 noPass++;
475 }
476 else
477 {
478 G4cout << "r is NOT less than 0.1. r = " << r_history.back() << G4endl;
479 }
480
481
482// VOV
483 for ( i = 0 ; i < N ; i++ )
484 {
485 first_ally [ i ] = 1.0 / history_grid [ N + i ];
486 second_ally [ i ] = vov_history [ N + i ];
487 }
488
489 pearson_r = calc_Pearson_r ( N , first_ally , second_ally );
490 t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) );
491
492 if ( t > 1.090546 )
493 {
494 G4cout << "VOV follows 1/std::sqrt(N)" << G4endl;
495 noPass++;
496 }
497 else
498 {
499 G4cout << "VOV does not follow 1/std::sqrt(N)" << G4endl;
500 }
501
502 if ( is_monotonically_decrease( second_ally ) == true )
503 {
504 G4cout << "VOV is monotonically decrease " << G4endl;
505 }
506 else
507 {
508 G4cout << "VOV is NOT monotonically decrease " << G4endl;
509 }
510
511// FOM
512
513 for ( i = 0 ; i < N ; i++ )
514 {
515 first_ally [ i ] = history_grid [ N + i ];
516 second_ally [ i ] = fom_history [ N + i ];
517 }
518
519 pearson_r = calc_Pearson_r ( N , first_ally , second_ally );
520 t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) );
521
522 if ( t < 0.429318 )
523 {
524 G4cout << "FOM distribution is RANDOM" << G4endl;
525 noPass++;
526 }
527 else
528 {
529 G4cout << "FOM distribution is not RANDOM" << G4endl;
530 }
531
532}
533
534
535
536G4double G4ConvergenceTester::calc_Pearson_r ( G4int N , std::vector<G4double> first_ally , std::vector<G4double> second_ally )
537{
538 G4double first_mean = 0.0;
539 G4double second_mean = 0.0;
540
541 G4int i;
542 for ( i = 0 ; i < N ; i++ )
543 {
544 first_mean += first_ally [ i ];
545 second_mean += second_ally [ i ];
546 }
547 first_mean = first_mean / N;
548 second_mean = second_mean / N;
549
550 G4double a = 0.0;
551 for ( i = 0 ; i < N ; i++ )
552 {
553 a += ( first_ally [ i ] - first_mean ) * ( second_ally [ i ] - second_mean );
554 }
555
556 G4double b1 = 0.0;
557 G4double b2 = 0.0;
558 for ( i = 0 ; i < N ; i++ )
559 {
560 b1 += ( first_ally [ i ] - first_mean ) * ( first_ally [ i ] - first_mean );
561 b2 += ( second_ally [ i ] - second_mean ) * ( second_ally [ i ] - second_mean );
562 }
563
564 G4double r = a / std::sqrt ( b1 * b2 );
565
566 return r;
567}
568
569
570
571G4bool G4ConvergenceTester::is_monotonically_decrease ( std::vector<G4double> ally )
572{
573
574 std::vector<G4double>::iterator it;
575 for ( it = ally.begin() ; it != ally.end() - 1 ; it++ )
576 {
577 if ( *it < *(it+1) ) return FALSE;
578 }
579
580 noPass++;
581 return TRUE;
582}
583
584
585
586//void G4ConvergenceTester::calc_slope_fit ( std::vector<G4double> largest_socres )
587void G4ConvergenceTester::calc_slope_fit ( std::vector<G4double> )
588{
589
590 // create PDF bins
591 G4double max = largest_scores.front();
592 G4int last = int ( largest_scores.size() );
593 G4double min = 0.0;
594 if ( largest_scores.back() != 0 )
595 {
596 min = largest_scores.back();
597 }
598 else
599 {
600 min = largest_scores[ last-1 ];
601 last = last - 1;
602 }
603
604 //G4cout << "largest " << max << G4endl;
605 //G4cout << "last " << min << G4endl;
606
607 if ( max*0.99 < min )
608 {
609 // upper limit is assumed to have been reached
610 slope = 10.0;
611 return;
612 }
613
614 std::vector < G4double > pdf_grid;
615
616 pdf_grid.resize( noBinOfPDF+1 ); // no grid = no bins + 1
617 pdf_grid[ 0 ] = max;
618 pdf_grid[ noBinOfPDF ] = min;
619 G4double log10_max = std::log10( max );
620 G4double log10_min = std::log10( min );
621 G4double log10_delta = log10_max - log10_min;
622 for ( G4int i = 1 ; i < noBinOfPDF ; i++ )
623 {
624 pdf_grid[i] = std::pow ( 10.0 , log10_max - log10_delta/10.0*(i) );
625 //G4cout << "pdf i " << i << " " << pdf_grid[i] << G4endl;
626 }
627
628 std::vector < G4double > pdf;
629 pdf.resize( noBinOfPDF );
630
631 for ( G4int j=0 ; j < last ; j ++ )
632 {
633 for ( G4int i = 0 ; i < 11 ; i++ )
634 {
635 if ( largest_scores[j] >= pdf_grid[i+1] )
636 {
637 pdf[i] += 1.0 / ( pdf_grid[i] - pdf_grid[i+1] ) / n;
638 //G4cout << "pdf " << j << " " << i << " " << largest_scores[j] << " " << G4endl;
639 break;
640 }
641 }
642 }
643
644 f_xi.resize( noBinOfPDF );
645 f_yi.resize( noBinOfPDF );
646 for ( G4int i = 0 ; i < noBinOfPDF ; i++ )
647 {
648 //G4cout << "pdf i " << i << " " << (pdf_grid[i]+pdf_grid[i+1])/2 << " " << pdf[i] << G4endl;
649 f_xi[i] = (pdf_grid[i]+pdf_grid[i+1])/2;
650 f_yi[i] = pdf[i];
651 }
652
653 // number of variables ( a and k )
654 minimizer = new G4SimplexDownhill<G4ConvergenceTester> ( this , 2 );
655 //G4double minimum = minimizer->GetMinimum();
656 std::vector<G4double> mp = minimizer->GetMinimumPoint();
657 G4double k = mp[1];
658
659 //G4cout << "SLOPE " << 1/mp[1]+1 << G4endl;
660 //G4cout << "SLOPE a " << mp[0] << G4endl;
661 //G4cout << "SLOPE k " << mp[1] << G4endl;
662 //G4cout << "SLOPE minimum " << minimizer->GetMinimum() << G4endl;
663
664 slope = 1/mp[1]+1;
665 if ( k < 1.0/9 ) // Please look Pareto distribution with "sigma=a" and "k"
666 {
667 slope = 10;
668 }
669 if ( slope > 10 )
670 {
671 slope = 10;
672 }
673}
674
675
676
677G4double G4ConvergenceTester::slope_fitting_function ( std::vector< G4double > x )
678{
679
680 G4double a = x[0];
681 G4double k = x[1];
682
683 if ( a <= 0 )
684 {
685 return 3.402823466e+38; // FLOAT_MAX
686 }
687 if ( k == 0 )
688 {
689 return 3.402823466e+38; // FLOAT_MAX
690 }
691
692// f_xi and f_yi is filled at "calc_slope_fit"
693
694 G4double y = 0.0;
695 G4int i;
696 for ( i = 0 ; i < int ( f_yi.size() ) ; i++ )
697 {
698 //if ( 1/a * ( 1 + k * f_xi [ i ] / a ) < 0 )
699 if ( ( 1 + k * f_xi [ i ] / a ) < 0 )
700 {
701 y +=3.402823466e+38; // FLOAT_MAX
702 }
703 else
704 {
705 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 ) );
706 }
707 }
708// G4cout << "y = " << y << G4endl;
709
710 return y;
711}
Note: See TracBrowser for help on using the repository browser.