Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/medical/dna/dnadamage2/src/PhysGeoImport.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 // Authors: J. Naoki D. Kondo (UCSF, US) : 10/10/2021
 27 //          J. Ramos-Mendez and B. Faddegon (UCSF, US)
 28 //
 29 /// \file PhysGeoImport.cc
 30 /// \brief Implementation of the plasmid load methods for the geometry
 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........oooOO0OOooo........oooOO0OOooo.....
 40 
 41 PhysGeoImport::PhysGeoImport() {}
 42 
 43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
 44 
 45 G4LogicalVolume* PhysGeoImport::CreateLogicVolumeXYZ(G4String fileName)
 46 {
 47   G4NistManager* man = G4NistManager::Instance();
 48   fEnvelopeWater = man->FindOrBuildMaterial("G4_WATER");
 49   fpWater =
 50     man->BuildMaterialWithNewDensity("G4_WATER_MODIFIED", "G4_WATER", 1.0 * g / cm / cm / cm);
 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 * angstrom;
 62   G4double PyRs = 2.9389169420478556 * angstrom;
 63   G4double PzRs = 2.9389169420478556 * angstrom;
 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 + "_solid";
 83   G4Box* box_solid =
 84     new G4Box(boxNameSolid, 0.5 * (fXMax - fXMin) + 0.5 * 3.4 * nm,
 85               0.5 * (fYMax - fYMin) + 0.5 * 3.4 * nm, 0.5 * (fZMax - fZMin) + 0.5 * 3.4 * nm);
 86 
 87   G4String boxNameLogic = fGeoName + "_logic";
 88   G4LogicalVolume* box_logic =
 89     new G4LogicalVolume(box_solid, fEnvelopeWater, boxNameLogic, 0, 0, 0);
 90 
 91   // Desoxyribose
 92 
 93   G4Ellipsoid* RSolidSugar = new G4Ellipsoid("sdeoxyribose", PxRs, PyRs, PzRs, -PzRs, .445 * PzRs);
 94   G4LogicalVolume* RSugar = new G4LogicalVolume(RSolidSugar, fpWater, "ldeoxyribose", 0, 0, 0);
 95   G4VisAttributes* MyVisAtt_Rs = new G4VisAttributes(G4Colour(G4Colour::Red()));
 96   MyVisAtt_Rs->SetForceSolid(true);
 97   RSugar->SetVisAttributes(MyVisAtt_Rs);
 98 
 99   // Phosphoric Acid
100 
101   G4Ellipsoid* YSolidSugar = new G4Ellipsoid("sphosphate", PxYs, PyYs, PzYs, -PzYs, .9 * angstrom);
102 
103   G4LogicalVolume* YSugar = new G4LogicalVolume(YSolidSugar, fpWater, "lphosphate", 0, 0, 0);
104   G4VisAttributes* MyVisAtt_Ys = new G4VisAttributes(G4Colour(G4Colour::Yellow()));
105   MyVisAtt_Ys->SetForceSolid(true);
106   YSugar->SetVisAttributes(MyVisAtt_Ys);
107 
108   // Base Pairs
109 
110   G4Ellipsoid* Base1a = new G4Ellipsoid("BasePair1a", PxBp, PyBp, PzBp, -PzBp, 1.15 * angstrom);
111   G4LogicalVolume* BaseP1a = new G4LogicalVolume(Base1a, fpWater, "BasePair1a", 0, 0, 0);
112   G4VisAttributes* MyVisAtt_Bp1a = new G4VisAttributes(G4Colour(G4Colour::Green()));
113   MyVisAtt_Bp1a->SetForceSolid(true);
114   BaseP1a->SetVisAttributes(MyVisAtt_Bp1a);
115 
116   G4Ellipsoid* Base1b = new G4Ellipsoid("BasePair1b", PxBp, PyBp, PzBp, 1.15 * angstrom, PzBp);
117   G4LogicalVolume* BaseP1b = new G4LogicalVolume(Base1b, fpWater, "BasePair1b", 0, 0, 0);
118 
119   G4VisAttributes* MyVisAtt_Bp1b = new G4VisAttributes(G4Colour(G4Colour::Green()));
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; vertex++) {
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 - xin) * (xfn - xin)) + ((yfn - yin) * (yfn - yin))));
136     G4double theta0 = std::atan2(xfn - xin, yfn - yin);
137     G4double lenght = std::sqrt(((xfn - xin) * (xfn - xin)) + ((yfn - yin) * (yfn - yin))
138                                 + ((zfn - zin) * (zfn - zin)));
139     G4double dl = 1.0 / (lenght / (3.4 * angstrom));
140     G4int nChain = (fVertexes[vertex] - fVertexes[vertex + 1]).mag() / (0.34 * nm);
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, (.5 * PzRs) + PzYs);
155       plus2.rotateX(-des1);
156       plus2.rotateZ(-posi);
157 
158       G4ThreeVector plus2alt = G4ThreeVector(0, 0, (.5 * PzRs) + PzYs);
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(x1, y1, z1);
167       G4ThreeVector position2i = G4ThreeVector(x2, y2, z2) + plus2;
168       G4ThreeVector position2ialt = G4ThreeVector(x2, y2, -z2) - plus2alt;
169       G4ThreeVector position3i = G4ThreeVector(x3, y3, z3);
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) + xin;
177       G4double y = dl * nseg * (yfn - yin) + yin;
178       G4double z = dl * nseg * (zfn - zin) + zin;
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 = G4RotationMatrix();
193       rotm1.rotateX(xrot1);
194       rotm1.rotateZ(-posi);
195       rotm1.rotateY(yrot1);
196       rotm1.rotateX(phi0);
197       rotm1.rotateZ(-theta0);
198       G4ThreeVector position1 = position1i + G4ThreeVector(x, y, z);
199       G4Transform3D transform1(rotm1, position1);
200 
201       G4double yrot1alt = theta + pi;
202       G4double xrot1alt = des1;
203       G4RotationMatrix rotm1alt = G4RotationMatrix();
204       rotm1alt.rotateX(xrot1alt);
205       rotm1alt.rotateZ(-posi);
206       rotm1alt.rotateY(yrot1alt);
207       rotm1alt.rotateX(phi0);
208       rotm1alt.rotateZ(-theta0);
209       G4ThreeVector position1alt = -position1i + G4ThreeVector(x, y, z);
210       G4Transform3D transform1alt(rotm1alt, position1alt);
211 
212       G4double yrot2 = theta;
213       G4double xrot2 = -des2;
214       G4RotationMatrix rotm2 = G4RotationMatrix();
215       rotm2.rotateX(xrot2);
216       rotm2.rotateY(yrot2 - bet1 + 0.8726646259971648);
217       rotm2.rotateX(phi0);
218       rotm2.rotateZ(-theta0);
219       G4ThreeVector position2 = position2i + G4ThreeVector(x, y, z);
220       G4Transform3D transform2(rotm2, position2);
221 
222       G4double yrot2alt = theta + pi;
223       G4double xrot2alt = des2;
224       G4RotationMatrix rotm2alt = G4RotationMatrix();
225       rotm2alt.rotateX(xrot2alt);
226       rotm2alt.rotateY(yrot2alt + bet1 - 0.8726646259971648);
227       rotm2alt.rotateX(phi0);
228       rotm2alt.rotateZ(-theta0);
229       G4ThreeVector position2alt = position2ialt + G4ThreeVector(x, y, z);
230       G4Transform3D transform2alt(rotm2alt, position2alt);
231 
232       G4double yrot3 = theta;
233       G4RotationMatrix rotm3 = G4RotationMatrix();
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 + G4ThreeVector(x, y, z);
240       G4Transform3D transform3(rotm3, position3);
241 
242       G4double yrot3alt = theta + pi;
243       G4RotationMatrix rotm3alt = G4RotationMatrix();
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 + G4ThreeVector(x, y, z);
250       G4Transform3D transform3alt(rotm3alt, position3alt);
251 
252       new G4PVPlacement(transform1, RSugar, "deoxyribose1", box_logic, false, index, false);
253 
254       new G4PVPlacement(transform1alt, RSugar, "deoxyribose2", box_logic, false, index, false);
255 
256       new G4PVPlacement(transform3, BaseP1a, "BasePair1", box_logic, false, index, false);
257 
258       new G4PVPlacement(transform3alt, BaseP1a, "BasePair2", box_logic, false, index, false);
259 
260       new G4PVPlacement(transform2, YSugar, "phosphate1", box_logic, false, index, false);
261 
262       new G4PVPlacement(transform2alt, YSugar, "phosphate2", box_logic, false, index, false);
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, 1});
272       fSampleDNADetails.push_back({-1, index, 2});
273 
274       index++;
275     }
276   }
277 
278   return box_logic;
279 }
280 
281 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
282 
283 void PhysGeoImport::ReadFile(G4String fileName)
284 {
285   G4double x, y, z;
286   fXMin = 1 * mm, fYMin = 1 * mm, fZMin = 1 * mm;
287   fXMax = -1 * mm, fYMax = -1 * mm, fZMax = -1 * mm;
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(fVertexes[i]);
311   }
312 
313   VertRed.push_back(fVertexes[0]);
314 
315   fVertexes = VertRed;
316 }
317 
318 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
319