radiobiology/src/RBEAccumulable.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 RBEAccumulable.cc
0027 /// \brief Implementation of the RadioBio::RBEAccumulable class
0028
0029 #include "RBEAccumulable.hh"
0030
0031 #include "G4ParticleDefinition.hh"
0032
0033 #include "Hit.hh"
0034 #include "Manager.hh"
0035 #include "RBE.hh"
0036 #include "VoxelizedSensitiveDetector.hh"
0037 #include <G4SystemOfUnits.hh>
0038
0039 #include <tuple>
0040
0041 namespace RadioBio
0042 {
0043
0044 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0045
0046 RBEAccumulable::RBEAccumulable() : VRadiobiologicalAccumulable("RBE") {}
0047
0048 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0049
0050 void RBEAccumulable::Merge(const G4VAccumulable& rhs)
0051 {
0052   if (GetVerboseLevel() > 1) {
0053     G4cout << "RBEAccumulable::Merge()" << G4endl;
0054   }
0055   const RBEAccumulable& other = dynamic_cast<const RBEAccumulable&>(rhs);
0056   fAlphaNumerator += other.fAlphaNumerator;
0057   fDenominator += other.fDenominator;
0058   fBetaNumerator += other.fBetaNumerator;
0059 }
0060
0061 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0062
0063 void RBEAccumulable::Reset()
0064 {
0065   if (GetVerboseLevel() > 0) {
0066     G4cout << "RBEAccumulable::Reset()" << G4endl;
0067   }
0068   if (fInitialized) {
0069     fAlphaNumerator = 0.0;
0070     fBetaNumerator = 0.0;
0071     fDenominator = 0.0;
0072   }
0073   else {
0074     Initialize();
0075   }
0076 }
0077
0078 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0079
0080 // to accumulate given the hit
0081 void RBEAccumulable::Accumulate(Hit* hit)
0082 {
0083   G4double kineticEnergy = hit->GetEkinMean();
0084   G4int A = hit->GetPartType()->GetAtomicMass();
0085   G4double energyDeposit = hit->GetDeltaE();
0086   G4double DX = hit->GetTrackLength();
0087   G4int Z = hit->GetPartType()->GetAtomicNumber();
0088   G4int i = hit->GetXindex();
0089   G4int j = hit->GetYindex();
0090   G4int k = hit->GetZindex();
0091
0092   // If A is zero, return immediately to avoid division by zero
0093   if (!A) return;
0094
0095   Accumulate(kineticEnergy / A, energyDeposit, DX, Z, i, j, k);
0096 }
0097
0098 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0099
0100 void RBEAccumulable::Accumulate(G4double E, G4double energyDeposit, G4double dX, G4int Z, G4int i,
0101                                 G4int j, G4int k)
0102 {
0103   if (!fInitialized) {
0104     G4Exception("RBEAccumulable::Accumulate", "NotInitialized", FatalException,
0105                 "Accumulable not initialized. Must be a programming error.");
0106   }
0107   if (GetVerboseLevel() > 2) {
0108     G4cout << "RBEAccumulable::Accumulate() in " << i << ", " << j << ", " << k << G4endl;
0109   }
0110   if (energyDeposit <= 0) {
0111     return;
0112   }
0113
0114   // Get the global voxel number for the given indexes
0115   size_t n = VoxelizedSensitiveDetector::GetInstance()->GetThisVoxelNumber(i, j, k);
0116
0117   // Calculate only for hadrons that traveled finite range and released finite energy.
0118   if ((Z >= 1) && (dX > 0) && (E > 0)) {
0119     RBE* rbe = dynamic_cast<RBE*>(Manager::GetInstance()->GetQuantity("RBE"));
0120     std::tuple<G4double, G4double> alpha_beta = rbe->GetHitAlphaAndBeta(E, Z);
0121     fDenominator[n] += energyDeposit;
0122     fAlphaNumerator[n] += std::get<0>(alpha_beta) * energyDeposit;
0123     fBetaNumerator[n] += std::sqrt(std::get<1>(alpha_beta)) * energyDeposit;
0124   }
0125 }
0126
0127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0128
0129 G4int RBEAccumulable::GetVerboseLevel() const
0130 {
0131   // Return same verbosity of RBE class
0132   return Manager::GetInstance()->GetQuantity("RBE")->GetVerboseLevel();
0133 }
0134
0135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0136
0137 void RBEAccumulable::Initialize()
0138 {
0139   if (GetVerboseLevel() > 0) {
0140     G4cout << "RBEAccumulable::Initialize(): " << G4endl;
0141   }
0142
0143   auto voxSensDet = VoxelizedSensitiveDetector::GetInstance();
0144
0145   fVoxelsAlongX = voxSensDet->GetVoxelNumberAlongX();
0146   fVoxelsAlongY = voxSensDet->GetVoxelNumberAlongY();
0147   fVoxelsAlongZ = voxSensDet->GetVoxelNumberAlongZ();
0148   fVoxels = fVoxelsAlongX * fVoxelsAlongY * fVoxelsAlongZ;
0149
0150   if (GetVerboseLevel() > 1) {
0151     G4cout << fVoxels << " voxels." << G4endl;
0152   }
0153
0154   fAlphaNumerator = array_type(0.0, fVoxels);
0155   fBetaNumerator = array_type(0.0, fVoxels);
0156   fDenominator = array_type(0.0, fVoxels);
0157   fInitialized = true;
0158 }
0159
0160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0161
0162 }  // namespace RadioBio