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