TestEm7/src/RunAction.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file electromagnetic/TestEm7/src/RunAction.cc
0027 /// \brief Implementation of the RunAction class
0028 //
0029 //
0030 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0031 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0032
0033 #include "RunAction.hh"
0034
0035 #include "DetectorConstruction.hh"
0036 #include "PhysicsList.hh"
0037 #include "PrimaryGeneratorAction.hh"
0038 #include "StepMax.hh"
0039
0040 #include "G4Run.hh"
0041 #include "G4RunManager.hh"
0042 #include "G4SystemOfUnits.hh"
0043 #include "G4UnitsTable.hh"
0044 #include "G4ios.hh"
0045 #include "Randomize.hh"
0046
0047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0048
0049 RunAction::RunAction(DetectorConstruction* det, PhysicsList* phys, PrimaryGeneratorAction* kin)
0050   : G4UserRunAction(),
0051     fAnalysisManager(0),
0052     fDetector(det),
0053     fPhysics(phys),
0054     fKinematic(kin),
0055     fTallyEdep(new G4double[kMaxTally]),
0056     fProjRange(0.),
0057     fProjRange2(0.),
0058     fEdeptot(0.),
0059     fEniel(0.),
0060     fNbPrimarySteps(0),
0061     fRange(0)
0062 {
0063   // Book predefined histograms
0064   BookHisto();
0065 }
0066
0067 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0068
0069 RunAction::~RunAction()
0070 {
0071   delete[] fTallyEdep;
0072 }
0073
0074 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0075
0076 void RunAction::BeginOfRunAction(const G4Run* aRun)
0077 {
0078   G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
0079
0080   if (!fAnalysisManager) {
0081     BookHisto();
0082   }
0083
0084   CLHEP::HepRandom::showEngineStatus();
0085
0086   // initialize projected range, tallies, Ebeam, and book histograms
0087   //
0088   fNbPrimarySteps = 0;
0089   fRange = 0;
0090   fProjRange = fProjRange2 = 0.;
0091   fEdeptot = fEniel = 0.;
0092   for (G4int j = 0; j < kMaxTally; ++j) {
0093     fTallyEdep[j] = 0.;
0094   }
0095   fKinematic->ResetEbeamCumul();
0096
0097   if (fAnalysisManager->IsActive()) {
0098     fAnalysisManager->OpenFile();
0099
0100     // histogram "1" is defined by the length of the target
0101     // zoomed histograms are defined by UI command
0102     G4double length = fDetector->GetAbsorSizeX();
0103     G4double stepMax = fPhysics->GetStepMaxProcess()->GetMaxStep();
0104     G4int nbBins = 100;
0105     if (stepMax < DBL_MAX) {
0106       G4int nb = (G4int)(0.5 + length / stepMax);
0107       nbBins = std::min(std::max(nbBins, nb), 10000);
0108     }
0109     fAnalysisManager->SetH1(1, nbBins, 0., length, "mm");
0110   }
0111 }
0112
0113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0114
0115 void RunAction::EndOfRunAction(const G4Run* aRun)
0116 {
0117   G4int nbofEvents = aRun->GetNumberOfEvent();
0118   if (nbofEvents == 0) return;
0119
0120   // run conditions
0121   //
0122   const G4Material* material = fDetector->GetAbsorMaterial();
0123   G4double density = material->GetDensity();
0124
0125   G4String particle = fKinematic->GetParticleGun()->GetParticleDefinition()->GetParticleName();
0126   G4double energy = fKinematic->GetParticleGun()->GetParticleEnergy();
0127   G4cout << "\n The run consists of " << nbofEvents << " " << particle << " of "
0128          << G4BestUnit(energy, "Energy") << " through "
0129          << G4BestUnit(fDetector->GetAbsorSizeX(), "Length") << " of " << material->GetName()
0130          << " (density: " << G4BestUnit(density, "Volumic Mass") << ")" << G4endl;
0131
0132   // compute projected range and straggling
0133   //
0134   if (fRange > 0) {
0135     fProjRange /= fRange;
0136     fProjRange2 /= fRange;
0137   }
0138   G4double rms = fProjRange2 - fProjRange * fProjRange;
0139   if (rms > 0.)
0140     rms = std::sqrt(rms);
0141   else
0142     rms = 0.;
0143
0144   G4double nstep = G4double(fNbPrimarySteps) / G4double(nbofEvents);
0145
0146   G4cout.precision(6);
0147   G4cout << "\n Projected Range= " << G4BestUnit(fProjRange, "Length")
0148          << "   rms= " << G4BestUnit(rms, "Length") << G4endl;
0149   G4cout << " Mean number of primary steps = " << nstep << G4endl;
0150
0151   // compute energy deposition and niel
0152   //
0153   fEdeptot /= nbofEvents;
0154   G4cout << " Total energy deposit= " << G4BestUnit(fEdeptot, "Energy") << G4endl;
0155   fEniel /= nbofEvents;
0156   G4cout << " niel energy deposit = " << G4BestUnit(fEniel, "Energy") << G4endl;
0157
0158   // print dose in tallies
0159   //
0160   G4int tallyNumber = fDetector->GetTallyNumber();
0161   if (tallyNumber > 0) {
0162     G4double Ebeam = fKinematic->GetEbeamCumul();
0163     G4cout << "\n---------------------------------------------------------\n";
0164     G4cout << " Cumulated Doses : \tEdep      \tEdep/Ebeam \tDose" << G4endl;
0165     for (G4int j = 0; j < tallyNumber; ++j) {
0166       G4double Edep = fTallyEdep[j], ratio = 100 * Edep / Ebeam;
0167       G4double tallyMass = fDetector->GetTallyMass(j);
0168       G4double Dose = Edep / tallyMass;
0169       G4cout << " tally " << j << ": \t \t" << G4BestUnit(Edep, "Energy") << "\t" << ratio
0170              << " % \t" << G4BestUnit(Dose, "Dose") << G4endl;
0171     }
0172     G4cout << "\n---------------------------------------------------------\n";
0173     G4cout << G4endl;
0174   }
0175
0176   if (fAnalysisManager->IsActive()) {
0177     // normalize histograms
0178     //
0179     for (G4int j = 1; j < 3; ++j) {
0180       G4double binWidth = fAnalysisManager->GetH1Width(j);
0181       G4double fac = (mm / MeV) / (nbofEvents * binWidth);
0182       fAnalysisManager->ScaleH1(j, fac);
0183     }
0184     fAnalysisManager->ScaleH1(3, 1. / nbofEvents);
0185
0186     // save histograms
0187     fAnalysisManager->Write();
0188     fAnalysisManager->CloseFile();
0189   }
0190
0191   // show Rndm status
0192   //
0193   CLHEP::HepRandom::showEngineStatus();
0194 }
0195
0196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0197
0198 void RunAction::BookHisto()
0199 {
0200   // Create or get analysis manager
0201   // The choice of analysis technology is done via selection of a namespace
0202   // in HistoManager.hh
0203   fAnalysisManager = G4AnalysisManager::Instance();
0204   fAnalysisManager->SetDefaultFileType("root");
0205   fAnalysisManager->SetFileName("testem7");
0206   fAnalysisManager->SetVerboseLevel(1);
0207   fAnalysisManager->SetActivation(true);  // enable inactivation of histograms
0208
0209   // Define histograms start values
0210   const G4int kMaxHisto = 4;
0211   const G4String id[] = {"h0", "h1", "h2", "h3"};
0212   const G4String title[] = {
0213     "dummy",  // 0
0214     "Edep (MeV/mm) along absorber ",  // 1
0215     "Edep (MeV/mm) along absorber zoomed",  // 2
0216     "projectile range"  // 3
0217   };
0218
0219   // Default values (to be reset via /analysis/h1/set command)
0220   G4int nbins = 100;
0221   G4double vmin = 0.;
0222   G4double vmax = 100.;
0223
0224   // Create all histograms as inactivated
0225   // as we have not yet set nbins, vmin, vmax
0226   for (G4int k = 0; k < kMaxHisto; ++k) {
0227     G4int ih = fAnalysisManager->CreateH1(id[k], title[k], nbins, vmin, vmax);
0228     G4bool activ = false;
0229     if (k == 1) activ = true;
0230     fAnalysisManager->SetH1Activation(ih, activ);
0231   }
0232 }
0233
0234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......