GammaTherapy/src/HistoManager.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 HistoManager.cc
0027 /// \brief Implementation of the HistoManager class
0028
0029 #include "HistoManager.hh"
0030
0031 #include "DetectorConstruction.hh"
0032
0033 #include "G4SystemOfUnits.hh"
0034
0035 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0036
0037 HistoManager::HistoManager()
0038 {
0039   fNBinsZ = 60;
0040   fNBinsR = 80;
0041   fNBinsE = 200;
0042
0043   fAbsorberZ = 300. * mm;
0044   fAbsorberR = 200. * mm;
0045   fScoreZ = 100. * mm;
0046   fMaxEnergy = 50. * MeV;
0047
0048   fStepZ = fStepR = fStepE = 0.0;
0049
0050   Book();
0051 }
0052
0053 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0054
0055 HistoManager::~HistoManager() {}
0056
0057 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0058
0059 void HistoManager::Book()
0060 {
0061   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0062
0063   // Create or get analysis manager
0064   analysisManager->SetDefaultFileType("root");
0065   analysisManager->SetVerboseLevel(1);
0066   analysisManager->SetActivation(true);  // enable inactivation of histograms
0067
0068   // Creating an 1-dimensional histograms in the root directory of the tree
0069   fHisto.assign(10, 0);
0070   int iHisto = 0;
0071   fHisto[iHisto] = analysisManager->CreateH1(
0072     "10", "Energy deposit at radius (mm) normalised on 1st channel", fNBinsR, 0., fAbsorberR / mm);
0073
0074   iHisto++;
0075   fHisto[iHisto] = analysisManager->CreateH1(
0076     "11", "Energy deposit at radius (mm) normalised to integral", fNBinsR, 0., fAbsorberR / mm);
0077
0078   iHisto++;
0079   fHisto[iHisto] = analysisManager->CreateH1(
0080     "12", "Energy deposit (MeV/kg/electron) at radius (mm)", fNBinsR, 0., fAbsorberR / mm);
0081
0082   iHisto++;
0083   fHisto[iHisto] = analysisManager->CreateH1("13", "Energy profile (MeV/kg/electron) over Z (mm)",
0084                                              fNBinsZ, 0., fAbsorberZ / mm);
0085
0086   iHisto++;
0087   fHisto[iHisto] =
0088     analysisManager->CreateH1("14", "Energy profile (MeV/kg/electron) over Z (mm) at Central Voxel",
0089                               fNBinsZ, 0., fAbsorberZ / mm);
0090
0091   iHisto++;
0092   fHisto[iHisto] = analysisManager->CreateH1("15", "Energy (MeV) of fGamma produced in the target",
0093                                              fNBinsE, 0., fMaxEnergy / MeV);
0094
0095   iHisto++;
0096   fHisto[iHisto] = analysisManager->CreateH1("16", "Energy (MeV) of fGamma before phantom", fNBinsE,
0097                                              0., fMaxEnergy / MeV);
0098
0099   iHisto++;
0100   fHisto[iHisto] = analysisManager->CreateH1("17", "Energy (MeV) of electrons produced in phantom",
0101                                              fNBinsE, 0., fMaxEnergy / MeV);
0102
0103   iHisto++;
0104   fHisto[iHisto] = analysisManager->CreateH1("18", "Energy (MeV) of electrons produced in target",
0105                                              fNBinsE, 0., fMaxEnergy / MeV);
0106
0107   iHisto++;
0108   fHisto[iHisto] = analysisManager->CreateH1(
0109     "19", "Gamma Energy Fluence (MeV/cm2) at radius(mm) in front of phantom", fNBinsR, 0.,
0110     fAbsorberR / mm);
0111
0112   // Create all histograms as inactivated
0113   // as we have not yet set nbins, vmin, vmax
0114   for (int i = 0; i < iHisto + 1; i++)
0115     analysisManager->SetH1Activation(i, false);
0116 }
0117
0118 void HistoManager::Update(DetectorConstruction* det, bool bForceActivation)
0119 {
0120   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0121
0122   if (bForceActivation) {
0123     for (int i = 0; i < (int)fHisto.size(); i++)
0124       analysisManager->SetH1Activation(fHisto[i], true);
0125     analysisManager->SetActivation(true);
0126   }
0127
0128   if (analysisManager->IsActive()) {
0129     // Check nBinsR / fAbsorberR histograms
0130     if (det->GetNumberDivR() != fNBinsR || std::fabs(det->GetAbsorberR() - fAbsorberR) > 0.01 * mm)
0131     {
0132       fNBinsR = det->GetNumberDivR();
0133       fAbsorberR = det->GetAbsorberR();
0134       std::vector<G4int> histoId{0, 1, 2, 9};
0135       for (auto v : histoId) {
0136         analysisManager->SetH1(fHisto[v], fNBinsR, 0., fAbsorberR / mm);
0137       }
0138     }
0139
0140     // Check nBinsZ / fAbsorberZ histograms
0141     if (det->GetNumberDivZ() != fNBinsZ || std::fabs(det->GetAbsorberZ() - fAbsorberZ) > 0.01 * mm)
0142     {
0143       fNBinsZ = det->GetNumberDivZ();
0144       fAbsorberZ = det->GetAbsorberZ();
0145       std::vector<G4int> histoId{3, 4};
0146       for (auto v : histoId) {
0147         analysisManager->SetH1(fHisto[v], fNBinsZ, 0., fAbsorberZ / mm);
0148       }
0149     }
0150
0151     // Check nBinsE / fAbsorberE histograms
0152     if (det->GetNumberDivE() != fNBinsE || std::fabs(det->GetMaxEnergy() - fMaxEnergy) > 0.01) {
0153       fNBinsE = det->GetNumberDivE();
0154       fMaxEnergy = det->GetMaxEnergy();
0155       std::vector<G4int> histoId{5, 6, 7, 8};
0156       for (auto v : histoId) {
0157         analysisManager->SetH1(fHisto[v], fNBinsE, 0., fMaxEnergy / MeV);
0158       }
0159     }
0160   }
0161 }
0162
0163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0164
0165 void HistoManager::DumpHistoParameters()
0166 {
0167   if (!G4AnalysisManager::Instance()->IsActive()) return;
0168
0169   for (int i = 0; i < (int)fHisto.size(); i++) {
0170     G4int histoId = fHisto[i];
0171     G4String title = G4AnalysisManager::Instance()->GetH1Title(histoId);
0172     G4int nbins = G4AnalysisManager::Instance()->GetH1Nbins(histoId);
0173     G4double xmin = G4AnalysisManager::Instance()->GetH1Xmin(histoId);
0174     G4double xmax = G4AnalysisManager::Instance()->GetH1Xmax(histoId);
0175     G4double width = G4AnalysisManager::Instance()->GetH1Width(histoId);
0176     G4cout << "Histogram parameters : " << i << " " << histoId << " : " << nbins << " ";
0177     G4cout << xmin << "/" << xmax << " " << width << G4endl;
0178   }
0179 }