| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 /// \file biasing/ReverseMC01/src/RMC01Analysi << 26 // $Id: RMC01AnalysisManager.cc,v 1.7 2010-11-11 14:39:42 ldesorgh Exp $
27 /// \brief Implementation of the RMC01Analysis << 27 // GEANT4 tag $Name: not supported by cvs2svn $
28 // <<
29 // 28 //
30 ////////////////////////////////////////////// 29 //////////////////////////////////////////////////////////////
31 // Class Name: RMC01AnalysisManage << 30 // Class Name: RMC01AnalysisManager
32 // Author: L. Desorgher << 31 // Author: L. Desorgher
33 // Organisation: SpaceIT GmbH << 32 // Organisation: SpaceIT GmbH
34 // Contract: ESA contract 21435/ << 33 // Contract: ESA contract 21435/08/NL/AT
35 // Customer: ESA/ESTEC << 34 // Customer: ESA/ESTEC
36 ////////////////////////////////////////////// 35 //////////////////////////////////////////////////////////////
37 <<
38 //....oooOO0OOooo........oooOO0OOooo........oo <<
39 //....oooOO0OOooo........oooOO0OOooo........oo <<
40 <<
41 #include "RMC01AnalysisManager.hh" 36 #include "RMC01AnalysisManager.hh"
42 <<
43 #include "RMC01AnalysisManagerMessenger.hh" <<
44 #include "RMC01SD.hh" <<
45 <<
46 #include "G4AdjointSimManager.hh" 37 #include "G4AdjointSimManager.hh"
47 #include "G4Electron.hh" <<
48 #include "G4Gamma.hh" <<
49 #include "G4PhysicalConstants.hh" <<
50 #include "G4Proton.hh" <<
51 #include "G4RunManager.hh" <<
52 #include "G4SDManager.hh" 38 #include "G4SDManager.hh"
53 #include "G4SystemOfUnits.hh" << 39 #include "RMC01SD.hh"
54 #include "G4THitsCollection.hh" 40 #include "G4THitsCollection.hh"
>> 41 #include "G4Electron.hh"
>> 42 #include "G4Proton.hh"
>> 43 #include "G4Gamma.hh"
>> 44 #include "Histo1DVar.hh"
>> 45 #include "Histo2DVar.hh"
55 #include "G4Timer.hh" 46 #include "G4Timer.hh"
>> 47 #include "G4RunManager.hh"
>> 48 #include "RMC01AnalysisManagerMessenger.hh"
>> 49 RMC01AnalysisManager* RMC01AnalysisManager::instance = 0;
56 50
57 RMC01AnalysisManager* RMC01AnalysisManager::fI << 51 ////////////////////////////////////////////////////////////////////////////////
58 << 52 //
59 //....oooOO0OOooo........oooOO0OOooo........oo << 53 RMC01AnalysisManager::RMC01AnalysisManager()
60 << 54 {
61 RMC01AnalysisManager::RMC01AnalysisManager() << 55
62 : fAccumulated_edep(0.), << 56 theMsg = new RMC01AnalysisManagerMessenger(this);
63 fAccumulated_edep2(0.), << 57
64 fMean_edep(0.), << 58 //----------------------
65 fError_mean_edep(0.), << 59 //Creation of histograms
66 fRelative_error(0.), << 60 //----------------------
67 fElapsed_time(0.), << 61
68 fPrecision_to_reach(0.), << 62 //Energy binning of the histograms : 60 log bins over [1keV-1GeV]
69 fStop_run_if_precision_reached(true), << 63
70 fNb_evt_modulo_for_convergence_test(5000), << 64 G4double bins[61];
71 fEdep_rmatrix_vs_electron_prim_energy(0), << 65 size_t nbin=60;
72 fElectron_current_rmatrix_vs_electron_prim << 66 G4double Emin=1.*keV;
73 fGamma_current_rmatrix_vs_electron_prim_en << 67 G4double Emax=1.*GeV;
74 fEdep_rmatrix_vs_gamma_prim_energy(0), << 68 for ( size_t i=0; i <= nbin; i++) {
75 fElectron_current_rmatrix_vs_gamma_prim_en << 69 G4double val_bin;
76 fGamma_current_rmatrix_vs_gamma_prim_energ << 70 val_bin=Emin * std::pow(10., i * std::log10(Emax/Emin)/nbin);
77 fEdep_rmatrix_vs_proton_prim_energy(0), << 71 G4double exp_10=4.-int(std::log10(val_bin));
78 fElectron_current_rmatrix_vs_proton_prim_e << 72 G4double factor =std::pow(10., exp_10);
79 fProton_current_rmatrix_vs_proton_prim_ene << 73 val_bin=int(factor*val_bin)/factor;
80 fGamma_current_rmatrix_vs_proton_prim_ener << 74 bins[i] = val_bin;
81 fFactoryOn(false), << 75
82 fPrimSpectrumType(EXPO), << 76 }
83 fAlpha_or_E0(.5 * MeV), << 77
84 fAmplitude_prim_spectrum(1.), << 78
85 fEmin_prim_spectrum(1. * keV), << 79 //Histograms for :
86 fEmax_prim_spectrum(20. * MeV), << 80 // 1)the forward simulation results
87 fAdjoint_sim_mode(true), << 81 // 2)the Reverse MC simulation results normalised to a user spectrum
88 fNb_evt_per_adj_evt(2) << 82 //-----------------------------------------------------------------------------------
89 { << 83
90 fMsg = new RMC01AnalysisManagerMessenger(thi << 84 edep_vs_prim_ekin = new Histo1DVar("edep_vs_prim_ekin",bins, nbin+1, LEFT);
91 << 85 electron_current = new Histo1DVar("electron_current",bins, nbin+1, LEFT);
>> 86 proton_current= new Histo1DVar("proton_current",bins, nbin+1, LEFT);
>> 87 gamma_current= new Histo1DVar("gamma_current",bins, nbin+1, LEFT);
>> 88
>> 89
>> 90
>> 91 //Answer matrices for the adjoint simulation only
>> 92 //-----------------------------------------------
>> 93
>> 94
>> 95 //answer matrices for external isotropic e- source
>> 96
>> 97 edep_answer_matrix_vs_electron_prim_energy = new Histo1DVar("edep_answer_matrix_vs_electron_prim_energy", bins, nbin+1, LEFT);
>> 98
>> 99 electron_current_answer_matrix_vs_electron_prim_energy = new Histo2DVar("electron_current_answer_matrix_vs_electron_prim_energy", bins, nbin+1, LEFT, bins, nbin+1, LEFT);
>> 100
>> 101 gamma_current_answer_matrix_vs_electron_prim_energy = new Histo2DVar("gamma_current_answer_matrix_vs_electron_prim_energy", bins, nbin+1, LEFT, bins, nbin+1, LEFT);
>> 102
>> 103
>> 104 //answer matrices for external isotropic gamma source
>> 105
>> 106 edep_answer_matrix_vs_gamma_prim_energy = new Histo1DVar("edep_answer_matrix_vs_gamma_prim_energy", bins, nbin+1, LEFT);
>> 107
>> 108 electron_current_answer_matrix_vs_gamma_prim_energy = new Histo2DVar("electron_current_answer_matrix_vs_gamma_prim_energy", bins, nbin+1, LEFT, bins, nbin+1, LEFT);
>> 109
>> 110 gamma_current_answer_matrix_vs_gamma_prim_energy = new Histo2DVar("gamma_current_answer_matrix_vs_gamma_prim_energy", bins, nbin+1, LEFT, bins, nbin+1, LEFT);
>> 111
>> 112
>> 113
>> 114 //answer matrices for external isotropic proton source
>> 115
>> 116 edep_answer_matrix_vs_proton_prim_energy = new Histo1DVar("edep_answer_matrix_vs_proton_prim_energy", bins, nbin+1, LEFT);
>> 117
>> 118 electron_current_answer_matrix_vs_proton_prim_energy = new Histo2DVar("electron_current_answer_matrix_vs_proton_prim_energy", bins, nbin+1, LEFT, bins, nbin+1, LEFT);
>> 119
>> 120 gamma_current_answer_matrix_vs_proton_prim_energy = new Histo2DVar("gamma_current_answer_matrix_vs_proton_prim_energy", bins, nbin+1, LEFT, bins, nbin+1, LEFT);
>> 121
>> 122 proton_current_answer_matrix_vs_proton_prim_energy = new Histo2DVar("proton_current_answer_matrix_vs_proton_prim_energy", bins, nbin+1, LEFT, bins, nbin+1, LEFT);
>> 123
>> 124
>> 125
>> 126 //--------------------------
>> 127 // Convergence test variables
>> 128 //--------------------------
>> 129
>> 130 precision_to_reach =0.;
>> 131 stop_run_if_precision_reached =true;
>> 132 nb_evt_modulo_for_convergence_test =5000;
>> 133
>> 134
>> 135
92 //------------- 136 //-------------
93 // Timer for convergence vector << 137 //Timer for convergence vector
94 //------------- 138 //-------------
>> 139
>> 140 theTimer = new G4Timer();
>> 141
>> 142
>> 143
>> 144
>> 145
95 146
96 fTimer = new G4Timer(); << 147
97 <<
98 //--------------------------------- 148 //---------------------------------
99 // Primary particle ID for normalisation of << 149 //Primary spectrum for normalisation of adjoint results
100 //--------------------------------- 150 //---------------------------------
101 << 151
102 fPrimPDG_ID = G4Electron::Electron()->GetPDG << 152 thePrimPDG_ID = G4Electron::Electron()->GetPDGEncoding();
103 << 153 alpha_or_E0=.5*MeV;
104 fFileName[0] = "sim"; << 154 amplitude_prim_spectrum =1.;
>> 155 emin_prim_spectrum=1.*keV;
>> 156 emax_prim_spectrum = 20.*MeV;
>> 157
>> 158
>> 159 }
>> 160 ////////////////////////////////////////////////////////////////////////////////
>> 161 //
>> 162 RMC01AnalysisManager::~RMC01AnalysisManager()
>> 163 {
>> 164 delete edep_vs_prim_ekin;
>> 165 delete electron_current;
>> 166 delete proton_current;
>> 167 delete gamma_current;
>> 168
>> 169 delete edep_answer_matrix_vs_electron_prim_energy;
>> 170 delete electron_current_answer_matrix_vs_electron_prim_energy;
>> 171 delete gamma_current_answer_matrix_vs_electron_prim_energy;
>> 172
>> 173 delete edep_answer_matrix_vs_gamma_prim_energy;
>> 174 delete electron_current_answer_matrix_vs_gamma_prim_energy;
>> 175 delete gamma_current_answer_matrix_vs_gamma_prim_energy;
>> 176
>> 177 delete edep_answer_matrix_vs_proton_prim_energy;
>> 178 delete electron_current_answer_matrix_vs_proton_prim_energy;
>> 179 delete proton_current_answer_matrix_vs_proton_prim_energy;
>> 180 delete gamma_current_answer_matrix_vs_proton_prim_energy;
>> 181
>> 182
105 } 183 }
106 << 184 ////////////////////////////////////////////////////////////////////////////////
107 //....oooOO0OOooo........oooOO0OOooo........oo << 185 //
108 <<
109 RMC01AnalysisManager::~RMC01AnalysisManager() <<
110 { <<
111 ; <<
112 } <<
113 <<
114 //....oooOO0OOooo........oooOO0OOooo........oo <<
115 <<
116 RMC01AnalysisManager* RMC01AnalysisManager::Ge 186 RMC01AnalysisManager* RMC01AnalysisManager::GetInstance()
117 { 187 {
118 if (fInstance == 0) fInstance = new RMC01Ana << 188 if (instance == 0) instance = new RMC01AnalysisManager;
119 return fInstance; << 189 return instance;
120 } 190 }
121 << 191 ////////////////////////////////////////////////////////////////////////////////
122 //....oooOO0OOooo........oooOO0OOooo........oo << 192 //
123 <<
124 void RMC01AnalysisManager::BeginOfRun(const G4 193 void RMC01AnalysisManager::BeginOfRun(const G4Run* aRun)
125 { << 194 { accumulated_edep =0.;
126 fAccumulated_edep = 0.; << 195 accumulated_edep2 =0.;
127 fAccumulated_edep2 = 0.; << 196 relative_error=1.;
128 fNentry = 0.0; << 197 mean_edep=0.;
129 fRelative_error = 1.; << 198 error_mean_edep=0.;
130 fMean_edep = 0.; << 199 adjoint_sim_mode =G4AdjointSimManager::GetInstance()->GetAdjointSimMode();
131 fError_mean_edep = 0.; << 200 if (adjoint_sim_mode){
132 fAdjoint_sim_mode = G4AdjointSimManager::Get << 201 nb_evt_per_adj_evt=aRun->GetNumberOfEventToBeProcessed()/G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
133 << 202 ConvergenceFileOutput.open("ConvergenceOfAdjointSimulationResults.txt", std::ios::out);
134 if (fAdjoint_sim_mode) { << 203 ConvergenceFileOutput<<"Normalised Edep[MeV]\terror[MeV]\tcomputing_time[s]"<<std::endl;
135 fNb_evt_per_adj_evt = aRun->GetNumberOfEve << 204 }
136 / G4AdjointSimManage << 205 else {
137 fConvergenceFileOutput.open("ConvergenceOf << 206 ConvergenceFileOutput.open("ConvergenceOfForwardSimulationResults.txt", std::ios::out);
138 fConvergenceFileOutput << "Normalised Edep << 207 ConvergenceFileOutput<<"Edep per event [MeV]\terror[MeV]\tcomputing_time[s]"<<std::endl;
139 } << 208 }
140 else { << 209 ConvergenceFileOutput.setf(std::ios::scientific);
141 fConvergenceFileOutput.open("ConvergenceOf << 210 ConvergenceFileOutput.precision(6);
142 fConvergenceFileOutput << "Edep per event << 211 ResetHistograms();
143 } << 212 theTimer->Start();
144 fConvergenceFileOutput.setf(std::ios::scient << 213 elapsed_time=0.;
145 fConvergenceFileOutput.precision(6); <<
146 <<
147 fTimer->Start(); <<
148 fElapsed_time = 0.; <<
149 <<
150 Book(); <<
151 } 214 }
152 << 215 ////////////////////////////////////////////////////////////////////////////////
153 //....oooOO0OOooo........oooOO0OOooo........oo << 216 //
154 <<
155 void RMC01AnalysisManager::EndOfRun(const G4Ru 217 void RMC01AnalysisManager::EndOfRun(const G4Run* aRun)
156 { << 218 { theTimer->Stop();
157 fTimer->Stop(); << 219
158 G4int nb_evt = aRun->GetNumberOfEvent(); << 220
159 G4double factor = 1. / nb_evt; << 221 if (!adjoint_sim_mode){
160 if (!fAdjoint_sim_mode) { << 222
161 G4cout << "Results of forward simulation!" << 223 G4cout<<"Results of forward simulation!"<<std::endl;
162 G4cout << "edep per event [MeV] = " << fMe << 224 G4cout<<"edep per event [MeV] = "<<mean_edep<<std::endl;
163 G4cout << "error[MeV] = " << fError_mean_e << 225 G4cout<<"error[MeV] = "<<error_mean_edep<<std::endl;
>> 226 G4int nb_evt=aRun->GetNumberOfEvent();
>> 227 WriteHisto(edep_vs_prim_ekin,1./nb_evt, G4String("Fwd_Edep_vs_EkinPrim.txt"),G4String("E1[MeV]\t\tE2[MeV]\t\tEdep[MeV]\terr_Edep[MeV]\n"));
>> 228 WriteHisto(electron_current,1./nb_evt, G4String("Fwd_ElectronCurrent.txt"),G4String("E1[MeV]\t\tE2[MeV]\t\tnb entering electron\terr\n"));
>> 229 WriteHisto(proton_current,1./nb_evt, G4String("Fwd_ProtonCurrent.txt"),G4String("E1[MeV]\t\tE2[MeV]\t\tnb entering proton\terr[]\n"));
>> 230 WriteHisto(gamma_current,1./nb_evt, G4String("Fwd_GammaCurrent.txt"),G4String("E1[MeV]\t\tE2[MeV]\t\tnb entering gamma\terr[]\n"));
164 } 231 }
165 << 232
166 else { 233 else {
167 G4cout << "Results of reverse/adjoint simu << 234 G4cout<<"Results of reverse/adjoint simulation!"<<std::endl;
168 G4cout << "normalised edep [MeV] = " << fM << 235 G4cout<<"normalised edep [MeV] = "<<mean_edep<<std::endl;
169 G4cout << "error[MeV] = " << fError_mean_e << 236 G4cout<<"error[MeV] = "<<error_mean_edep<<std::endl;
170 factor = 1. * G4AdjointSimManager::GetInst << 237
171 / aRun->GetNumberOfEvent(); << 238
172 } << 239 G4double factor=1.*G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun()*nb_evt_per_adj_evt/aRun->GetNumberOfEvent();
173 Save(factor); << 240
174 fConvergenceFileOutput.close(); << 241
>> 242
>> 243 WriteHisto(edep_vs_prim_ekin,factor, G4String("Adj_Edep_vs_EkinPrim.txt"),G4String("E1[MeV]\t\tE2[MeV]\t\tEdep[MeV]\terr_Edep[MeV]\n"));
>> 244 WriteHisto(electron_current,factor, G4String("Adj_ElectronCurrent.txt"),G4String("E1[MeV]\t\tE2[MeV]\t\tnb entering electron\terr\n"));
>> 245 WriteHisto(proton_current,factor, G4String("Adj_ProtonCurrent.txt"),G4String("E1[MeV]\t\tE2[MeV]\t\tnb entering proton\terr[]\n"));
>> 246 WriteHisto(gamma_current,factor, G4String("Adj_GammaCurrent.txt"),G4String("E1[MeV]\t\tE2[MeV]\t\tnb entering gamma\terr[]\n"));
>> 247
>> 248
>> 249 WriteHisto(edep_answer_matrix_vs_electron_prim_energy,factor,
>> 250 G4String("Adj_Edep_vs_EkinPrimElectron_Answer.txt"),
>> 251 G4String("E1[MeV]\t\tE2[MeV]\t\tEdep Efficiency[MeV*cm2*MeV*str]\terr_Edep[MeV*cm2*MeV*str]\n"));
>> 252
>> 253 WriteHisto(electron_current_answer_matrix_vs_electron_prim_energy,factor,
>> 254 G4String("Adj_ElectronCurrent_vs_EkinPrimElectron_Answer.txt"),
>> 255 G4String("Eprim1[MeV]\t\tEprim2[MeV]\t\tEsec1[MeV]\t\tEsec2[MeV]\t Current Efficiency[cm2*MeV*str]\terr[cm2*MeV*str]\n"));
>> 256
>> 257 WriteHisto(gamma_current_answer_matrix_vs_electron_prim_energy,factor,
>> 258 G4String("Adj_GammaCurrent_vs_EkinPrimElectron_Answer.txt"),
>> 259 G4String("Eprim1[MeV]\t\tEprim2[MeV]\t\tEsec1[MeV]\t\tEsec2[MeV]\t Current Efficiency[cm2*MeV*str]\terr[cm2*MeV*str]\n"));
>> 260
>> 261
>> 262
>> 263 WriteHisto(edep_answer_matrix_vs_gamma_prim_energy,factor,
>> 264 G4String("Adj_Edep_vs_EkinPrimGamma_Answer.txt"),
>> 265 G4String("E1[MeV]\t\tE2[MeV]\t\tEdep Efficiency[MeV*cm2*MeV*str]\terr_Edep[MeV*cm2*MeV*str]\n"));
>> 266
>> 267 WriteHisto(electron_current_answer_matrix_vs_gamma_prim_energy,factor,
>> 268 G4String("Adj_ElectronCurrent_vs_EkinPrimGamma_Answer.txt"),
>> 269 G4String("Eprim1[MeV]\t\tEprim2[MeV]\t\tEsec1[MeV]\t\tEsec2[MeV]\t Current Efficiency[cm2*MeV*str]\terr[cm2*MeV*str]\n"));
>> 270
>> 271 WriteHisto(gamma_current_answer_matrix_vs_gamma_prim_energy,factor,
>> 272 G4String("Adj_GammaCurrent_vs_EkinPrimGamma_Answer.txt"),
>> 273 G4String("Eprim1[MeV]\t\tEprim2[MeV]\t\tEsec1[MeV]\t\tEsec2[MeV]\t Current Efficiency[cm2*MeV*str]\terr[cm2*MeV*str]\n"));
>> 274
>> 275
>> 276
>> 277 WriteHisto(edep_answer_matrix_vs_proton_prim_energy,factor,
>> 278 G4String("Adj_Edep_vs_EkinPrimProton_Answer.txt"),
>> 279 G4String("E1[MeV]\t\tE2[MeV]\t\tEdep Efficiency[MeV*cm2*MeV*str]\terr_Edep[MeV*cm2*MeV*str]\n"));
>> 280
>> 281 WriteHisto(electron_current_answer_matrix_vs_proton_prim_energy,factor,
>> 282 G4String("Adj_ElectronCurrent_vs_EkinPrimProton_Answer.txt"),
>> 283 G4String("Eprim1[MeV]\t\tEprim2[MeV]\t\tEsec1[MeV]\t\tEsec2[MeV]\t Current Efficiency[cm2*MeV*str]\terr[cm2*MeV*str]\n"));
>> 284
>> 285 WriteHisto(gamma_current_answer_matrix_vs_proton_prim_energy,factor,
>> 286 G4String("Adj_GammaCurrent_vs_EkinPrimProton_Answer.txt"),
>> 287 G4String("Eprim1[MeV]\t\tEprim2[MeV]\t\tEsec1[MeV]\t\tEsec2[MeV]\t Current Efficiency[cm2*MeV*str]\terr[cm2*MeV*str]\n"));
>> 288
>> 289 WriteHisto(proton_current_answer_matrix_vs_proton_prim_energy,factor,
>> 290 G4String("Adj_ProtonCurrent_vs_EkinPrimProton_Answer.txt"),
>> 291 G4String("Eprim1[MeV]\t\tEprim2[MeV]\t\tEsec1[MeV]\t\tEsec2[MeV]\t Current Efficiency[cm2*MeV*str]\terr[cm2*MeV*str]\n"));
>> 292 }
>> 293 ConvergenceFileOutput.close();
>> 294 }
>> 295 ////////////////////////////////////////////////////////////////////////////////
>> 296 //
>> 297 void RMC01AnalysisManager::BeginOfEvent(const G4Event* )
>> 298 {;
175 } 299 }
176 << 300 ////////////////////////////////////////////////////////////////////////////////
177 //....oooOO0OOooo........oooOO0OOooo........oo << 301 //
178 <<
179 void RMC01AnalysisManager::BeginOfEvent(const <<
180 { <<
181 ; <<
182 } <<
183 <<
184 //....oooOO0OOooo........oooOO0OOooo........oo <<
185 <<
186 void RMC01AnalysisManager::EndOfEvent(const G4 302 void RMC01AnalysisManager::EndOfEvent(const G4Event* anEvent)
187 { << 303 {
188 if (fAdjoint_sim_mode) <<
189 EndOfEventForAdjointSimulation(anEvent); <<
190 else <<
191 EndOfEventForForwardSimulation(anEvent); <<
192 <<
193 // Test convergence. The error is already co <<
194 //-------------------------------------- <<
195 G4int nb_event = anEvent->GetEventID() + 1; <<
196 // G4double factor=1.; <<
197 if (fAdjoint_sim_mode) { <<
198 G4double n_adj_evt = nb_event / fNb_evt_pe <<
199 // nb_event/fNb_evt_per_adj_evt; <<
200 if (n_adj_evt * fNb_evt_per_adj_evt == nb_ <<
201 nb_event = static_cast<G4int>(n_adj_evt) <<
202 } <<
203 else <<
204 nb_event = 0; <<
205 } <<
206 <<
207 if (nb_event > 100 && fStop_run_if_precision <<
208 G4cout << fPrecision_to_reach * 100. << "% <<
209 fTimer->Stop(); <<
210 fElapsed_time += fTimer->GetRealElapsed(); <<
211 fConvergenceFileOutput << fMean_edep << '\ <<
212 << std::endl; <<
213 G4RunManager::GetRunManager()->AbortRun(tr <<
214 } <<
215 304
216 if (nb_event > 0 && nb_event % fNb_evt_modul <<
217 fTimer->Stop(); <<
218 fElapsed_time += fTimer->GetRealElapsed(); <<
219 fTimer->Start(); <<
220 fConvergenceFileOutput << fMean_edep << '\ <<
221 << std::endl; <<
222 } <<
223 } <<
224 <<
225 //....oooOO0OOooo........oooOO0OOooo........oo <<
226 <<
227 void RMC01AnalysisManager::EndOfEventForForwar <<
228 { <<
229 G4SDManager* SDman = G4SDManager::GetSDMpoin <<
230 G4HCofThisEvent* HCE = anEvent->GetHCofThisE <<
231 RMC01DoubleWithWeightHitsCollection* edepCol <<
232 (RMC01DoubleWithWeightHitsCollection*)(HCE <<
233 305
234 RMC01DoubleWithWeightHitsCollection* electro <<
235 (RMC01DoubleWithWeightHitsCollection*)(HCE <<
236 306
237 RMC01DoubleWithWeightHitsCollection* protonC << 307 if (adjoint_sim_mode) EndOfEventForAdjointSimulation(anEvent);
238 (RMC01DoubleWithWeightHitsCollection*)(HCE << 308 else EndOfEventForForwardSimulation(anEvent);
239 << 309
240 RMC01DoubleWithWeightHitsCollection* gammaCu << 310
241 (RMC01DoubleWithWeightHitsCollection*)(HCE << 311 //Test convergence. The error is already computed
242 << 312 //--------------------------------------
243 // Total energy deposited in Event << 313 G4int nb_event=anEvent->GetEventID()+1;
244 //------------------------------- << 314 //G4double factor=1.;
245 G4double totEdep = 0; << 315 if (adjoint_sim_mode) {
246 std::size_t i; << 316 G4double n_adj_evt= nb_event/nb_evt_per_adj_evt;
247 for (i = 0; i < edepCollection->entries(); + << 317 // nb_event/nb_evt_per_adj_evt;
248 totEdep += (*edepCollection)[i]->GetValue( << 318 if (n_adj_evt*nb_evt_per_adj_evt == nb_event) {
249 << 319 nb_event =static_cast<G4int>(n_adj_evt);
250 if (totEdep > 0.) { << 320 //factor=1.*G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
251 fAccumulated_edep += totEdep; << 321 }
252 fAccumulated_edep2 += totEdep * totEdep; << 322 else nb_event=0;
253 fNentry += 1.0; << 323
254 G4PrimaryParticle* thePrimary = anEvent->G << 324 }
255 G4double E0 = thePrimary->GetG4code()->Get << 325
256 G4double P = thePrimary->GetMomentum().mag << 326
257 G4double prim_ekin = std::sqrt(E0 * E0 + P << 327
258 fEdep_vs_prim_ekin->fill(prim_ekin, totEde << 328 if (nb_event>100 && stop_run_if_precision_reached && precision_to_reach >relative_error) {
259 } << 329 G4cout<<precision_to_reach*100.<<"% Precision reached!"<<std::endl;
260 ComputeMeanEdepAndError(anEvent, fMean_edep, << 330 theTimer->Stop();
261 if (fError_mean_edep > 0) fRelative_error = << 331 elapsed_time+=theTimer->GetRealElapsed();
262 << 332 ConvergenceFileOutput<<mean_edep<<'\t'<<error_mean_edep<<'\t'<<elapsed_time<<std::endl;
263 // Particle current on sensitive cylinder << 333 G4RunManager::GetRunManager()->AbortRun(true);
264 //------------------------------------- << 334 }
265 << 335
266 for (i = 0; i < electronCurrentCollection->e << 336
267 G4double ekin = (*electronCurrentCollectio << 337 if (nb_event>0 && nb_event % nb_evt_modulo_for_convergence_test == 0) {
268 G4double weight = (*electronCurrentCollect << 338 theTimer->Stop();
269 fElectron_current->fill(ekin, weight); << 339 elapsed_time+=theTimer->GetRealElapsed();
270 } << 340 theTimer->Start();
271 << 341 ConvergenceFileOutput<<mean_edep<<'\t'<<error_mean_edep<<'\t'<<elapsed_time<<std::endl;
272 for (i = 0; i < protonCurrentCollection->ent << 342
273 G4double ekin = (*protonCurrentCollection) << 343 }
274 G4double weight = (*protonCurrentCollectio << 344
275 fProton_current->fill(ekin, weight); << 345
276 } << 346
277 << 347
278 for (i = 0; i < gammaCurrentCollection->entr << 348
279 G4double ekin = (*gammaCurrentCollection)[ << 349
280 G4double weight = (*gammaCurrentCollection << 350 }
281 fGamma_current->fill(ekin, weight); << 351 ////////////////////////////////////////////////////////////////////////////////
282 } << 352 //
283 } << 353 void RMC01AnalysisManager::EndOfEventForForwardSimulation(const G4Event* anEvent)
284 << 354 {
285 //....oooOO0OOooo........oooOO0OOooo........oo << 355
286 << 356 G4SDManager* SDman = G4SDManager::GetSDMpointer();
287 void RMC01AnalysisManager::EndOfEventForAdjoin << 357 G4HCofThisEvent* HCE = anEvent->GetHCofThisEvent();
288 { << 358 RMC01DoubleWithWeightHitsCollection* EdepCollection = (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("edep")));
289 // Output from Sensitive volume computed dur << 359 RMC01DoubleWithWeightHitsCollection* ElectronCurrentCollection = (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("current_electron")));
>> 360 RMC01DoubleWithWeightHitsCollection* ProtonCurrentCollection = (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("current_proton")));
>> 361 RMC01DoubleWithWeightHitsCollection* GammaCurrentCollection = (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("current_gamma")));
>> 362
>> 363 //Total energy deposited in Event
>> 364 //-------------------------------
>> 365 G4double totEdep=0;
>> 366 G4int i;
>> 367 for (i=0;i<EdepCollection->entries();i++) totEdep+=(*EdepCollection)[i]->GetValue()*(*EdepCollection)[i]->GetWeight();
>> 368
>> 369 if (totEdep>0.){
>> 370 accumulated_edep +=totEdep ;
>> 371 accumulated_edep2 +=totEdep*totEdep;
>> 372 G4PrimaryParticle* thePrimary=anEvent->GetPrimaryVertex()->GetPrimary();
>> 373 G4double E0= thePrimary->GetG4code()->GetPDGMass();
>> 374 G4double P=thePrimary->GetMomentum().mag();
>> 375 G4double prim_ekin =std::sqrt(E0*E0+P*P)-E0;
>> 376 edep_vs_prim_ekin->fill(prim_ekin,totEdep);
>> 377 }
>> 378 ComputeMeanEdepAndError(anEvent,mean_edep,error_mean_edep);
>> 379 if (error_mean_edep>0) relative_error= error_mean_edep/mean_edep;
>> 380
>> 381
>> 382
>> 383 //Particle current on sensitive cylinder
>> 384 //-------------------------------------
>> 385
>> 386 for (i=0;i<ElectronCurrentCollection->entries();i++) {
>> 387 G4double ekin =(*ElectronCurrentCollection)[i]->GetValue();
>> 388 G4double weight=(*ElectronCurrentCollection)[i]->GetWeight();
>> 389 electron_current->fill(ekin,weight);
>> 390 }
>> 391
>> 392 for (i=0;i<ProtonCurrentCollection->entries();i++) {
>> 393 G4double ekin =(*ProtonCurrentCollection)[i]->GetValue();
>> 394 G4double weight=(*ProtonCurrentCollection)[i]->GetWeight();
>> 395 proton_current->fill(ekin,weight);
>> 396 }
>> 397
>> 398 for (i=0;i<GammaCurrentCollection->entries();i++) {
>> 399 G4double ekin =(*GammaCurrentCollection)[i]->GetValue();
>> 400 G4double weight=(*GammaCurrentCollection)[i]->GetWeight();
>> 401 gamma_current->fill(ekin,weight);
>> 402 }
>> 403
>> 404
>> 405 }
>> 406 ////////////////////////////////////////////////////////////////////////////////
>> 407 //
>> 408 void RMC01AnalysisManager::EndOfEventForAdjointSimulation(const G4Event* anEvent)
>> 409 {
>> 410 //Output from Sensitive volume computed during the forward tracking phase
290 //------------------------------------------ 411 //-----------------------------------------------------------------------
>> 412
291 G4SDManager* SDman = G4SDManager::GetSDMpoin 413 G4SDManager* SDman = G4SDManager::GetSDMpointer();
292 G4HCofThisEvent* HCE = anEvent->GetHCofThisE 414 G4HCofThisEvent* HCE = anEvent->GetHCofThisEvent();
293 RMC01DoubleWithWeightHitsCollection* edepCol << 415 RMC01DoubleWithWeightHitsCollection* EdepCollection = (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("edep")));
294 (RMC01DoubleWithWeightHitsCollection*)(HCE << 416 RMC01DoubleWithWeightHitsCollection* ElectronCurrentCollection = (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("current_electron")));
295 << 417 RMC01DoubleWithWeightHitsCollection* ProtonCurrentCollection = (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("current_proton")));
296 RMC01DoubleWithWeightHitsCollection* electro << 418 RMC01DoubleWithWeightHitsCollection* GammaCurrentCollection = (RMC01DoubleWithWeightHitsCollection*)(HCE->GetHC(SDman->GetCollectionID("current_gamma")));
297 (RMC01DoubleWithWeightHitsCollection*)(HCE << 419
298 << 420 //Output from adjoint tracking phase
299 RMC01DoubleWithWeightHitsCollection* protonC <<
300 (RMC01DoubleWithWeightHitsCollection*)(HCE <<
301 <<
302 RMC01DoubleWithWeightHitsCollection* gammaCu <<
303 (RMC01DoubleWithWeightHitsCollection*)(HCE <<
304 <<
305 // Computation of total energy deposited in <<
306 //------------------------------- <<
307 G4double totEdep = 0; <<
308 std::size_t i; <<
309 for (i = 0; i < edepCollection->entries(); + <<
310 totEdep += (*edepCollection)[i]->GetValue( <<
311 <<
312 // Output from adjoint tracking phase <<
313 //------------------------------------------ 421 //----------------------------------------------------------------------------
314 << 422
315 G4AdjointSimManager* theAdjointSimManager = 423 G4AdjointSimManager* theAdjointSimManager = G4AdjointSimManager::GetInstance();
316 << 424 G4int pdg_nb =theAdjointSimManager->GetFwdParticlePDGEncodingAtEndOfLastAdjointTrack();
317 size_t nb_adj_track = theAdjointSimManager-> << 425 G4double prim_ekin=theAdjointSimManager->GetEkinAtEndOfLastAdjointTrack();
318 G4double total_normalised_weight = 0.; << 426 G4double adj_weight=theAdjointSimManager->GetWeightAtEndOfLastAdjointTrack();
319 << 427
320 // We need to loop over the adjoint tracks t << 428
321 // surface. << 429 //Factor of normalisation to user selected primary spectrum (power law or exponential)
322 for (std::size_t j = 0; j < nb_adj_track; ++ << 430 //------------------------------------------------------------------------------------
323 G4int pdg_nb = theAdjointSimManager->GetFw << 431
324 G4double prim_ekin = theAdjointSimManager- << 432 G4double normalised_weight = 0.;
325 G4double adj_weight = theAdjointSimManager << 433 if (pdg_nb== thePrimPDG_ID && prim_ekin>= emin_prim_spectrum && prim_ekin<= emax_prim_spectrum)
326 << 434 normalised_weight = adj_weight*PrimDiffAndDirectionalFluxForAdjointSim(prim_ekin);
327 // Factor of normalisation to user defined << 435
328 //---------------------------------------- << 436
329 G4double normalised_weight = 0.; << 437 //Answer matrices
330 if (pdg_nb == fPrimPDG_ID && prim_ekin >= << 438 //-------------
331 && prim_ekin <= fEmax_prim_spectrum) << 439 Histo1DVar* edep_answer_matrix =0;
332 normalised_weight = adj_weight * PrimDif << 440 Histo2DVar* electron_current_answer_matrix =0;
333 total_normalised_weight += normalised_weig << 441 Histo2DVar* gamma_current_answer_matrix =0;
334 << 442 Histo2DVar* proton_current_answer_matrix =0;
335 // Answer matrices << 443
336 //------------- << 444
337 G4H1* edep_rmatrix = 0; << 445 if (pdg_nb == G4Electron::Electron()->GetPDGEncoding()){ //electron answer matrices
338 G4H2* electron_current_rmatrix = 0; << 446 edep_answer_matrix = edep_answer_matrix_vs_electron_prim_energy;
339 G4H2* gamma_current_rmatrix = 0; << 447 electron_current_answer_matrix = electron_current_answer_matrix_vs_electron_prim_energy;
340 G4H2* proton_current_rmatrix = 0; << 448 gamma_current_answer_matrix = gamma_current_answer_matrix_vs_electron_prim_energy;
341 << 449 }
342 if (pdg_nb == G4Electron::Electron()->GetP << 450 else if (pdg_nb == G4Gamma::Gamma()->GetPDGEncoding()){ //electron answer matrices
343 edep_rmatrix = fEdep_rmatrix_vs_electron << 451 edep_answer_matrix = edep_answer_matrix_vs_gamma_prim_energy;
344 electron_current_rmatrix = fElectron_cur << 452 electron_current_answer_matrix = electron_current_answer_matrix_vs_gamma_prim_energy;
345 gamma_current_rmatrix = fGamma_current_r << 453 gamma_current_answer_matrix = gamma_current_answer_matrix_vs_gamma_prim_energy;
346 } << 454 }
347 else if (pdg_nb == G4Gamma::Gamma()->GetPD << 455 else if (pdg_nb == G4Proton::Proton()->GetPDGEncoding()){ //electron answer matrices
348 // gammma answer matrices << 456 edep_answer_matrix = edep_answer_matrix_vs_proton_prim_energy;
349 edep_rmatrix = fEdep_rmatrix_vs_gamma_pr << 457 electron_current_answer_matrix = electron_current_answer_matrix_vs_proton_prim_energy;
350 electron_current_rmatrix = fElectron_cur << 458 gamma_current_answer_matrix = gamma_current_answer_matrix_vs_proton_prim_energy;
351 gamma_current_rmatrix = fGamma_current_r << 459 proton_current_answer_matrix = proton_current_answer_matrix_vs_proton_prim_energy;
352 } << 460 }
353 else if (pdg_nb == G4Proton::Proton()->Get << 461
354 // proton answer matrices << 462
355 edep_rmatrix = fEdep_rmatrix_vs_proton_p << 463
356 electron_current_rmatrix = fElectron_cur << 464 //Registering of total energy deposited in Event
357 gamma_current_rmatrix = fGamma_current_r <<
358 proton_current_rmatrix = fProton_current <<
359 } <<
360 // Register histo edep vs prim ekin <<
361 //---------------------------------- <<
362 if (normalised_weight > 0) fEdep_vs_prim_e <<
363 // Registering answer matrix <<
364 //--------------------------- <<
365 edep_rmatrix->fill(prim_ekin, totEdep * ad <<
366 <<
367 // Registering of current of particles on <<
368 //---------------------------------------- <<
369 <<
370 for (i = 0; i < electronCurrentCollection- <<
371 G4double ekin = (*electronCurrentCollect <<
372 G4double weight = (*electronCurrentColle <<
373 fElectron_current->fill(ekin, weight * n <<
374 electron_current_rmatrix->fill(prim_ekin <<
375 } <<
376 for (i = 0; i < protonCurrentCollection->e <<
377 G4double ekin = (*protonCurrentCollectio <<
378 G4double weight = (*protonCurrentCollect <<
379 fProton_current->fill(ekin, weight * nor <<
380 proton_current_rmatrix->fill(prim_ekin, <<
381 } <<
382 for (i = 0; i < gammaCurrentCollection->en <<
383 G4double ekin = (*gammaCurrentCollection <<
384 G4double weight = (*gammaCurrentCollecti <<
385 fGamma_current->fill(ekin, weight * norm <<
386 gamma_current_rmatrix->fill(prim_ekin, e <<
387 } <<
388 } <<
389 <<
390 // Registering of total energy deposited in <<
391 //------------------------------- 465 //-------------------------------
392 G4bool new_mean_computed = false; << 466 G4double totEdep=0;
393 if (totEdep > 0.) { << 467 G4int i;
394 if (total_normalised_weight > 0.) { << 468 for (i=0;i<EdepCollection->entries();i++) totEdep+=(*EdepCollection)[i]->GetValue()*(*EdepCollection)[i]->GetWeight();
395 G4double edep = totEdep * total_normalis << 469
396 << 470 G4bool new_mean_computed=false;
397 // Check if the edep is not wrongly too << 471 if (totEdep>0.){
398 //-------------------------------------- << 472 if (normalised_weight>0.){
399 G4double new_mean(0.0), new_error(0.0); << 473 G4double edep=totEdep* normalised_weight;
400 fAccumulated_edep += edep; << 474
401 fAccumulated_edep2 += edep * edep; << 475 //Check if the edep is not wrongly too high
402 fNentry += 1.0; << 476 //-----------------------------------------
403 ComputeMeanEdepAndError(anEvent, new_mea << 477 G4double new_mean , new_error;
404 G4double new_relative_error = 1.; << 478
405 if (new_error > 0) new_relative_error = << 479 accumulated_edep +=edep;
406 if (fRelative_error < 0.10 && new_relati << 480 accumulated_edep2 +=edep*edep;
407 G4cout << "Potential wrong adjoint wei << 481
408 G4cout << "The results of this event w << 482 ComputeMeanEdepAndError(anEvent,new_mean,new_error);
409 G4cout << "previous mean edep [MeV] " << 483 G4double new_relative_error = 1.;
410 G4cout << "previous relative error " < << 484 if ( new_error >0) new_relative_error = new_error/ new_mean;
411 G4cout << "new rejected mean edep [MeV << 485 if (relative_error <0.10 && new_relative_error>1.5*relative_error) { // rejected
412 G4cout << "new rejected relative error << 486 G4cout<<"Potential wrong adjoint weight!"<<std::endl;
413 fAccumulated_edep -= edep; << 487 G4cout<<"The results of this event will not be registered!"<<std::endl;
414 fAccumulated_edep2 -= edep * edep; << 488 G4cout<<"previous mean edep [MeV] "<< mean_edep<<std::endl;
415 fNentry -= 1.0; << 489 G4cout<<"previous relative error "<< relative_error<<std::endl;
416 return; << 490 G4cout<<"new rejected mean edep [MeV] "<< new_mean<<std::endl;
417 } << 491 G4cout<<"new rejected relative error "<< new_relative_error<<std::endl;
418 else { // accepted << 492 accumulated_edep -=edep;
419 fMean_edep = new_mean; << 493 accumulated_edep2 -=edep*edep;
420 fError_mean_edep = new_error; << 494 return;
421 fRelative_error = new_relative_error; << 495 }
422 new_mean_computed = true; << 496 else { //accepted
423 } << 497 mean_edep = new_mean;
424 } << 498 error_mean_edep = new_error;
425 << 499 relative_error =new_relative_error;
426 if (!new_mean_computed) { << 500 new_mean_computed=true;
427 ComputeMeanEdepAndError(anEvent, fMean_e << 501 }
428 fRelative_error = (fMean_edep > 0.0) ? f << 502 edep_vs_prim_ekin->fill(prim_ekin,edep);
429 } << 503 }
430 } << 504
>> 505 // Registering answer matrix
>> 506 //---------------------------
>> 507
>> 508 edep_answer_matrix->fill(prim_ekin,totEdep*adj_weight/cm2);
>> 509
>> 510
>> 511 }
>> 512 if (!new_mean_computed){
>> 513 ComputeMeanEdepAndError(anEvent,mean_edep,error_mean_edep);
>> 514 if (error_mean_edep>0) relative_error= error_mean_edep/mean_edep;
>> 515 }
>> 516
>> 517
>> 518 //Registering of current of particles on the sensitive volume
>> 519 //------------------------------------------------------------
>> 520
>> 521 for (i=0;i<ElectronCurrentCollection->entries();i++) {
>> 522 G4double ekin =(*ElectronCurrentCollection)[i]->GetValue();
>> 523 G4double weight=(*ElectronCurrentCollection)[i]->GetWeight();
>> 524 electron_current->fill(ekin,weight*normalised_weight);
>> 525 electron_current_answer_matrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
>> 526 }
>> 527
>> 528 for (i=0;i<ProtonCurrentCollection->entries();i++) {
>> 529 G4double ekin =(*ProtonCurrentCollection)[i]->GetValue();
>> 530 G4double weight=(*ProtonCurrentCollection)[i]->GetWeight();
>> 531 proton_current->fill(ekin,weight*normalised_weight);
>> 532 proton_current_answer_matrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
>> 533 }
>> 534
>> 535 for (i=0;i<GammaCurrentCollection->entries();i++) {
>> 536 G4double ekin =(*GammaCurrentCollection)[i]->GetValue();
>> 537 G4double weight=(*GammaCurrentCollection)[i]->GetWeight();
>> 538 gamma_current->fill(ekin,weight*normalised_weight);
>> 539 gamma_current_answer_matrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
>> 540
>> 541 }
>> 542
>> 543
431 } 544 }
432 << 545 ////////////////////////////////////////////////////////////////////////////////
433 //....oooOO0OOooo........oooOO0OOooo........oo << 546 //
434 << 547 G4double RMC01AnalysisManager::PrimDiffAndDirectionalFluxForAdjointSim(G4double prim_energy)
435 G4double RMC01AnalysisManager::PrimDiffAndDirF << 548 {
436 { << 549 G4double flux=amplitude_prim_spectrum;
437 G4double flux = fAmplitude_prim_spectrum; << 550 if ( the_prim_spectrum_type ==EXPO) flux*=std::exp(-prim_energy/alpha_or_E0);
438 if (fPrimSpectrumType == EXPO) << 551 else flux*=std::pow(prim_energy, -alpha_or_E0);
439 flux *= std::exp(-prim_energy / fAlpha_or_ <<
440 else <<
441 flux *= std::pow(prim_energy, -fAlpha_or_E <<
442 return flux; 552 return flux;
443 } 553 }
444 << 554 ////////////////////////////////////////////////////////////////////////////////
445 //....oooOO0OOooo........oooOO0OOooo........oo << 555 //
446 /* << 556 void RMC01AnalysisManager::WriteHisto(Histo1DVar* anHisto, G4double scaling_factor, G4String fileName, G4String header_lines)
447 void RMC01AnalysisManager::WriteHisto(G4H1* a <<
448 G4double scaling_factor, G4String <<
449 { std::fstream FileOutput(fileName, std::ios:: 557 { std::fstream FileOutput(fileName, std::ios::out);
450 FileOutput<<header_lines; 558 FileOutput<<header_lines;
>> 559
451 FileOutput.setf(std::ios::scientific); 560 FileOutput.setf(std::ios::scientific);
452 FileOutput.precision(6); 561 FileOutput.precision(6);
453 562
454 for (G4int i =0;i<G4int(anHisto->axis().bins << 563 size_t nxbins = anHisto->part.total_bins();
455 FileOutput<<anHisto->axis().bin_lower_ << 564 for (size_t i =0;i<nxbins;i++) {
456 <<'\t'<<anHisto->axis().bin_uppe << 565 FileOutput<<anHisto->part.get_bin_position(i)
457 <<'\t'<<anHisto->bin_height(i)*s << 566 <<'\t'<<anHisto->part.get_bin_position(i+1)
458 <<'\t'<<anHisto->bin_error(i)*sc << 567 <<'\t'<<anHisto->get_bin_value(int(i))*scaling_factor
>> 568 <<'\t'<<anHisto->get_bin_error(int(i))*scaling_factor<<std::endl;
459 } 569 }
460 } 570 }
461 << 571 ////////////////////////////////////////////////////////////////////////////////
462 //....oooOO0OOooo........oooOO0OOooo........oo << 572 //
463 << 573 void RMC01AnalysisManager::WriteHisto(Histo2DVar* anHisto, G4double scaling_factor, G4String fileName, G4String header_lines)
464 void RMC01AnalysisManager::WriteHisto(G4H2* a <<
465 G4double scaling_factor, G4String <<
466 { std::fstream FileOutput(fileName, std::ios:: 574 { std::fstream FileOutput(fileName, std::ios::out);
467 FileOutput<<header_lines; 575 FileOutput<<header_lines;
468 << 576
469 FileOutput.setf(std::ios::scientific); 577 FileOutput.setf(std::ios::scientific);
470 FileOutput.precision(6); 578 FileOutput.precision(6);
471 << 579
472 for (G4int i =0;i<G4int(anHisto->axis_x().bi << 580 size_t nxbins = anHisto->x_part.total_bins();
473 for (G4int j =0;j<G4int(anHisto->axis_y(). << 581 size_t nybins = anHisto->y_part.total_bins();
474 FileOutput<<anHisto->axis_x().bin_lower << 582 for (size_t i =0;i<nxbins;i++) {
475 <<'\t'<<anHisto->axis_x() << 583 for (size_t j =0;j<nybins;j++) {
476 <<'\t'<<anHisto->axis_y().b << 584 FileOutput<<anHisto->x_part.get_bin_position(i)
477 <<'\t'<<anHisto->axis_y().b << 585 <<'\t'<<anHisto->x_part.get_bin_position(i+1)
478 <<'\t'<<anHisto->bin_height << 586 <<'\t'<<anHisto->y_part.get_bin_position(j)
479 <<'\t'<<anHisto->bin_error(i,j << 587 <<'\t'<<anHisto->y_part.get_bin_position(j+1)
480 << 588 <<'\t'<<anHisto->get_bin_value(int(i),int(j))*scaling_factor
481 } << 589 <<'\t'<<anHisto->get_bin_error(int(i),int(j))*scaling_factor<<std::endl;
>> 590 }
482 } 591 }
483 } 592 }
484 */ << 593 ////////////////////////////////////////////////////////////////////////////////
485 //....oooOO0OOooo........oooOO0OOooo........oo << 594 //
486 << 595 void RMC01AnalysisManager::ResetHistograms()
487 void RMC01AnalysisManager::ComputeMeanEdepAndE << 596 { edep_vs_prim_ekin->reset();
488 << 597 electron_current->reset();
489 { << 598 proton_current->reset();
490 G4int nb_event = anEvent->GetEventID() + 1; << 599 gamma_current->reset();
491 G4double factor = 1.; << 600
492 if (fAdjoint_sim_mode) { << 601 edep_answer_matrix_vs_electron_prim_energy->reset();
493 nb_event /= fNb_evt_per_adj_evt; << 602 electron_current_answer_matrix_vs_electron_prim_energy->reset();
494 factor = 1. * G4AdjointSimManager::GetInst << 603 gamma_current_answer_matrix_vs_electron_prim_energy->reset();
495 } << 604
496 << 605 edep_answer_matrix_vs_gamma_prim_energy->reset();
497 // VI: error computation now is based on num << 606 electron_current_answer_matrix_vs_gamma_prim_energy->reset();
498 // number of events << 607 gamma_current_answer_matrix_vs_gamma_prim_energy->reset();
499 // LD: This is wrong! With the use of fNentr << 608
500 // correctly normalised. The mean and th << 609 edep_answer_matrix_vs_proton_prim_energy->reset();
501 // with nb_event. The old computation ha << 610 electron_current_answer_matrix_vs_proton_prim_energy->reset();
502 // VI: OK, but let computations be double << 611 proton_current_answer_matrix_vs_proton_prim_energy->reset();
503 if (nb_event > 0) { << 612 gamma_current_answer_matrix_vs_proton_prim_energy->reset();
504 G4double norm = 1.0 / (G4double)nb_event; <<
505 mean = fAccumulated_edep * norm; <<
506 G4double mean_x2 = fAccumulated_edep2 * no <<
507 G4double zz = mean_x2 - mean * mean; <<
508 /* <<
509 G4cout << "Nevt= " << nb_event << " mean <<
510 << " mean_x2= " << mean_x2 << " <<
511 << zz << G4endl; <<
512 */ <<
513 error = factor * std::sqrt(std::max(zz, 0. <<
514 mean *= factor; <<
515 } <<
516 else { <<
517 mean = 0; <<
518 error = 0; <<
519 } <<
520 // G4cout << "Aend: " << mean << " " << erro <<
521 } 613 }
522 << 614 ////////////////////////////////////////////////////////////////////////////////
523 //....oooOO0OOooo........oooOO0OOooo........oo << 615 //
524 << 616 void RMC01AnalysisManager::ComputeMeanEdepAndError(const G4Event* anEvent,G4double& mean,G4double& error)
525 void RMC01AnalysisManager::SetPrimaryExpSpectr << 617 {
526 << 618 G4int nb_event=anEvent->GetEventID()+1;
527 << 619 G4double factor=1.;
528 { << 620 if (adjoint_sim_mode) {
529 fPrimSpectrumType = EXPO; << 621 nb_event /=nb_evt_per_adj_evt;
530 if (particle_name == "e-") << 622 factor=1.*G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
531 fPrimPDG_ID = G4Electron::Electron()->GetP << 623 }
532 else if (particle_name == "gamma") << 624
533 fPrimPDG_ID = G4Gamma::Gamma()->GetPDGEnco << 625 //error computation
534 else if (particle_name == "proton") << 626 if (nb_event>1) {
535 fPrimPDG_ID = G4Proton::Proton()->GetPDGEn << 627 mean = accumulated_edep/nb_event;
536 else { << 628 G4double mean_x2 =accumulated_edep2/nb_event;
537 G4cout << "The particle that you did selec << 629 error = factor*std::sqrt(mean_x2-mean*mean)/std::sqrt(G4double(nb_event));
538 << "list for primary [e-, gamma, pr << 630 mean *=factor;
539 return; << 631 } else {
>> 632 mean=0;
>> 633 error=0;
540 } 634 }
541 fAlpha_or_E0 = E0; <<
542 fAmplitude_prim_spectrum = <<
543 omni_fluence / E0 / (std::exp(-Emin / E0) <<
544 fEmin_prim_spectrum = Emin; <<
545 fEmax_prim_spectrum = Emax; <<
546 } 635 }
547 << 636 ////////////////////////////////////////////////////////////////////////////////
548 //....oooOO0OOooo........oooOO0OOooo........oo << 637 //
549 << 638 void RMC01AnalysisManager::SetPrimaryExponentialSpectrumForAdjointSim(const G4String& particle_name, G4double omni_fluence,G4double E0, G4double Emin,G4double Emax)
550 void RMC01AnalysisManager::SetPrimaryPowerLawS << 639 { the_prim_spectrum_type = EXPO;
551 << 640 if (particle_name == "e-" ) thePrimPDG_ID = G4Electron::Electron()->GetPDGEncoding();
552 << 641 else if (particle_name == "gamma") thePrimPDG_ID = G4Gamma::Gamma()->GetPDGEncoding();
553 << 642 else if (particle_name == "proton") thePrimPDG_ID = G4Proton::Proton()->GetPDGEncoding();
554 { <<
555 fPrimSpectrumType = POWER; <<
556 if (particle_name == "e-") <<
557 fPrimPDG_ID = G4Electron::Electron()->GetP <<
558 else if (particle_name == "gamma") <<
559 fPrimPDG_ID = G4Gamma::Gamma()->GetPDGEnco <<
560 else if (particle_name == "proton") <<
561 fPrimPDG_ID = G4Proton::Proton()->GetPDGEn <<
562 else { <<
563 G4cout << "The particle that you did selec <<
564 << " list for primary [e-, gamma, p <<
565 return; <<
566 } <<
567 <<
568 if (alpha == 1.) { <<
569 fAmplitude_prim_spectrum = omni_fluence / <<
570 } <<
571 else { 643 else {
572 G4double p = 1. - alpha; << 644 G4cout<<"The particle that you did select is not in the candidate list for primary [e-, gamma, proton]!"<<G4endl;
573 fAmplitude_prim_spectrum = omni_fluence / << 645 return;
574 } << 646 }
575 << 647
576 fAlpha_or_E0 = alpha; << 648
577 fEmin_prim_spectrum = Emin; << 649 alpha_or_E0 = E0 ;
578 fEmax_prim_spectrum = Emax; << 650 amplitude_prim_spectrum = omni_fluence/E0/(std::exp(-Emin/E0)-std::exp(-Emax/E0))/4./pi;
>> 651 emin_prim_spectrum = Emin ;
>> 652 emax_prim_spectrum = Emax;
579 } 653 }
580 << 654 ////////////////////////////////////////////////////////////////////////////////
581 //....oooOO0OOooo........oooOO0OOooo........oo << 655 //
582 << 656 void RMC01AnalysisManager::SetPrimaryPowerLawSpectrumForAdjointSim(const G4String& particle_name, G4double omni_fluence,G4double alpha, G4double Emin,G4double Emax)
583 void RMC01AnalysisManager::Book() << 657 { the_prim_spectrum_type =POWER;
584 { << 658 if (particle_name == "e-" ) thePrimPDG_ID = G4Electron::Electron()->GetPDGEncoding();
585 //---------------------- << 659 else if (particle_name == "gamma") thePrimPDG_ID = G4Gamma::Gamma()->GetPDGEncoding();
586 // Creation of histograms << 660 else if (particle_name == "proton") thePrimPDG_ID = G4Proton::Proton()->GetPDGEncoding();
587 //---------------------- << 661 else {
588 << 662 G4cout<<"The particle that you did select is not in the candidate list for primary [e-, gamma, proton]!"<<G4endl;
589 // Energy binning of the histograms : 60 log << 663 return;
590 << 664 }
591 G4double emin = 1. * keV; << 665
592 G4double emax = 1. * GeV; << 666
593 << 667 if (alpha ==1.) {
594 // file_name << 668 amplitude_prim_spectrum = omni_fluence/std::log(Emax/Emin)/4./pi;
595 fFileName[0] = "forward_sim"; << 669 }
596 if (fAdjoint_sim_mode) fFileName[0] = "adjoi << 670 else {
597 << 671 G4double p=1.-alpha;
598 // Histo manager << 672 amplitude_prim_spectrum = omni_fluence/p/(std::pow(Emax,p)-std::pow(Emin,p))/4./pi;
599 G4AnalysisManager* theHistoManager = G4Analy << 673
600 theHistoManager->SetDefaultFileType("root"); << 674 }
601 G4String extension = theHistoManager->GetFil << 675
602 fFileName[1] = fFileName[0] + "." + extensio << 676 alpha_or_E0 = alpha ;
603 theHistoManager->SetFirstHistoId(1); << 677 emin_prim_spectrum = Emin ;
604 << 678 emax_prim_spectrum = Emax;
605 G4bool fileOpen = theHistoManager->OpenFile( <<
606 if (!fileOpen) { <<
607 G4cout << "\n---> RMC01AnalysisManager::Bo <<
608 return; <<
609 } <<
610 <<
611 // Create directories <<
612 theHistoManager->SetHistoDirectoryName("hist <<
613 <<
614 // Histograms for : <<
615 // 1)the forward simulation results <<
616 // 2)the Reverse MC simulation resu <<
617 //------------------------------------------ <<
618 <<
619 G4int idHisto = <<
620 theHistoManager->CreateH1(G4String("Edep_v <<
621 60, emin, emax, <<
622 fEdep_vs_prim_ekin = theHistoManager->GetH1( <<
623 <<
624 idHisto = theHistoManager->CreateH1(G4String <<
625 emax, "n <<
626 <<
627 fElectron_current = theHistoManager->GetH1(i <<
628 <<
629 idHisto = theHistoManager->CreateH1(G4String <<
630 emax, "n <<
631 fProton_current = theHistoManager->GetH1(idH <<
632 <<
633 idHisto = theHistoManager->CreateH1(G4String <<
634 "none", <<
635 fGamma_current = theHistoManager->GetH1(idHi <<
636 <<
637 // Response matrices for the adjoint simulat <<
638 //------------------------------------------ <<
639 if (fAdjoint_sim_mode) { <<
640 // Response matrices for external isotropi <<
641 //---------------------------------------- <<
642 <<
643 idHisto = theHistoManager->CreateH1(G4Stri <<
644 G4Stri <<
645 emax, <<
646 fEdep_rmatrix_vs_electron_prim_energy = th <<
647 <<
648 idHisto = theHistoManager->CreateH2( <<
649 G4String("Electron_current_rmatrix_vs_el <<
650 G4String("electron current RM vs e- pri <<
651 "none", "none", "none", G4String("log"), <<
652 <<
653 fElectron_current_rmatrix_vs_electron_prim <<
654 <<
655 idHisto = theHistoManager->CreateH2(G4Stri <<
656 G4Stri <<
657 emin, <<
658 G4Stri <<
659 <<
660 fGamma_current_rmatrix_vs_electron_prim_en <<
661 <<
662 // Response matrices for external isotropi <<
663 <<
664 idHisto = theHistoManager->CreateH1(G4Stri <<
665 G4Stri <<
666 emax, <<
667 fEdep_rmatrix_vs_gamma_prim_energy = theHi <<
668 <<
669 idHisto = theHistoManager->CreateH2(G4Stri <<
670 G4Stri <<
671 60, em <<
672 "none" <<
673 <<
674 fElectron_current_rmatrix_vs_gamma_prim_en <<
675 <<
676 idHisto = theHistoManager->CreateH2(G4Stri <<
677 G4Stri <<
678 emin, <<
679 G4Stri <<
680 <<
681 fGamma_current_rmatrix_vs_gamma_prim_energ <<
682 <<
683 // Response matrices for external isotropi <<
684 idHisto = theHistoManager->CreateH1(G4Stri <<
685 G4Stri <<
686 emax, <<
687 fEdep_rmatrix_vs_proton_prim_energy = theH <<
688 <<
689 idHisto = theHistoManager->CreateH2(G4Stri <<
690 G4Stri <<
691 60, em <<
692 "none" <<
693 <<
694 fElectron_current_rmatrix_vs_proton_prim_e <<
695 <<
696 idHisto = theHistoManager->CreateH2(G4Stri <<
697 G4Stri <<
698 emin, <<
699 G4Stri <<
700 <<
701 fGamma_current_rmatrix_vs_proton_prim_ener <<
702 <<
703 idHisto = theHistoManager->CreateH2(G4Stri <<
704 G4Stri <<
705 emin, <<
706 G4Stri <<
707 <<
708 fProton_current_rmatrix_vs_proton_prim_ene <<
709 } <<
710 fFactoryOn = true; <<
711 G4cout << "\n----> Histogram Tree is opened <<
712 } <<
713 <<
714 //....oooOO0OOooo........oooOO0OOooo........oo <<
715 <<
716 void RMC01AnalysisManager::Save(G4double scali <<
717 { <<
718 if (fFactoryOn) { <<
719 G4AnalysisManager* theHistoManager = G4Ana <<
720 // scaling of results <<
721 //----------------- <<
722 <<
723 for (G4int ind = 1; ind <= theHistoManager <<
724 theHistoManager->SetH1Ascii(ind, true); <<
725 theHistoManager->ScaleH1(ind, scaling_fa <<
726 } <<
727 for (G4int ind = 1; ind <= theHistoManager <<
728 theHistoManager->ScaleH2(ind, scaling_fa <<
729 <<
730 theHistoManager->Write(); <<
731 theHistoManager->CloseFile(); <<
732 G4cout << "\n----> Histogram Tree is saved <<
733 <<
734 theHistoManager->Clear(); <<
735 fFactoryOn = false; <<
736 } <<
737 } 679 }
738 680