| 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 PhysGeoImport.cc
28 /// \brief Implementation of the PhysGeoImport
29
30 #include "PhysGeoImport.hh"
31 #include "Randomize.hh"
32 #include "G4Sphere.hh"
33 #include "G4PhysicalConstants.hh"
34
35 //....oooOO0OOooo........oooOO0OOooo........oo
36
37 PhysGeoImport::PhysGeoImport() :
38 fFactor(1.),
39 fGeoName("")
40 {
41 DefineMaterial();
42 }
43
44 //....oooOO0OOooo........oooOO0OOooo........oo
45
46 PhysGeoImport::PhysGeoImport(G4bool isVisu) :
47 fIsVisu(isVisu),
48 fFactor(1.),
49 fGeoName("")
50 {
51 if (fIsVisu) G4cout<<" For Checking Visual
52 DefineMaterial();
53 }
54
55 //....oooOO0OOooo........oooOO0OOooo........oo
56
57 G4LogicalVolume* PhysGeoImport::CreateLogicVol
58 {
59 // The idea is to return a pointer to the
60 // Create a temporary material map to assi
61
62 std::map<G4String, G4Material*> materials;
63
64 // Loop on each material to build the map
65 for(G4int i=0, ie=fMaterialVect.size(); i<
66 {
67 G4String matName("");
68
69 if(fMaterialVect[i]->GetName()=="adeni
70 else if(fMaterialVect[i]->GetName()=="
71 else if(fMaterialVect[i]->GetName()=="
72 else if(fMaterialVect[i]->GetName()=="
73 else if(fMaterialVect[i]->GetName()=="
74 else if(fMaterialVect[i]->GetName()=="
75 else if(fMaterialVect[i]->GetName()=="
76 else if(fMaterialVect[i]->GetName()=="
77 else if(fMaterialVect[i]->GetName()=="
78 else
79 {
80 matName = fMaterialVect[i]->GetNam
81 }
82
83 materials[matName] = fMaterialVect[i];
84 }
85
86
87 // Parse the input file
88
89 voxelName = ParseFile(fileName);
90
91 // Create the volumes
92
93 G4String boxNameSolid = fGeoName+"_solid";
94 G4Box* box_solid = new G4Box(boxNameSolid,
95
96 G4String boxNameLogic = fGeoName+"_logic";
97 G4LogicalVolume* box_logic = new G4Logical
98
99 // sort the molecules
100 std::sort(fMolecules.begin(), fMolecules.e
101
102 // Loop on all the parsed molecules
103 for(G4int i=0, ie=fMolecules.size(); i<ie;
104 {
105 // Retrieve general molecule informati
106 //
107 G4String name = fMolecules[i].fName;
108 G4String materialName = "water";// fMo
109 G4double radius = fMolecules[i].fRadiu
110 G4double waterRadius = fMolecules[i].f
111 G4ThreeVector moleculePosition = fMole
112 G4int copyNum = fMolecules[i].fCopyNum
113
114 // Water hydration shell volume part
115
116 G4Orb* moleculeWater_solid = 0;
117 G4VSolid* moleculeWaterCut_solid = 0;
118 G4LogicalVolume* moleculeWater_logic =
119
120 // If water radius != 0 then we have a
121 G4double tol = 0.0001;
122 if(waterRadius > (0 + tol)*nm)
123 {
124 G4Material* mat = 0;
125
126 mat=materials[materialName];
127
128 G4String nameWaterSolid = name+"_w
129 moleculeWater_solid = new G4Orb(na
130 moleculeWaterCut_solid = CreateCut
131
132 G4String nameWaterLogic = name+"_w
133 moleculeWater_logic = new G4Logica
134
135 G4String nameWaterPhys = name+"_wa
136 new G4PVPlacement(0, moleculePosit
137 }
138
139 // Dna volume part
140
141 G4Orb* molecule_solid = 0;
142 G4VSolid* moleculeCut_solid = 0;
143 G4LogicalVolume* molecule_logic = 0;
144
145 G4String nameSolid = fMolecules[i].fNa
146 molecule_solid = new G4Orb(nameSolid,
147 moleculeCut_solid = CreateCutSolid(mol
148
149 G4String nameLogic = name+"_logic";
150 molecule_logic = new G4LogicalVolume(m
151
152 // If there was a water hydration shel
153 // placed within it and, thus, its rel
154 if(waterRadius > (0 + tol)*nm)
155 {
156 G4ThreeVector position(0.,0.,0.);
157 G4String namePhys = name+"_phys";
158 new G4PVPlacement(0, position, mol
159 }
160 // If not, coordinates are those of th
161 else
162 {
163 G4ThreeVector position = moleculeP
164 G4String namePhys = name+"_phys";
165 new G4PVPlacement(0, position, mol
166 }
167 }
168
169 // Clear the containers
170 fMolecules.clear();
171 fRadiusMap.clear();
172 fWaterRadiusMap.clear();
173
174 return box_logic;
175 }
176
177 //....oooOO0OOooo........oooOO0OOooo........oo
178
179 G4String PhysGeoImport::ParseFile(const G4Stri
180 {
181 G4cout<<"Start parsing of "<<fileName<<G4e
182
183 // Clear the containers
184 fMolecules.clear();
185 fRadiusMap.clear();
186 fWaterRadiusMap.clear();
187
188 // Setup the input stream
189 std::ifstream file(fileName.c_str());
190
191 // Check if the file was correctly opened
192 if(!file.is_open())
193 {
194 // Geant4 exception
195 G4String msg = fileName+" could not be
196 G4Exception("PhysGeoImport::ParseFile(
197 }
198
199 // Define the line string variable
200 G4String line;
201 G4int preBpCpNo = -1;
202 G4int preHistoneCpNo =-1;
203 G4int NoBp = 0;
204 G4int NoHistone = 0;
205 // Read the file line per line
206 while(std::getline(file, line) )
207 {
208 // Check the line to determine if it i
209 if(line.empty()) continue; // skip the
210
211 // Data string stream
212 std::istringstream issLine(line);
213
214 // String to determine the first lette
215 G4String firstItem;
216
217 // Put the first letter/word within th
218 issLine >> firstItem;
219
220 // Check first letter to determine if
221 if(firstItem=="#") continue; // skip t
222
223 // Use the file
224 else if(firstItem=="_Name")
225 {
226 G4String name;
227 issLine >> name;
228
229 fGeoName = name;
230 }
231 else if(firstItem=="_Size")
232 {
233 G4double size;
234 issLine >> size;
235 size *= fFactor*nm;
236
237 fSize = size;
238 }
239 else if(firstItem == "_Version")
240 {
241
242 }
243 else if(firstItem=="_Number")
244 {
245
246 }
247 else if(firstItem=="_Radius")
248 {
249 G4String name;
250 issLine >> name;
251
252 G4double radius;
253 issLine >> radius;
254 radius *= fFactor*nm;
255
256 G4double waterRadius;
257 issLine >> waterRadius;
258 waterRadius *= fFactor*nm;
259
260 fRadiusMap[name] = radius;
261 fWaterRadiusMap[name] = waterRadiu
262 }
263 else if(firstItem=="_pl")
264 {
265 G4String name;
266 issLine >> name;
267
268 G4String material;
269 issLine >> material;
270
271 G4int strand;
272 issLine >> strand;
273
274 G4int copyNumber;
275 issLine >> copyNumber;
276
277 if (name == "histone") {
278 if (copyNumber != preHistoneCp
279 NoHistone ++;
280 preHistoneCpNo = copyNumbe
281 }
282 } else {
283 if (copyNumber != preBpCpNo) {
284 NoBp++;
285 preBpCpNo = copyNumber;
286 }
287 }
288
289 G4double x;
290 issLine >> x;
291 x *= fFactor*nm;
292
293 G4double y;
294 issLine >> y;
295 y *= fFactor*nm;
296
297 G4double z;
298 issLine >> z;
299 z *= fFactor*nm;
300
301 Molecule molecule(name, copyNumber
302
303 fMolecules.push_back(molecule);
304 }
305 else
306 {
307 // Geant4 exception
308 G4String msg = firstItem+" is not
309 G4Exception("PhysGeoImport::ParseF
310 }
311 }
312
313 // Close the file once the reading is done
314 file.close();
315
316 G4cout<<"End parsing of "<<fileName<<G4end
317 fVoxelNbHistoneMap.insert({fGeoName,NoHist
318 fVoxelNbBpMap.insert({fGeoName,NoBp});
319 return fGeoName;
320 }
321
322 //....oooOO0OOooo........oooOO0OOooo........oo
323
324 G4VSolid* PhysGeoImport::CreateCutSolid(G4Orb
325
326
327
328 {
329 // The idea behing this method is to cut o
330 // If a reference and a target volumes are
331 // The reference is already selected when
332
333 // Use the tiny space to differentiate the
334 G4double tinySpace = 0.001*fFactor*nm;
335
336 // Cutted solid to be returned
337 G4SubtractionSolid* solidCut(NULL);
338
339 // Some flags
340 G4bool isCutted = false;
341 G4bool isOurVol = false;
342
343 // Radius of the molecule to cut
344 G4double radiusRef;
345 if(molRef.fRadiusWater==0) radiusRef = mol
346 else radiusRef = molRef.fRadiusWater;
347
348 // Reference volume position
349 G4ThreeVector posRef = molRef.fPosition;
350 // Before looping on all the volumes we ch
351 //reference volume overlaps with its conta
352 if(std::abs(posRef.x() ) + radiusRef > fSi
353 || std::abs(posRef.y() ) + radiusR
354 || std::abs(posRef.z() ) + radiusR
355 {
356 // If we enter here, then the referenc
357 // Box used to cut
358 G4Box* solidBox = new G4Box("solid_box
359 G4ThreeVector posBox;
360
361 // Create a dummy rotation matrix
362 G4RotationMatrix *rotMat = new G4Rotat
363 rotMat->rotate(0, G4ThreeVector(0,0,1)
364
365 // To choose the cut direction
366
367 // Up
368 if(std::abs( posRef.y() + radiusRef )
369 {
370 posBox = -posRef + G4ThreeVector(0
371
372 // If the volume is cutted for the
373 if(!isCutted)
374 {
375 solidCut = new G4SubtractionSo
376 isCutted = true;
377 }
378 // For the other times
379 else solidCut = new G4SubtractionS
380 }
381
382 // Down
383 if(std::abs( posRef.y() - radiusRef )
384 {
385 posBox = -posRef + G4ThreeVector(0
386
387 // If the volume is cutted for the
388 if(!isCutted)
389 {
390 solidCut = new G4SubtractionSo
391 isCutted = true;
392 }
393 // For the other times
394 else solidCut = new G4SubtractionS
395 }
396
397 // Left
398 if(std::abs( posRef.x() + radiusRef )
399 {
400 posBox = -posRef + G4ThreeVector(f
401
402 // If the volume is cutted for the
403 if(!isCutted)
404 {
405 solidCut = new G4SubtractionSo
406 isCutted = true;
407 }
408 // For the other times
409 else solidCut = new G4SubtractionS
410 }
411
412 // Right
413 if(std::abs( posRef.x() - radiusRef )
414 {
415 posBox = -posRef + G4ThreeVector(-
416
417 // If the volume is cutted for the
418 if(!isCutted)
419 {
420 solidCut = new G4SubtractionSo
421 isCutted = true;
422 }
423 // For the other times
424 else solidCut = new G4SubtractionS
425 }
426
427 // Forward
428 if(std::abs( posRef.z() + radiusRef )
429 {
430 posBox = -posRef + G4ThreeVector(0
431
432 // If the volume is cutted for the
433 if(!isCutted)
434 {
435 solidCut = new G4SubtractionSo
436 isCutted = true;
437 }
438 // For the other times
439 else solidCut = new G4SubtractionS
440 }
441
442 // Backward
443 if(std::abs( posRef.z() - radiusRef )
444 {
445 posBox = -posRef + G4ThreeVector(0
446
447 // If the volume is cutted for the
448 if(!isCutted)
449 {
450 solidCut = new G4SubtractionSo
451 isCutted = true;
452 }
453 // For the other times
454 else solidCut = new G4SubtractionS
455 }
456 }
457
458 // Look the other volumes of the voxel
459 // Loop on all the target volumes (other v
460 for(G4int i=0, ie=molList.size(); i<ie; ++
461 {
462 G4ThreeVector posTar = molList[i].fPos
463
464 G4double rTar = posRef.z();
465 G4double zTar = posTar.z();
466
467 if(zTar>rTar+20*fFactor*nm)
468 {
469 break;
470 }
471 else if(zTar<rTar-20*fFactor*nm)
472 {
473 continue;
474 }
475
476 // Retrieve current target sphere info
477 G4double radiusTar;
478 if(molList[i].fRadiusWater==0) radiusT
479 else radiusTar = molList[i].fRadiusWat
480
481 // Compute the distance reference-targ
482 G4double distance = std::abs( (posRef
483
484 // Use the distance to check if the cu
485 // This can only happen once per loop.
486 if(distance==0 && !isOurVol)
487 {
488 // Target volume is also reference
489
490 // Set the flag
491 isOurVol = true;
492
493 // Next iteration
494 continue;
495 }
496 // If the condition is correct more th
497 else if(distance == 0 && isOurVol)
498 {
499 G4ExceptionDescription desmsg;
500 desmsg << "DetectorConstruction::C
501 G4Exception("ChemGeoImport::BuildG
502 }
503
504 // If the volumes are differents then
505 // close enough to overlap and, thus,
506 else if(distance <= radiusRef+radiusTa
507 {
508 // Volumes are close enough, there
509
510 // Box used to cut
511 G4Box* solidBox = new G4Box("solid
512
513 // This part is tricky.
514 // The goal is to calculate the po
515
516 // diff vector to from ref to tar
517 G4ThreeVector diff = posTar - posR
518
519 // Find the intersection point and
520 G4double d = (std::pow(radiusRef,2
521 (2*distance) + solidBo
522
523 // If we are in another volume we
524 //differentiate without ambiguitie
525 // (may not be necessary)
526 if(in) d -= 2*tinySpace;
527
528 // Position of the box in order to
529 // "* ( diff/diff.getR() )" is nec
530 G4ThreeVector pos = d *( diff/diff
531
532 // Get the rotation angles because
533 // to give the right "cut angle".
534 G4double phi = std::acos(pos.getZ(
535 G4double theta = std::acos( pos.ge
536
537 if(pos.getY()<0) theta = -theta;
538
539 G4ThreeVector rotAxisForPhi(1*fFac
540 rotAxisForPhi.rotateZ(theta+pi/2);
541
542 // Create the rotation matrix
543 G4RotationMatrix *rotMat = new G4R
544 rotMat->rotate(-phi, rotAxisForPhi
545
546 // Rotate it again
547 G4ThreeVector rotZAxis(0.,0.,1*fFa
548 rotMat->rotate(theta, rotZAxis);
549
550 // If the volume is cutted for the
551 if(!isCutted) solidCut = new G4Sub
552 solidBox, rotMat, pos);
553
554 // For the other times
555 else solidCut = new G4SubtractionS
556
557 // Set the cut flag
558 isCutted = true;
559 }
560 }
561
562 // If there was at least one cut then we r
563 if(isCutted) return solidCut;
564
565 // Otherwise, we return the original volum
566 else return solidOrbRef;
567
568 }
569
570 //....oooOO0OOooo........oooOO0OOooo........oo
571
572 void PhysGeoImport::DefineMaterial()
573 {
574 G4NistManager * man = G4NistManager::Insta
575 fpWater = man->FindOrBuildMaterial("G4_WAT
576 fMaterialVect.push_back(fpWater);
577 fVacuum = man->FindOrBuildMaterial("G4_Gal
578 fMaterialVect.push_back(fVacuum);
579
580 G4double z, a, density;
581 G4String name, symbol;
582 G4int nComponents, nAtoms;
583
584 a = 12.0107*g/mole;
585 G4Element* elC = new G4Element(name="Carbo
586 a = 1.00794*g/mole;
587 G4Element* elH = new G4Element(name="Hydr
588 a = 15.9994*g/mole;
589 G4Element* elO = new G4Element(name="Oxyg
590 a = 14.00674*g/mole;
591 G4Element* elN = new G4Element(name="Nitr
592 a = 30.973762*g/mole;
593 G4Element* elP = new G4Element(name="Phos
594
595 density = 1.*g/cm3;
596
597 fTHF = new G4Material("THF", density, nCom
598 fTHF->AddElement(elC, nAtoms=4);
599 fTHF->AddElement(elH, nAtoms=8);
600 fTHF->AddElement(elO, nAtoms=1);
601
602 fPY = new G4Material("PY", density, nCompo
603 fPY->AddElement(elC, nAtoms=4);
604 fPY->AddElement(elH, nAtoms=4);
605 fPY->AddElement(elN, nAtoms=2);
606
607 fPU = new G4Material("PU", density, nCompo
608 fPU->AddElement(elC, nAtoms=5);
609 fPU->AddElement(elH, nAtoms=4);
610 fPU->AddElement(elN, nAtoms=4);
611
612 fTMP = new G4Material("TMP", density,4);
613 fTMP->AddElement(elC, nAtoms=3);
614 fTMP->AddElement(elH, nAtoms=6);
615 fTMP->AddElement(elP, nAtoms=1);
616 fTMP->AddElement(elO, nAtoms=4);
617
618 fMaterialVect.push_back(fTHF);
619 fMaterialVect.push_back(fPY);
620 fMaterialVect.push_back(fPU);
621 fMaterialVect.push_back(fTMP);
622
623 // Mixture creation
624 fSugarMixt = new G4Material("sugarMixt", d
625 fSugarMixt->AddMaterial(fTHF,0.5);
626 fSugarMixt->AddMaterial(fTMP,0.5);
627 fMaterialVect.push_back(fSugarMixt);
628
629 // Real DNA materials creation
630 //
631 // Density calculation:
632 //
633 // d=m/V and M=m/n <=> m = M*n
634 // <=> d = M*n/V
635 // <=> d = (M*nb)/(Na*V)
636
637 G4double MBackbone = (5.*12.0104 + 10.*1.0
638 G4double VBackbone = 1.823912/20. * std::p
639 G4double densityBaTHF = (MBackbone/(g/mole
640
641 fDeoxyribose = new G4Material("backbone_TH
642 fDeoxyribose->AddElement(elC, nAtoms=5);
643 fDeoxyribose->AddElement(elH, nAtoms=10);
644 fDeoxyribose->AddElement(elO, nAtoms=4);
645 fMaterialVect.push_back(fDeoxyribose);
646
647 MBackbone = (3.*1.00794 + 1.*30.973762 + 4
648 VBackbone = 1.19114/20. * std::pow(1.e-9,
649 G4double densityBaTMP = (MBackbone/(g/mole
650
651 fPhosphate = new G4Material("backbone_TMP"
652 fPhosphate->AddElement(elH, nAtoms=3);
653 fPhosphate->AddElement(elP, nAtoms=1);
654 fPhosphate->AddElement(elO, nAtoms=4);
655 fMaterialVect.push_back(fPhosphate);
656
657 G4double MCytosine = (4.*12.0104 + 5.*1.00
658 G4double VBase = 1.8527205/20. * std::pow(
659 G4double densityCy = (MCytosine/(g/mole) *
660
661 fCytosine_PY = new G4Material("cytosine_PY
662 fCytosine_PY->AddElement(elC, nAtoms=4);
663 fCytosine_PY->AddElement(elH, nAtoms=5);
664 fCytosine_PY->AddElement(elN, nAtoms=3);
665 fCytosine_PY->AddElement(elO, nAtoms=1);
666 fMaterialVect.push_back(fCytosine_PY);
667
668 G4double MThy = (5.*12.0104 + 6.*1.00794 +
669 VBase = 1.8527205/20. * std::pow(1.e-9, 3.
670 densityCy = (MThy/(g/mole) * 1. )/(6.022e2
671
672 fThymine_PY = new G4Material("thymine_PY",
673 fThymine_PY->AddElement(elC, nAtoms=5);
674 fThymine_PY->AddElement(elH, nAtoms=6);
675 fThymine_PY->AddElement(elN, nAtoms=2);
676 fThymine_PY->AddElement(elO, nAtoms=2);
677 fMaterialVect.push_back(fThymine_PY);
678
679 G4double MAdenine = (5.*12.0104 + 5.*1.007
680 VBase = 1.8527205/20. * std::pow(1.e-9, 3.
681 G4double densityAd = (MAdenine/(g/mole) *
682
683 fAdenine_PU = new G4Material("adenine_PU",
684 fAdenine_PU->AddElement(elC, nAtoms=5);
685 fAdenine_PU->AddElement(elH, nAtoms=5);
686 fAdenine_PU->AddElement(elN, nAtoms=5);
687 fMaterialVect.push_back(fAdenine_PU);
688
689 MAdenine = (5.*12.0104 + 5.*1.00794 + 5.*1
690 VBase = 1.8527205/20. * std::pow(1.e-9, 3.
691 densityAd = (MAdenine/(g/mole) * 1. )/(6.0
692
693 fGuanine_PU = new G4Material("guanine_PU",
694 fGuanine_PU->AddElement(elC, nAtoms=5);
695 fGuanine_PU->AddElement(elH, nAtoms=5);
696 fGuanine_PU->AddElement(elN, nAtoms=5);
697 fGuanine_PU->AddElement(elO, nAtoms=1);
698 fMaterialVect.push_back(fGuanine_PU);
699
700 fHomogeneous_dna = new G4Material("homogen
701 fHomogeneous_dna->AddMaterial(fpWater, 0.3
702 fHomogeneous_dna->AddMaterial(fDeoxyribose
703 fHomogeneous_dna->AddMaterial(fPhosphate,
704 fHomogeneous_dna->AddMaterial(fAdenine_PU,
705 fHomogeneous_dna->AddMaterial(fGuanine_PU,
706 fHomogeneous_dna->AddMaterial(fCytosine_PY
707 fHomogeneous_dna->AddMaterial(fThymine_PY,
708
709 fMaterialVect.push_back(fHomogeneous_dna);
710 }
711
712 //....oooOO0OOooo........oooOO0OOooo........oo
713
714 G4LogicalVolume* PhysGeoImport::CreateNucleusL
715 {
716 // Here, we parse the nucleus file and bui
717
718 G4cout<<"Start the nucleus creation..."<<G
719
720 G4VSolid* nucleusSolid = 0;
721 G4LogicalVolume* nucleusLogic = 0;
722
723 G4bool isNucleusBuild (false);
724
725 // Open the world file
726 std::ifstream nucleusFile;
727 nucleusFile.open(fileName.c_str());
728
729 // Check if the file was correctly opened
730 if(!nucleusFile.is_open())
731 {
732 // Geant4 exception
733 G4String msg = fileName+" could not be
734 G4Exception("PhysGeoImport::ParseFile(
735 }
736
737 // Define the line string variable
738 G4String line;
739
740 // Read the file line per line
741 while(std::getline(nucleusFile, line) && !
742 {
743 // Check the line to determine if it i
744 if(line.empty()) continue; // skip the
745
746 // Data string stream
747 std::istringstream issLine(line);
748
749 // String to determine the first lette
750 G4String firstItem;
751
752 // Put the first letter/word within th
753 issLine >> firstItem;
754
755 // Check first letter to determine if
756 if(firstItem=="#") continue; // skip t
757
758 // Use the file
759
760 else if(firstItem == "_Name")
761 {
762 issLine >> fNucleusName;
763 }
764
765 else if(firstItem=="_Type")
766 {
767 issLine >> fNucleusType;
768
769 if (fNucleusType == "Ellipsoid")
770 {
771 G4float semiX, semiY, semiZ;
772
773 issLine >> semiX >> semiY >> s
774
775 semiX *= fFactor*nm;
776 semiY *= fFactor*nm;
777 semiZ *= fFactor*nm;
778
779 fNucleusData["SemiX"] = semiX;
780 fNucleusData["SemiY"] = semiY;
781 fNucleusData["SemiZ"] = semiZ;
782
783 nucleusSolid = new G4Ellipsoid
784 nucleusLogic = new G4LogicalVo
785 G4double r3 = semiX*semiY*semi
786 fNucleusVolume = 4.0*pi*r3/3.0
787 } else if (fNucleusType == "Ellipt
788 G4float semiX, semiY, semiZ;
789
790 issLine >> semiX >> semiY >> s
791
792 semiX *= fFactor*nm;
793 semiY *= fFactor*nm;
794 semiZ *= fFactor*nm;
795
796 fNucleusData["SemiX"] = semiX;
797 fNucleusData["SemiY"] = semiY;
798 fNucleusData["SemiZ"] = semiZ;
799
800 nucleusSolid = new G4Elliptica
801 nucleusLogic = new G4LogicalVo
802 G4double r3 = 2*semiX*semiY*se
803 fNucleusVolume = pi*r3;
804 } else if (fNucleusType == "Spheri
805 G4float radius;
806 issLine >> radius;
807 radius *= fFactor*nm;
808 fNucleusData["SemiX"] = radius
809 fNucleusData["SemiY"] = radius
810 fNucleusData["SemiZ"] = radius
811 nucleusSolid = new G4Orb(fNucl
812 nucleusLogic = new G4LogicalVo
813 G4double r3 = radius*radius*ra
814 fNucleusVolume = 4.0*pi*r3/3.0
815 } else {
816 G4String msg =fNucleusType+" i
817 G4Exception("PhysGeoImport::Cr
818 "Geo_InputFile_Unknown_Nucleus
819 }
820
821 isNucleusBuild = true;
822 }
823 }
824
825 nucleusFile.close();
826
827 if(!isNucleusBuild)
828 {
829 G4String msg = "Nucleus data were not
830 G4Exception("PhysGeoImport::CreateNucl
831 "Geo_InputFile_Unknown_Nucleus", Fatal
832 }
833
834 return nucleusLogic;
835 }
836
837 //....oooOO0OOooo........oooOO0OOooo........oo
838
839 std::vector<Voxel>* PhysGeoImport::CreateVoxel
840 {
841 fTotalNbBpPlacedInGeo = 0;
842 fTotalNbHistonePlacedInGeo = 0;
843 G4cout<<"Start the voxel data generation..
844
845 std::vector<Voxel>* voxels = new std::vect
846
847 // Clear any previously loaded world
848 voxels->clear();
849
850 voxels->reserve(1000000);
851
852 // Open the world file
853 std::ifstream nucleusFile;
854 nucleusFile.open(fileName.c_str());
855
856 // Check if the file was correctly opened
857 if(!nucleusFile.is_open())
858 {
859 // Geant4 exception
860 G4String msg = fileName+" could not be
861 G4Exception("PhysGeoImport::ParseFile(
862 }
863
864 // Initialise the copy number to 0
865 G4int voxelCopyNumber = 0;
866
867 // Define the line string variable
868 G4String line;
869
870 auto whatchromatintype = [] (Voxel::VoxelT
871 if (voxt == Voxel::Straight || voxt ==
872 voxt == Voxel::Up || voxt == Voxel
873 else if (voxt == Voxel::Straight2 || v
874 voxt == Voxel::Up2 || voxt == Voxe
875 else return ChromatinType::fUnspecifie
876 };
877 // Read the file line per line
878 while(std::getline(nucleusFile, line) )
879 {
880 // Check the line to determine if it i
881 if(line.empty()) continue; // skip the
882
883 // Data string stream
884 std::istringstream issLine(line);
885
886 // String to determine the first lette
887 G4String firstItem;
888
889 // Put the first letter/word within th
890 issLine >> firstItem;
891
892 // Check first letter to determine if
893 if(firstItem=="#") continue; // skip t
894
895 // Use the file
896
897 else if(firstItem == "_pl")
898 {
899 // Set the flags
900
901 G4String name;
902 issLine >> name;
903
904 G4int chromoCpN;
905 issLine >> chromoCpN;
906
907 G4int domainCpN;
908 issLine >> domainCpN;
909
910 G4double x;
911 issLine >> x;
912 x *= fFactor*nm;
913
914 G4double y;
915 issLine >> y;
916 y *= fFactor*nm;
917
918 G4double z;
919 issLine >> z;
920 z *= fFactor*nm;
921
922 G4double xx;
923 issLine >> xx;
924
925 G4double yx;
926 issLine >> yx;
927
928 G4double zx;
929 issLine >> zx;
930
931 G4double xy;
932 issLine >> xy;
933
934 G4double yy;
935 issLine >> yy;
936
937 G4double zy;
938 issLine >> zy;
939
940 G4double xz;
941 issLine >> xz;
942
943 G4double yz;
944 issLine >> yz;
945
946 G4double zz;
947 issLine >> zz;
948
949
950 G4RotationMatrix* rot = new G4Rota
951 G4ThreeVector(yx,yy,yz),G4ThreeVec
952
953 Voxel::VoxelType type = Voxel::Oth
954
955 if(name=="voxelStraight" || name==
956 {
957 type = Voxel::Straight;
958 }
959 else if(name=="voxelUp" || name=="
960 {
961 type = Voxel::Up;
962 }
963 else if(name=="voxelDown" || name=
964 {
965 type = Voxel::Down;
966 }
967 else if(name=="voxelRight" || name
968 {
969 type = Voxel::Right;
970 }
971 else if(name=="voxelLeft" || name=
972 {
973 type = Voxel::Left;
974 }
975 else if(name=="voxelStraight2" ||
976 {
977 type = Voxel::Straight2;
978 }
979 else if(name=="voxelUp2" || name==
980 {
981 type = Voxel::Up2;
982 }
983 else if(name=="voxelDown2" || name
984 {
985 type = Voxel::Down2;
986 }
987 else if(name=="voxelRight2" || nam
988 {
989 type = Voxel::Right2;
990 }
991 else if(name=="voxelLeft2" || name
992 {
993 type = Voxel::Left2;
994 }
995 else
996 {
997 G4ExceptionDescription msg ;
998 msg << "The voxel named "<<nam
999 G4Exception("PhysGeoImport::Cr
1000 }
1001
1002 voxels->push_back(Voxel(voxelCopy
1003 G4ThreeVector(x,y,z), rot) );
1004 fTotalNbBpPlacedInGeo += fVoxelNb
1005 fTotalNbHistonePlacedInGeo += fVo
1006 voxelCopyNumber++;
1007 auto chromatintype = whatchromat
1008 if (fChromatinTypeCount.find(chro
1009 fChromatinTypeCount.insert({c
1010 else fChromatinTypeCount[chromati
1011 }
1012 }
1013
1014 nucleusFile.close();
1015
1016 G4cout<<"End the voxel data generation"<<
1017
1018 return voxels;
1019 }
1020
1021 //....oooOO0OOooo........oooOO0OOooo........o