| 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: S. Meylan and C. Villagrasa (IRSN,
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........oo
47 // Molecule struct def
48 struct DNAParser::Molecule
49 {
50 Molecule(std::string name, int copyNumber,
51 double waterRadius, const std::st
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) cons
70 };
71
72 //....oooOO0OOooo........oooOO0OOooo........oo
73 DNAParser::DNAParser() : fSize(0), fGeoName(""
74 {
75 EnumParser();
76 }
77
78 //....oooOO0OOooo........oooOO0OOooo........oo
79
80 DNAParser::~DNAParser() = default;
81
82 //....oooOO0OOooo........oooOO0OOooo........oo
83
84 G4LogicalVolume* DNAParser::CreateLogicVolume(
85 {
86 G4NistManager* pMan = G4NistManager::Instanc
87 fpWater = pMan->FindOrBuildMaterial("G4_WATE
88
89 G4String boxNameSolid = fGeoName + "_solid";
90
91 G4Box* pBoxSolid = new G4Box(boxNameSolid, f
92
93 G4String boxNameLogic = fGeoName + "_logic";
94
95 G4LogicalVolume* pBoxLogic = new G4LogicalVo
96
97 std::sort(fMolecules.begin(), fMolecules.end
98
99 for (int i = 0, ie = fMolecules.size(); i <
100 G4String name = fMolecules[i].fName;
101 G4double radius = fMolecules[i].fRadius;
102 G4double waterRadius = fMolecules[i].fRadi
103 G4ThreeVector moleculePosition = fMolecule
104
105 int copyNum = fMolecules[i].fCopyNumber;
106
107 G4Orb* pMoleculeWaterSolid = nullptr;
108 G4VSolid* pMoleculeWaterCutSolid = nullptr
109 G4LogicalVolume* pMoleculeWaterLogic = nul
110 G4VPhysicalVolume* pPhysicalVolumeWater =
111
112 G4double tolerence = 1e-4 * nm;
113
114 if (waterRadius > tolerence) {
115 G4String nameWaterSolid = name + "_water
116 pMoleculeWaterSolid = new G4Orb(nameWate
117 pMoleculeWaterCutSolid =
118 CreateCutSolid(pMoleculeWaterSolid, fM
119
120 G4String nameWaterLogic = name + "_water
121
122 pMoleculeWaterLogic = new G4LogicalVolum
123 G4String nameWaterPhys = name + "_water"
124 pPhysicalVolumeWater = new G4PVPlacement
125
126
127 auto it = fEnumMap.find(nameWaterPhys);
128 if (it != fEnumMap.end()) {
129 fGeometryMap.insert(std::make_pair(pPh
130 }
131 else {
132 G4ExceptionDescription exceptionDescri
133 exceptionDescription << nameWaterPhys
134 G4Exception("DNAParser::CreateLogicVol
135 exceptionDescription);
136 }
137 }
138 // Dna volume part
139
140 G4Orb* pMoleculeSolid = nullptr;
141 G4VSolid* pMoleculeCutSolid = nullptr;
142 G4LogicalVolume* pMoleculeLogic = nullptr;
143 G4VPhysicalVolume* pPhysicalVolumeMolecule
144
145 G4String nameSolid = fMolecules[i].fName +
146 pMoleculeSolid = new G4Orb(nameSolid, radi
147
148 pMoleculeCutSolid = CreateCutSolid(pMolecu
149
150 G4String nameLogic = name + "_logic";
151
152 pMoleculeLogic = new G4LogicalVolume(pMole
153 if (waterRadius > tolerence) {
154 G4ThreeVector position(0, 0, 0);
155 std::string namePhys = name;
156 pPhysicalVolumeMolecule = new G4PVPlacem
157
158
159 auto it = fEnumMap.find(namePhys);
160 if (it != fEnumMap.end()) {
161 fGeometryMap.insert(std::make_pair(pPh
162 }
163 else {
164 G4ExceptionDescription exceptionDescri
165 exceptionDescription << namePhys << "
166 G4Exception("DNAParser::CreateLogicVol
167 exceptionDescription);
168 }
169 }
170 else {
171 G4ThreeVector position = moleculePositio
172 G4String namePhys = name;
173 pPhysicalVolumeMolecule =
174 new G4PVPlacement(0, position, pMolecu
175
176 auto it = fEnumMap.find(namePhys);
177
178 if (it != fEnumMap.end()) {
179 fGeometryMap.insert(std::make_pair(pPh
180 }
181 else {
182 G4ExceptionDescription exceptionDescri
183 exceptionDescription << namePhys << "
184 G4Exception("DNAParser::CreateLogicVol
185 exceptionDescription);
186 }
187 }
188 }
189 fMolecules.clear();
190 fRadiusMap.clear();
191 fWaterRadiusMap.clear();
192 return pBoxLogic;
193 }
194
195 //....oooOO0OOooo........oooOO0OOooo........oo
196
197 void DNAParser::ParseFile(const std::string& f
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 exceptionDescriptio
207 exceptionDescription << fileName << "could
208 G4Exception("DNAParser::ParseFile()", "DNA
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, G4Th
265 fWaterRadiusMap[name],
266 fMolecules.push_back(molecule);
267
268 auto itt = fEnumMap.find(name);
269
270 if (itt != fEnumMap.end()) {
271 if (itt->second != DNAVolumeType::phos
272 if (itt->second == DNAVolumeType::de
273 || itt->second == DNAVolumeType:
274 {
275 name = "Deoxyribose";
276 }
277 if (itt->second == DNAVolumeType::ba
278 name = "Adenine";
279 }
280 if (itt->second == DNAVolumeType::ba
281 name = "Thymine";
282 }
283 if (itt->second == DNAVolumeType::ba
284 name = "Cytosine";
285 }
286 if (itt->second == DNAVolumeType::ba
287 name = "Guanine";
288 }
289 if (itt->second == DNAVolumeType::hi
290 name = "Histone";
291 }
292
293 fpGun->AddMolecule(name, G4ThreeVect
294 }
295 }
296 else {
297 G4ExceptionDescription exceptionDescri
298 exceptionDescription << name << " coul
299 G4Exception("DNAParser::ParseFile()",
300 }
301 }
302 }
303 file.close();
304 }
305
306 //....oooOO0OOooo........oooOO0OOooo........oo
307
308 G4VSolid* DNAParser::CreateCutSolid(G4Orb* pSo
309 std::vecto
310 {
311 // The idea behing this method is to cut ove
312 // one of them (the reference) and checking
313 // If a reference and a target volumes are c
314 // The reference is already selected when we
315 // Use the tiny space to differentiate the f
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
333 || std::abs(posRef.y()) + radiusRef > fS
334 || std::abs(posRef.z()) + radiusRef > fS
335 {
336 G4Box* solidBox = new G4Box("solid_box_for
337 G4ThreeVector posBox;
338 G4RotationMatrix rotMat;
339 rotMat.rotate(0, G4ThreeVector(0, 0, 1));
340
341 if (std::abs(posRef.y() + radiusRef) > fSi
342 posBox = -posRef + G4ThreeVector(0, fSiz
343 if (!isCutted) {
344 pSolidCut = new G4SubtractionSolid("so
345 isCutted = true;
346 }
347 else {
348 pSolidCut = new G4SubtractionSolid("so
349 }
350 }
351 // Down
352 if (std::abs(posRef.y() - radiusRef) > fSi
353 posBox = -posRef + G4ThreeVector(0, -fSi
354
355 if (!isCutted) {
356 pSolidCut = new G4SubtractionSolid("so
357 isCutted = true;
358 }
359 else {
360 pSolidCut = new G4SubtractionSolid("so
361 }
362 }
363 // Left
364 if (std::abs(posRef.x() + radiusRef) > fSi
365 posBox = -posRef + G4ThreeVector(fSize,
366 if (!isCutted) {
367 pSolidCut = new G4SubtractionSolid("so
368 isCutted = true;
369 }
370 else {
371 pSolidCut = new G4SubtractionSolid("so
372 }
373 }
374
375 // Right
376 if (std::abs(posRef.x() - radiusRef) > fSi
377 posBox = -posRef + G4ThreeVector(-fSize,
378 if (!isCutted) {
379 pSolidCut = new G4SubtractionSolid("so
380 isCutted = true;
381 }
382 else {
383 pSolidCut = new G4SubtractionSolid("so
384 }
385 }
386 // Forward
387 if (std::abs(posRef.z() + radiusRef) > fSi
388 posBox = -posRef + G4ThreeVector(0, 0, f
389 if (!isCutted) {
390 pSolidCut = new G4SubtractionSolid("so
391 isCutted = true;
392 }
393 else {
394 pSolidCut = new G4SubtractionSolid("so
395 }
396 }
397
398 // Backward
399 if (std::abs(posRef.z() - radiusRef) > fSi
400 posBox = -posRef + G4ThreeVector(0, 0, -
401 if (!isCutted) {
402 pSolidCut = new G4SubtractionSolid("so
403 isCutted = true;
404 }
405 else {
406 pSolidCut = new G4SubtractionSolid("so
407 }
408 }
409 }
410
411 for (int i = 0, ie = molList.size(); i < ie;
412 G4ThreeVector posTar = molList[i].fPositio
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 - pos
432
433 if (distance == 0 && !isOurVol) {
434 isOurVol = true;
435 continue;
436 }
437 else if (distance == 0 && isOurVol) {
438 G4ExceptionDescription exceptionDescript
439 exceptionDescription << "CreateCutSolid:
440 << "are placed at t
441 G4Exception("DNAParser::CreateCutSolid()
442 exceptionDescription);
443 }
444 else if (distance <= radiusRef + radiusTar
445 G4Box* solidBox = new G4Box("solid_box_f
446 G4ThreeVector diff = posTar - posRef;
447 G4double d =
448 (std::pow(radiusRef, 2) - std::pow(rad
449 + solidBox->GetZHalfLength() - tinySpa
450 if (in) {
451 d -= 2 * tinySpace;
452 }
453 G4ThreeVector pos = d * (diff / diff.get
454 G4double phi = std::acos(pos.getZ() / po
455 G4double theta = std::acos(pos.getX() /
456
457 if (pos.getY() < 0) {
458 theta = -theta;
459 }
460
461 G4ThreeVector rotAxisForPhi(1 * nm, 0.,
462 rotAxisForPhi.rotateZ(theta + CLHEP::pi
463
464 G4RotationMatrix* rotMat = new G4Rotatio
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("so
472 }
473 else {
474 pSolidCut = new G4SubtractionSolid("so
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........oo
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"] = deoxyribose
503 fEnumMap["phosphate1_water"] = phosphate1_wa
504 fEnumMap["deoxyribose2_water"] = deoxyribose
505 fEnumMap["phosphate2_water"] = phosphate2_wa
506
507 fEnumMap["base_adenine_water"] = base_adenin
508 fEnumMap["base_guanine_water"] = base_guanin
509 fEnumMap["base_cytosine_water"] = base_cytos
510 fEnumMap["base_thymine_water"] = base_thymin
511
512 fEnumMap["histone"] = histone;
513 fEnumMap["physWorld"] = physWorld;
514 fEnumMap["VoxelStraight"] = VoxelStraight;
515 }
516