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