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