00001
00002
00003
00004
00005
00006 #include "RunAction.hh"
00007 #include "DetectorConstruction.hh"
00008
00009 #include "G4Run.hh"
00010 #include "G4RunManager.hh"
00011 #include "G4UnitsTable.hh"
00012
00013 #include <iomanip>
00014
00015
00016
00017 RunAction::RunAction(DetectorConstruction* det)
00018 :m_detector(det)
00019 { }
00020
00021
00022
00023 RunAction::~RunAction()
00024 { }
00025
00026
00027
00028 void RunAction::BeginOfRunAction(const G4Run* run)
00029 {
00030 G4cout << "### Run " << run->GetRunID() << " start." << G4endl;
00031
00032
00033
00034 m_sumEdep = m_sumEdep2 = 0.;
00035 m_nbEntry = 0;
00036
00037
00038
00039 m_nbSlices = 100;
00040 m_astronautHeight = m_detector->GetAstronautHeight();
00041 m_sliceWidth = m_astronautHeight/m_nbSlices;
00042 m_depthDose.resize(m_nbSlices, 0.);
00043 }
00044
00045
00046
00047 void RunAction::SumDepthDose(G4double zLocal, G4double edep)
00048 {
00049
00050 if (zLocal > m_astronautHeight) {
00051 G4cout << "\n --> warning from RunAction::SumDepthDose() : "
00052 << " zLocal beyond astronaut height : " << zLocal/mm << G4endl;
00053 return;
00054 }
00055
00056
00057 G4int k = int(zLocal/m_sliceWidth);
00058
00059
00060 m_depthDose[k] += edep;
00061 }
00062
00063
00064
00065 void RunAction::EndOfRunAction(const G4Run*)
00066 {
00067 G4double meanEd = 0., meanEd2 = 0., rmsEd = 0., error_rel_sum = 0.;
00068 if (m_nbEntry) {
00069 meanEd = m_sumEdep/m_nbEntry; meanEd2 = m_sumEdep2/m_nbEntry;
00070 G4double variance = meanEd2 - meanEd*meanEd;
00071 if (variance > 0.) rmsEd = std::sqrt(variance);
00072 error_rel_sum = rmsEd/(meanEd*std::sqrt(double(m_nbEntry)));
00073 }
00074
00075
00076 std::ios::fmtflags mode = G4cout.flags();
00077 G4cout.setf(std::ios::fixed,std::ios::floatfield);
00078 G4int prec = G4cout.precision(3);
00079
00080
00081
00082 G4cout
00083 << "\n ---> Nb of non-empty events = " << m_nbEntry
00084 << "\n ---> Mean Edep per event = " << G4BestUnit(meanEd,"Energy")
00085 << "\n ---> Total Energy deposited = " << G4BestUnit(m_sumEdep,"Energy")
00086 << " +- " << G4BestUnit(m_sumEdep*error_rel_sum,"Energy")
00087 << " (-> " << 100*error_rel_sum << " %)"
00088 << G4endl;
00089
00090
00091
00092 G4double mass = m_detector->GetAstronaut()->GetLogicalVolume()->GetMass();
00093
00094
00095
00096 G4double dose = m_sumEdep/mass;
00097 G4cout.setf(mode,std::ios::floatfield);
00098 G4cout.precision(5);
00099 G4cout << " ---> Dose = " << G4BestUnit(dose,"Dose")
00100 << " +- " << G4BestUnit(dose*error_rel_sum,"Dose") << "\n"
00101 << G4endl;
00102
00103
00104 G4cout.setf(mode,std::ios::floatfield);
00105 G4cout.precision(prec);
00106
00107
00108
00109 for (G4int k=0; k<m_nbSlices; k++) m_depthDose[k] /= mass;
00110
00111
00112
00113 WriteDepthDose("depthDose");
00114
00115
00116 G4cout.setf(mode,std::ios::floatfield);
00117 G4cout.precision(prec);
00118 }
00119
00120
00121
00122 #include <fstream>
00123
00124 void RunAction::WriteDepthDose(G4String name)
00125 {
00126 G4String fileName = name + ".ascii";
00127 std::ofstream File(fileName, std::ios::out);
00128
00129 std::ios::fmtflags mode = File.flags();
00130 File.setf( std::ios::scientific, std::ios::floatfield );
00131 G4int prec = File.precision(4);
00132
00133 File << " Longitudinal depth dose distribution \n "
00134 << "\n zLocal (mm)\t m_depthDose(gray) \n" << G4endl;
00135
00136 G4double zLocal;
00137 for (G4int k=0; k<m_nbSlices; k++) {
00138 zLocal = (k+0.5)*m_sliceWidth;
00139 File << " " << zLocal/mm
00140 << "\t " << m_depthDose[k]/gray << G4endl;
00141 }
00142
00143
00144 File.setf(mode,std::ios::floatfield);
00145 File.precision(prec);
00146 }
00147
00148
00149