| 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 }