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 DetectorConstruction.cc 28 /// \brief Implementation of the DetectorConst 29 30 #include "DetectorConstruction.hh" 31 #include "Analysis.hh" 32 #include "DetectorConstructionMessenger.hh" 33 34 #include "G4SystemOfUnits.hh" 35 #include "G4Region.hh" 36 #include "G4ProductionCuts.hh" 37 #include "G4UserLimits.hh" 38 #include "G4NistManager.hh" 39 #include "G4RunManager.hh" 40 #include "G4UnionSolid.hh" 41 #include "G4SubtractionSolid.hh" 42 #include "G4Filesystem.hh" 43 44 RunningMode gRunMode = RunningMode::Phys; 45 46 //....oooOO0OOooo........oooOO0OOooo........oo 47 48 DetectorConstruction::DetectorConstruction(G4d 49 G4VUserDetectorConstruction(), fFactor(fac 50 { 51 fDetectorMessenger = new DetectorConstruct 52 fWorldBoxSizeX = fWorldBoxSizeY = fWorldBo 53 if (gRunMode == RunningMode::Chem) fChemGe 54 } 55 56 //....oooOO0OOooo........oooOO0OOooo........oo 57 58 G4VPhysicalVolume* DetectorConstruction::Const 59 60 { 61 G4NistManager * man = G4NistManager::Insta 62 fWater = man->FindOrBuildMaterial("G4_WATE 63 if (gRunMode == RunningMode::Phys) Constru 64 else if (gRunMode == RunningMode::Chem) Co 65 else { 66 // only a world volume 67 G4ExceptionDescription msg; 68 msg <<"Only world volume is constructe 69 <<" Make sure you choose the corre 70 <<" Ignore this message if you int 71 G4Exception("DetectorConstruction::Con 72 "", JustWarning, msg); 73 fSolidWorld = new G4Box("solidWorld",f 74 fLogicWorld = new G4LogicalVolume(fSol 75 fPhysWorld = new G4PVPlacement(0,G4Thr 76 fLogicWorld, 0, false, false); 77 } 78 return fPhysWorld; 79 } 80 81 //....oooOO0OOooo........oooOO0OOooo........oo 82 83 G4VPhysicalVolume * DetectorConstruction::Cons 84 { 85 PhysGeoImport geo(fBVisu); 86 geo.SetFactor(fFactor); 87 88 std::map<G4String, G4LogicalVolume*> voxel 89 // Create an empty logical nucleus 90 91 G4LogicalVolume* logicNucleus = geo.Create 92 if (!fUsingUserDefinedSizesForWorld) { 93 G4double scale = 2.0; 94 fWorldBoxSizeX = scale*geo.GetNucleusS 95 fWorldBoxSizeY = scale*geo.GetNucleusS 96 fWorldBoxSizeZ = geo.GetNucleusSizeDat 97 } 98 99 G4cout<<"===>>World box sizes: SemiX = "<< 100 <<G4BestUnit(fWorldBoxSizeY,"Length" 101 <<G4BestUnit(fWorldBoxSizeZ,"Length" 102 fSolidWorld = new G4Box("solidWorld",fWorl 103 fLogicWorld = new G4LogicalVolume(fSolidWo 104 fPhysWorld = new G4PVPlacement(0,G4ThreeVe 105 // then place cell nucleus in world 106 new G4PVPlacement(0, G4ThreeVector(), logi 107 108 G4LogicalVolume* logicStraightVoxel{nullpt 109 *logicLeftVoxel{nullpt 110 G4LogicalVolume* logicStraightVoxel2{nullp 111 , *logicLeftVoxel2{nul 112 113 for (const auto & entry : fVoxelDefFilesLi 114 G4String name=""; 115 auto ptemp = geo.CreateLogicVolume(ent 116 if (name=="voxelStraight" || name=="Vo 117 logicStraightVoxel = ptemp; 118 voxelMap["VoxelStraight"] = logicS 119 } 120 if (name=="voxelUp" || name=="VoxelUp" 121 logicUpVoxel = ptemp; 122 voxelMap["VoxelUp"] = logicUpVoxel 123 } 124 if (name=="voxelDown" || name=="VoxelD 125 logicDownVoxel = ptemp; 126 voxelMap["VoxelDown"] = logicDownV 127 } 128 if (name=="voxelRight" || name=="Voxel 129 logicRightVoxel = ptemp; 130 voxelMap["VoxelRight"] = logicRigh 131 } 132 if (name=="voxelLeft" || name=="VoxelL 133 logicLeftVoxel = ptemp; 134 voxelMap["VoxelLeft"] = logicLeftV 135 } 136 if (name=="voxelStraight2" || name=="V 137 logicStraightVoxel2 = ptemp; 138 voxelMap["VoxelStraight2"] = logic 139 } 140 if (name=="voxelUp2" || name=="VoxelUp 141 logicUpVoxel2 = ptemp; 142 voxelMap["VoxelUp2"] = logicUpVoxe 143 } 144 if (name=="voxelDown2" || name=="Voxel 145 logicDownVoxel2 = ptemp; 146 voxelMap["VoxelDown2"] = logicDown 147 } 148 if (name=="voxelRight2" || name=="Voxe 149 logicRightVoxel2 = ptemp; 150 voxelMap["VoxelRight2"] = logicRig 151 } 152 if (name=="voxelLeft2" || name=="Voxel 153 logicLeftVoxel2 = ptemp; 154 voxelMap["VoxelLeft2"] = logicLeft 155 } 156 } 157 158 // Create the voxel data table 159 std::vector<Voxel>* voxelTable = geo.Creat 160 const G4int nucleusSize = voxelTable->size 161 G4cout<<"================================= 162 G4cout<<"=====> Number of Histones in each 163 for (auto [key, value] : geo.GetVoxelNbHis 164 G4cout<<key<<": \t\t\t"<<value<<G4endl 165 } 166 G4cout<<"=====> Number of Basepairs in eac 167 for (auto [key, value] : geo.GetVoxelNbBpM 168 G4cout<<key<<": \t\t\t"<<value<<G4endl 169 } 170 G4cout <<"=====> Total Number of Histones 171 <<geo.GetTotalNbHistonePlacedInGeo 172 G4cout <<"=====> Total Number of Basepair 173 <<geo.GetTotalNbBpPlacedInGeo()<<G 174 G4cout <<"=====> Number of each chromatin 175 for (auto [key, value] : geo.GetChromatinT 176 G4String chromatinname = "Heterochroma 177 if (key == ChromatinType::fEuchromatin 178 G4cout<<chromatinname<<": \t\t\t"<<val 179 } 180 G4cout<<"================================= 181 182 // Create the voxel parameterisation 183 if (voxelMap.size() > 0) { 184 // The following dummy declarations ar 185 // in VoxelParameterisation to prevent 186 G4double prevZpos = -fWorldBoxSizeZ; 187 for (auto it = voxelMap.begin(); it != 188 G4double posX = fWorldBoxSizeX - g 189 G4double posY = fWorldBoxSizeY - g 190 G4double posZ = prevZpos + geo.Get 191 new G4PVPlacement(0, G4ThreeVector 192 it->second, "xx", 193 } 194 // End dummy declaration 195 G4int nthreads=0; 196 #ifdef G4MULTITHREADED 197 nthreads = G4RunManager::GetRunManager 198 #endif 199 if ( (nthreads >1) && (voxelMap.size() 200 G4ExceptionDescription msg; 201 msg <<"Number of thread is "<<nthr 202 <<"\nThere will be no DNA cons 203 G4Exception("DetectorConstruction: 204 "RunningMode", JustWar 205 } else { 206 G4VPVParameterisation* voxelParam 207 new G4PVParameterised("VoxelParam" 208 logicNucleus, kUnd 209 } 210 } else { 211 G4ExceptionDescription msg; 212 msg <<"It seems that voxel-definition 213 <<" There will be no DNA constitue 214 G4Exception("DetectorConstruction::Con 215 "Geo_InputFile_NoFile", Ju 216 } 217 Analysis::GetAnalysis()->RecordVoxelDefFil 218 Analysis::GetAnalysis()->SetTotalNbBpPlace 219 Analysis::GetAnalysis()->SetTotalNbHistone 220 Analysis::GetAnalysis()->SetNucleusVolume( 221 Analysis::GetAnalysis()->SetNucleusMassDen 222 logicNucleus->GetMaterial()->GetDensit 223 return fPhysWorld; 224 } 225 226 //....oooOO0OOooo........oooOO0OOooo........oo 227 228 G4VPhysicalVolume *DetectorConstruction::Const 229 { 230 if (fWorldBoxSizeX < fVoxelHalfSizeXYZ) { 231 fWorldBoxSizeX = fWorldBoxSizeY = fWor 232 } 233 fSolidWorld = new G4Box("solidWorld", fWor 234 fLogicWorld = new G4LogicalVolume(fSolidWo 235 fPhysWorld = new G4PVPlacement(0,G4ThreeVe 236 237 if (fVoxelHalfSizeXYZ>0) { 238 G4cout<<"=====> Construct voxel ...... 239 G4Box* solidVoxel = new G4Box("solidVo 240 G4LogicalVolume* logicVoxel = new G4Lo 241 new G4PVPlacement(0, G4ThreeVector(), 242 } 243 return fPhysWorld; 244 } 245 246 //....oooOO0OOooo........oooOO0OOooo........oo 247 248 void DetectorConstruction::SetCellDefFilePath( 249 { 250 const G4fs::path thisP{std::string(finput) 251 G4fs::file_status fst = G4fs::file_status{ 252 auto isExist = G4fs::status_known(fst) ? G 253 if (! isExist) { 254 G4String msg = "File " + finput + "doe 255 G4Exception("DetectorConstruction::Set 256 "Geo_InputFileNotOpened", FatalExcepti 257 } else { 258 fCellDefFilePath = finput; 259 Analysis::GetAnalysis()->RecordCellDef 260 } 261 } 262 263 //....oooOO0OOooo........oooOO0OOooo........oo 264 265 void DetectorConstruction::AddVoxelDefFile(con 266 { 267 const G4fs::path thisP{std::string(finput) 268 G4fs::file_status fst = G4fs::file_status{ 269 auto isExist = G4fs::status_known(fst) ? G 270 if (! isExist) { 271 G4String msg = "File " + finput + "doe 272 G4Exception("DetectorConstruction::Add 273 "Geo_InputFileNotOpened", FatalExcepti 274 } else { 275 fVoxelDefFilesList.insert(finput); 276 } 277 } 278 279 //....oooOO0OOooo........oooOO0OOooo........oo 280 281 void DetectorConstruction::SetWorldBoxSizes(G4 282 { 283 G4double sizex = v.getX(); 284 G4double sizey = v.getY(); 285 G4double sizez = v.getZ(); 286 if (sizex > 0. && sizey > 0. && sizez > 0. 287 fWorldBoxSizeX = sizex; 288 fWorldBoxSizeY = sizey; 289 fWorldBoxSizeZ = sizez; 290 fUsingUserDefinedSizesForWorld = true; 291 } else { 292 G4ExceptionDescription msg ; 293 msg << " Check your setting? One world 294 G4Exception("DetectorConstruction::Set 295 "Geo_WorldSizes", FatalException, msg) 296 } 297 } 298 299 //....oooOO0OOooo........oooOO0OOooo........oo 300 301 void DetectorConstruction::ParseGeoFileForChem 302 { 303 fChemGeoImport->ParseFiles(fn); 304 fVoxelHalfSizeXYZ = fChemGeoImport->GetSiz 305 }