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

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