src/RunAction.cc

Go to the documentation of this file.

//
// $Id: RunAction.cc 168 2008-09-16 13:23:52Z maire $
// 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

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

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

#ifdef G4ANALYSIS_USE
#include "AIDA/AIDA.h"
#endif

#include <iomanip>

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

RunAction::RunAction(DetectorConstruction* det)
:m_detector(det)
{
  BookHistoFile( "histograms", "root");
}

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

RunAction::~RunAction()
{ 
  SaveHistoFile();
}

//....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.);
}

//....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; 
}

//....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;
  
  // write m_depthDose array on a file
  //
  ////WriteDepthDose("m_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);
  
  std::ios::fmtflags mode = File.flags();  
  File.setf( std::ios::scientific, std::ios::floatfield );
  G4int prec = File.precision(4);
      
  File << "  Longitudinal depth dose distribution \n " 
       << "\n  zLocal (mm)\t  m_depthDose(gray) \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;
  }
  
  // restaure default formats
  File.setf(mode,std::ios::floatfield);
  File.precision(prec);         
}

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

void RunAction::BookHistoFile(G4String fileName, G4String type)
{          
  histoFile = fileName + "." + type;
  m_af = 0; m_tree = 0;
  m_histo[0] = m_histo[1] = 0;
  
#ifdef G4ANALYSIS_USE
   
  // check allowed type
  if ((type != "root") && (type != "hbook")) {
    G4cout << "\n----> Error in BookHisto. Type = " << type
           << " not recognized" <<  G4endl;
    return;
  }
  
 // Creating the analysis factory
 m_af = AIDA_createAnalysisFactory();
 
 if (m_af) {
   // Creating the m_tree factory
   AIDA::ITreeFactory* tf = m_af->createTreeFactory();
 
   // Creating a m_tree mapped to an hbook file.
   G4bool readOnly  = false;
   G4bool createNew = true;
   G4String options =  "export=root";
   m_tree = tf->create(histoFile,type,readOnly,createNew, options);
   delete tf;
   
   if (m_tree) {
     // Creating a histogram factory
     AIDA::IHistogramFactory* hf = m_af->createHistogramFactory(*m_tree);

     // Creating the histograms
     m_histo[0]=hf->createHistogram1D
                   ("1","total energy deposit in astronaut",100,0.,500*MeV);
                   
     delete hf;
     G4cout << "\n----> Histogram m_tree is opened in " << histoFile << G4endl;
   }
 }
#endif  
}

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

void RunAction::SaveHistoFile()
{
#ifdef G4ANALYSIS_USE
  m_tree->commit();       // Writing the histograms to the file
  m_tree->close();        // and closing the m_tree (and the file)
  G4cout << "\n----> Histogram m_tree is saved in " << histoFile << G4endl;
      
  delete m_tree;
  delete m_af;
#endif
}

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

---