| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Authors: J. Naoki D. Kondo (UCSF, US) : 10/
27 // J. Ramos-Mendez and B. Faddegon (U
28 //
29 /// \file PhysGeoImport.cc
30 /// \brief Implementation of the plasmid load
31
32 #include "PhysGeoImport.hh"
33
34 #include "G4DNAChemistryManager.hh"
35 #include "G4Ellipsoid.hh"
36 #include "G4VPhysicalVolume.hh"
37 #include "Randomize.hh"
38
39 //....oooOO0OOooo........oooOO0OOooo........oo
40
41 PhysGeoImport::PhysGeoImport() {}
42
43 //....oooOO0OOooo........oooOO0OOooo........oo
44
45 G4LogicalVolume* PhysGeoImport::CreateLogicVol
46 {
47 G4NistManager* man = G4NistManager::Instance
48 fEnvelopeWater = man->FindOrBuildMaterial("G
49 fpWater =
50 man->BuildMaterialWithNewDensity("G4_WATER
51
52 ReadFile(fileName);
53
54 G4double des1 = 1.1344640137963142;
55 G4double des2 = des1 + (CLHEP::pi * .5);
56 G4double ang = 0.6283185307179586;
57 G4double bet1 = 0.6283185307179586 * 2;
58 G4double posi = 1.0471975511965976;
59 G4double sep = .1 * angstrom;
60 // Geometries Sizes
61 G4double PxRs = 2.9389169420478556 * angstro
62 G4double PyRs = 2.9389169420478556 * angstro
63 G4double PzRs = 2.9389169420478556 * angstro
64 G4double PxYs = 2.7 * angstrom;
65 G4double PyYs = 2.7 * angstrom;
66 G4double PzYs = 2.7 * angstrom;
67 G4double PxBp = 2.45 * angstrom;
68 G4double PyBp = 2.45 * angstrom;
69 G4double PzBp = 2.45 * angstrom;
70
71 G4double xin = -170 * angstrom;
72 G4double yin = -170 * angstrom;
73 G4double zin = -170 * angstrom;
74 G4double xfn = 170 * angstrom;
75 G4double yfn = 170 * angstrom;
76 G4double zfn = 170 * angstrom;
77
78 G4int nVertex = fVertexes.size();
79
80 // Envelope
81
82 std::string boxNameSolid = fGeoName + "_soli
83 G4Box* box_solid =
84 new G4Box(boxNameSolid, 0.5 * (fXMax - fXM
85 0.5 * (fYMax - fYMin) + 0.5 * 3.
86
87 G4String boxNameLogic = fGeoName + "_logic";
88 G4LogicalVolume* box_logic =
89 new G4LogicalVolume(box_solid, fEnvelopeWa
90
91 // Desoxyribose
92
93 G4Ellipsoid* RSolidSugar = new G4Ellipsoid("
94 G4LogicalVolume* RSugar = new G4LogicalVolum
95 G4VisAttributes* MyVisAtt_Rs = new G4VisAttr
96 MyVisAtt_Rs->SetForceSolid(true);
97 RSugar->SetVisAttributes(MyVisAtt_Rs);
98
99 // Phosphoric Acid
100
101 G4Ellipsoid* YSolidSugar = new G4Ellipsoid("
102
103 G4LogicalVolume* YSugar = new G4LogicalVolum
104 G4VisAttributes* MyVisAtt_Ys = new G4VisAttr
105 MyVisAtt_Ys->SetForceSolid(true);
106 YSugar->SetVisAttributes(MyVisAtt_Ys);
107
108 // Base Pairs
109
110 G4Ellipsoid* Base1a = new G4Ellipsoid("BaseP
111 G4LogicalVolume* BaseP1a = new G4LogicalVolu
112 G4VisAttributes* MyVisAtt_Bp1a = new G4VisAt
113 MyVisAtt_Bp1a->SetForceSolid(true);
114 BaseP1a->SetVisAttributes(MyVisAtt_Bp1a);
115
116 G4Ellipsoid* Base1b = new G4Ellipsoid("BaseP
117 G4LogicalVolume* BaseP1b = new G4LogicalVolu
118
119 G4VisAttributes* MyVisAtt_Bp1b = new G4VisAt
120 MyVisAtt_Bp1b->SetForceSolid(true);
121 BaseP1b->SetVisAttributes(MyVisAtt_Bp1b);
122
123 G4int index = 0;
124 G4double cAngle = 0;
125 G4double pi = CLHEP::pi;
126 for (int vertex = 0; vertex < nVertex - 1; v
127 xin = fVertexes[vertex][0] - fOffsetX;
128 yin = fVertexes[vertex][1] - fOffsetY;
129 zin = fVertexes[vertex][2] - fOffsetZ;
130 xfn = fVertexes[vertex + 1][0] - fOffsetX;
131 yfn = fVertexes[vertex + 1][1] - fOffsetY;
132 zfn = fVertexes[vertex + 1][2] - fOffsetZ;
133
134 G4double phi0 =
135 std::atan2(zfn - zin, std::sqrt(((xfn -
136 G4double theta0 = std::atan2(xfn - xin, yf
137 G4double lenght = std::sqrt(((xfn - xin) *
138 + ((zfn - zin)
139 G4double dl = 1.0 / (lenght / (3.4 * angst
140 G4int nChain = (fVertexes[vertex] - fVerte
141
142 for (G4int nseg = 0; nseg < nChain; nseg++
143 cAngle += ang;
144
145 G4double theta = cAngle;
146 G4double x1 = 0;
147 G4double y1 = 0;
148 G4double z1 = ((2 * PzBp) + PzRs + sep);
149
150 G4double x2 = 0;
151 G4double y2 = 0;
152 G4double z2 = ((2 * PzBp) + PzRs + sep);
153
154 G4ThreeVector plus2 = G4ThreeVector(0, 0
155 plus2.rotateX(-des1);
156 plus2.rotateZ(-posi);
157
158 G4ThreeVector plus2alt = G4ThreeVector(0
159 plus2alt.rotateX(-des1);
160 plus2alt.rotateZ(posi);
161
162 G4double x3 = 0;
163 G4double y3 = 0;
164 G4double z3 = PzBp + sep;
165
166 G4ThreeVector position1i = G4ThreeVector
167 G4ThreeVector position2i = G4ThreeVector
168 G4ThreeVector position2ialt = G4ThreeVec
169 G4ThreeVector position3i = G4ThreeVector
170
171 position1i.rotateY(theta);
172 position2i.rotateY(theta);
173 position2ialt.rotateY(theta);
174 position3i.rotateY(theta);
175
176 G4double x = dl * nseg * (xfn - xin) + x
177 G4double y = dl * nseg * (yfn - yin) + y
178 G4double z = dl * nseg * (zfn - zin) + z
179
180 position1i.rotateX(phi0);
181 position2i.rotateX(phi0);
182 position2ialt.rotateX(phi0);
183 position3i.rotateX(phi0);
184
185 position1i.rotateZ(-theta0);
186 position2i.rotateZ(-theta0);
187 position2ialt.rotateZ(-theta0);
188 position3i.rotateZ(-theta0);
189
190 G4double yrot1 = theta;
191 G4double xrot1 = -des1;
192 G4RotationMatrix rotm1 = G4RotationMatri
193 rotm1.rotateX(xrot1);
194 rotm1.rotateZ(-posi);
195 rotm1.rotateY(yrot1);
196 rotm1.rotateX(phi0);
197 rotm1.rotateZ(-theta0);
198 G4ThreeVector position1 = position1i + G
199 G4Transform3D transform1(rotm1, position
200
201 G4double yrot1alt = theta + pi;
202 G4double xrot1alt = des1;
203 G4RotationMatrix rotm1alt = G4RotationMa
204 rotm1alt.rotateX(xrot1alt);
205 rotm1alt.rotateZ(-posi);
206 rotm1alt.rotateY(yrot1alt);
207 rotm1alt.rotateX(phi0);
208 rotm1alt.rotateZ(-theta0);
209 G4ThreeVector position1alt = -position1i
210 G4Transform3D transform1alt(rotm1alt, po
211
212 G4double yrot2 = theta;
213 G4double xrot2 = -des2;
214 G4RotationMatrix rotm2 = G4RotationMatri
215 rotm2.rotateX(xrot2);
216 rotm2.rotateY(yrot2 - bet1 + 0.872664625
217 rotm2.rotateX(phi0);
218 rotm2.rotateZ(-theta0);
219 G4ThreeVector position2 = position2i + G
220 G4Transform3D transform2(rotm2, position
221
222 G4double yrot2alt = theta + pi;
223 G4double xrot2alt = des2;
224 G4RotationMatrix rotm2alt = G4RotationMa
225 rotm2alt.rotateX(xrot2alt);
226 rotm2alt.rotateY(yrot2alt + bet1 - 0.872
227 rotm2alt.rotateX(phi0);
228 rotm2alt.rotateZ(-theta0);
229 G4ThreeVector position2alt = position2ia
230 G4Transform3D transform2alt(rotm2alt, po
231
232 G4double yrot3 = theta;
233 G4RotationMatrix rotm3 = G4RotationMatri
234 rotm3.rotateX(-pi / 2);
235 rotm3.rotateZ(-ang);
236 rotm3.rotateY(yrot3);
237 rotm3.rotateX(phi0);
238 rotm3.rotateZ(-theta0);
239 G4ThreeVector position3 = position3i + G
240 G4Transform3D transform3(rotm3, position
241
242 G4double yrot3alt = theta + pi;
243 G4RotationMatrix rotm3alt = G4RotationMa
244 rotm3alt.rotateX(pi / 2);
245 rotm3alt.rotateZ(-ang);
246 rotm3alt.rotateY(yrot3alt);
247 rotm3alt.rotateX(phi0);
248 rotm3alt.rotateZ(-theta0);
249 G4ThreeVector position3alt = -position3i
250 G4Transform3D transform3alt(rotm3alt, po
251
252 new G4PVPlacement(transform1, RSugar, "d
253
254 new G4PVPlacement(transform1alt, RSugar,
255
256 new G4PVPlacement(transform3, BaseP1a, "
257
258 new G4PVPlacement(transform3alt, BaseP1a
259
260 new G4PVPlacement(transform2, YSugar, "p
261
262 new G4PVPlacement(transform2alt, YSugar,
263
264 G4ThreeVector Deoxy1 = position1;
265 G4ThreeVector Deoxy2 = position1alt;
266
267 fSampleDNAPositions.push_back(Deoxy1);
268 fSampleDNAPositions.push_back(Deoxy2);
269 fSampleDNANames.push_back("Deoxyribose")
270 fSampleDNANames.push_back("Deoxyribose")
271 fSampleDNADetails.push_back({-1, index,
272 fSampleDNADetails.push_back({-1, index,
273
274 index++;
275 }
276 }
277
278 return box_logic;
279 }
280
281 //....oooOO0OOooo........oooOO0OOooo........oo
282
283 void PhysGeoImport::ReadFile(G4String fileName
284 {
285 G4double x, y, z;
286 fXMin = 1 * mm, fYMin = 1 * mm, fZMin = 1 *
287 fXMax = -1 * mm, fYMax = -1 * mm, fZMax = -1
288 std::ifstream plasmidFile(fileName);
289 while (true) {
290 plasmidFile >> x >> y >> z;
291 if (!plasmidFile.good()) break;
292 x *= nm;
293 y *= nm;
294 z *= nm;
295 fVertexes.push_back(G4ThreeVector(x, y, z)
296 if (fXMin > x) fXMin = x;
297 if (fXMax < x) fXMax = x;
298 if (fYMin > y) fYMin = y;
299 if (fYMax < y) fYMax = y;
300 if (fZMin > z) fZMin = z;
301 if (fZMax < z) fZMax = z;
302 }
303 plasmidFile.close();
304 fOffsetX = (fXMin + fXMax) * 0.5;
305 fOffsetY = (fYMin + fYMax) * 0.5;
306 fOffsetZ = (fZMin + fZMax) * 0.5;
307
308 std::vector<G4ThreeVector> VertRed;
309 for (size_t i = 0; i < fVertexes.size(); i++
310 if (i % 15 == 0) VertRed.push_back(fVertex
311 }
312
313 VertRed.push_back(fVertexes[0]);
314
315 fVertexes = VertRed;
316 }
317
318 //....oooOO0OOooo........oooOO0OOooo........oo
319