Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Static Public Member Functions
RMC01AnalysisManager Class Reference

#include <RMC01AnalysisManager.hh>

Public Member Functions

 ~RMC01AnalysisManager ()
 
void BeginOfRun (const G4Run *)
 
void EndOfRun (const G4Run *)
 
void BeginOfEvent (const G4Event *)
 
void EndOfEvent (const G4Event *)
 
void SetPrimaryExpSpectrumForAdjointSim (const G4String &particle_name, G4double fluence, G4double E0, G4double Emin, G4double Emax)
 
void SetPrimaryPowerLawSpectrumForAdjointSim (const G4String &particle_name, G4double fluence, G4double alpha, G4double Emin, G4double Emax)
 
void SetPrecision (G4double precision)
 
void book ()
 
void save (G4double scaling_factor)
 

Static Public Member Functions

static RMC01AnalysisManager * GetInstance ()
 

Detailed Description

Definition at line 77 of file RMC01AnalysisManager.hh.

Constructor & Destructor Documentation

RMC01AnalysisManager::~RMC01AnalysisManager ( )

Definition at line 103 of file RMC01AnalysisManager.cc.

104 {;
105 }

Member Function Documentation

void RMC01AnalysisManager::BeginOfEvent ( const G4Event *  )

Definition at line 177 of file RMC01AnalysisManager.cc.

Referenced by RMC01EventAction::BeginOfEventAction().

178 { ;
179 }
void RMC01AnalysisManager::BeginOfRun ( const G4Run *  aRun)

Definition at line 117 of file RMC01AnalysisManager.cc.

References book(), G4AdjointSimManager::GetAdjointSimMode(), G4AdjointSimManager::GetInstance(), G4AdjointSimManager::GetNbEvtOfLastRun(), G4Run::GetNumberOfEventToBeProcessed(), and G4Timer::Start().

Referenced by RMC01RunAction::BeginOfRunAction().

118 {
119 
120 
121  fAccumulated_edep =0.;
122  fAccumulated_edep2 =0.;
123  fRelative_error=1.;
124  fMean_edep=0.;
125  fError_mean_edep=0.;
126  fAdjoint_sim_mode =G4AdjointSimManager::GetInstance()->GetAdjointSimMode();
127 
128  if (fAdjoint_sim_mode){
129  fNb_evt_per_adj_evt=aRun->GetNumberOfEventToBeProcessed()/
130  G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
131  fConvergenceFileOutput.open("ConvergenceOfAdjointSimulationResults.txt",
132  std::ios::out);
133  fConvergenceFileOutput<<
134  "Normalised Edep[MeV]\terror[MeV]\tcomputing_time[s]"<<std::endl;
135  }
136  else {
137  fConvergenceFileOutput.open("ConvergenceOfForwardSimulationResults.txt",
138  std::ios::out);
139  fConvergenceFileOutput<<
140  "Edep per event [MeV]\terror[MeV]\tcomputing_time[s]"
141  <<std::endl;
142  }
143  fConvergenceFileOutput.setf(std::ios::scientific);
144  fConvergenceFileOutput.precision(6);
145 
146  fTimer->Start();
147  fElapsed_time=0.;
148 
149  book();
150 }
G4AdjointSimManager::GetInstance
static G4AdjointSimManager * GetInstance()
Definition: G4AdjointSimManager.cc:129
G4Run::GetNumberOfEventToBeProcessed
G4int GetNumberOfEventToBeProcessed() const
Definition: G4Run.hh:83
G4Timer::Start
void Start()
void RMC01AnalysisManager::book ( )

Definition at line 580 of file RMC01AnalysisManager.cc.

References G4VAnalysisManager::CreateH1(), G4cout, G4endl, G4VAnalysisManager::GetFileType(), G4RootAnalysisManager::GetH1(), G4RootAnalysisManager::GetH2(), python.hepunit::GeV, python.hepunit::keV, G4VAnalysisManager::OpenFile(), G4VAnalysisManager::SetFirstHistoId(), and G4VAnalysisManager::SetHistoDirectoryName().

Referenced by BeginOfRun().

581 {
582  //----------------------
583  //Creation of histograms
584  //----------------------
585 
586  //Energy binning of the histograms : 60 log bins over [1keV-1GeV]
587 
588  G4double emin=1.*keV;
589  G4double emax=1.*GeV;
590 
591  //file_name
592  fFileName[0]="forward_sim";
593  if (fAdjoint_sim_mode) fFileName[0]="adjoint_sim";
594 
595  //Histo manager
596  G4AnalysisManager* theHistoManager = G4AnalysisManager::Instance();
597  G4String extension = theHistoManager->GetFileType();
598  fFileName[1] = fFileName[0] + "." + extension;
599  theHistoManager->SetFirstHistoId(1);
600 
601  G4bool fileOpen = theHistoManager->OpenFile(fFileName[0]);
602  if (!fileOpen) {
603  G4cout << "\n---> RMC01AnalysisManager::book(): cannot open " << fFileName[1]
604  << G4endl;
605  return;
606  }
607 
608  // Create directories
609  theHistoManager->SetHistoDirectoryName("histo");
610 
611  //Histograms for :
612  // 1)the forward simulation results
613  // 2)the Reverse MC simulation results normalised to a user spectrum
614  //------------------------------------------------------------------------
615 
616 
617  G4int idHisto =
618  theHistoManager->CreateH1(G4String("Edep_vs_prim_ekin"),
619  G4String("edep vs e- primary energy"),60,emin,emax,
620  "none","none",G4String("log"));
621  fEdep_vs_prim_ekin = theHistoManager->GetH1(idHisto);
622 
623  idHisto = theHistoManager->CreateH1(G4String("elecron_current"),
624  G4String("electron"),60,emin,emax,
625  "none","none",G4String("log"));
626 
627  fElectron_current = theHistoManager->GetH1(idHisto);
628 
629  idHisto= theHistoManager->CreateH1(G4String("proton_current"),
630  G4String("proton"),60,emin,emax,
631  "none","none",G4String("log"));
632  fProton_current=theHistoManager->GetH1(idHisto);
633 
634 
635  idHisto= theHistoManager->CreateH1(G4String("gamma_current"),
636  G4String("gamma"),60,emin,emax,
637  "none","none",G4String("log"));
638  fGamma_current=theHistoManager->GetH1(idHisto);
639 
640 
641  //Response matrices for the adjoint simulation only
642  //-----------------------------------------------
643  if (fAdjoint_sim_mode){
644  //Response matrices for external isotropic e- source
645  //--------------------------------------------------
646 
647  idHisto =
648  theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_electron_prim_energy"),
649  G4String("electron RM vs e- primary energy"),60,emin,emax,
650  "none","none",G4String("log"));
651  fEdep_rmatrix_vs_electron_prim_energy = theHistoManager->GetH1(idHisto);
652 
653  idHisto =
654  theHistoManager->
655  CreateH2(G4String("Electron_current_rmatrix_vs_electron_prim_energy"),
656  G4String("electron current RM vs e- primary energy"),
657  60,emin,emax,60,emin,emax,
658  "none","none","none","none",G4String("log"),G4String("log"));
659 
660  fElectron_current_rmatrix_vs_electron_prim_energy =
661  theHistoManager->GetH2(idHisto);
662 
663  idHisto =
664  theHistoManager->
665  CreateH2(G4String("Gamma_current_rmatrix_vs_electron_prim_energy"),
666  G4String("gamma current RM vs e- primary energy"),
667  60,emin,emax,60,emin,emax,
668  "none","none","none","none",G4String("log"),G4String("log"));
669 
670 
671  fGamma_current_rmatrix_vs_electron_prim_energy =
672  theHistoManager->GetH2(idHisto);
673 
674 
675  //Response matrices for external isotropic gamma source
676 
677  idHisto =
678  theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_gamma_prim_energy"),
679  G4String("electron RM vs gamma primary energy"),60,emin,emax,
680  "none","none",G4String("log"));
681  fEdep_rmatrix_vs_gamma_prim_energy = theHistoManager->GetH1(idHisto);
682 
683  idHisto =
684  theHistoManager->
685  CreateH2(G4String("Electron_current_rmatrix_vs_gamma_prim_energy"),
686  G4String("electron current RM vs gamma primary energy"),
687  60,emin,emax,60,emin,emax,
688  "none","none","none","none",G4String("log"),G4String("log"));
689 
690  fElectron_current_rmatrix_vs_gamma_prim_energy =
691  theHistoManager->GetH2(idHisto);
692 
693  idHisto =
694  theHistoManager->
695  CreateH2(G4String("Gamma_current_rmatrix_vs_gamma_prim_energy"),
696  G4String("gamma current RM vs gamma primary energy"),
697  60,emin,emax,60,emin,emax,
698  "none","none","none","none",G4String("log"),G4String("log"));
699 
700  fGamma_current_rmatrix_vs_gamma_prim_energy =
701  theHistoManager->GetH2(idHisto);
702 
703 
704 
705  //Response matrices for external isotropic proton source
706  idHisto =
707  theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_proton_prim_energy"),
708  G4String("electron RM vs proton primary energy"),60,emin,emax,
709  "none","none",G4String("log"));
710  fEdep_rmatrix_vs_proton_prim_energy = theHistoManager->GetH1(idHisto);
711 
712  idHisto =
713  theHistoManager->
714  CreateH2(G4String("Electron_current_rmatrix_vs_proton_prim_energy"),
715  G4String("electron current RM vs proton primary energy"),
716  60,emin,emax,60,emin,emax,
717  "none","none","none","none",G4String("log"),G4String("log"));
718 
719  fElectron_current_rmatrix_vs_proton_prim_energy =
720  theHistoManager->GetH2(idHisto);
721 
722  idHisto =
723  theHistoManager->
724  CreateH2(G4String("Gamma_current_rmatrix_vs_proton_prim_energy"),
725  G4String("gamma current RM vs proton primary energy"),
726  60,emin,emax,60,emin,emax,
727  "none","none","none","none",G4String("log"),G4String("log"));
728 
729  fGamma_current_rmatrix_vs_proton_prim_energy =
730  theHistoManager->GetH2(idHisto);
731 
732  idHisto =
733  theHistoManager->
734  CreateH2(G4String("Proton_current_rmatrix_vs_proton_prim_energy"),
735  G4String("proton current RM vs proton primary energy"),
736  60,emin,emax,60,emin,emax,
737  "none","none","none","none",G4String("log"),G4String("log"));
738 
739  fProton_current_rmatrix_vs_proton_prim_energy =
740  theHistoManager->GetH2(idHisto);
741  }
742  fFactoryOn = true;
743  G4cout << "\n----> Histogram Tree is opened in " << fFileName[1] << G4endl;
744 }
G4VAnalysisManager::SetHistoDirectoryName
G4bool SetHistoDirectoryName(const G4String &dirName)
Definition: G4VAnalysisManager.cc:185
G4VAnalysisManager::SetFirstHistoId
G4bool SetFirstHistoId(G4int firstId)
Definition: G4VAnalysisManager.cc:401
G4VAnalysisManager::CreateH1
G4int CreateH1(const G4String &name, const G4String &title, G4int nbins, G4double xmin, G4double xmax, const G4String &unitName="none", const G4String &fcnName="none", const G4String &binSchemeName="linear")
Definition: G4VAnalysisManager.cc:215
G4int
int G4int
Definition: G4Types.hh:78
G4cout
G4GLOB_DLL std::ostream G4cout
G4bool
bool G4bool
Definition: G4Types.hh:79
G4VAnalysisManager::GetFileType
G4String GetFileType() const
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4RootAnalysisManager::GetH1
tools::histo::h1d * GetH1(G4int id, G4bool warn=true, G4bool onlyIfActive=true) const
G4RootAnalysisManager::GetH2
tools::histo::h2d * GetH2(G4int id, G4bool warn=true, G4bool onlyIfActive=true) const
void RMC01AnalysisManager::EndOfEvent ( const G4Event *  anEvent)

Definition at line 183 of file RMC01AnalysisManager.cc.

References G4RunManager::AbortRun(), G4cout, G4Event::GetEventID(), G4Timer::GetRealElapsed(), G4RunManager::GetRunManager(), G4Timer::Start(), and G4Timer::Stop().

Referenced by RMC01EventAction::EndOfEventAction().

184 {
185  if (fAdjoint_sim_mode) EndOfEventForAdjointSimulation(anEvent);
186  else EndOfEventForForwardSimulation(anEvent);
187 
188  //Test convergence. The error is already computed
189  //--------------------------------------
190  G4int nb_event=anEvent->GetEventID()+1;
191  //G4double factor=1.;
192  if (fAdjoint_sim_mode) {
193  G4double n_adj_evt= nb_event/fNb_evt_per_adj_evt;
194  // nb_event/fNb_evt_per_adj_evt;
195  if (n_adj_evt*fNb_evt_per_adj_evt == nb_event) {
196  nb_event =static_cast<G4int>(n_adj_evt);
197  }
198  else nb_event=0;
199 
200  }
201 
202  if (nb_event>100 && fStop_run_if_precision_reached &&
203  fPrecision_to_reach >fRelative_error) {
204  G4cout<<fPrecision_to_reach*100.<<"% Precision reached!"<<std::endl;
205  fTimer->Stop();
206  fElapsed_time+=fTimer->GetRealElapsed();
207  fConvergenceFileOutput<<fMean_edep<<'\t'<<fError_mean_edep
208  <<'\t'<<fElapsed_time<<std::endl;
209  G4RunManager::GetRunManager()->AbortRun(true);
210  }
211 
212  if (nb_event>0 && nb_event % fNb_evt_modulo_for_convergence_test == 0) {
213  fTimer->Stop();
214  fElapsed_time+=fTimer->GetRealElapsed();
215  fTimer->Start();
216  fConvergenceFileOutput<<fMean_edep<<'\t'<<fError_mean_edep<<'\t'
217  <<fElapsed_time<<std::endl;
218  }
219 }
G4RunManager::AbortRun
virtual void AbortRun(G4bool softAbort=false)
Definition: G4RunManager.cc:551
G4int
int G4int
Definition: G4Types.hh:78
G4Event::GetEventID
G4int GetEventID() const
Definition: G4Event.hh:140
G4cout
G4GLOB_DLL std::ostream G4cout
G4RunManager::GetRunManager
static G4RunManager * GetRunManager()
Definition: G4RunManager.cc:74
G4Timer::GetRealElapsed
G4double GetRealElapsed() const
Definition: G4Timer.cc:107
G4Timer::Stop
void Stop()
G4Timer::Start
void Start()
G4double
double G4double
Definition: G4Types.hh:76
void RMC01AnalysisManager::EndOfRun ( const G4Run *  aRun)

Definition at line 154 of file RMC01AnalysisManager.cc.

References G4cout, G4AdjointSimManager::GetInstance(), G4AdjointSimManager::GetNbEvtOfLastRun(), G4Run::GetNumberOfEvent(), save(), and G4Timer::Stop().

Referenced by RMC01RunAction::EndOfRunAction().

155 { fTimer->Stop();
156  G4int nb_evt=aRun->GetNumberOfEvent();
157  G4double factor =1./ nb_evt;
158  if (!fAdjoint_sim_mode){
159  G4cout<<"Results of forward simulation!"<<std::endl;
160  G4cout<<"edep per event [MeV] = "<<fMean_edep<<std::endl;
161  G4cout<<"error[MeV] = "<<fError_mean_edep<<std::endl;
162  }
163 
164  else {
165  G4cout<<"Results of reverse/adjoint simulation!"<<std::endl;
166  G4cout<<"normalised edep [MeV] = "<<fMean_edep<<std::endl;
167  G4cout<<"error[MeV] = "<<fError_mean_edep<<std::endl;
168  factor=1.*G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun()
169  *fNb_evt_per_adj_evt/aRun->GetNumberOfEvent();
170  }
171  save(factor);
172  fConvergenceFileOutput.close();
173 }
G4AdjointSimManager::GetInstance
static G4AdjointSimManager * GetInstance()
Definition: G4AdjointSimManager.cc:129
G4int
int G4int
Definition: G4Types.hh:78
RMC01AnalysisManager::save
void save(G4double scaling_factor)
Definition: RMC01AnalysisManager.cc:748
G4cout
G4GLOB_DLL std::ostream G4cout
G4Run::GetNumberOfEvent
G4int GetNumberOfEvent() const
Definition: G4Run.hh:79
G4Timer::Stop
void Stop()
G4double
double G4double
Definition: G4Types.hh:76
RMC01AnalysisManager * RMC01AnalysisManager::GetInstance ( void  )
static

Definition at line 109 of file RMC01AnalysisManager.cc.

Referenced by RMC01EventAction::BeginOfEventAction(), RMC01EventAction::EndOfEventAction(), and RMC01RunAction::RMC01RunAction().

110 {
111  if (fInstance == 0) fInstance = new RMC01AnalysisManager;
112  return fInstance;
113 }
void RMC01AnalysisManager::save ( G4double  scaling_factor)

Definition at line 748 of file RMC01AnalysisManager.cc.

References G4VAnalysisManager::CloseFile(), G4cout, G4endl, G4VAnalysisManager::GetNofH1s(), G4VAnalysisManager::GetNofH2s(), G4VAnalysisManager::ScaleH1(), G4VAnalysisManager::ScaleH2(), G4VAnalysisManager::SetH1Ascii(), and G4VAnalysisManager::Write().

Referenced by EndOfRun().

749 { if (fFactoryOn) {
750  G4AnalysisManager* theHistoManager = G4AnalysisManager::Instance();
751  //scaling of results
752  //-----------------
753 
754  for (int ind=1; ind<=theHistoManager->GetNofH1s();ind++){
755  theHistoManager->SetH1Ascii(ind,true);
756  theHistoManager->ScaleH1(ind,scaling_factor);
757  }
758  for (int ind=1; ind<=theHistoManager->GetNofH2s();ind++)
759  theHistoManager->ScaleH2(ind,scaling_factor);
760 
761  theHistoManager->Write();
762  theHistoManager->CloseFile();
763  G4cout << "\n----> Histogram Tree is saved in " << fFileName[1] << G4endl;
764 
765  delete G4AnalysisManager::Instance();
766  fFactoryOn = false;
767  }
768 }
G4VAnalysisManager::ScaleH1
G4bool ScaleH1(G4int id, G4double factor)
Definition: G4VAnalysisManager.cc:332
G4VAnalysisManager::ScaleH2
G4bool ScaleH2(G4int id, G4double factor)
Definition: G4VAnalysisManager.cc:338
G4VAnalysisManager::GetNofH1s
G4int GetNofH1s() const
Definition: G4VAnalysisManager.cc:487
G4VAnalysisManager::SetH1Ascii
void SetH1Ascii(G4int id, G4bool ascii)
Definition: G4VAnalysisManager.cc:523
G4cout
G4GLOB_DLL std::ostream G4cout
G4endl
#define G4endl
Definition: G4ios.hh:61
G4VAnalysisManager::GetNofH2s
G4int GetNofH2s() const
Definition: G4VAnalysisManager.cc:493
void RMC01AnalysisManager::SetPrecision ( G4double  precision)
inline

Definition at line 97 of file RMC01AnalysisManager.hh.

Referenced by RMC01AnalysisManagerMessenger::SetNewValue().

97  {
98  fPrecision_to_reach =precision/100.;};
void RMC01AnalysisManager::SetPrimaryExpSpectrumForAdjointSim ( const G4String &  particle_name,
G4double  fluence,
G4double  E0,
G4double  Emin,
G4double  Emax 
)

Definition at line 521 of file RMC01AnalysisManager.cc.

References G4Electron::Electron(), EXPO, G4cout, G4endl, G4Gamma::Gamma(), G4ParticleDefinition::GetPDGEncoding(), python.hepunit::pi, and G4Proton::Proton().

Referenced by RMC01AnalysisManagerMessenger::SetNewValue().

525 { fPrimSpectrumType = EXPO;
526  if (particle_name == "e-" ) fPrimPDG_ID =
527  G4Electron::Electron()->GetPDGEncoding();
528  else if (particle_name == "gamma") fPrimPDG_ID =
529  G4Gamma::Gamma()->GetPDGEncoding();
530  else if (particle_name == "proton") fPrimPDG_ID =
531  G4Proton::Proton()->GetPDGEncoding();
532  else {
533  G4cout<<"The particle that you did select is not in the candidate "<<
534  "list for primary [e-, gamma, proton]!"<<G4endl;
535  return;
536  }
537  fAlpha_or_E0 = E0 ;
538  fAmplitude_prim_spectrum = omni_fluence/E0/
539  (std::exp(-Emin/E0)-std::exp(-Emax/E0))/4./pi;
540  fEmin_prim_spectrum = Emin ;
541  fEmax_prim_spectrum = Emax;
542 }
EXPO
Definition: RMC01AnalysisManager.hh:71
G4cout
G4GLOB_DLL std::ostream G4cout
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4endl
#define G4endl
Definition: G4ios.hh:61
void RMC01AnalysisManager::SetPrimaryPowerLawSpectrumForAdjointSim ( const G4String &  particle_name,
G4double  fluence,
G4double  alpha,
G4double  Emin,
G4double  Emax 
)

Definition at line 546 of file RMC01AnalysisManager.cc.

References G4Electron::Electron(), G4cout, G4endl, G4Gamma::Gamma(), G4ParticleDefinition::GetPDGEncoding(), python.hepunit::pi, POWER, and G4Proton::Proton().

Referenced by RMC01AnalysisManagerMessenger::SetNewValue().

549 { fPrimSpectrumType =POWER;
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 
562 
563  if (alpha ==1.) {
564  fAmplitude_prim_spectrum = omni_fluence/std::log(Emax/Emin)/4./pi;
565  }
566  else {
567  G4double p=1.-alpha;
568  fAmplitude_prim_spectrum = omni_fluence/p/(std::pow(Emax,p)
569  -std::pow(Emin,p))/4./pi;
570  }
571 
572  fAlpha_or_E0 = alpha ;
573  fEmin_prim_spectrum = Emin ;
574  fEmax_prim_spectrum = Emax;
575 
576 }
POWER
Definition: RMC01AnalysisManager.hh:71
p
const char * p
Definition: xmltok.h:285
G4cout
G4GLOB_DLL std::ostream G4cout
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
The documentation for this class was generated from the following files:
Generated on Wed Apr 30 2014 15:57:30 for Geant4.10 by   doxygen 1.8.7