Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // 27 /// \file PhysGeoImport.cc 28 /// \brief Implementation of the PhysGeoImport class 29 30 #include "PhysGeoImport.hh" 31 #include "Randomize.hh" 32 #include "G4Sphere.hh" 33 #include "G4PhysicalConstants.hh" 34 35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 36 37 PhysGeoImport::PhysGeoImport() : 38 fFactor(1.), 39 fGeoName("") 40 { 41 DefineMaterial(); 42 } 43 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 45 46 PhysGeoImport::PhysGeoImport(G4bool isVisu) : 47 fIsVisu(isVisu), 48 fFactor(1.), 49 fGeoName("") 50 { 51 if (fIsVisu) G4cout<<" For Checking Visualization!"<<G4endl; 52 DefineMaterial(); 53 } 54 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 56 57 G4LogicalVolume* PhysGeoImport::CreateLogicVolume(const G4String& fileName,G4String& voxelName) 58 { 59 // The idea is to return a pointer to the mother volume of the input file as output. 60 // Create a temporary material map to assign materials 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<ie; ++i) 66 { 67 G4String matName(""); 68 69 if(fMaterialVect[i]->GetName()=="adenine_PU") matName = "adenine"; 70 else if(fMaterialVect[i]->GetName()=="thymine_PY") matName = "thymine"; 71 else if(fMaterialVect[i]->GetName()=="guanine_PU") matName = "guanine"; 72 else if(fMaterialVect[i]->GetName()=="cytosine_PY") matName = "cytosine"; 73 else if(fMaterialVect[i]->GetName()=="backbone_THF") matName = "deoxyribose"; 74 else if(fMaterialVect[i]->GetName()=="backbone_TMP") matName = "phosphate"; 75 else if(fMaterialVect[i]->GetName()=="G4_WATER") matName = "water"; 76 else if(fMaterialVect[i]->GetName()=="homogeneous_dna") matName = "homogeneous_dna"; 77 else if(fMaterialVect[i]->GetName()=="G4_Galactic") matName = "vacuum"; 78 else 79 { 80 matName = fMaterialVect[i]->GetName(); 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, fSize/2, fSize/2, fSize/2); 95 96 G4String boxNameLogic = fGeoName+"_logic"; 97 G4LogicalVolume* box_logic = new G4LogicalVolume(box_solid, fpWater, boxNameLogic); 98 99 // sort the molecules 100 std::sort(fMolecules.begin(), fMolecules.end() ); 101 102 // Loop on all the parsed molecules 103 for(G4int i=0, ie=fMolecules.size(); i<ie; ++i) 104 { 105 // Retrieve general molecule informations 106 // 107 G4String name = fMolecules[i].fName; 108 G4String materialName = "water";// fMolecules[i].fMaterial; 109 G4double radius = fMolecules[i].fRadius; 110 G4double waterRadius = fMolecules[i].fRadiusWater; 111 G4ThreeVector moleculePosition = fMolecules[i].fPosition; 112 G4int copyNum = fMolecules[i].fCopyNumber; 113 114 // Water hydration shell volume part 115 116 G4Orb* moleculeWater_solid = 0; 117 G4VSolid* moleculeWaterCut_solid = 0; 118 G4LogicalVolume* moleculeWater_logic = 0; 119 120 // If water radius != 0 then we have a water hydration shell 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+"_water_solid"; 129 moleculeWater_solid = new G4Orb(nameWaterSolid, waterRadius); 130 moleculeWaterCut_solid = CreateCutSolid(moleculeWater_solid, fMolecules[i], fMolecules, false); 131 132 G4String nameWaterLogic = name+"_water_logic"; 133 moleculeWater_logic = new G4LogicalVolume(moleculeWaterCut_solid, mat, nameWaterLogic); 134 135 G4String nameWaterPhys = name+"_water_phys"; 136 new G4PVPlacement(0, moleculePosition, moleculeWater_logic, nameWaterPhys, box_logic, false, copyNum); 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].fName+"_solid"; 146 molecule_solid = new G4Orb(nameSolid, radius); 147 moleculeCut_solid = CreateCutSolid(molecule_solid, fMolecules[i], fMolecules, true); 148 149 G4String nameLogic = name+"_logic"; 150 molecule_logic = new G4LogicalVolume(moleculeCut_solid, materials[materialName], nameLogic); 151 152 // If there was a water hydration shell volume then the current dna volume is 153 // placed within it and, thus, its relative coordinates are 0,0,0. 154 if(waterRadius > (0 + tol)*nm) 155 { 156 G4ThreeVector position(0.,0.,0.); 157 G4String namePhys = name+"_phys"; 158 new G4PVPlacement(0, position, molecule_logic, namePhys, moleculeWater_logic, false, copyNum); 159 } 160 // If not, coordinates are those of the molecule 161 else 162 { 163 G4ThreeVector position = moleculePosition; 164 G4String namePhys = name+"_phys"; 165 new G4PVPlacement(0, position, molecule_logic, namePhys, box_logic, false, copyNum); 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........oooOO0OOooo........oooOO0OOooo...... 178 179 G4String PhysGeoImport::ParseFile(const G4String& fileName) 180 { 181 G4cout<<"Start parsing of "<<fileName<<G4endl; 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 opened"; 196 G4Exception("PhysGeoImport::ParseFile()", "Geo_InputFileNotOpened", FatalException, msg); 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 is empty 209 if(line.empty()) continue; // skip the line if it is empty 210 211 // Data string stream 212 std::istringstream issLine(line); 213 214 // String to determine the first letter/word 215 G4String firstItem; 216 217 // Put the first letter/word within the string 218 issLine >> firstItem; 219 220 // Check first letter to determine if the line is data or comment 221 if(firstItem=="#") continue; // skip the line if it is comment 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] = waterRadius; 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 != preHistoneCpNo ) { 279 NoHistone ++; 280 preHistoneCpNo = copyNumber; 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, G4ThreeVector(x, y, z), fRadiusMap[name], fWaterRadiusMap[name], material, strand); 302 303 fMolecules.push_back(molecule); 304 } 305 else 306 { 307 // Geant4 exception 308 G4String msg = firstItem+" is not defined in the parser. Check the input file: "+fileName+"."; 309 G4Exception("PhysGeoImport::ParseFile()", "Geo_WrongParse", FatalException, msg); 310 } 311 } 312 313 // Close the file once the reading is done 314 file.close(); 315 316 G4cout<<"End parsing of "<<fileName<<G4endl; 317 fVoxelNbHistoneMap.insert({fGeoName,NoHistone}); 318 fVoxelNbBpMap.insert({fGeoName,NoBp}); 319 return fGeoName; 320 } 321 322 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 323 324 G4VSolid* PhysGeoImport::CreateCutSolid(G4Orb *solidOrbRef, 325 Molecule &molRef, 326 std::vector<Molecule> &molList, 327 G4bool in) 328 { 329 // The idea behing this method is to cut overlap volumes by selecting one of them (the reference) and checking all the other volumes (the targets). 330 // If a reference and a target volumes are close enough to overlap they will be cut. 331 // The reference is already selected when we enter this method. 332 333 // Use the tiny space to differentiate the frontiers (may not be necessary) 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 = molRef.fRadius; 346 else radiusRef = molRef.fRadiusWater; 347 348 // Reference volume position 349 G4ThreeVector posRef = molRef.fPosition; 350 // Before looping on all the volumes we check if the current 351 //reference volume overlaps with its container voxel boundaries. 352 if(std::abs(posRef.x() ) + radiusRef > fSize/2 // along x 353 || std::abs(posRef.y() ) + radiusRef > fSize/2 // along y 354 || std::abs(posRef.z() ) + radiusRef > fSize/2) // along z 355 { 356 // If we enter here, then the reference volume overlaps with the boundaries of the container voxel 357 // Box used to cut 358 G4Box* solidBox = new G4Box("solid_box_for_cut", fSize/2, fSize/2, fSize/2); 359 G4ThreeVector posBox; 360 361 // Create a dummy rotation matrix 362 G4RotationMatrix *rotMat = new G4RotationMatrix; 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 ) > fSize/2 ) 369 { 370 posBox = -posRef + G4ThreeVector(0,fSize,0); 371 372 // If the volume is cutted for the first time 373 if(!isCutted) 374 { 375 solidCut = new G4SubtractionSolid("solidCut", solidOrbRef, solidBox, rotMat, posBox); 376 isCutted = true; 377 } 378 // For the other times 379 else solidCut = new G4SubtractionSolid("solidCut", solidCut, solidBox, rotMat, posBox); 380 } 381 382 // Down 383 if(std::abs( posRef.y() - radiusRef ) > fSize/2 ) 384 { 385 posBox = -posRef + G4ThreeVector(0,-fSize,0); 386 387 // If the volume is cutted for the first time 388 if(!isCutted) 389 { 390 solidCut = new G4SubtractionSolid("solidCut", solidOrbRef,solidBox,rotMat, posBox); 391 isCutted = true; 392 } 393 // For the other times 394 else solidCut = new G4SubtractionSolid("solidCut", solidCut, solidBox, rotMat, posBox); 395 } 396 397 // Left 398 if(std::abs( posRef.x() + radiusRef ) > fSize/2 ) 399 { 400 posBox = -posRef + G4ThreeVector(fSize,0,0); 401 402 // If the volume is cutted for the first time 403 if(!isCutted) 404 { 405 solidCut = new G4SubtractionSolid("solidCut", solidOrbRef, solidBox, rotMat, posBox); 406 isCutted = true; 407 } 408 // For the other times 409 else solidCut = new G4SubtractionSolid("solidCut", solidCut, solidBox, rotMat, posBox); 410 } 411 412 // Right 413 if(std::abs( posRef.x() - radiusRef ) > fSize/2 ) 414 { 415 posBox = -posRef + G4ThreeVector(-fSize,0,0); 416 417 // If the volume is cutted for the first time 418 if(!isCutted) 419 { 420 solidCut = new G4SubtractionSolid("solidCut", solidOrbRef, solidBox, rotMat, posBox); 421 isCutted = true; 422 } 423 // For the other times 424 else solidCut = new G4SubtractionSolid("solidCut", solidCut, solidBox, rotMat, posBox); 425 } 426 427 // Forward 428 if(std::abs( posRef.z() + radiusRef ) > fSize/2 ) 429 { 430 posBox = -posRef + G4ThreeVector(0,0,fSize); 431 432 // If the volume is cutted for the first time 433 if(!isCutted) 434 { 435 solidCut = new G4SubtractionSolid("solidCut", solidOrbRef, solidBox, rotMat, posBox); 436 isCutted = true; 437 } 438 // For the other times 439 else solidCut = new G4SubtractionSolid("solidCut", solidCut, solidBox, rotMat, posBox); 440 } 441 442 // Backward 443 if(std::abs( posRef.z() - radiusRef ) > fSize/2 ) 444 { 445 posBox = -posRef + G4ThreeVector(0,0,-fSize); 446 447 // If the volume is cutted for the first time 448 if(!isCutted) 449 { 450 solidCut = new G4SubtractionSolid("solidCut", solidOrbRef, solidBox, rotMat, posBox); 451 isCutted = true; 452 } 453 // For the other times 454 else solidCut = new G4SubtractionSolid("solidCut", solidCut, solidBox, rotMat, posBox); 455 } 456 } 457 458 // Look the other volumes of the voxel 459 // Loop on all the target volumes (other volumes with potential overlaps) 460 for(G4int i=0, ie=molList.size(); i<ie; ++i) 461 { 462 G4ThreeVector posTar = molList[i].fPosition; 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 informations 477 G4double radiusTar; 478 if(molList[i].fRadiusWater==0) radiusTar = molList[i].fRadius; 479 else radiusTar = molList[i].fRadiusWater; 480 481 // Compute the distance reference-target 482 G4double distance = std::abs( (posRef - posTar).getR() ); 483 484 // Use the distance to check if the current target is also the reference. 485 // This can only happen once per loop. 486 if(distance==0 && !isOurVol) 487 { 488 // Target volume is also reference volume. 489 490 // Set the flag 491 isOurVol = true; 492 493 // Next iteration 494 continue; 495 } 496 // If the condition is correct more than one time then there is a mistake somewhere. 497 else if(distance == 0 && isOurVol) 498 { 499 G4ExceptionDescription desmsg; 500 desmsg << "DetectorConstruction::CreateCutSolid: Two volumes are placed at the same position."; 501 G4Exception("ChemGeoImport::BuildGeometry", "", FatalException, desmsg); 502 } 503 504 // If the volumes are differents then we want to know if they are 505 // close enough to overlap and, thus, to intiate a cut. 506 else if(distance <= radiusRef+radiusTar) 507 { 508 // Volumes are close enough, there will be a cut 509 510 // Box used to cut 511 G4Box* solidBox = new G4Box("solid_box_for_cut", 2*radiusTar, 2*radiusTar, 2*radiusTar); 512 513 // This part is tricky. 514 // The goal is to calculate the position of the intersection center 515 516 // diff vector to from ref to tar 517 G4ThreeVector diff = posTar - posRef; 518 519 // Find the intersection point and add to it half the length of the box used to cut 520 G4double d = (std::pow(radiusRef,2)-std::pow(radiusTar,2)+std::pow(distance,2) ) / 521 (2*distance) + solidBox->GetZHalfLength() - tinySpace; 522 523 // If we are in another volume we choose to double the tinySpace to 524 //differentiate without ambiguities the inner and outer volume frontiers. 525 // (may not be necessary) 526 if(in) d -= 2*tinySpace; 527 528 // Position of the box in order to achieve the cut. 529 // "* ( diff/diff.getR() )" is necessary to get a vector in the right direction as output 530 G4ThreeVector pos = d *( diff/diff.getR() ); 531 532 // Get the rotation angles because the box used to cut needs to be rotated 533 // to give the right "cut angle". 534 G4double phi = std::acos(pos.getZ()/pos.getR()); 535 G4double theta = std::acos( pos.getX() / ( pos.getR()*std::cos(pi/2.-phi) ) ); 536 537 if(pos.getY()<0) theta = -theta; 538 539 G4ThreeVector rotAxisForPhi(1*fFactor*nm,0.,0.); 540 rotAxisForPhi.rotateZ(theta+pi/2); 541 542 // Create the rotation matrix 543 G4RotationMatrix *rotMat = new G4RotationMatrix; 544 rotMat->rotate(-phi, rotAxisForPhi); 545 546 // Rotate it again 547 G4ThreeVector rotZAxis(0.,0.,1*fFactor*nm); 548 rotMat->rotate(theta, rotZAxis); 549 550 // If the volume is cutted for the first time 551 if(!isCutted) solidCut = new G4SubtractionSolid("solidCut", solidOrbRef, 552 solidBox, rotMat, pos); 553 554 // For the other times 555 else solidCut = new G4SubtractionSolid("solidCut", solidCut, solidBox, rotMat, pos); 556 557 // Set the cut flag 558 isCutted = true; 559 } 560 } 561 562 // If there was at least one cut then we return the cutted volume 563 if(isCutted) return solidCut; 564 565 // Otherwise, we return the original volume 566 else return solidOrbRef; 567 568 } 569 570 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 571 572 void PhysGeoImport::DefineMaterial() 573 { 574 G4NistManager * man = G4NistManager::Instance(); 575 fpWater = man->FindOrBuildMaterial("G4_WATER"); 576 fMaterialVect.push_back(fpWater); 577 fVacuum = man->FindOrBuildMaterial("G4_Galactic"); 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="Carbon", symbol="C", z=6., a); 586 a = 1.00794*g/mole; 587 G4Element* elH = new G4Element(name="Hydrogen",symbol="H" , z= 1., a); 588 a = 15.9994*g/mole; 589 G4Element* elO = new G4Element(name="Oxygen" ,symbol="O" , z= 8., a); 590 a = 14.00674*g/mole; 591 G4Element* elN = new G4Element(name="Nitrogen" ,symbol="N" , z= 7., a); 592 a = 30.973762*g/mole; 593 G4Element* elP = new G4Element(name="Phosphorus" ,symbol="P" , z= 15., a); 594 595 density = 1.*g/cm3; 596 597 fTHF = new G4Material("THF", density, nComponents=3); 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, nComponents=3); 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, nComponents=3); 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", density, nComponents=2); 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.00794 + 4.*15.9994)*g/mole; 638 G4double VBackbone = 1.823912/20. * std::pow(1.e-9, 3.) *m3; 639 G4double densityBaTHF = (MBackbone/(g/mole) * 1. )/(6.022e23*VBackbone/cm3) *g/cm3; // g/cm3 640 641 fDeoxyribose = new G4Material("backbone_THF", densityBaTHF, nComponents=3); 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.*15.9994)*g/mole; 648 VBackbone = 1.19114/20. * std::pow(1.e-9, 3.) *m3; 649 G4double densityBaTMP = (MBackbone/(g/mole) * 1. )/(6.022e23*VBackbone/cm3) *g/cm3; // g/cm3 650 651 fPhosphate = new G4Material("backbone_TMP", densityBaTMP,3); 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.00794 + 3.*14.0067 + 1.*15.9994)*g/mole; 658 G4double VBase = 1.8527205/20. * std::pow(1.e-9, 3.) *m3; 659 G4double densityCy = (MCytosine/(g/mole) * 1. )/(6.022e23*VBase/cm3) *g/cm3; // g/cm3 660 661 fCytosine_PY = new G4Material("cytosine_PY", densityCy, nComponents=4); 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 + 2.*14.0067 + 2.*15.9994)*g/mole; 669 VBase = 1.8527205/20. * std::pow(1.e-9, 3.) *m3; 670 densityCy = (MThy/(g/mole) * 1. )/(6.022e23*VBase/cm3) *g/cm3; // g/cm3 671 672 fThymine_PY = new G4Material("thymine_PY", densityCy, nComponents=4); 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.00794 + 5.*14.0067)*g/mole; 680 VBase = 1.8527205/20. * std::pow(1.e-9, 3.) *m3; 681 G4double densityAd = (MAdenine/(g/mole) * 1. )/(6.022e23*VBase/cm3) *g/cm3; // g/cm3 682 683 fAdenine_PU = new G4Material("adenine_PU", densityAd, nComponents=3); 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.*14.0067 + 1.*15.999)*g/mole; 690 VBase = 1.8527205/20. * std::pow(1.e-9, 3.) *m3; 691 densityAd = (MAdenine/(g/mole) * 1. )/(6.022e23*VBase/cm3) *g/cm3; // g/cm3 692 693 fGuanine_PU = new G4Material("guanine_PU", densityAd, nComponents=4); 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("homogeneous_dna", 19.23e-21/(12.30781*1.e-21) *g/cm3, nComponents=7);//21 701 fHomogeneous_dna->AddMaterial(fpWater, 0.37); 702 fHomogeneous_dna->AddMaterial(fDeoxyribose, 0.23); 703 fHomogeneous_dna->AddMaterial(fPhosphate, 0.17); 704 fHomogeneous_dna->AddMaterial(fAdenine_PU, 0.06); 705 fHomogeneous_dna->AddMaterial(fGuanine_PU, 0.07); 706 fHomogeneous_dna->AddMaterial(fCytosine_PY, 0.05); 707 fHomogeneous_dna->AddMaterial(fThymine_PY, 0.05); 708 709 fMaterialVect.push_back(fHomogeneous_dna); 710 } 711 712 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 713 714 G4LogicalVolume* PhysGeoImport::CreateNucleusLogicVolume(const G4String& fileName) 715 { 716 // Here, we parse the nucleus file and build a logical Geant4 volume from it. 717 718 G4cout<<"Start the nucleus creation..."<<G4endl; 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 opened"; 734 G4Exception("PhysGeoImport::ParseFile()", "Geo_InputFileNotOpened", FatalException, msg); 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) && !isNucleusBuild) 742 { 743 // Check the line to determine if it is empty 744 if(line.empty()) continue; // skip the line if it is empty 745 746 // Data string stream 747 std::istringstream issLine(line); 748 749 // String to determine the first letter/word 750 G4String firstItem; 751 752 // Put the first letter/word within the string 753 issLine >> firstItem; 754 755 // Check first letter to determine if the line is data or comment 756 if(firstItem=="#") continue; // skip the line if it is comment 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 >> semiZ; 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(fNucleusName, semiX, semiY, semiZ); 784 nucleusLogic = new G4LogicalVolume(nucleusSolid, fpWater, "nucleus_logic"); 785 G4double r3 = semiX*semiY*semiZ/m3; 786 fNucleusVolume = 4.0*pi*r3/3.0; 787 } else if (fNucleusType == "EllipticCylinder") { 788 G4float semiX, semiY, semiZ; 789 790 issLine >> semiX >> semiY >> semiZ; 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 G4EllipticalTube(fNucleusName, semiX, semiY, semiZ); 801 nucleusLogic = new G4LogicalVolume(nucleusSolid, fpWater, "nucleus_logic"); 802 G4double r3 = 2*semiX*semiY*semiZ/m3; // multiplied by 2 cuz semiZ is hafl of height 803 fNucleusVolume = pi*r3; 804 } else if (fNucleusType == "Spherical") { 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(fNucleusName,radius); 812 nucleusLogic = new G4LogicalVolume(nucleusSolid, fpWater, "nucleus_logic"); 813 G4double r3 = radius*radius*radius/m3; 814 fNucleusVolume = 4.0*pi*r3/3.0; 815 } else { 816 G4String msg =fNucleusType+" is not registered."; 817 G4Exception("PhysGeoImport::CreateNucleusLogicVolume()", 818 "Geo_InputFile_Unknown_Nucleus", FatalException, msg); 819 } 820 821 isNucleusBuild = true; 822 } 823 } 824 825 nucleusFile.close(); 826 827 if(!isNucleusBuild) 828 { 829 G4String msg = "Nucleus data were not found for "+fNucleusType; 830 G4Exception("PhysGeoImport::CreateNucleusLogicVolume()", 831 "Geo_InputFile_Unknown_Nucleus", FatalException, msg); 832 } 833 834 return nucleusLogic; 835 } 836 837 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 838 839 std::vector<Voxel>* PhysGeoImport::CreateVoxelsData(const G4String& fileName) 840 { 841 fTotalNbBpPlacedInGeo = 0; 842 fTotalNbHistonePlacedInGeo = 0; 843 G4cout<<"Start the voxel data generation..."<<G4endl; 844 845 std::vector<Voxel>* voxels = new std::vector<Voxel>; 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 opened"; 861 G4Exception("PhysGeoImport::ParseFile()", "Geo_InputFileNotOpened", FatalException, msg); 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::VoxelType voxt) -> ChromatinType { 871 if (voxt == Voxel::Straight || voxt == Voxel::Left || voxt == Voxel::Right || 872 voxt == Voxel::Up || voxt == Voxel::Down) return ChromatinType::fHeterochromatin; 873 else if (voxt == Voxel::Straight2 || voxt == Voxel::Left2 || voxt == Voxel::Right2 || 874 voxt == Voxel::Up2 || voxt == Voxel::Down2) return ChromatinType::fEuchromatin; 875 else return ChromatinType::fUnspecified; 876 }; 877 // Read the file line per line 878 while(std::getline(nucleusFile, line) ) 879 { 880 // Check the line to determine if it is empty 881 if(line.empty()) continue; // skip the line if it is empty 882 883 // Data string stream 884 std::istringstream issLine(line); 885 886 // String to determine the first letter/word 887 G4String firstItem; 888 889 // Put the first letter/word within the string 890 issLine >> firstItem; 891 892 // Check first letter to determine if the line is data or comment 893 if(firstItem=="#") continue; // skip the line if it is comment 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 G4RotationMatrix(G4ThreeVector(xx,xy,xz), 951 G4ThreeVector(yx,yy,yz),G4ThreeVector(zx,zy,zz)); 952 953 Voxel::VoxelType type = Voxel::Other; 954 955 if(name=="voxelStraight" || name=="VoxelStraight") 956 { 957 type = Voxel::Straight; 958 } 959 else if(name=="voxelUp" || name=="VoxelUp") 960 { 961 type = Voxel::Up; 962 } 963 else if(name=="voxelDown" || name=="VoxelDown") 964 { 965 type = Voxel::Down; 966 } 967 else if(name=="voxelRight" || name=="VoxelRight") 968 { 969 type = Voxel::Right; 970 } 971 else if(name=="voxelLeft" || name=="VoxelLeft") 972 { 973 type = Voxel::Left; 974 } 975 else if(name=="voxelStraight2" || name=="VoxelStraight2") 976 { 977 type = Voxel::Straight2; 978 } 979 else if(name=="voxelUp2" || name=="VoxelUp2") 980 { 981 type = Voxel::Up2; 982 } 983 else if(name=="voxelDown2" || name=="VoxelDown2") 984 { 985 type = Voxel::Down2; 986 } 987 else if(name=="voxelRight2" || name=="VoxelRight2") 988 { 989 type = Voxel::Right2; 990 } 991 else if(name=="voxelLeft2" || name=="VoxelLeft2") 992 { 993 type = Voxel::Left2; 994 } 995 else 996 { 997 G4ExceptionDescription msg ; 998 msg << "The voxel named "<<name<<" is not registered in the parser"; 999 G4Exception("PhysGeoImport::CreateVoxelsData", "", FatalException, msg, ""); 1000 } 1001 1002 voxels->push_back(Voxel(voxelCopyNumber, chromoCpN, domainCpN, type, 1003 G4ThreeVector(x,y,z), rot) ); 1004 fTotalNbBpPlacedInGeo += fVoxelNbBpMap[name]; 1005 fTotalNbHistonePlacedInGeo += fVoxelNbHistoneMap[name]; 1006 voxelCopyNumber++; 1007 auto chromatintype = whatchromatintype(type); 1008 if (fChromatinTypeCount.find(chromatintype) == fChromatinTypeCount.end()) 1009 fChromatinTypeCount.insert({chromatintype,1}); 1010 else fChromatinTypeCount[chromatintype]++; 1011 } 1012 } 1013 1014 nucleusFile.close(); 1015 1016 G4cout<<"End the voxel data generation"<<G4endl; 1017 1018 return voxels; 1019 } 1020 1021 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......