| Geant4 Cross Reference |
1 #!/usr/bin/env python
2 # coding: utf-8
3
4 # To run, change the links and definitions in
5 # type in the terminal containing the file
6 # $python3 molecularDNAsurvival.py
7
8 # Authors:
9 # K. Chatzipapas (*), D. Sakata, H. N. Tran, M
10 # (*) contact: k.chatzipapas@yahoo.com
11
12 # Authors' note : This is a very early version
13 # Needs to optimized for each
14 # version will be released in
15
16
17 # Define filenames
18 outputFile = "molecularDNAsurvival.txt"
19 iRootFile = "../molecular-dna.root"
20 print("\nInput file:",iRootFile,"\n")
21
22 # Define cell parameters
23 r3 = 7100*1e-09 * 2500*1e-09 * 7100*1e-09 # a
24 NBP = 6405886128 # bp / Length of the DNA chai
25 mass = 997 * 4 * 3.141592 * r3 / 3 # water
26
27 # Integration parameters
28 t0 = 0.0 # Starting time point (dimensionles
29 t1 = 336 # Final time point (dimensionless)
30 dt = 0.001 # Intergration time step (dimension
31 numpoints = int( (t1 - t0) / dt )
32
33 # Survival model Parameters
34 # THESE PARAMETERS HAVE TO BE DEFINED FOR PROP
35 # AS DESCRIBED IN (Stewart RD. Radiat Res. 200
36 maxDose = 9
37 DR1 = 720 #Gy/h SF-dose
38 bp = 1 #Gbp
39 #DR2 = 600. #Gy/h SF-time
40 #DR3 = 12000. #Gy/h FAR
41
42 # usual value of gamma 0.25 (changes the end v
43 gamma = 0.19 #other published values: 0.18932
44 # Lamb1 the end of the curve Lamb2 is connecte
45 Lamb1 = 0.10 #other published values: 0.01259
46 Lamb2 = 1 #other published values: 33062.9
47 #
48 Beta1 = 0 #other published values: 0.01932
49 Beta2 = 0
50 # Defines the curvature as well as the end val
51 Eta = 0.00005 #other published values: #7.5
52 # name of the cell (used in the output)
53 cell = "test"
54
55 # ONLY ADVANCED USERS AFTER THIS POINT
56 ##############################################
57 import ROOT as root
58 import numpy as np
59 import matplotlib.pyplot as plt
60 from scipy.integrate import ode
61 import threading
62
63 # Definition of the survival model
64 def SF_system( T_, DR_, Y_, Sig1_, Sig2_, gamm
65
66 global T
67 global DR
68 global Y
69 global Sig1
70 global Sig2
71 global gamma
72 global Lamb1
73 global Lamb2
74 global Beta1
75 global Beta2
76 global Eta
77
78 T = T_
79 DR = DR_
80 Y = Y_
81 Sig1 = Sig1_
82 Sig2 = Sig2_
83 gamma = gamma_
84 Lamb1 = Lamb1_
85 Lamb2 = Lamb2_
86 Beta1 = Beta1_
87 Beta2 = Beta2_
88 Eta = Eta_
89
90
91 return T, DR, Y, Sig1, Sig2, gamma, Lamb1,
92
93 def SF_function( t, x, p ):
94 if(t<=T):
95 #print("t<=T")
96 dx[0] = DR*Y*Sig1 - Lamb1*x[0] - Eta*x
97 dx[1] = DR*Y*Sig2 - Lamb2*x[1] - Eta*x
98 dx[2] = Beta1*Lamb1*x[0] + Beta2*Lamb2
99 dx[3] = (1.-Beta1)*Lamb1*x[0] + (1.-Be
100 else:
101 dx[0] = - Lamb1*x[0] - Eta*x[0]*(x[0]
102 dx[1] = - Lamb2*x[1] - Eta*x[1]*(x[0]
103 dx[2] = Beta1*Lamb1*x[0] + Beta2*Lamb2
104 dx[3] = (1.-Beta1)*Lamb1*x[0] + (1.-Be
105 return dx
106
107
108 DSBBPID = []
109 repDSB = 0
110 irrDSB = 0
111 acc_edep = 0
112 dose = 0
113 i = 0
114 eVtoJ = 1.60218e-19
115
116
117 # Analyze root files
118 f = root.TFile(iRootFile)
119
120 fTree = f.Get("tuples/damage")
121 for e in fTree:
122 if (e.TypeClassification=="DSB" or e.TypeC
123 DSBBPID.append((e.Event, e.BasePair))
124 i += 1
125
126 gTree = f.Get("tuples/classification")
127 for e in gTree:
128 repDSB += e.DSB
129 irrDSB += e.DSBp + 2*e.DSBpp
130
131 ffTree = f.Get("tuples/chromosome_hits")
132 for ev in ffTree:
133 acc_edep += (ev.e_chromosome_kev + ev.e_dn
134
135
136 # Calculate the total number of DSBs
137 totalDSB = repDSB + irrDSB
138 # Calculate the absorbed dose
139 dose = acc_edep * eVtoJ / mass # Dose in Gy
140
141 #print ("r3=", r3)
142 #print("Dose :", dose, "\nTotal DSB/Gy :", tot
143
144 NBP_of_cell = 4682000000 # This is a constan
145 NBPcell = NBP/NBP_of_cell
146 totalDSBperGyCell = totalDSB/NBPcell/dose
147 #totalDSBperGyCell = totalDSB/dose
148 print ("Dose :", dose, "\nTotal DSB/Gy/Cell :"
149
150 kFactor = 1 # This factor is used to normaliz
151 S1 = repDSB/NBPcell/dose/kFactor
152 S2 = irrDSB/NBPcell/dose/kFactor
153 print("Repairable DSBs/Gy/cell :", S1, "\nIrre
154
155
156 # In[ ]:
157
158
159 T = 0
160 DR = 0
161 Y = 0
162 kozi = False
163 sum = [0]*(maxDose+1)
164 outputfile = [""] * (maxDose+1)
165
166 # Create the time samples for the output of th
167 # I use a large number of points, only because
168 # a plot of the solution that looks nice.
169 #t = [(t1-t0) * float(i) / (numpoints - 1) for
170 t = np.linspace(t0, t1, numpoints)
171
172 dx = np.zeros(4)
173 x = np.zeros(4)
174 dsbs = 0.0
175 Sig1 = S1
176 Sig2 = S2
177
178 filename = "txt/molecularDNAsurvival_"+cell+"
179 outputfile1 = open(filename,"w")
180 outputfile1.write("Dose \t Survival Fraction\n
181
182 for j in range(maxDose+1):
183 #j=1
184 sol = []
185 time = []
186 dose = j
187 expEndTime = dose/DR1 #h
188
189 #////double StopTime = (expEndTime+ n/min
190
191 arguments = SF_system( expEndTime, DR1,bp,
192 p = [ T, DR, Y, Sig1, Sig2, gamma, Lamb1,
193
194 Stepper = "dopri5" # dopri5,vode,dop853
195
196 if (j==0):
197 solver = ode(SF_function).set_integrat
198 print("#*#*#*#*#*# Ignore this warning
199 else:
200 solver = ode(SF_function).set_integrat
201 solver.set_initial_value( x, t0 ).set_f_pa
202 #print(solver.t,solver.y)
203
204 # Repeatedly call the `integrate` method t
205 # solution to time t+dt, and save the solu
206 while solver.t < t1: # and solver.successf
207 solver.integrate(solver.t+dt)
208
209 sum[j] = solver.y[2]
210 print("Survival probability of cell ("+cel
211 outputfile1.write(str(dose) +"\t"+ str(np.
212
213 c = np.linspace(0, maxDose, num=maxDose+1)
214 sum=np.array(sum)
215 outputfile1.close()
216
217 print("\nCalculation concluded. Check the fina
218
219 plt.yscale('log')
220 plt.xlim([0, 8])
221 plt.ylim([1e-2, 1])
222 plt.plot(c,np.exp(-sum[:]))
223 plt.show()
224
225 print("\nThank you for using molecularDNAsurvi
226