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