src/RunAction.cc

Go to the documentation of this file.

//
// $Id: RunAction.cc 222 2008-11-21 08:12:11Z jacquem $
// 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

#include "RunAction.hh"
#include "DetectorConstruction.hh"

#include "G4Run.hh"
#include "G4RunManager.hh"
#include "G4UnitsTable.hh"

#include <iomanip>

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

RunAction::RunAction(DetectorConstruction* det)
:m_detector(det)
{ }

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

RunAction::~RunAction()
{ }

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void RunAction::BeginOfRunAction(const G4Run* run)
{  
  G4cout << "### Run " << run->GetRunID() << " start." << G4endl;
  
  //initialize total energy deposit
  //
  m_sumEdep = m_sumEdep2 = 0.;
  m_nbEntry = 0;
  
  //construct vector of m_depthDose distribution
  //
  m_nbSlices = 100;
  m_astronautHeight = m_detector->GetAstronautHeight();
  m_sliceWidth = m_astronautHeight/m_nbSlices;
  m_depthDose.resize(m_nbSlices, 0.);
  m_depthDoseSD   .resize(m_nbSlices, 0.);
  m_depthTotalEkin.resize(m_nbSlices, 0.);
  m_depthNbProton .resize(m_nbSlices, 0 );
  
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void RunAction::SumDepthDose(G4double zLocal, G4double edep)
{
  // check coherence
  if (zLocal > m_astronautHeight) {
    G4cout << "\n --> warning from RunAction::SumDepthDose() : "
           << " zLocal beyond astronaut height : " << zLocal/mm << G4endl;
    return;
  }
  
  // compute slice number               
  G4int k = int(zLocal/m_sliceWidth);
  
  // sum edep
  m_depthDose[k] += edep; 
}


void RunAction::SumDepthDoseSD(G4double zLocal, G4double edep)
{
  // check coherence
  if (zLocal > m_astronautHeight) {
    G4cout << "\n --> warning from RunAction::SumDepthDose() : "
           << " zLocal beyond astronaut height : " << zLocal/mm << G4endl;
    return;
  }
  
  // compute slice number               
  G4int k = int(zLocal/m_sliceWidth);
  
  // sum edep
  m_depthDoseSD[k] += edep; 
}


void RunAction::SumDepthDoseSD(G4int k, G4double edep, G4int nbProton, G4double totalEkin)
{
  // check coherence
  if (k >= m_nbSlices) {
    G4cout << "\n --> warning from RunAction::SumDepthDoseSD() : "
           << " k beyond nb slices : " << k << G4endl;
    return;
  }
  
  // sum edep, number of proton, and total proton eKin:
  m_depthDoseSD[k]    += edep; 
  m_depthNbProton[k]  += nbProton;
  m_depthTotalEkin[k] += totalEkin;
}


G4int RunAction::GetSliceNumber(G4double zLocal) const
{
  // compute slice number               
  G4int k = int(zLocal/m_sliceWidth);
  
  // check coherence
  if (k >= m_nbSlices) {
    G4cout << "\n --> warning from RunAction::GetSliceNumber() : "
           << " k beyond nb slices : " << k << " returning -1" << G4endl;
    k = -1;
  }
  
  return k;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void RunAction::EndOfRunAction(const G4Run*)
{ 
  G4double meanEd = 0., meanEd2 = 0., rmsEd = 0., error_rel_sum = 0.;
  if (m_nbEntry) { 
    meanEd = m_sumEdep/m_nbEntry; meanEd2 = m_sumEdep2/m_nbEntry;
    G4double variance = meanEd2 - meanEd*meanEd;
    if (variance > 0.) rmsEd = std::sqrt(variance);
    error_rel_sum = rmsEd/(meanEd*std::sqrt(double(m_nbEntry)));
  }
   
  // choose printing format
  std::ios::fmtflags mode = G4cout.flags();
  G4cout.setf(std::ios::fixed,std::ios::floatfield);
  G4int prec = G4cout.precision(3);   
  
  //print total energy deposit
  //
  G4cout
    << "\n ---> Nb of non-empty events = " << m_nbEntry
    << "\n ---> Mean Edep per event    = " << G4BestUnit(meanEd,"Energy")     
    << "\n ---> Total Energy deposited = " << G4BestUnit(m_sumEdep,"Energy")
    << " +- " << G4BestUnit(m_sumEdep*error_rel_sum,"Energy")
    << "  (-> " << 100*error_rel_sum << " %)" 
    << G4endl;
    
  //get mass of astronaut
  //
  G4double mass = m_detector->GetAstronaut()->GetLogicalVolume()->GetMass();
  
  //compute dose
  //
  G4double dose = m_sumEdep/mass;
  G4cout.setf(mode,std::ios::floatfield);
  G4cout.precision(5);     
  G4cout << " ---> Dose = " << G4BestUnit(dose,"Dose") 
         << " +- " << G4BestUnit(dose*error_rel_sum,"Dose") << "\n" 
         << G4endl;
  
  // restaure default formats
  G4cout.setf(mode,std::ios::floatfield);
  G4cout.precision(prec);
        
  // convert m_depthDose array from energy to dose
  //
  for (G4int k=0; k<m_nbSlices; k++)  m_depthDose[k] /= mass;
  //
  for (G4int k=0; k<m_nbSlices; k++)  m_depthDoseSD[k] /= mass;


  
  // write m_depthDose array on a file
  //
  WriteDepthDose("depthDose");
  
  // restaure default formats
  G4cout.setf(mode,std::ios::floatfield);
  G4cout.precision(prec);    
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

#include <fstream>

void RunAction::WriteDepthDose(G4String name)
{ 
  G4String fileName = name + ".ascii";
  std::ofstream File(fileName, std::ios::out);
  
  G4String fileNameSD = name + "SD"+ ".ascii";
  std::ofstream FileSD(fileNameSD, std::ios::out);
  
  std::ios::fmtflags mode = File.flags();  
  File.setf( std::ios::scientific, std::ios::floatfield );
  G4int prec = File.precision(4);

  std::ios::fmtflags modeSD = FileSD.flags();  
  FileSD.setf( std::ios::scientific, std::ios::floatfield );
  G4int precSD = FileSD.precision(4);
      
  File   << "  Longitudinal depth dose distribution \n " 
         << "\n  zLocal (mm)\t  m_depthDose(gray) \n" << G4endl;

  FileSD << "  Longitudinal depth dose distribution \n " 
         << "\n  zLocal (mm)\t  Dose(gray) \t #proton \t <eKin> (MeV) \n" << G4endl;
           
  G4double zLocal;    
  for (G4int k=0; k<m_nbSlices; k++) {
     zLocal = (k+0.5)*m_sliceWidth;
     File << "  "   << zLocal/mm 
          << "\t  " << m_depthDose[k]/gray << G4endl;
  }

  // -- write data from SD :
  for (G4int k=0; k<m_nbSlices; k++) {
      zLocal = (k+0.5)*m_sliceWidth;
      G4double meanEkin(-1.*MeV);
      if (m_depthNbProton[k]>0) meanEkin = m_depthTotalEkin[k]/m_depthNbProton[k];
      FileSD << "  "   << zLocal/mm 
             << "\t  " << m_depthDoseSD[k]/gray
             << "\t  " << m_depthNbProton[k]
             << "\t  " << meanEkin / MeV
             << G4endl;
  }
  
  // restaure default formats
  File.setf(mode,std::ios::floatfield);
  File.precision(prec);

  FileSD.setf(modeSD,std::ios::floatfield);
  FileSD.precision(precSD);
     
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

---