File indexing completed on 2026-03-31 07:50:27
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033 #include "DNAParser.hh"
0034
0035 #include "G4Box.hh"
0036 #include "G4DNAChemistryManager.hh"
0037 #include "G4LogicalVolume.hh"
0038 #include "G4Material.hh"
0039 #include "G4NistManager.hh"
0040 #include "G4Orb.hh"
0041 #include "G4PVPlacement.hh"
0042 #include "G4SubtractionSolid.hh"
0043 #include "G4SystemOfUnits.hh"
0044 #include "G4VPhysicalVolume.hh"
0045 #include "G4VSolid.hh"
0046
0047 #include <fstream>
0048
0049
0050
0051 struct DNAParser::Molecule
0052 {
0053 Molecule(std::string name, int copyNumber, const G4ThreeVector& position, double radius,
0054 double waterRadius, const std::string& material, int strand)
0055 : fName(name),
0056 fMaterial(material),
0057 fCopyNumber(copyNumber),
0058 fPosition(position),
0059 fRadius(radius),
0060 fRadiusWater(waterRadius),
0061 fStrand(strand)
0062 {}
0063
0064 G4String fName;
0065 G4String fMaterial;
0066 G4int fCopyNumber;
0067 G4ThreeVector fPosition;
0068 G4double fRadius;
0069 G4double fRadiusWater;
0070 G4int fStrand;
0071
0072 G4bool operator<(const Molecule& rhs) const { return (fPosition.z() < rhs.fPosition.z()); }
0073 };
0074
0075
0076 DNAParser::DNAParser() : fSize(0), fGeoName(""), fpWater(nullptr), fpGun(new G4MoleculeGun())
0077 {
0078 EnumParser();
0079 }
0080
0081
0082
0083 DNAParser::~DNAParser() = default;
0084
0085
0086
0087 G4LogicalVolume* DNAParser::CreateLogicVolume()
0088 {
0089 G4NistManager* pMan = G4NistManager::Instance();
0090 fpWater = pMan->FindOrBuildMaterial("G4_WATER");
0091
0092 G4String boxNameSolid = fGeoName + "_solid";
0093
0094 G4Box* pBoxSolid = new G4Box(boxNameSolid, fSize / 2, fSize / 2, fSize / 2);
0095
0096 G4String boxNameLogic = fGeoName + "_logic";
0097
0098 G4LogicalVolume* pBoxLogic = new G4LogicalVolume(pBoxSolid, fpWater, boxNameLogic);
0099
0100 std::sort(fMolecules.begin(), fMolecules.end());
0101
0102 for (int i = 0, ie = fMolecules.size(); i < ie; ++i) {
0103 G4String name = fMolecules[i].fName;
0104 G4double radius = fMolecules[i].fRadius;
0105 G4double waterRadius = fMolecules[i].fRadiusWater;
0106 G4ThreeVector moleculePosition = fMolecules[i].fPosition;
0107
0108 int copyNum = fMolecules[i].fCopyNumber;
0109
0110 G4Orb* pMoleculeWaterSolid = nullptr;
0111 G4VSolid* pMoleculeWaterCutSolid = nullptr;
0112 G4LogicalVolume* pMoleculeWaterLogic = nullptr;
0113 G4VPhysicalVolume* pPhysicalVolumeWater = nullptr;
0114
0115 G4double tolerence = 1e-4 * nm;
0116
0117 if (waterRadius > tolerence) {
0118 G4String nameWaterSolid = name + "_water_solid";
0119 pMoleculeWaterSolid = new G4Orb(nameWaterSolid, waterRadius);
0120 pMoleculeWaterCutSolid =
0121 CreateCutSolid(pMoleculeWaterSolid, fMolecules[i], fMolecules, false);
0122
0123 G4String nameWaterLogic = name + "_water_logic";
0124
0125 pMoleculeWaterLogic = new G4LogicalVolume(pMoleculeWaterCutSolid, fpWater, nameWaterLogic);
0126 G4String nameWaterPhys = name + "_water";
0127 pPhysicalVolumeWater = new G4PVPlacement(0, moleculePosition, pMoleculeWaterLogic,
0128 nameWaterPhys, pBoxLogic, false, copyNum);
0129
0130 auto it = fEnumMap.find(nameWaterPhys);
0131 if (it != fEnumMap.end()) {
0132 fGeometryMap.insert(std::make_pair(pPhysicalVolumeWater, it->second));
0133 }
0134 else {
0135 G4ExceptionDescription exceptionDescription;
0136 exceptionDescription << nameWaterPhys << " could not be exsited";
0137 G4Exception("DNAParser::CreateLogicVolume()", "DNAParser001", FatalException,
0138 exceptionDescription);
0139 }
0140 }
0141
0142
0143 G4Orb* pMoleculeSolid = nullptr;
0144 G4VSolid* pMoleculeCutSolid = nullptr;
0145 G4LogicalVolume* pMoleculeLogic = nullptr;
0146 G4VPhysicalVolume* pPhysicalVolumeMolecule = nullptr;
0147
0148 G4String nameSolid = fMolecules[i].fName + "_solid";
0149 pMoleculeSolid = new G4Orb(nameSolid, radius);
0150
0151 pMoleculeCutSolid = CreateCutSolid(pMoleculeSolid, fMolecules[i], fMolecules, true);
0152
0153 G4String nameLogic = name + "_logic";
0154
0155 pMoleculeLogic = new G4LogicalVolume(pMoleculeCutSolid, fpWater, nameLogic);
0156 if (waterRadius > tolerence) {
0157 G4ThreeVector position(0, 0, 0);
0158 std::string namePhys = name;
0159 pPhysicalVolumeMolecule = new G4PVPlacement(0, position, pMoleculeLogic, namePhys,
0160 pMoleculeWaterLogic, false, copyNum);
0161
0162 auto it = fEnumMap.find(namePhys);
0163 if (it != fEnumMap.end()) {
0164 fGeometryMap.insert(std::make_pair(pPhysicalVolumeMolecule, it->second));
0165 }
0166 else {
0167 G4ExceptionDescription exceptionDescription;
0168 exceptionDescription << namePhys << " could not be exsited";
0169 G4Exception("DNAParser::CreateLogicVolume()", "DNAParser002", FatalException,
0170 exceptionDescription);
0171 }
0172 }
0173 else {
0174 G4ThreeVector position = moleculePosition;
0175 G4String namePhys = name;
0176 pPhysicalVolumeMolecule =
0177 new G4PVPlacement(0, position, pMoleculeLogic, namePhys, pBoxLogic, false, copyNum);
0178
0179 auto it = fEnumMap.find(namePhys);
0180
0181 if (it != fEnumMap.end()) {
0182 fGeometryMap.insert(std::make_pair(pPhysicalVolumeMolecule, it->second));
0183 }
0184 else {
0185 G4ExceptionDescription exceptionDescription;
0186 exceptionDescription << namePhys << " could not be exsited";
0187 G4Exception("DNAParser::CreateLogicVolume()", "DNAParser003", FatalException,
0188 exceptionDescription);
0189 }
0190 }
0191 }
0192 fMolecules.clear();
0193 fRadiusMap.clear();
0194 fWaterRadiusMap.clear();
0195 return pBoxLogic;
0196 }
0197
0198
0199
0200 void DNAParser::ParseFile(const std::string& fileName)
0201 {
0202 fMolecules.clear();
0203 fRadiusMap.clear();
0204 fWaterRadiusMap.clear();
0205
0206 std::ifstream file(fileName.c_str());
0207
0208 if (!file.is_open()) {
0209 G4ExceptionDescription exceptionDescription;
0210 exceptionDescription << fileName << "could not be opened";
0211 G4Exception("DNAParser::ParseFile()", "DNAParser002", FatalException, exceptionDescription);
0212 }
0213 std::string line;
0214 while (std::getline(file, line)) {
0215 if (line.empty()) continue;
0216 std::istringstream issLine(line);
0217 std::string firstItem;
0218 issLine >> firstItem;
0219
0220 if (firstItem == "_Name") {
0221 G4String name = "";
0222 issLine >> name;
0223 fGeoName = name;
0224 }
0225 if (firstItem == "_Size") {
0226 G4double size;
0227 issLine >> size;
0228 size *= nm;
0229 fSize = size;
0230 }
0231 if (firstItem == "_Radius") {
0232 G4String name;
0233 issLine >> name;
0234 G4double radius;
0235 issLine >> radius;
0236 radius *= nm;
0237 G4double waterRadius;
0238 issLine >> waterRadius;
0239 waterRadius *= nm;
0240 fRadiusMap[name] = radius;
0241 fWaterRadiusMap[name] = waterRadius;
0242 }
0243 if (firstItem == "_pl") {
0244 std::string name;
0245 issLine >> name;
0246 G4String material;
0247 issLine >> material;
0248
0249 G4int strand;
0250 issLine >> strand;
0251
0252 G4int copyNumber;
0253 issLine >> copyNumber;
0254
0255 G4double x;
0256 issLine >> x;
0257 x *= nm;
0258
0259 G4double y;
0260 issLine >> y;
0261 y *= nm;
0262
0263 G4double z;
0264 issLine >> z;
0265 z *= nm;
0266
0267 Molecule molecule(name, copyNumber, G4ThreeVector(x, y, z), fRadiusMap[name],
0268 fWaterRadiusMap[name], material, strand);
0269 fMolecules.push_back(molecule);
0270
0271 auto itt = fEnumMap.find(name);
0272
0273 if (itt != fEnumMap.end()) {
0274 if (itt->second != DNAVolumeType::phosphate1 && itt->second != DNAVolumeType::phosphate2) {
0275 if (itt->second == DNAVolumeType::deoxyribose1
0276 || itt->second == DNAVolumeType::deoxyribose2)
0277 {
0278 name = "Deoxyribose";
0279 }
0280 if (itt->second == DNAVolumeType::base_adenine) {
0281 name = "Adenine";
0282 }
0283 if (itt->second == DNAVolumeType::base_thymine) {
0284 name = "Thymine";
0285 }
0286 if (itt->second == DNAVolumeType::base_cytosine) {
0287 name = "Cytosine";
0288 }
0289 if (itt->second == DNAVolumeType::base_guanine) {
0290 name = "Guanine";
0291 }
0292 if (itt->second == DNAVolumeType::histone) {
0293 name = "Histone";
0294 }
0295
0296 fpGun->AddMolecule(name, G4ThreeVector(x, y, z), 1 * picosecond);
0297 }
0298 }
0299 else {
0300 G4ExceptionDescription exceptionDescription;
0301 exceptionDescription << name << " could not be exsited";
0302 G4Exception("DNAParser::ParseFile()", "DNAParser005", FatalException, exceptionDescription);
0303 }
0304 }
0305 }
0306 file.close();
0307 }
0308
0309
0310
0311 G4VSolid* DNAParser::CreateCutSolid(G4Orb* pSolidOrbRef, Molecule& molRef,
0312 std::vector<Molecule>& molList, G4bool in)
0313 {
0314
0315
0316
0317
0318
0319 G4double tinySpace = 0.001 * nm;
0320
0321 G4SubtractionSolid* pSolidCut(NULL);
0322 G4bool isCutted = false;
0323 G4bool isOurVol = false;
0324 G4double radiusRef;
0325
0326 if (molRef.fRadiusWater == 0) {
0327 radiusRef = molRef.fRadius;
0328 }
0329 else {
0330 radiusRef = molRef.fRadiusWater;
0331 }
0332
0333 G4ThreeVector posRef = molRef.fPosition;
0334
0335 if (std::abs(posRef.x()) + radiusRef > fSize / 2
0336 || std::abs(posRef.y()) + radiusRef > fSize / 2
0337 || std::abs(posRef.z()) + radiusRef > fSize / 2)
0338 {
0339 G4Box* solidBox = new G4Box("solid_box_for_cut", fSize / 2, fSize / 2, fSize / 2);
0340 G4ThreeVector posBox;
0341 G4RotationMatrix rotMat;
0342 rotMat.rotate(0, G4ThreeVector(0, 0, 1));
0343
0344 if (std::abs(posRef.y() + radiusRef) > fSize / 2) {
0345 posBox = -posRef + G4ThreeVector(0, fSize, 0);
0346 if (!isCutted) {
0347 pSolidCut = new G4SubtractionSolid("solidCut", pSolidOrbRef, solidBox, &rotMat, posBox);
0348 isCutted = true;
0349 }
0350 else {
0351 pSolidCut = new G4SubtractionSolid("solidCut", pSolidCut, solidBox, &rotMat, posBox);
0352 }
0353 }
0354
0355 if (std::abs(posRef.y() - radiusRef) > fSize / 2) {
0356 posBox = -posRef + G4ThreeVector(0, -fSize, 0);
0357
0358 if (!isCutted) {
0359 pSolidCut = new G4SubtractionSolid("solidCut", pSolidOrbRef, solidBox, &rotMat, posBox);
0360 isCutted = true;
0361 }
0362 else {
0363 pSolidCut = new G4SubtractionSolid("solidCut", pSolidCut, solidBox, &rotMat, posBox);
0364 }
0365 }
0366
0367 if (std::abs(posRef.x() + radiusRef) > fSize / 2) {
0368 posBox = -posRef + G4ThreeVector(fSize, 0, 0);
0369 if (!isCutted) {
0370 pSolidCut = new G4SubtractionSolid("solidCut", pSolidOrbRef, solidBox, &rotMat, posBox);
0371 isCutted = true;
0372 }
0373 else {
0374 pSolidCut = new G4SubtractionSolid("solidCut", pSolidCut, solidBox, &rotMat, posBox);
0375 }
0376 }
0377
0378
0379 if (std::abs(posRef.x() - radiusRef) > fSize / 2) {
0380 posBox = -posRef + G4ThreeVector(-fSize, 0, 0);
0381 if (!isCutted) {
0382 pSolidCut = new G4SubtractionSolid("solidCut", pSolidOrbRef, solidBox, &rotMat, posBox);
0383 isCutted = true;
0384 }
0385 else {
0386 pSolidCut = new G4SubtractionSolid("solidCut", pSolidCut, solidBox, &rotMat, posBox);
0387 }
0388 }
0389
0390 if (std::abs(posRef.z() + radiusRef) > fSize / 2) {
0391 posBox = -posRef + G4ThreeVector(0, 0, fSize);
0392 if (!isCutted) {
0393 pSolidCut = new G4SubtractionSolid("solidCut", pSolidOrbRef, solidBox, &rotMat, posBox);
0394 isCutted = true;
0395 }
0396 else {
0397 pSolidCut = new G4SubtractionSolid("solidCut", pSolidCut, solidBox, &rotMat, posBox);
0398 }
0399 }
0400
0401
0402 if (std::abs(posRef.z() - radiusRef) > fSize / 2) {
0403 posBox = -posRef + G4ThreeVector(0, 0, -fSize);
0404 if (!isCutted) {
0405 pSolidCut = new G4SubtractionSolid("solidCut", pSolidOrbRef, solidBox, &rotMat, posBox);
0406 isCutted = true;
0407 }
0408 else {
0409 pSolidCut = new G4SubtractionSolid("solidCut", pSolidCut, solidBox, &rotMat, posBox);
0410 }
0411 }
0412 }
0413
0414 for (int i = 0, ie = molList.size(); i < ie; ++i) {
0415 G4ThreeVector posTar = molList[i].fPosition;
0416 G4double rTar = posRef.z();
0417 G4double zTar = posTar.z();
0418
0419 if (zTar > rTar + 20 * nm) {
0420 break;
0421 }
0422 else if (zTar < rTar - 20 * nm) {
0423 continue;
0424 }
0425
0426 G4double radiusTar;
0427 if (molList[i].fRadiusWater == 0) {
0428 radiusTar = molList[i].fRadius;
0429 }
0430 else {
0431 radiusTar = molList[i].fRadiusWater;
0432 }
0433
0434 G4double distance = std::abs((posRef - posTar).getR());
0435
0436 if (distance == 0 && !isOurVol) {
0437 isOurVol = true;
0438 continue;
0439 }
0440 else if (distance == 0 && isOurVol) {
0441 G4ExceptionDescription exceptionDescription;
0442 exceptionDescription << "CreateCutSolid: Two volumes "
0443 << "are placed at the same position.";
0444 G4Exception("DNAParser::CreateCutSolid()", "DNAParser004", FatalException,
0445 exceptionDescription);
0446 }
0447 else if (distance <= radiusRef + radiusTar) {
0448 G4Box* solidBox = new G4Box("solid_box_for_cut", 2 * radiusTar, 2 * radiusTar, 2 * radiusTar);
0449 G4ThreeVector diff = posTar - posRef;
0450 G4double d =
0451 (std::pow(radiusRef, 2) - std::pow(radiusTar, 2) + std::pow(distance, 2)) / (2 * distance)
0452 + solidBox->GetZHalfLength() - tinySpace;
0453 if (in) {
0454 d -= 2 * tinySpace;
0455 }
0456 G4ThreeVector pos = d * (diff / diff.getR());
0457 G4double phi = std::acos(pos.getZ() / pos.getR());
0458 G4double theta = std::acos(pos.getX() / (pos.getR() * std::cos(CLHEP::pi / 2. - phi)));
0459
0460 if (pos.getY() < 0) {
0461 theta = -theta;
0462 }
0463
0464 G4ThreeVector rotAxisForPhi(1 * nm, 0., 0.);
0465 rotAxisForPhi.rotateZ(theta + CLHEP::pi / 2.);
0466
0467 G4RotationMatrix* rotMat = new G4RotationMatrix;
0468 rotMat->rotate(-phi, rotAxisForPhi);
0469
0470 G4ThreeVector rotZAxis(0., 0., 1 * nm);
0471 rotMat->rotate(theta, rotZAxis);
0472
0473 if (!isCutted) {
0474 pSolidCut = new G4SubtractionSolid("solidCut", pSolidOrbRef, solidBox, rotMat, pos);
0475 }
0476 else {
0477 pSolidCut = new G4SubtractionSolid("solidCut", pSolidCut, solidBox, rotMat, pos);
0478 }
0479 isCutted = true;
0480 }
0481 }
0482
0483 if (isCutted) {
0484 return pSolidCut;
0485 }
0486 else {
0487 return pSolidOrbRef;
0488 }
0489 }
0490
0491
0492
0493 void DNAParser::EnumParser()
0494 {
0495 fEnumMap["deoxyribose1"] = deoxyribose1;
0496 fEnumMap["phosphate1"] = phosphate1;
0497 fEnumMap["deoxyribose2"] = deoxyribose2;
0498 fEnumMap["phosphate2"] = phosphate2;
0499
0500 fEnumMap["base_cytosine"] = base_cytosine;
0501 fEnumMap["base_guanine"] = base_guanine;
0502 fEnumMap["base_thymine"] = base_thymine;
0503 fEnumMap["base_adenine"] = base_adenine;
0504
0505 fEnumMap["deoxyribose1_water"] = deoxyribose1_water;
0506 fEnumMap["phosphate1_water"] = phosphate1_water;
0507 fEnumMap["deoxyribose2_water"] = deoxyribose2_water;
0508 fEnumMap["phosphate2_water"] = phosphate2_water;
0509
0510 fEnumMap["base_adenine_water"] = base_adenine_water;
0511 fEnumMap["base_guanine_water"] = base_guanine_water;
0512 fEnumMap["base_cytosine_water"] = base_cytosine_water;
0513 fEnumMap["base_thymine_water"] = base_thymine_water;
0514
0515 fEnumMap["histone"] = histone;
0516 fEnumMap["physWorld"] = physWorld;
0517 fEnumMap["VoxelStraight"] = VoxelStraight;
0518 }
