| Geant4 Cross Reference |
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: S. Meylan and C. Villagrasa (IRSN, France)
27 // Updated: H. Tran (IRSN), France: 20/12/2018
28 //
29
30 #include "DNAParser.hh"
31
32 #include "G4Box.hh"
33 #include "G4DNAChemistryManager.hh"
34 #include "G4LogicalVolume.hh"
35 #include "G4Material.hh"
36 #include "G4NistManager.hh"
37 #include "G4Orb.hh"
38 #include "G4PVPlacement.hh"
39 #include "G4SubtractionSolid.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4VPhysicalVolume.hh"
42 #include "G4VSolid.hh"
43
44 #include <fstream>
45
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47 // Molecule struct def
48 struct DNAParser::Molecule
49 {
50 Molecule(std::string name, int copyNumber, const G4ThreeVector& position, double radius,
51 double waterRadius, const std::string& material, int strand)
52 : fName(name),
53 fMaterial(material),
54 fCopyNumber(copyNumber),
55 fPosition(position),
56 fRadius(radius),
57 fRadiusWater(waterRadius),
58 fStrand(strand)
59 {}
60
61 G4String fName;
62 G4String fMaterial;
63 G4int fCopyNumber;
64 G4ThreeVector fPosition;
65 G4double fRadius;
66 G4double fRadiusWater;
67 G4int fStrand;
68
69 G4bool operator<(const Molecule& rhs) const { return (fPosition.z() < rhs.fPosition.z()); }
70 };
71
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73 DNAParser::DNAParser() : fSize(0), fGeoName(""), fpWater(nullptr), fpGun(new G4MoleculeGun())
74 {
75 EnumParser();
76 }
77
78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
79
80 DNAParser::~DNAParser() = default;
81
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
83
84 G4LogicalVolume* DNAParser::CreateLogicVolume()
85 {
86 G4NistManager* pMan = G4NistManager::Instance();
87 fpWater = pMan->FindOrBuildMaterial("G4_WATER");
88
89 G4String boxNameSolid = fGeoName + "_solid";
90
91 G4Box* pBoxSolid = new G4Box(boxNameSolid, fSize / 2, fSize / 2, fSize / 2);
92
93 G4String boxNameLogic = fGeoName + "_logic";
94
95 G4LogicalVolume* pBoxLogic = new G4LogicalVolume(pBoxSolid, fpWater, boxNameLogic);
96
97 std::sort(fMolecules.begin(), fMolecules.end());
98
99 for (int i = 0, ie = fMolecules.size(); i < ie; ++i) {
100 G4String name = fMolecules[i].fName;
101 G4double radius = fMolecules[i].fRadius;
102 G4double waterRadius = fMolecules[i].fRadiusWater;
103 G4ThreeVector moleculePosition = fMolecules[i].fPosition;
104
105 int copyNum = fMolecules[i].fCopyNumber;
106
107 G4Orb* pMoleculeWaterSolid = nullptr;
108 G4VSolid* pMoleculeWaterCutSolid = nullptr;
109 G4LogicalVolume* pMoleculeWaterLogic = nullptr;
110 G4VPhysicalVolume* pPhysicalVolumeWater = nullptr;
111
112 G4double tolerence = 1e-4 * nm;
113
114 if (waterRadius > tolerence) {
115 G4String nameWaterSolid = name + "_water_solid";
116 pMoleculeWaterSolid = new G4Orb(nameWaterSolid, waterRadius);
117 pMoleculeWaterCutSolid =
118 CreateCutSolid(pMoleculeWaterSolid, fMolecules[i], fMolecules, false);
119
120 G4String nameWaterLogic = name + "_water_logic";
121
122 pMoleculeWaterLogic = new G4LogicalVolume(pMoleculeWaterCutSolid, fpWater, nameWaterLogic);
123 G4String nameWaterPhys = name + "_water";
124 pPhysicalVolumeWater = new G4PVPlacement(0, moleculePosition, pMoleculeWaterLogic,
125 nameWaterPhys, pBoxLogic, false, copyNum);
126
127 auto it = fEnumMap.find(nameWaterPhys);
128 if (it != fEnumMap.end()) {
129 fGeometryMap.insert(std::make_pair(pPhysicalVolumeWater, it->second));
130 }
131 else {
132 G4ExceptionDescription exceptionDescription;
133 exceptionDescription << nameWaterPhys << " could not be exsited";
134 G4Exception("DNAParser::CreateLogicVolume()", "DNAParser001", FatalException,
135 exceptionDescription);
136 }
137 }
138 // Dna volume part
139
140 G4Orb* pMoleculeSolid = nullptr;
141 G4VSolid* pMoleculeCutSolid = nullptr;
142 G4LogicalVolume* pMoleculeLogic = nullptr;
143 G4VPhysicalVolume* pPhysicalVolumeMolecule = nullptr;
144
145 G4String nameSolid = fMolecules[i].fName + "_solid";
146 pMoleculeSolid = new G4Orb(nameSolid, radius);
147
148 pMoleculeCutSolid = CreateCutSolid(pMoleculeSolid, fMolecules[i], fMolecules, true);
149
150 G4String nameLogic = name + "_logic";
151
152 pMoleculeLogic = new G4LogicalVolume(pMoleculeCutSolid, fpWater, nameLogic);
153 if (waterRadius > tolerence) {
154 G4ThreeVector position(0, 0, 0);
155 std::string namePhys = name;
156 pPhysicalVolumeMolecule = new G4PVPlacement(0, position, pMoleculeLogic, namePhys,
157 pMoleculeWaterLogic, false, copyNum);
158
159 auto it = fEnumMap.find(namePhys);
160 if (it != fEnumMap.end()) {
161 fGeometryMap.insert(std::make_pair(pPhysicalVolumeMolecule, it->second));
162 }
163 else {
164 G4ExceptionDescription exceptionDescription;
165 exceptionDescription << namePhys << " could not be exsited";
166 G4Exception("DNAParser::CreateLogicVolume()", "DNAParser002", FatalException,
167 exceptionDescription);
168 }
169 }
170 else {
171 G4ThreeVector position = moleculePosition;
172 G4String namePhys = name;
173 pPhysicalVolumeMolecule =
174 new G4PVPlacement(0, position, pMoleculeLogic, namePhys, pBoxLogic, false, copyNum);
175
176 auto it = fEnumMap.find(namePhys);
177
178 if (it != fEnumMap.end()) {
179 fGeometryMap.insert(std::make_pair(pPhysicalVolumeMolecule, it->second));
180 }
181 else {
182 G4ExceptionDescription exceptionDescription;
183 exceptionDescription << namePhys << " could not be exsited";
184 G4Exception("DNAParser::CreateLogicVolume()", "DNAParser003", FatalException,
185 exceptionDescription);
186 }
187 }
188 }
189 fMolecules.clear();
190 fRadiusMap.clear();
191 fWaterRadiusMap.clear();
192 return pBoxLogic;
193 }
194
195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
196
197 void DNAParser::ParseFile(const std::string& fileName)
198 {
199 fMolecules.clear();
200 fRadiusMap.clear();
201 fWaterRadiusMap.clear();
202
203 std::ifstream file(fileName.c_str());
204
205 if (!file.is_open()) {
206 G4ExceptionDescription exceptionDescription;
207 exceptionDescription << fileName << "could not be opened";
208 G4Exception("DNAParser::ParseFile()", "DNAParser002", FatalException, exceptionDescription);
209 }
210 std::string line;
211 while (std::getline(file, line)) {
212 if (line.empty()) continue;
213 std::istringstream issLine(line);
214 std::string firstItem;
215 issLine >> firstItem;
216
217 if (firstItem == "_Name") {
218 G4String name = "";
219 issLine >> name;
220 fGeoName = name;
221 }
222 if (firstItem == "_Size") {
223 G4double size;
224 issLine >> size;
225 size *= nm;
226 fSize = size;
227 }
228 if (firstItem == "_Radius") {
229 G4String name;
230 issLine >> name;
231 G4double radius;
232 issLine >> radius;
233 radius *= nm;
234 G4double waterRadius;
235 issLine >> waterRadius;
236 waterRadius *= nm;
237 fRadiusMap[name] = radius;
238 fWaterRadiusMap[name] = waterRadius;
239 }
240 if (firstItem == "_pl") {
241 std::string name;
242 issLine >> name;
243 G4String material;
244 issLine >> material;
245
246 G4int strand;
247 issLine >> strand;
248
249 G4int copyNumber;
250 issLine >> copyNumber;
251
252 G4double x;
253 issLine >> x;
254 x *= nm;
255
256 G4double y;
257 issLine >> y;
258 y *= nm;
259
260 G4double z;
261 issLine >> z;
262 z *= nm;
263
264 Molecule molecule(name, copyNumber, G4ThreeVector(x, y, z), fRadiusMap[name],
265 fWaterRadiusMap[name], material, strand);
266 fMolecules.push_back(molecule);
267
268 auto itt = fEnumMap.find(name);
269
270 if (itt != fEnumMap.end()) {
271 if (itt->second != DNAVolumeType::phosphate1 && itt->second != DNAVolumeType::phosphate2) {
272 if (itt->second == DNAVolumeType::deoxyribose1
273 || itt->second == DNAVolumeType::deoxyribose2)
274 {
275 name = "Deoxyribose";
276 }
277 if (itt->second == DNAVolumeType::base_adenine) {
278 name = "Adenine";
279 }
280 if (itt->second == DNAVolumeType::base_thymine) {
281 name = "Thymine";
282 }
283 if (itt->second == DNAVolumeType::base_cytosine) {
284 name = "Cytosine";
285 }
286 if (itt->second == DNAVolumeType::base_guanine) {
287 name = "Guanine";
288 }
289 if (itt->second == DNAVolumeType::histone) {
290 name = "Histone";
291 }
292
293 fpGun->AddMolecule(name, G4ThreeVector(x, y, z), 1 * picosecond);
294 }
295 }
296 else {
297 G4ExceptionDescription exceptionDescription;
298 exceptionDescription << name << " could not be exsited";
299 G4Exception("DNAParser::ParseFile()", "DNAParser005", FatalException, exceptionDescription);
300 }
301 }
302 }
303 file.close();
304 }
305
306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
307
308 G4VSolid* DNAParser::CreateCutSolid(G4Orb* pSolidOrbRef, Molecule& molRef,
309 std::vector<Molecule>& molList, G4bool in)
310 {
311 // The idea behing this method is to cut overlap volumes by selecting
312 // one of them (the reference) and checking all the other volumes (the targets).
313 // If a reference and a target volumes are close enough to overlap they will be cut.
314 // The reference is already selected when we enter this method.
315 // Use the tiny space to differentiate the frontiers (may not be necessary)
316 G4double tinySpace = 0.001 * nm;
317
318 G4SubtractionSolid* pSolidCut(NULL);
319 G4bool isCutted = false;
320 G4bool isOurVol = false;
321 G4double radiusRef;
322
323 if (molRef.fRadiusWater == 0) {
324 radiusRef = molRef.fRadius;
325 }
326 else {
327 radiusRef = molRef.fRadiusWater;
328 }
329
330 G4ThreeVector posRef = molRef.fPosition;
331
332 if (std::abs(posRef.x()) + radiusRef > fSize / 2 // along x
333 || std::abs(posRef.y()) + radiusRef > fSize / 2 // along y
334 || std::abs(posRef.z()) + radiusRef > fSize / 2) // along z
335 {
336 G4Box* solidBox = new G4Box("solid_box_for_cut", fSize / 2, fSize / 2, fSize / 2);
337 G4ThreeVector posBox;
338 G4RotationMatrix rotMat;
339 rotMat.rotate(0, G4ThreeVector(0, 0, 1));
340
341 if (std::abs(posRef.y() + radiusRef) > fSize / 2) {
342 posBox = -posRef + G4ThreeVector(0, fSize, 0);
343 if (!isCutted) {
344 pSolidCut = new G4SubtractionSolid("solidCut", pSolidOrbRef, solidBox, &rotMat, posBox);
345 isCutted = true;
346 }
347 else {
348 pSolidCut = new G4SubtractionSolid("solidCut", pSolidCut, solidBox, &rotMat, posBox);
349 }
350 }
351 // Down
352 if (std::abs(posRef.y() - radiusRef) > fSize / 2) {
353 posBox = -posRef + G4ThreeVector(0, -fSize, 0);
354
355 if (!isCutted) {
356 pSolidCut = new G4SubtractionSolid("solidCut", pSolidOrbRef, solidBox, &rotMat, posBox);
357 isCutted = true;
358 }
359 else {
360 pSolidCut = new G4SubtractionSolid("solidCut", pSolidCut, solidBox, &rotMat, posBox);
361 }
362 }
363 // Left
364 if (std::abs(posRef.x() + radiusRef) > fSize / 2) {
365 posBox = -posRef + G4ThreeVector(fSize, 0, 0);
366 if (!isCutted) {
367 pSolidCut = new G4SubtractionSolid("solidCut", pSolidOrbRef, solidBox, &rotMat, posBox);
368 isCutted = true;
369 }
370 else {
371 pSolidCut = new G4SubtractionSolid("solidCut", pSolidCut, solidBox, &rotMat, posBox);
372 }
373 }
374
375 // Right
376 if (std::abs(posRef.x() - radiusRef) > fSize / 2) {
377 posBox = -posRef + G4ThreeVector(-fSize, 0, 0);
378 if (!isCutted) {
379 pSolidCut = new G4SubtractionSolid("solidCut", pSolidOrbRef, solidBox, &rotMat, posBox);
380 isCutted = true;
381 }
382 else {
383 pSolidCut = new G4SubtractionSolid("solidCut", pSolidCut, solidBox, &rotMat, posBox);
384 }
385 }
386 // Forward
387 if (std::abs(posRef.z() + radiusRef) > fSize / 2) {
388 posBox = -posRef + G4ThreeVector(0, 0, fSize);
389 if (!isCutted) {
390 pSolidCut = new G4SubtractionSolid("solidCut", pSolidOrbRef, solidBox, &rotMat, posBox);
391 isCutted = true;
392 }
393 else {
394 pSolidCut = new G4SubtractionSolid("solidCut", pSolidCut, solidBox, &rotMat, posBox);
395 }
396 }
397
398 // Backward
399 if (std::abs(posRef.z() - radiusRef) > fSize / 2) {
400 posBox = -posRef + G4ThreeVector(0, 0, -fSize);
401 if (!isCutted) {
402 pSolidCut = new G4SubtractionSolid("solidCut", pSolidOrbRef, solidBox, &rotMat, posBox);
403 isCutted = true;
404 }
405 else {
406 pSolidCut = new G4SubtractionSolid("solidCut", pSolidCut, solidBox, &rotMat, posBox);
407 }
408 }
409 }
410
411 for (int i = 0, ie = molList.size(); i < ie; ++i) {
412 G4ThreeVector posTar = molList[i].fPosition;
413 G4double rTar = posRef.z();
414 G4double zTar = posTar.z();
415
416 if (zTar > rTar + 20 * nm) {
417 break;
418 }
419 else if (zTar < rTar - 20 * nm) {
420 continue;
421 }
422
423 G4double radiusTar;
424 if (molList[i].fRadiusWater == 0) {
425 radiusTar = molList[i].fRadius;
426 }
427 else {
428 radiusTar = molList[i].fRadiusWater;
429 }
430
431 G4double distance = std::abs((posRef - posTar).getR());
432
433 if (distance == 0 && !isOurVol) {
434 isOurVol = true;
435 continue;
436 }
437 else if (distance == 0 && isOurVol) {
438 G4ExceptionDescription exceptionDescription;
439 exceptionDescription << "CreateCutSolid: Two volumes "
440 << "are placed at the same position.";
441 G4Exception("DNAParser::CreateCutSolid()", "DNAParser004", FatalException,
442 exceptionDescription);
443 }
444 else if (distance <= radiusRef + radiusTar) {
445 G4Box* solidBox = new G4Box("solid_box_for_cut", 2 * radiusTar, 2 * radiusTar, 2 * radiusTar);
446 G4ThreeVector diff = posTar - posRef;
447 G4double d =
448 (std::pow(radiusRef, 2) - std::pow(radiusTar, 2) + std::pow(distance, 2)) / (2 * distance)
449 + solidBox->GetZHalfLength() - tinySpace;
450 if (in) {
451 d -= 2 * tinySpace;
452 }
453 G4ThreeVector pos = d * (diff / diff.getR());
454 G4double phi = std::acos(pos.getZ() / pos.getR());
455 G4double theta = std::acos(pos.getX() / (pos.getR() * std::cos(CLHEP::pi / 2. - phi)));
456
457 if (pos.getY() < 0) {
458 theta = -theta;
459 }
460
461 G4ThreeVector rotAxisForPhi(1 * nm, 0., 0.);
462 rotAxisForPhi.rotateZ(theta + CLHEP::pi / 2.);
463
464 G4RotationMatrix* rotMat = new G4RotationMatrix;
465 rotMat->rotate(-phi, rotAxisForPhi);
466
467 G4ThreeVector rotZAxis(0., 0., 1 * nm);
468 rotMat->rotate(theta, rotZAxis);
469
470 if (!isCutted) {
471 pSolidCut = new G4SubtractionSolid("solidCut", pSolidOrbRef, solidBox, rotMat, pos);
472 }
473 else {
474 pSolidCut = new G4SubtractionSolid("solidCut", pSolidCut, solidBox, rotMat, pos);
475 }
476 isCutted = true;
477 }
478 }
479
480 if (isCutted) {
481 return pSolidCut;
482 }
483 else {
484 return pSolidOrbRef;
485 }
486 }
487
488 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
489
490 void DNAParser::EnumParser()
491 {
492 fEnumMap["deoxyribose1"] = deoxyribose1;
493 fEnumMap["phosphate1"] = phosphate1;
494 fEnumMap["deoxyribose2"] = deoxyribose2;
495 fEnumMap["phosphate2"] = phosphate2;
496
497 fEnumMap["base_cytosine"] = base_cytosine;
498 fEnumMap["base_guanine"] = base_guanine;
499 fEnumMap["base_thymine"] = base_thymine;
500 fEnumMap["base_adenine"] = base_adenine;
501
502 fEnumMap["deoxyribose1_water"] = deoxyribose1_water;
503 fEnumMap["phosphate1_water"] = phosphate1_water;
504 fEnumMap["deoxyribose2_water"] = deoxyribose2_water;
505 fEnumMap["phosphate2_water"] = phosphate2_water;
506
507 fEnumMap["base_adenine_water"] = base_adenine_water;
508 fEnumMap["base_guanine_water"] = base_guanine_water;
509 fEnumMap["base_cytosine_water"] = base_cytosine_water;
510 fEnumMap["base_thymine_water"] = base_thymine_water;
511
512 fEnumMap["histone"] = histone;
513 fEnumMap["physWorld"] = physWorld;
514 fEnumMap["VoxelStraight"] = VoxelStraight;
515 }
516