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