RMC01AnalysisManager Class Reference

#include <RMC01AnalysisManager.hh>

Collaboration diagram for RMC01AnalysisManager:

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 ()

Private Member Functions
	RMC01AnalysisManager ()
void	EndOfEventForForwardSimulation (const G4Event *anEvent)
void	EndOfEventForAdjointSimulation (const G4Event *anEvent)
G4double	PrimDiffAndDirFluxForAdjointSim (G4double prim_energy)
void	ComputeMeanEdepAndError (const G4Event *anEvent, G4double &mean, G4double &error)

Private Attributes
RMC01AnalysisManagerMessenger *	fMsg
G4AnaH1 *	fEdep_vs_prim_ekin
G4AnaH1 *	fElectron_current
G4AnaH1 *	fProton_current
G4AnaH1 *	fGamma_current
G4double	fAccumulated_edep
G4double	fAccumulated_edep2
G4double	fNentry
G4double	fMean_edep
G4double	fError_mean_edep
G4double	fRelative_error
G4double	fElapsed_time
G4double	fPrecision_to_reach
G4bool	fStop_run_if_precision_reached
G4int	fNb_evt_modulo_for_convergence_test
G4AnaH1 *	fEdep_rmatrix_vs_electron_prim_energy
G4AnaH2 *	fElectron_current_rmatrix_vs_electron_prim_energy
G4AnaH2 *	fGamma_current_rmatrix_vs_electron_prim_energy
G4AnaH1 *	fEdep_rmatrix_vs_gamma_prim_energy
G4AnaH2 *	fElectron_current_rmatrix_vs_gamma_prim_energy
G4AnaH2 *	fGamma_current_rmatrix_vs_gamma_prim_energy
G4AnaH1 *	fEdep_rmatrix_vs_proton_prim_energy
G4AnaH2 *	fElectron_current_rmatrix_vs_proton_prim_energy
G4AnaH2 *	fProton_current_rmatrix_vs_proton_prim_energy
G4AnaH2 *	fGamma_current_rmatrix_vs_proton_prim_energy
G4String	fFileName [2]
G4bool	fFactoryOn
PRIM_SPECTRUM_TYPE	fPrimSpectrumType
G4int	fPrimPDG_ID
G4double	fAlpha_or_E0
G4double	fAmplitude_prim_spectrum
G4double	fEmin_prim_spectrum
G4double	fEmax_prim_spectrum
G4bool	fAdjoint_sim_mode
G4int	fNb_evt_per_adj_evt
G4Timer *	fTimer
std::fstream	fConvergenceFileOutput

Static Private Attributes
static RMC01AnalysisManager *	fInstance = 0

Detailed Description

Definition at line 77 of file RMC01AnalysisManager.hh.

Constructor & Destructor Documentation

~RMC01AnalysisManager()

RMC01AnalysisManager::~RMC01AnalysisManager

(

)

Definition at line 103 of file RMC01AnalysisManager.cc.

{;
}

RMC01AnalysisManager()

RMC01AnalysisManager::RMC01AnalysisManager ( )

private

Definition at line 61 of file RMC01AnalysisManager.cc.

 :fAccumulated_edep(0.), fAccumulated_edep2(0.), fMean_edep(0.),
  fError_mean_edep(0.), fRelative_error(0.), fElapsed_time(0.),
  fPrecision_to_reach(0.),fStop_run_if_precision_reached(true),
  fNb_evt_modulo_for_convergence_test(5000),
  fEdep_rmatrix_vs_electron_prim_energy(0),
  fElectron_current_rmatrix_vs_electron_prim_energy(0),
  fGamma_current_rmatrix_vs_electron_prim_energy(0),
  fEdep_rmatrix_vs_gamma_prim_energy(0),
  fElectron_current_rmatrix_vs_gamma_prim_energy(0),
  fGamma_current_rmatrix_vs_gamma_prim_energy(0),
  fEdep_rmatrix_vs_proton_prim_energy(0),
  fElectron_current_rmatrix_vs_proton_prim_energy(0),
  fProton_current_rmatrix_vs_proton_prim_energy(0),
  fGamma_current_rmatrix_vs_proton_prim_energy(0),
  fFactoryOn(false),
  fPrimSpectrumType(EXPO),
  fAlpha_or_E0(.5*MeV),fAmplitude_prim_spectrum (1.),
  fEmin_prim_spectrum(1.*keV),fEmax_prim_spectrum (20.*MeV),
  fAdjoint_sim_mode(true),fNb_evt_per_adj_evt(2)
{ 
  
  fMsg = new RMC01AnalysisManagerMessenger(this);

  //-------------
  //Timer for convergence vector
  //-------------
  
  fTimer = new G4Timer();

  //---------------------------------
  //Primary particle ID for normalisation of adjoint results
  //---------------------------------
  
  fPrimPDG_ID = G4Electron::Electron()->GetPDGEncoding();
  
  fFileName[0] = "sim";

}

Here is the call graph for this function:

Here is the caller graph for this function:

Member Function Documentation

BeginOfEvent()

void RMC01AnalysisManager::BeginOfEvent

(

const G4Event *

)

Definition at line 178 of file RMC01AnalysisManager.cc.

{ ;
}

Here is the caller graph for this function:

BeginOfRun()

void RMC01AnalysisManager::BeginOfRun

(

const G4Run *

aRun

)

Definition at line 117 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:

book()

void RMC01AnalysisManager::book

(

)

Definition at line 589 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:

ComputeMeanEdepAndError()

void RMC01AnalysisManager::ComputeMeanEdepAndError ( const G4Event *  anEvent, G4double &  mean, G4double &  error  )

private

Definition at line 499 of file RMC01AnalysisManager.cc.

{  
   G4int nb_event=anEvent->GetEventID()+1;
   G4double factor=1.;
   if (fAdjoint_sim_mode) {
      nb_event /=fNb_evt_per_adj_evt;
      factor=1.*G4AdjointSimManager::GetInstance()->GetNbEvtOfLastRun();
   }
   
   // VI: error computation now is based on number of entries and not 
   //     number of events
   if (fNentry > 1.0) {
      mean = fAccumulated_edep/fNentry;
      G4double mean_x2 = fAccumulated_edep2/fNentry;
      /*
      G4cout << "Nevt= " << nb_event <<  " mean= " << mean 
             << "  mean_x2= " <<  mean_x2 << " x2 - x*x= " 
             << mean_x2-mean*mean << G4endl;
      */
      error = factor*std::sqrt(mean_x2-mean*mean)/std::sqrt(fNentry);
      mean *=factor;
   }
   else {
      mean=0;
      error=0;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

EndOfEvent()

void RMC01AnalysisManager::EndOfEvent

(

const G4Event *

anEvent

)

Definition at line 184 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:

EndOfEventForAdjointSimulation()

void RMC01AnalysisManager::EndOfEventForAdjointSimulation ( const G4Event *  anEvent )

private

Definition at line 292 of file RMC01AnalysisManager.cc.

{  
  //Output from Sensitive volume computed during the forward tracking phase
  //-----------------------------------------------------------------------
  G4SDManager* SDman = G4SDManager::GetSDMpointer();
  G4HCofThisEvent* HCE = anEvent->GetHCofThisEvent();
  RMC01DoubleWithWeightHitsCollection* edepCollection =
                  (RMC01DoubleWithWeightHitsCollection*)(
                    HCE->GetHC(SDman->GetCollectionID("edep")));

  RMC01DoubleWithWeightHitsCollection* electronCurrentCollection =
                  (RMC01DoubleWithWeightHitsCollection*)(
                 HCE->GetHC(SDman->GetCollectionID("current_electron")));

  RMC01DoubleWithWeightHitsCollection* protonCurrentCollection =
                  (RMC01DoubleWithWeightHitsCollection*)(
                 HCE->GetHC(SDman->GetCollectionID("current_proton")));

  RMC01DoubleWithWeightHitsCollection* gammaCurrentCollection =
                  (RMC01DoubleWithWeightHitsCollection*)(
                 HCE->GetHC(SDman->GetCollectionID("current_gamma")));
  
  //Output from adjoint tracking phase
  //----------------------------------------------------------------------------
  
  G4AdjointSimManager* theAdjointSimManager =
                   G4AdjointSimManager::GetInstance();
  G4int pdg_nb =theAdjointSimManager
           ->GetFwdParticlePDGEncodingAtEndOfLastAdjointTrack();
  G4double prim_ekin=theAdjointSimManager->GetEkinAtEndOfLastAdjointTrack();
  G4double adj_weight=theAdjointSimManager->GetWeightAtEndOfLastAdjointTrack();
 
  
  //Factor of normalisation to user defined prim spectrum (power law or exp)
  //---------------------------------------------------------------------------
  
  G4double normalised_weight = 0.;
  if (pdg_nb== fPrimPDG_ID && prim_ekin>= fEmin_prim_spectrum
                                         && prim_ekin<= fEmax_prim_spectrum)
  normalised_weight =
                adj_weight*PrimDiffAndDirFluxForAdjointSim(prim_ekin);
  
  //Answer matrices
  //-------------
   G4AnaH1* edep_rmatrix =0;
   G4AnaH2* electron_current_rmatrix =0;
   G4AnaH2* gamma_current_rmatrix =0;
   G4AnaH2* proton_current_rmatrix =0;

   if (pdg_nb == G4Electron::Electron()->GetPDGEncoding()){ //e- answer matrices
     edep_rmatrix = fEdep_rmatrix_vs_electron_prim_energy;

     electron_current_rmatrix =
                          fElectron_current_rmatrix_vs_electron_prim_energy;

     gamma_current_rmatrix = fGamma_current_rmatrix_vs_electron_prim_energy;
   }
   else if (pdg_nb == G4Gamma::Gamma()->GetPDGEncoding()){
     //gammma answer matrices
     edep_rmatrix = fEdep_rmatrix_vs_gamma_prim_energy;
     electron_current_rmatrix = fElectron_current_rmatrix_vs_gamma_prim_energy;
     gamma_current_rmatrix = fGamma_current_rmatrix_vs_gamma_prim_energy;
   }
   else if (pdg_nb == G4Proton::Proton()->GetPDGEncoding()){
     //proton answer matrices
     edep_rmatrix = fEdep_rmatrix_vs_proton_prim_energy;
     electron_current_rmatrix = fElectron_current_rmatrix_vs_proton_prim_energy;
     gamma_current_rmatrix = fGamma_current_rmatrix_vs_proton_prim_energy;
     proton_current_rmatrix = fProton_current_rmatrix_vs_proton_prim_energy;
   }
  
  //Registering of total energy deposited in Event
  //-------------------------------
   G4double totEdep=0; 
   G4int i;
   for (i=0;i<edepCollection->entries();i++)
           totEdep+=(*edepCollection)[i]->GetValue()*
                                            (*edepCollection)[i]->GetWeight();
   
   G4bool new_mean_computed=false;
   if (totEdep>0.){
     if (normalised_weight>0.){
       G4double edep=totEdep* normalised_weight;
                
      //Check if the edep is not wrongly too high
      //-----------------------------------------
      G4double new_mean , new_error;
                
      fAccumulated_edep +=edep;
      fAccumulated_edep2 +=edep*edep;
      fNentry += 1.0;
      ComputeMeanEdepAndError(anEvent,new_mean,new_error);
      G4double new_relative_error = 1.;
      if ( new_error >0) new_relative_error = new_error/ new_mean;
      if (fRelative_error <0.10 && new_relative_error>1.5*fRelative_error) {
         G4cout<<"Potential wrong adjoint weight!"<<std::endl;
         G4cout<<"The results of this event will not be registered!"
                                                            <<std::endl;
         G4cout<<"previous mean edep [MeV] "<< fMean_edep<<std::endl;
         G4cout<<"previous relative error "<< fRelative_error<<std::endl;
         G4cout<<"new rejected mean edep [MeV] "<< new_mean<<std::endl;
         G4cout<<"new rejected relative error "<< new_relative_error
                                                              <<std::endl;
         fAccumulated_edep -=edep;
         fAccumulated_edep2 -=edep*edep;
         fNentry -= 1.0;
         return; 
                }
                else { //accepted
                        fMean_edep = new_mean;
                        fError_mean_edep = new_error;
                        fRelative_error =new_relative_error;
                        new_mean_computed=true;
                }         
                fEdep_vs_prim_ekin->fill(prim_ekin,edep); 
        }
        
        // Registering answer matrix
        //---------------------------
        
        edep_rmatrix->fill(prim_ekin,totEdep*adj_weight/cm2);
   }
   if (!new_mean_computed){
     ComputeMeanEdepAndError(anEvent,fMean_edep,fError_mean_edep);
     if (fError_mean_edep>0) fRelative_error= fError_mean_edep/fMean_edep;
   }         
    
   
  //Registering of current of particles on the sensitive volume
  //------------------------------------------------------------
   
   for (i=0;i<electronCurrentCollection->entries();i++) {
     G4double ekin =(*electronCurrentCollection)[i]->GetValue();
     G4double weight=(*electronCurrentCollection)[i]->GetWeight();
     fElectron_current->fill(ekin,weight*normalised_weight);
     electron_current_rmatrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
   }
   
   for (i=0;i<protonCurrentCollection->entries();i++) {
     G4double ekin =(*protonCurrentCollection)[i]->GetValue();
     G4double weight=(*protonCurrentCollection)[i]->GetWeight();
     fProton_current->fill(ekin,weight*normalised_weight);
     proton_current_rmatrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
  }        
   
   for (i=0;i<gammaCurrentCollection->entries();i++) {
     G4double ekin =(*gammaCurrentCollection)[i]->GetValue();
     G4double weight=(*gammaCurrentCollection)[i]->GetWeight();
     fGamma_current->fill(ekin,weight*normalised_weight);
     gamma_current_rmatrix->fill(prim_ekin,ekin,weight*adj_weight/cm2);
   }
}

Here is the call graph for this function:

Here is the caller graph for this function:

EndOfEventForForwardSimulation()

void RMC01AnalysisManager::EndOfEventForForwardSimulation ( const G4Event *  anEvent )

private

Definition at line 224 of file RMC01AnalysisManager.cc.

{  
   
   G4SDManager* SDman = G4SDManager::GetSDMpointer();
   G4HCofThisEvent* HCE = anEvent->GetHCofThisEvent();
   RMC01DoubleWithWeightHitsCollection* edepCollection =
         (RMC01DoubleWithWeightHitsCollection*)
              (HCE->GetHC(SDman->GetCollectionID("edep")));

   RMC01DoubleWithWeightHitsCollection* electronCurrentCollection =
             (RMC01DoubleWithWeightHitsCollection*)
                (HCE->GetHC(SDman->GetCollectionID("current_electron")));

   RMC01DoubleWithWeightHitsCollection* protonCurrentCollection =
             (RMC01DoubleWithWeightHitsCollection*)
                   (HCE->GetHC(SDman->GetCollectionID("current_proton")));

   RMC01DoubleWithWeightHitsCollection* gammaCurrentCollection =
             (RMC01DoubleWithWeightHitsCollection*)
                     (HCE->GetHC(SDman->GetCollectionID("current_gamma")));
   
   //Total energy deposited in Event
   //-------------------------------
   G4double totEdep=0; 
   G4int i;
   for (i=0;i<edepCollection->entries();i++)
        totEdep+=(*edepCollection)[i]->GetValue()
                  *(*edepCollection)[i]->GetWeight();
   
   if (totEdep>0.){
           fAccumulated_edep +=totEdep ;
           fAccumulated_edep2 +=totEdep*totEdep;
           fNentry += 1.0;
           G4PrimaryParticle* thePrimary=anEvent->GetPrimaryVertex()->GetPrimary();
           G4double E0= thePrimary->GetG4code()->GetPDGMass();
           G4double P=thePrimary->GetMomentum().mag();
           G4double prim_ekin =std::sqrt(E0*E0+P*P)-E0;
           fEdep_vs_prim_ekin->fill(prim_ekin,totEdep);
   } 
   ComputeMeanEdepAndError(anEvent,fMean_edep,fError_mean_edep);
   if (fError_mean_edep>0) fRelative_error= fError_mean_edep/fMean_edep;
                   
   //Particle current on sensitive cylinder
   //-------------------------------------
   
   for (i=0;i<electronCurrentCollection->entries();i++) {
           G4double ekin =(*electronCurrentCollection)[i]->GetValue();
        G4double weight=(*electronCurrentCollection)[i]->GetWeight();
        fElectron_current->fill(ekin,weight);
   }
   
   for (i=0;i<protonCurrentCollection->entries();i++) {
           G4double ekin =(*protonCurrentCollection)[i]->GetValue();
        G4double weight=(*protonCurrentCollection)[i]->GetWeight();
        fProton_current->fill(ekin,weight);
   }        
   
   for (i=0;i<gammaCurrentCollection->entries();i++) {
           G4double ekin =(*gammaCurrentCollection)[i]->GetValue();
        G4double weight=(*gammaCurrentCollection)[i]->GetWeight();
        fGamma_current->fill(ekin,weight);
   }
   
}

Here is the call graph for this function:

Here is the caller graph for this function:

EndOfRun()

void RMC01AnalysisManager::EndOfRun

(

const G4Run *

aRun

)

Definition at line 155 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:

GetInstance()

RMC01AnalysisManager * RMC01AnalysisManager::GetInstance ( void  )

static

Definition at line 109 of file RMC01AnalysisManager.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

PrimDiffAndDirFluxForAdjointSim()

G4double RMC01AnalysisManager::PrimDiffAndDirFluxForAdjointSim ( G4double  prim_energy )

private

Definition at line 448 of file RMC01AnalysisManager.cc.

{ 
  G4double flux=fAmplitude_prim_spectrum;
  if ( fPrimSpectrumType ==EXPO)      flux*=std::exp(-prim_energy/fAlpha_or_E0);
  else flux*=std::pow(prim_energy, -fAlpha_or_E0);
  return flux;
}

Here is the caller graph for this function:

save()

void RMC01AnalysisManager::save

(

G4double

scaling_factor

)

Definition at line 757 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:

SetPrecision()

void RMC01AnalysisManager::SetPrecision ( G4double  precision )

inline

Definition at line 97 of file RMC01AnalysisManager.hh.

                                               {
                                   fPrecision_to_reach =precision/100.;};

Here is the call graph for this function:

Here is the caller graph for this function:

SetPrimaryExpSpectrumForAdjointSim()

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

Definition at line 530 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:

SetPrimaryPowerLawSpectrumForAdjointSim()

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

Definition at line 555 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:

Member Data Documentation

fAccumulated_edep

G4double RMC01AnalysisManager::fAccumulated_edep

private

Definition at line 138 of file RMC01AnalysisManager.hh.

fAccumulated_edep2

G4double RMC01AnalysisManager::fAccumulated_edep2

private

Definition at line 139 of file RMC01AnalysisManager.hh.

fAdjoint_sim_mode

G4bool RMC01AnalysisManager::fAdjoint_sim_mode

private

Definition at line 179 of file RMC01AnalysisManager.hh.

fAlpha_or_E0

G4double RMC01AnalysisManager::fAlpha_or_E0

private

Definition at line 175 of file RMC01AnalysisManager.hh.

fAmplitude_prim_spectrum

G4double RMC01AnalysisManager::fAmplitude_prim_spectrum

private

Definition at line 176 of file RMC01AnalysisManager.hh.

fConvergenceFileOutput

std::fstream RMC01AnalysisManager::fConvergenceFileOutput

private

Definition at line 185 of file RMC01AnalysisManager.hh.

fEdep_rmatrix_vs_electron_prim_energy

G4AnaH1* RMC01AnalysisManager::fEdep_rmatrix_vs_electron_prim_energy

private

Definition at line 152 of file RMC01AnalysisManager.hh.

fEdep_rmatrix_vs_gamma_prim_energy

G4AnaH1* RMC01AnalysisManager::fEdep_rmatrix_vs_gamma_prim_energy

private

Definition at line 156 of file RMC01AnalysisManager.hh.

fEdep_rmatrix_vs_proton_prim_energy

G4AnaH1* RMC01AnalysisManager::fEdep_rmatrix_vs_proton_prim_energy

private

Definition at line 160 of file RMC01AnalysisManager.hh.

fEdep_vs_prim_ekin

G4AnaH1* RMC01AnalysisManager::fEdep_vs_prim_ekin

private

Definition at line 126 of file RMC01AnalysisManager.hh.

fElapsed_time

G4double RMC01AnalysisManager::fElapsed_time

private

Definition at line 144 of file RMC01AnalysisManager.hh.

fElectron_current

G4AnaH1* RMC01AnalysisManager::fElectron_current

private

Definition at line 127 of file RMC01AnalysisManager.hh.

fElectron_current_rmatrix_vs_electron_prim_energy

G4AnaH2* RMC01AnalysisManager::fElectron_current_rmatrix_vs_electron_prim_energy

private

Definition at line 153 of file RMC01AnalysisManager.hh.

fElectron_current_rmatrix_vs_gamma_prim_energy

G4AnaH2* RMC01AnalysisManager::fElectron_current_rmatrix_vs_gamma_prim_energy

private

Definition at line 157 of file RMC01AnalysisManager.hh.

fElectron_current_rmatrix_vs_proton_prim_energy

G4AnaH2* RMC01AnalysisManager::fElectron_current_rmatrix_vs_proton_prim_energy

private

Definition at line 161 of file RMC01AnalysisManager.hh.

fEmax_prim_spectrum

G4double RMC01AnalysisManager::fEmax_prim_spectrum

private

Definition at line 178 of file RMC01AnalysisManager.hh.

fEmin_prim_spectrum

G4double RMC01AnalysisManager::fEmin_prim_spectrum

private

Definition at line 177 of file RMC01AnalysisManager.hh.

fError_mean_edep

G4double RMC01AnalysisManager::fError_mean_edep

private

Definition at line 142 of file RMC01AnalysisManager.hh.

fFactoryOn

G4bool RMC01AnalysisManager::fFactoryOn

private

Definition at line 166 of file RMC01AnalysisManager.hh.

fFileName

G4String RMC01AnalysisManager::fFileName[2]

private

Definition at line 165 of file RMC01AnalysisManager.hh.

fGamma_current

G4AnaH1* RMC01AnalysisManager::fGamma_current

private

Definition at line 129 of file RMC01AnalysisManager.hh.

fGamma_current_rmatrix_vs_electron_prim_energy

G4AnaH2* RMC01AnalysisManager::fGamma_current_rmatrix_vs_electron_prim_energy

private

Definition at line 154 of file RMC01AnalysisManager.hh.

fGamma_current_rmatrix_vs_gamma_prim_energy

G4AnaH2* RMC01AnalysisManager::fGamma_current_rmatrix_vs_gamma_prim_energy

private

Definition at line 158 of file RMC01AnalysisManager.hh.

fGamma_current_rmatrix_vs_proton_prim_energy

G4AnaH2* RMC01AnalysisManager::fGamma_current_rmatrix_vs_proton_prim_energy

private

Definition at line 163 of file RMC01AnalysisManager.hh.

fInstance

RMC01AnalysisManager * RMC01AnalysisManager::fInstance = 0

staticprivate

Definition at line 106 of file RMC01AnalysisManager.hh.

fMean_edep

G4double RMC01AnalysisManager::fMean_edep

private

Definition at line 141 of file RMC01AnalysisManager.hh.

fMsg

RMC01AnalysisManagerMessenger* RMC01AnalysisManager::fMsg

private

Definition at line 122 of file RMC01AnalysisManager.hh.

fNb_evt_modulo_for_convergence_test

G4int RMC01AnalysisManager::fNb_evt_modulo_for_convergence_test

private

Definition at line 147 of file RMC01AnalysisManager.hh.

fNb_evt_per_adj_evt

G4int RMC01AnalysisManager::fNb_evt_per_adj_evt

private

Definition at line 180 of file RMC01AnalysisManager.hh.

fNentry

G4double RMC01AnalysisManager::fNentry

private

Definition at line 140 of file RMC01AnalysisManager.hh.

fPrecision_to_reach

G4double RMC01AnalysisManager::fPrecision_to_reach

private

Definition at line 145 of file RMC01AnalysisManager.hh.

fPrimPDG_ID

G4int RMC01AnalysisManager::fPrimPDG_ID

private

Definition at line 174 of file RMC01AnalysisManager.hh.

fPrimSpectrumType

PRIM_SPECTRUM_TYPE RMC01AnalysisManager::fPrimSpectrumType

private

Definition at line 173 of file RMC01AnalysisManager.hh.

fProton_current

G4AnaH1* RMC01AnalysisManager::fProton_current

private

Definition at line 128 of file RMC01AnalysisManager.hh.

fProton_current_rmatrix_vs_proton_prim_energy

G4AnaH2* RMC01AnalysisManager::fProton_current_rmatrix_vs_proton_prim_energy

private

Definition at line 162 of file RMC01AnalysisManager.hh.

fRelative_error

G4double RMC01AnalysisManager::fRelative_error

private

Definition at line 143 of file RMC01AnalysisManager.hh.

fStop_run_if_precision_reached

G4bool RMC01AnalysisManager::fStop_run_if_precision_reached

private

Definition at line 146 of file RMC01AnalysisManager.hh.

fTimer

G4Timer* RMC01AnalysisManager::fTimer

private

Definition at line 184 of file RMC01AnalysisManager.hh.

---

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

• RMC01AnalysisManager.hh
• RMC01AnalysisManager.cc