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