| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 /// \file radiobiology/src/RBEAccumulable.cc
28 /// \brief Implementation of the RadioBio::RBEAccumulable class
29
30 #include "RBEAccumulable.hh"
31
32 #include "G4ParticleDefinition.hh"
33
34 #include "Hit.hh"
35 #include "Manager.hh"
36 #include "RBE.hh"
37 #include "VoxelizedSensitiveDetector.hh"
38 #include <G4SystemOfUnits.hh>
39
40 #include <tuple>
41
42 namespace RadioBio
43 {
44
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
46
47 RBEAccumulable::RBEAccumulable() : VRadiobiologicalAccumulable("RBE") {}
48
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50
51 void RBEAccumulable::Merge(const G4VAccumulable& rhs)
52 {
53 if (GetVerboseLevel() > 1) {
54 G4cout << "RBEAccumulable::Merge()" << G4endl;
55 }
56 const RBEAccumulable& other = dynamic_cast<const RBEAccumulable&>(rhs);
57 fAlphaNumerator += other.fAlphaNumerator;
58 fDenominator += other.fDenominator;
59 fBetaNumerator += other.fBetaNumerator;
60 }
61
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
63
64 void RBEAccumulable::Reset()
65 {
66 if (GetVerboseLevel() > 0) {
67 G4cout << "RBEAccumulable::Reset()" << G4endl;
68 }
69 if (fInitialized) {
70 fAlphaNumerator = 0.0;
71 fBetaNumerator = 0.0;
72 fDenominator = 0.0;
73 }
74 else {
75 Initialize();
76 }
77 }
78
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
80
81 // to accumulate given the hit
82 void RBEAccumulable::Accumulate(Hit* hit)
83 {
84 G4double kineticEnergy = hit->GetEkinMean();
85 G4int A = hit->GetPartType()->GetAtomicMass();
86 G4double energyDeposit = hit->GetDeltaE();
87 G4double DX = hit->GetTrackLength();
88 G4int Z = hit->GetPartType()->GetAtomicNumber();
89 G4int i = hit->GetXindex();
90 G4int j = hit->GetYindex();
91 G4int k = hit->GetZindex();
92
93 // If A is zero, return immediately to avoid division by zero
94 if (!A) return;
95
96 Accumulate(kineticEnergy / A, energyDeposit, DX, Z, i, j, k);
97 }
98
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100
101 void RBEAccumulable::Accumulate(G4double E, G4double energyDeposit, G4double dX, G4int Z, G4int i,
102 G4int j, G4int k)
103 {
104 if (!fInitialized) {
105 G4Exception("RBEAccumulable::Accumulate", "NotInitialized", FatalException,
106 "Accumulable not initialized. Must be a programming error.");
107 }
108 if (GetVerboseLevel() > 2) {
109 G4cout << "RBEAccumulable::Accumulate() in " << i << ", " << j << ", " << k << G4endl;
110 }
111 if (energyDeposit <= 0) {
112 return;
113 }
114
115 // Get the global voxel number for the given indexes
116 size_t n = VoxelizedSensitiveDetector::GetInstance()->GetThisVoxelNumber(i, j, k);
117
118 // Calculate only for hadrons that traveled finite range and released finite energy.
119 if ((Z >= 1) && (dX > 0) && (E > 0)) {
120 RBE* rbe = dynamic_cast<RBE*>(Manager::GetInstance()->GetQuantity("RBE"));
121 std::tuple<G4double, G4double> alpha_beta = rbe->GetHitAlphaAndBeta(E, Z);
122 fDenominator[n] += energyDeposit;
123 fAlphaNumerator[n] += std::get<0>(alpha_beta) * energyDeposit;
124 fBetaNumerator[n] += std::sqrt(std::get<1>(alpha_beta)) * energyDeposit;
125 }
126 }
127
128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
129
130 G4int RBEAccumulable::GetVerboseLevel() const
131 {
132 // Return same verbosity of RBE class
133 return Manager::GetInstance()->GetQuantity("RBE")->GetVerboseLevel();
134 }
135
136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
137
138 void RBEAccumulable::Initialize()
139 {
140 if (GetVerboseLevel() > 0) {
141 G4cout << "RBEAccumulable::Initialize(): " << G4endl;
142 }
143
144 auto voxSensDet = VoxelizedSensitiveDetector::GetInstance();
145
146 fVoxelsAlongX = voxSensDet->GetVoxelNumberAlongX();
147 fVoxelsAlongY = voxSensDet->GetVoxelNumberAlongY();
148 fVoxelsAlongZ = voxSensDet->GetVoxelNumberAlongZ();
149 fVoxels = fVoxelsAlongX * fVoxelsAlongY * fVoxelsAlongZ;
150
151 if (GetVerboseLevel() > 1) {
152 G4cout << fVoxels << " voxels." << G4endl;
153 }
154
155 fAlphaNumerator = array_type(0.0, fVoxels);
156 fBetaNumerator = array_type(0.0, fVoxels);
157 fDenominator = array_type(0.0, fVoxels);
158 fInitialized = true;
159 }
160
161 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
162
163 } // namespace RadioBio
164