| 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 this script and
8 # type in the terminal containing the file
9 # $python3 createSDD.py
10
11 # Author contact, Konstantinos Chatzipapas, chatzipa@cenbg.in2p3.fr, 14/03/2022
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 information exists in the simulation definition file
35
36 # Define if it is a simulation of a single-track irradiation (0), a delivered dose (1) or a fluence (2)
37 doseFluence = 1
38 amountDoseFluence = 0
39
40 # Define dose rate
41 doseRate = float(0)
42
43 # Define the number of chromosomes and the length in MBp
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), S (3), G2 (4) and M (5), together with percentage, e.g [3,0.7] indicates a cell 70% of the way through S phase)
51 cellCyclePhase = "0, 0.0"
52
53 # Define DNA structure : whole nucleus (0), a heterochromatin region (1),
54 # euchromatin region (2),
55 # a mixed (heterochromatin and euchromatin) region (3),
56 # single DNA fiber (4), DNA wrapped around a single histone (5),
57 # DNA plasmid (6) or a simple circular (7) or straight (8) DNA section.
58 dnaStructure = 4
59 nakedWet = 1 # naked (0), wet (1)
60
61 # Define if the real experiment was in-vitro (0) or in-vivo (1)
62 vitroVivo = 0
63
64 # Define the proliferation status, quiescent (0) or proliferating (1)
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) or including chemistry (1)
73 bpNm = 0 # BPs (0) or in nm (1)
74 bpThres = 10.0 # the distance in BPs or nm between backbone lesions that are considered DSBs
75 baseLesions = -1 # This value then determines the distance (in BP or nm) beyond the outer backbone damages where base damages are also stored in the same site (float).
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])/1000
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])/1000
166 targetY = float(spl[4])/1000
167 targetZ = float(spl[5])/1000
168 r3 = float(spl[3])*1e-09*float(spl[4])*1e-09*float(spl[5])*1e-09
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*float(spl[4])*1e-06*float(spl[5])*1e-06
175 #print ("r3=", r3)
176 #print(targetType,targetShape,targetRadius,targetX,targetY,targetZ)
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/directDamageLower":
183 directDamageL = spl[1]
184 #print(directDamageL)
185 if spl[0] =="/dnadamage/directDamageUpper":
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 Chatzipapas, chatzipa@cenbg.in2p3.fr, "
204 "14/03/2022;\n")
205 ofile.write("***Important information*********************************************\n"
206 "To provide some extra information on the quantification of DSB, "
207 "DSB+ and DSB++, the last column of data section, includes the values "
208 "of 4, 5 that correspond to DSB+ and DSB++ respectively;\n"
209 "*********************************************************************\n")
210 #Adjust description
211 ofile.write("Simulation Details, "+title+". DNA damages from direct and indirect effects;\n")
212 ofile.write("Source, Monoenergetic "+sourceShape+" "+particle+" "+sourceTropus+"tropic,")
213 ofile.write(" centered at "+sourceX+" "+sourceY+" "+sourceZ+" "+sourceUnits)
214 ofile.write(" of a cell nucleus, containing a free DNA segment. Energy: ")
215 ofile.write(str(energy)+" "+energyUnits+";\n")
216 ######
217 ofile.write("Source type, "+str(sourceType)+";\n") # Needs improvement
218 ofile.write("Incident particles, "+str(incidentParticles)+";\n")
219 ofile.write("Mean particle energy, "+str(energy)+" MeV;\n")
220 ofile.write("Energy distribution, "+energyDistribution+";\n") # Needs improvement
221 ofile.write("Particle fraction, "+str(particleFraction)+";\n") # Needs improvement
222
223 f = root.TFile(iRootFile)
224
225 eVtoJ = 1.60218e-19
226
227 mass = 997 * 4 * 3.141592 * r3 / 3 # density * 4/3 * pi * r3
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.e_dna_kev) *1e3 # eV
236
237
238 # Calculate the absorbed dose
239 amountDoseFluence = acc_edep * eVtoJ / mass # Dose in Gy
240
241
242
243 ofile.write("Dose or fluence, "+str(doseFluence)+", "+str(amountDoseFluence)+";\n") # Needs improvement
244 ofile.write("Dose rate, "+str(doseRate)+";\n")
245 ofile.write("Irradiation target, Simple free DNA fragment in a "+targetShape+" with ")
246 if targetShape == "sphere":
247 ofile.write("radius "+str(targetRadius)+" um;\n")
248 else: #targetShape == "ellipse":
249 ofile.write("dimensions "+str(targetX)+" "+str(targetY)+" "+str(targetZ)+" um;\n")
250 ofile.write("Volumes, 0,"+str(worldSize)+","+str(worldSize)+","+str(worldSize)+","+str(-worldSize)+","+str(-worldSize)+","+str(-worldSize)+",")
251 if targetShape == "sphere":
252 ofile.write(" 1,"+str(targetRadius)+","+str(targetRadius)+","+str(targetRadius)+";\n")
253 elif targetShape == "ellipse":
254 ofile.write(" 1,"+str(targetX)+","+str(targetY)+","+str(targetZ)+";\n")
255 else:
256 ofile.write(" 0,"+str(targetX)+","+str(targetY)+","+str(targetZ)+","+str(+targetX)+","+str(+targetY)+","+str(+targetZ)+";\n")
257 ofile.write("Chromosome sizes, "+str(nbChromo)+", "+str(chromoLength)+";\n")
258 ofile.write("DNA Density, "+str(dnaDensity)+";\n")
259 ofile.write("Cell Cycle Phase, "+cellCyclePhase+";\n")
260 ofile.write("DNA Structure, "+str(dnaStructure)+", "+str(nakedWet)+";\n")
261 ofile.write("In vitro / in vivo, "+str(vitroVivo)+";\n")
262 ofile.write("Proliferation status, "+str(proliferation)+";\n")
263 ofile.write("Microenvironment, "+str(temperature)+", "+str(oxygen)+";\n")
264 ofile.write("Damage definition, "+str(directIndirect)+", "+str(bpNm)+", "+str(bpThres)+", "+str(baseLesions)+", "+str(directDamageL)+";\n")
265 ofile.write("Time, "+str(time)+" "+timeUnit+";\n")
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, "+str(nEntries)+", "+str(number)+";\n")
277 ofile.write("Data entries, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0;\n")
278 ofile.write("Data was generated in minimal output format and used as an example. ")
279 ofile.write("Please modify description to match different conditions;\n")
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 == "DSBd"):
296 DSB = 1
297 elif (entry.SourceClassification == "DSBi"):
298 DSB = 2
299 elif (entry.SourceClassification == "DSBh" or entry.SourceClassification == "DSBm"):
300 DSB = 3
301 #print (entry.TypeClassification, entry.SourceClassification)
302 elif (entry.TypeClassification == "DSB+"):
303 DSB = 4
304 elif (entry.TypeClassification == "DSB++"):
305 DSB = 5
306 #if (entry.SourceClassification == "DSBh" or entry.SourceClassification == "DSBm"):
307 #print (entry.TypeClassification, entry.SourceClassification)
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_um)+", "+str(entry.Position_y_um)+", "+str(entry.Position_z_um)+"; ")
316 ofile.write(str(entry.BaseDamage)+", "+str(entry.StrandDamage)+", "+str(DSB)+";\n")
317 Primaryflag = False
318 Eventflag = False
319 else:
320 if Eventflag:
321 ofile.write("1, "+str(entry.Event)+"; ")
322 ofile.write(str(entry.Position_x_um)+", "+str(entry.Position_y_um)+", "+str(entry.Position_z_um)+"; ")
323 ofile.write(str(entry.BaseDamage)+", "+str(entry.StrandDamage)+", "+str(DSB)+";\n")
324 Eventflag = False
325 else:
326 ofile.write("0, "+str(entry.Event)+"; ")
327 ofile.write(str(entry.Position_x_um)+", "+str(entry.Position_y_um)+", "+str(entry.Position_z_um)+"; ")
328 ofile.write(str(entry.BaseDamage)+", "+str(entry.StrandDamage)+", "+str(DSB)+";\n")
329
330
331 # In[8]:
332
333 print("Output file: ", outputFile)
334
335
336