| Geant4 Cross Reference |
1 #!/usr/bin/env python
2 # coding: utf-8
3
4 # In[1]:
5
6
7 # To run, change the links and definitions in
8 # type in the terminal containing the file
9 # $python3 createSDD.py
10
11 # Author contact, Konstantinos Chatzipapas, ch
12
13
14 # In[2]:
15
16
17 import ROOT as root
18 import numpy as np
19
20
21 # In[3]:
22
23
24 # Define filenames
25 outputFile = "SDD_minFormatMolecularDNA.txt"
26 iMacFile = "human_cell.mac"
27 iRootFile = "molecular-dna.root"
28
29
30 # In[4]:
31
32
33 # Definitions
34 # They may be adapted in the code, if such inf
35
36 # Define if it is a simulation of a single-tra
37 doseFluence = 1
38 amountDoseFluence = 0
39
40 # Define dose rate
41 doseRate = float(0)
42
43 # Define the number of chromosomes and the len
44 nbChromo = 1
45 chromoLength = 6405.886128
46
47 # Define DNA density (BPs/um3)
48 dnaDensity = 1.2*1e07
49
50 # Define cell cycle phase (use G0 (1), G1 (2),
51 cellCyclePhase = "0, 0.0"
52
53 # Define DNA structure : whole nucleus (0), a
54 # euchromatin region (
55 # a mixed (heterochrom
56 # single DNA fiber (4)
57 # DNA plasmid (6) or a
58 dnaStructure = 4
59 nakedWet = 1 # naked (0), wet (1)
60
61 # Define if the real experiment was in-vitro (
62 vitroVivo = 0
63
64 # Define the proliferation status, quiescent (
65 proliferation = 0
66
67 # Define the microenvironment
68 temperature = 27 # degrees
69 oxygen = 0.0 # molarity
70
71 # Damage definition
72 directIndirect = 1 # direct effects only (0)
73 bpNm = 0 # BPs (0) or in nm (1)
74 bpThres = 10.0 # the distance in BPs or
75 baseLesions = -1 # This value then determi
76
77 # Default, if nothing exists in mac file
78 sourceTropus = "iso"
79
80
81 # In[5]:
82
83
84 # Initialize several parameters
85 title = "No Title Included"
86 incidentParticles = " "
87 particle = " "
88 energy = " "
89 energyUnits = " "
90 sourceShape = " "
91 sourceTropus = " "
92 sourceX = "0"
93 sourceY = "0"
94 sourceZ = "0"
95 sourceUnits = "nm"
96 targetType = " "
97 targetShape = " "
98 targetRadius = 0
99 targetX = 0
100 targetY = 0
101 targetZ = 0
102 r3 = 0
103 worldSize = 0
104 directDamageL = " "
105 directDamageU = " "
106 time = " "
107 timeUnit = " "
108
109
110
111 # In[6]:
112
113
114 with open(iMacFile, 'r') as infile:
115 for line in infile:
116 spl = line.split(" ")
117 if spl[0] == "###":
118 title = spl[1]
119 #print(title)
120 if spl[0] =="/run/beamOn":
121 incidentParticles = int(spl[1])
122 #print(incidentParticles)
123 if spl[0] =="/gps/particle":
124 spl2 = spl[1].split("\n")
125 particle = spl2[0]
126 #print(particle)
127 if spl[0] =="/gps/energy":
128 sourceType = 1
129 energyDistribution = "M, 0"
130 particleFraction = float(1)
131 energy = float(spl[1])
132 spl2 = spl[2].split("\n")
133 energyUnits = spl2[0]
134 #print(energy,energyUnits)
135 if spl[0] =="/gps/pos/shape":
136 spl2 = spl[1].split("\n")
137 sourceShape = spl2[0]
138 #print(sourceShape)
139 if spl[0] =="/gps/ang/type":
140 spl2 = spl[1].split("\n")
141 sourceTropus = spl2[0]
142 #print("type",sourceTropus)
143 if spl[0] =="/gps/pos/centre":
144 sourceX = spl[1]
145 sourceY = spl[2]
146 sourceZ = spl[3]
147 spl2 = spl[4].split("\n")
148 sourceUnits = spl2[0]
149 if spl[0] =="/chromosome/add":
150 targetType = spl[1]
151 targetShape = spl[2]
152 if targetShape == "sphere":
153 #print (spl)
154 if spl[4]=="nm\n":
155 targetRadius = float(spl[3
156 r3 = float(spl[3])*1e-09
157 #print ("r3=", r3)
158 elif spl[4]=="um\n":
159 targetRadius = float(spl[3
160 r3 = float(spl[3])*1e-06
161 #print ("r3=", r3)
162 if targetShape == "ellipse":
163 #print (spl)
164 if spl[9]=="nm":
165 targetX = float(spl[3])/10
166 targetY = float(spl[4])/10
167 targetZ = float(spl[5])/10
168 r3 = float(spl[3])*1e-09*f
169 #print ("r3=", r3)
170 elif spl[9]=="um":
171 targetX = float(spl[3])
172 targetY = float(spl[4])
173 targetZ = float(spl[5])
174 r3 = float(spl[3])*1e-06*f
175 #print ("r3=", r3)
176 #print(targetType,targetShape,targ
177 if spl[0] =="/world/worldSize":
178 worldSize = float(spl[1])/2
179 if spl[2] =="nm\n":
180 worldSize = worldSize/1000
181 #print(worldSize)
182 if spl[0] =="/dnadamage/directDamageLo
183 directDamageL = spl[1]
184 #print(directDamageL)
185 if spl[0] =="/dnadamage/directDamageUp
186 directDamageU = spl[1]
187 #print(directDamageU)
188 if spl[0] =="/scheduler/endTime":
189 time = spl[1]
190 spl2 = spl[2].split("\n")
191 timeUnit = spl2[0]
192 #print(time, timeUnit)
193
194
195 # In[7]:
196
197
198 # Creating SDD file
199 with open(outputFile, 'w') as ofile:
200 # Starting with header part of SDD file
201 ofile.write("\nSDD version, SDDv1.0;\n")
202 ofile.write("Software, MolecularDNA;\n")
203 ofile.write("Author contact, Konstantinos
204 "14/03/2022;\n")
205 ofile.write("***Important information*****
206 "To provide some extra informa
207 "DSB+ and DSB++, the last colu
208 "of 4, 5 that correspond to DS
209 "*****************************
210 #Adjust description
211 ofile.write("Simulation Details, "+title+"
212 ofile.write("Source, Monoenergetic "+sourc
213 ofile.write(" centered at "+sourceX+" "+so
214 ofile.write(" of a cell nucleus, containin
215 ofile.write(str(energy)+" "+energyUnits+";
216 ######
217 ofile.write("Source type, "+str(sourceType
218 ofile.write("Incident particles, "+str(inc
219 ofile.write("Mean particle energy, "+str(e
220 ofile.write("Energy distribution, "+energy
221 ofile.write("Particle fraction, "+str(part
222
223 f = root.TFile(iRootFile)
224
225 eVtoJ = 1.60218e-19
226
227 mass = 997 * 4 * 3.141592 * r3 / 3 # d
228
229 acc_edep = 0
230 dose = 0
231 nbEntries = 0
232 ffTree = f.Get("tuples/chromosome_hits")
233 nbEntries += ffTree.GetEntries()
234 for ev in ffTree:
235 acc_edep += (ev.e_chromosome_kev + ev.
236
237
238 # Calculate the absorbed dose
239 amountDoseFluence = acc_edep * eVtoJ / mas
240
241
242
243 ofile.write("Dose or fluence, "+str(doseFl
244 ofile.write("Dose rate, "+str(doseRate)+";
245 ofile.write("Irradiation target, Simple fr
246 if targetShape == "sphere":
247 ofile.write("radius "+str(targetRadius
248 else: #targetShape == "ellipse":
249 ofile.write("dimensions "+str(targetX)
250 ofile.write("Volumes, 0,"+str(worldSize)+"
251 if targetShape == "sphere":
252 ofile.write(" 1,"+str(targetRadius)+",
253 elif targetShape == "ellipse":
254 ofile.write(" 1,"+str(targetX)+","+str
255 else:
256 ofile.write(" 0,"+str(targetX)+","+str
257 ofile.write("Chromosome sizes, "+str(nbChr
258 ofile.write("DNA Density, "+str(dnaDensity
259 ofile.write("Cell Cycle Phase, "+cellCycle
260 ofile.write("DNA Structure, "+str(dnaStruc
261 ofile.write("In vitro / in vivo, "+str(vit
262 ofile.write("Proliferation status, "+str(p
263 ofile.write("Microenvironment, "+str(tempe
264 ofile.write("Damage definition, "+str(dire
265 ofile.write("Time, "+str(time)+" "+timeUni
266
267 number = 0
268 gTree = f.Get("tuples/primary_source")
269 number += gTree.GetEntries()
270
271 nEntries = 0
272 fTree = f.Get("tuples/damage")
273 nEntries += fTree.GetEntries()
274 #print (nEntries)
275
276 ofile.write("Damage and primary count, "+s
277 ofile.write("Data entries, 1, 1, 0, 0, 0,
278 ofile.write("Data was generated in minimal
279 ofile.write("Please modify description to
280 ofile.write("\n***EndOfHeader***;\n\n")
281
282 ##############################################
283 # Writing the data part of SDD file
284 currentEvent = 0
285 Primaryflag = True
286 Eventflag = False
287
288 for entry in fTree:
289 if entry.Event!=currentEvent:
290 Eventflag = True
291 currentEvent = entry.Event
292 #print(entry.SourceClassification)
293 DSB = 0
294 if (entry.TypeClassification == "DSB")
295 if (entry.SourceClassification ==
296 DSB = 1
297 elif (entry.SourceClassification =
298 DSB = 2
299 elif (entry.SourceClassification =
300 DSB = 3
301 #print (entry.TypeClassificati
302 elif (entry.TypeClassification == "DSB
303 DSB = 4
304 elif (entry.TypeClassification == "DSB
305 DSB = 5
306 #if (entry.SourceClassification ==
307 #print (entry.TypeClassificati
308 else:
309 DSB = 0
310
311 #print(DSB)
312 if Primaryflag:
313 #print("True")
314 ofile.write("2, "+str(entry.Event)
315 ofile.write(str(entry.Position_x_u
316 ofile.write(str(entry.BaseDamage)+
317 Primaryflag = False
318 Eventflag = False
319 else:
320 if Eventflag:
321 ofile.write("1, "+str(entry.Ev
322 ofile.write(str(entry.Position
323 ofile.write(str(entry.BaseDama
324 Eventflag = False
325 else:
326 ofile.write("0, "+str(entry.Ev
327 ofile.write(str(entry.Position
328 ofile.write(str(entry.BaseDama
329
330
331 # In[8]:
332
333 print("Output file: ", outputFile)
334
335
336