| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 /// \file biasing/ReverseMC01/src/RMC01Analysi
27 /// \brief Implementation of the RMC01Analysis
28 //
29 //
30 //////////////////////////////////////////////
31 // Class Name: RMC01AnalysisManage
32 // Author: L. Desorgher
33 // Organisation: SpaceIT GmbH
34 // Contract: ESA contract 21435/
35 // Customer: ESA/ESTEC
36 //////////////////////////////////////////////
37
38 //....oooOO0OOooo........oooOO0OOooo........oo
39 //....oooOO0OOooo........oooOO0OOooo........oo
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::fI
58
59 //....oooOO0OOooo........oooOO0OOooo........oo
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
73 fGamma_current_rmatrix_vs_electron_prim_en
74 fEdep_rmatrix_vs_gamma_prim_energy(0),
75 fElectron_current_rmatrix_vs_gamma_prim_en
76 fGamma_current_rmatrix_vs_gamma_prim_energ
77 fEdep_rmatrix_vs_proton_prim_energy(0),
78 fElectron_current_rmatrix_vs_proton_prim_e
79 fProton_current_rmatrix_vs_proton_prim_ene
80 fGamma_current_rmatrix_vs_proton_prim_ener
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(thi
91
92 //-------------
93 // Timer for convergence vector
94 //-------------
95
96 fTimer = new G4Timer();
97
98 //---------------------------------
99 // Primary particle ID for normalisation of
100 //---------------------------------
101
102 fPrimPDG_ID = G4Electron::Electron()->GetPDG
103
104 fFileName[0] = "sim";
105 }
106
107 //....oooOO0OOooo........oooOO0OOooo........oo
108
109 RMC01AnalysisManager::~RMC01AnalysisManager()
110 {
111 ;
112 }
113
114 //....oooOO0OOooo........oooOO0OOooo........oo
115
116 RMC01AnalysisManager* RMC01AnalysisManager::Ge
117 {
118 if (fInstance == 0) fInstance = new RMC01Ana
119 return fInstance;
120 }
121
122 //....oooOO0OOooo........oooOO0OOooo........oo
123
124 void RMC01AnalysisManager::BeginOfRun(const G4
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::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
147 fTimer->Start();
148 fElapsed_time = 0.;
149
150 Book();
151 }
152
153 //....oooOO0OOooo........oooOO0OOooo........oo
154
155 void RMC01AnalysisManager::EndOfRun(const G4Ru
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!"
162 G4cout << "edep per event [MeV] = " << fMe
163 G4cout << "error[MeV] = " << fError_mean_e
164 }
165
166 else {
167 G4cout << "Results of reverse/adjoint simu
168 G4cout << "normalised edep [MeV] = " << fM
169 G4cout << "error[MeV] = " << fError_mean_e
170 factor = 1. * G4AdjointSimManager::GetInst
171 / aRun->GetNumberOfEvent();
172 }
173 Save(factor);
174 fConvergenceFileOutput.close();
175 }
176
177 //....oooOO0OOooo........oooOO0OOooo........oo
178
179 void RMC01AnalysisManager::BeginOfEvent(const
180 {
181 ;
182 }
183
184 //....oooOO0OOooo........oooOO0OOooo........oo
185
186 void RMC01AnalysisManager::EndOfEvent(const G4
187 {
188 if (fAdjoint_sim_mode)
189 EndOfEventForAdjointSimulation(anEvent);
190 else
191 EndOfEventForForwardSimulation(anEvent);
192
193 // Test convergence. The error is already co
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_pe
199 // nb_event/fNb_evt_per_adj_evt;
200 if (n_adj_evt * fNb_evt_per_adj_evt == nb_
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
208 G4cout << fPrecision_to_reach * 100. << "%
209 fTimer->Stop();
210 fElapsed_time += fTimer->GetRealElapsed();
211 fConvergenceFileOutput << fMean_edep << '\
212 << std::endl;
213 G4RunManager::GetRunManager()->AbortRun(tr
214 }
215
216 if (nb_event > 0 && nb_event % fNb_evt_modul
217 fTimer->Stop();
218 fElapsed_time += fTimer->GetRealElapsed();
219 fTimer->Start();
220 fConvergenceFileOutput << fMean_edep << '\
221 << std::endl;
222 }
223 }
224
225 //....oooOO0OOooo........oooOO0OOooo........oo
226
227 void RMC01AnalysisManager::EndOfEventForForwar
228 {
229 G4SDManager* SDman = G4SDManager::GetSDMpoin
230 G4HCofThisEvent* HCE = anEvent->GetHCofThisE
231 RMC01DoubleWithWeightHitsCollection* edepCol
232 (RMC01DoubleWithWeightHitsCollection*)(HCE
233
234 RMC01DoubleWithWeightHitsCollection* electro
235 (RMC01DoubleWithWeightHitsCollection*)(HCE
236
237 RMC01DoubleWithWeightHitsCollection* protonC
238 (RMC01DoubleWithWeightHitsCollection*)(HCE
239
240 RMC01DoubleWithWeightHitsCollection* gammaCu
241 (RMC01DoubleWithWeightHitsCollection*)(HCE
242
243 // Total energy deposited in Event
244 //-------------------------------
245 G4double totEdep = 0;
246 std::size_t i;
247 for (i = 0; i < edepCollection->entries(); +
248 totEdep += (*edepCollection)[i]->GetValue(
249
250 if (totEdep > 0.) {
251 fAccumulated_edep += totEdep;
252 fAccumulated_edep2 += totEdep * totEdep;
253 fNentry += 1.0;
254 G4PrimaryParticle* thePrimary = anEvent->G
255 G4double E0 = thePrimary->GetG4code()->Get
256 G4double P = thePrimary->GetMomentum().mag
257 G4double prim_ekin = std::sqrt(E0 * E0 + P
258 fEdep_vs_prim_ekin->fill(prim_ekin, totEde
259 }
260 ComputeMeanEdepAndError(anEvent, fMean_edep,
261 if (fError_mean_edep > 0) fRelative_error =
262
263 // Particle current on sensitive cylinder
264 //-------------------------------------
265
266 for (i = 0; i < electronCurrentCollection->e
267 G4double ekin = (*electronCurrentCollectio
268 G4double weight = (*electronCurrentCollect
269 fElectron_current->fill(ekin, weight);
270 }
271
272 for (i = 0; i < protonCurrentCollection->ent
273 G4double ekin = (*protonCurrentCollection)
274 G4double weight = (*protonCurrentCollectio
275 fProton_current->fill(ekin, weight);
276 }
277
278 for (i = 0; i < gammaCurrentCollection->entr
279 G4double ekin = (*gammaCurrentCollection)[
280 G4double weight = (*gammaCurrentCollection
281 fGamma_current->fill(ekin, weight);
282 }
283 }
284
285 //....oooOO0OOooo........oooOO0OOooo........oo
286
287 void RMC01AnalysisManager::EndOfEventForAdjoin
288 {
289 // Output from Sensitive volume computed dur
290 //------------------------------------------
291 G4SDManager* SDman = G4SDManager::GetSDMpoin
292 G4HCofThisEvent* HCE = anEvent->GetHCofThisE
293 RMC01DoubleWithWeightHitsCollection* edepCol
294 (RMC01DoubleWithWeightHitsCollection*)(HCE
295
296 RMC01DoubleWithWeightHitsCollection* electro
297 (RMC01DoubleWithWeightHitsCollection*)(HCE
298
299 RMC01DoubleWithWeightHitsCollection* protonC
300 (RMC01DoubleWithWeightHitsCollection*)(HCE
301
302 RMC01DoubleWithWeightHitsCollection* gammaCu
303 (RMC01DoubleWithWeightHitsCollection*)(HCE
304
305 // Computation of total energy deposited in
306 //-------------------------------
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
315 G4AdjointSimManager* theAdjointSimManager =
316
317 size_t nb_adj_track = theAdjointSimManager->
318 G4double total_normalised_weight = 0.;
319
320 // We need to loop over the adjoint tracks t
321 // surface.
322 for (std::size_t j = 0; j < nb_adj_track; ++
323 G4int pdg_nb = theAdjointSimManager->GetFw
324 G4double prim_ekin = theAdjointSimManager-
325 G4double adj_weight = theAdjointSimManager
326
327 // Factor of normalisation to user defined
328 //----------------------------------------
329 G4double normalised_weight = 0.;
330 if (pdg_nb == fPrimPDG_ID && prim_ekin >=
331 && prim_ekin <= fEmax_prim_spectrum)
332 normalised_weight = adj_weight * PrimDif
333 total_normalised_weight += normalised_weig
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()->GetP
343 edep_rmatrix = fEdep_rmatrix_vs_electron
344 electron_current_rmatrix = fElectron_cur
345 gamma_current_rmatrix = fGamma_current_r
346 }
347 else if (pdg_nb == G4Gamma::Gamma()->GetPD
348 // gammma answer matrices
349 edep_rmatrix = fEdep_rmatrix_vs_gamma_pr
350 electron_current_rmatrix = fElectron_cur
351 gamma_current_rmatrix = fGamma_current_r
352 }
353 else if (pdg_nb == G4Proton::Proton()->Get
354 // proton answer matrices
355 edep_rmatrix = fEdep_rmatrix_vs_proton_p
356 electron_current_rmatrix = fElectron_cur
357 gamma_current_rmatrix = fGamma_current_r
358 proton_current_rmatrix = fProton_current
359 }
360 // Register histo edep vs prim ekin
361 //----------------------------------
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 //-------------------------------
392 G4bool new_mean_computed = false;
393 if (totEdep > 0.) {
394 if (total_normalised_weight > 0.) {
395 G4double edep = totEdep * total_normalis
396
397 // Check if the edep is not wrongly too
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_mea
404 G4double new_relative_error = 1.;
405 if (new_error > 0) new_relative_error =
406 if (fRelative_error < 0.10 && new_relati
407 G4cout << "Potential wrong adjoint wei
408 G4cout << "The results of this event w
409 G4cout << "previous mean edep [MeV] "
410 G4cout << "previous relative error " <
411 G4cout << "new rejected mean edep [MeV
412 G4cout << "new rejected relative error
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_e
428 fRelative_error = (fMean_edep > 0.0) ? f
429 }
430 }
431 }
432
433 //....oooOO0OOooo........oooOO0OOooo........oo
434
435 G4double RMC01AnalysisManager::PrimDiffAndDirF
436 {
437 G4double flux = fAmplitude_prim_spectrum;
438 if (fPrimSpectrumType == EXPO)
439 flux *= std::exp(-prim_energy / fAlpha_or_
440 else
441 flux *= std::pow(prim_energy, -fAlpha_or_E
442 return flux;
443 }
444
445 //....oooOO0OOooo........oooOO0OOooo........oo
446 /*
447 void RMC01AnalysisManager::WriteHisto(G4H1* a
448 G4double scaling_factor, G4String
449 { std::fstream FileOutput(fileName, std::ios::
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
455 FileOutput<<anHisto->axis().bin_lower_
456 <<'\t'<<anHisto->axis().bin_uppe
457 <<'\t'<<anHisto->bin_height(i)*s
458 <<'\t'<<anHisto->bin_error(i)*sc
459 }
460 }
461
462 //....oooOO0OOooo........oooOO0OOooo........oo
463
464 void RMC01AnalysisManager::WriteHisto(G4H2* a
465 G4double scaling_factor, G4String
466 { std::fstream FileOutput(fileName, std::ios::
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().bi
473 for (G4int j =0;j<G4int(anHisto->axis_y().
474 FileOutput<<anHisto->axis_x().bin_lower
475 <<'\t'<<anHisto->axis_x()
476 <<'\t'<<anHisto->axis_y().b
477 <<'\t'<<anHisto->axis_y().b
478 <<'\t'<<anHisto->bin_height
479 <<'\t'<<anHisto->bin_error(i,j
480
481 }
482 }
483 }
484 */
485 //....oooOO0OOooo........oooOO0OOooo........oo
486
487 void RMC01AnalysisManager::ComputeMeanEdepAndE
488
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::GetInst
495 }
496
497 // VI: error computation now is based on num
498 // number of events
499 // LD: This is wrong! With the use of fNentr
500 // correctly normalised. The mean and th
501 // with nb_event. The old computation ha
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 * no
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 }
520 // G4cout << "Aend: " << mean << " " << erro
521 }
522
523 //....oooOO0OOooo........oooOO0OOooo........oo
524
525 void RMC01AnalysisManager::SetPrimaryExpSpectr
526
527
528 {
529 fPrimSpectrumType = EXPO;
530 if (particle_name == "e-")
531 fPrimPDG_ID = G4Electron::Electron()->GetP
532 else if (particle_name == "gamma")
533 fPrimPDG_ID = G4Gamma::Gamma()->GetPDGEnco
534 else if (particle_name == "proton")
535 fPrimPDG_ID = G4Proton::Proton()->GetPDGEn
536 else {
537 G4cout << "The particle that you did selec
538 << "list for primary [e-, gamma, pr
539 return;
540 }
541 fAlpha_or_E0 = E0;
542 fAmplitude_prim_spectrum =
543 omni_fluence / E0 / (std::exp(-Emin / E0)
544 fEmin_prim_spectrum = Emin;
545 fEmax_prim_spectrum = Emax;
546 }
547
548 //....oooOO0OOooo........oooOO0OOooo........oo
549
550 void RMC01AnalysisManager::SetPrimaryPowerLawS
551
552
553
554 {
555 fPrimSpectrumType = POWER;
556 if (particle_name == "e-")
557 fPrimPDG_ID = G4Electron::Electron()->GetP
558 else if (particle_name == "gamma")
559 fPrimPDG_ID = G4Gamma::Gamma()->GetPDGEnco
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 {
572 G4double p = 1. - alpha;
573 fAmplitude_prim_spectrum = omni_fluence /
574 }
575
576 fAlpha_or_E0 = alpha;
577 fEmin_prim_spectrum = Emin;
578 fEmax_prim_spectrum = Emax;
579 }
580
581 //....oooOO0OOooo........oooOO0OOooo........oo
582
583 void RMC01AnalysisManager::Book()
584 {
585 //----------------------
586 // Creation of histograms
587 //----------------------
588
589 // Energy binning of the histograms : 60 log
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
611 // Create directories
612 theHistoManager->SetHistoDirectoryName("hist
613
614 // Histograms for :
615 // 1)the forward simulation results
616 // 2)the Reverse MC simulation resu
617 //------------------------------------------
618
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
633 idHisto = theHistoManager->CreateH1(G4String
634 "none",
635 fGamma_current = theHistoManager->GetH1(idHi
636
637 // Response matrices for the adjoint simulat
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
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 }
710 fFactoryOn = true;
711 G4cout << "\n----> Histogram Tree is opened
712 }
713
714 //....oooOO0OOooo........oooOO0OOooo........oo
715
716 void RMC01AnalysisManager::Save(G4double scali
717 {
718 if (fFactoryOn) {
719 G4AnalysisManager* theHistoManager = G4Ana
720 // scaling of results
721 //-----------------
722
723 for (G4int ind = 1; ind <= theHistoManager
724 theHistoManager->SetH1Ascii(ind, true);
725 theHistoManager->ScaleH1(ind, scaling_fa
726 }
727 for (G4int ind = 1; ind <= theHistoManager
728 theHistoManager->ScaleH2(ind, scaling_fa
729
730 theHistoManager->Write();
731 theHistoManager->CloseFile();
732 G4cout << "\n----> Histogram Tree is saved
733
734 theHistoManager->Clear();
735 fFactoryOn = false;
736 }
737 }
738