G4ConvergenceTester.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 //
00028 // Convergence Tests for Monte Carlo results.
00029 //
00030 // Reference
00031 // MCNP(TM) -A General Monte Carlo N-Particle Transport Code
00032 // Version 4B
00033 // Judith F. Briesmeister, Editor
00034 // LA-12625-M, Issued: March 1997, UC 705 and UC 700
00035 // CHAPTER 2. GEOMETRY, DATA, PHYSICS, AND MATHEMATICS
00036 //        VI. ESTIMATION OF THE MONTE CARLO PRECISION
00037 //
00038 // Positives numbers are assumed for inputs
00039 //
00040 // Koi, Tatsumi (SLAC/SCCS)
00041 //
00042 
00043 #include "G4ConvergenceTester.hh"
00044 
00045 G4ConvergenceTester::G4ConvergenceTester()
00046  : n(0), sum(0.), mean(0.), var(0.), sd(0.), r(0.), efficiency(0.),
00047    r2eff(0.), r2int(0.), shift(0.), vov(0.), fom(0.), largest(0.),
00048    largest_score_happened(0), mean_1(0.), var_1(0.), sd_1(0.), r_1(0.),
00049    shift_1(0.), vov_1(0.), fom_1(0.), noBinOfHistory(16), slope(0.),
00050    noBinOfPDF(10), minimizer(0), noPass(0), noTotal(8), statsAreUpdated(true)
00051 {
00052    nonzero_histories.clear();
00053    largest_scores.clear();
00054    largest_scores.push_back( 0.0 );
00055 
00056    history_grid.resize( noBinOfHistory , 0 );
00057    mean_history.resize( noBinOfHistory , 0.0 );
00058    var_history.resize( noBinOfHistory , 0.0 );
00059    sd_history.resize( noBinOfHistory , 0.0 );
00060    r_history.resize( noBinOfHistory , 0.0 );
00061    vov_history.resize( noBinOfHistory , 0.0 );
00062    fom_history.resize( noBinOfHistory , 0.0 );
00063    shift_history.resize( noBinOfHistory , 0.0 );
00064    e_history.resize( noBinOfHistory , 0.0 );
00065    r2eff_history.resize( noBinOfHistory , 0.0 );
00066    r2int_history.resize( noBinOfHistory , 0.0 );
00067 
00068    timer = new G4Timer();
00069    timer->Start();
00070    cpu_time.clear();
00071    cpu_time.push_back( 0.0 );
00072 }
00073 
00074 
00075 
00076 G4ConvergenceTester::~G4ConvergenceTester()
00077 {
00078    delete timer;
00079 }
00080 
00081 
00082 
00083 void G4ConvergenceTester::AddScore( G4double x )
00084 { 
00085 
00086    //G4cout << x << G4endl;
00087 
00088    timer->Stop();
00089    cpu_time.push_back( timer->GetSystemElapsed() + timer->GetUserElapsed() );
00090 
00091    if ( x == 0.0 ) 
00092    { 
00093    }
00094    else
00095    {
00096        nonzero_histories.insert( std::pair< G4int , G4double > ( n , x ) );
00097        if ( x > largest_scores.back() )  
00098        { 
00099 //        Following serch  should become faster if begin from bottom.
00100           std::vector< G4double >::iterator it; 
00101           for ( it = largest_scores.begin() ; it != largest_scores.end() ; it++ )
00102           { 
00103              if ( x > *it ) 
00104              { 
00105                 largest_scores.insert( it , x );
00106                 break;
00107              }
00108           }
00109 
00110           if ( largest_scores.size() > 201 )
00111           {
00112              largest_scores.pop_back();
00113           } 
00114           //G4cout << largest_scores.size() << " " << largest_scores.front() << " " << largest_scores.back() << G4endl;
00115        }
00116        sum += x; 
00117    }
00118 
00119     // Data has been added so statistics have not been updated to new values
00120    statsAreUpdated = false;
00121    n++;
00122    return; 
00123 }
00124 
00125 
00126 
00127 void G4ConvergenceTester::calStat()
00128 {
00129     
00130 
00131    efficiency = double( nonzero_histories.size() ) / n; 
00132 
00133    mean = sum / n;
00134    
00135    G4double sum_x2 = 0.0; 
00136    var = 0.0;
00137    shift = 0.0;
00138    vov = 0.0;
00139 
00140    G4double xi; 
00141    std::map< G4int , G4double >::iterator it;
00142    for ( it = nonzero_histories.begin() ; it != nonzero_histories.end() ; it++ )
00143    {
00144       xi = it->second;
00145       sum_x2 += xi * xi;
00146       var += ( xi - mean ) * ( xi - mean );
00147       shift += ( xi - mean ) * ( xi - mean ) * ( xi - mean );
00148       vov += ( xi - mean ) * ( xi - mean ) * ( xi - mean ) * ( xi - mean );
00149    }
00150 
00151    var += ( n - nonzero_histories.size() ) * mean * mean;
00152    shift += ( n - nonzero_histories.size() ) * mean * mean * mean * ( -1 );
00153    vov += ( n - nonzero_histories.size() ) * mean * mean * mean * mean;
00154 
00155    vov = vov / ( var * var ) - 1.0 / n; 
00156 
00157    var = var/(n-1);
00158 
00159    sd = std::sqrt ( var );
00160 
00161    r = sd / mean / std::sqrt ( G4double(n) ); 
00162 
00163    r2eff = ( 1 - efficiency ) / ( efficiency * n );
00164    r2int = sum_x2 / ( sum * sum ) - 1 / ( efficiency * n );
00165    
00166    shift = shift / ( 2 * var * n );
00167    
00168    fom =  1 / (r*r) / cpu_time.back(); 
00169 
00170 // Find Largest History 
00171    //G4double largest = 0.0;
00172    largest = 0.0;
00173    largest_score_happened = 0;
00174    G4double spend_time_of_largest = 0.0;
00175    for ( it = nonzero_histories.begin() ; it != nonzero_histories.end() ; it++ )
00176    {
00177       if ( std::abs ( it->second ) > largest )
00178       {
00179          largest = it->second;
00180          largest_score_happened = it->first;
00181          spend_time_of_largest = cpu_time [ it->first+1 ] - cpu_time [ it->first ];
00182       }
00183    }
00184 
00185    mean_1 = 0.0;
00186    var_1 = 0.0;
00187    shift_1 = 0.0;
00188    vov_1 = 0.0;
00189    sd_1 = 0.0;
00190    r_1 = 0.0;
00191    vov_1 = 0.0;
00192 
00193 //   G4cout << "The largest history  = " << largest << G4endl;
00194 
00195    mean_1 = ( sum + largest ) / ( n + 1 );
00196 
00197    for ( it = nonzero_histories.begin() ; it != nonzero_histories.end() ; it++ )
00198    {
00199       xi = it->second;
00200       var_1 += ( xi - mean_1 ) * ( xi - mean_1 );
00201       shift_1 += ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 );
00202       vov_1 += ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 );
00203    }
00204    xi = largest;
00205    var_1 += ( xi - mean_1 ) * ( xi - mean_1 );
00206    shift_1 += ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 );
00207    vov_1 += ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 );
00208 
00209    var_1 += ( n - nonzero_histories.size() ) * mean_1 * mean_1;
00210    shift_1 += ( n - nonzero_histories.size() ) * mean_1 * mean_1 * mean_1 * ( -1 );
00211    vov_1 += ( n - nonzero_histories.size() ) * mean_1 * mean_1 * mean_1 * mean_1;
00212 
00213    vov_1 = vov_1 / ( var_1 * var_1 ) - 1.0 / ( n + 1 ); 
00214 
00215    var_1 = var_1 / n ;
00216 
00217    sd_1 = std::sqrt ( var_1 );
00218 
00219    r_1 = sd_1 / mean_1 / std::sqrt ( G4double(n + 1) ); 
00220 
00221    shift_1 = shift_1 / ( 2 * var_1 * ( n + 1 ) );
00222 
00223    fom_1 = 1 / ( r * r ) / ( cpu_time.back() + spend_time_of_largest );
00224 
00225    if ( nonzero_histories.size() < 500 ) 
00226    {
00227       G4cout << "Number of non zero history too small to calculate SLOPE" << G4endl;
00228    }
00229    else
00230    {
00231       G4int i = int ( nonzero_histories.size() );
00232 
00233       // 5% criterion
00234       G4int j = int ( i * 0.05 );
00235       while ( int( largest_scores.size() ) > j ) 
00236       {
00237          largest_scores.pop_back();
00238       }
00239       calc_slope_fit( largest_scores );   
00240    }
00241 
00242    calc_grid_point_of_history();
00243    calc_stat_history();
00244    
00245    // statistics have been calculated and this function does not need
00246    // to be called again until data has been added
00247    statsAreUpdated = true;
00248 }
00249 
00250 
00251 
00252 void G4ConvergenceTester::calc_grid_point_of_history()
00253 {
00254 
00255 // histroy_grid [ 0,,,15 ] 
00256 // history_grid [0] 1/16 ,,, history_grid [15] 16/16
00257 // if number of event is x then history_grid [15] become x-1. 
00258 // 16 -> noBinOfHisotry
00259 
00260    G4int i;
00261    for ( i = 1 ; i <= noBinOfHistory ; i++ )
00262    { 
00263       history_grid [ i-1 ] = int ( n / ( double( noBinOfHistory ) ) * i - 0.1 );
00264       //G4cout << "history_grid " << i-1  << " " << history_grid [ i-1 ] << G4endl;
00265    }
00266 
00267 }
00268 
00269 
00270 
00271 void G4ConvergenceTester::calc_stat_history()
00272 {
00273 //   G4cout << "i/16  till_ith  mean  var  sd  r  vov  fom  shift  e  r2eff  r2int" << G4endl;
00274 
00275    G4int i;
00276    for ( i = 1 ; i <=  noBinOfHistory  ; i++ )
00277    {
00278 
00279       G4int ith = history_grid [ i-1 ];
00280 
00281       G4int nonzero_till_ith = 0;
00282       G4double xi;
00283       G4double mean_till_ith = 0.0;  
00284       std::map< G4int , G4double >::iterator it;
00285 
00286       for ( it = nonzero_histories.begin() ; it !=nonzero_histories.end() ; it++ )
00287       {
00288          if ( it->first <= ith )
00289          {
00290             xi = it->second;
00291             mean_till_ith += xi;
00292             nonzero_till_ith++; 
00293          }
00294       }
00295 
00296       mean_till_ith = mean_till_ith / ( ith+1 ); 
00297       mean_history [ i-1 ] = mean_till_ith;
00298 
00299       G4double sum_x2_till_ith = 0.0;
00300       G4double var_till_ith = 0.0;
00301       G4double vov_till_ith = 0.0;
00302       G4double shift_till_ith = 0.0;
00303   
00304       for ( it = nonzero_histories.begin() ; it !=nonzero_histories.end() ; it++ )
00305       {
00306          if ( it->first <= ith )
00307          {
00308          xi = it->second;
00309          sum_x2_till_ith += xi * xi; 
00310          var_till_ith += ( xi - mean_till_ith ) * ( xi - mean_till_ith );
00311          shift_till_ith += ( xi - mean_till_ith ) * ( xi - mean_till_ith ) * ( xi - mean_till_ith );
00312          vov_till_ith += ( xi - mean_till_ith ) * ( xi - mean_till_ith ) * ( xi - mean_till_ith ) * ( xi - mean_till_ith );
00313          }
00314       }
00315 
00316       var_till_ith += ( (ith+1) - nonzero_till_ith ) * mean_till_ith * mean_till_ith;
00317 
00318       vov_till_ith += ( (ith+1) - nonzero_till_ith ) * mean_till_ith * mean_till_ith * mean_till_ith * mean_till_ith ;
00319       vov_till_ith = vov_till_ith / ( var_till_ith * var_till_ith ) - 1.0 / (ith+1); 
00320       vov_history [ i-1 ] = vov_till_ith;
00321 
00322       var_till_ith = var_till_ith / ( ith+1 - 1 );
00323       var_history [ i-1 ] = var_till_ith;
00324       sd_history [ i-1 ] = std::sqrt( var_till_ith );
00325       r_history [ i-1 ] = std::sqrt( var_till_ith ) / mean_till_ith / std::sqrt ( 1.0*(ith+1) );
00326 
00327       fom_history [ i-1 ] = 1 / ( r_history [ i-1 ] *  r_history [ i-1 ] ) / cpu_time [ ith ];
00328    
00329       shift_till_ith += ( (ith+1) - nonzero_till_ith ) * mean_till_ith * mean_till_ith * mean_till_ith * ( -1 );
00330       shift_till_ith = shift_till_ith / ( 2 * var_till_ith * (ith+1) );
00331       shift_history [ i-1 ] = shift_till_ith;
00332 
00333       e_history [ i-1 ] = 1.0*nonzero_till_ith / (ith+1);
00334       r2eff_history [ i-1 ] = ( 1 - e_history [ i-1 ] ) / (  e_history [ i-1 ] * (ith+1) );
00335 
00336       G4double sum_till_ith =  mean_till_ith * (ith+1); 
00337       r2int_history [ i-1 ] = ( sum_x2_till_ith ) / ( sum_till_ith * sum_till_ith ) - 1 / ( e_history [ i-1 ] * (ith+1) );
00338 
00339    }
00340 
00341 }
00342 
00343 
00344 
00345 void G4ConvergenceTester::ShowResult(std::ostream& out)
00346 {
00347     // if data has been added since the last computation of the statistical values (not statsAreUpdated)
00348     // call calStat to recompute the statistical values
00349    if(!statsAreUpdated) { calStat(); }
00350 
00351    out << "EFFICIENCY = " << efficiency << G4endl;
00352    out << "MEAN = " << mean << G4endl;
00353    out << "VAR = " << var << G4endl;
00354    out << "SD = " << sd << G4endl;
00355    out << "R = "<< r << G4endl;
00356    out << "SHIFT = "<< shift << G4endl;
00357    out << "VOV = "<< vov << G4endl;
00358    out << "FOM = "<< fom << G4endl;
00359 
00360    out << "THE LARGEST SCORE = " << largest << " and it happend at " << largest_score_happened << "th event" << G4endl;
00361    out << "Affected Mean = " << mean_1 << " and its ratio to orignal is " << mean_1/mean << G4endl;
00362    out << "Affected VAR = " << var_1 << " and its ratio to orignal is " << var_1/var << G4endl;
00363    out << "Affected R = " << r_1 << " and its ratio to orignal is " << r_1/r << G4endl;
00364    out << "Affected SHIFT = " << shift_1 << " and its ratio to orignal is " << shift_1/shift << G4endl;
00365    out << "Affected FOM = " << fom_1 << " and its ratio to orignal is " << fom_1/fom << G4endl;
00366 
00367    check_stat_history(out);
00368 
00369 // check SLOPE and output result
00370    if ( slope >= 3 )
00371    {    
00372       noPass++;
00373       out << "SLOPE is large enough" << G4endl; 
00374    }
00375    else
00376    {
00377       out << "SLOPE is not large enough" << G4endl; 
00378    }
00379 
00380    out << "This result passes " << noPass << " / "<< noTotal << " Convergence Test." << G4endl; 
00381    out << G4endl;
00382 
00383 }
00384 
00385 void G4ConvergenceTester::ShowHistory(std::ostream& out)
00386 {
00387    out << "i/" << noBinOfHistory << " till_ith  mean  var  sd  r  vov  fom  shift  e  r2eff  r2int" << G4endl;
00388    for ( G4int i = 1 ; i <=  noBinOfHistory  ; i++ )
00389    {
00390       out << i << " " 
00391              << history_grid [ i-1 ] << " "
00392              << mean_history [ i-1 ] << " "  
00393              << var_history [ i-1 ] << " "  
00394              << sd_history [ i-1 ] << " "  
00395              << r_history [ i-1 ] << " "  
00396              << vov_history [ i-1 ] << " "  
00397              << fom_history [ i-1 ] << " " 
00398              << shift_history [ i-1 ] << " "  
00399              << e_history [ i-1 ] << " " 
00400              << r2eff_history [ i-1 ] << " " 
00401              << r2int_history [ i-1 ] << " " 
00402              << G4endl;
00403    }
00404 }
00405 
00406 void G4ConvergenceTester::check_stat_history(std::ostream& out)
00407 {
00408 
00409 // 1 sigma rejection for null hypothesis 
00410 
00411    std::vector<G4double> first_ally;
00412    std::vector<G4double> second_ally;
00413 
00414 // use 2nd half of hisories
00415    G4int N = mean_history.size() / 2;
00416    G4int i;
00417  
00418    G4double pearson_r;
00419    G4double t;
00420 
00421    first_ally.resize( N );
00422    second_ally.resize( N );
00423 
00424 // Mean
00425 
00426    for ( i = 0 ; i < N ; i++ ) 
00427    {
00428       first_ally [ i ] = history_grid [ N + i ];
00429       second_ally [ i ] = mean_history [ N + i ];
00430    }
00431 
00432    pearson_r = calc_Pearson_r ( N , first_ally , second_ally );
00433    t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) );
00434 
00435    if ( t < 0.429318 ) // Student t of (Degree of freedom = N-2 )  
00436    {    
00437       out << "MEAN distribution is  RANDOM" << G4endl; 
00438       noPass++;
00439    }
00440    else
00441    {
00442       out << "MEAN distribution is not RANDOM" << G4endl; 
00443    }
00444 
00445 
00446 // R
00447 
00448    for ( i = 0 ; i < N ; i++ ) 
00449    {
00450       first_ally [ i ] = 1.0 / std::sqrt ( G4double(history_grid [ N + i ]) );
00451       second_ally [ i ] = r_history [ N + i ];
00452    }
00453 
00454    pearson_r = calc_Pearson_r ( N , first_ally , second_ally );
00455    t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) );
00456       
00457    if ( t > 1.090546 )
00458    {    
00459       out << "r follows 1/std::sqrt(N)" << G4endl; 
00460       noPass++;
00461    }
00462    else
00463    {
00464       out << "r does not follow 1/std::sqrt(N)" << G4endl; 
00465    }
00466 
00467    if (  is_monotonically_decrease( second_ally ) == true ) 
00468    {
00469       out << "r is monotonically decrease " << G4endl;
00470    }
00471    else
00472    {
00473       out << "r is NOT monotonically decrease " << G4endl;
00474    }
00475 
00476    if ( r_history.back() < 0.1 )  
00477    {
00478       out << "r is less than 0.1. r = " <<  r_history.back() << G4endl;
00479       noPass++;
00480    }
00481    else
00482    {
00483       out << "r is NOT less than 0.1. r = " <<  r_history.back() << G4endl;
00484    }
00485 
00486 
00487 // VOV
00488    for ( i = 0 ; i < N ; i++ ) 
00489    {
00490       first_ally [ i ] = 1.0 / history_grid [ N + i ];
00491       second_ally [ i ] = vov_history [ N + i ];
00492    }
00493 
00494    pearson_r = calc_Pearson_r ( N , first_ally , second_ally );
00495    t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) );
00496       
00497    if ( t > 1.090546 )
00498    {    
00499       out << "VOV follows 1/std::sqrt(N)" << G4endl; 
00500       noPass++;
00501    }
00502    else
00503    {
00504       out << "VOV does not follow 1/std::sqrt(N)" << G4endl; 
00505    }
00506 
00507    if ( is_monotonically_decrease( second_ally ) == true )
00508    {
00509       out << "VOV is monotonically decrease " << G4endl;
00510    }
00511    else
00512    {
00513       out << "VOV is NOT monotonically decrease " << G4endl;
00514    }
00515 
00516 // FOM
00517 
00518    for ( i = 0 ; i < N ; i++ ) 
00519    {
00520       first_ally [ i ] = history_grid [ N + i ];
00521       second_ally [ i ] = fom_history [ N + i ];
00522    }
00523 
00524    pearson_r = calc_Pearson_r ( N , first_ally , second_ally );
00525    t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) );
00526 
00527    if ( t < 0.429318 )
00528    {    
00529       out << "FOM distribution is RANDOM" << G4endl; 
00530       noPass++;
00531    }
00532    else
00533    {
00534       out << "FOM distribution is not RANDOM" << G4endl; 
00535    }
00536 
00537 }
00538 
00539 
00540 
00541 G4double G4ConvergenceTester::calc_Pearson_r ( G4int N , std::vector<G4double> first_ally , std::vector<G4double> second_ally )
00542 {
00543    G4double first_mean = 0.0;  
00544    G4double second_mean = 0.0;
00545 
00546    G4int i;
00547    for ( i = 0 ; i < N ; i++ )
00548    {
00549       first_mean += first_ally [ i ]; 
00550       second_mean += second_ally [ i ]; 
00551    }
00552    first_mean = first_mean / N;
00553    second_mean = second_mean / N;
00554    
00555    G4double a = 0.0; 
00556    for ( i = 0 ; i < N ; i++ )
00557    {
00558       a += ( first_ally [ i ] - first_mean ) * ( second_ally [ i ] - second_mean );
00559    }
00560 
00561    G4double b1 = 0.0; 
00562    G4double b2 = 0.0;
00563    for ( i = 0 ; i < N ; i++ )
00564    {
00565       b1 += ( first_ally [ i ] - first_mean ) * ( first_ally [ i ] - first_mean );
00566       b2 += ( second_ally [ i ] - second_mean ) * ( second_ally [ i ] - second_mean );
00567    }
00568    
00569    G4double rds = a / std::sqrt ( b1 * b2 );  
00570 
00571    return rds; 
00572 }
00573 
00574 
00575 
00576 G4bool G4ConvergenceTester::is_monotonically_decrease ( std::vector<G4double> ally )
00577 {
00578 
00579    std::vector<G4double>::iterator it;
00580    for ( it = ally.begin() ; it != ally.end() - 1 ; it++ )
00581    {
00582       if ( *it < *(it+1) ) return FALSE;
00583    }
00584 
00585    noPass++;
00586    return TRUE;
00587 }
00588 
00589 
00590 
00591 //void G4ConvergenceTester::calc_slope_fit ( std::vector<G4double> largest_socres )
00592 void G4ConvergenceTester::calc_slope_fit ( std::vector<G4double> )
00593 {
00594 
00595    // create PDF bins 
00596    G4double max = largest_scores.front();
00597    G4int last = int ( largest_scores.size() );
00598    G4double min = 0.0;
00599    if (  largest_scores.back() !=  0 ) 
00600    {
00601       min = largest_scores.back();
00602    }
00603    else
00604    {
00605       min = largest_scores[ last-1 ];
00606       last = last - 1;
00607    }
00608     
00609    //G4cout << "largest " << max << G4endl;
00610    //G4cout << "last  " << min << G4endl;
00611 
00612    if ( max*0.99 < min )  
00613    {
00614       // upper limit is assumed to have been reached
00615       slope = 10.0;
00616       return;
00617    }
00618 
00619    std::vector < G4double >  pdf_grid;
00620 
00621    pdf_grid.resize( noBinOfPDF+1 );   // no grid  = no bins + 1
00622    pdf_grid[ 0 ] = max; 
00623    pdf_grid[ noBinOfPDF ] = min; 
00624    G4double log10_max = std::log10( max );
00625    G4double log10_min = std::log10( min );
00626    G4double log10_delta = log10_max - log10_min;
00627    for ( G4int i = 1 ; i < noBinOfPDF ; i++ )
00628    {
00629       pdf_grid[i] = std::pow ( 10.0 , log10_max - log10_delta/10.0*(i) );    
00630       //G4cout << "pdf i " << i << " " << pdf_grid[i] << G4endl;
00631    }
00632    
00633    std::vector < G4double >  pdf;
00634    pdf.resize( noBinOfPDF ); 
00635 
00636    for ( G4int j=0 ; j < last ; j ++ )
00637    {
00638       for ( G4int i = 0 ; i < 11 ; i++ )
00639       {
00640          if ( largest_scores[j] >= pdf_grid[i+1] )  
00641          {
00642             pdf[i] += 1.0 / ( pdf_grid[i] - pdf_grid[i+1] ) / n;
00643             //G4cout << "pdf " << j << " " << i << " " <<  largest_scores[j]  << " " << G4endl;
00644             break;
00645          }
00646       }
00647    }
00648 
00649    f_xi.resize( noBinOfPDF );
00650    f_yi.resize( noBinOfPDF );
00651    for ( G4int i = 0 ; i < noBinOfPDF ; i++ )
00652    {
00653       //G4cout << "pdf i " << i << " " <<  (pdf_grid[i]+pdf_grid[i+1])/2 << " " << pdf[i] << G4endl;
00654       f_xi[i] = (pdf_grid[i]+pdf_grid[i+1])/2;
00655       f_yi[i] = pdf[i];
00656    }
00657 
00658    //                                                  number of variables ( a and k )  
00659    minimizer = new G4SimplexDownhill<G4ConvergenceTester> ( this , 2 ); 
00660    //G4double minimum =  minimizer->GetMinimum();
00661    std::vector<G4double> mp = minimizer->GetMinimumPoint();
00662    G4double k = mp[1];
00663 
00664    //G4cout << "SLOPE " << 1/mp[1]+1 << G4endl;
00665    //G4cout << "SLOPE  a " << mp[0] << G4endl;
00666    //G4cout << "SLOPE  k " << mp[1] << G4endl;
00667    //G4cout << "SLOPE  minimum " << minimizer->GetMinimum() << G4endl;
00668 
00669    slope = 1/mp[1]+1;
00670    if ( k < 1.0/9 )  // Please look Pareto distribution with "sigma=a" and "k" 
00671    {
00672       slope = 10;
00673    } 
00674    if ( slope > 10 ) 
00675    {
00676       slope = 10; 
00677    }
00678 }
00679 
00680 
00681 
00682 G4double G4ConvergenceTester::slope_fitting_function ( std::vector< G4double > x )
00683 {
00684 
00685    G4double a = x[0];
00686    G4double k = x[1];
00687 
00688    if ( a <= 0 ) 
00689    {
00690       return 3.402823466e+38;  // FLOAT_MAX
00691    } 
00692    if ( k == 0 ) 
00693    {
00694       return 3.402823466e+38;  // FLOAT_MAX
00695    }
00696 
00697 // f_xi and f_yi is filled at "calc_slope_fit"
00698 
00699    G4double y = 0.0;
00700    G4int i;
00701    for ( i = 0 ; i < int ( f_yi.size() ) ; i++ )
00702    {
00703       //if ( 1/a * ( 1 + k * f_xi [ i ] / a ) < 0 )
00704       if ( ( 1 + k * f_xi [ i ] / a ) < 0 )
00705       {
00706          y +=3.402823466e+38;  // FLOAT_MAX 
00707       }
00708       else 
00709       {
00710          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 ) );
00711       }
00712    }
00713 //   G4cout << "y = " << y << G4endl;
00714 
00715    return y;
00716 }

Generated on Mon May 27 17:47:57 2013 for Geant4 by  doxygen 1.4.7