Warning, file /geant4/examples/advanced/dna/dsbandrepair/src/DetectorConstruction.cc was not indexed
or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030 #include "DetectorConstruction.hh"
0031 #include "Analysis.hh"
0032 #include "DetectorConstructionMessenger.hh"
0033
0034 #include "G4SystemOfUnits.hh"
0035 #include "G4Region.hh"
0036 #include "G4ProductionCuts.hh"
0037 #include "G4UserLimits.hh"
0038 #include "G4NistManager.hh"
0039 #include "G4RunManager.hh"
0040 #include "G4UnionSolid.hh"
0041 #include "G4SubtractionSolid.hh"
0042 #include "G4Filesystem.hh"
0043
0044 RunningMode gRunMode = RunningMode::Phys;
0045
0046
0047
0048 DetectorConstruction::DetectorConstruction(G4double factor, G4int verbose, G4bool isVisu) :
0049 G4VUserDetectorConstruction(), fFactor(factor), fVerbose(verbose), fBVisu(isVisu)
0050 {
0051 fDetectorMessenger = new DetectorConstructionMessenger(this);
0052 fWorldBoxSizeX = fWorldBoxSizeY = fWorldBoxSizeZ = 1*nm;
0053 if (gRunMode == RunningMode::Chem) fChemGeoImport = std::make_unique<ChemGeoImport>();
0054 }
0055
0056
0057
0058 G4VPhysicalVolume* DetectorConstruction::Construct()
0059
0060 {
0061 G4NistManager * man = G4NistManager::Instance();
0062 fWater = man->FindOrBuildMaterial("G4_WATER");
0063 if (gRunMode == RunningMode::Phys) ConstructFullCellNucleusGeo();
0064 else if (gRunMode == RunningMode::Chem) ConstructVoxelGeo();
0065 else {
0066
0067 G4ExceptionDescription msg;
0068 msg <<"Only world volume is constructed."
0069 <<" Make sure you choose the correct running mode."
0070 <<" Ignore this message if you intentionally test the code.";
0071 G4Exception("DetectorConstruction::Construct()",
0072 "", JustWarning, msg);
0073 fSolidWorld = new G4Box("solidWorld",fWorldBoxSizeX, fWorldBoxSizeY, fWorldBoxSizeZ);
0074 fLogicWorld = new G4LogicalVolume(fSolidWorld, fWater, "logicWorld");
0075 fPhysWorld = new G4PVPlacement(0,G4ThreeVector(),"physWorld",
0076 fLogicWorld, 0, false, false);
0077 }
0078 return fPhysWorld;
0079 }
0080
0081
0082
0083 G4VPhysicalVolume * DetectorConstruction::ConstructFullCellNucleusGeo()
0084 {
0085 PhysGeoImport geo(fBVisu);
0086 geo.SetFactor(fFactor);
0087
0088 std::map<G4String, G4LogicalVolume*> voxelMap;
0089
0090
0091 G4LogicalVolume* logicNucleus = geo.CreateNucleusLogicVolume(fCellDefFilePath);
0092 if (!fUsingUserDefinedSizesForWorld) {
0093 G4double scale = 2.0;
0094 fWorldBoxSizeX = scale*geo.GetNucleusSizeData()["SemiX"];
0095 fWorldBoxSizeY = scale*geo.GetNucleusSizeData()["SemiY"];
0096 fWorldBoxSizeZ = geo.GetNucleusSizeData()["SemiZ"];
0097 }
0098
0099 G4cout<<"===>>World box sizes: SemiX = "<<G4BestUnit(fWorldBoxSizeX,"Length")<<" , SemiY = "
0100 <<G4BestUnit(fWorldBoxSizeY,"Length")<<", SemiZ = "
0101 <<G4BestUnit(fWorldBoxSizeZ,"Length")<<"<<===="<<G4endl;
0102 fSolidWorld = new G4Box("solidWorld",fWorldBoxSizeX, fWorldBoxSizeY, fWorldBoxSizeZ);
0103 fLogicWorld = new G4LogicalVolume(fSolidWorld, fWater, "logicWorld");
0104 fPhysWorld = new G4PVPlacement(0,G4ThreeVector(),"physWorld", fLogicWorld, 0, false, false);
0105
0106 new G4PVPlacement(0, G4ThreeVector(), logicNucleus, "nucleus_pl", fLogicWorld, false, false);
0107
0108 G4LogicalVolume* logicStraightVoxel{nullptr}, *logicUpVoxel{nullptr}, *logicDownVoxel{nullptr},
0109 *logicLeftVoxel{nullptr},*logicRightVoxel{nullptr};
0110 G4LogicalVolume* logicStraightVoxel2{nullptr},*logicUpVoxel2{nullptr},*logicDownVoxel2{nullptr}
0111 , *logicLeftVoxel2{nullptr},*logicRightVoxel2{nullptr};
0112
0113 for (const auto & entry : fVoxelDefFilesList) {
0114 G4String name="";
0115 auto ptemp = geo.CreateLogicVolume(entry,name);
0116 if (name=="voxelStraight" || name=="VoxelStraight") {
0117 logicStraightVoxel = ptemp;
0118 voxelMap["VoxelStraight"] = logicStraightVoxel;
0119 }
0120 if (name=="voxelUp" || name=="VoxelUp") {
0121 logicUpVoxel = ptemp;
0122 voxelMap["VoxelUp"] = logicUpVoxel;
0123 }
0124 if (name=="voxelDown" || name=="VoxelDown") {
0125 logicDownVoxel = ptemp;
0126 voxelMap["VoxelDown"] = logicDownVoxel;
0127 }
0128 if (name=="voxelRight" || name=="VoxelRight") {
0129 logicRightVoxel = ptemp;
0130 voxelMap["VoxelRight"] = logicRightVoxel;
0131 }
0132 if (name=="voxelLeft" || name=="VoxelLeft") {
0133 logicLeftVoxel = ptemp;
0134 voxelMap["VoxelLeft"] = logicLeftVoxel;
0135 }
0136 if (name=="voxelStraight2" || name=="VoxelStraight2") {
0137 logicStraightVoxel2 = ptemp;
0138 voxelMap["VoxelStraight2"] = logicStraightVoxel2;
0139 }
0140 if (name=="voxelUp2" || name=="VoxelUp2") {
0141 logicUpVoxel2 = ptemp;
0142 voxelMap["VoxelUp2"] = logicUpVoxel2;
0143 }
0144 if (name=="voxelDown2" || name=="VoxelDown2") {
0145 logicDownVoxel2 = ptemp;
0146 voxelMap["VoxelDown2"] = logicDownVoxel2;
0147 }
0148 if (name=="voxelRight2" || name=="VoxelRight2") {
0149 logicRightVoxel2 = ptemp;
0150 voxelMap["VoxelRight2"] = logicRightVoxel2;
0151 }
0152 if (name=="voxelLeft2" || name=="VoxelLeft2") {
0153 logicLeftVoxel2 = ptemp;
0154 voxelMap["VoxelLeft2"] = logicLeftVoxel2;
0155 }
0156 }
0157
0158
0159 std::vector<Voxel>* voxelTable = geo.CreateVoxelsData(fCellDefFilePath);
0160 const G4int nucleusSize = voxelTable->size();
0161 G4cout<<"============================================================================"<<G4endl;
0162 G4cout<<"=====> Number of Histones in each voxel: "<<G4endl;
0163 for (auto [key, value] : geo.GetVoxelNbHistoneMap()) {
0164 G4cout<<key<<": \t\t\t"<<value<<G4endl;
0165 }
0166 G4cout<<"=====> Number of Basepairs in each voxel: "<<G4endl;
0167 for (auto [key, value] : geo.GetVoxelNbBpMap()) {
0168 G4cout<<key<<": \t\t\t"<<value<<G4endl;
0169 }
0170 G4cout <<"=====> Total Number of Histones placed in geometry: \t"
0171 <<geo.GetTotalNbHistonePlacedInGeo()<<G4endl;
0172 G4cout <<"=====> Total Number of Basepairs placed in geometry: \t"
0173 <<geo.GetTotalNbBpPlacedInGeo()<<G4endl;
0174 G4cout <<"=====> Number of each chromatin type placed in geometry: "<<G4endl;
0175 for (auto [key, value] : geo.GetChromatinTypeCountMap()) {
0176 G4String chromatinname = "Heterochromatin";
0177 if (key == ChromatinType::fEuchromatin) chromatinname = "Euchromatin";
0178 G4cout<<chromatinname<<": \t\t\t"<<value<<G4endl;
0179 }
0180 G4cout<<"============================================================================"<<G4endl;
0181
0182
0183 if (voxelMap.size() > 0) {
0184
0185
0186 G4double prevZpos = -fWorldBoxSizeZ;
0187 for (auto it = voxelMap.begin(); it != voxelMap.end(); it++) {
0188 G4double posX = fWorldBoxSizeX - geo.GetVoxelFullSize();
0189 G4double posY = fWorldBoxSizeY - geo.GetVoxelFullSize();
0190 G4double posZ = prevZpos + geo.GetVoxelFullSize();
0191 new G4PVPlacement(0, G4ThreeVector(posX,posY,posZ),
0192 it->second, "xx", fLogicWorld, false, 0);
0193 }
0194
0195 G4int nthreads=0;
0196 #ifdef G4MULTITHREADED
0197 nthreads = G4RunManager::GetRunManager()->GetNumberOfThreads();
0198 #endif
0199 if ( (nthreads >1) && (voxelMap.size() >1)) {
0200 G4ExceptionDescription msg;
0201 msg <<"Number of thread is "<<nthreads<<" > 1; Thus, dsbandrepair will run in testing mode. "
0202 <<"\nThere will be no DNA constituents placed inside Cell Nucleus!!!";
0203 G4Exception("DetectorConstruction::ConstructFullCellNucleusGeo",
0204 "RunningMode", JustWarning, msg);
0205 } else {
0206 G4VPVParameterisation* voxelParam = new VoxelParameterisation(voxelMap, voxelTable);
0207 new G4PVParameterised("VoxelParam", voxelMap.begin()->second,
0208 logicNucleus, kUndefined, nucleusSize, voxelParam);
0209 }
0210 } else {
0211 G4ExceptionDescription msg;
0212 msg <<"It seems that voxel-definition files are not provided."
0213 <<" There will be no DNA constituents placed inside Cell Nucleus!!!";
0214 G4Exception("DetectorConstruction::ConstructFullCellNucleusGeo",
0215 "Geo_InputFile_NoFile", JustWarning, msg);
0216 }
0217 Analysis::GetAnalysis()->RecordVoxelDefFilesList(fVoxelDefFilesList);
0218 Analysis::GetAnalysis()->SetTotalNbBpPlacedInGeo(geo.GetTotalNbBpPlacedInGeo());
0219 Analysis::GetAnalysis()->SetTotalNbHistonePlacedInGeo(geo.GetTotalNbHistonePlacedInGeo());
0220 Analysis::GetAnalysis()->SetNucleusVolume(geo.GetNucleusVolume());
0221 Analysis::GetAnalysis()->SetNucleusMassDensity(
0222 logicNucleus->GetMaterial()->GetDensity()/(kg/m3));
0223 return fPhysWorld;
0224 }
0225
0226
0227
0228 G4VPhysicalVolume *DetectorConstruction::ConstructVoxelGeo()
0229 {
0230 if (fWorldBoxSizeX < fVoxelHalfSizeXYZ) {
0231 fWorldBoxSizeX = fWorldBoxSizeY = fWorldBoxSizeZ = fVoxelHalfSizeXYZ*1.1;
0232 }
0233 fSolidWorld = new G4Box("solidWorld", fWorldBoxSizeX, fWorldBoxSizeY, fWorldBoxSizeZ);
0234 fLogicWorld = new G4LogicalVolume(fSolidWorld, fWater, "logicWorld");
0235 fPhysWorld = new G4PVPlacement(0,G4ThreeVector(),"physWorld", fLogicWorld, 0, false, 0);
0236
0237 if (fVoxelHalfSizeXYZ>0) {
0238 G4cout<<"=====> Construct voxel ........"<<G4endl;
0239 G4Box* solidVoxel = new G4Box("solidVoxel", fVoxelHalfSizeXYZ, fVoxelHalfSizeXYZ, fVoxelHalfSizeXYZ );
0240 G4LogicalVolume* logicVoxel = new G4LogicalVolume(solidVoxel, fWater, "logicVoxel");
0241 new G4PVPlacement(0, G4ThreeVector(), logicVoxel, "physVoxel", fLogicWorld, false, 0);
0242 }
0243 return fPhysWorld;
0244 }
0245
0246
0247
0248 void DetectorConstruction::SetCellDefFilePath(const G4String finput)
0249 {
0250 const G4fs::path thisP{std::string(finput)};
0251 G4fs::file_status fst = G4fs::file_status{};
0252 auto isExist = G4fs::status_known(fst) ? G4fs::exists(fst) : G4fs::exists(thisP);
0253 if (! isExist) {
0254 G4String msg = "File " + finput + "does not exist !!! ";
0255 G4Exception("DetectorConstruction::SetNucleusDefFilePath()",
0256 "Geo_InputFileNotOpened", FatalException, msg);
0257 } else {
0258 fCellDefFilePath = finput;
0259 Analysis::GetAnalysis()->RecordCellDefFiliePath(fCellDefFilePath);
0260 }
0261 }
0262
0263
0264
0265 void DetectorConstruction::AddVoxelDefFile(const G4String finput)
0266 {
0267 const G4fs::path thisP{std::string(finput)};
0268 G4fs::file_status fst = G4fs::file_status{};
0269 auto isExist = G4fs::status_known(fst) ? G4fs::exists(fst) : G4fs::exists(thisP);
0270 if (! isExist) {
0271 G4String msg = "File " + finput + "does not exist !!! ";
0272 G4Exception("DetectorConstruction::AddVoxelDefFile()",
0273 "Geo_InputFileNotOpened", FatalException, msg);
0274 } else {
0275 fVoxelDefFilesList.insert(finput);
0276 }
0277 }
0278
0279
0280
0281 void DetectorConstruction::SetWorldBoxSizes(G4ThreeVector v)
0282 {
0283 G4double sizex = v.getX();
0284 G4double sizey = v.getY();
0285 G4double sizez = v.getZ();
0286 if (sizex > 0. && sizey > 0. && sizez > 0.) {
0287 fWorldBoxSizeX = sizex;
0288 fWorldBoxSizeY = sizey;
0289 fWorldBoxSizeZ = sizez;
0290 fUsingUserDefinedSizesForWorld = true;
0291 } else {
0292 G4ExceptionDescription msg ;
0293 msg << " Check your setting? One world dimensions is <= 0 !!! ";
0294 G4Exception("DetectorConstruction::SetWorldBoxSizes",
0295 "Geo_WorldSizes", FatalException, msg);
0296 }
0297 }
0298
0299
0300
0301 void DetectorConstruction::ParseGeoFileForChemMode(const G4String fn)
0302 {
0303 fChemGeoImport->ParseFiles(fn);
0304 fVoxelHalfSizeXYZ = fChemGeoImport->GetSize()/2.;
0305 }
