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