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 102 of file RMC01AnalysisManager.cc.

{;
}

Member Function Documentation

void RMC01AnalysisManager::BeginOfEvent

(

const G4Event *

)

Definition at line 176 of file RMC01AnalysisManager.cc.

{ ;
}

Here is the caller graph for this function:

void RMC01AnalysisManager::BeginOfRun

(

const G4Run *

aRun

)

Definition at line 116 of file RMC01AnalysisManager.cc.

{

   fAccumulated_edep =0.;
   fAccumulated_edep2 =0.;
   fNentry = 0.0;
   fRelative_error=1.;
   fMean_edep=0.;
   fError_mean_edep=0.;
   fAdjoint_sim_mode =G4AdjointSimManager::GetInstance()->GetAdjointSimMode();

   if (fAdjoint_sim_mode){
           fNb_evt_per_adj_evt=aRun->GetNumberOfEventToBeProcessed()/
                       G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
    fConvergenceFileOutput.open("ConvergenceOfAdjointSimulationResults.txt",
                                                                std::ios::out);
    fConvergenceFileOutput<<
           "Normalised Edep[MeV]\terror[MeV]\tcomputing_time[s]"<<std::endl;
   }
   else {
     fConvergenceFileOutput.open("ConvergenceOfForwardSimulationResults.txt",
                                                               std::ios::out);
     fConvergenceFileOutput<<
         "Edep per event [MeV]\terror[MeV]\tcomputing_time[s]"
                                                                 <<std::endl;
   }
   fConvergenceFileOutput.setf(std::ios::scientific);
   fConvergenceFileOutput.precision(6);         

   fTimer->Start();
   fElapsed_time=0.;

   Book();
}

Here is the call graph for this function:

Here is the caller graph for this function:

void RMC01AnalysisManager::Book

(

)

Definition at line 601 of file RMC01AnalysisManager.cc.

{
  //----------------------
  //Creation of histograms
  //----------------------

  //Energy binning of the histograms : 60 log bins over [1keV-1GeV]

  G4double emin=1.*keV;
  G4double emax=1.*GeV;

  //file_name
  fFileName[0]="forward_sim";
  if (fAdjoint_sim_mode) fFileName[0]="adjoint_sim";

  //Histo manager
   G4AnalysisManager* theHistoManager = G4AnalysisManager::Instance();
   G4String extension = theHistoManager->GetFileType();
   fFileName[1] = fFileName[0] + "." + extension;
   theHistoManager->SetFirstHistoId(1);

   G4bool fileOpen = theHistoManager->OpenFile(fFileName[0]);
   if (!fileOpen) {
     G4cout << "\n---> RMC01AnalysisManager::Book(): cannot open "
             << fFileName[1]
              << G4endl;
     return;
   }

    // Create directories
   theHistoManager->SetHistoDirectoryName("histo");

  //Histograms for :
  //        1)the forward simulation results
  //        2)the Reverse MC  simulation results normalised to a user spectrum
  //------------------------------------------------------------------------

   G4int idHisto =
      theHistoManager->CreateH1(G4String("Edep_vs_prim_ekin"),
      G4String("edep vs e- primary energy"),60,emin,emax,
      "none","none",G4String("log"));
  fEdep_vs_prim_ekin = theHistoManager->GetH1(idHisto);

  idHisto = theHistoManager->CreateH1(G4String("elecron_current"),
        G4String("electron"),60,emin,emax,
        "none","none",G4String("log"));

  fElectron_current  =  theHistoManager->GetH1(idHisto);

  idHisto= theHistoManager->CreateH1(G4String("proton_current"),
        G4String("proton"),60,emin,emax,
        "none","none",G4String("log"));
  fProton_current=theHistoManager->GetH1(idHisto);

  idHisto= theHistoManager->CreateH1(G4String("gamma_current"),
          G4String("gamma"),60,emin,emax,
          "none","none",G4String("log"));
  fGamma_current=theHistoManager->GetH1(idHisto);

  //Response matrices for the adjoint simulation only
  //-----------------------------------------------
  if (fAdjoint_sim_mode){
  //Response matrices for external isotropic e- source
  //--------------------------------------------------

   idHisto =
    theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_electron_prim_energy"),
    G4String("electron RM vs e- primary energy"),60,emin,emax,
    "none","none",G4String("log"));
   fEdep_rmatrix_vs_electron_prim_energy = theHistoManager->GetH1(idHisto);

   idHisto =
      theHistoManager->
         CreateH2(G4String("Electron_current_rmatrix_vs_electron_prim_energy"),
         G4String("electron current  RM vs e- primary energy"),
         60,emin,emax,60,emin,emax,
         "none","none","none","none",G4String("log"),G4String("log"));

   fElectron_current_rmatrix_vs_electron_prim_energy =
                                             theHistoManager->GetH2(idHisto);

   idHisto =
        theHistoManager->
           CreateH2(G4String("Gamma_current_rmatrix_vs_electron_prim_energy"),
           G4String("gamma current  RM vs e- primary energy"),
           60,emin,emax,60,emin,emax,
           "none","none","none","none",G4String("log"),G4String("log"));

   fGamma_current_rmatrix_vs_electron_prim_energy =
                                           theHistoManager->GetH2(idHisto);

  //Response matrices for external isotropic gamma source

   idHisto =
    theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_gamma_prim_energy"),
         G4String("electron RM vs gamma primary energy"),60,emin,emax,
        "none","none",G4String("log"));
   fEdep_rmatrix_vs_gamma_prim_energy = theHistoManager->GetH1(idHisto);

   idHisto =
        theHistoManager->
          CreateH2(G4String("Electron_current_rmatrix_vs_gamma_prim_energy"),
          G4String("electron current  RM vs gamma primary energy"),
          60,emin,emax,60,emin,emax,
        "none","none","none","none",G4String("log"),G4String("log"));

   fElectron_current_rmatrix_vs_gamma_prim_energy =
                                               theHistoManager->GetH2(idHisto);

   idHisto =
          theHistoManager->
             CreateH2(G4String("Gamma_current_rmatrix_vs_gamma_prim_energy"),
             G4String("gamma current  RM vs gamma primary energy"),
             60,emin,emax,60,emin,emax,
             "none","none","none","none",G4String("log"),G4String("log"));

   fGamma_current_rmatrix_vs_gamma_prim_energy =
                                             theHistoManager->GetH2(idHisto);

  //Response matrices for external isotropic proton source
   idHisto =
     theHistoManager->CreateH1(G4String("Edep_rmatrix_vs_proton_prim_energy"),
         G4String("electron RM vs proton primary energy"),60,emin,emax,
         "none","none",G4String("log"));
   fEdep_rmatrix_vs_proton_prim_energy = theHistoManager->GetH1(idHisto);

   idHisto =
         theHistoManager->
           CreateH2(G4String("Electron_current_rmatrix_vs_proton_prim_energy"),
           G4String("electron current  RM vs proton primary energy"),
           60,emin,emax,60,emin,emax,
         "none","none","none","none",G4String("log"),G4String("log"));

    fElectron_current_rmatrix_vs_proton_prim_energy =
                                                theHistoManager->GetH2(idHisto);

   idHisto =
           theHistoManager->
              CreateH2(G4String("Gamma_current_rmatrix_vs_proton_prim_energy"),
              G4String("gamma current  RM vs proton primary energy"),
              60,emin,emax,60,emin,emax,
              "none","none","none","none",G4String("log"),G4String("log"));

   fGamma_current_rmatrix_vs_proton_prim_energy =
                                              theHistoManager->GetH2(idHisto);

   idHisto =
              theHistoManager->
              CreateH2(G4String("Proton_current_rmatrix_vs_proton_prim_energy"),
               G4String("proton current  RM vs proton primary energy"),
               60,emin,emax,60,emin,emax,
               "none","none","none","none",G4String("log"),G4String("log"));

      fProton_current_rmatrix_vs_proton_prim_energy =
                                                theHistoManager->GetH2(idHisto);
  }
  fFactoryOn = true;
  G4cout << "\n----> Histogram Tree is opened in " << fFileName[1] << G4endl;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void RMC01AnalysisManager::EndOfEvent

(

const G4Event *

anEvent

)

Definition at line 182 of file RMC01AnalysisManager.cc.

{  
   if (fAdjoint_sim_mode) EndOfEventForAdjointSimulation(anEvent);
   else EndOfEventForForwardSimulation(anEvent);

   //Test convergence. The error is already computed
   //--------------------------------------
   G4int nb_event=anEvent->GetEventID()+1;
   //G4double factor=1.;
   if (fAdjoint_sim_mode) {
           G4double  n_adj_evt= nb_event/fNb_evt_per_adj_evt;
        // nb_event/fNb_evt_per_adj_evt;
        if (n_adj_evt*fNb_evt_per_adj_evt == nb_event) {
                nb_event =static_cast<G4int>(n_adj_evt);
        }        
        else nb_event=0;
   }
   
   if (nb_event>100 && fStop_run_if_precision_reached &&
                                      fPrecision_to_reach >fRelative_error) {
      G4cout<<fPrecision_to_reach*100.<<"%  Precision reached!"<<std::endl;
      fTimer->Stop();
      fElapsed_time+=fTimer->GetRealElapsed();
      fConvergenceFileOutput<<fMean_edep<<'\t'<<fError_mean_edep
                                         <<'\t'<<fElapsed_time<<std::endl;
      G4RunManager::GetRunManager()->AbortRun(true);
   }
   
   if (nb_event>0 && nb_event % fNb_evt_modulo_for_convergence_test == 0) {
      fTimer->Stop();
      fElapsed_time+=fTimer->GetRealElapsed();
      fTimer->Start();
      fConvergenceFileOutput<<fMean_edep<<'\t'<<fError_mean_edep<<'\t'
                                                   <<fElapsed_time<<std::endl;
   }
}   

Here is the call graph for this function:

Here is the caller graph for this function:

void RMC01AnalysisManager::EndOfRun

(

const G4Run *

aRun

)

Definition at line 153 of file RMC01AnalysisManager.cc.

{ fTimer->Stop();
  G4int nb_evt=aRun->GetNumberOfEvent();
  G4double factor =1./ nb_evt;
  if (!fAdjoint_sim_mode){
   G4cout<<"Results of forward simulation!"<<std::endl;
   G4cout<<"edep per event [MeV] = "<<fMean_edep<<std::endl;
   G4cout<<"error[MeV] = "<<fError_mean_edep<<std::endl;
  }
  
  else {
   G4cout<<"Results of reverse/adjoint simulation!"<<std::endl;
   G4cout<<"normalised edep [MeV] = "<<fMean_edep<<std::endl;
   G4cout<<"error[MeV] = "<<fError_mean_edep<<std::endl;
   factor=1.*G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun()
                                 *fNb_evt_per_adj_evt/aRun->GetNumberOfEvent();
  }
  Save(factor);
  fConvergenceFileOutput.close();
}

Here is the call graph for this function:

Here is the caller graph for this function:

RMC01AnalysisManager * RMC01AnalysisManager::GetInstance ( void  )

static

Definition at line 108 of file RMC01AnalysisManager.cc.

{
  if (fInstance == 0) fInstance = new RMC01AnalysisManager;
  return fInstance;
}

Here is the caller graph for this function:

void RMC01AnalysisManager::Save

(

G4double

scaling_factor

)

Definition at line 763 of file RMC01AnalysisManager.cc.

{ if (fFactoryOn) {
    G4AnalysisManager* theHistoManager = G4AnalysisManager::Instance();
    //scaling of results
    //-----------------

    for (int ind=1; ind<=theHistoManager->GetNofH1s();ind++){
       theHistoManager->SetH1Ascii(ind,true);
       theHistoManager->ScaleH1(ind,scaling_factor);
    }
    for (int ind=1; ind<=theHistoManager->GetNofH2s();ind++)
                        theHistoManager->ScaleH2(ind,scaling_factor);

    theHistoManager->Write();
    theHistoManager->CloseFile();
    G4cout << "\n----> Histogram Tree is saved in " << fFileName[1] << G4endl;

    delete G4AnalysisManager::Instance();
    fFactoryOn = false;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void RMC01AnalysisManager::SetPrecision ( G4double  precision )

inline

Definition at line 97 of file RMC01AnalysisManager.hh.

                                               {
                                   fPrecision_to_reach =precision/100.;};

Here is the caller graph for this function:

void RMC01AnalysisManager::SetPrimaryExpSpectrumForAdjointSim	(	const G4String &	particle_name,
		G4double	fluence,
		G4double	E0,
		G4double	Emin,
		G4double	Emax
	)

Definition at line 542 of file RMC01AnalysisManager.cc.

{ fPrimSpectrumType = EXPO;
  if (particle_name == "e-" ) fPrimPDG_ID =
                      G4Electron::Electron()->GetPDGEncoding();
  else if (particle_name == "gamma") fPrimPDG_ID =
                                       G4Gamma::Gamma()->GetPDGEncoding();
  else if (particle_name == "proton") fPrimPDG_ID =
                                      G4Proton::Proton()->GetPDGEncoding();
  else {
   G4cout<<"The particle that you did select is not in the candidate "<<
       "list for primary [e-, gamma, proton]!"<<G4endl;
          return;
  }        
  fAlpha_or_E0 = E0 ;
  fAmplitude_prim_spectrum = omni_fluence/E0/
                              (std::exp(-Emin/E0)-std::exp(-Emax/E0))/4./pi;
  fEmin_prim_spectrum = Emin ;
  fEmax_prim_spectrum = Emax;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void RMC01AnalysisManager::SetPrimaryPowerLawSpectrumForAdjointSim	(	const G4String &	particle_name,
		G4double	fluence,
		G4double	alpha,
		G4double	Emin,
		G4double	Emax
	)

Definition at line 567 of file RMC01AnalysisManager.cc.

{ fPrimSpectrumType  =POWER;
  if (particle_name == "e-" ) fPrimPDG_ID =
                              G4Electron::Electron()->GetPDGEncoding();
  else if (particle_name == "gamma") fPrimPDG_ID =
                              G4Gamma::Gamma()->GetPDGEncoding();
  else if (particle_name == "proton") fPrimPDG_ID =
                              G4Proton::Proton()->GetPDGEncoding();
  else {
          G4cout<<"The particle that you did select is not in the candidate"<<
          " list for primary [e-, gamma, proton]!"<<G4endl;
          return;
  }        
  

 if (alpha ==1.) {
         fAmplitude_prim_spectrum = omni_fluence/std::log(Emax/Emin)/4./pi;
 }
 else {
         G4double p=1.-alpha;
         fAmplitude_prim_spectrum = omni_fluence/p/(std::pow(Emax,p)
                                                   -std::pow(Emin,p))/4./pi;
 }

  fAlpha_or_E0 = alpha ;
  fEmin_prim_spectrum = Emin ;
  fEmax_prim_spectrum = Emax;

}

Here is the call graph for this function:

Here is the caller graph for this function:

---

The documentation for this class was generated from the following files:

• source/geant4.10.03.p03/examples/extended/biasing/ReverseMC01/include/RMC01AnalysisManager.hh
• source/geant4.10.03.p03/examples/extended/biasing/ReverseMC01/src/RMC01AnalysisManager.cc