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

Last change on this file since 1350 was 1350, checked in by garnier, 13 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.