repairmodels/src/TLKModel.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 TLKModel.cc
0028 /// \brief Implementation of the TLKModel class
0029
0030 #include "TLKModel.hh"
0031
0032 #include "ClassifiedDamage.hh"
0033 #include "Damage.hh"
0034 #include "DamageClassifier.hh"
0035 #include "ODESolver.hh"
0036 #include <cmath>
0037 #include <iostream>
0038 #include <fstream>
0039 #include <functional>
0040
0041 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0042
0043 TLKModel::TLKModel(double pLambda1,double pLambda2, double pBeta1, double pBeta2,double pEta):
0044 fLambda1(pLambda1),
0045 fLambda2(pLambda2),
0046 fBeta1(pBeta1),
0047 fBeta2(pBeta2),
0048 fEta(pEta)
0049 {
0050     fSingleDSBYield = 0;
0051     fComplexDSBYield = 0;
0052     fBpForDSB = 10;
0053     fStartTime = 0.0;
0054     fStopTime = 480.0;
0055     fStepTime = 0.048;
0056 }
0057
0058 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0059
0060 std::vector<double> TLKModel::TLK_odes_system(double t,std::vector<double> y){
0061
0062     std::vector<double> dxdt;
0063     double dxdt1 = -fLambda1*y[0]-fEta*y[0]*(y[0]+y[1]);
0064     double dxdt2 = -fLambda2*y[1]-fEta*y[1]*(y[0]+y[1]);
0065     double dxdt3 = fBeta1*fLambda1*y[0]+fBeta2*fLambda2*y[1]+0.25*fEta*(y[0]+y[1])*(y[0]+y[1]);
0066     dxdt.push_back(dxdt1);
0067     dxdt.push_back(dxdt2);
0068     dxdt.push_back(dxdt3);
0069     return dxdt;
0070 }
0071
0072 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0073
0074 void TLKModel::ComputeAndSetDamageInput(std::vector<Damage> vecDamage)
0075 {
0076     DamageClassifier damClass;
0077     auto classifiedDamage = damClass.MakeCluster(vecDamage,fBpForDSB,false);
0078
0079     fComplexDSBYield = damClass.GetNumComplexDSB(classifiedDamage);
0080     fSingleDSBYield = damClass.GetNumDSB(classifiedDamage)-fComplexDSBYield;
0081
0082
0083     fComplexDSBYield = fComplexDSBYield/fDose;
0084     fSingleDSBYield = fSingleDSBYield/fDose;
0085
0086 }
0087
0088 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0089
0090 double TLKModel::ComputeSF(double pDose)
0091 {
0092     std::vector<double> y( 3 );
0093
0094     y[0] = fSingleDSBYield*pDose;
0095     y[1] = fComplexDSBYield*pDose;
0096     y[2] = 0;
0097     std::function<std::vector<double>(double,std::vector<double>)>
0098     func = [this] (double t,std::vector<double> y) -> std::vector<double> {
0099         return TLK_odes_system(t,y);
0100     };
0101     ODESolver odeSolver;
0102     odeSolver.RungeKutta4(func,y,fStartTime,fStopTime,fStepTime);
0103     return std::exp(-y[2]);
0104 }
0105
0106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0107
0108 void TLKModel::CalculateRepair(double pDoseMax, double pDeltaDose)
0109 {
0110     fSFCurve.clear();
0111
0112     for(double dose=0.;dose<=pDoseMax;dose+=pDeltaDose)
0113     {
0114         fSFCurve.push_back(std::make_pair(dose,ComputeSF(dose)));
0115     }
0116
0117 }
0118
0119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0120
0121 void TLKModel::WriteOutput(std::string pFileName)
0122 {
0123     std::fstream file;
0124     file.open(pFileName.c_str(), std::ios_base::out);
0125     //Header part
0126     file <<"#============================================= TLK  MODEL =============================================#\n";
0127     file << "                                  TLK Model, CalculateRepair with:\n";
0128     file << "#Single DSB  = " << fSingleDSBYield   << " (DSB/Gy)                     " << "#Complex DSB    = " << fComplexDSBYield << " (DSB/Gy)\n";
0129     file << "#Lambda1  = " << fLambda1     << "                     " << "#Lambda2    = " << fLambda2 << "\n";
0130     file << "#Beta1    = " << fBeta1       << "                  " << "#Beta2      = " << fBeta2 << "\n";
0131     file << "#Eta     = " << fEta         << "\n";
0132     file <<"#========================================================================================================#\n";
0133
0134     file << "Dose (Gy)\tSF\n";
0135     //End Header part
0136     for(int i=0;i<fSFCurve.size();i++)
0137     {
0138         file << fSFCurve[i].first << "\t" << fSFCurve[i].second << "\n";
0139     }
0140
0141     file.close();
0142 }
0143
0144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......