dsbandrepair/src/DetectorConstruction.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 /// \file DetectorConstruction.cc
0028 /// \brief Implementation of the DetectorConstruction class
0029
0030 #include "DetectorConstruction.hh"
0031 #include "InformationKeeper.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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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         // only a world volume
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0082
0083 G4VPhysicalVolume * DetectorConstruction::ConstructFullCellNucleusGeo()
0084 {
0085     PhysGeoImport geo(fBVisu);
0086     geo.SetFactor(fFactor);
0087
0088     std::map<G4String, G4LogicalVolume*> voxelMap;
0089     // Create an empty logical nucleus
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     // then place cell nucleus in world
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     // Create the voxel data table
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     // Create the voxel parameterisation
0183     if (voxelMap.size() > 0) {
0184         // The following dummy declarations are nescessary for SetLogicalVolum()
0185         // in VoxelParameterisation to prevent from coredump
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         // End dummy declaration
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     InformationKeeper::Instance()->RecordVoxelDefFilesList(fVoxelDefFilesList);
0218     InformationKeeper::Instance()->SetTotalNbBpPlacedInGeo(geo.GetTotalNbBpPlacedInGeo());
0219     InformationKeeper::Instance()->SetTotalNbHistonePlacedInGeo(geo.GetTotalNbHistonePlacedInGeo());
0220     InformationKeeper::Instance()->SetNucleusVolume(geo.GetNucleusVolume());
0221     InformationKeeper::Instance()->SetNucleusMassDensity(
0222         logicNucleus->GetMaterial()->GetDensity()/(kg/m3)); // density in kg/m3;
0223     return fPhysWorld;
0224 }
0225
0226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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         InformationKeeper::Instance()->RecordCellDefFiliePath(fCellDefFilePath);
0260     }
0261 }
0262
0263 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0300
0301 void DetectorConstruction::ParseGeoFileForChemMode(const G4String fn)
0302 {
0303     fChemGeoImport->ParseFiles(fn);
0304     fVoxelHalfSizeXYZ = fChemGeoImport->GetSize()/2.;
0305 }