| Geant4 Cross Reference |
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 // G4ConvergenceTester
27 //
28 // Class description:
29 //
30 // Convergence Tests for Monte Carlo results.
31 //
32 // Reference
33 // MCNP(TM) -A General Monte Carlo N-Particle Transport Code
34 // Version 4B
35 // Judith F. Briesmeister, Editor
36 // LA-12625-M, Issued: March 1997, UC 705 and UC 700
37 // CHAPTER 2. GEOMETRY, DATA, PHYSICS, AND MATHEMATICS
38 // VI. ESTIMATION OF THE MONTE CARLO PRECISION
39 //
40 // Positive numbers are assumed for input values
41
42 // Author: Tatsumi Koi (SLAC/SCCS)
43 // --------------------------------------------------------------------
44 #ifndef G4ConvergenceTester_hh
45 #define G4ConvergenceTester_hh 1
46
47 #include "G4SimplexDownhill.hh"
48 #include "G4Timer.hh"
49 #include "globals.hh"
50
51 #include <map>
52 #include <vector>
53
54 class G4ConvergenceTester
55 {
56 public:
57
58 G4ConvergenceTester(const G4String& theName = "NONAME");
59 ~G4ConvergenceTester();
60 G4ConvergenceTester(G4double);
61
62 void AddScore(G4double);
63
64 inline G4ConvergenceTester& operator+=(G4double val)
65 {
66 this->AddScore(val);
67 return *this;
68 }
69
70 void ShowHistory(std::ostream& out = G4cout);
71 void ShowResult(std::ostream& out = G4cout);
72 // Default to G4cout but can be redirected to another ostream
73
74 inline G4double GetValueOfMinimizingFunction(std::vector<G4double> x)
75 {
76 return slope_fitting_function(std::move(x));
77 }
78
79 void ComputeStatistics() { calStat(); }
80 // Explicitly calculate statistics
81
82 // All accessors check to make sure value is current before returning
83
84 inline G4double GetMean() { CheckIsUpdated(); return mean; }
85 inline G4double GetStandardDeviation() { CheckIsUpdated(); return sd; }
86 inline G4double GetVariance() { CheckIsUpdated(); return var; }
87 inline G4double GetR() { CheckIsUpdated(); return r; }
88 inline G4double GetEfficiency() { CheckIsUpdated(); return efficiency; }
89 inline G4double GetR2eff() { CheckIsUpdated(); return r2eff; }
90 inline G4double GetR2int() { CheckIsUpdated(); return r2int; }
91 inline G4double GetShift() { CheckIsUpdated(); return shift; }
92 inline G4double GetVOV() { CheckIsUpdated(); return vov; }
93 inline G4double GetFOM() { CheckIsUpdated(); return fom; }
94
95 private:
96
97 void calStat();
98 // Boolean value of 'statsAreUpdated' is set to TRUE at end of calStat
99 // and set to FALSE at end of AddScore().
100 // NOTE: A thread lock needs to be put in AddScore() so calStat()
101 // is not executed in one thread while AddScore() is adding data
102
103 inline void CheckIsUpdated()
104 {
105 if(!statsAreUpdated) { calStat(); }
106 }
107
108 void calc_grid_point_of_history();
109 void calc_stat_history();
110 void check_stat_history(std::ostream& out = G4cout);
111 G4double calc_Pearson_r(G4int,std::vector<G4double>,std::vector<G4double>);
112 G4bool is_monotonically_decrease(const std::vector<G4double>&);
113 void calc_slope_fit(const std::vector<G4double>&);
114 G4double slope_fitting_function(std::vector<G4double>);
115
116 private:
117
118 G4String name;
119 std::map<G4int, G4double> nonzero_histories; // (ith-history, score value)
120 G4int n = 0; // number of history
121 G4double sum = 0.0; // sum of scores;
122
123 G4Timer* timer = nullptr;
124 std::vector<G4double> cpu_time;
125
126 G4double mean = 0.0;
127 G4double var = 0.0;
128 G4double sd = 0.0;
129 G4double r = 0.0; // relative err sd/mean/sqrt(n)
130 G4double efficiency = 0.0; // rate of non zero score
131 G4double r2eff = 0.0;
132 G4double r2int = 0.0;
133 G4double shift = 0.0;
134 G4double vov = 0.0;
135 G4double fom = 0.0;
136
137 G4double largest = 0.0;
138 G4int largest_score_happened = 0;
139
140 G4double mean_1 = 0.0;
141 G4double var_1 = 0.0;
142 G4double sd_1 = 0.0;
143 G4double r_1 = 0.0; // relative err sd/mean/sqrt(n)
144 G4double shift_1 = 0.0;
145 G4double vov_1 = 0.0;
146 G4double fom_1 = 0.0;
147
148 G4int noBinOfHistory = 16;
149 std::vector<G4int> history_grid;
150 std::vector<G4double> mean_history;
151 std::vector<G4double> var_history;
152 std::vector<G4double> sd_history;
153 std::vector<G4double> r_history;
154 std::vector<G4double> vov_history;
155 std::vector<G4double> fom_history;
156 std::vector<G4double> shift_history;
157 std::vector<G4double> e_history;
158 std::vector<G4double> r2eff_history;
159 std::vector<G4double> r2int_history;
160
161 G4double slope = 0.0;
162 std::vector<G4double> largest_scores;
163 std::vector<G4double> f_xi;
164 std::vector<G4double> f_yi;
165 G4int noBinOfPDF = 10;
166 G4SimplexDownhill<G4ConvergenceTester>* minimizer = nullptr;
167
168 G4int noPass = 0;
169 G4int noTotal = 8; // Total number of tests
170
171 G4bool statsAreUpdated = true;
172 G4bool showHistory = true;
173 G4bool calcSLOPE = true;
174 };
175
176 #endif
177