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