Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-02-23 09:22:01 

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 //
0027 /// file: DNAGeometry.cc
0028 /// brief:
0029 /*
0030  * This class builds the DNA geometry. It interacts in a non-standard way
0031  * with the DNAWorld class, to create a physical DNA geometry
0032  * in a parallel world.
0033  *
0034  */
0035 #include "DNAGeometry.hh"
0036 
0037 #include "DNAGeometryMessenger.hh"
0038 #include "OctreeNode.hh"
0039 #include "PlacementVolumeInfo.hh"
0040 #include "VirtualChromosome.hh"
0041 
0042 #include "G4Box.hh"
0043 #include "G4Color.hh"
0044 #include "G4Ellipsoid.hh"
0045 #include "G4NistManager.hh"
0046 #include "G4PVParameterised.hh"
0047 #include "G4PVPlacement.hh"
0048 #include "G4PhysicalConstants.hh"
0049 #include "G4VisAttributes.hh"
0050 
0051 #include <cmath>
0052 #include <fstream>
0053 
0054 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0055 
0056 G4bool compareLVByName::operator()(const G4LogicalVolume* lhs, const G4LogicalVolume* rhs) const
0057 {
0058   return lhs->GetName() < rhs->GetName();
0059 }
0060 
0061 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0062 
0063 DNAGeometry::DNAGeometry()
0064 {
0065   fpMessenger = new DNAGeometryMessenger(this);
0066   fpChromosomeMapper = new ChromosomeMapper();
0067 
0068   // Create Damage Model and add some default values
0069   fpDamageModel = new DamageModel();
0070   G4NistManager* man = G4NistManager::Instance();
0071   fpSugarMaterial = man->FindOrBuildMaterial("G4_DNA_DEOXYRIBOSE");
0072   fpMediumMaterial = man->FindOrBuildMaterial("G4_WATER");
0073   fpPhosphateMaterial = man->FindOrBuildMaterial("G4_DNA_PHOSPHATE");
0074   fpGuanineMaterial = man->FindOrBuildMaterial("G4_DNA_GUANINE");
0075   fpAdenineMaterial = man->FindOrBuildMaterial("G4_DNA_ADENINE");
0076   fpThymineMaterial = man->FindOrBuildMaterial("G4_DNA_THYMINE");
0077   fpCytosineMaterial = man->FindOrBuildMaterial("G4_DNA_CYTOSINE");
0078   fpHistoneMaterial = man->FindOrBuildMaterial("G4_WATER");
0079   fpDNAPhysicsWorld = new DNAWorld();
0080   fpDrawCellAttrs = new G4VisAttributes();
0081   fpDrawCellAttrs->SetColor(G4Color(0, 0, 1, 0.2));
0082   fpDrawCellAttrs->SetForceSolid(true);
0083   // set default molecule sizes
0084   fMoleculeSizes["phosphate"] = G4ThreeVector(2.282354, 2.282354, 2.282354) * angstrom;
0085   fMoleculeSizes["sugar"] = G4ThreeVector(2.632140, 2.632140, 2.632140) * angstrom;
0086   fMoleculeSizes["guanine"] = G4ThreeVector(3.631503, 3.799953, 1.887288) * angstrom;
0087   fMoleculeSizes["cytosine"] = G4ThreeVector(3.597341, 3.066331, 1.779361) * angstrom;
0088   fMoleculeSizes["adenine"] = G4ThreeVector(3.430711, 3.743504, 1.931958) * angstrom;
0089   fMoleculeSizes["thymine"] = G4ThreeVector(4.205943, 3.040448, 2.003359) * angstrom;
0090   fMoleculeSizes["histone"] = G4ThreeVector(25, 25, 25) * angstrom;
0091 }
0092 
0093 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0094 
0095 DNAGeometry::~DNAGeometry()
0096 {
0097   delete fpMessenger;
0098 }
0099 
0100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0101 
0102 void DNAGeometry::BuildDNA(G4LogicalVolume* vol)
0103 {
0104   // TODO: Add Assertion tests to make sure files have been loaded
0105   fpDNAVolumeChem = vol;
0106 
0107   // TODO: Make a complete copy of vol and it's physical volume
0108   auto* volAsEllipsoid = dynamic_cast<G4Ellipsoid*>(vol->GetSolid());
0109   auto* volAsBox = dynamic_cast<G4Box*>(vol->GetSolid());
0110 
0111   G4VSolid* cloneSolid;
0112   if (volAsEllipsoid != nullptr) {
0113     cloneSolid = new G4Ellipsoid(*((G4Ellipsoid*)vol->GetSolid()));
0114   }
0115   else if (volAsBox != nullptr) {
0116     cloneSolid = new G4Box(*((G4Box*)vol->GetSolid()));
0117   }
0118   else {
0119     G4ExceptionDescription errmsg;
0120     errmsg << "An invalid physical volume shape has been used to hold the DNA volume" << G4endl;
0121     G4Exception("MolecularDnaGeometry", "ERR_INVALID_DNA_SHAPE", FatalException, errmsg);
0122     return;
0123   }
0124 
0125   fpDNAVolumePhys = new G4LogicalVolume(cloneSolid, vol->GetMaterial(),
0126                                         "DNAPhysLV");  // dousatsu
0127   FillParameterisedSpace();
0128   fpDNAPhysicsWorld->SetDNAVolumePointer(fpDNAVolumePhys);
0129 }
0130 
0131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0132 
0133 void DNAGeometry::FillParameterisedSpace()
0134 {
0135   // Create map for each logical volume with its parametrisation
0136   // LV, placement_index, copy number, position, rotation
0137   std::vector<placement> placementVector;
0138   std::map<G4LogicalVolume*, G4int> currentCopyNumber;
0139   G4int numberOfChains = 0;
0140 
0141   // read in and create the LV's specified in the macro
0142   std::map<G4String, G4LogicalVolume*> namesToLVs;
0143   for (const auto& it : fVoxelNames) {
0144     G4String thisVoxelName = it.first;
0145     G4String thisVoxelFile = it.second;
0146 
0147     auto voxel = LoadVoxelVolume(thisVoxelName, thisVoxelFile);
0148     G4LogicalVolume* lv = voxel.first;
0149     fInfoMap[lv] = voxel.second;
0150     namesToLVs[thisVoxelName] = lv;
0151     currentCopyNumber[lv] = 0;
0152     numberOfChains = this->GetPVInfo(lv)->GetNumberOfChains();
0153     if (!(numberOfChains == 1 || numberOfChains == 4 || numberOfChains == 8)) {
0154       G4cout << "Geometry contains volumes with none of 1/4/8"
0155              << " chains. This can cause errors" << G4endl;
0156     }
0157   }
0158 
0159   // Go through the voxel types to build a placement vector
0160   for (unsigned ii = 0; ii < fVoxelTypes.size(); ii++) {
0161     G4LogicalVolume* lv = namesToLVs[fVoxelTypes[ii]];
0162 
0163     placement thisPlacement = std::make_tuple(lv, fVoxelIndices[ii], currentCopyNumber[lv]++,
0164                                               fVoxelPositions[ii], fVoxelRotations[ii]);
0165     placementVector.push_back(thisPlacement);
0166   }
0167 
0168   // Do each placement.
0169   std::map<G4String, int64_t> histonesPerChromosome;  // dousatsu
0170   std::map<G4String, int64_t> basePairsPerChromosome;  // dousatsu
0171   // std::map<G4String, long long> basePairsPerChromosome;//ORG
0172   for (auto it = placementVector.begin(); it != placementVector.end(); it++) {
0173     G4LogicalVolume* thisLogical = std::get<0>(*it);
0174     if (thisLogical == nullptr) {
0175       G4Exception("DNAGeometry::FillParameterisedSpace", "ERR_LOG_VOL_EXISTS_FALSE", FatalException,
0176                   "At least one of the placement names"
0177                   " specified doesn't exist");
0178     }
0179     G4int thisIndex = std::get<1>(*it);
0180     G4int thisCopyNo = std::get<2>(*it);
0181 
0182     G4ThreeVector thisPosition = std::get<3>(*it);
0183     G4ThreeVector thisRotation = std::get<4>(*it);
0184     std::stringstream ss;
0185     ss << thisLogical->GetName() << "-" << thisIndex;
0186 
0187     auto inv = new G4RotationMatrix();
0188     inv->rotateX(thisRotation.getX());
0189     inv->rotateY(thisRotation.getY());
0190     inv->rotateZ(thisRotation.getZ());
0191     auto pRot = new G4RotationMatrix(inv->inverse());
0192     delete inv;
0193 
0194     G4bool isPhysicallyPlaced = false;
0195     VirtualChromosome* thisChromosome = this->GetChromosomeMapper()->GetChromosome(thisPosition);
0196     if (thisChromosome != nullptr) {
0197       isPhysicallyPlaced = true;
0198       // Place for Physics
0199       new G4PVPlacement(pRot, thisPosition, thisLogical, ss.str(), fpDNAVolumePhys, false,
0200                         thisCopyNo, this->GetOverlaps());
0201       // Place for Chemistry
0202       new G4PVPlacement(pRot, thisPosition, this->GetMatchingChemVolume(thisLogical), ss.str(),
0203                         fpDNAVolumeChem, false, thisCopyNo, this->GetOverlaps());
0204 
0205       // dousatsu==============
0206       if (histonesPerChromosome.find(thisChromosome->GetName()) == histonesPerChromosome.end()) {
0207         histonesPerChromosome[thisChromosome->GetName()] = 0;
0208       }
0209       histonesPerChromosome[thisChromosome->GetName()] +=
0210         this->GetPVInfo(thisLogical)->GetTotalHistones();
0211       // dousatsu==============
0212 
0213       if (basePairsPerChromosome.find(thisChromosome->GetName()) == basePairsPerChromosome.end()) {
0214         basePairsPerChromosome[thisChromosome->GetName()] = 0;
0215       }
0216       basePairsPerChromosome[thisChromosome->GetName()] +=
0217         this->GetPVInfo(thisLogical)->GetTotalBasePairs();
0218     }
0219 
0220     // add the placement to the local register, even if it isn't in a
0221     // chromosome, as this tracks base pair index number.
0222     if (numberOfChains == 4 || numberOfChains == 8) {
0223       AddFourChainPlacement(it, placementVector.begin(), isPhysicallyPlaced);
0224     }
0225     else if (numberOfChains == 1) {
0226       AddSingleChainPlacement(it, placementVector.begin(), isPhysicallyPlaced);
0227     }
0228     else {
0229       G4ExceptionDescription errmsg;
0230       errmsg << "Having none of 1/4/8 chains in a placement is not "
0231              << "supported, the application will crash" << G4endl;
0232       G4Exception("DNAGeometry::FillParameterisedSpace", "Number of chains not supported",
0233                   FatalException, errmsg);
0234     }
0235   }
0236 
0237   // dousatsu------------------------------
0238   G4cout << "Histones placed per chromosome:" << G4endl;
0239   G4cout << "Chromosome               Histones" << G4endl;
0240   G4cout << "-----------------------------------" << G4endl;
0241   for (auto& it : histonesPerChromosome) {
0242     G4cout << it.first.c_str() << " : " << it.second << G4endl;
0243   }
0244   // dousatsu------------------------------
0245 
0246   G4cout << "Base Pairs placed per chromosome:" << G4endl;
0247   G4cout << "Chromosome               Base Pairs" << G4endl;
0248   G4cout << "-----------------------------------" << G4endl;
0249   for (auto& it : basePairsPerChromosome) {
0250     G4cout << it.first.c_str() << " : " << it.second << G4endl;
0251   }
0252   G4cout << "-------------------------------" << G4endl;
0253 
0254   if (this->GetVerbosity() > 0) {
0255     this->PrintOctreeStats();
0256   }
0257 }
0258 
0259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0260 
0261 void DNAGeometry::AddSingleChainPlacement(const std::vector<placement>::iterator it,
0262                                           const std::vector<placement>::iterator,
0263                                           G4bool isPhysicallyPlaced)
0264 {
0265   G4LogicalVolume* thisLogical = std::get<0>(*it);
0266   G4bool twist = fVoxelTwist[thisLogical->GetName()];
0267   AddNewPlacement(thisLogical, {{0, 1, 2, 3, 4, 5, 6, 7}}, twist, isPhysicallyPlaced);
0268 }
0269 
0270 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0271 
0272 void DNAGeometry::AddFourChainPlacement(const std::vector<placement>::iterator it,
0273                                         const std::vector<placement>::iterator begin,
0274                                         G4bool isPhysicallyPlaced)
0275 {
0276   G4LogicalVolume* thisLogical = std::get<0>(*it);
0277   G4int thisIndex = std::get<1>(*it);
0278   // Use info to build the placement transforms.
0279   std::array<int, 8> global_chain{};
0280   G4bool twist = fVoxelTwist[thisLogical->GetName()];
0281   if (it == begin) {
0282     global_chain = {{0, 1, 2, 3, 4, 5, 6, 7}};
0283   }
0284   else {
0285     // work out global chain
0286     placement previous = *(it - 1);
0287     G4ThreeVector oldRotation = std::get<4>(previous);
0288     G4ThreeVector thisRotation = std::get<4>(*it);
0289 
0290     G4RotationMatrix rotnew;
0291     rotnew.rotateX(thisRotation.getX());
0292     rotnew.rotateY(thisRotation.getY());
0293     rotnew.rotateZ(thisRotation.getZ());
0294 
0295     rotnew.rectify();
0296 
0297     G4RotationMatrix rotold;
0298     rotold.rotateX(oldRotation.getX());
0299     rotold.rotateY(oldRotation.getY());
0300     rotold.rotateZ(oldRotation.getZ());
0301 
0302     G4ThreeVector zdiff = rotold.colZ() - rotnew.colZ();
0303     if (zdiff.mag2() > 0.1) {
0304       rotold = rotold.rotate(halfpi, rotold.colY());
0305     }
0306     rotold.rectify();
0307     // rotnew.colz == rotold.colz now.
0308     // There exists now a rotation around rotnew.colZ()
0309     // that transforms rotold to rotnew.
0310     G4RotationMatrix relative = rotnew * rotold.inverse();
0311     relative.rectify();
0312     G4ThreeVector axis = relative.getAxis();
0313     G4double angle = relative.getDelta();
0314     // get the sign of the rotation, and correct for quadrant.
0315     G4int sign = (axis.dot(rotnew.colZ()) >= 0) ? 1 : -1;
0316     if (sign < 0) {
0317       angle = 2 * pi - angle;
0318     }
0319     G4int zrots = ((G4int)(2.05 * angle / pi)) % 4;  // ORG
0320 
0321     global_chain = {{(this->GetGlobalChain(thisIndex - 1, 0) + zrots) % 4,
0322                      (this->GetGlobalChain(thisIndex - 1, 1) + zrots) % 4,
0323                      (this->GetGlobalChain(thisIndex - 1, 2) + zrots) % 4,
0324                      (this->GetGlobalChain(thisIndex - 1, 3) + zrots) % 4,
0325                      4 + (this->GetGlobalChain(thisIndex - 1, 4) + zrots) % 4,
0326                      4 + (this->GetGlobalChain(thisIndex - 1, 5) + zrots) % 4,
0327                      4 + (this->GetGlobalChain(thisIndex - 1, 6) + zrots) % 4,
0328                      4 + (this->GetGlobalChain(thisIndex - 1, 7) + zrots) % 4}};
0329     twist = false;
0330   }
0331   this->AddNewPlacement(thisLogical, global_chain, twist, isPhysicallyPlaced);
0332 }
0333 
0334 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0335 
0336 void DNAGeometry::AddNewPlacement(const G4LogicalVolume* lv, std::array<int, 8> global_chain,
0337                                   G4bool twist, G4bool isPhysicallyPlaced)
0338 {
0339   placement_transform pt;
0340   std::array<int64_t, 8> start{};  // dousatsu
0341   std::array<int64_t, 8> end{};  // dousatsu
0342   // std::array<long long, 8> start;//ORG
0343   // std::array<long long, 8> end;//ORG
0344   if (fPlacementTransformations.empty()) {
0345     start = {{0, 0, 0, 0, 0, 0, 0, 0}};
0346     end = {{this->GetPVInfo(lv)->GetPairsOnChain(global_chain[0]),
0347             this->GetPVInfo(lv)->GetPairsOnChain(global_chain[1]),
0348             this->GetPVInfo(lv)->GetPairsOnChain(global_chain[2]),
0349             this->GetPVInfo(lv)->GetPairsOnChain(global_chain[3])}};
0350   }
0351   else {
0352     placement_transform previous = fPlacementTransformations.back();
0353     std::array<int64_t, 8> previousEnd = std::get<2>(previous);  // dousatsu
0354     // std::array<long long, 8> previousEnd = std::get<2>(previous);//ORG
0355     start = {{previousEnd[0], previousEnd[1], previousEnd[2], previousEnd[3], previousEnd[4],
0356               previousEnd[5], previousEnd[6], previousEnd[7]}};
0357     end = {{previousEnd[0] + this->GetPVInfo(lv)->GetPairsOnChain(global_chain[0]),
0358             previousEnd[1] + this->GetPVInfo(lv)->GetPairsOnChain(global_chain[1]),
0359             previousEnd[2] + this->GetPVInfo(lv)->GetPairsOnChain(global_chain[2]),
0360             previousEnd[3] + this->GetPVInfo(lv)->GetPairsOnChain(global_chain[3]),
0361             previousEnd[4] + this->GetPVInfo(lv)->GetPairsOnChain(global_chain[4]),
0362             previousEnd[5] + this->GetPVInfo(lv)->GetPairsOnChain(global_chain[5]),
0363             previousEnd[6] + this->GetPVInfo(lv)->GetPairsOnChain(global_chain[6]),
0364             previousEnd[7] + this->GetPVInfo(lv)->GetPairsOnChain(global_chain[7])}};
0365   }
0366   pt = std::make_tuple(global_chain, start, end, twist, isPhysicallyPlaced);
0367   fPlacementTransformations.push_back(pt);
0368 }
0369 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0370 
0371 std::pair<G4LogicalVolume*, PlacementVolumeInfo*>
0372 DNAGeometry::LoadVoxelVolume(const G4String& volumeName, const G4String& filename)
0373 {
0374   if (this->GetVerbosity() > 0) {
0375     G4cout << "Loading Voxel: " << volumeName << G4endl;
0376   }
0377 
0378   auto thisInfo = new PlacementVolumeInfo();
0379 
0380   G4double vxdim = 0.5 * fVoxelSideLength.getX();
0381   G4double vydim = 0.5 * fVoxelSideLength.getY();
0382   G4double vzdim = 0.5 * fVoxelSideLength.getZ();
0383 
0384   // Make Volumes, and put them into the maps for phys/chem comparison
0385   auto physicsPhysical = new G4Box(volumeName, vxdim, vydim, vzdim);
0386   auto physicsLogical = new G4LogicalVolume(physicsPhysical, fpMediumMaterial, volumeName);
0387   auto chemPhysical = new G4Box(volumeName, vxdim, vydim, vzdim);
0388   auto chemLogical = new G4LogicalVolume(chemPhysical, fpMediumMaterial, volumeName);
0389 
0390   if (this->GetDrawCellVolumes()) {
0391     chemLogical->SetVisAttributes(fpDrawCellAttrs);
0392   }
0393 
0394   fChemToPhys[chemLogical] = physicsLogical;
0395   fPhysToChem[physicsLogical] = chemLogical;
0396 
0397   physicsLogical->SetSmartless(this->GetSmartless());
0398 
0399   auto thisOctree = new OctreeNode(G4ThreeVector(0, 0, 0), G4ThreeVector(vxdim, vydim, vzdim), 1);
0400   auto thisHistoneOctree =
0401     new OctreeNode(G4ThreeVector(0, 0, 0), G4ThreeVector(vxdim, vydim, vzdim), 1);
0402 
0403   // open and load file
0404   if (!utility::Path_exists(filename)) {
0405     G4ExceptionDescription errmsg;
0406     errmsg << "The file: " << filename << " could not be found." << G4endl;
0407     G4Exception("DNAGeometry::LoadVoxelVolume", "File Not Found", FatalException, errmsg);
0408   }
0409   std::ifstream infile(filename, std::ifstream::in);
0410   G4String currentline;
0411   G4String lower_name;
0412   G4ThreeVector pos;
0413   G4ThreeVector rot;
0414   G4ThreeVector mol_size;
0415   std::vector<molecule_t> molecules;
0416   std::vector<molecule_t> histones;
0417   G4bool file_contains_size = false;
0418   G4int line_fields = 0;
0419   G4int line_number = 0;
0420   G4int uncommented_line_number = 0;
0421 
0422   if (infile.is_open()) {
0423     while (getline(infile, currentline)) {
0424       if (currentline[0] != '#') {
0425         std::vector<G4String> line = utility::Split(currentline, ' ');
0426         // validation
0427         if (uncommented_line_number == 0) {
0428           line_fields = line.size();
0429           file_contains_size = (line_fields == 14);
0430         }
0431         if ((line_fields != 10) && (line_fields != 14)) {
0432           G4ExceptionDescription errmsg;
0433           errmsg << filename << " should have 10 or 14 fields, only got " << line_fields << G4endl;
0434           G4Exception("DNAGeometry::LoadVoxelVolume", "BadFileFormat", FatalErrorInArgument,
0435                       errmsg);
0436         }
0437         if (line_fields != (G4int)line.size()) {
0438           G4ExceptionDescription errmsg;
0439           errmsg << "Unexpected number of fields on line " << line_number << " of file " << filename
0440                  << G4endl << "Expected " << line_fields << " fields, but got " << line.size()
0441                  << G4endl;
0442           G4Exception("DNAGeometry::LoadVoxelVolume", "BadFileFormat", FatalErrorInArgument,
0443                       errmsg);
0444         }
0445         molecule_t thisMolecule;
0446 
0447         thisMolecule.fname = (G4String)line.at(0);
0448         thisMolecule.fshape = (G4String)line.at(1);
0449         thisMolecule.fchain_id = (G4int)std::stoi(line.at(2));
0450         thisMolecule.fstrand_id = (G4int)std::stoi(line.at(3));
0451         thisMolecule.fbase_idx = (int64_t)std::stoll(line.at(4));
0452 
0453         lower_name = (G4String)line.at(0);
0454         G4StrUtil::to_lower(lower_name);
0455         // lower_name.toLower();
0456 
0457         G4int size_offset = (file_contains_size) ? 3 : 0;
0458         if (file_contains_size) {
0459           mol_size =
0460             G4ThreeVector(std::stod(line.at(5)) * angstrom, std::stod(line.at(6)) * angstrom,
0461                           std::stod(line.at(7)) * angstrom);
0462         }
0463         else {
0464           if (fMoleculeSizes.count(lower_name) == 0) {
0465             G4ExceptionDescription errmsg;
0466             errmsg << "A molecule size has not been defined for " << lower_name << G4endl;
0467             G4Exception("DNAGeometry::LoadVoxelVolume", "MoleculeNotDefined", FatalErrorInArgument,
0468                         errmsg);
0469           }
0470           mol_size = fMoleculeSizes.at(lower_name);
0471         }
0472         pos = G4ThreeVector(std::stod(line.at(5 + size_offset)) * angstrom,
0473                             std::stod(line.at(6 + size_offset)) * angstrom,
0474                             std::stod(line.at(7 + size_offset)) * angstrom);
0475 
0476         rot =
0477           G4ThreeVector(std::stod(line.at(8 + size_offset)), std::stod(line.at(9 + size_offset)),
0478                         std::stod(line.at(10 + size_offset)));
0479 
0480         thisMolecule.fsize = mol_size;
0481         thisMolecule.fposition = pos;
0482         thisMolecule.frotation = rot;
0483 
0484         if (lower_name == "histone") {
0485           histones.push_back(thisMolecule);
0486         }
0487         else {
0488           molecules.push_back(thisMolecule);
0489         }
0490         uncommented_line_number++;
0491       }
0492       line_number++;
0493     }
0494     infile.close();
0495   }
0496   else {
0497     G4ExceptionDescription errmsg;
0498     errmsg << "Voxel Position file read error" << G4endl << "Name: " << volumeName << G4endl
0499            << "Filename: " << filename << G4endl << "Format should be: "
0500            << "# NAME SHAPE CHAIN_ID STRAND_ID BP_INDEX SIZE_X SIZE_Y "
0501            << "SIZE_Z POS_X POS_Y POS_Z ROT_X ROT_Y ROT_Z" << G4endl << "or" << G4endl
0502            << "# NAME SHAPE CHAIN_ID STRAND_ID BP_INDEX POS_X POS_Y POS_Z "
0503               "ROT_X ROT_Y ROT_Z"
0504            << G4endl
0505            << "Separator is space, lines with '#' as the first  character are "
0506               "commented"
0507            << G4endl;
0508 
0509     G4Exception("DNAGeometry::LoadVoxelVolume", "BadFile", FatalErrorInArgument, errmsg);
0510   }
0511 
0512   // We can place the elements now, specify an order that they'll be placed
0513   // for the union solid cuts.
0514   G4VPhysicalVolume* pv_placement = nullptr;
0515 
0516   // Place Histones
0517   G4int his_arr_size = histones.size();
0518   for (G4int ii = 0; ii != his_arr_size; ++ii) {
0519     molecule_t thisMolecule = histones[ii];
0520     if (this->GetVerbosity() > 4) {
0521       G4cout << "Placing molecule " << ii << ": " << thisMolecule.fname << G4endl;
0522     }
0523     // Identify number of histones in voxel volume
0524     pv_placement = PlaceHistone(physicsLogical, thisMolecule);
0525     thisHistoneOctree->AddPhysicalVolume(pv_placement);
0526   }
0527   thisInfo->SetNHistones(his_arr_size);
0528   thisInfo->SetHistoneOctree(thisHistoneOctree);
0529 
0530   // Place Molecules
0531   G4int mol_arr_size = molecules.size();
0532   for (G4int ii = 0; ii != mol_arr_size; ii++) {
0533     molecule_t thisMolecule = molecules[ii];
0534     if (this->GetVerbosity() > 4) {
0535       G4cout << "Placing molecule " << ii << ": " << thisMolecule.fname << G4endl;
0536     }
0537     // Identify number of molecules on chain
0538     if (thisInfo->GetPairsOnChain(thisMolecule.fchain_id) < (thisMolecule.fbase_idx + 1)) {
0539       thisInfo->SetPairsOnChain(thisMolecule.fchain_id, thisMolecule.fbase_idx + 1);
0540     }
0541 
0542     // Place Volume
0543     if (thisMolecule.fname == "Sugar") {
0544       molecule_t phosphate = molecules[ii - 1];
0545       // G4int phosph_idx = (thisMolecule.strand_id == 0) ? ii - 1 : ii - 1;
0546       G4int phosph_idx = ii - 1;
0547       if ((phosph_idx >= 0) && (phosph_idx < (G4int)mol_arr_size)) {
0548         phosphate = molecules[phosph_idx];
0549         if (phosphate.fchain_id != thisMolecule.fchain_id) {
0550           // There is no next sugar, just the voxel wall
0551           phosphate = molecule_t();
0552           phosphate.fname = "wall";
0553         }
0554       }
0555       else {
0556         // There is no next sugar, just the voxel wall
0557         phosphate = molecule_t();
0558         phosphate.fname = "wall";
0559       }
0560       pv_placement = PlaceSugar(physicsLogical, thisMolecule, phosphate);
0561     }
0562     else if (thisMolecule.fname == "Phosphate") {
0563       molecule_t sugar;
0564       G4int sugar_idx = (thisMolecule.fstrand_id == 0) ? ii + 7 : ii - 5;
0565       if ((sugar_idx >= 0) && (sugar_idx < mol_arr_size)) {
0566         sugar = molecules[sugar_idx];
0567         if (sugar.fchain_id != thisMolecule.fchain_id) {
0568           // There is no next sugar, just the voxel wall
0569           sugar = molecule_t();
0570           sugar.fname = "wall";
0571         }
0572       }
0573       else {
0574         // There is no next sugar, just the voxel wall
0575         sugar = molecule_t();
0576         sugar.fname = "wall";
0577       }
0578       pv_placement = PlacePhosphate(physicsLogical, thisMolecule, sugar);
0579     }
0580     else if ((thisMolecule.fname == "Guanine") || (thisMolecule.fname == "Adenine")
0581              || (thisMolecule.fname == "Cytosine") || (thisMolecule.fname == "Thymine"))
0582     {
0583       molecule_t sugar = molecules[ii - 1];
0584       molecule_t oppBasePair;
0585       molecule_t nextBasePair;
0586 
0587       G4int nextBPidx;
0588       G4bool endOfChain = false;
0589       G4bool nextBPexists;
0590       if (thisMolecule.fstrand_id == 0) {
0591         oppBasePair = molecules[ii + 3];
0592         nextBPexists = (ii + 6 < mol_arr_size);
0593         if (nextBPexists) {
0594           endOfChain = (thisMolecule.fchain_id != molecules[ii + 6].fchain_id);
0595         }
0596         nextBPidx = (nextBPexists && !endOfChain) ? ii + 6 : ii - 6;
0597       }
0598       else {
0599         oppBasePair = molecules[ii - 3];
0600         nextBPexists = (ii - 6 >= 0);
0601         if (nextBPexists) {
0602           endOfChain = (thisMolecule.fchain_id != molecules[ii - 6].fchain_id);
0603         }
0604         nextBPidx = (nextBPexists && !endOfChain) ? ii - 6 : ii + 6;
0605       }
0606       nextBasePair = molecules[nextBPidx];
0607       pv_placement = PlaceBase(physicsLogical, thisMolecule, oppBasePair, nextBasePair, sugar);
0608     }
0609     else {
0610       G4ExceptionDescription errmsg;
0611       errmsg << "Molecule name: " << thisMolecule.fname << " is invalid." << G4endl << "Aborting"
0612              << G4endl;
0613       G4Exception("DNAGeometry::LoadVoxelVolume", "Invalid molecule name", FatalException, errmsg);
0614     }
0615     if (this->GetOverlaps()) {
0616       if (pv_placement->CheckOverlaps()) {
0617         // NOTE: if there are overlaps, the simulation will probably
0618         // crash. This is the desired behaviour (Fail Loud)
0619         G4cout << "Overlap, translation: " << pv_placement->GetFrameTranslation() << G4endl;
0620         G4cout << "Overlap, rotationX: " << pv_placement->GetFrameRotation()->colX() << G4endl;
0621         G4cout << "Overlap, rotationY: " << pv_placement->GetFrameRotation()->colY() << G4endl;
0622         G4cout << "Overlap, rotationZ: " << pv_placement->GetFrameRotation()->colZ() << G4endl;
0623       }
0624     }
0625     if (this->GetDrawCellVolumes()) {
0626       // do not draw DNA, as we are drawing the cell volume
0627       pv_placement->GetLogicalVolume()->SetVisAttributes(G4VisAttributes::GetInvisible());
0628     }
0629     thisOctree->AddPhysicalVolume(pv_placement);
0630   }
0631   thisInfo->SetOctree(thisOctree);
0632   return std::make_pair(physicsLogical, thisInfo);
0633 }
0634 
0635 void DNAGeometry::FillVoxelVectors()
0636 {
0637   G4String filename = fFractalCurveFile;
0638   if (!utility::Path_exists(filename)) {
0639     G4ExceptionDescription errmsg;
0640     errmsg << "The file: " << filename << " could not be found." << G4endl;
0641     G4Exception("DNAGeometry::FillVoxelVectors", "File Not Found", FatalException, errmsg);
0642   }
0643   std::ifstream infile(filename, std::ifstream::in);
0644   G4String currentline;
0645   G4ThreeVector pos;
0646   G4ThreeVector rot;
0647 
0648   fVoxelTypes.clear();
0649   fVoxelIndices.clear();
0650   fVoxelPositions.clear();
0651   fVoxelRotations.clear();
0652 
0653   if (infile.is_open()) {
0654     while (getline(infile, currentline)) {
0655       if (currentline[0] != '#') {
0656         try {
0657           G4int pi_or_one = fAnglesAsPi ? pi : 1;
0658           std::vector<G4String> line = utility::Split(currentline, ' ');
0659           fVoxelIndices.push_back(std::stoi(line.at(0)));
0660           fVoxelTypes.push_back((G4String)line.at(1));
0661           pos = G4ThreeVector(std::stod(line.at(2)) * fFractalScaling.getX(),
0662                               std::stod(line.at(3)) * fFractalScaling.getY(),
0663                               std::stod(line.at(4)) * fFractalScaling.getZ());
0664           rot = G4ThreeVector(std::stod(line.at(5)), std::stod(line.at(6)), std::stod(line.at(7)));
0665           fVoxelPositions.push_back(pos);
0666           fVoxelRotations.push_back(rot * pi_or_one);
0667         }
0668         catch (int exception) {
0669           G4cout << "Voxel Position file read error" << G4endl;
0670           G4cout << "Got Line: " << currentline << G4endl;
0671           G4cout << "Format should be: "
0672                  << "IDX KIND POS_X POS_Y POS_Z EUL_PSI EUL_THETA"
0673                  << " EUL_PHI" << G4endl;
0674           G4cout << "Separator is space, lines with '#' as the first "
0675                  << "character are commented" << G4endl;
0676           std::exit(1);
0677         }
0678       }
0679     }
0680     infile.close();
0681   }
0682   else {
0683     G4cout << "Voxel Position file read error on open" << G4endl;
0684     G4cout << "Format should be: "
0685            << "IDX KIND POS_X POS_Y POS_Z EUL_PSI EUL_THETA EUL_PHI" << G4endl;
0686     G4cout << "Separator is space, lines with '#' as the first "
0687            << "character are commented" << G4endl;
0688   }
0689 }
0690 
0691 OctreeNode* DNAGeometry::GetTopOctreeNode(G4LogicalVolume* lv) const
0692 {
0693   const PlacementVolumeInfo* info = GetPVInfo(lv);
0694   return (info == nullptr) ? nullptr : info->GetOctree();
0695 }
0696 
0697 const PlacementVolumeInfo* DNAGeometry::GetPVInfo(const G4LogicalVolume* lv) const
0698 {
0699   if (fInfoMap.find(lv) != fInfoMap.end()) {
0700     return fInfoMap.at(lv);
0701   }
0702   else if (fChemToPhys.find(lv) != fChemToPhys.end()) {
0703     return fInfoMap.at(fChemToPhys.at(lv));
0704   }
0705   else {
0706     return nullptr;
0707   }
0708 }
0709 
0710 void DNAGeometry::PrintOctreeStats()
0711 {
0712   G4cout << "Name         Vols     EndNodes  Max(vols/node)" << G4endl;
0713   G4cout << "----------------------------------------------" << G4endl;
0714   OctreeNode* currentOctree;
0715   for (auto& it : fInfoMap) {
0716     char buffer[80];
0717     currentOctree = it.second->GetOctree();
0718     G4int vols = currentOctree->GetContents().size();
0719     const char* lvname = it.first->GetName().c_str();
0720     G4int mvols = currentOctree->GetMaxContents();
0721     G4int nodes = currentOctree->GetNumberOfTerminalNodes();
0722     std::snprintf(buffer, 80, "%-12s %-7i  %-7i   %-9i", lvname, vols, nodes, mvols);
0723     G4cout << buffer << G4endl;
0724   }
0725   G4cout << "-------------------------------------------" << G4endl;
0726 }
0727 
0728 ////////////////////////////////////////////////////////////////////////////////
0729 ////////////////////////////////////////////////////////////////////////////////
0730 ////////////////////////////////////////////////////////////////////////////////
0731 
0732 G4VPhysicalVolume* DNAGeometry::PlacePhosphate(G4LogicalVolume* physicsLogical,
0733                                                const molecule_t& thisMolecule,
0734                                                const molecule_t& sugar)
0735 {
0736   // Define Vis Attributes
0737   auto vattrs = new G4VisAttributes(G4Color::Yellow());
0738   vattrs->SetForceSolid(true);
0739 
0740   std::stringstream ss;
0741   ss << thisMolecule.fname << "-" << thisMolecule.fchain_id << "-" << thisMolecule.fstrand_id << "-"
0742      << thisMolecule.fbase_idx;
0743   G4String name = ss.str();
0744   // Initial Rotation of the sphere
0745   G4RotationMatrix rot;
0746   rot.rotateX(thisMolecule.frotation.getX());
0747   rot.rotateY(thisMolecule.frotation.getY());
0748   rot.rotateZ(thisMolecule.frotation.getZ());
0749 
0750   G4VPhysicalVolume* pv_placement;
0751   G4RotationMatrix* rot_pointer;
0752   if (sugar.fname == "wall") {
0753     if (this->GetVerbosity() >= 2) {
0754       G4cout << "Wall placement" << G4endl;
0755     }
0756     G4ThreeVector abspos(std::abs(thisMolecule.fposition.getX()),
0757                          std::abs(thisMolecule.fposition.getY()),
0758                          std::abs(thisMolecule.fposition.getZ()));
0759     G4ThreeVector overlaps = abspos + thisMolecule.fsize - fVoxelSideLength;
0760 
0761     G4Ellipsoid* mol;
0762     // NOTE: will work to cut molecules at boundary only when the chain is
0763     // leaving the box.
0764     G4double z_cut = -1 * utility::Min(overlaps);
0765     if (z_cut < 0)  // all molecules fit
0766     {
0767       mol = new G4Ellipsoid(name, thisMolecule.fsize.getX(), thisMolecule.fsize.getY(),
0768                             thisMolecule.fsize.getZ());
0769     }
0770     else  // need to cut molecule to fit in box
0771     {
0772       mol = new G4Ellipsoid(name, thisMolecule.fsize.getX(), thisMolecule.fsize.getY(),
0773                             thisMolecule.fsize.getZ(), -thisMolecule.fsize.getZ(), z_cut);
0774     }
0775 
0776     rot_pointer = new G4RotationMatrix(rot.inverse());
0777     auto molLogic = new G4LogicalVolume(mol, fpPhosphateMaterial, name);
0778 
0779     molLogic->SetVisAttributes(vattrs);
0780     pv_placement = new G4PVPlacement(rot_pointer, thisMolecule.fposition, molLogic, name,
0781                                      physicsLogical, false, 0, this->GetOverlaps());
0782   }
0783   else {
0784     // Case where there is a sugar molecule that we would intersect with
0785     G4ThreeVector z_direction = sugar.fposition - thisMolecule.fposition;
0786     G4double separation = z_direction.mag();
0787     G4double z_cut = separation - sugar.fsize.getX();
0788 
0789     G4ThreeVector local_z = G4ThreeVector(0, 0, 1);
0790     G4ThreeVector z_new = z_direction / z_direction.mag();
0791 
0792     G4ThreeVector ortho = local_z.cross(z_new);
0793     G4double dp = local_z.dot(z_new);
0794     G4double angle = -std::atan2(ortho.mag(), dp);
0795 
0796     G4RotationMatrix second_rot(ortho, angle);
0797     rot_pointer = new G4RotationMatrix(second_rot);
0798 
0799     auto mol = new G4Ellipsoid(name, thisMolecule.fsize.getX(), thisMolecule.fsize.getY(),
0800                                thisMolecule.fsize.getZ(), -thisMolecule.fsize.getZ(), z_cut);
0801 
0802     auto molLogic = new G4LogicalVolume(mol, fpPhosphateMaterial, name);
0803 
0804     molLogic->SetVisAttributes(vattrs);
0805     pv_placement = new G4PVPlacement(rot_pointer, thisMolecule.fposition, molLogic, name,
0806                                      physicsLogical, false, 0, this->GetOverlaps());
0807   }
0808 
0809   return pv_placement;
0810 }
0811 
0812 G4VPhysicalVolume* DNAGeometry::PlaceSugar(G4LogicalVolume* physicsLogical,
0813                                            const molecule_t& thisMolecule,
0814                                            const molecule_t& phosphate)
0815 {
0816   // Define Vis Attributes
0817   auto vattrs = new G4VisAttributes(G4Color::Red());
0818   vattrs->SetForceSolid(true);
0819 
0820   // Name:
0821   std::stringstream ss;
0822   ss << thisMolecule.fname << "-" << thisMolecule.fchain_id << "-" << thisMolecule.fstrand_id << "-"
0823      << thisMolecule.fbase_idx;
0824   G4String name = ss.str();
0825   // Initial Rotation of the sphere
0826   G4RotationMatrix rot;
0827   rot.rotateX(thisMolecule.frotation.getX());
0828   rot.rotateY(thisMolecule.frotation.getY());
0829   rot.rotateZ(thisMolecule.frotation.getZ());
0830 
0831   G4VPhysicalVolume* pv_placement;
0832   G4RotationMatrix* rot_pointer;
0833   if (phosphate.fname == "wall") {
0834     if (this->GetVerbosity() >= 2) {
0835       G4cout << "Wall placement" << G4endl;
0836     }
0837     G4ThreeVector abspos(std::abs(thisMolecule.fposition.getX()),
0838                          std::abs(thisMolecule.fposition.getY()),
0839                          std::abs(thisMolecule.fposition.getZ()));
0840     G4ThreeVector overlaps = abspos + thisMolecule.fsize - fVoxelSideLength;
0841 
0842     G4Ellipsoid* mol;
0843     // NOTE: will work to cut molecules at boundary only when the chain is
0844     // leaving the box.
0845     G4double z_cut = -1 * utility::Min(overlaps);
0846     if (z_cut < 0)  // all molecules fit
0847     {
0848       mol = new G4Ellipsoid(name, thisMolecule.fsize.getX(), thisMolecule.fsize.getY(),
0849                             thisMolecule.fsize.getZ());
0850     }
0851     else  // need to cut molecule to fit in box
0852     {
0853       mol = new G4Ellipsoid(name, thisMolecule.fsize.getX(), thisMolecule.fsize.getY(),
0854                             thisMolecule.fsize.getZ(), -thisMolecule.fsize.getZ(), z_cut);
0855     }
0856 
0857     rot_pointer = new G4RotationMatrix(rot.inverse());
0858     auto molLogic = new G4LogicalVolume(mol, fpSugarMaterial, name);
0859 
0860     molLogic->SetVisAttributes(vattrs);
0861     pv_placement = new G4PVPlacement(rot_pointer, thisMolecule.fposition, molLogic, name,
0862                                      physicsLogical, false, 0, this->GetOverlaps());
0863   }
0864   else {
0865     // Get information about the phsophate
0866     G4ThreeVector z_direction = phosphate.fposition - thisMolecule.fposition;
0867     G4double separation = z_direction.mag();
0868     G4double z_cut = separation - phosphate.fsize.getX();
0869 
0870     G4ThreeVector local_z = G4ThreeVector(0, 0, 1);
0871     G4ThreeVector z_new = z_direction / z_direction.mag();
0872 
0873     G4ThreeVector ortho = local_z.cross(z_new);
0874     G4double dp = local_z.dot(z_new);
0875     G4double angle = -std::atan2(ortho.mag(), dp);
0876 
0877     G4RotationMatrix first_rot(ortho, angle);
0878     rot_pointer = new G4RotationMatrix(first_rot);
0879 
0880     auto mol = new G4Ellipsoid(name, thisMolecule.fsize.getX(), thisMolecule.fsize.getY(),
0881                                thisMolecule.fsize.getZ(), -thisMolecule.fsize.getZ(), z_cut);
0882 
0883     auto molLogic = new G4LogicalVolume(mol, fpSugarMaterial, name);
0884 
0885     molLogic->SetVisAttributes(vattrs);
0886 
0887     pv_placement = new G4PVPlacement(rot_pointer, thisMolecule.fposition, molLogic, name,
0888                                      physicsLogical, false, 0, this->GetOverlaps());
0889   }
0890   return pv_placement;
0891 }
0892 
0893 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0894 
0895 G4VPhysicalVolume* DNAGeometry::PlaceBase(G4LogicalVolume* physicsLogical,
0896                                           const molecule_t& thisMolecule,
0897                                           const molecule_t& oppBasePair,
0898                                           const molecule_t& nextBasePair, const molecule_t& sugar)
0899 {
0900   std::stringstream ss;
0901   ss << thisMolecule.fname << "-" << thisMolecule.fchain_id << "-" << thisMolecule.fstrand_id << "-"
0902      << thisMolecule.fbase_idx;
0903   G4String name = ss.str();
0904   // Initial Rotation of the sphere
0905   G4RotationMatrix rot;
0906   rot.rotateX(thisMolecule.frotation.getX());
0907   rot.rotateY(thisMolecule.frotation.getY());
0908   rot.rotateZ(thisMolecule.frotation.getZ());
0909 
0910   G4ThreeVector sugar_dirn = sugar.fposition - thisMolecule.fposition;
0911   G4ThreeVector oppBPdirn = oppBasePair.fposition - thisMolecule.fposition;
0912   G4ThreeVector nextBPdirn = nextBasePair.fposition - thisMolecule.fposition;
0913   // Do a cut for the sugar, do a scaling for the base pair
0914   G4double sugar_cut = sugar_dirn.mag() - sugar.fsize.getX();
0915 
0916   G4double base_scaling = oppBPdirn.mag() - oppBasePair.fsize.getY();
0917   base_scaling = std::min(1., base_scaling / thisMolecule.fsize.getY());
0918 
0919   // Find the new constrained values
0920   G4double thisShapeX = base_scaling * thisMolecule.fsize.getX();
0921   G4double thisShapeY = base_scaling * thisMolecule.fsize.getY();
0922   G4double thisShapeZ = base_scaling * thisMolecule.fsize.getZ();
0923   thisShapeZ = std::min(thisShapeZ, 1.7 * angstrom);
0924 
0925   // first rotation to turn the ellipse to run along sugar_dirn
0926   G4ThreeVector z_old = G4ThreeVector(0, 0, 1);
0927   G4ThreeVector z_new = sugar_dirn / sugar_dirn.mag();
0928 
0929   G4ThreeVector ortho = z_old.cross(z_new);
0930   G4double angle = std::acos(z_new.dot(z_old));
0931 
0932   // Now find out quadrant of angle
0933   G4RotationMatrix first_rot;
0934   G4RotationMatrix first_rot_1(ortho, angle);
0935   G4RotationMatrix first_rot_2(ortho, angle + halfpi);
0936   G4RotationMatrix first_rot_3(ortho, angle + pi);
0937   G4RotationMatrix first_rot_4(ortho, angle + 3 * halfpi);
0938 
0939   if (first_rot_1(z_old).dot(z_new) > 0.99) {
0940     first_rot = first_rot_1;
0941   }
0942   else if (first_rot_2(z_old).dot(z_new) > 0.99) {
0943     first_rot = first_rot_2;
0944   }
0945   else if (first_rot_3(z_old).dot(z_new) > 0.99) {
0946     first_rot = first_rot_3;
0947   }
0948   else if (first_rot_4(z_old).dot(z_new) > 0.99) {
0949     first_rot = first_rot_4;
0950   }
0951   else {
0952     G4cout << "Could not identify first rotation" << G4endl;
0953   }
0954 
0955   // second rotation
0956   // Get y vector aligned with next_base_pair_direction through a rotation
0957   // around sugar_dirn
0958   G4ThreeVector y = first_rot(G4ThreeVector(0, 1, 0));
0959   G4ThreeVector up = rot(G4ThreeVector(0, 0, 1));
0960 
0961   G4double interval = 0.1;  // precision for finding theta
0962   G4double theta = -2 * interval;
0963   G4double prev1 = -DBL_MAX;
0964   G4double prev2 = -DBL_MAX;
0965   G4double val = -DBL_MAX;
0966   G4RotationMatrix z_rotation(z_new, theta);
0967   while (theta <= twopi + interval) {
0968     prev2 = prev1;
0969     prev1 = val;
0970     theta += interval;
0971     z_rotation.rotate(interval, z_new);
0972     val = z_rotation(y).dot(up);
0973     if ((prev1 < val) && (prev1 < prev2)) {
0974       // prev1 contains a local minimum
0975       theta -= interval;
0976       break;
0977     }
0978   }
0979   G4RotationMatrix second_rot(first_rot.inverse()(z_new), theta);
0980 
0981   auto rot_pointer = new G4RotationMatrix(second_rot.inverse() * first_rot.inverse());
0982   auto mol = new G4Ellipsoid(name, thisShapeX, thisShapeZ, thisShapeY, -thisShapeY, sugar_cut);
0983 
0984   G4LogicalVolume* molLogic = nullptr;
0985   G4Color color;
0986   if (thisMolecule.fname == "Guanine") {
0987     color = G4Color::Green();
0988     molLogic = new G4LogicalVolume(mol, fpGuanineMaterial, name);
0989   }
0990   else if (thisMolecule.fname == "Adenine") {
0991     color = G4Color::Blue();
0992     molLogic = new G4LogicalVolume(mol, fpAdenineMaterial, name);
0993   }
0994   else if (thisMolecule.fname == "Cytosine") {
0995     color = G4Color::Cyan();
0996     molLogic = new G4LogicalVolume(mol, fpCytosineMaterial, name);
0997   }
0998   else if (thisMolecule.fname == "Thymine") {
0999     color = G4Color::Magenta();
1000     molLogic = new G4LogicalVolume(mol, fpThymineMaterial, name);
1001   }
1002   else {
1003     G4ExceptionDescription errmsg;
1004     errmsg << "Molecule name: " << thisMolecule.fname << " is invalid." << G4endl << "Aborting"
1005            << G4endl;
1006     G4Exception("DNAGeometry::PlaceBase", "Invalid molecule name", FatalException, errmsg);
1007   }
1008 
1009   // Define Vis Attributes
1010   auto vattrs = new G4VisAttributes(color);
1011   vattrs->SetForceSolid(true);
1012   molLogic->SetVisAttributes(vattrs);
1013 
1014   G4VPhysicalVolume* pv_placement =
1015     new G4PVPlacement(rot_pointer, thisMolecule.fposition, molLogic, name, physicsLogical, false, 0,
1016                       this->GetOverlaps());
1017 
1018   return pv_placement;
1019 }
1020 
1021 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1022 
1023 G4VPhysicalVolume* DNAGeometry::PlaceHistone(G4LogicalVolume* physicsLogical,
1024                                              const molecule_t& thisMolecule)
1025 {
1026   G4VPhysicalVolume* pv_placement;
1027 
1028   // Define Vis Attributes
1029   auto vattrs = new G4VisAttributes(G4Color::Gray());
1030   vattrs->SetColor(G4Colour(0, 0, 1, 0.2));
1031   vattrs->SetForceSolid(true);
1032 
1033   // Name:
1034   std::stringstream ss;
1035   ss << thisMolecule.fname << "-" << thisMolecule.fchain_id << "-" << thisMolecule.fstrand_id << "-"
1036      << thisMolecule.fbase_idx;
1037   G4String name = ss.str();
1038 
1039   //// Does the Histone fit in the voxel (assumes square voxels)
1040   // G4double max_coord = std::max(std::abs(thisMolecule.position.getX()),
1041   //                              std::abs(thisMolecule.position.getY()))
1042   // max_coord = std::max(max_coord, std::abs(thisMolecule.position.getZ()));
1043   //// if overlap
1044   // G4double radius = utility::max(thisMolecule.size);
1045   // G4double radiusx=0,radiusx=0,radiusx=0;
1046   // if ((max_coord + radius) > fVoxelSideLength)
1047   //{
1048   //    radius = 0.9*(fVoxelSideLength - max_coord)
1049   //    G4cout<<" Warning : Histone overlap with mother voxel volume "<<G4endl;
1050   //}
1051 
1052   // Initial Rotation of the sphere
1053   G4RotationMatrix rot;
1054   rot.rotateX(thisMolecule.frotation.getX());
1055   rot.rotateY(thisMolecule.frotation.getY());
1056   rot.rotateZ(thisMolecule.frotation.getZ());
1057 
1058   G4RotationMatrix* rot_pointer;
1059   auto mol = new G4Ellipsoid(name, thisMolecule.fsize.getX(), thisMolecule.fsize.getY(),
1060                              thisMolecule.fsize.getZ());
1061 
1062   fHistoneInteractionRadius = thisMolecule.fsize.mag();
1063 
1064   rot_pointer = new G4RotationMatrix(rot.inverse());
1065 
1066   auto molLogic = new G4LogicalVolume(mol, fpHistoneMaterial, name);
1067   molLogic->SetVisAttributes(vattrs);
1068 
1069   pv_placement = new G4PVPlacement(rot_pointer, thisMolecule.fposition, molLogic, name,
1070                                    physicsLogical, false, 0, this->GetOverlaps());
1071 
1072   return pv_placement;
1073 }
1074 // dousatsu --------------------------------------------------------------------
1075 
1076 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1077 
1078 void DNAGeometry::FindNearbyMolecules(const G4LogicalVolume* lv, const G4ThreeVector& localPosition,
1079                                       std::vector<G4VPhysicalVolume*>& daughters_pv_out,
1080                                       double searchRange)
1081 {
1082   const PlacementVolumeInfo* pvInfo = GetPVInfo(lv);
1083   if (pvInfo == nullptr) {
1084     return;
1085   }
1086   // never search more than the radical kill distance
1087   searchRange = std::min(searchRange, fRadicalKillDistance);
1088   OctreeNode* octreeNode = pvInfo->GetOctree();
1089   octreeNode->SearchOctree(localPosition, daughters_pv_out, searchRange);
1090 
1091   //     G4double r = std::pow(std::pow(localPosition.getX(), 2) +
1092   //     std::pow(localPosition.getY(), 2), 0.5)/nm; if (r > 6)
1093   //     {
1094   //       G4cout << "candidates found in radius " << searchRange/nm << " nm : "
1095   //              << daughters_pv_out.size() << G4endl;
1096   //       G4cout << "Radius from center of point: "
1097   //              << r << " nm" << G4endl;
1098   //       G4cout << "----------------------------------------------------" <<
1099   //       G4endl;
1100   //     }
1101 }
1102 
1103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1104 
1105 G4bool DNAGeometry::IsInsideHistone(const G4LogicalVolume* lv,
1106                                     const G4ThreeVector& localPosition) const
1107 {
1108   const PlacementVolumeInfo* pvInfo = GetPVInfo(lv);
1109   if (pvInfo == nullptr) {
1110     G4cout << " Warning : no PV for IsInsideHistone " << G4endl;
1111     return false;
1112   }
1113 
1114   if (!fUseHistoneScav) {
1115     return false;
1116   }
1117 
1118   OctreeNode* octreeNode = pvInfo->GetHistoneOctree();
1119   const std::vector<G4VPhysicalVolume*> result =
1120     octreeNode->SearchOctree(localPosition, fHistoneInteractionRadius);
1121   return !result.empty();
1122 }
1123 
1124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1125 
1126 void DNAGeometry::BasePairIndexTest()
1127 {
1128   G4cout << "------------------------" << G4endl;
1129   G4cout << "Begin BasePairIndexTest" << G4endl;
1130   G4bool pass = true;
1131   for (G4int ii = 0; ii != GetNumberOfChains(); ii++) {
1132     G4cout << "Testing Chain " << ii << G4endl;
1133     std::vector<int64_t> test_vector;  // dousatsu
1134     // std::vector<long long> test_vector;//ORG
1135     int64_t start_idx;
1136     int64_t end_idx;  // dousatsu
1137     // long long start_idx, end_idx;//ORG
1138     for (G4int jj = 0; jj != (G4int)fPlacementTransformations.size(); jj++) {
1139       start_idx = GetStartIdx(jj, ii);
1140       end_idx = GetEndIdx(jj, ii);
1141 
1142       test_vector.push_back(start_idx);
1143       test_vector.push_back(end_idx);
1144     }
1145     for (auto it = test_vector.begin(); it != test_vector.end(); it++) {
1146       if (it == test_vector.begin())  // exception because start of vector
1147       {
1148         if (it + 1 != test_vector.end()) {
1149           if (*it > *(it + 1)) {
1150             G4ExceptionDescription errmsg;
1151             errmsg << "BasePairIndexTest fails at start of test. "
1152                    << "The index " << *(it + 1) << " occurs after " << (*it) << "." << G4endl;
1153             G4Exception("DNAGeometry::BasePairIndexTest", "ERR_TEST_FAIL", FatalException, errmsg);
1154           }
1155         }
1156       }
1157       else if (it == test_vector.end() - 1)  // exception because vec end
1158       {
1159         if (*it < *(it - 1)) {
1160           G4ExceptionDescription errmsg;
1161           errmsg << "BasePairIndexTest fails. "
1162                  << "The index " << *(it - 1) << " occurs before " << (*it) << "." << G4endl;
1163           G4Exception("DNAGeometry::BasePairIndexTest", "ERR_TEST_FAIL", FatalException, errmsg);
1164         }
1165       }
1166       else  // default case
1167       {
1168         if (*it < *(it - 1)) {
1169           G4ExceptionDescription errmsg;
1170           errmsg << "BasePairIndexTest fails. "
1171                  << "The index " << *(it - 1) << " occurs before " << (*it) << "." << G4endl;
1172           G4Exception("DNAGeometry::BasePairIndexTest", "ERR_TEST_FAIL", FatalException, errmsg);
1173         }
1174         if (*it > *(it + 1)) {
1175           G4ExceptionDescription errmsg;
1176           errmsg << "BasePairIndexTest fails. "
1177                  << "The index " << *(it + 1) << " occurs after " << (*it) << "." << G4endl;
1178           G4Exception("DNAGeometry::BasePairIndexTest", "ERR_TEST_FAIL", FatalException, errmsg);
1179         }
1180       }
1181     }
1182   }
1183   if (pass) {
1184     G4cout << "BasePairIndexTest passed" << G4endl;
1185     G4cout << "------------------------" << G4endl;
1186   }
1187 }
1188 
1189 /**
1190  * Global unique id is an integer with the following bases
1191  * moleculeID: base ::molecule::Count
1192  * PlacementIdx: base this->GetNumberOfPlacements()
1193  * StrandID: base 2
1194  * ChainID: base GetNumberOfChains()
1195  * BasePairIdx: base GetMaxBPIdx()
1196  */
1197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1198 
1199 int64_t DNAGeometry::GetGlobalUniqueID(G4VPhysicalVolume* dnavol, const G4VTouchable* touch) const
1200 {
1201   const G4String& dnaName = dnavol->GetName();
1202   std::array<G4String, 4> pv_arr = utility::Get_four_elements(dnaName, '-');
1203   molecule mol = utility::GetMoleculeEnum(pv_arr.at(0));
1204   G4int chainIdx = std::stoi(pv_arr.at(1));
1205   G4int strandIdx = std::stoi(pv_arr.at(2));
1206   int64_t baseIdx = std::stoll(pv_arr.at(3));  // dousatsu
1207   // G4int baseIdx   = std::stoi(pv_arr.at(3));//ORG
1208 
1209   const G4String& placementName = touch->GetVolume()->GetName();
1210   G4int placeIdx = std::stoi(utility::Get_seperated_element(placementName, '-', 1));
1211 
1212   G4int chains = this->GetNumberOfChains();  // dousatsu
1213   G4int placements = this->GetNumberOfPlacements();  // dousatsu
1214   G4int molecules = ::molecule::Count;  // dousatsu
1215   // long long chains = this->GetNumberOfChains();//ORG
1216   // long long placements = this->GetNumberOfPlacements();//ORG
1217   // long long molecules  = ::molecule::Count;//ORG
1218 
1219   chainIdx = this->GetGlobalChain(placeIdx, chainIdx);
1220   // long long globalPair = baseIdx + this->GetStartIdx(placeIdx, chainIdx);
1221   if (this->GetStrandsFlipped(placeIdx)) {
1222     strandIdx = (strandIdx + 1) % 2;
1223   }
1224 
1225   int64_t val1 = mol;
1226   int64_t val2 = molecules * placeIdx;
1227   int64_t val3 = molecules * placements * strandIdx;
1228   int64_t val4 = molecules * placements * 2 * chainIdx;
1229   int64_t val5 = molecules * placements * 2 * chains * baseIdx;
1230   int64_t sum = val1 + val2 + val3 + val4 + val5;
1231   if (sum < 0) {
1232     G4ExceptionDescription errmsg;
1233     errmsg << "v1: " << val1 << G4endl << "v2: " << val2 << G4endl << "v3: " << val3 << G4endl
1234            << "v4: " << val4 << G4endl << "v5: " << val5 << G4endl << "sum: " << sum << G4endl
1235            << "placeIdx: " << placeIdx << G4endl << "placements: " << placements << G4endl
1236            << "strandIdx: " << strandIdx << G4endl << "chainIdx: " << chainIdx << G4endl
1237            << "chains: " << chains << G4endl << "baseIdx: " << baseIdx << G4endl << "mol: " << mol
1238            << G4endl;
1239     G4Exception("DNAGeometry::GetGlobalUniqueID", "ERR_BAD_ID", FatalException, errmsg);
1240   }
1241   return sum;
1242 }
1243 
1244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1245 
1246 // copy of the GetGlobalUniqueID except with integer inputs rather than strings
1247 int64_t DNAGeometry::GetGlobalUniqueIDTest(int64_t mol, int64_t placeIdx, int64_t chainIdx,
1248                                            int64_t strandIdx, int64_t baseIdx) const
1249 {
1250   int64_t chains = this->GetNumberOfChains();  // dousatsu
1251   int64_t placements = this->GetNumberOfPlacements();  // dousatsu
1252   int64_t molecules = ::molecule::Count;  // dousatsu
1253   // long long chains = this->GetNumberOfChains();//ORG
1254   // long long placements = this->GetNumberOfPlacements();//ORG
1255   // long long molecules  = ::molecule::Count;//ORG
1256 
1257   chainIdx = this->GetGlobalChain(placeIdx, chainIdx);
1258   if (this->GetStrandsFlipped(placeIdx)) {
1259     strandIdx = (strandIdx + 1) % 2;
1260   }
1261 
1262   int64_t val1 = mol;
1263   int64_t val2 = molecules * placeIdx;
1264   int64_t val3 = molecules * placements * strandIdx;
1265   int64_t val4 = molecules * placements * 2 * chainIdx;
1266   int64_t val5 = molecules * placements * 2 * chains * baseIdx;
1267   int64_t sum = val1 + val2 + val3 + val4 + val5;
1268   if (sum < 0) {
1269     G4ExceptionDescription errmsg;
1270     errmsg << "v1: " << val1 << G4endl << "v2: " << val2 << G4endl << "v3: " << val3 << G4endl
1271            << "v4: " << val4 << G4endl << "v5: " << val5 << G4endl << "sum: " << sum << G4endl
1272            << " molecules : " << molecules << G4endl << "placeIdx: " << placeIdx << G4endl
1273            << "placements: " << placements << G4endl << "strandIdx: " << strandIdx << G4endl
1274            << "chainIdx: " << chainIdx << G4endl << "chains: " << chains << G4endl
1275            << "baseIdx: " << baseIdx << G4endl << "mol: " << mol << G4endl;
1276     G4Exception("DNAGeometry::GetGlobalUniqueIDTest",  // dousatsu
1277                 "ERR_BAD_ID", FatalException, errmsg);
1278   }
1279   return sum;
1280 }
1281 
1282 void DNAGeometry::UniqueIDTest()
1283 {
1284   G4cout << "Unique ID Test starting." << G4endl;
1285   G4int molecules = ::molecule::Count;
1286 
1287   G4int chains = this->GetNumberOfChains();
1288   G4int placements = this->GetNumberOfPlacements();
1289   G4int basepairs = 100000;  // Maximum number of basepairs to consider per block
1290 
1291   G4bool pass = true;
1292   G4int counter = 0;
1293   for (G4int pment = 0; pment != (placements); ++pment) {
1294     for (G4int chain = 0; chain != (chains); ++chain) {
1295       for (G4int mol = 0; mol != (molecules); ++mol) {
1296         // repeat 10 times along the strand/chain/molecule
1297         for (int64_t ii = 0; ii != 10; ++ii) {
1298           counter++;
1299           G4int strand = G4UniformRand() > 0.5 ? 0 : 1;
1300           G4int basepair = (G4int)basepairs * G4UniformRand();  // dousatsu
1301           int64_t unique_id =  // dousatsu
1302                                // long long unique_id = //ORG
1303             this->GetGlobalUniqueIDTest(mol, pment, chain, strand, basepair);
1304 
1305           G4int real_strand = (this->GetStrandsFlipped(pment)) ?  // ORG
1306                                 (strand + 1) % 2
1307                                                                : strand;
1308           G4int real_chain = this->GetGlobalChain(pment, chain);  // ORG
1309           int64_t real_bp = basepair + this->GetStartIdx(pment, real_chain);  // dousatsu
1310           // long long real_bp = basepair + this->GetStartIdx(pment,
1311           // real_chain);//ORG
1312 
1313           G4int obs_mol = this->GetMoleculeFromUniqueID(unique_id);
1314           G4int obs_placement = this->GetPlacementIndexFromUniqueID(unique_id);
1315           G4int obs_strand = this->GetStrandIDFromUniqueID(unique_id);
1316           G4int obs_chain = this->GetChainIDFromUniqueID(unique_id);
1317           int64_t obs_basepair = this->GetBasePairFromUniqueID(unique_id);  // dousatsu
1318           // long long obs_basepair = this->GetBasePairFromUniqueID
1319           // (unique_id);//ORG
1320 
1321           if ((pment != obs_placement) || (real_strand != obs_strand) || (real_chain != obs_chain)
1322               || (obs_mol != mol) || (real_bp != obs_basepair))
1323           {
1324             pass = false;
1325             G4cerr << "Unique ID test failed at ID: " << unique_id << G4endl;
1326             G4cerr << "real_placement: " << pment << " observed: " << obs_placement << G4endl;
1327             G4cerr << "real_strand: " << real_strand << " observed: " << obs_strand << G4endl;
1328             G4cerr << "real_chain: " << real_chain << " observed: " << obs_chain << G4endl;
1329             G4cerr << "real_mol: " << mol << " observed: " << obs_mol << G4endl;
1330             G4cerr << "real_bp: " << real_bp << " observed: " << obs_basepair << G4endl;
1331           }
1332         }
1333       }
1334     }
1335   }
1336   if (!pass) {
1337     G4Exception("DNAGeometry::UniqueIDTest", "ERR_TEST_FAIL", FatalException,
1338                 "Algorithm to generate unique ID's failed");
1339   }
1340   G4cout << "Unique ID Test completed, after running " << counter << " tests" << G4endl;
1341 }
1342 
1343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1344 
1345 G4Material* DNAGeometry::GetMaterialFromUniqueID(int64_t idx) const
1346 {
1347   molecule mol = GetMoleculeFromUniqueID(idx);
1348   switch (mol) {
1349     case ::molecule::SUGAR:
1350       return fpSugarMaterial;
1351     case ::molecule::PHOSPHATE:
1352       return fpPhosphateMaterial;
1353     case ::molecule::GUANINE:
1354       return fpGuanineMaterial;
1355     case ::molecule::ADENINE:
1356       return fpAdenineMaterial;
1357     case ::molecule::THYMINE:
1358       return fpThymineMaterial;
1359     case ::molecule::CYTOSINE:
1360       return fpCytosineMaterial;
1361     case ::molecule::UNSPECIFIED:
1362       G4cout << "DNAGeometry::GetMaterialFromUniqueID: "
1363              << "Encountered an unspecified material." << G4endl;
1364       return fpMediumMaterial;
1365     default:
1366       G4cout << "DNAGeometry::GetMaterialFromUniqueID: "
1367              << "Encountered a bad material enumerator." << G4endl;
1368       return fpMediumMaterial;
1369   }
1370 }
1371 
1372 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1373 
1374 G4int DNAGeometry::GetPlacementIndexFromUniqueID(int64_t idx) const  // dousatsu
1375 {
1376   // remove molecule
1377   idx -= idx % (int64_t)::molecule::Count;
1378   // remove everything above molecule
1379   idx = idx % (int64_t)(::molecule::Count * this->GetNumberOfPlacements());
1380   return (G4int)idx / (::molecule::Count);
1381 }
1382 
1383 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1384 
1385 G4int DNAGeometry::GetStrandIDFromUniqueID(int64_t idx) const  // dousatsu
1386 {
1387   // remove molecule
1388   idx -= (int64_t)(idx % ::molecule::Count);
1389   // remove placement
1390   idx -= (int64_t)(idx % (::molecule::Count * this->GetNumberOfPlacements()));
1391   // remove everything above strand
1392   idx = (int64_t)(idx % (::molecule::Count * this->GetNumberOfPlacements() * 2));
1393   // recover strand
1394   return (G4int)idx / (::molecule::Count * this->GetNumberOfPlacements());
1395 }
1396 
1397 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1398 
1399 G4int DNAGeometry::GetChainIDFromUniqueID(int64_t idx) const  // dousatsu
1400 {
1401   // remove molecule
1402   idx -= idx % (int64_t)::molecule::Count;
1403   // remove placement
1404   idx -= idx % (int64_t)(::molecule::Count * this->GetNumberOfPlacements());
1405   // remove strand
1406   idx -= idx % (int64_t)(::molecule::Count * this->GetNumberOfPlacements() * 2);
1407   // remove everything above chain
1408   idx =
1409     idx
1410     % (int64_t)(::molecule::Count * this->GetNumberOfPlacements() * 2 * this->GetNumberOfChains());
1411   // recover chain
1412   G4int chain = (G4int)idx / (int64_t)(::molecule::Count * this->GetNumberOfPlacements() * 2);
1413   return chain;
1414 }
1415 
1416 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1417 
1418 int64_t DNAGeometry::GetBasePairFromUniqueID(int64_t idx) const
1419 {
1420   // remove molecule
1421   idx -= idx % ::molecule::Count;
1422 
1423   // remove placement
1424   int64_t placement_ = idx % (int64_t)(::molecule::Count * this->GetNumberOfPlacements());
1425   idx -= placement_;
1426 
1427   // find placement for later
1428   placement_ = placement_ / (int64_t)(::molecule::Count);
1429 
1430   // remove strand
1431   idx -= idx % (int64_t)(::molecule::Count * this->GetNumberOfPlacements() * 2);
1432 
1433   // remove chain
1434   int64_t chain =
1435     idx
1436     % (int64_t)(::molecule::Count * this->GetNumberOfPlacements() * 2 * this->GetNumberOfChains());
1437   idx -= chain;
1438 
1439   // find chain for later
1440   chain = chain / (int64_t)(::molecule::Count * this->GetNumberOfPlacements() * 2);
1441 
1442   // only base pairs left
1443   int64_t bp =
1444     idx
1445     / (int64_t)(::molecule::Count * this->GetNumberOfPlacements() * 2 * this->GetNumberOfChains());
1446 
1447   // use chain and placement to get BP offset
1448   bp = bp + this->GetStartIdx(placement_, chain);
1449   return bp;
1450 }
1451 
1452 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1453 
1454 int64_t DNAGeometry::GetMaxBPIdx() const  // dousatsu
1455 {
1456   int64_t mx = 0;  // dousatsu
1457   for (auto val : std::get<2>(*fPlacementTransformations.end())) {
1458     if (val > mx) {
1459       mx = val;
1460     }
1461   }
1462   return mx;
1463 }
1464 
1465 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1466 
1467 G4int DNAGeometry::GetLocalChain(G4int vol_idx, G4int global_chain) const
1468 {
1469   auto chains = std::get<0>(fPlacementTransformations[vol_idx]);
1470   G4int local = 0;
1471   for (auto chain : chains) {
1472     if (global_chain == chain) {
1473       return local;
1474     }
1475     local++;
1476   }
1477   return -1;
1478 }
1479 
1480 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1481 
1482 G4LogicalVolume* DNAGeometry::GetMatchingPhysVolume(G4LogicalVolume* chem) const
1483 {
1484   auto it = fChemToPhys.find(chem);
1485   if (it != fChemToPhys.end()) {
1486     return it->second;
1487   }
1488   else if (fPhysToChem.find(chem) != fPhysToChem.end()) {
1489     G4cout << "A phyics volume was used to find a chem volume" << G4endl;
1490     return chem;
1491   }
1492   else {
1493     return nullptr;
1494   }
1495 }
1496 
1497 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1498 
1499 G4LogicalVolume* DNAGeometry::GetMatchingChemVolume(G4LogicalVolume* physics) const
1500 {
1501   auto it = fPhysToChem.find(physics);
1502   if (it != fPhysToChem.end()) {
1503     return it->second;
1504   }
1505   else if (fChemToPhys.find(physics) != fChemToPhys.end()) {
1506     G4cout << "A chem volume was used to find a physics volume" << G4endl;
1507     return physics;
1508   }
1509   else {
1510     return nullptr;
1511   }
1512 }
1513 
1514 void DNAGeometry::ChromosomeTest()
1515 {
1516   fpChromosomeMapper->Test();
1517 }
1518 
1519 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1520 
1521 void DNAGeometry::AddChangeMoleculeSize(const G4String& name, const G4ThreeVector& size)
1522 {
1523   G4String lower_name = name;
1524   // lower_name.toLower();
1525   G4StrUtil::to_lower(lower_name);
1526   if (!(fEnableCustomMoleculeSizes)) {
1527     G4cerr << "Custom Molecule Sizes need to be enabled first" << G4endl;
1528   }
1529   else {
1530     if (fMoleculeSizes.count(lower_name) == 0) {
1531       // add molecule
1532       fMoleculeSizes[lower_name] = size;
1533     }
1534     else {
1535       // error message for replacement
1536       G4cout << lower_name << " is already defined in molecule size, I will update" << G4endl
1537              << " this molecule's size from" << fMoleculeSizes[lower_name] << G4endl << " to "
1538              << size << G4endl;
1539       fMoleculeSizes[lower_name] = size;
1540     }
1541   }
1542 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!