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