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