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