| 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 DetectorConstruction.cc
28 /// \brief Implementation of the DetectorConstruction class
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........oooOO0OOooo........oooOO0OOooo......
47
48 DetectorConstruction::DetectorConstruction(G4double factor, G4int verbose, G4bool isVisu) :
49 G4VUserDetectorConstruction(), fFactor(factor), fVerbose(verbose), fBVisu(isVisu)
50 {
51 fDetectorMessenger = new DetectorConstructionMessenger(this);
52 fWorldBoxSizeX = fWorldBoxSizeY = fWorldBoxSizeZ = 1*nm;
53 if (gRunMode == RunningMode::Chem) fChemGeoImport = std::make_unique<ChemGeoImport>();
54 }
55
56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57
58 G4VPhysicalVolume* DetectorConstruction::Construct()
59
60 {
61 G4NistManager * man = G4NistManager::Instance();
62 fWater = man->FindOrBuildMaterial("G4_WATER");
63 if (gRunMode == RunningMode::Phys) ConstructFullCellNucleusGeo();
64 else if (gRunMode == RunningMode::Chem) ConstructVoxelGeo();
65 else {
66 // only a world volume
67 G4ExceptionDescription msg;
68 msg <<"Only world volume is constructed."
69 <<" Make sure you choose the correct running mode."
70 <<" Ignore this message if you intentionally test the code.";
71 G4Exception("DetectorConstruction::Construct()",
72 "", JustWarning, msg);
73 fSolidWorld = new G4Box("solidWorld",fWorldBoxSizeX, fWorldBoxSizeY, fWorldBoxSizeZ);
74 fLogicWorld = new G4LogicalVolume(fSolidWorld, fWater, "logicWorld");
75 fPhysWorld = new G4PVPlacement(0,G4ThreeVector(),"physWorld",
76 fLogicWorld, 0, false, false);
77 }
78 return fPhysWorld;
79 }
80
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82
83 G4VPhysicalVolume * DetectorConstruction::ConstructFullCellNucleusGeo()
84 {
85 PhysGeoImport geo(fBVisu);
86 geo.SetFactor(fFactor);
87
88 std::map<G4String, G4LogicalVolume*> voxelMap;
89 // Create an empty logical nucleus
90
91 G4LogicalVolume* logicNucleus = geo.CreateNucleusLogicVolume(fCellDefFilePath);
92 if (!fUsingUserDefinedSizesForWorld) {
93 G4double scale = 2.0;
94 fWorldBoxSizeX = scale*geo.GetNucleusSizeData()["SemiX"];
95 fWorldBoxSizeY = scale*geo.GetNucleusSizeData()["SemiY"];
96 fWorldBoxSizeZ = geo.GetNucleusSizeData()["SemiZ"];
97 }
98
99 G4cout<<"===>>World box sizes: SemiX = "<<G4BestUnit(fWorldBoxSizeX,"Length")<<" , SemiY = "
100 <<G4BestUnit(fWorldBoxSizeY,"Length")<<", SemiZ = "
101 <<G4BestUnit(fWorldBoxSizeZ,"Length")<<"<<===="<<G4endl;
102 fSolidWorld = new G4Box("solidWorld",fWorldBoxSizeX, fWorldBoxSizeY, fWorldBoxSizeZ);
103 fLogicWorld = new G4LogicalVolume(fSolidWorld, fWater, "logicWorld");
104 fPhysWorld = new G4PVPlacement(0,G4ThreeVector(),"physWorld", fLogicWorld, 0, false, false);
105 // then place cell nucleus in world
106 new G4PVPlacement(0, G4ThreeVector(), logicNucleus, "nucleus_pl", fLogicWorld, false, false);
107
108 G4LogicalVolume* logicStraightVoxel{nullptr}, *logicUpVoxel{nullptr}, *logicDownVoxel{nullptr},
109 *logicLeftVoxel{nullptr},*logicRightVoxel{nullptr};
110 G4LogicalVolume* logicStraightVoxel2{nullptr},*logicUpVoxel2{nullptr},*logicDownVoxel2{nullptr}
111 , *logicLeftVoxel2{nullptr},*logicRightVoxel2{nullptr};
112
113 for (const auto & entry : fVoxelDefFilesList) {
114 G4String name="";
115 auto ptemp = geo.CreateLogicVolume(entry,name);
116 if (name=="voxelStraight" || name=="VoxelStraight") {
117 logicStraightVoxel = ptemp;
118 voxelMap["VoxelStraight"] = logicStraightVoxel;
119 }
120 if (name=="voxelUp" || name=="VoxelUp") {
121 logicUpVoxel = ptemp;
122 voxelMap["VoxelUp"] = logicUpVoxel;
123 }
124 if (name=="voxelDown" || name=="VoxelDown") {
125 logicDownVoxel = ptemp;
126 voxelMap["VoxelDown"] = logicDownVoxel;
127 }
128 if (name=="voxelRight" || name=="VoxelRight") {
129 logicRightVoxel = ptemp;
130 voxelMap["VoxelRight"] = logicRightVoxel;
131 }
132 if (name=="voxelLeft" || name=="VoxelLeft") {
133 logicLeftVoxel = ptemp;
134 voxelMap["VoxelLeft"] = logicLeftVoxel;
135 }
136 if (name=="voxelStraight2" || name=="VoxelStraight2") {
137 logicStraightVoxel2 = ptemp;
138 voxelMap["VoxelStraight2"] = logicStraightVoxel2;
139 }
140 if (name=="voxelUp2" || name=="VoxelUp2") {
141 logicUpVoxel2 = ptemp;
142 voxelMap["VoxelUp2"] = logicUpVoxel2;
143 }
144 if (name=="voxelDown2" || name=="VoxelDown2") {
145 logicDownVoxel2 = ptemp;
146 voxelMap["VoxelDown2"] = logicDownVoxel2;
147 }
148 if (name=="voxelRight2" || name=="VoxelRight2") {
149 logicRightVoxel2 = ptemp;
150 voxelMap["VoxelRight2"] = logicRightVoxel2;
151 }
152 if (name=="voxelLeft2" || name=="VoxelLeft2") {
153 logicLeftVoxel2 = ptemp;
154 voxelMap["VoxelLeft2"] = logicLeftVoxel2;
155 }
156 }
157
158 // Create the voxel data table
159 std::vector<Voxel>* voxelTable = geo.CreateVoxelsData(fCellDefFilePath);
160 const G4int nucleusSize = voxelTable->size();
161 G4cout<<"============================================================================"<<G4endl;
162 G4cout<<"=====> Number of Histones in each voxel: "<<G4endl;
163 for (auto [key, value] : geo.GetVoxelNbHistoneMap()) {
164 G4cout<<key<<": \t\t\t"<<value<<G4endl;
165 }
166 G4cout<<"=====> Number of Basepairs in each voxel: "<<G4endl;
167 for (auto [key, value] : geo.GetVoxelNbBpMap()) {
168 G4cout<<key<<": \t\t\t"<<value<<G4endl;
169 }
170 G4cout <<"=====> Total Number of Histones placed in geometry: \t"
171 <<geo.GetTotalNbHistonePlacedInGeo()<<G4endl;
172 G4cout <<"=====> Total Number of Basepairs placed in geometry: \t"
173 <<geo.GetTotalNbBpPlacedInGeo()<<G4endl;
174 G4cout <<"=====> Number of each chromatin type placed in geometry: "<<G4endl;
175 for (auto [key, value] : geo.GetChromatinTypeCountMap()) {
176 G4String chromatinname = "Heterochromatin";
177 if (key == ChromatinType::fEuchromatin) chromatinname = "Euchromatin";
178 G4cout<<chromatinname<<": \t\t\t"<<value<<G4endl;
179 }
180 G4cout<<"============================================================================"<<G4endl;
181
182 // Create the voxel parameterisation
183 if (voxelMap.size() > 0) {
184 // The following dummy declarations are nescessary for SetLogicalVolum()
185 // in VoxelParameterisation to prevent from coredump
186 G4double prevZpos = -fWorldBoxSizeZ;
187 for (auto it = voxelMap.begin(); it != voxelMap.end(); it++) {
188 G4double posX = fWorldBoxSizeX - geo.GetVoxelFullSize();
189 G4double posY = fWorldBoxSizeY - geo.GetVoxelFullSize();
190 G4double posZ = prevZpos + geo.GetVoxelFullSize();
191 new G4PVPlacement(0, G4ThreeVector(posX,posY,posZ),
192 it->second, "xx", fLogicWorld, false, 0);
193 }
194 // End dummy declaration
195 G4int nthreads=0;
196 #ifdef G4MULTITHREADED
197 nthreads = G4RunManager::GetRunManager()->GetNumberOfThreads();
198 #endif
199 if ( (nthreads >1) && (voxelMap.size() >1)) {
200 G4ExceptionDescription msg;
201 msg <<"Number of thread is "<<nthreads<<" > 1; Thus, dsbandrepair will run in testing mode. "
202 <<"\nThere will be no DNA constituents placed inside Cell Nucleus!!!";
203 G4Exception("DetectorConstruction::ConstructFullCellNucleusGeo",
204 "RunningMode", JustWarning, msg);
205 } else {
206 G4VPVParameterisation* voxelParam = new VoxelParameterisation(voxelMap, voxelTable);
207 new G4PVParameterised("VoxelParam", voxelMap.begin()->second,
208 logicNucleus, kUndefined, nucleusSize, voxelParam);
209 }
210 } else {
211 G4ExceptionDescription msg;
212 msg <<"It seems that voxel-definition files are not provided."
213 <<" There will be no DNA constituents placed inside Cell Nucleus!!!";
214 G4Exception("DetectorConstruction::ConstructFullCellNucleusGeo",
215 "Geo_InputFile_NoFile", JustWarning, msg);
216 }
217 Analysis::GetAnalysis()->RecordVoxelDefFilesList(fVoxelDefFilesList);
218 Analysis::GetAnalysis()->SetTotalNbBpPlacedInGeo(geo.GetTotalNbBpPlacedInGeo());
219 Analysis::GetAnalysis()->SetTotalNbHistonePlacedInGeo(geo.GetTotalNbHistonePlacedInGeo());
220 Analysis::GetAnalysis()->SetNucleusVolume(geo.GetNucleusVolume());
221 Analysis::GetAnalysis()->SetNucleusMassDensity(
222 logicNucleus->GetMaterial()->GetDensity()/(kg/m3)); // density in kg/m3;
223 return fPhysWorld;
224 }
225
226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
227
228 G4VPhysicalVolume *DetectorConstruction::ConstructVoxelGeo()
229 {
230 if (fWorldBoxSizeX < fVoxelHalfSizeXYZ) {
231 fWorldBoxSizeX = fWorldBoxSizeY = fWorldBoxSizeZ = fVoxelHalfSizeXYZ*1.1;
232 }
233 fSolidWorld = new G4Box("solidWorld", fWorldBoxSizeX, fWorldBoxSizeY, fWorldBoxSizeZ);
234 fLogicWorld = new G4LogicalVolume(fSolidWorld, fWater, "logicWorld");
235 fPhysWorld = new G4PVPlacement(0,G4ThreeVector(),"physWorld", fLogicWorld, 0, false, 0);
236
237 if (fVoxelHalfSizeXYZ>0) {
238 G4cout<<"=====> Construct voxel ........"<<G4endl;
239 G4Box* solidVoxel = new G4Box("solidVoxel", fVoxelHalfSizeXYZ, fVoxelHalfSizeXYZ, fVoxelHalfSizeXYZ );
240 G4LogicalVolume* logicVoxel = new G4LogicalVolume(solidVoxel, fWater, "logicVoxel");
241 new G4PVPlacement(0, G4ThreeVector(), logicVoxel, "physVoxel", fLogicWorld, false, 0);
242 }
243 return fPhysWorld;
244 }
245
246 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
247
248 void DetectorConstruction::SetCellDefFilePath(const G4String finput)
249 {
250 const G4fs::path thisP{std::string(finput)};
251 G4fs::file_status fst = G4fs::file_status{};
252 auto isExist = G4fs::status_known(fst) ? G4fs::exists(fst) : G4fs::exists(thisP);
253 if (! isExist) {
254 G4String msg = "File " + finput + "does not exist !!! ";
255 G4Exception("DetectorConstruction::SetNucleusDefFilePath()",
256 "Geo_InputFileNotOpened", FatalException, msg);
257 } else {
258 fCellDefFilePath = finput;
259 Analysis::GetAnalysis()->RecordCellDefFiliePath(fCellDefFilePath);
260 }
261 }
262
263 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
264
265 void DetectorConstruction::AddVoxelDefFile(const G4String finput)
266 {
267 const G4fs::path thisP{std::string(finput)};
268 G4fs::file_status fst = G4fs::file_status{};
269 auto isExist = G4fs::status_known(fst) ? G4fs::exists(fst) : G4fs::exists(thisP);
270 if (! isExist) {
271 G4String msg = "File " + finput + "does not exist !!! ";
272 G4Exception("DetectorConstruction::AddVoxelDefFile()",
273 "Geo_InputFileNotOpened", FatalException, msg);
274 } else {
275 fVoxelDefFilesList.insert(finput);
276 }
277 }
278
279 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
280
281 void DetectorConstruction::SetWorldBoxSizes(G4ThreeVector v)
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 dimensions is <= 0 !!! ";
294 G4Exception("DetectorConstruction::SetWorldBoxSizes",
295 "Geo_WorldSizes", FatalException, msg);
296 }
297 }
298
299 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
300
301 void DetectorConstruction::ParseGeoFileForChemMode(const G4String fn)
302 {
303 fChemGeoImport->ParseFiles(fn);
304 fVoxelHalfSizeXYZ = fChemGeoImport->GetSize()/2.;
305 }