source: trunk/examples/extended/analysis/N03Con/src/G4ConvergenceTester.cc @ 1292

Last change on this file since 1292 was 807, checked in by garnier, 16 years ago

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