| 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 G4int 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 G4int 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 (size_t j=0;j<nb_adj_track;j++) {
323 G4int pdg_nb = theAdjointSimManager->GetFw << 333 G4int pdg_nb =theAdjointSimManager
324 G4double prim_ekin = theAdjointSimManager- << 334 ->GetFwdParticlePDGEncodingAtEndOfLastAdjointTrack(j);
325 G4double adj_weight = theAdjointSimManager << 335 G4double prim_ekin=theAdjointSimManager
326 << 336 ->GetEkinAtEndOfLastAdjointTrack(j);
327 // Factor of normalisation to user defined << 337 G4double adj_weight=theAdjointSimManager
>> 338 ->GetWeightAtEndOfLastAdjointTrack(j);
>> 339
>> 340
>> 341 //Factor of normalisation to user defined prim spectrum (power law or exp)
328 //---------------------------------------- 342 //------------------------------------------------------------------------
329 G4double normalised_weight = 0.; 343 G4double normalised_weight = 0.;
330 if (pdg_nb == fPrimPDG_ID && prim_ekin >= << 344 if (pdg_nb== fPrimPDG_ID && prim_ekin>= fEmin_prim_spectrum
331 && prim_ekin <= fEmax_prim_spectrum) << 345 && prim_ekin<= fEmax_prim_spectrum)
332 normalised_weight = adj_weight * PrimDif << 346 normalised_weight =
>> 347 adj_weight*PrimDiffAndDirFluxForAdjointSim(prim_ekin);
333 total_normalised_weight += normalised_weig 348 total_normalised_weight += normalised_weight;
334 << 349
335 // Answer matrices << 350 //Answer matrices
336 //------------- 351 //-------------
337 G4H1* edep_rmatrix = 0; << 352 G4AnaH1* edep_rmatrix =0;
338 G4H2* electron_current_rmatrix = 0; << 353 G4AnaH2* electron_current_rmatrix =0;
339 G4H2* gamma_current_rmatrix = 0; << 354 G4AnaH2* gamma_current_rmatrix =0;
340 G4H2* proton_current_rmatrix = 0; << 355 G4AnaH2* proton_current_rmatrix =0;
341 356
342 if (pdg_nb == G4Electron::Electron()->GetP << 357 if (pdg_nb == G4Electron::Electron()->GetPDGEncoding()){ //e- matrices
343 edep_rmatrix = fEdep_rmatrix_vs_electron 358 edep_rmatrix = fEdep_rmatrix_vs_electron_prim_energy;
344 electron_current_rmatrix = fElectron_cur << 359 electron_current_rmatrix =
>> 360 fElectron_current_rmatrix_vs_electron_prim_energy;
345 gamma_current_rmatrix = fGamma_current_r 361 gamma_current_rmatrix = fGamma_current_rmatrix_vs_electron_prim_energy;
346 } 362 }
347 else if (pdg_nb == G4Gamma::Gamma()->GetPD << 363 else if (pdg_nb == G4Gamma::Gamma()->GetPDGEncoding()){
348 // gammma answer matrices << 364 //gammma answer matrices
349 edep_rmatrix = fEdep_rmatrix_vs_gamma_pr 365 edep_rmatrix = fEdep_rmatrix_vs_gamma_prim_energy;
350 electron_current_rmatrix = fElectron_cur 366 electron_current_rmatrix = fElectron_current_rmatrix_vs_gamma_prim_energy;
351 gamma_current_rmatrix = fGamma_current_r 367 gamma_current_rmatrix = fGamma_current_rmatrix_vs_gamma_prim_energy;
352 } 368 }
353 else if (pdg_nb == G4Proton::Proton()->Get << 369 else if (pdg_nb == G4Proton::Proton()->GetPDGEncoding()){
354 // proton answer matrices << 370 //proton answer matrices
355 edep_rmatrix = fEdep_rmatrix_vs_proton_p 371 edep_rmatrix = fEdep_rmatrix_vs_proton_prim_energy;
356 electron_current_rmatrix = fElectron_cur << 372 electron_current_rmatrix =
>> 373 fElectron_current_rmatrix_vs_proton_prim_energy;
357 gamma_current_rmatrix = fGamma_current_r 374 gamma_current_rmatrix = fGamma_current_rmatrix_vs_proton_prim_energy;
358 proton_current_rmatrix = fProton_current 375 proton_current_rmatrix = fProton_current_rmatrix_vs_proton_prim_energy;
359 } 376 }
360 // Register histo edep vs prim ekin << 377 //Register histo edep vs prim ekin
361 //---------------------------------- 378 //----------------------------------
362 if (normalised_weight > 0) fEdep_vs_prim_e << 379 if (normalised_weight>0) fEdep_vs_prim_ekin
>> 380 ->fill(prim_ekin,totEdep*normalised_weight);
363 // Registering answer matrix 381 // Registering answer matrix
364 //--------------------------- 382 //---------------------------
365 edep_rmatrix->fill(prim_ekin, totEdep * ad << 383 edep_rmatrix->fill(prim_ekin,totEdep*adj_weight/cm2);
366 << 384
367 // Registering of current of particles on << 385 //Registering of current of particles on the sensitive volume
368 //---------------------------------------- << 386 //------------------------------------------------------------
369 << 387
370 for (i = 0; i < electronCurrentCollection- << 388 for (i=0;i<electronCurrentCollection->entries();i++) {
371 G4double ekin = (*electronCurrentCollect << 389 G4double ekin =(*electronCurrentCollection)[i]->GetValue();
372 G4double weight = (*electronCurrentColle << 390 G4double weight=(*electronCurrentCollection)[i]->GetWeight();
373 fElectron_current->fill(ekin, weight * n << 391 fElectron_current->fill(ekin,weight*normalised_weight);
374 electron_current_rmatrix->fill(prim_ekin << 392 electron_current_rmatrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
375 } << 393 }
376 for (i = 0; i < protonCurrentCollection->e << 394 for (i=0;i<protonCurrentCollection->entries();i++) {
377 G4double ekin = (*protonCurrentCollectio << 395 G4double ekin =(*protonCurrentCollection)[i]->GetValue();
378 G4double weight = (*protonCurrentCollect << 396 G4double weight=(*protonCurrentCollection)[i]->GetWeight();
379 fProton_current->fill(ekin, weight * nor << 397 fProton_current->fill(ekin,weight*normalised_weight);
380 proton_current_rmatrix->fill(prim_ekin, << 398 proton_current_rmatrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
381 } << 399 }
382 for (i = 0; i < gammaCurrentCollection->en << 400 for (i=0;i<gammaCurrentCollection->entries();i++) {
383 G4double ekin = (*gammaCurrentCollection << 401 G4double ekin =(*gammaCurrentCollection)[i]->GetValue();
384 G4double weight = (*gammaCurrentCollecti << 402 G4double weight=(*gammaCurrentCollection)[i]->GetWeight();
385 fGamma_current->fill(ekin, weight * norm << 403 fGamma_current->fill(ekin,weight*normalised_weight);
386 gamma_current_rmatrix->fill(prim_ekin, e << 404 gamma_current_rmatrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
387 } << 405 }
388 } 406 }
389 407
390 // Registering of total energy deposited in << 408 //Registering of total energy deposited in Event
391 //------------------------------- 409 //-------------------------------
392 G4bool new_mean_computed = false; << 410 G4bool new_mean_computed=false;
393 if (totEdep > 0.) { << 411 if (totEdep>0.){
394 if (total_normalised_weight > 0.) { << 412 if (total_normalised_weight>0.){
395 G4double edep = totEdep * total_normalis << 413 G4double edep=totEdep* total_normalised_weight;
396 << 414
397 // Check if the edep is not wrongly too << 415 //Check if the edep is not wrongly too high
398 //-------------------------------------- << 416 //-----------------------------------------
399 G4double new_mean(0.0), new_error(0.0); << 417 G4double new_mean , new_error;
400 fAccumulated_edep += edep; << 418 fAccumulated_edep +=edep;
401 fAccumulated_edep2 += edep * edep; << 419 fAccumulated_edep2 +=edep*edep;
402 fNentry += 1.0; << 420 fNentry += 1.0;
403 ComputeMeanEdepAndError(anEvent, new_mea << 421 ComputeMeanEdepAndError(anEvent,new_mean,new_error);
404 G4double new_relative_error = 1.; << 422 G4double new_relative_error = 1.;
405 if (new_error > 0) new_relative_error = << 423 if ( new_error >0) new_relative_error = new_error/ new_mean;
406 if (fRelative_error < 0.10 && new_relati << 424 if (fRelative_error <0.10 && new_relative_error>1.5*fRelative_error) {
407 G4cout << "Potential wrong adjoint wei << 425 G4cout<<"Potential wrong adjoint weight!"<<std::endl;
408 G4cout << "The results of this event w << 426 G4cout<<"The results of this event will not be registered!"
409 G4cout << "previous mean edep [MeV] " << 427 <<std::endl;
410 G4cout << "previous relative error " < << 428 G4cout<<"previous mean edep [MeV] "<< fMean_edep<<std::endl;
411 G4cout << "new rejected mean edep [MeV << 429 G4cout<<"previous relative error "<< fRelative_error<<std::endl;
412 G4cout << "new rejected relative error << 430 G4cout<<"new rejected mean edep [MeV] "<< new_mean<<std::endl;
413 fAccumulated_edep -= edep; << 431 G4cout<<"new rejected relative error "<< new_relative_error
414 fAccumulated_edep2 -= edep * edep; << 432 <<std::endl;
415 fNentry -= 1.0; << 433 fAccumulated_edep -=edep;
416 return; << 434 fAccumulated_edep2 -=edep*edep;
417 } << 435 fNentry -= 1.0;
418 else { // accepted << 436 return;
419 fMean_edep = new_mean; << 437 }
420 fError_mean_edep = new_error; << 438 else { //accepted
421 fRelative_error = new_relative_error; << 439 fMean_edep = new_mean;
422 new_mean_computed = true; << 440 fError_mean_edep = new_error;
423 } << 441 fRelative_error =new_relative_error;
>> 442 new_mean_computed=true;
>> 443 }
>> 444
>> 445 }
>> 446 if (!new_mean_computed){
>> 447 ComputeMeanEdepAndError(anEvent,fMean_edep,fError_mean_edep);
>> 448 if (fError_mean_edep>0) fRelative_error= fError_mean_edep/fMean_edep;
424 } 449 }
425 << 450 }
426 if (!new_mean_computed) { <<
427 ComputeMeanEdepAndError(anEvent, fMean_e <<
428 fRelative_error = (fMean_edep > 0.0) ? f <<
429 } <<
430 } <<
431 } 451 }
432 452
433 //....oooOO0OOooo........oooOO0OOooo........oo 453 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
434 454
435 G4double RMC01AnalysisManager::PrimDiffAndDirF << 455 G4double RMC01AnalysisManager::PrimDiffAndDirFluxForAdjointSim(
436 { << 456 G4double prim_energy)
437 G4double flux = fAmplitude_prim_spectrum; << 457 {
438 if (fPrimSpectrumType == EXPO) << 458 G4double flux=fAmplitude_prim_spectrum;
439 flux *= std::exp(-prim_energy / fAlpha_or_ << 459 if ( fPrimSpectrumType ==EXPO) flux*=std::exp(-prim_energy/fAlpha_or_E0);
440 else << 460 else flux*=std::pow(prim_energy, -fAlpha_or_E0);
441 flux *= std::pow(prim_energy, -fAlpha_or_E <<
442 return flux; 461 return flux;
443 } 462 }
444 463
445 //....oooOO0OOooo........oooOO0OOooo........oo 464 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
446 /* 465 /*
447 void RMC01AnalysisManager::WriteHisto(G4H1* a << 466 void RMC01AnalysisManager::WriteHisto(G4AnaH1* anHisto,
448 G4double scaling_factor, G4String 467 G4double scaling_factor, G4String fileName, G4String header_lines)
449 { std::fstream FileOutput(fileName, std::ios:: 468 { std::fstream FileOutput(fileName, std::ios::out);
450 FileOutput<<header_lines; 469 FileOutput<<header_lines;
451 FileOutput.setf(std::ios::scientific); 470 FileOutput.setf(std::ios::scientific);
452 FileOutput.precision(6); 471 FileOutput.precision(6);
453 472
454 for (G4int i =0;i<G4int(anHisto->axis().bins << 473 for (G4int i =0;i<G4int(anHisto->axis().bins());i++) {
455 FileOutput<<anHisto->axis().bin_lower_ 474 FileOutput<<anHisto->axis().bin_lower_edge(i)
456 <<'\t'<<anHisto->axis().bin_uppe 475 <<'\t'<<anHisto->axis().bin_upper_edge(i)
457 <<'\t'<<anHisto->bin_height(i)*s 476 <<'\t'<<anHisto->bin_height(i)*scaling_factor
458 <<'\t'<<anHisto->bin_error(i)*sc 477 <<'\t'<<anHisto->bin_error(i)*scaling_factor<<std::endl;
459 } 478 }
460 } 479 }
461 480
462 //....oooOO0OOooo........oooOO0OOooo........oo 481 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
463 482
464 void RMC01AnalysisManager::WriteHisto(G4H2* a << 483 void RMC01AnalysisManager::WriteHisto(G4AnaH2* anHisto,
465 G4double scaling_factor, G4String 484 G4double scaling_factor, G4String fileName, G4String header_lines)
466 { std::fstream FileOutput(fileName, std::ios:: 485 { std::fstream FileOutput(fileName, std::ios::out);
467 FileOutput<<header_lines; 486 FileOutput<<header_lines;
468 << 487
469 FileOutput.setf(std::ios::scientific); 488 FileOutput.setf(std::ios::scientific);
470 FileOutput.precision(6); 489 FileOutput.precision(6);
471 490
472 for (G4int i =0;i<G4int(anHisto->axis_x().bi << 491 for (G4int i =0;i<G4int(anHisto->axis_x().bins());i++) {
473 for (G4int j =0;j<G4int(anHisto->axis_y(). << 492 for (G4int j =0;j<G4int(anHisto->axis_y().bins());j++) {
474 FileOutput<<anHisto->axis_x().bin_lower 493 FileOutput<<anHisto->axis_x().bin_lower_edge(i)
475 <<'\t'<<anHisto->axis_x() 494 <<'\t'<<anHisto->axis_x().bin_upper_edge(i)
476 <<'\t'<<anHisto->axis_y().b 495 <<'\t'<<anHisto->axis_y().bin_lower_edge(i)
477 <<'\t'<<anHisto->axis_y().b 496 <<'\t'<<anHisto->axis_y().bin_upper_edge(i)
478 <<'\t'<<anHisto->bin_height 497 <<'\t'<<anHisto->bin_height(i,j)*scaling_factor
479 <<'\t'<<anHisto->bin_error(i,j 498 <<'\t'<<anHisto->bin_error(i,j)*scaling_factor
480 499 <<std::endl;
481 } 500 }
482 } 501 }
483 } 502 }
484 */ 503 */
485 //....oooOO0OOooo........oooOO0OOooo........oo 504 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
486 505
487 void RMC01AnalysisManager::ComputeMeanEdepAndE << 506 void RMC01AnalysisManager::ComputeMeanEdepAndError(
488 << 507 const G4Event* anEvent,G4double& mean,G4double& error)
489 { << 508 {
490 G4int nb_event = anEvent->GetEventID() + 1; << 509 G4int nb_event=anEvent->GetEventID()+1;
491 G4double factor = 1.; << 510 G4double factor=1.;
492 if (fAdjoint_sim_mode) { << 511 if (fAdjoint_sim_mode) {
493 nb_event /= fNb_evt_per_adj_evt; << 512 nb_event /=fNb_evt_per_adj_evt;
494 factor = 1. * G4AdjointSimManager::GetInst << 513 factor=1.*G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
495 } << 514 }
496 << 515
497 // VI: error computation now is based on num << 516 // VI: error computation now is based on number of entries and not
498 // number of events << 517 // number of events
499 // LD: This is wrong! With the use of fNentr << 518 // LD: This is wrong! With the use of fNentry the results were no longer
500 // correctly normalised. The mean and th << 519 // correctly normalised. The mean and the error should be computed
501 // with nb_event. The old computation ha << 520 // with nb_event. The old computation has been reset.
502 // VI: OK, but let computations be double << 521 G4float nb_event_float = G4float(nb_event);
503 if (nb_event > 0) { << 522 if (nb_event_float >1.) {
504 G4double norm = 1.0 / (G4double)nb_event; << 523 mean = fAccumulated_edep/nb_event_float;
505 mean = fAccumulated_edep * norm; << 524 G4double mean_x2 = fAccumulated_edep2/nb_event_float;
506 G4double mean_x2 = fAccumulated_edep2 * no << 525 /*
507 G4double zz = mean_x2 - mean * mean; << 526 G4cout << "Nevt= " << nb_event << " mean= " << mean
508 /* << 527 << " mean_x2= " << mean_x2 << " x2 - x*x= "
509 G4cout << "Nevt= " << nb_event << " mean << 528 << mean_x2-mean*mean << G4endl;
510 << " mean_x2= " << mean_x2 << " << 529 */
511 << zz << G4endl; << 530 error = factor*std::sqrt(mean_x2-mean*mean)/std::sqrt(nb_event_float);
512 */ << 531 mean *=factor;
513 error = factor * std::sqrt(std::max(zz, 0. << 532 }
514 mean *= factor; << 533 else {
515 } << 534 mean=0;
516 else { << 535 error=0;
517 mean = 0; <<
518 error = 0; <<
519 } 536 }
520 // G4cout << "Aend: " << mean << " " << erro <<
521 } 537 }
522 538
523 //....oooOO0OOooo........oooOO0OOooo........oo 539 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
524 540
525 void RMC01AnalysisManager::SetPrimaryExpSpectr << 541 void RMC01AnalysisManager::SetPrimaryExpSpectrumForAdjointSim(
526 << 542 const G4String& particle_name, G4double omni_fluence,
527 << 543 G4double E0, G4double Emin,
528 { << 544 G4double Emax)
529 fPrimSpectrumType = EXPO; << 545 { fPrimSpectrumType = EXPO;
530 if (particle_name == "e-") << 546 if (particle_name == "e-" ) fPrimPDG_ID =
531 fPrimPDG_ID = G4Electron::Electron()->GetP << 547 G4Electron::Electron()->GetPDGEncoding();
532 else if (particle_name == "gamma") << 548 else if (particle_name == "gamma") fPrimPDG_ID =
533 fPrimPDG_ID = G4Gamma::Gamma()->GetPDGEnco << 549 G4Gamma::Gamma()->GetPDGEncoding();
534 else if (particle_name == "proton") << 550 else if (particle_name == "proton") fPrimPDG_ID =
535 fPrimPDG_ID = G4Proton::Proton()->GetPDGEn << 551 G4Proton::Proton()->GetPDGEncoding();
536 else { 552 else {
537 G4cout << "The particle that you did selec << 553 G4cout<<"The particle that you did select is not in the candidate "<<
538 << "list for primary [e-, gamma, pr << 554 "list for primary [e-, gamma, proton]!"<<G4endl;
539 return; << 555 return;
540 } << 556 }
541 fAlpha_or_E0 = E0; << 557 fAlpha_or_E0 = E0 ;
542 fAmplitude_prim_spectrum = << 558 fAmplitude_prim_spectrum = omni_fluence/E0/
543 omni_fluence / E0 / (std::exp(-Emin / E0) << 559 (std::exp(-Emin/E0)-std::exp(-Emax/E0))/4./pi;
544 fEmin_prim_spectrum = Emin; << 560 fEmin_prim_spectrum = Emin ;
545 fEmax_prim_spectrum = Emax; 561 fEmax_prim_spectrum = Emax;
546 } 562 }
547 563
548 //....oooOO0OOooo........oooOO0OOooo........oo 564 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
549 565
550 void RMC01AnalysisManager::SetPrimaryPowerLawS << 566 void RMC01AnalysisManager::SetPrimaryPowerLawSpectrumForAdjointSim(
551 << 567 const G4String& particle_name, G4double omni_fluence,
552 << 568 G4double alpha, G4double Emin,G4double Emax)
553 << 569 { fPrimSpectrumType =POWER;
554 { << 570 if (particle_name == "e-" ) fPrimPDG_ID =
555 fPrimSpectrumType = POWER; << 571 G4Electron::Electron()->GetPDGEncoding();
556 if (particle_name == "e-") << 572 else if (particle_name == "gamma") fPrimPDG_ID =
557 fPrimPDG_ID = G4Electron::Electron()->GetP << 573 G4Gamma::Gamma()->GetPDGEncoding();
558 else if (particle_name == "gamma") << 574 else if (particle_name == "proton") fPrimPDG_ID =
559 fPrimPDG_ID = G4Gamma::Gamma()->GetPDGEnco << 575 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 { 576 else {
572 G4double p = 1. - alpha; << 577 G4cout<<"The particle that you did select is not in the candidate"<<
573 fAmplitude_prim_spectrum = omni_fluence / << 578 " list for primary [e-, gamma, proton]!"<<G4endl;
574 } << 579 return;
>> 580 }
>> 581
>> 582
>> 583 if (alpha ==1.) {
>> 584 fAmplitude_prim_spectrum = omni_fluence/std::log(Emax/Emin)/4./pi;
>> 585 }
>> 586 else {
>> 587 G4double p=1.-alpha;
>> 588 fAmplitude_prim_spectrum = omni_fluence/p/(std::pow(Emax,p)
>> 589 -std::pow(Emin,p))/4./pi;
>> 590 }
575 591
576 fAlpha_or_E0 = alpha; << 592 fAlpha_or_E0 = alpha ;
577 fEmin_prim_spectrum = Emin; << 593 fEmin_prim_spectrum = Emin ;
578 fEmax_prim_spectrum = Emax; 594 fEmax_prim_spectrum = Emax;
>> 595
579 } 596 }
580 597
581 //....oooOO0OOooo........oooOO0OOooo........oo 598 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
582 599
583 void RMC01AnalysisManager::Book() 600 void RMC01AnalysisManager::Book()
584 { 601 {
585 //---------------------- 602 //----------------------
586 // Creation of histograms << 603 //Creation of histograms
587 //---------------------- 604 //----------------------
588 605
589 // Energy binning of the histograms : 60 log << 606 //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 607
611 // Create directories << 608 G4double emin=1.*keV;
612 theHistoManager->SetHistoDirectoryName("hist << 609 G4double emax=1.*GeV;
613 610
614 // Histograms for : << 611 //file_name
615 // 1)the forward simulation results << 612 fFileName[0]="forward_sim";
616 // 2)the Reverse MC simulation resu << 613 if (fAdjoint_sim_mode) fFileName[0]="adjoint_sim";
>> 614
>> 615 //Histo manager
>> 616 G4AnalysisManager* theHistoManager = G4AnalysisManager::Instance();
>> 617 G4String extension = theHistoManager->GetFileType();
>> 618 fFileName[1] = fFileName[0] + "." + extension;
>> 619 theHistoManager->SetFirstHistoId(1);
>> 620
>> 621 G4bool fileOpen = theHistoManager->OpenFile(fFileName[0]);
>> 622 if (!fileOpen) {
>> 623 G4cout << "\n---> RMC01AnalysisManager::Book(): cannot open "
>> 624 << fFileName[1]
>> 625 << G4endl;
>> 626 return;
>> 627 }
>> 628
>> 629 // Create directories
>> 630 theHistoManager->SetHistoDirectoryName("histo");
>> 631
>> 632 //Histograms for :
>> 633 // 1)the forward simulation results
>> 634 // 2)the Reverse MC simulation results normalised to a user spectrum
617 //------------------------------------------ 635 //------------------------------------------------------------------------
618 636
619 G4int idHisto = << 637 G4int idHisto =
620 theHistoManager->CreateH1(G4String("Edep_v << 638 theHistoManager->CreateH1(G4String("Edep_vs_prim_ekin"),
621 60, emin, emax, << 639 G4String("edep vs e- primary energy"),60,emin,emax,
>> 640 "none","none",G4String("log"));
622 fEdep_vs_prim_ekin = theHistoManager->GetH1( 641 fEdep_vs_prim_ekin = theHistoManager->GetH1(idHisto);
623 642
624 idHisto = theHistoManager->CreateH1(G4String << 643 idHisto = theHistoManager->CreateH1(G4String("elecron_current"),
625 emax, "n << 644 G4String("electron"),60,emin,emax,
>> 645 "none","none",G4String("log"));
>> 646
>> 647 fElectron_current = theHistoManager->GetH1(idHisto);
>> 648
>> 649 idHisto= theHistoManager->CreateH1(G4String("proton_current"),
>> 650 G4String("proton"),60,emin,emax,
>> 651 "none","none",G4String("log"));
>> 652 fProton_current=theHistoManager->GetH1(idHisto);
>> 653
>> 654 idHisto= theHistoManager->CreateH1(G4String("gamma_current"),
>> 655 G4String("gamma"),60,emin,emax,
>> 656 "none","none",G4String("log"));
>> 657 fGamma_current=theHistoManager->GetH1(idHisto);
626 658
627 fElectron_current = theHistoManager->GetH1(i << 659 //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 //------------------------------------------ 660 //-----------------------------------------------
639 if (fAdjoint_sim_mode) { << 661 if (fAdjoint_sim_mode){
640 // Response matrices for external isotropi << 662 //Response matrices for external isotropic e- source
641 //---------------------------------------- << 663 //--------------------------------------------------
642 << 664
643 idHisto = theHistoManager->CreateH1(G4Stri << 665 idHisto =
644 G4Stri << 666 theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_electron_prim_energy"),
645 emax, << 667 G4String("electron RM vs e- primary energy"),60,emin,emax,
646 fEdep_rmatrix_vs_electron_prim_energy = th << 668 "none","none",G4String("log"));
647 << 669 fEdep_rmatrix_vs_electron_prim_energy = theHistoManager->GetH1(idHisto);
648 idHisto = theHistoManager->CreateH2( << 670
649 G4String("Electron_current_rmatrix_vs_el << 671 idHisto =
650 G4String("electron current RM vs e- pri << 672 theHistoManager->
651 "none", "none", "none", G4String("log"), << 673 CreateH2(G4String("Electron_current_rmatrix_vs_electron_prim_energy"),
652 << 674 G4String("electron current RM vs e- primary energy"),
653 fElectron_current_rmatrix_vs_electron_prim << 675 60,emin,emax,60,emin,emax,
654 << 676 "none","none","none","none",G4String("log"),G4String("log"));
655 idHisto = theHistoManager->CreateH2(G4Stri << 677
656 G4Stri << 678 fElectron_current_rmatrix_vs_electron_prim_energy =
657 emin, << 679 theHistoManager->GetH2(idHisto);
658 G4Stri << 680
>> 681 idHisto =
>> 682 theHistoManager->
>> 683 CreateH2(G4String("Gamma_current_rmatrix_vs_electron_prim_energy"),
>> 684 G4String("gamma current RM vs e- primary energy"),
>> 685 60,emin,emax,60,emin,emax,
>> 686 "none","none","none","none",G4String("log"),G4String("log"));
>> 687
>> 688 fGamma_current_rmatrix_vs_electron_prim_energy =
>> 689 theHistoManager->GetH2(idHisto);
>> 690
>> 691 //Response matrices for external isotropic gamma source
>> 692
>> 693 idHisto =
>> 694 theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_gamma_prim_energy"),
>> 695 G4String("electron RM vs gamma primary energy"),60,emin,emax,
>> 696 "none","none",G4String("log"));
>> 697 fEdep_rmatrix_vs_gamma_prim_energy = theHistoManager->GetH1(idHisto);
>> 698
>> 699 idHisto =
>> 700 theHistoManager->
>> 701 CreateH2(G4String("Electron_current_rmatrix_vs_gamma_prim_energy"),
>> 702 G4String("electron current RM vs gamma primary energy"),
>> 703 60,emin,emax,60,emin,emax,
>> 704 "none","none","none","none",G4String("log"),G4String("log"));
>> 705
>> 706 fElectron_current_rmatrix_vs_gamma_prim_energy =
>> 707 theHistoManager->GetH2(idHisto);
>> 708
>> 709 idHisto =
>> 710 theHistoManager->
>> 711 CreateH2(G4String("Gamma_current_rmatrix_vs_gamma_prim_energy"),
>> 712 G4String("gamma current RM vs gamma primary energy"),
>> 713 60,emin,emax,60,emin,emax,
>> 714 "none","none","none","none",G4String("log"),G4String("log"));
>> 715
>> 716 fGamma_current_rmatrix_vs_gamma_prim_energy =
>> 717 theHistoManager->GetH2(idHisto);
>> 718
>> 719 //Response matrices for external isotropic proton source
>> 720 idHisto =
>> 721 theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_proton_prim_energy"),
>> 722 G4String("electron RM vs proton primary energy"),60,emin,emax,
>> 723 "none","none",G4String("log"));
>> 724 fEdep_rmatrix_vs_proton_prim_energy = theHistoManager->GetH1(idHisto);
>> 725
>> 726 idHisto =
>> 727 theHistoManager->
>> 728 CreateH2(G4String("Electron_current_rmatrix_vs_proton_prim_energy"),
>> 729 G4String("electron current RM vs proton primary energy"),
>> 730 60,emin,emax,60,emin,emax,
>> 731 "none","none","none","none",G4String("log"),G4String("log"));
>> 732
>> 733 fElectron_current_rmatrix_vs_proton_prim_energy =
>> 734 theHistoManager->GetH2(idHisto);
>> 735
>> 736 idHisto =
>> 737 theHistoManager->
>> 738 CreateH2(G4String("Gamma_current_rmatrix_vs_proton_prim_energy"),
>> 739 G4String("gamma current RM vs proton primary energy"),
>> 740 60,emin,emax,60,emin,emax,
>> 741 "none","none","none","none",G4String("log"),G4String("log"));
>> 742
>> 743 fGamma_current_rmatrix_vs_proton_prim_energy =
>> 744 theHistoManager->GetH2(idHisto);
>> 745
>> 746 idHisto =
>> 747 theHistoManager->
>> 748 CreateH2(G4String("Proton_current_rmatrix_vs_proton_prim_energy"),
>> 749 G4String("proton current RM vs proton primary energy"),
>> 750 60,emin,emax,60,emin,emax,
>> 751 "none","none","none","none",G4String("log"),G4String("log"));
659 752
660 fGamma_current_rmatrix_vs_electron_prim_en << 753 fProton_current_rmatrix_vs_proton_prim_energy =
661 << 754 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 } 755 }
710 fFactoryOn = true; 756 fFactoryOn = true;
711 G4cout << "\n----> Histogram Tree is opened 757 G4cout << "\n----> Histogram Tree is opened in " << fFileName[1] << G4endl;
712 } 758 }
713 759
714 //....oooOO0OOooo........oooOO0OOooo........oo 760 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
715 761
716 void RMC01AnalysisManager::Save(G4double scali 762 void RMC01AnalysisManager::Save(G4double scaling_factor)
717 { << 763 { if (fFactoryOn) {
718 if (fFactoryOn) { <<
719 G4AnalysisManager* theHistoManager = G4Ana 764 G4AnalysisManager* theHistoManager = G4AnalysisManager::Instance();
720 // scaling of results << 765 //scaling of results
721 //----------------- 766 //-----------------
722 767
723 for (G4int ind = 1; ind <= theHistoManager << 768 for (int ind=1; ind<=theHistoManager->GetNofH1s();ind++){
724 theHistoManager->SetH1Ascii(ind, true); << 769 theHistoManager->SetH1Ascii(ind,true);
725 theHistoManager->ScaleH1(ind, scaling_fa << 770 theHistoManager->ScaleH1(ind,scaling_factor);
726 } 771 }
727 for (G4int ind = 1; ind <= theHistoManager << 772 for (int ind=1; ind<=theHistoManager->GetNofH2s();ind++)
728 theHistoManager->ScaleH2(ind, scaling_fa << 773 theHistoManager->ScaleH2(ind,scaling_factor);
729 774
730 theHistoManager->Write(); 775 theHistoManager->Write();
731 theHistoManager->CloseFile(); 776 theHistoManager->CloseFile();
732 G4cout << "\n----> Histogram Tree is saved 777 G4cout << "\n----> Histogram Tree is saved in " << fFileName[1] << G4endl;
733 778
734 theHistoManager->Clear(); << 779 delete G4AnalysisManager::Instance();
735 fFactoryOn = false; 780 fFactoryOn = false;
736 } 781 }
737 } 782 }
738 783