| 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 TLKModel.cc
28 /// \brief Implementation of the TLKModel class
29
30 #include "TLKModel.hh"
31
32 #include "ClassifiedDamage.hh"
33 #include "Damage.hh"
34 #include "DamageClassifier.hh"
35 #include "ODESolver.hh"
36 #include <cmath>
37 #include <iostream>
38 #include <fstream>
39 #include <functional>
40
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
42
43 TLKModel::TLKModel(double pLambda1,double pLambda2, double pBeta1, double pBeta2,double pEta):
44 fLambda1(pLambda1),
45 fLambda2(pLambda2),
46 fBeta1(pBeta1),
47 fBeta2(pBeta2),
48 fEta(pEta)
49 {
50 fSingleDSBYield = 0;
51 fComplexDSBYield = 0;
52 fBpForDSB = 10;
53 fStartTime = 0.0;
54 fStopTime = 480.0;
55 fStepTime = 0.048;
56 }
57
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
59
60 std::vector<double> TLKModel::TLK_odes_system(double t,std::vector<double> y){
61
62 std::vector<double> dxdt;
63 double dxdt1 = -fLambda1*y[0]-fEta*y[0]*(y[0]+y[1]);
64 double dxdt2 = -fLambda2*y[1]-fEta*y[1]*(y[0]+y[1]);
65 double dxdt3 = fBeta1*fLambda1*y[0]+fBeta2*fLambda2*y[1]+0.25*fEta*(y[0]+y[1])*(y[0]+y[1]);
66 dxdt.push_back(dxdt1);
67 dxdt.push_back(dxdt2);
68 dxdt.push_back(dxdt3);
69 return dxdt;
70 }
71
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73
74 void TLKModel::ComputeAndSetDamageInput(std::vector<Damage> vecDamage)
75 {
76 DamageClassifier damClass;
77 auto classifiedDamage = damClass.MakeCluster(vecDamage,fBpForDSB,false);
78
79 fComplexDSBYield = damClass.GetNumComplexDSB(classifiedDamage);
80 fSingleDSBYield = damClass.GetNumDSB(classifiedDamage)-fComplexDSBYield;
81
82
83 fComplexDSBYield = fComplexDSBYield/fDose;
84 fSingleDSBYield = fSingleDSBYield/fDose;
85
86 }
87
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
89
90 double TLKModel::ComputeSF(double pDose)
91 {
92 std::vector<double> y( 3 );
93
94 y[0] = fSingleDSBYield*pDose;
95 y[1] = fComplexDSBYield*pDose;
96 y[2] = 0;
97 std::function<std::vector<double>(double,std::vector<double>)>
98 func = [this] (double t,std::vector<double> y) -> std::vector<double> {
99 return TLK_odes_system(t,y);
100 };
101 ODESolver odeSolver;
102 odeSolver.RungeKutta4(func,y,fStartTime,fStopTime,fStepTime);
103 return std::exp(-y[2]);
104 }
105
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107
108 void TLKModel::CalculateRepair(double pDoseMax, double pDeltaDose)
109 {
110 fSFCurve.clear();
111
112 for(double dose=0.;dose<=pDoseMax;dose+=pDeltaDose)
113 {
114 fSFCurve.push_back(std::make_pair(dose,ComputeSF(dose)));
115 }
116
117 }
118
119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
120
121 void TLKModel::WriteOutput(std::string pFileName)
122 {
123 std::fstream file;
124 file.open(pFileName.c_str(), std::ios_base::out);
125 //Header part
126 file <<"#============================================= TLK MODEL =============================================#\n";
127 file << " TLK Model, CalculateRepair with:\n";
128 file << "#Single DSB = " << fSingleDSBYield << " (DSB/Gy) " << "#Complex DSB = " << fComplexDSBYield << " (DSB/Gy)\n";
129 file << "#Lambda1 = " << fLambda1 << " " << "#Lambda2 = " << fLambda2 << "\n";
130 file << "#Beta1 = " << fBeta1 << " " << "#Beta2 = " << fBeta2 << "\n";
131 file << "#Eta = " << fEta << "\n";
132 file <<"#========================================================================================================#\n";
133
134 file << "Dose (Gy)\tSF\n";
135 //End Header part
136 for(int i=0;i<fSFCurve.size();i++)
137 {
138 file << fSFCurve[i].first << "\t" << fSFCurve[i].second << "\n";
139 }
140
141 file.close();
142 }
143
144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......