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