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 #include "G4Analyser.hh"
00033 #include <cmath>
00034 #include <iomanip>
00035
00036 G4Analyser::G4Analyser()
00037 : verboseLevel(0), eventNumber(0.0), averageMultiplicity(0.0),
00038 averageProtonNumber(0.0), averageNeutronNumber(0.0),
00039 averagePionNumber(0.0), averageNucleonKinEnergy(0.0),
00040 averageProtonKinEnergy(0.0), averageNeutronKinEnergy(0.0),
00041 averagePionKinEnergy(0.0), averageExitationEnergy(0.0),
00042 averageOutgoingNuclei(0.0), fissy_prob(0.0), averagePionPl(0.0),
00043 averagePionMin(0.0), averagePion0(0.0), averageA(0.0), averageZ(0.0),
00044 inel_csec(0.0), withNuclei(false) {
00045 if (verboseLevel > 3) {
00046 G4cout << " >>> G4Analyser::G4Analyser" << G4endl;
00047 }
00048 }
00049
00050 void G4Analyser::setInelCsec(G4double csec, G4bool withn) {
00051
00052 if (verboseLevel > 3) {
00053 G4cout << " >>> G4Analyser::setInelCsec" << G4endl;
00054 }
00055
00056 inel_csec = csec;
00057 withNuclei = withn;
00058
00059 if (verboseLevel > 3) {
00060 G4cout << " total inelastic " << inel_csec << G4endl;
00061 }
00062 }
00063
00064 void G4Analyser::setWatchers(const std::vector<G4NuclWatcher>& watchers) {
00065
00066 if (verboseLevel > 3) {
00067 G4cout << " >>> G4Analyser::setWatchers" << G4endl;
00068 }
00069
00070 ana_watchers = watchers;
00071 if (verboseLevel > 3) {
00072 G4cout << " watchers set " << watchers.size() << G4endl;
00073 }
00074 }
00075
00076 void G4Analyser::try_watchers(G4int a, G4int z, G4bool if_nucl) {
00077
00078 if (verboseLevel > 3) {
00079 G4cout << " >>> G4Analyser::try_watchers" << G4endl;
00080 }
00081
00082 for (G4int iw = 0; iw < G4int(ana_watchers.size()); iw++) {
00083
00084 if (if_nucl) {
00085
00086 if (ana_watchers[iw].look_forNuclei()) ana_watchers[iw].watch(a, z);
00087
00088 } else {
00089
00090 if (!ana_watchers[iw].look_forNuclei()) ana_watchers[iw].watch(a, z);
00091 }
00092 }
00093 }
00094
00095
00096 void G4Analyser::analyse(const G4CollisionOutput& output) {
00097
00098 if (verboseLevel > 3) {
00099 G4cout << " >>> G4Analyser::analyse" << G4endl;
00100 }
00101
00102 if (withNuclei) {
00103 const std::vector<G4InuclNuclei>& nucleus = output.getOutgoingNuclei();
00104
00105
00106 if (nucleus.size() > 0) {
00107 G4int nbig = 0;
00108 averageOutgoingNuclei += nucleus.size();
00109
00110 for (G4int in = 0; in < G4int(nucleus.size()); in++) {
00111 averageExitationEnergy += nucleus[in].getExitationEnergy();
00112
00113 G4int a = nucleus[in].getA();
00114 G4int z = nucleus[in].getZ();
00115
00116 if (in == 0) {
00117 averageA += a;
00118 averageZ += z;
00119 };
00120
00121 if (a > 10) nbig++;
00122 try_watchers(a, z, true);
00123 };
00124
00125 if (nbig > 1) fissy_prob += 1.0;
00126 eventNumber += 1.0;
00127 const std::vector<G4InuclElementaryParticle>& particles =
00128 output.getOutgoingParticles();
00129 averageMultiplicity += particles.size();
00130
00131 for (G4int i = 0; i < G4int(particles.size()); i++) {
00132 G4int ap = 0;
00133 G4int zp = 0;
00134
00135 if (particles[i].nucleon()) {
00136 averageNucleonKinEnergy += particles[i].getKineticEnergy();
00137
00138 if (particles[i].type() == 1) {
00139 zp = 1;
00140 ap = 1;
00141 averageProtonNumber += 1.0;
00142 averageProtonKinEnergy += particles[i].getKineticEnergy();
00143
00144 } else {
00145 ap = 1;
00146 zp = 0;
00147 averageNeutronNumber += 1.0;
00148 averageNeutronKinEnergy += particles[i].getKineticEnergy();
00149 };
00150
00151 } else if (particles[i].pion()) {
00152 averagePionKinEnergy += particles[i].getKineticEnergy();
00153 averagePionNumber += 1.0;
00154 ap = 0;
00155
00156 if (particles[i].type() == 3) {
00157 zp = 1;
00158 averagePionPl += 1.0;
00159
00160 } else if (particles[i].type() == 5) {
00161 zp = -1;
00162 averagePionMin += 1.0;
00163
00164 } else if (particles[i].type() == 7) {
00165 zp = 0;
00166 averagePion0 += 1.0;
00167 };
00168 };
00169 try_watchers(ap, zp, false);
00170 };
00171 };
00172
00173 } else {
00174 eventNumber += 1.0;
00175 const std::vector<G4InuclElementaryParticle>& particles =
00176 output.getOutgoingParticles();
00177 averageMultiplicity += particles.size();
00178
00179 for (G4int i = 0; i < G4int(particles.size()); i++) {
00180
00181 if (particles[i].nucleon()) {
00182 averageNucleonKinEnergy += particles[i].getKineticEnergy();
00183
00184 if (particles[i].type() == 1) {
00185 averageProtonNumber += 1.0;
00186 averageProtonKinEnergy += particles[i].getKineticEnergy();
00187
00188 } else {
00189 averageNeutronNumber += 1.0;
00190 averageNeutronKinEnergy += particles[i].getKineticEnergy();
00191 }
00192
00193 } else if (particles[i].pion()) {
00194 averagePionKinEnergy += particles[i].getKineticEnergy();
00195 averagePionNumber += 1.0;
00196 }
00197 }
00198 }
00199 }
00200
00201 void G4Analyser::printResultsSimple() {
00202
00203 if (verboseLevel > 3) {
00204 G4cout << " >>> G4Analyser::printResultsSimple" << G4endl;
00205 }
00206
00207 G4cout << " Number of events " << G4int(eventNumber + 0.1) << G4endl
00208 << " average multiplicity " << averageMultiplicity / eventNumber << G4endl
00209 << " average proton number " << averageProtonNumber / eventNumber << G4endl
00210 << " average neutron number " << averageNeutronNumber / eventNumber << G4endl
00211 << " average nucleon Ekin " << averageNucleonKinEnergy /
00212 (averageProtonNumber + averageNeutronNumber) << G4endl
00213 << " average proton Ekin " << averageProtonKinEnergy / (averageProtonNumber +
00214 1.0e-10) << G4endl
00215 << " average neutron Ekin " << averageNeutronKinEnergy / (averageNeutronNumber +
00216 1.0e-10) << G4endl
00217 << " average pion number " << averagePionNumber / eventNumber << G4endl
00218 << " average pion Ekin " << averagePionKinEnergy / (averagePionNumber +
00219 1.0e-10) << G4endl;
00220 if (withNuclei) {
00221 G4cout
00222 << " average Exitation Energy " <<
00223 averageExitationEnergy / averageOutgoingNuclei << G4endl
00224 << " average num of fragments " << averageOutgoingNuclei / eventNumber << G4endl;
00225 G4cout << " fission prob. " << fissy_prob / eventNumber << " c.sec " <<
00226 inel_csec * fissy_prob / eventNumber << G4endl;
00227 }
00228 }
00229
00230 void G4Analyser::printResults() {
00231
00232 if (verboseLevel > 3) {
00233 G4cout << " >>> G4Analyser::printResults" << G4endl;
00234 }
00235
00236 G4cout << " Number of events " << G4int(eventNumber + 0.1) << G4endl
00237 << " average multiplicity " << averageMultiplicity / eventNumber << G4endl
00238 << " average proton number " << averageProtonNumber / eventNumber << G4endl
00239 << " average neutron number " << averageNeutronNumber / eventNumber << G4endl
00240 << " average nucleon Ekin " << averageNucleonKinEnergy /
00241 (averageProtonNumber + averageNeutronNumber) << G4endl
00242 << " average proton Ekin " << averageProtonKinEnergy / (averageProtonNumber +
00243 1.0e-10) << G4endl
00244 << " average neutron Ekin " << averageNeutronKinEnergy / (averageNeutronNumber +
00245 1.0e-10) << G4endl
00246 << " average pion number " << averagePionNumber / eventNumber << G4endl
00247 << " average pion Ekin " << averagePionKinEnergy / (averagePionNumber +
00248 1.0e-10) << G4endl
00249 << " average pi+ " << averagePionPl / eventNumber << G4endl
00250 << " average pi- " << averagePionMin / eventNumber << G4endl
00251 << " average pi0 " << averagePion0 / eventNumber << G4endl;
00252
00253 if (withNuclei) {
00254 G4cout
00255 << " average A " << averageA / eventNumber << G4endl
00256 << " average Z " << averageZ / eventNumber << G4endl
00257 << " average Exitation Energy " <<
00258 averageExitationEnergy / averageOutgoingNuclei << G4endl
00259 << " average num of fragments " << averageOutgoingNuclei / eventNumber << G4endl;
00260 G4cout << " fission prob. " << fissy_prob / eventNumber << " c.sec " <<
00261 inel_csec * fissy_prob / eventNumber << G4endl;
00262 handleWatcherStatistics();
00263 }
00264 }
00265
00266 void G4Analyser::handleWatcherStatistics() {
00267
00268 if (verboseLevel > 3) {
00269 G4cout << " >>> G4Analyser::handleWatcherStatistics" << G4endl;
00270 }
00271
00272
00273
00274 if (verboseLevel > 3) {
00275 G4cout << " >>>Izotop analysis:" << G4endl;
00276 };
00277
00278 G4double fgr = 0.0;
00279 G4double averat = 0.0;
00280 G4double ave_err = 0.0;
00281 G4double gl_chsq = 0.0;
00282 G4double tot_exper = 0.0;
00283 G4double tot_exper_err = 0.0;
00284 G4double tot_inucl = 0.0;
00285 G4double tot_inucl_err = 0.0;
00286 G4double checked = 0.0;
00287
00288 for (G4int iw = 0; iw < G4int(ana_watchers.size()); iw++) {
00289 ana_watchers[iw].setInuclCs(inel_csec, G4int(eventNumber));
00290 ana_watchers[iw].print();
00291
00292 if (ana_watchers[iw].to_check()) {
00293 std::pair<G4double, G4double> rat_err = ana_watchers[iw].getAverageRatio();
00294 averat += rat_err.first;
00295 ave_err += rat_err.second;
00296 gl_chsq += ana_watchers[iw].getChsq();
00297 std::pair<G4double, G4double> cs_err = ana_watchers[iw].getExpCs();
00298 tot_exper += cs_err.first;
00299 tot_exper_err += cs_err.second;
00300 std::pair<G4double, G4double> inucl_cs_err = ana_watchers[iw].getInuclCs();
00301 tot_inucl += inucl_cs_err.first;
00302 tot_inucl_err += inucl_cs_err.second;
00303 G4double iz_checked = ana_watchers[iw].getNmatched();
00304
00305 if (iz_checked > 0.0) {
00306 fgr += ana_watchers[iw].getLhood();
00307 checked += iz_checked;
00308 };
00309 };
00310 };
00311
00312 if (checked > 0.0) {
00313 gl_chsq = std::sqrt(gl_chsq) / checked;
00314 averat /= checked;
00315 ave_err /= checked;
00316 fgr = std::pow(10.0, std::sqrt(fgr / checked));
00317 };
00318
00319 if (verboseLevel > 3) {
00320 G4cout << " total exper c.s. " << tot_exper << " err " << tot_exper_err <<
00321 " tot inucl c.s. " << tot_inucl << " err " << tot_inucl_err << G4endl;
00322 G4cout << " checked total " << checked << " lhood " << fgr << G4endl
00323 << " average ratio " << averat << " err " << ave_err << G4endl
00324 << " global chsq " << gl_chsq << G4endl;
00325 }
00326 }
00327
00328 void G4Analyser::printResultsNtuple() {
00329
00330 if (verboseLevel > 3) {
00331 G4cout << " >>> G4Analyser::printResultsNtuple" << G4endl;
00332 }
00333
00334
00335
00336 G4cout <<
00337 std::setw(15) << int(eventNumber + 0.1) <<
00338 std::setw(15) << averageMultiplicity / eventNumber <<
00339 std::setw(15) << averageProtonNumber / eventNumber <<
00340 std::setw(15) << averageNeutronNumber / eventNumber << " " <<
00341 std::setw(15) << averageNucleonKinEnergy / (averageProtonNumber + averageNeutronNumber) << " " <<
00342 std::setw(15) << averageProtonKinEnergy / (averageProtonNumber + 1.0e-10) << " " <<
00343 std::setw(15) << averageNeutronKinEnergy / (averageNeutronNumber + 1.0e-10) << " " <<
00344 std::setw(15) << averagePionNumber / eventNumber << " " <<
00345 std::setw(15) << averagePionKinEnergy / (averagePionNumber + 1.0e-10) << G4endl;
00346 }
