| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 /// file: DNAGeometry.cc
28 /// brief:
29 /*
30 * This class builds the DNA geometry. It interacts in a non-standard way
31 * with the DNAWorld class, to create a physical DNA geometry
32 * in a parallel world.
33 *
34 */
35 #include "DNAGeometry.hh"
36
37 #include "DNAGeometryMessenger.hh"
38 #include "OctreeNode.hh"
39 #include "PlacementVolumeInfo.hh"
40 #include "VirtualChromosome.hh"
41
42 #include "G4Box.hh"
43 #include "G4Color.hh"
44 #include "G4Ellipsoid.hh"
45 #include "G4NistManager.hh"
46 #include "G4PVParameterised.hh"
47 #include "G4PVPlacement.hh"
48 #include "G4PhysicalConstants.hh"
49 #include "G4VisAttributes.hh"
50
51 #include <cmath>
52 #include <fstream>
53
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55
56 G4bool compareLVByName::operator()(const G4LogicalVolume* lhs, const G4LogicalVolume* rhs) const
57 {
58 return lhs->GetName() < rhs->GetName();
59 }
60
61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
62
63 DNAGeometry::DNAGeometry()
64 {
65 fpMessenger = new DNAGeometryMessenger(this);
66 fpChromosomeMapper = new ChromosomeMapper();
67
68 // Create Damage Model and add some default values
69 fpDamageModel = new DamageModel();
70 G4NistManager* man = G4NistManager::Instance();
71 fpSugarMaterial = man->FindOrBuildMaterial("G4_DNA_DEOXYRIBOSE");
72 fpMediumMaterial = man->FindOrBuildMaterial("G4_WATER");
73 fpPhosphateMaterial = man->FindOrBuildMaterial("G4_DNA_PHOSPHATE");
74 fpGuanineMaterial = man->FindOrBuildMaterial("G4_DNA_GUANINE");
75 fpAdenineMaterial = man->FindOrBuildMaterial("G4_DNA_ADENINE");
76 fpThymineMaterial = man->FindOrBuildMaterial("G4_DNA_THYMINE");
77 fpCytosineMaterial = man->FindOrBuildMaterial("G4_DNA_CYTOSINE");
78 fpHistoneMaterial = man->FindOrBuildMaterial("G4_WATER");
79 fpDNAPhysicsWorld = new DNAWorld();
80 fpDrawCellAttrs = new G4VisAttributes();
81 fpDrawCellAttrs->SetColor(G4Color(0, 0, 1, 0.2));
82 fpDrawCellAttrs->SetForceSolid(true);
83 // set default molecule sizes
84 fMoleculeSizes["phosphate"] = G4ThreeVector(2.282354, 2.282354, 2.282354) * angstrom;
85 fMoleculeSizes["sugar"] = G4ThreeVector(2.632140, 2.632140, 2.632140) * angstrom;
86 fMoleculeSizes["guanine"] = G4ThreeVector(3.631503, 3.799953, 1.887288) * angstrom;
87 fMoleculeSizes["cytosine"] = G4ThreeVector(3.597341, 3.066331, 1.779361) * angstrom;
88 fMoleculeSizes["adenine"] = G4ThreeVector(3.430711, 3.743504, 1.931958) * angstrom;
89 fMoleculeSizes["thymine"] = G4ThreeVector(4.205943, 3.040448, 2.003359) * angstrom;
90 fMoleculeSizes["histone"] = G4ThreeVector(25, 25, 25) * angstrom;
91 }
92
93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
94
95 DNAGeometry::~DNAGeometry()
96 {
97 delete fpMessenger;
98 }
99
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101
102 void DNAGeometry::BuildDNA(G4LogicalVolume* vol)
103 {
104 // TODO: Add Assertion tests to make sure files have been loaded
105 fpDNAVolumeChem = vol;
106
107 // TODO: Make a complete copy of vol and it's physical volume
108 auto* volAsEllipsoid = dynamic_cast<G4Ellipsoid*>(vol->GetSolid());
109 auto* volAsBox = dynamic_cast<G4Box*>(vol->GetSolid());
110
111 G4VSolid* cloneSolid;
112 if (volAsEllipsoid != nullptr) {
113 cloneSolid = new G4Ellipsoid(*((G4Ellipsoid*)vol->GetSolid()));
114 }
115 else if (volAsBox != nullptr) {
116 cloneSolid = new G4Box(*((G4Box*)vol->GetSolid()));
117 }
118 else {
119 G4ExceptionDescription errmsg;
120 errmsg << "An invalid physical volume shape has been used to hold the DNA volume" << G4endl;
121 G4Exception("MolecularDnaGeometry", "ERR_INVALID_DNA_SHAPE", FatalException, errmsg);
122 return;
123 }
124
125 fpDNAVolumePhys = new G4LogicalVolume(cloneSolid, vol->GetMaterial(),
126 "DNAPhysLV"); // dousatsu
127 FillParameterisedSpace();
128 fpDNAPhysicsWorld->SetDNAVolumePointer(fpDNAVolumePhys);
129 }
130
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132
133 void DNAGeometry::FillParameterisedSpace()
134 {
135 // Create map for each logical volume with its parametrisation
136 // LV, placement_index, copy number, position, rotation
137 std::vector<placement> placementVector;
138 std::map<G4LogicalVolume*, G4int> currentCopyNumber;
139 G4int numberOfChains = 0;
140
141 // read in and create the LV's specified in the macro
142 std::map<G4String, G4LogicalVolume*> namesToLVs;
143 for (const auto& it : fVoxelNames) {
144 G4String thisVoxelName = it.first;
145 G4String thisVoxelFile = it.second;
146
147 auto voxel = LoadVoxelVolume(thisVoxelName, thisVoxelFile);
148 G4LogicalVolume* lv = voxel.first;
149 fInfoMap[lv] = voxel.second;
150 namesToLVs[thisVoxelName] = lv;
151 currentCopyNumber[lv] = 0;
152 numberOfChains = this->GetPVInfo(lv)->GetNumberOfChains();
153 if (!(numberOfChains == 1 || numberOfChains == 4 || numberOfChains == 8)) {
154 G4cout << "Geometry contains volumes with none of 1/4/8"
155 << " chains. This can cause errors" << G4endl;
156 }
157 }
158
159 // Go through the voxel types to build a placement vector
160 for (unsigned ii = 0; ii < fVoxelTypes.size(); ii++) {
161 G4LogicalVolume* lv = namesToLVs[fVoxelTypes[ii]];
162
163 placement thisPlacement = std::make_tuple(lv, fVoxelIndices[ii], currentCopyNumber[lv]++,
164 fVoxelPositions[ii], fVoxelRotations[ii]);
165 placementVector.push_back(thisPlacement);
166 }
167
168 // Do each placement.
169 std::map<G4String, int64_t> histonesPerChromosome; // dousatsu
170 std::map<G4String, int64_t> basePairsPerChromosome; // dousatsu
171 // std::map<G4String, long long> basePairsPerChromosome;//ORG
172 for (auto it = placementVector.begin(); it != placementVector.end(); it++) {
173 G4LogicalVolume* thisLogical = std::get<0>(*it);
174 if (thisLogical == nullptr) {
175 G4Exception("DNAGeometry::FillParameterisedSpace", "ERR_LOG_VOL_EXISTS_FALSE", FatalException,
176 "At least one of the placement names"
177 " specified doesn't exist");
178 }
179 G4int thisIndex = std::get<1>(*it);
180 G4int thisCopyNo = std::get<2>(*it);
181
182 G4ThreeVector thisPosition = std::get<3>(*it);
183 G4ThreeVector thisRotation = std::get<4>(*it);
184 std::stringstream ss;
185 ss << thisLogical->GetName() << "-" << thisIndex;
186
187 auto inv = new G4RotationMatrix();
188 inv->rotateX(thisRotation.getX());
189 inv->rotateY(thisRotation.getY());
190 inv->rotateZ(thisRotation.getZ());
191 auto pRot = new G4RotationMatrix(inv->inverse());
192 delete inv;
193
194 G4bool isPhysicallyPlaced = false;
195 VirtualChromosome* thisChromosome = this->GetChromosomeMapper()->GetChromosome(thisPosition);
196 if (thisChromosome != nullptr) {
197 isPhysicallyPlaced = true;
198 // Place for Physics
199 new G4PVPlacement(pRot, thisPosition, thisLogical, ss.str(), fpDNAVolumePhys, false,
200 thisCopyNo, this->GetOverlaps());
201 // Place for Chemistry
202 new G4PVPlacement(pRot, thisPosition, this->GetMatchingChemVolume(thisLogical), ss.str(),
203 fpDNAVolumeChem, false, thisCopyNo, this->GetOverlaps());
204
205 // dousatsu==============
206 if (histonesPerChromosome.find(thisChromosome->GetName()) == histonesPerChromosome.end()) {
207 histonesPerChromosome[thisChromosome->GetName()] = 0;
208 }
209 histonesPerChromosome[thisChromosome->GetName()] +=
210 this->GetPVInfo(thisLogical)->GetTotalHistones();
211 // dousatsu==============
212
213 if (basePairsPerChromosome.find(thisChromosome->GetName()) == basePairsPerChromosome.end()) {
214 basePairsPerChromosome[thisChromosome->GetName()] = 0;
215 }
216 basePairsPerChromosome[thisChromosome->GetName()] +=
217 this->GetPVInfo(thisLogical)->GetTotalBasePairs();
218 }
219
220 // add the placement to the local register, even if it isn't in a
221 // chromosome, as this tracks base pair index number.
222 if (numberOfChains == 4 || numberOfChains == 8) {
223 AddFourChainPlacement(it, placementVector.begin(), isPhysicallyPlaced);
224 }
225 else if (numberOfChains == 1) {
226 AddSingleChainPlacement(it, placementVector.begin(), isPhysicallyPlaced);
227 }
228 else {
229 G4ExceptionDescription errmsg;
230 errmsg << "Having none of 1/4/8 chains in a placement is not "
231 << "supported, the application will crash" << G4endl;
232 G4Exception("DNAGeometry::FillParameterisedSpace", "Number of chains not supported",
233 FatalException, errmsg);
234 }
235 }
236
237 // dousatsu------------------------------
238 G4cout << "Histones placed per chromosome:" << G4endl;
239 G4cout << "Chromosome Histones" << G4endl;
240 G4cout << "-----------------------------------" << G4endl;
241 for (auto& it : histonesPerChromosome) {
242 G4cout << it.first.c_str() << " : " << it.second << G4endl;
243 }
244 // dousatsu------------------------------
245
246 G4cout << "Base Pairs placed per chromosome:" << G4endl;
247 G4cout << "Chromosome Base Pairs" << G4endl;
248 G4cout << "-----------------------------------" << G4endl;
249 for (auto& it : basePairsPerChromosome) {
250 G4cout << it.first.c_str() << " : " << it.second << G4endl;
251 }
252 G4cout << "-------------------------------" << G4endl;
253
254 if (this->GetVerbosity() > 0) {
255 this->PrintOctreeStats();
256 }
257 }
258
259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
260
261 void DNAGeometry::AddSingleChainPlacement(const std::vector<placement>::iterator it,
262 const std::vector<placement>::iterator,
263 G4bool isPhysicallyPlaced)
264 {
265 G4LogicalVolume* thisLogical = std::get<0>(*it);
266 G4bool twist = fVoxelTwist[thisLogical->GetName()];
267 AddNewPlacement(thisLogical, {{0, 1, 2, 3, 4, 5, 6, 7}}, twist, isPhysicallyPlaced);
268 }
269
270 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
271
272 void DNAGeometry::AddFourChainPlacement(const std::vector<placement>::iterator it,
273 const std::vector<placement>::iterator begin,
274 G4bool isPhysicallyPlaced)
275 {
276 G4LogicalVolume* thisLogical = std::get<0>(*it);
277 G4int thisIndex = std::get<1>(*it);
278 // Use info to build the placement transforms.
279 std::array<int, 8> global_chain{};
280 G4bool twist = fVoxelTwist[thisLogical->GetName()];
281 if (it == begin) {
282 global_chain = {{0, 1, 2, 3, 4, 5, 6, 7}};
283 }
284 else {
285 // work out global chain
286 placement previous = *(it - 1);
287 G4ThreeVector oldRotation = std::get<4>(previous);
288 G4ThreeVector thisRotation = std::get<4>(*it);
289
290 G4RotationMatrix rotnew;
291 rotnew.rotateX(thisRotation.getX());
292 rotnew.rotateY(thisRotation.getY());
293 rotnew.rotateZ(thisRotation.getZ());
294
295 rotnew.rectify();
296
297 G4RotationMatrix rotold;
298 rotold.rotateX(oldRotation.getX());
299 rotold.rotateY(oldRotation.getY());
300 rotold.rotateZ(oldRotation.getZ());
301
302 G4ThreeVector zdiff = rotold.colZ() - rotnew.colZ();
303 if (zdiff.mag2() > 0.1) {
304 rotold = rotold.rotate(halfpi, rotold.colY());
305 }
306 rotold.rectify();
307 // rotnew.colz == rotold.colz now.
308 // There exists now a rotation around rotnew.colZ()
309 // that transforms rotold to rotnew.
310 G4RotationMatrix relative = rotnew * rotold.inverse();
311 relative.rectify();
312 G4ThreeVector axis = relative.getAxis();
313 G4double angle = relative.getDelta();
314 // get the sign of the rotation, and correct for quadrant.
315 G4int sign = (axis.dot(rotnew.colZ()) >= 0) ? 1 : -1;
316 if (sign < 0) {
317 angle = 2 * pi - angle;
318 }
319 G4int zrots = ((G4int)(2.05 * angle / pi)) % 4; // ORG
320
321 global_chain = {{(this->GetGlobalChain(thisIndex - 1, 0) + zrots) % 4,
322 (this->GetGlobalChain(thisIndex - 1, 1) + zrots) % 4,
323 (this->GetGlobalChain(thisIndex - 1, 2) + zrots) % 4,
324 (this->GetGlobalChain(thisIndex - 1, 3) + zrots) % 4,
325 4 + (this->GetGlobalChain(thisIndex - 1, 4) + zrots) % 4,
326 4 + (this->GetGlobalChain(thisIndex - 1, 5) + zrots) % 4,
327 4 + (this->GetGlobalChain(thisIndex - 1, 6) + zrots) % 4,
328 4 + (this->GetGlobalChain(thisIndex - 1, 7) + zrots) % 4}};
329 twist = false;
330 }
331 this->AddNewPlacement(thisLogical, global_chain, twist, isPhysicallyPlaced);
332 }
333
334 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
335
336 void DNAGeometry::AddNewPlacement(const G4LogicalVolume* lv, std::array<int, 8> global_chain,
337 G4bool twist, G4bool isPhysicallyPlaced)
338 {
339 placement_transform pt;
340 std::array<int64_t, 8> start{}; // dousatsu
341 std::array<int64_t, 8> end{}; // dousatsu
342 // std::array<long long, 8> start;//ORG
343 // std::array<long long, 8> end;//ORG
344 if (fPlacementTransformations.empty()) {
345 start = {{0, 0, 0, 0, 0, 0, 0, 0}};
346 end = {{this->GetPVInfo(lv)->GetPairsOnChain(global_chain[0]),
347 this->GetPVInfo(lv)->GetPairsOnChain(global_chain[1]),
348 this->GetPVInfo(lv)->GetPairsOnChain(global_chain[2]),
349 this->GetPVInfo(lv)->GetPairsOnChain(global_chain[3])}};
350 }
351 else {
352 placement_transform previous = fPlacementTransformations.back();
353 std::array<int64_t, 8> previousEnd = std::get<2>(previous); // dousatsu
354 // std::array<long long, 8> previousEnd = std::get<2>(previous);//ORG
355 start = {{previousEnd[0], previousEnd[1], previousEnd[2], previousEnd[3], previousEnd[4],
356 previousEnd[5], previousEnd[6], previousEnd[7]}};
357 end = {{previousEnd[0] + this->GetPVInfo(lv)->GetPairsOnChain(global_chain[0]),
358 previousEnd[1] + this->GetPVInfo(lv)->GetPairsOnChain(global_chain[1]),
359 previousEnd[2] + this->GetPVInfo(lv)->GetPairsOnChain(global_chain[2]),
360 previousEnd[3] + this->GetPVInfo(lv)->GetPairsOnChain(global_chain[3]),
361 previousEnd[4] + this->GetPVInfo(lv)->GetPairsOnChain(global_chain[4]),
362 previousEnd[5] + this->GetPVInfo(lv)->GetPairsOnChain(global_chain[5]),
363 previousEnd[6] + this->GetPVInfo(lv)->GetPairsOnChain(global_chain[6]),
364 previousEnd[7] + this->GetPVInfo(lv)->GetPairsOnChain(global_chain[7])}};
365 }
366 pt = std::make_tuple(global_chain, start, end, twist, isPhysicallyPlaced);
367 fPlacementTransformations.push_back(pt);
368 }
369 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
370
371 std::pair<G4LogicalVolume*, PlacementVolumeInfo*>
372 DNAGeometry::LoadVoxelVolume(const G4String& volumeName, const G4String& filename)
373 {
374 if (this->GetVerbosity() > 0) {
375 G4cout << "Loading Voxel: " << volumeName << G4endl;
376 }
377
378 auto thisInfo = new PlacementVolumeInfo();
379
380 G4double vxdim = 0.5 * fVoxelSideLength.getX();
381 G4double vydim = 0.5 * fVoxelSideLength.getY();
382 G4double vzdim = 0.5 * fVoxelSideLength.getZ();
383
384 // Make Volumes, and put them into the maps for phys/chem comparison
385 auto physicsPhysical = new G4Box(volumeName, vxdim, vydim, vzdim);
386 auto physicsLogical = new G4LogicalVolume(physicsPhysical, fpMediumMaterial, volumeName);
387 auto chemPhysical = new G4Box(volumeName, vxdim, vydim, vzdim);
388 auto chemLogical = new G4LogicalVolume(chemPhysical, fpMediumMaterial, volumeName);
389
390 if (this->GetDrawCellVolumes()) {
391 chemLogical->SetVisAttributes(fpDrawCellAttrs);
392 }
393
394 fChemToPhys[chemLogical] = physicsLogical;
395 fPhysToChem[physicsLogical] = chemLogical;
396
397 physicsLogical->SetSmartless(this->GetSmartless());
398
399 auto thisOctree = new OctreeNode(G4ThreeVector(0, 0, 0), G4ThreeVector(vxdim, vydim, vzdim), 1);
400 auto thisHistoneOctree =
401 new OctreeNode(G4ThreeVector(0, 0, 0), G4ThreeVector(vxdim, vydim, vzdim), 1);
402
403 // open and load file
404 if (!utility::Path_exists(filename)) {
405 G4ExceptionDescription errmsg;
406 errmsg << "The file: " << filename << " could not be found." << G4endl;
407 G4Exception("DNAGeometry::LoadVoxelVolume", "File Not Found", FatalException, errmsg);
408 }
409 std::ifstream infile(filename, std::ifstream::in);
410 G4String currentline;
411 G4String lower_name;
412 G4ThreeVector pos;
413 G4ThreeVector rot;
414 G4ThreeVector mol_size;
415 std::vector<molecule_t> molecules;
416 std::vector<molecule_t> histones;
417 G4bool file_contains_size = false;
418 G4int line_fields = 0;
419 G4int line_number = 0;
420 G4int uncommented_line_number = 0;
421
422 if (infile.is_open()) {
423 while (getline(infile, currentline)) {
424 if (currentline[0] != '#') {
425 std::vector<G4String> line = utility::Split(currentline, ' ');
426 // validation
427 if (uncommented_line_number == 0) {
428 line_fields = line.size();
429 file_contains_size = (line_fields == 14);
430 }
431 if ((line_fields != 10) && (line_fields != 14)) {
432 G4ExceptionDescription errmsg;
433 errmsg << filename << " should have 10 or 14 fields, only got " << line_fields << G4endl;
434 G4Exception("DNAGeometry::LoadVoxelVolume", "BadFileFormat", FatalErrorInArgument,
435 errmsg);
436 }
437 if (line_fields != (G4int)line.size()) {
438 G4ExceptionDescription errmsg;
439 errmsg << "Unexpected number of fields on line " << line_number << " of file " << filename
440 << G4endl << "Expected " << line_fields << " fields, but got " << line.size()
441 << G4endl;
442 G4Exception("DNAGeometry::LoadVoxelVolume", "BadFileFormat", FatalErrorInArgument,
443 errmsg);
444 }
445 molecule_t thisMolecule;
446
447 thisMolecule.fname = (G4String)line.at(0);
448 thisMolecule.fshape = (G4String)line.at(1);
449 thisMolecule.fchain_id = (G4int)std::stoi(line.at(2));
450 thisMolecule.fstrand_id = (G4int)std::stoi(line.at(3));
451 thisMolecule.fbase_idx = (int64_t)std::stoll(line.at(4));
452
453 lower_name = (G4String)line.at(0);
454 G4StrUtil::to_lower(lower_name);
455 // lower_name.toLower();
456
457 G4int size_offset = (file_contains_size) ? 3 : 0;
458 if (file_contains_size) {
459 mol_size =
460 G4ThreeVector(std::stod(line.at(5)) * angstrom, std::stod(line.at(6)) * angstrom,
461 std::stod(line.at(7)) * angstrom);
462 }
463 else {
464 if (fMoleculeSizes.count(lower_name) == 0) {
465 G4ExceptionDescription errmsg;
466 errmsg << "A molecule size has not been defined for " << lower_name << G4endl;
467 G4Exception("DNAGeometry::LoadVoxelVolume", "MoleculeNotDefined", FatalErrorInArgument,
468 errmsg);
469 }
470 mol_size = fMoleculeSizes.at(lower_name);
471 }
472 pos = G4ThreeVector(std::stod(line.at(5 + size_offset)) * angstrom,
473 std::stod(line.at(6 + size_offset)) * angstrom,
474 std::stod(line.at(7 + size_offset)) * angstrom);
475
476 rot =
477 G4ThreeVector(std::stod(line.at(8 + size_offset)), std::stod(line.at(9 + size_offset)),
478 std::stod(line.at(10 + size_offset)));
479
480 thisMolecule.fsize = mol_size;
481 thisMolecule.fposition = pos;
482 thisMolecule.frotation = rot;
483
484 if (lower_name == "histone") {
485 histones.push_back(thisMolecule);
486 }
487 else {
488 molecules.push_back(thisMolecule);
489 }
490 uncommented_line_number++;
491 }
492 line_number++;
493 }
494 infile.close();
495 }
496 else {
497 G4ExceptionDescription errmsg;
498 errmsg << "Voxel Position file read error" << G4endl << "Name: " << volumeName << G4endl
499 << "Filename: " << filename << G4endl << "Format should be: "
500 << "# NAME SHAPE CHAIN_ID STRAND_ID BP_INDEX SIZE_X SIZE_Y "
501 << "SIZE_Z POS_X POS_Y POS_Z ROT_X ROT_Y ROT_Z" << G4endl << "or" << G4endl
502 << "# NAME SHAPE CHAIN_ID STRAND_ID BP_INDEX POS_X POS_Y POS_Z "
503 "ROT_X ROT_Y ROT_Z"
504 << G4endl
505 << "Separator is space, lines with '#' as the first character are "
506 "commented"
507 << G4endl;
508
509 G4Exception("DNAGeometry::LoadVoxelVolume", "BadFile", FatalErrorInArgument, errmsg);
510 }
511
512 // We can place the elements now, specify an order that they'll be placed
513 // for the union solid cuts.
514 G4VPhysicalVolume* pv_placement = nullptr;
515
516 // Place Histones
517 G4int his_arr_size = histones.size();
518 for (G4int ii = 0; ii != his_arr_size; ++ii) {
519 molecule_t thisMolecule = histones[ii];
520 if (this->GetVerbosity() > 4) {
521 G4cout << "Placing molecule " << ii << ": " << thisMolecule.fname << G4endl;
522 }
523 // Identify number of histones in voxel volume
524 pv_placement = PlaceHistone(physicsLogical, thisMolecule);
525 thisHistoneOctree->AddPhysicalVolume(pv_placement);
526 }
527 thisInfo->SetNHistones(his_arr_size);
528 thisInfo->SetHistoneOctree(thisHistoneOctree);
529
530 // Place Molecules
531 G4int mol_arr_size = molecules.size();
532 for (G4int ii = 0; ii != mol_arr_size; ii++) {
533 molecule_t thisMolecule = molecules[ii];
534 if (this->GetVerbosity() > 4) {
535 G4cout << "Placing molecule " << ii << ": " << thisMolecule.fname << G4endl;
536 }
537 // Identify number of molecules on chain
538 if (thisInfo->GetPairsOnChain(thisMolecule.fchain_id) < (thisMolecule.fbase_idx + 1)) {
539 thisInfo->SetPairsOnChain(thisMolecule.fchain_id, thisMolecule.fbase_idx + 1);
540 }
541
542 // Place Volume
543 if (thisMolecule.fname == "Sugar") {
544 molecule_t phosphate = molecules[ii - 1];
545 // G4int phosph_idx = (thisMolecule.strand_id == 0) ? ii - 1 : ii - 1;
546 G4int phosph_idx = ii - 1;
547 if ((phosph_idx >= 0) && (phosph_idx < (G4int)mol_arr_size)) {
548 phosphate = molecules[phosph_idx];
549 if (phosphate.fchain_id != thisMolecule.fchain_id) {
550 // There is no next sugar, just the voxel wall
551 phosphate = molecule_t();
552 phosphate.fname = "wall";
553 }
554 }
555 else {
556 // There is no next sugar, just the voxel wall
557 phosphate = molecule_t();
558 phosphate.fname = "wall";
559 }
560 pv_placement = PlaceSugar(physicsLogical, thisMolecule, phosphate);
561 }
562 else if (thisMolecule.fname == "Phosphate") {
563 molecule_t sugar;
564 G4int sugar_idx = (thisMolecule.fstrand_id == 0) ? ii + 7 : ii - 5;
565 if ((sugar_idx >= 0) && (sugar_idx < mol_arr_size)) {
566 sugar = molecules[sugar_idx];
567 if (sugar.fchain_id != thisMolecule.fchain_id) {
568 // There is no next sugar, just the voxel wall
569 sugar = molecule_t();
570 sugar.fname = "wall";
571 }
572 }
573 else {
574 // There is no next sugar, just the voxel wall
575 sugar = molecule_t();
576 sugar.fname = "wall";
577 }
578 pv_placement = PlacePhosphate(physicsLogical, thisMolecule, sugar);
579 }
580 else if ((thisMolecule.fname == "Guanine") || (thisMolecule.fname == "Adenine")
581 || (thisMolecule.fname == "Cytosine") || (thisMolecule.fname == "Thymine"))
582 {
583 molecule_t sugar = molecules[ii - 1];
584 molecule_t oppBasePair;
585 molecule_t nextBasePair;
586
587 G4int nextBPidx;
588 G4bool endOfChain = false;
589 G4bool nextBPexists;
590 if (thisMolecule.fstrand_id == 0) {
591 oppBasePair = molecules[ii + 3];
592 nextBPexists = (ii + 6 < mol_arr_size);
593 if (nextBPexists) {
594 endOfChain = (thisMolecule.fchain_id != molecules[ii + 6].fchain_id);
595 }
596 nextBPidx = (nextBPexists && !endOfChain) ? ii + 6 : ii - 6;
597 }
598 else {
599 oppBasePair = molecules[ii - 3];
600 nextBPexists = (ii - 6 >= 0);
601 if (nextBPexists) {
602 endOfChain = (thisMolecule.fchain_id != molecules[ii - 6].fchain_id);
603 }
604 nextBPidx = (nextBPexists && !endOfChain) ? ii - 6 : ii + 6;
605 }
606 nextBasePair = molecules[nextBPidx];
607 pv_placement = PlaceBase(physicsLogical, thisMolecule, oppBasePair, nextBasePair, sugar);
608 }
609 else {
610 G4ExceptionDescription errmsg;
611 errmsg << "Molecule name: " << thisMolecule.fname << " is invalid." << G4endl << "Aborting"
612 << G4endl;
613 G4Exception("DNAGeometry::LoadVoxelVolume", "Invalid molecule name", FatalException, errmsg);
614 }
615 if (this->GetOverlaps()) {
616 if (pv_placement->CheckOverlaps()) {
617 // NOTE: if there are overlaps, the simulation will probably
618 // crash. This is the desired behaviour (Fail Loud)
619 G4cout << "Overlap, translation: " << pv_placement->GetFrameTranslation() << G4endl;
620 G4cout << "Overlap, rotationX: " << pv_placement->GetFrameRotation()->colX() << G4endl;
621 G4cout << "Overlap, rotationY: " << pv_placement->GetFrameRotation()->colY() << G4endl;
622 G4cout << "Overlap, rotationZ: " << pv_placement->GetFrameRotation()->colZ() << G4endl;
623 }
624 }
625 if (this->GetDrawCellVolumes()) {
626 // do not draw DNA, as we are drawing the cell volume
627 pv_placement->GetLogicalVolume()->SetVisAttributes(G4VisAttributes::GetInvisible());
628 }
629 thisOctree->AddPhysicalVolume(pv_placement);
630 }
631 thisInfo->SetOctree(thisOctree);
632 return std::make_pair(physicsLogical, thisInfo);
633 }
634
635 void DNAGeometry::FillVoxelVectors()
636 {
637 G4String filename = fFractalCurveFile;
638 if (!utility::Path_exists(filename)) {
639 G4ExceptionDescription errmsg;
640 errmsg << "The file: " << filename << " could not be found." << G4endl;
641 G4Exception("DNAGeometry::FillVoxelVectors", "File Not Found", FatalException, errmsg);
642 }
643 std::ifstream infile(filename, std::ifstream::in);
644 G4String currentline;
645 G4ThreeVector pos;
646 G4ThreeVector rot;
647
648 fVoxelTypes.clear();
649 fVoxelIndices.clear();
650 fVoxelPositions.clear();
651 fVoxelRotations.clear();
652
653 if (infile.is_open()) {
654 while (getline(infile, currentline)) {
655 if (currentline[0] != '#') {
656 try {
657 G4int pi_or_one = fAnglesAsPi ? pi : 1;
658 std::vector<G4String> line = utility::Split(currentline, ' ');
659 fVoxelIndices.push_back(std::stoi(line.at(0)));
660 fVoxelTypes.push_back((G4String)line.at(1));
661 pos = G4ThreeVector(std::stod(line.at(2)) * fFractalScaling.getX(),
662 std::stod(line.at(3)) * fFractalScaling.getY(),
663 std::stod(line.at(4)) * fFractalScaling.getZ());
664 rot = G4ThreeVector(std::stod(line.at(5)), std::stod(line.at(6)), std::stod(line.at(7)));
665 fVoxelPositions.push_back(pos);
666 fVoxelRotations.push_back(rot * pi_or_one);
667 }
668 catch (int exception) {
669 G4cout << "Voxel Position file read error" << G4endl;
670 G4cout << "Got Line: " << currentline << G4endl;
671 G4cout << "Format should be: "
672 << "IDX KIND POS_X POS_Y POS_Z EUL_PSI EUL_THETA"
673 << " EUL_PHI" << G4endl;
674 G4cout << "Separator is space, lines with '#' as the first "
675 << "character are commented" << G4endl;
676 std::exit(1);
677 }
678 }
679 }
680 infile.close();
681 }
682 else {
683 G4cout << "Voxel Position file read error on open" << G4endl;
684 G4cout << "Format should be: "
685 << "IDX KIND POS_X POS_Y POS_Z EUL_PSI EUL_THETA EUL_PHI" << G4endl;
686 G4cout << "Separator is space, lines with '#' as the first "
687 << "character are commented" << G4endl;
688 }
689 }
690
691 OctreeNode* DNAGeometry::GetTopOctreeNode(G4LogicalVolume* lv) const
692 {
693 const PlacementVolumeInfo* info = GetPVInfo(lv);
694 return (info == nullptr) ? nullptr : info->GetOctree();
695 }
696
697 const PlacementVolumeInfo* DNAGeometry::GetPVInfo(const G4LogicalVolume* lv) const
698 {
699 if (fInfoMap.find(lv) != fInfoMap.end()) {
700 return fInfoMap.at(lv);
701 }
702 else if (fChemToPhys.find(lv) != fChemToPhys.end()) {
703 return fInfoMap.at(fChemToPhys.at(lv));
704 }
705 else {
706 return nullptr;
707 }
708 }
709
710 void DNAGeometry::PrintOctreeStats()
711 {
712 G4cout << "Name Vols EndNodes Max(vols/node)" << G4endl;
713 G4cout << "----------------------------------------------" << G4endl;
714 OctreeNode* currentOctree;
715 for (auto& it : fInfoMap) {
716 char buffer[80];
717 currentOctree = it.second->GetOctree();
718 G4int vols = currentOctree->GetContents().size();
719 const char* lvname = it.first->GetName().c_str();
720 G4int mvols = currentOctree->GetMaxContents();
721 G4int nodes = currentOctree->GetNumberOfTerminalNodes();
722 std::snprintf(buffer, 80, "%-12s %-7i %-7i %-9i", lvname, vols, nodes, mvols);
723 G4cout << buffer << G4endl;
724 }
725 G4cout << "-------------------------------------------" << G4endl;
726 }
727
728 ////////////////////////////////////////////////////////////////////////////////
729 ////////////////////////////////////////////////////////////////////////////////
730 ////////////////////////////////////////////////////////////////////////////////
731
732 G4VPhysicalVolume* DNAGeometry::PlacePhosphate(G4LogicalVolume* physicsLogical,
733 const molecule_t& thisMolecule,
734 const molecule_t& sugar)
735 {
736 // Define Vis Attributes
737 auto vattrs = new G4VisAttributes(G4Color::Yellow());
738 vattrs->SetForceSolid(true);
739
740 std::stringstream ss;
741 ss << thisMolecule.fname << "-" << thisMolecule.fchain_id << "-" << thisMolecule.fstrand_id << "-"
742 << thisMolecule.fbase_idx;
743 G4String name = ss.str();
744 // Initial Rotation of the sphere
745 G4RotationMatrix rot;
746 rot.rotateX(thisMolecule.frotation.getX());
747 rot.rotateY(thisMolecule.frotation.getY());
748 rot.rotateZ(thisMolecule.frotation.getZ());
749
750 G4VPhysicalVolume* pv_placement;
751 G4RotationMatrix* rot_pointer;
752 if (sugar.fname == "wall") {
753 if (this->GetVerbosity() >= 2) {
754 G4cout << "Wall placement" << G4endl;
755 }
756 G4ThreeVector abspos(std::abs(thisMolecule.fposition.getX()),
757 std::abs(thisMolecule.fposition.getY()),
758 std::abs(thisMolecule.fposition.getZ()));
759 G4ThreeVector overlaps = abspos + thisMolecule.fsize - fVoxelSideLength;
760
761 G4Ellipsoid* mol;
762 // NOTE: will work to cut molecules at boundary only when the chain is
763 // leaving the box.
764 G4double z_cut = -1 * utility::Min(overlaps);
765 if (z_cut < 0) // all molecules fit
766 {
767 mol = new G4Ellipsoid(name, thisMolecule.fsize.getX(), thisMolecule.fsize.getY(),
768 thisMolecule.fsize.getZ());
769 }
770 else // need to cut molecule to fit in box
771 {
772 mol = new G4Ellipsoid(name, thisMolecule.fsize.getX(), thisMolecule.fsize.getY(),
773 thisMolecule.fsize.getZ(), -thisMolecule.fsize.getZ(), z_cut);
774 }
775
776 rot_pointer = new G4RotationMatrix(rot.inverse());
777 auto molLogic = new G4LogicalVolume(mol, fpPhosphateMaterial, name);
778
779 molLogic->SetVisAttributes(vattrs);
780 pv_placement = new G4PVPlacement(rot_pointer, thisMolecule.fposition, molLogic, name,
781 physicsLogical, false, 0, this->GetOverlaps());
782 }
783 else {
784 // Case where there is a sugar molecule that we would intersect with
785 G4ThreeVector z_direction = sugar.fposition - thisMolecule.fposition;
786 G4double separation = z_direction.mag();
787 G4double z_cut = separation - sugar.fsize.getX();
788
789 G4ThreeVector local_z = G4ThreeVector(0, 0, 1);
790 G4ThreeVector z_new = z_direction / z_direction.mag();
791
792 G4ThreeVector ortho = local_z.cross(z_new);
793 G4double dp = local_z.dot(z_new);
794 G4double angle = -std::atan2(ortho.mag(), dp);
795
796 G4RotationMatrix second_rot(ortho, angle);
797 rot_pointer = new G4RotationMatrix(second_rot);
798
799 auto mol = new G4Ellipsoid(name, thisMolecule.fsize.getX(), thisMolecule.fsize.getY(),
800 thisMolecule.fsize.getZ(), -thisMolecule.fsize.getZ(), z_cut);
801
802 auto molLogic = new G4LogicalVolume(mol, fpPhosphateMaterial, name);
803
804 molLogic->SetVisAttributes(vattrs);
805 pv_placement = new G4PVPlacement(rot_pointer, thisMolecule.fposition, molLogic, name,
806 physicsLogical, false, 0, this->GetOverlaps());
807 }
808
809 return pv_placement;
810 }
811
812 G4VPhysicalVolume* DNAGeometry::PlaceSugar(G4LogicalVolume* physicsLogical,
813 const molecule_t& thisMolecule,
814 const molecule_t& phosphate)
815 {
816 // Define Vis Attributes
817 auto vattrs = new G4VisAttributes(G4Color::Red());
818 vattrs->SetForceSolid(true);
819
820 // Name:
821 std::stringstream ss;
822 ss << thisMolecule.fname << "-" << thisMolecule.fchain_id << "-" << thisMolecule.fstrand_id << "-"
823 << thisMolecule.fbase_idx;
824 G4String name = ss.str();
825 // Initial Rotation of the sphere
826 G4RotationMatrix rot;
827 rot.rotateX(thisMolecule.frotation.getX());
828 rot.rotateY(thisMolecule.frotation.getY());
829 rot.rotateZ(thisMolecule.frotation.getZ());
830
831 G4VPhysicalVolume* pv_placement;
832 G4RotationMatrix* rot_pointer;
833 if (phosphate.fname == "wall") {
834 if (this->GetVerbosity() >= 2) {
835 G4cout << "Wall placement" << G4endl;
836 }
837 G4ThreeVector abspos(std::abs(thisMolecule.fposition.getX()),
838 std::abs(thisMolecule.fposition.getY()),
839 std::abs(thisMolecule.fposition.getZ()));
840 G4ThreeVector overlaps = abspos + thisMolecule.fsize - fVoxelSideLength;
841
842 G4Ellipsoid* mol;
843 // NOTE: will work to cut molecules at boundary only when the chain is
844 // leaving the box.
845 G4double z_cut = -1 * utility::Min(overlaps);
846 if (z_cut < 0) // all molecules fit
847 {
848 mol = new G4Ellipsoid(name, thisMolecule.fsize.getX(), thisMolecule.fsize.getY(),
849 thisMolecule.fsize.getZ());
850 }
851 else // need to cut molecule to fit in box
852 {
853 mol = new G4Ellipsoid(name, thisMolecule.fsize.getX(), thisMolecule.fsize.getY(),
854 thisMolecule.fsize.getZ(), -thisMolecule.fsize.getZ(), z_cut);
855 }
856
857 rot_pointer = new G4RotationMatrix(rot.inverse());
858 auto molLogic = new G4LogicalVolume(mol, fpSugarMaterial, name);
859
860 molLogic->SetVisAttributes(vattrs);
861 pv_placement = new G4PVPlacement(rot_pointer, thisMolecule.fposition, molLogic, name,
862 physicsLogical, false, 0, this->GetOverlaps());
863 }
864 else {
865 // Get information about the phsophate
866 G4ThreeVector z_direction = phosphate.fposition - thisMolecule.fposition;
867 G4double separation = z_direction.mag();
868 G4double z_cut = separation - phosphate.fsize.getX();
869
870 G4ThreeVector local_z = G4ThreeVector(0, 0, 1);
871 G4ThreeVector z_new = z_direction / z_direction.mag();
872
873 G4ThreeVector ortho = local_z.cross(z_new);
874 G4double dp = local_z.dot(z_new);
875 G4double angle = -std::atan2(ortho.mag(), dp);
876
877 G4RotationMatrix first_rot(ortho, angle);
878 rot_pointer = new G4RotationMatrix(first_rot);
879
880 auto mol = new G4Ellipsoid(name, thisMolecule.fsize.getX(), thisMolecule.fsize.getY(),
881 thisMolecule.fsize.getZ(), -thisMolecule.fsize.getZ(), z_cut);
882
883 auto molLogic = new G4LogicalVolume(mol, fpSugarMaterial, name);
884
885 molLogic->SetVisAttributes(vattrs);
886
887 pv_placement = new G4PVPlacement(rot_pointer, thisMolecule.fposition, molLogic, name,
888 physicsLogical, false, 0, this->GetOverlaps());
889 }
890 return pv_placement;
891 }
892
893 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
894
895 G4VPhysicalVolume* DNAGeometry::PlaceBase(G4LogicalVolume* physicsLogical,
896 const molecule_t& thisMolecule,
897 const molecule_t& oppBasePair,
898 const molecule_t& nextBasePair, const molecule_t& sugar)
899 {
900 std::stringstream ss;
901 ss << thisMolecule.fname << "-" << thisMolecule.fchain_id << "-" << thisMolecule.fstrand_id << "-"
902 << thisMolecule.fbase_idx;
903 G4String name = ss.str();
904 // Initial Rotation of the sphere
905 G4RotationMatrix rot;
906 rot.rotateX(thisMolecule.frotation.getX());
907 rot.rotateY(thisMolecule.frotation.getY());
908 rot.rotateZ(thisMolecule.frotation.getZ());
909
910 G4ThreeVector sugar_dirn = sugar.fposition - thisMolecule.fposition;
911 G4ThreeVector oppBPdirn = oppBasePair.fposition - thisMolecule.fposition;
912 G4ThreeVector nextBPdirn = nextBasePair.fposition - thisMolecule.fposition;
913 // Do a cut for the sugar, do a scaling for the base pair
914 G4double sugar_cut = sugar_dirn.mag() - sugar.fsize.getX();
915
916 G4double base_scaling = oppBPdirn.mag() - oppBasePair.fsize.getY();
917 base_scaling = std::min(1., base_scaling / thisMolecule.fsize.getY());
918
919 // Find the new constrained values
920 G4double thisShapeX = base_scaling * thisMolecule.fsize.getX();
921 G4double thisShapeY = base_scaling * thisMolecule.fsize.getY();
922 G4double thisShapeZ = base_scaling * thisMolecule.fsize.getZ();
923 thisShapeZ = std::min(thisShapeZ, 1.7 * angstrom);
924
925 // first rotation to turn the ellipse to run along sugar_dirn
926 G4ThreeVector z_old = G4ThreeVector(0, 0, 1);
927 G4ThreeVector z_new = sugar_dirn / sugar_dirn.mag();
928
929 G4ThreeVector ortho = z_old.cross(z_new);
930 G4double angle = std::acos(z_new.dot(z_old));
931
932 // Now find out quadrant of angle
933 G4RotationMatrix first_rot;
934 G4RotationMatrix first_rot_1(ortho, angle);
935 G4RotationMatrix first_rot_2(ortho, angle + halfpi);
936 G4RotationMatrix first_rot_3(ortho, angle + pi);
937 G4RotationMatrix first_rot_4(ortho, angle + 3 * halfpi);
938
939 if (first_rot_1(z_old).dot(z_new) > 0.99) {
940 first_rot = first_rot_1;
941 }
942 else if (first_rot_2(z_old).dot(z_new) > 0.99) {
943 first_rot = first_rot_2;
944 }
945 else if (first_rot_3(z_old).dot(z_new) > 0.99) {
946 first_rot = first_rot_3;
947 }
948 else if (first_rot_4(z_old).dot(z_new) > 0.99) {
949 first_rot = first_rot_4;
950 }
951 else {
952 G4cout << "Could not identify first rotation" << G4endl;
953 }
954
955 // second rotation
956 // Get y vector aligned with next_base_pair_direction through a rotation
957 // around sugar_dirn
958 G4ThreeVector y = first_rot(G4ThreeVector(0, 1, 0));
959 G4ThreeVector up = rot(G4ThreeVector(0, 0, 1));
960
961 G4double interval = 0.1; // precision for finding theta
962 G4double theta = -2 * interval;
963 G4double prev1 = -DBL_MAX;
964 G4double prev2 = -DBL_MAX;
965 G4double val = -DBL_MAX;
966 G4RotationMatrix z_rotation(z_new, theta);
967 while (theta <= twopi + interval) {
968 prev2 = prev1;
969 prev1 = val;
970 theta += interval;
971 z_rotation.rotate(interval, z_new);
972 val = z_rotation(y).dot(up);
973 if ((prev1 < val) && (prev1 < prev2)) {
974 // prev1 contains a local minimum
975 theta -= interval;
976 break;
977 }
978 }
979 G4RotationMatrix second_rot(first_rot.inverse()(z_new), theta);
980
981 auto rot_pointer = new G4RotationMatrix(second_rot.inverse() * first_rot.inverse());
982 auto mol = new G4Ellipsoid(name, thisShapeX, thisShapeZ, thisShapeY, -thisShapeY, sugar_cut);
983
984 G4LogicalVolume* molLogic = nullptr;
985 G4Color color;
986 if (thisMolecule.fname == "Guanine") {
987 color = G4Color::Green();
988 molLogic = new G4LogicalVolume(mol, fpGuanineMaterial, name);
989 }
990 else if (thisMolecule.fname == "Adenine") {
991 color = G4Color::Blue();
992 molLogic = new G4LogicalVolume(mol, fpAdenineMaterial, name);
993 }
994 else if (thisMolecule.fname == "Cytosine") {
995 color = G4Color::Cyan();
996 molLogic = new G4LogicalVolume(mol, fpCytosineMaterial, name);
997 }
998 else if (thisMolecule.fname == "Thymine") {
999 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 }