| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 /// file: DNAGeometry.cc
28 /// brief:
29 /*
30 * This class builds the DNA geometry. It inte
31 * with the DNAWorld class, to create a physic
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........oo
55
56 G4bool compareLVByName::operator()(const G4Log
57 {
58 return lhs->GetName() < rhs->GetName();
59 }
60
61 //....oooOO0OOooo........oooOO0OOooo........oo
62
63 DNAGeometry::DNAGeometry()
64 {
65 fpMessenger = new DNAGeometryMessenger(this)
66 fpChromosomeMapper = new ChromosomeMapper();
67
68 // Create Damage Model and add some default
69 fpDamageModel = new DamageModel();
70 G4NistManager* man = G4NistManager::Instance
71 fpSugarMaterial = man->FindOrBuildMaterial("
72 fpMediumMaterial = man->FindOrBuildMaterial(
73 fpPhosphateMaterial = man->FindOrBuildMateri
74 fpGuanineMaterial = man->FindOrBuildMaterial
75 fpAdenineMaterial = man->FindOrBuildMaterial
76 fpThymineMaterial = man->FindOrBuildMaterial
77 fpCytosineMaterial = man->FindOrBuildMateria
78 fpHistoneMaterial = man->FindOrBuildMaterial
79 fpDNAPhysicsWorld = new DNAWorld();
80 fpDrawCellAttrs = new G4VisAttributes();
81 fpDrawCellAttrs->SetColor(G4Color(0, 0, 1, 0
82 fpDrawCellAttrs->SetForceSolid(true);
83 // set default molecule sizes
84 fMoleculeSizes["phosphate"] = G4ThreeVector(
85 fMoleculeSizes["sugar"] = G4ThreeVector(2.63
86 fMoleculeSizes["guanine"] = G4ThreeVector(3.
87 fMoleculeSizes["cytosine"] = G4ThreeVector(3
88 fMoleculeSizes["adenine"] = G4ThreeVector(3.
89 fMoleculeSizes["thymine"] = G4ThreeVector(4.
90 fMoleculeSizes["histone"] = G4ThreeVector(25
91 }
92
93 //....oooOO0OOooo........oooOO0OOooo........oo
94
95 DNAGeometry::~DNAGeometry()
96 {
97 delete fpMessenger;
98 }
99
100 //....oooOO0OOooo........oooOO0OOooo........oo
101
102 void DNAGeometry::BuildDNA(G4LogicalVolume* vo
103 {
104 // TODO: Add Assertion tests to make sure fi
105 fpDNAVolumeChem = vol;
106
107 // TODO: Make a complete copy of vol and it'
108 auto* volAsEllipsoid = dynamic_cast<G4Ellips
109 auto* volAsBox = dynamic_cast<G4Box*>(vol->G
110
111 G4VSolid* cloneSolid;
112 if (volAsEllipsoid != nullptr) {
113 cloneSolid = new G4Ellipsoid(*((G4Ellipsoi
114 }
115 else if (volAsBox != nullptr) {
116 cloneSolid = new G4Box(*((G4Box*)vol->GetS
117 }
118 else {
119 G4ExceptionDescription errmsg;
120 errmsg << "An invalid physical volume shap
121 G4Exception("MolecularDnaGeometry", "ERR_I
122 return;
123 }
124
125 fpDNAVolumePhys = new G4LogicalVolume(cloneS
126 "DNAPh
127 FillParameterisedSpace();
128 fpDNAPhysicsWorld->SetDNAVolumePointer(fpDNA
129 }
130
131 //....oooOO0OOooo........oooOO0OOooo........oo
132
133 void DNAGeometry::FillParameterisedSpace()
134 {
135 // Create map for each logical volume with i
136 // LV, placement_index, copy number, positio
137 std::vector<placement> placementVector;
138 std::map<G4LogicalVolume*, G4int> currentCop
139 G4int numberOfChains = 0;
140
141 // read in and create the LV's specified in
142 std::map<G4String, G4LogicalVolume*> namesTo
143 for (const auto& it : fVoxelNames) {
144 G4String thisVoxelName = it.first;
145 G4String thisVoxelFile = it.second;
146
147 auto voxel = LoadVoxelVolume(thisVoxelName
148 G4LogicalVolume* lv = voxel.first;
149 fInfoMap[lv] = voxel.second;
150 namesToLVs[thisVoxelName] = lv;
151 currentCopyNumber[lv] = 0;
152 numberOfChains = this->GetPVInfo(lv)->GetN
153 if (!(numberOfChains == 1 || numberOfChain
154 G4cout << "Geometry contains volumes wit
155 << " chains. This can cause error
156 }
157 }
158
159 // Go through the voxel types to build a pla
160 for (unsigned ii = 0; ii < fVoxelTypes.size(
161 G4LogicalVolume* lv = namesToLVs[fVoxelTyp
162
163 placement thisPlacement = std::make_tuple(
164
165 placementVector.push_back(thisPlacement);
166 }
167
168 // Do each placement.
169 std::map<G4String, int64_t> histonesPerChrom
170 std::map<G4String, int64_t> basePairsPerChro
171 // std::map<G4String, long long> basePairsPe
172 for (auto it = placementVector.begin(); it !
173 G4LogicalVolume* thisLogical = std::get<0>
174 if (thisLogical == nullptr) {
175 G4Exception("DNAGeometry::FillParameteri
176 "At least one of the placeme
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>(*
183 G4ThreeVector thisRotation = std::get<4>(*
184 std::stringstream ss;
185 ss << thisLogical->GetName() << "-" << thi
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->inve
192 delete inv;
193
194 G4bool isPhysicallyPlaced = false;
195 VirtualChromosome* thisChromosome = this->
196 if (thisChromosome != nullptr) {
197 isPhysicallyPlaced = true;
198 // Place for Physics
199 new G4PVPlacement(pRot, thisPosition, th
200 thisCopyNo, this->GetO
201 // Place for Chemistry
202 new G4PVPlacement(pRot, thisPosition, th
203 fpDNAVolumeChem, false
204
205 // dousatsu==============
206 if (histonesPerChromosome.find(thisChrom
207 histonesPerChromosome[thisChromosome->
208 }
209 histonesPerChromosome[thisChromosome->Ge
210 this->GetPVInfo(thisLogical)->GetTotal
211 // dousatsu==============
212
213 if (basePairsPerChromosome.find(thisChro
214 basePairsPerChromosome[thisChromosome-
215 }
216 basePairsPerChromosome[thisChromosome->G
217 this->GetPVInfo(thisLogical)->GetTotal
218 }
219
220 // add the placement to the local register
221 // chromosome, as this tracks base pair in
222 if (numberOfChains == 4 || numberOfChains
223 AddFourChainPlacement(it, placementVecto
224 }
225 else if (numberOfChains == 1) {
226 AddSingleChainPlacement(it, placementVec
227 }
228 else {
229 G4ExceptionDescription errmsg;
230 errmsg << "Having none of 1/4/8 chains i
231 << "supported, the application wi
232 G4Exception("DNAGeometry::FillParameteri
233 FatalException, errmsg);
234 }
235 }
236
237 // dousatsu------------------------------
238 G4cout << "Histones placed per chromosome:"
239 G4cout << "Chromosome Histones
240 G4cout << "---------------------------------
241 for (auto& it : histonesPerChromosome) {
242 G4cout << it.first.c_str() << " : " << it.
243 }
244 // dousatsu------------------------------
245
246 G4cout << "Base Pairs placed per chromosome:
247 G4cout << "Chromosome Base Pai
248 G4cout << "---------------------------------
249 for (auto& it : basePairsPerChromosome) {
250 G4cout << it.first.c_str() << " : " << it.
251 }
252 G4cout << "-------------------------------"
253
254 if (this->GetVerbosity() > 0) {
255 this->PrintOctreeStats();
256 }
257 }
258
259 //....oooOO0OOooo........oooOO0OOooo........oo
260
261 void DNAGeometry::AddSingleChainPlacement(cons
262 cons
263 G4bo
264 {
265 G4LogicalVolume* thisLogical = std::get<0>(*
266 G4bool twist = fVoxelTwist[thisLogical->GetN
267 AddNewPlacement(thisLogical, {{0, 1, 2, 3, 4
268 }
269
270 //....oooOO0OOooo........oooOO0OOooo........oo
271
272 void DNAGeometry::AddFourChainPlacement(const
273 const
274 G4bool
275 {
276 G4LogicalVolume* thisLogical = std::get<0>(*
277 G4int thisIndex = std::get<1>(*it);
278 // Use info to build the placement transform
279 std::array<int, 8> global_chain{};
280 G4bool twist = fVoxelTwist[thisLogical->GetN
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>(pr
288 G4ThreeVector thisRotation = std::get<4>(*
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() - rotn
303 if (zdiff.mag2() > 0.1) {
304 rotold = rotold.rotate(halfpi, rotold.co
305 }
306 rotold.rectify();
307 // rotnew.colz == rotold.colz now.
308 // There exists now a rotation around rotn
309 // that transforms rotold to rotnew.
310 G4RotationMatrix relative = rotnew * rotol
311 relative.rectify();
312 G4ThreeVector axis = relative.getAxis();
313 G4double angle = relative.getDelta();
314 // get the sign of the rotation, and corre
315 G4int sign = (axis.dot(rotnew.colZ()) >= 0
316 if (sign < 0) {
317 angle = 2 * pi - angle;
318 }
319 G4int zrots = ((G4int)(2.05 * angle / pi))
320
321 global_chain = {{(this->GetGlobalChain(thi
322 (this->GetGlobalChain(thi
323 (this->GetGlobalChain(thi
324 (this->GetGlobalChain(thi
325 4 + (this->GetGlobalChain
326 4 + (this->GetGlobalChain
327 4 + (this->GetGlobalChain
328 4 + (this->GetGlobalChain
329 twist = false;
330 }
331 this->AddNewPlacement(thisLogical, global_ch
332 }
333
334 //....oooOO0OOooo........oooOO0OOooo........oo
335
336 void DNAGeometry::AddNewPlacement(const G4Logi
337 G4bool twist
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)->GetPairsOnCha
347 this->GetPVInfo(lv)->GetPairsOnCha
348 this->GetPVInfo(lv)->GetPairsOnCha
349 this->GetPVInfo(lv)->GetPairsOnCha
350 }
351 else {
352 placement_transform previous = fPlacementT
353 std::array<int64_t, 8> previousEnd = std::
354 // std::array<long long, 8> previousEnd =
355 start = {{previousEnd[0], previousEnd[1],
356 previousEnd[5], previousEnd[6],
357 end = {{previousEnd[0] + this->GetPVInfo(l
358 previousEnd[1] + this->GetPVInfo(l
359 previousEnd[2] + this->GetPVInfo(l
360 previousEnd[3] + this->GetPVInfo(l
361 previousEnd[4] + this->GetPVInfo(l
362 previousEnd[5] + this->GetPVInfo(l
363 previousEnd[6] + this->GetPVInfo(l
364 previousEnd[7] + this->GetPVInfo(l
365 }
366 pt = std::make_tuple(global_chain, start, en
367 fPlacementTransformations.push_back(pt);
368 }
369 //....oooOO0OOooo........oooOO0OOooo........oo
370
371 std::pair<G4LogicalVolume*, PlacementVolumeInf
372 DNAGeometry::LoadVoxelVolume(const G4String& v
373 {
374 if (this->GetVerbosity() > 0) {
375 G4cout << "Loading Voxel: " << volumeName
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
385 auto physicsPhysical = new G4Box(volumeName,
386 auto physicsLogical = new G4LogicalVolume(ph
387 auto chemPhysical = new G4Box(volumeName, vx
388 auto chemLogical = new G4LogicalVolume(chemP
389
390 if (this->GetDrawCellVolumes()) {
391 chemLogical->SetVisAttributes(fpDrawCellAt
392 }
393
394 fChemToPhys[chemLogical] = physicsLogical;
395 fPhysToChem[physicsLogical] = chemLogical;
396
397 physicsLogical->SetSmartless(this->GetSmartl
398
399 auto thisOctree = new OctreeNode(G4ThreeVect
400 auto thisHistoneOctree =
401 new OctreeNode(G4ThreeVector(0, 0, 0), G4T
402
403 // open and load file
404 if (!utility::Path_exists(filename)) {
405 G4ExceptionDescription errmsg;
406 errmsg << "The file: " << filename << " co
407 G4Exception("DNAGeometry::LoadVoxelVolume"
408 }
409 std::ifstream infile(filename, std::ifstream
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::
426 // validation
427 if (uncommented_line_number == 0) {
428 line_fields = line.size();
429 file_contains_size = (line_fields ==
430 }
431 if ((line_fields != 10) && (line_field
432 G4ExceptionDescription errmsg;
433 errmsg << filename << " should have
434 G4Exception("DNAGeometry::LoadVoxelV
435 errmsg);
436 }
437 if (line_fields != (G4int)line.size())
438 G4ExceptionDescription errmsg;
439 errmsg << "Unexpected number of fiel
440 << G4endl << "Expected " << l
441 << G4endl;
442 G4Exception("DNAGeometry::LoadVoxelV
443 errmsg);
444 }
445 molecule_t thisMolecule;
446
447 thisMolecule.fname = (G4String)line.at
448 thisMolecule.fshape = (G4String)line.a
449 thisMolecule.fchain_id = (G4int)std::s
450 thisMolecule.fstrand_id = (G4int)std::
451 thisMolecule.fbase_idx = (int64_t)std:
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_siz
458 if (file_contains_size) {
459 mol_size =
460 G4ThreeVector(std::stod(line.at(5)
461 std::stod(line.at(7)
462 }
463 else {
464 if (fMoleculeSizes.count(lower_name)
465 G4ExceptionDescription errmsg;
466 errmsg << "A molecule size has not
467 G4Exception("DNAGeometry::LoadVoxe
468 errmsg);
469 }
470 mol_size = fMoleculeSizes.at(lower_n
471 }
472 pos = G4ThreeVector(std::stod(line.at(
473 std::stod(line.at(
474 std::stod(line.at(
475
476 rot =
477 G4ThreeVector(std::stod(line.at(8 +
478 std::stod(line.at(10 +
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"
499 << "Filename: " << filename << G4en
500 << "# NAME SHAPE CHAIN_ID STRAND_ID
501 << "SIZE_Z POS_X POS_Y POS_Z ROT_X
502 << "# NAME SHAPE CHAIN_ID STRAND_ID
503 "ROT_X ROT_Y ROT_Z"
504 << G4endl
505 << "Separator is space, lines with
506 "commented"
507 << G4endl;
508
509 G4Exception("DNAGeometry::LoadVoxelVolume"
510 }
511
512 // We can place the elements now, specify an
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 << "
522 }
523 // Identify number of histones in voxel vo
524 pv_placement = PlaceHistone(physicsLogical
525 thisHistoneOctree->AddPhysicalVolume(pv_pl
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 << "
536 }
537 // Identify number of molecules on chain
538 if (thisInfo->GetPairsOnChain(thisMolecule
539 thisInfo->SetPairsOnChain(thisMolecule.f
540 }
541
542 // Place Volume
543 if (thisMolecule.fname == "Sugar") {
544 molecule_t phosphate = molecules[ii - 1]
545 // G4int phosph_idx = (thisMolecule.stra
546 G4int phosph_idx = ii - 1;
547 if ((phosph_idx >= 0) && (phosph_idx < (
548 phosphate = molecules[phosph_idx];
549 if (phosphate.fchain_id != thisMolecul
550 // There is no next sugar, just the
551 phosphate = molecule_t();
552 phosphate.fname = "wall";
553 }
554 }
555 else {
556 // There is no next sugar, just the vo
557 phosphate = molecule_t();
558 phosphate.fname = "wall";
559 }
560 pv_placement = PlaceSugar(physicsLogical
561 }
562 else if (thisMolecule.fname == "Phosphate"
563 molecule_t sugar;
564 G4int sugar_idx = (thisMolecule.fstrand_
565 if ((sugar_idx >= 0) && (sugar_idx < mol
566 sugar = molecules[sugar_idx];
567 if (sugar.fchain_id != thisMolecule.fc
568 // There is no next sugar, just the
569 sugar = molecule_t();
570 sugar.fname = "wall";
571 }
572 }
573 else {
574 // There is no next sugar, just the vo
575 sugar = molecule_t();
576 sugar.fname = "wall";
577 }
578 pv_placement = PlacePhosphate(physicsLog
579 }
580 else if ((thisMolecule.fname == "Guanine")
581 || (thisMolecule.fname == "Cytosi
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
595 }
596 nextBPidx = (nextBPexists && !endOfCha
597 }
598 else {
599 oppBasePair = molecules[ii - 3];
600 nextBPexists = (ii - 6 >= 0);
601 if (nextBPexists) {
602 endOfChain = (thisMolecule.fchain_id
603 }
604 nextBPidx = (nextBPexists && !endOfCha
605 }
606 nextBasePair = molecules[nextBPidx];
607 pv_placement = PlaceBase(physicsLogical,
608 }
609 else {
610 G4ExceptionDescription errmsg;
611 errmsg << "Molecule name: " << thisMolec
612 << G4endl;
613 G4Exception("DNAGeometry::LoadVoxelVolum
614 }
615 if (this->GetOverlaps()) {
616 if (pv_placement->CheckOverlaps()) {
617 // NOTE: if there are overlaps, the si
618 // crash. This is the desired behaviou
619 G4cout << "Overlap, translation: " <<
620 G4cout << "Overlap, rotationX: " << pv
621 G4cout << "Overlap, rotationY: " << pv
622 G4cout << "Overlap, rotationZ: " << pv
623 }
624 }
625 if (this->GetDrawCellVolumes()) {
626 // do not draw DNA, as we are drawing th
627 pv_placement->GetLogicalVolume()->SetVis
628 }
629 thisOctree->AddPhysicalVolume(pv_placement
630 }
631 thisInfo->SetOctree(thisOctree);
632 return std::make_pair(physicsLogical, thisIn
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 << " co
641 G4Exception("DNAGeometry::FillVoxelVectors
642 }
643 std::ifstream infile(filename, std::ifstream
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 :
658 std::vector<G4String> line = utility
659 fVoxelIndices.push_back(std::stoi(li
660 fVoxelTypes.push_back((G4String)line
661 pos = G4ThreeVector(std::stod(line.a
662 std::stod(line.a
663 std::stod(line.a
664 rot = G4ThreeVector(std::stod(line.a
665 fVoxelPositions.push_back(pos);
666 fVoxelRotations.push_back(rot * pi_o
667 }
668 catch (int exception) {
669 G4cout << "Voxel Position file read
670 G4cout << "Got Line: " << currentlin
671 G4cout << "Format should be: "
672 << "IDX KIND POS_X POS_Y POS_
673 << " EUL_PHI" << G4endl;
674 G4cout << "Separator is space, lines
675 << "character are commented"
676 std::exit(1);
677 }
678 }
679 }
680 infile.close();
681 }
682 else {
683 G4cout << "Voxel Position file read error
684 G4cout << "Format should be: "
685 << "IDX KIND POS_X POS_Y POS_Z EUL_
686 G4cout << "Separator is space, lines with
687 << "character are commented" << G4e
688 }
689 }
690
691 OctreeNode* DNAGeometry::GetTopOctreeNode(G4Lo
692 {
693 const PlacementVolumeInfo* info = GetPVInfo(
694 return (info == nullptr) ? nullptr : info->G
695 }
696
697 const PlacementVolumeInfo* DNAGeometry::GetPVI
698 {
699 if (fInfoMap.find(lv) != fInfoMap.end()) {
700 return fInfoMap.at(lv);
701 }
702 else if (fChemToPhys.find(lv) != fChemToPhys
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 M
713 G4cout << "---------------------------------
714 OctreeNode* currentOctree;
715 for (auto& it : fInfoMap) {
716 char buffer[80];
717 currentOctree = it.second->GetOctree();
718 G4int vols = currentOctree->GetContents().
719 const char* lvname = it.first->GetName().c
720 G4int mvols = currentOctree->GetMaxContent
721 G4int nodes = currentOctree->GetNumberOfTe
722 std::snprintf(buffer, 80, "%-12s %-7i %-7
723 G4cout << buffer << G4endl;
724 }
725 G4cout << "---------------------------------
726 }
727
728 //////////////////////////////////////////////
729 //////////////////////////////////////////////
730 //////////////////////////////////////////////
731
732 G4VPhysicalVolume* DNAGeometry::PlacePhosphate
733
734
735 {
736 // Define Vis Attributes
737 auto vattrs = new G4VisAttributes(G4Color::Y
738 vattrs->SetForceSolid(true);
739
740 std::stringstream ss;
741 ss << thisMolecule.fname << "-" << thisMolec
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
757 std::abs(thisMolecule
758 std::abs(thisMolecule
759 G4ThreeVector overlaps = abspos + thisMole
760
761 G4Ellipsoid* mol;
762 // NOTE: will work to cut molecules at bou
763 // leaving the box.
764 G4double z_cut = -1 * utility::Min(overlap
765 if (z_cut < 0) // all molecules fit
766 {
767 mol = new G4Ellipsoid(name, thisMolecule
768 thisMolecule.fsize
769 }
770 else // need to cut molecule to fit in bo
771 {
772 mol = new G4Ellipsoid(name, thisMolecule
773 thisMolecule.fsize
774 }
775
776 rot_pointer = new G4RotationMatrix(rot.inv
777 auto molLogic = new G4LogicalVolume(mol, f
778
779 molLogic->SetVisAttributes(vattrs);
780 pv_placement = new G4PVPlacement(rot_point
781 physicsLo
782 }
783 else {
784 // Case where there is a sugar molecule th
785 G4ThreeVector z_direction = sugar.fpositio
786 G4double separation = z_direction.mag();
787 G4double z_cut = separation - sugar.fsize.
788
789 G4ThreeVector local_z = G4ThreeVector(0, 0
790 G4ThreeVector z_new = z_direction / z_dire
791
792 G4ThreeVector ortho = local_z.cross(z_new)
793 G4double dp = local_z.dot(z_new);
794 G4double angle = -std::atan2(ortho.mag(),
795
796 G4RotationMatrix second_rot(ortho, angle);
797 rot_pointer = new G4RotationMatrix(second_
798
799 auto mol = new G4Ellipsoid(name, thisMolec
800 thisMolecule.fs
801
802 auto molLogic = new G4LogicalVolume(mol, f
803
804 molLogic->SetVisAttributes(vattrs);
805 pv_placement = new G4PVPlacement(rot_point
806 physicsLo
807 }
808
809 return pv_placement;
810 }
811
812 G4VPhysicalVolume* DNAGeometry::PlaceSugar(G4L
813 con
814 con
815 {
816 // Define Vis Attributes
817 auto vattrs = new G4VisAttributes(G4Color::R
818 vattrs->SetForceSolid(true);
819
820 // Name:
821 std::stringstream ss;
822 ss << thisMolecule.fname << "-" << thisMolec
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
838 std::abs(thisMolecule
839 std::abs(thisMolecule
840 G4ThreeVector overlaps = abspos + thisMole
841
842 G4Ellipsoid* mol;
843 // NOTE: will work to cut molecules at bou
844 // leaving the box.
845 G4double z_cut = -1 * utility::Min(overlap
846 if (z_cut < 0) // all molecules fit
847 {
848 mol = new G4Ellipsoid(name, thisMolecule
849 thisMolecule.fsize
850 }
851 else // need to cut molecule to fit in bo
852 {
853 mol = new G4Ellipsoid(name, thisMolecule
854 thisMolecule.fsize
855 }
856
857 rot_pointer = new G4RotationMatrix(rot.inv
858 auto molLogic = new G4LogicalVolume(mol, f
859
860 molLogic->SetVisAttributes(vattrs);
861 pv_placement = new G4PVPlacement(rot_point
862 physicsLo
863 }
864 else {
865 // Get information about the phsophate
866 G4ThreeVector z_direction = phosphate.fpos
867 G4double separation = z_direction.mag();
868 G4double z_cut = separation - phosphate.fs
869
870 G4ThreeVector local_z = G4ThreeVector(0, 0
871 G4ThreeVector z_new = z_direction / z_dire
872
873 G4ThreeVector ortho = local_z.cross(z_new)
874 G4double dp = local_z.dot(z_new);
875 G4double angle = -std::atan2(ortho.mag(),
876
877 G4RotationMatrix first_rot(ortho, angle);
878 rot_pointer = new G4RotationMatrix(first_r
879
880 auto mol = new G4Ellipsoid(name, thisMolec
881 thisMolecule.fs
882
883 auto molLogic = new G4LogicalVolume(mol, f
884
885 molLogic->SetVisAttributes(vattrs);
886
887 pv_placement = new G4PVPlacement(rot_point
888 physicsLo
889 }
890 return pv_placement;
891 }
892
893 //....oooOO0OOooo........oooOO0OOooo........oo
894
895 G4VPhysicalVolume* DNAGeometry::PlaceBase(G4Lo
896 cons
897 cons
898 cons
899 {
900 std::stringstream ss;
901 ss << thisMolecule.fname << "-" << thisMolec
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 -
911 G4ThreeVector oppBPdirn = oppBasePair.fposit
912 G4ThreeVector nextBPdirn = nextBasePair.fpos
913 // Do a cut for the sugar, do a scaling for
914 G4double sugar_cut = sugar_dirn.mag() - suga
915
916 G4double base_scaling = oppBPdirn.mag() - op
917 base_scaling = std::min(1., base_scaling / t
918
919 // Find the new constrained values
920 G4double thisShapeX = base_scaling * thisMol
921 G4double thisShapeY = base_scaling * thisMol
922 G4double thisShapeZ = base_scaling * thisMol
923 thisShapeZ = std::min(thisShapeZ, 1.7 * angs
924
925 // first rotation to turn the ellipse to run
926 G4ThreeVector z_old = G4ThreeVector(0, 0, 1)
927 G4ThreeVector z_new = sugar_dirn / sugar_dir
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 +
936 G4RotationMatrix first_rot_3(ortho, angle +
937 G4RotationMatrix first_rot_4(ortho, angle +
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.9
943 first_rot = first_rot_2;
944 }
945 else if (first_rot_3(z_old).dot(z_new) > 0.9
946 first_rot = first_rot_3;
947 }
948 else if (first_rot_4(z_old).dot(z_new) > 0.9
949 first_rot = first_rot_4;
950 }
951 else {
952 G4cout << "Could not identify first rotati
953 }
954
955 // second rotation
956 // Get y vector aligned with next_base_pair_
957 // around sugar_dirn
958 G4ThreeVector y = first_rot(G4ThreeVector(0,
959 G4ThreeVector up = rot(G4ThreeVector(0, 0, 1
960
961 G4double interval = 0.1; // precision for f
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.invers
980
981 auto rot_pointer = new G4RotationMatrix(seco
982 auto mol = new G4Ellipsoid(name, thisShapeX,
983
984 G4LogicalVolume* molLogic = nullptr;
985 G4Color color;
986 if (thisMolecule.fname == "Guanine") {
987 color = G4Color::Green();
988 molLogic = new G4LogicalVolume(mol, fpGuan
989 }
990 else if (thisMolecule.fname == "Adenine") {
991 color = G4Color::Blue();
992 molLogic = new G4LogicalVolume(mol, fpAden
993 }
994 else if (thisMolecule.fname == "Cytosine") {
995 color = G4Color::Cyan();
996 molLogic = new G4LogicalVolume(mol, fpCyto
997 }
998 else if (thisMolecule.fname == "Thymine") {
999 color = G4Color::Magenta();
1000 molLogic = new G4LogicalVolume(mol, fpThy
1001 }
1002 else {
1003 G4ExceptionDescription errmsg;
1004 errmsg << "Molecule name: " << thisMolecu
1005 << G4endl;
1006 G4Exception("DNAGeometry::PlaceBase", "In
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, thisMolecu
1016 this->GetOverlaps());
1017
1018 return pv_placement;
1019 }
1020
1021 //....oooOO0OOooo........oooOO0OOooo........o
1022
1023 G4VPhysicalVolume* DNAGeometry::PlaceHistone(
1024
1025 {
1026 G4VPhysicalVolume* pv_placement;
1027
1028 // Define Vis Attributes
1029 auto vattrs = new G4VisAttributes(G4Color::
1030 vattrs->SetColor(G4Colour(0, 0, 1, 0.2));
1031 vattrs->SetForceSolid(true);
1032
1033 // Name:
1034 std::stringstream ss;
1035 ss << thisMolecule.fname << "-" << thisMole
1036 << thisMolecule.fbase_idx;
1037 G4String name = ss.str();
1038
1039 //// Does the Histone fit in the voxel (ass
1040 // G4double max_coord = std::max(std::abs(t
1041 // std::abs(th
1042 // max_coord = std::max(max_coord, std::abs
1043 //// if overlap
1044 // G4double radius = utility::max(thisMolec
1045 // G4double radiusx=0,radiusx=0,radiusx=0;
1046 // if ((max_coord + radius) > fVoxelSideLen
1047 //{
1048 // radius = 0.9*(fVoxelSideLength - max_
1049 // G4cout<<" Warning : Histone overlap w
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, thisMolecu
1060 thisMolecule.fsi
1061
1062 fHistoneInteractionRadius = thisMolecule.fs
1063
1064 rot_pointer = new G4RotationMatrix(rot.inve
1065
1066 auto molLogic = new G4LogicalVolume(mol, fp
1067 molLogic->SetVisAttributes(vattrs);
1068
1069 pv_placement = new G4PVPlacement(rot_pointe
1070 physicsLog
1071
1072 return pv_placement;
1073 }
1074 // dousatsu ---------------------------------
1075
1076 //....oooOO0OOooo........oooOO0OOooo........o
1077
1078 void DNAGeometry::FindNearbyMolecules(const G
1079 std::ve
1080 double
1081 {
1082 const PlacementVolumeInfo* pvInfo = GetPVIn
1083 if (pvInfo == nullptr) {
1084 return;
1085 }
1086 // never search more than the radical kill
1087 searchRange = std::min(searchRange, fRadica
1088 OctreeNode* octreeNode = pvInfo->GetOctree(
1089 octreeNode->SearchOctree(localPosition, dau
1090
1091 // G4double r = std::pow(std::pow(local
1092 // std::pow(localPosition.getY(), 2), 0
1093 // {
1094 // G4cout << "candidates found in rad
1095 // << daughters_pv_out.size()
1096 // G4cout << "Radius from center of p
1097 // << r << " nm" << G4endl;
1098 // G4cout << "-----------------------
1099 // G4endl;
1100 // }
1101 }
1102
1103 //....oooOO0OOooo........oooOO0OOooo........o
1104
1105 G4bool DNAGeometry::IsInsideHistone(const G4L
1106 const G4T
1107 {
1108 const PlacementVolumeInfo* pvInfo = GetPVIn
1109 if (pvInfo == nullptr) {
1110 G4cout << " Warning : no PV for IsInsideH
1111 return false;
1112 }
1113
1114 if (!fUseHistoneScav) {
1115 return false;
1116 }
1117
1118 OctreeNode* octreeNode = pvInfo->GetHistone
1119 const std::vector<G4VPhysicalVolume*> resul
1120 octreeNode->SearchOctree(localPosition, f
1121 return !result.empty();
1122 }
1123
1124 //....oooOO0OOooo........oooOO0OOooo........o
1125
1126 void DNAGeometry::BasePairIndexTest()
1127 {
1128 G4cout << "------------------------" << G4e
1129 G4cout << "Begin BasePairIndexTest" << G4en
1130 G4bool pass = true;
1131 for (G4int ii = 0; ii != GetNumberOfChains(
1132 G4cout << "Testing Chain " << ii << G4end
1133 std::vector<int64_t> test_vector; // dou
1134 // std::vector<long long> test_vector;//O
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)fPlacemen
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 !=
1146 if (it == test_vector.begin()) // exce
1147 {
1148 if (it + 1 != test_vector.end()) {
1149 if (*it > *(it + 1)) {
1150 G4ExceptionDescription errmsg;
1151 errmsg << "BasePairIndexTest fail
1152 << "The index " << *(it +
1153 G4Exception("DNAGeometry::BasePai
1154 }
1155 }
1156 }
1157 else if (it == test_vector.end() - 1)
1158 {
1159 if (*it < *(it - 1)) {
1160 G4ExceptionDescription errmsg;
1161 errmsg << "BasePairIndexTest fails.
1162 << "The index " << *(it - 1)
1163 G4Exception("DNAGeometry::BasePairI
1164 }
1165 }
1166 else // default case
1167 {
1168 if (*it < *(it - 1)) {
1169 G4ExceptionDescription errmsg;
1170 errmsg << "BasePairIndexTest fails.
1171 << "The index " << *(it - 1)
1172 G4Exception("DNAGeometry::BasePairI
1173 }
1174 if (*it > *(it + 1)) {
1175 G4ExceptionDescription errmsg;
1176 errmsg << "BasePairIndexTest fails.
1177 << "The index " << *(it + 1)
1178 G4Exception("DNAGeometry::BasePairI
1179 }
1180 }
1181 }
1182 }
1183 if (pass) {
1184 G4cout << "BasePairIndexTest passed" << G
1185 G4cout << "------------------------" << G
1186 }
1187 }
1188
1189 /**
1190 * Global unique id is an integer with the fo
1191 * moleculeID: base ::molecule::Count
1192 * PlacementIdx: base this->GetNumberOfPlacem
1193 * StrandID: base 2
1194 * ChainID: base GetNumberOfChains()
1195 * BasePairIdx: base GetMaxBPIdx()
1196 */
1197 //....oooOO0OOooo........oooOO0OOooo........o
1198
1199 int64_t DNAGeometry::GetGlobalUniqueID(G4VPhy
1200 {
1201 const G4String& dnaName = dnavol->GetName()
1202 std::array<G4String, 4> pv_arr = utility::G
1203 molecule mol = utility::GetMoleculeEnum(pv_
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));
1207 // G4int baseIdx = std::stoi(pv_arr.at(3)
1208
1209 const G4String& placementName = touch->GetV
1210 G4int placeIdx = std::stoi(utility::Get_sep
1211
1212 G4int chains = this->GetNumberOfChains();
1213 G4int placements = this->GetNumberOfPlaceme
1214 G4int molecules = ::molecule::Count; // do
1215 // long long chains = this->GetNumberOfChai
1216 // long long placements = this->GetNumberOf
1217 // long long molecules = ::molecule::Count
1218
1219 chainIdx = this->GetGlobalChain(placeIdx, c
1220 // long long globalPair = baseIdx + this->G
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 * str
1228 int64_t val4 = molecules * placements * 2 *
1229 int64_t val5 = molecules * placements * 2 *
1230 int64_t sum = val1 + val2 + val3 + val4 + v
1231 if (sum < 0) {
1232 G4ExceptionDescription errmsg;
1233 errmsg << "v1: " << val1 << G4endl << "v2
1234 << "v4: " << val4 << G4endl << "v5
1235 << "placeIdx: " << placeIdx << G4e
1236 << "strandIdx: " << strandIdx << G
1237 << "chains: " << chains << G4endl
1238 << G4endl;
1239 G4Exception("DNAGeometry::GetGlobalUnique
1240 }
1241 return sum;
1242 }
1243
1244 //....oooOO0OOooo........oooOO0OOooo........o
1245
1246 // copy of the GetGlobalUniqueID except with
1247 int64_t DNAGeometry::GetGlobalUniqueIDTest(in
1248 in
1249 {
1250 int64_t chains = this->GetNumberOfChains();
1251 int64_t placements = this->GetNumberOfPlace
1252 int64_t molecules = ::molecule::Count; //
1253 // long long chains = this->GetNumberOfChai
1254 // long long placements = this->GetNumberOf
1255 // long long molecules = ::molecule::Count
1256
1257 chainIdx = this->GetGlobalChain(placeIdx, c
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 * str
1265 int64_t val4 = molecules * placements * 2 *
1266 int64_t val5 = molecules * placements * 2 *
1267 int64_t sum = val1 + val2 + val3 + val4 + v
1268 if (sum < 0) {
1269 G4ExceptionDescription errmsg;
1270 errmsg << "v1: " << val1 << G4endl << "v2
1271 << "v4: " << val4 << G4endl << "v5
1272 << " molecules : " << molecules <<
1273 << "placements: " << placements <<
1274 << "chainIdx: " << chainIdx << G4e
1275 << "baseIdx: " << baseIdx << G4end
1276 G4Exception("DNAGeometry::GetGlobalUnique
1277 "ERR_BAD_ID", FatalException,
1278 }
1279 return sum;
1280 }
1281
1282 void DNAGeometry::UniqueIDTest()
1283 {
1284 G4cout << "Unique ID Test starting." << G4e
1285 G4int molecules = ::molecule::Count;
1286
1287 G4int chains = this->GetNumberOfChains();
1288 G4int placements = this->GetNumberOfPlaceme
1289 G4int basepairs = 100000; // Maximum numbe
1290
1291 G4bool pass = true;
1292 G4int counter = 0;
1293 for (G4int pment = 0; pment != (placements)
1294 for (G4int chain = 0; chain != (chains);
1295 for (G4int mol = 0; mol != (molecules);
1296 // repeat 10 times along the strand/c
1297 for (int64_t ii = 0; ii != 10; ++ii)
1298 counter++;
1299 G4int strand = G4UniformRand() > 0.
1300 G4int basepair = (G4int)basepairs *
1301 int64_t unique_id = // dousatsu
1302 // long long u
1303 this->GetGlobalUniqueIDTest(mol,
1304
1305 G4int real_strand = (this->GetStran
1306 (strand + 1)
1307
1308 G4int real_chain = this->GetGlobalC
1309 int64_t real_bp = basepair + this->
1310 // long long real_bp = basepair + t
1311 // real_chain);//ORG
1312
1313 G4int obs_mol = this->GetMoleculeFr
1314 G4int obs_placement = this->GetPlac
1315 G4int obs_strand = this->GetStrandI
1316 G4int obs_chain = this->GetChainIDF
1317 int64_t obs_basepair = this->GetBas
1318 // long long obs_basepair = this->G
1319 // (unique_id);//ORG
1320
1321 if ((pment != obs_placement) || (re
1322 || (obs_mol != mol) || (real_bp
1323 {
1324 pass = false;
1325 G4cerr << "Unique ID test failed
1326 G4cerr << "real_placement: " << p
1327 G4cerr << "real_strand: " << real
1328 G4cerr << "real_chain: " << real_
1329 G4cerr << "real_mol: " << mol <<
1330 G4cerr << "real_bp: " << real_bp
1331 }
1332 }
1333 }
1334 }
1335 }
1336 if (!pass) {
1337 G4Exception("DNAGeometry::UniqueIDTest",
1338 "Algorithm to generate unique
1339 }
1340 G4cout << "Unique ID Test completed, after
1341 }
1342
1343 //....oooOO0OOooo........oooOO0OOooo........o
1344
1345 G4Material* DNAGeometry::GetMaterialFromUniqu
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::GetMaterialFrom
1363 << "Encountered an unspecified m
1364 return fpMediumMaterial;
1365 default:
1366 G4cout << "DNAGeometry::GetMaterialFrom
1367 << "Encountered a bad material e
1368 return fpMediumMaterial;
1369 }
1370 }
1371
1372 //....oooOO0OOooo........oooOO0OOooo........o
1373
1374 G4int DNAGeometry::GetPlacementIndexFromUniqu
1375 {
1376 // remove molecule
1377 idx -= idx % (int64_t)::molecule::Count;
1378 // remove everything above molecule
1379 idx = idx % (int64_t)(::molecule::Count * t
1380 return (G4int)idx / (::molecule::Count);
1381 }
1382
1383 //....oooOO0OOooo........oooOO0OOooo........o
1384
1385 G4int DNAGeometry::GetStrandIDFromUniqueID(in
1386 {
1387 // remove molecule
1388 idx -= (int64_t)(idx % ::molecule::Count);
1389 // remove placement
1390 idx -= (int64_t)(idx % (::molecule::Count *
1391 // remove everything above strand
1392 idx = (int64_t)(idx % (::molecule::Count *
1393 // recover strand
1394 return (G4int)idx / (::molecule::Count * th
1395 }
1396
1397 //....oooOO0OOooo........oooOO0OOooo........o
1398
1399 G4int DNAGeometry::GetChainIDFromUniqueID(int
1400 {
1401 // remove molecule
1402 idx -= idx % (int64_t)::molecule::Count;
1403 // remove placement
1404 idx -= idx % (int64_t)(::molecule::Count *
1405 // remove strand
1406 idx -= idx % (int64_t)(::molecule::Count *
1407 // remove everything above chain
1408 idx =
1409 idx
1410 % (int64_t)(::molecule::Count * this->Get
1411 // recover chain
1412 G4int chain = (G4int)idx / (int64_t)(::mole
1413 return chain;
1414 }
1415
1416 //....oooOO0OOooo........oooOO0OOooo........o
1417
1418 int64_t DNAGeometry::GetBasePairFromUniqueID(
1419 {
1420 // remove molecule
1421 idx -= idx % ::molecule::Count;
1422
1423 // remove placement
1424 int64_t placement_ = idx % (int64_t)(::mole
1425 idx -= placement_;
1426
1427 // find placement for later
1428 placement_ = placement_ / (int64_t)(::molec
1429
1430 // remove strand
1431 idx -= idx % (int64_t)(::molecule::Count *
1432
1433 // remove chain
1434 int64_t chain =
1435 idx
1436 % (int64_t)(::molecule::Count * this->Get
1437 idx -= chain;
1438
1439 // find chain for later
1440 chain = chain / (int64_t)(::molecule::Count
1441
1442 // only base pairs left
1443 int64_t bp =
1444 idx
1445 / (int64_t)(::molecule::Count * this->Get
1446
1447 // use chain and placement to get BP offset
1448 bp = bp + this->GetStartIdx(placement_, cha
1449 return bp;
1450 }
1451
1452 //....oooOO0OOooo........oooOO0OOooo........o
1453
1454 int64_t DNAGeometry::GetMaxBPIdx() const //
1455 {
1456 int64_t mx = 0; // dousatsu
1457 for (auto val : std::get<2>(*fPlacementTran
1458 if (val > mx) {
1459 mx = val;
1460 }
1461 }
1462 return mx;
1463 }
1464
1465 //....oooOO0OOooo........oooOO0OOooo........o
1466
1467 G4int DNAGeometry::GetLocalChain(G4int vol_id
1468 {
1469 auto chains = std::get<0>(fPlacementTransfo
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........o
1481
1482 G4LogicalVolume* DNAGeometry::GetMatchingPhys
1483 {
1484 auto it = fChemToPhys.find(chem);
1485 if (it != fChemToPhys.end()) {
1486 return it->second;
1487 }
1488 else if (fPhysToChem.find(chem) != fPhysToC
1489 G4cout << "A phyics volume was used to fi
1490 return chem;
1491 }
1492 else {
1493 return nullptr;
1494 }
1495 }
1496
1497 //....oooOO0OOooo........oooOO0OOooo........o
1498
1499 G4LogicalVolume* DNAGeometry::GetMatchingChem
1500 {
1501 auto it = fPhysToChem.find(physics);
1502 if (it != fPhysToChem.end()) {
1503 return it->second;
1504 }
1505 else if (fChemToPhys.find(physics) != fChem
1506 G4cout << "A chem volume was used to find
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........o
1520
1521 void DNAGeometry::AddChangeMoleculeSize(const
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
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 de
1537 << " this molecule's size from"
1538 << size << G4endl;
1539 fMoleculeSizes[lower_name] = size;
1540 }
1541 }
1542 }