microtrack/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 RunAction.cc
0027 /// \brief Implementation of the RunAction class
0028
0029 #include "RunAction.hh"
0030
0031 #include "G4AnalysisManager.hh"
0032 #include "G4RunManager.hh"
0033 #include "G4SystemOfUnits.hh"
0034 #include "G4UnitsTable.hh"
0035 #include "globals.hh"
0036
0037 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0038
0039 RunAction::RunAction()
0040 {
0041   // Create analysis manager
0042   auto analysisManager = G4AnalysisManager::Instance();
0043
0044   analysisManager->SetVerboseLevel(1);
0045   analysisManager->SetNtupleMerging(true);
0046
0047   // Define binning for each sort of histogram
0048   //
0049   //  Kinetic energy at the entrance and exit histograms
0050   const G4double kinEmin = 0.;
0051   const G4double kinEmax = 300.;
0052   const G4double kinEbinWidth = 0.1;
0053   const G4int kinEbins =
0054     static_cast<G4int>((kinEmax - kinEmin) / kinEbinWidth);
0055
0056   // Logarithmic binning for energy histograms
0057   const G4int minLog10E = -4.;
0058   const G4int maxLog10E = +4.;
0059   const G4int nBinsLog10E = (maxLog10E - minLog10E) * 50;
0060
0061   G4double binsLog10E[nBinsLog10E + 1];
0062   G4double binWidthLog10E =
0063     static_cast<G4double>((maxLog10E - minLog10E)) / nBinsLog10E;
0064
0065   for (G4int ii = 0; ii <= nBinsLog10E; ii++) {
0066     binsLog10E[ii] = std::pow(10., binWidthLog10E * ii + minLog10E);
0067   }
0068
0069   std::vector<G4double> vecBinsLog10E(binsLog10E, binsLog10E + nBinsLog10E + 1);
0070
0071   // Logarithmic binning for weighted numbers histograms
0072
0073   const G4int minLog10W = 0;
0074   const G4int maxLog10W = +6;
0075   const G4int nBinsLog10W = (maxLog10W - minLog10W) * 50;
0076
0077   G4double binsLog10W[nBinsLog10W + 1];
0078   G4double binWidthLog10W =
0079     static_cast<G4double>((maxLog10W - minLog10W)) / nBinsLog10W;
0080
0081   for (G4int ii = 0; ii <= nBinsLog10W; ii++) {
0082     binsLog10W[ii] = std::pow(10., binWidthLog10W * ii + minLog10W);
0083   }
0084
0085   std::vector<G4double> vecBinsLog10W(binsLog10W, binsLog10W + nBinsLog10W + 1);
0086
0087   // Linear binning for counter histograms
0088   G4int minCount = 0;
0089   G4int maxCount = 20000;
0090   G4int nBinsCount = (maxCount - minCount) / 100;
0091
0092   // Create histograms
0093   //
0094
0095   // Logarithmic binning  Energy histograms
0096   analysisManager->CreateH1(
0097     "fe",
0098     "Energy imparted per event [keV] (log binning)",
0099     vecBinsLog10E);
0100
0101   analysisManager->CreateH1(
0102     "efe",
0103     "Weighted energy imparted per event [keV] (log binning)",
0104     vecBinsLog10E);
0105
0106   analysisManager->CreateH1(
0107     "e2fe",
0108     "Squared-weighted energy imparted per event [keV] (log binning)",
0109     vecBinsLog10E);
0110
0111   // Logarithmic binning for lineal energy histograms
0112   analysisManager->CreateH1("fy", "Lineal energy [keV/um] (log binning)",
0113                             vecBinsLog10E);
0114
0115   analysisManager->CreateH1(
0116     "yfy",
0117     "Dose-weighted lineal energy [keV/um] (log binning)", vecBinsLog10E);
0118
0119   analysisManager->CreateH1(
0120     "y2fy",
0121     "Squared-weighted lineal energy [keV/um] (log binning)", vecBinsLog10E);
0122
0123   // Logarithmic binning for specific energy histograms
0124   analysisManager->CreateH1(
0125     "fz",
0126     "Single-event specific energy [Gy] (log binning)",
0127     vecBinsLog10E);
0128
0129   analysisManager->CreateH1(
0130     "zfz",
0131     "Dose-weighted single-event specific energy [Gy] (log binning)",
0132     vecBinsLog10E);
0133
0134   analysisManager->CreateH1(
0135     "z2fz",
0136     "Squared-weighted single-event specific energy [Gy] (log binning)",
0137     vecBinsLog10E);
0138
0139   // Counters histograms
0140   analysisManager->CreateH1(
0141     "Nsel", "Number of selectable hits per event",
0142     nBinsCount, minCount, maxCount);
0143
0144   analysisManager->CreateH1(
0145     "Nsite", "Number of hits in site", nBinsCount,
0146     minCount, maxCount);
0147
0148   analysisManager->CreateH1("Nint",
0149     "Number of selectable hits in site", nBinsCount,
0150     minCount, maxCount);
0151
0152   // Kinetic energy at the entrance and exit histograms
0153   analysisManager->CreateH1("KinE_in", "Kinetic energy at the entrance [MeV]",
0154                             kinEbins, kinEmin, kinEmax);
0155   analysisManager->CreateH1("KinE_out", "Kinetic energy at the exit [MeV]",
0156                             kinEbins, kinEmin, kinEmax);
0157
0158   // 2D histogram for Nsite vs weighted Edep
0159   analysisManager->CreateH2(
0160     "Nsite_vs_e",
0161     "Number of hits in site vs energy imparted [keV] (log-log)",
0162     vecBinsLog10E, vecBinsLog10W);
0163 }
0164
0165
0166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0167
0168 RunAction::~RunAction() = default;
0169
0170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0171
0172 void RunAction::BeginOfRunAction(const G4Run* /*aRun*/)
0173 {
0174   // inform the runManager to save random number seed
0175   // G4RunManager::GetRunManager()->SetRandomNumberStore(true);
0176
0177   // Get analysis manager
0178   auto analysisManager = G4AnalysisManager::Instance();
0179
0180   // Open an output file
0181   // The file extension will set the choice of the output format
0182   G4String fileName = "microtrack.root";
0183   // Other formats supported: .csv, .hdf5, .xml
0184   // G4String fileName = "microtrack.csv";
0185   // G4String fileName = "microtrack.hdf5";
0186   // G4String fileName = "microtrack.xml";
0187   analysisManager->OpenFile(fileName);
0188   G4cout << "Using " << analysisManager->GetType() << G4endl;
0189 }
0190
0191
0192 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0193
0194 void RunAction::EndOfRunAction(const G4Run* /*aRun*/)
0195 {
0196   auto analysisManager = G4AnalysisManager::Instance();
0197
0198   if (IsMaster() && analysisManager->GetH1(1)) {
0199
0200     G4cout << G4endl;
0201
0202     G4cout << "----> print histogram statistics for the entire run: "
0203            << G4endl << G4endl;
0204
0205     // print histogram statistics
0206     //
0207
0208     G4cout << "  Single-event energy imparted:\n"
0209            << "  ----------------------------"
0210            << G4endl << G4endl;
0211
0212     G4cout << "    Frequency-mean: \\varepsilon_{1,F} = "
0213            << analysisManager->GetH1(0)->mean() << " keV "
0214            << " (rms = " << analysisManager->GetH1(0)->rms() << " keV)"
0215            << G4endl;
0216
0217     G4cout << "    Dose-mean: \\varepsilon_{1,D} = "
0218            << analysisManager->GetH1(1)->mean() << " keV "
0219            << " (rms = " << analysisManager->GetH1(1)->rms() << " keV)"
0220            << G4endl;
0221
0222     G4cout << "    Squared-weighted histogram: mean = "
0223            << analysisManager->GetH1(2)->mean() << " keV "
0224            << " (rms = " << analysisManager->GetH1(2)->rms() << " keV)"
0225            << G4endl;
0226
0227     G4cout << G4endl
0228            << "  Lineal energy:\n"
0229            << "  -------------"
0230            << G4endl << G4endl;
0231
0232     G4cout << "    Frequency-mean: y_F = "
0233            << analysisManager->GetH1(3)->mean() << " keV/um "
0234            << " (rms = " << analysisManager->GetH1(3)->rms() << " keV/um)"
0235            << G4endl;
0236
0237     G4cout << "    Dose-mean: y_D = "
0238            << analysisManager->GetH1(4)->mean() << " keV/um "
0239            << " (rms = " << analysisManager->GetH1(4)->rms() << " keV/um)"
0240            << G4endl;
0241
0242     G4cout << "    Squared-weighted histogram: mean = "
0243            << analysisManager->GetH1(5)->mean() << " keV/um "
0244            << " (rms = " << analysisManager->GetH1(5)->rms() << " keV/um)"
0245            << G4endl;
0246
0247     G4cout << G4endl
0248            << "  Single-event specific energy:\n"
0249            << "  ----------------------------"
0250            << G4endl << G4endl;
0251
0252     G4cout << "    Frequency-mean: z_{1,F} = "
0253            << analysisManager->GetH1(6)->mean() << " Gy "
0254            << " (rms = " << analysisManager->GetH1(6)->rms() << " Gy)"
0255            << G4endl;
0256
0257     G4cout << "    Dose-mean: z_{1,D} = "
0258            << analysisManager->GetH1(7)->mean() << " Gy "
0259            << " (rms = " << analysisManager->GetH1(7)->rms() << " Gy)"
0260            << G4endl;
0261
0262     G4cout << "    Squared-weighted histogram: mean = "
0263            << analysisManager->GetH1(8)->mean() << " Gy"
0264            << " (rms = " << analysisManager->GetH1(8)->rms() << " Gy)"
0265            << G4endl;
0266
0267     G4cout << G4endl
0268            << "  Number of hits per event:\n"
0269            << "  ------------------------"
0270            << G4endl << G4endl;
0271
0272     G4cout << "    Eligible for site random placement, N_{sel}: mean = "
0273            << analysisManager->GetH1(9)->mean()
0274            << " (rms = " << analysisManager->GetH1(9)->rms() << ")"
0275            << G4endl;
0276
0277     G4cout << "    Within the site, N_{site}: mean = "
0278            << analysisManager->GetH1(10)->mean()
0279            << " (rms = " << analysisManager->GetH1(10)->rms() << ")"
0280            << G4endl;
0281
0282     G4cout << "    Within the site and eligible for site random placement,"
0283            << " N_{int}: mean = " << analysisManager->GetH1(11)->mean()
0284            << " (rms = " << analysisManager->GetH1(11)->rms() << ")"
0285            << G4endl;
0286
0287     G4cout << G4endl
0288            << "  Kinetic energy of the primary particle:\n"
0289            << "  --------------------------------------"
0290            << G4endl << G4endl;
0291
0292     G4cout << "    At entrance of SDbox, T_{in}: mean = "
0293            << G4BestUnit(analysisManager->GetH1(12)->mean(), "Energy")
0294            << " (rms = "
0295            << G4BestUnit(analysisManager->GetH1(12)->rms(), "Energy") << ")"
0296            << G4endl;
0297
0298     G4cout << "    At exit of SDbox, T_{out}: mean = "
0299            << G4BestUnit(analysisManager->GetH1(13)->mean(), "Energy")
0300            << " (rms = "
0301            << G4BestUnit(analysisManager->GetH1(13)->rms(), "Energy") << ")"
0302            << G4endl;
0303
0304     G4cout << G4endl;
0305   }
0306
0307   // save histograms & ntuple
0308   //
0309   analysisManager->Write();
0310   analysisManager->CloseFile();
0311
0312 }
0313
0314 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......