dnadamage1/src/DNAParser.cc

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