00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #include "G4NuclWatcher.hh"
00034 #include "globals.hh"
00035
00036 #include <algorithm>
00037 #include <vector>
00038 #include <cmath>
00039
00040 G4NuclWatcher::G4NuclWatcher(G4int z,
00041 const std::vector<G4double>& expa,
00042 const std::vector<G4double>& expcs,
00043 const std::vector<G4double>& experr,
00044 G4bool check,
00045 G4bool nucl)
00046 : nuclz(z), izotop_chsq(0.), average_ratio(0.), aver_rat_err(0.),
00047 aver_lhood(0.), aver_matched(0.), exper_as(expa), exper_cs(expcs),
00048 exper_err(experr), checkable(check), nucleable(nucl) {}
00049
00050 void G4NuclWatcher::watch(G4int a, G4int z) {
00051 const G4double small = 0.001;
00052
00053 if (std::abs(z-nuclz) >= small) return;
00054
00055 G4bool here = false;
00056 G4int simulatedAsSize = simulated_as.size();
00057 for (G4int i = 0; i<simulatedAsSize && !here; i++) {
00058 if (std::abs(simulated_as[i] - a) < small) {
00059 simulated_cs[i] += 1.0;
00060 here = true;
00061 }
00062 }
00063
00064 if (!here) {
00065 simulated_as.push_back(a);
00066 simulated_cs.push_back(1.0);
00067 }
00068 }
00069
00070 void G4NuclWatcher::setInuclCs(G4double csec, G4int nev) {
00071 G4int simulatedAsSize = simulated_as.size();
00072 for(G4int i = 0; i < simulatedAsSize ; i++) {
00073 double err = std::sqrt(simulated_cs[i]) / simulated_cs[i];
00074
00075 simulated_prob.push_back(simulated_cs[i] / nev);
00076 simulated_cs[i] *= csec / nev;
00077 simulated_errors.push_back(simulated_cs[i] * err);
00078 }
00079 }
00080
00081 std::pair<G4double, G4double> G4NuclWatcher::getExpCs() const {
00082 G4double cs = 0.0;
00083 G4double err = 0.0;
00084
00085 G4int experAsSize = exper_as.size();
00086 for(G4int iz = 0; iz < experAsSize; iz++) {
00087 cs += exper_cs[iz];
00088 err += exper_err[iz];
00089 }
00090
00091 return std::pair<G4double, G4double>(cs, err);
00092 }
00093
00094 std::pair<G4double, G4double> G4NuclWatcher::getInuclCs() const {
00095 G4double cs = 0.0;
00096 G4double err = 0.0;
00097 G4int simulatedAsSize = simulated_as.size();
00098 for(G4int iz = 0; iz < simulatedAsSize; iz++) {
00099 cs += simulated_cs[iz];
00100 err += simulated_errors[iz];
00101 }
00102
00103 return std::pair<G4double, G4double>(cs, err);
00104 }
00105
00106 void G4NuclWatcher::print() {
00107 const G4double small = 0.001;
00108
00109 G4cout << "\n ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ "
00110 << "\n **** izotop Z **** " << nuclz << G4endl;
00111
00112 izotop_chsq = 0.0;
00113 average_ratio = 0.0;
00114 aver_rat_err = 0.0;
00115 G4double exp_cs = 0.0;
00116 G4double exp_cs_err = 0.0;
00117 G4double inucl_cs = 0.0;
00118 G4double inucl_cs_err = 0.0;
00119 std::vector<G4bool> not_used(simulated_cs.size(), true);
00120 G4int nmatched = exper_as.size();
00121 G4int nused = simulated_cs.size();
00122 G4double lhood = 0.0;
00123 G4int experAsSize = exper_as.size();
00124
00125 for (G4int iz = 0; iz < experAsSize; iz++) {
00126 G4double a = exper_as[iz];
00127
00128 exp_cs += exper_cs[iz];
00129 exp_cs_err += exper_err[iz];
00130
00131 G4bool found = false;
00132 G4int simulatedAsSize = simulated_as.size();
00133 for (G4int i = 0; i<simulatedAsSize && !found; i++) {
00134 if (std::fabs(simulated_as[i] - a) < small) {
00135 G4double rat = simulated_cs[i] / exper_cs[iz];
00136
00137 lhood += std::log10(rat) * std::log10(rat);
00138
00139 G4double rat_err
00140 = std::sqrt(simulated_errors[i]*simulated_errors[i] +
00141 exper_err[iz]*exper_err[iz] * rat*rat) / exper_cs[iz];
00142 average_ratio += rat;
00143 aver_rat_err += rat_err;
00144
00145 G4cout << " A " << a << " exp.cs " << exper_cs[iz] << " err "
00146 << exper_err[iz] << G4endl << " sim. cs " << simulated_cs[i]
00147 << " err " << simulated_errors[i] << G4endl
00148 << " ratio " << rat << " err " << rat_err << G4endl
00149 << " simulated production rate " << simulated_prob[i] << G4endl;
00150
00151 not_used[i] = false;
00152 izotop_chsq += (rat - 1.0) * (rat - 1.0) / rat_err / rat_err;
00153 found = true;
00154 nused--;
00155 }
00156 }
00157
00158 if (found) nmatched--;
00159 else
00160 G4cout << " not found exper.: A " << a << " exp.cs " << exper_cs[iz]
00161 << " err " << exper_err[iz] << G4endl;
00162 }
00163
00164 G4cout << " not found in simulations " << nmatched << G4endl
00165 << " not found in exper: " << nused << G4endl;
00166
00167 G4int simulatedAsSize = simulated_as.size();
00168 for(G4int i = 0; i < simulatedAsSize; i++) {
00169 inucl_cs += simulated_cs[i];
00170 inucl_cs_err += simulated_errors[i];
00171
00172 if(not_used[i])
00173 G4cout << " extra simul.: A " << simulated_as[i] << " sim. cs "
00174 << simulated_cs[i] << " err " << simulated_errors[i] << G4endl;
00175
00176 G4cout << " simulated production rate " << simulated_prob[i] << G4endl;
00177 }
00178
00179 G4int matched = exper_as.size() - nmatched;
00180
00181 if (matched > 0) {
00182 aver_lhood = lhood;
00183 aver_matched = matched;
00184 lhood = std::pow(10.0, std::sqrt(lhood/matched));
00185
00186 G4cout << " matched " << matched << " CHSQ " << std::sqrt(izotop_chsq) / matched
00187 << G4endl
00188 << " raw chsq " << izotop_chsq << G4endl
00189 << " average ratio " << average_ratio / matched
00190 << " err " << aver_rat_err / matched << G4endl
00191 << " lhood " << lhood << G4endl;
00192 }
00193 else {
00194 izotop_chsq = 0.0;
00195 aver_lhood = 0.0;
00196 }
00197
00198 G4cout << " exper. cs " << exp_cs << " err " << exp_cs_err << G4endl
00199 << " inucl. cs " << inucl_cs << " err " << inucl_cs_err << G4endl
00200 << " ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ "
00201 << G4endl;
00202 }