| 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 // This example is provided by the Geant4-DNA collaboration
27 // Any report or published results obtained using the Geant4-DNA software
28 // shall cite the following Geant4-DNA collaboration publication:
29 // Med. Phys. 37 (2010) 4692-4708
30 // and papers
31 // M. Batmunkh et al. J Radiat Res Appl Sci 8 (2015) 498-507
32 // O. Belov et al. Physica Medica 32 (2016) 1510-1520
33 // The Geant4-DNA web site is available at http://geant4-dna.org
34 //
35 // -------------------------------------------------------------------
36 // November 2016
37 // -------------------------------------------------------------------
38 //
39 /// \file DetectorConstruction.cc
40 /// \brief Implementation of the DetectorConstruction class
41
42 #include "DetectorConstruction.hh"
43
44 #include "G4Colour.hh"
45 #include "G4PhysicalConstants.hh"
46 #include "G4ProductionCuts.hh"
47 #include "G4Region.hh"
48 #include "G4RotationMatrix.hh"
49 #include "G4SystemOfUnits.hh"
50 #include "G4VisAttributes.hh"
51
52 #include <algorithm>
53 #include <iostream>
54
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56
57 DetectorConstruction::DetectorConstruction()
58 {
59 fNeuronLoadParamz = new NeuronLoadDataFile();
60 DefineMaterials();
61 }
62
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64
65 DetectorConstruction::~DetectorConstruction()
66 {
67 delete fNeuronLoadParamz;
68 }
69
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71
72 void DetectorConstruction::DefineMaterials()
73 {
74 // Water is defined from NIST material database
75 G4NistManager* man = G4NistManager::Instance();
76 fpWaterMaterial = man->FindOrBuildMaterial("G4_WATER");
77 fpWorldMaterial = man->FindOrBuildMaterial("G4_Galactic");
78 }
79
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
81
82 G4VPhysicalVolume* DetectorConstruction::Construct()
83 {
84 G4cout << " ---- Begin of Neuron Construction! ---- "
85 << "\n"
86 << " ==========================================================" << G4endl;
87
88 // ===============================================
89 // WORLD VOLUME - filled by default material
90 // ===============================================
91
92 // Dimensions of world volume are calculated as overall dimensions of neuron!
93
94 G4double worldSizeX;
95 worldSizeX = fNeuronLoadParamz->GetdiagnlLength() * um;
96
97 if (worldSizeX <= 0.0) {
98 worldSizeX = 10. * cm;
99 }
100
101 G4double worldSizeY = worldSizeX;
102 G4double worldSizeZ = worldSizeX;
103 G4cout << " Side length of word volume is calculated : " << worldSizeX / um << " um" << G4endl;
104 G4VSolid* worldS = new G4Box("World", worldSizeX / 2, worldSizeY / 2, worldSizeZ / 2);
105
106 G4LogicalVolume* worldLV = new G4LogicalVolume(worldS, fpWorldMaterial, "World");
107
108 // Visualization attributes
109 G4VisAttributes* worldVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5, 0.1)); // Gray
110 worldVisAtt->SetVisibility(true);
111 worldLV->SetVisAttributes(worldVisAtt);
112
113 G4VPhysicalVolume* worldPV = new G4PVPlacement(0, // no rotation
114 G4ThreeVector(), // at (0,0,0)
115 "World", // its name
116 worldLV, // its logical volume
117 0, // its mother volume
118 false, // no boolean operation
119 0, // copy number
120 true); // checking overlaps forced to YES
121
122 // ===============================================
123 // HOMOGENEOUS MEDIUM - LARGE SPHERE VOLUME
124 // ===============================================
125
126 // Radius of water sphere calculated as overall dimensions of neuron.
127 fTotMassMedium = fNeuronLoadParamz->GetTotMassMedium();
128 fTotSurfMedium = fNeuronLoadParamz->GetTotSurfMedium();
129 G4double RadiusMedium = fNeuronLoadParamz->GetdiagnlLength() * um / 2.;
130 G4cout << " Radius of homogeneous medium is calculated : " << RadiusMedium / um << " um"
131 << G4endl;
132 G4VSolid* mediumS = new G4Orb("Medium", RadiusMedium);
133
134 G4LogicalVolume* mediumLV = new G4LogicalVolume(mediumS, fpWaterMaterial, "Medium");
135
136 // Visualization attributes
137 G4VisAttributes* mediumVisAtt = new G4VisAttributes(G4Colour(0.0, 1.0, 0.0, 0.02)); // Green
138 // mediumVisAtt->SetForceSolid(true);
139 // mediumVisAtt->SetForceWireframe (true);
140 mediumVisAtt->SetForceLineSegmentsPerCircle(180);
141 mediumVisAtt->SetVisibility(true);
142 mediumLV->SetVisAttributes(mediumVisAtt);
143
144 G4VPhysicalVolume* mediumPV =
145 new G4PVPlacement(0, G4ThreeVector(), "Medium", mediumLV, worldPV, false, 0, fCheckOverlaps);
146
147 // ===============================================
148 // TARGET - BOUNDING SLICE including NEURON
149 // ===============================================
150
151 // Dimensions of bounding slice volume defined as overall measure of neuron
152
153 G4double TargetSizeX = fNeuronLoadParamz->GetwidthB() * um;
154 G4double TargetSizeY = fNeuronLoadParamz->GetheightB() * um;
155 G4double TargetSizeZ = fNeuronLoadParamz->GetdepthB() * um;
156 fTotMassSlice = fNeuronLoadParamz->GetTotMassSlice();
157 G4cout << " Overall dimensions (um) of neuron morphology : "
158 << "\n"
159 << '\t' << " width = " << TargetSizeX / um << " height = " << TargetSizeY / um
160 << " depth = " << TargetSizeZ / um << G4endl;
161
162 G4cout << " Volume (um3), surface (um2) and mass (ug) of Bounding Slice are"
163 << " calculated : "
164 << "\n"
165 << '\t' << fNeuronLoadParamz->GetTotVolSlice() / std::pow(um, 3) << "; " << '\t'
166 << fNeuronLoadParamz->GetTotSurfSlice() / (um * um) << "; " << '\t'
167 << fNeuronLoadParamz->GetTotMassSlice() * 1e6 / g << "\n " << G4endl;
168
169 G4Box* boundingS =
170 new G4Box("BoundingSlice", TargetSizeX / 2., TargetSizeY / 2., TargetSizeZ / 2.);
171
172 G4LogicalVolume* boundingLV = new G4LogicalVolume(boundingS, fpWaterMaterial, "BoundingSlice");
173
174 // Visualization attributes with opacity!
175 G4VisAttributes* TargetVisAtt = new G4VisAttributes(G4Colour(1.0, 1.0, 1.0, 0.1));
176 TargetVisAtt->SetForceSolid(true);
177 TargetVisAtt->SetVisibility(true);
178 boundingLV->SetVisAttributes(TargetVisAtt);
179 new G4PVPlacement(0, G4ThreeVector(), "BoundingSlice", boundingLV, mediumPV, false, 0,
180 fCheckOverlaps);
181
182 // ===============================================
183 // NEURON MORPHOLOGY
184 // ===============================================
185
186 G4cout << " Volume (um3), surface (um2) and mass(ug) of Neuron "
187 << "are calculated : "
188 << "\n"
189 << '\t' << fNeuronLoadParamz->GetTotVolNeuron() / std::pow(um, 3) << "; " << '\t'
190 << fNeuronLoadParamz->GetTotSurfNeuron() / (um * um) << "; " << '\t'
191 << fNeuronLoadParamz->GetTotMassNeuron() * 1e6 / g << G4endl;
192 fTotMassNeuron = fNeuronLoadParamz->GetTotMassNeuron();
193 G4cout << " Total number of compartments into Neuron : " << fNeuronLoadParamz->GetnbNeuroncomp()
194 << G4endl;
195 G4cout << " Shift values (um) for Neuron translation are calculated : "
196 << "\n"
197 << '\t' << " shiftX = " << fNeuronLoadParamz->GetshiftX()
198 << " shiftY = " << fNeuronLoadParamz->GetshiftY()
199 << " shiftZ = " << fNeuronLoadParamz->GetshiftZ() << G4endl;
200
201 // Soma in Violet with opacity // 0.85,0.44,0.84
202 fSomaColour = new G4VisAttributes;
203 fSomaColour->SetColour(G4Colour(G4Colour(22 / 255., 200 / 255., 30 / 255.)));
204 fSomaColour->SetForceSolid(true); // true
205 fSomaColour->SetVisibility(true);
206
207 // Dendrites in Dark-Blue
208 fDendColour = new G4VisAttributes;
209 fDendColour->SetColour(G4Colour(G4Colour(0.0, 0.0, 0.5)));
210 fDendColour->SetForceSolid(true);
211 // fDendColour->SetVisibility(true);
212
213 // Axon in Maroon
214 fAxonColour = new G4VisAttributes;
215 fAxonColour->SetColour(G4Colour(G4Colour(0.5, 0.0, 0.0)));
216 fAxonColour->SetForceSolid(true);
217 fAxonColour->SetVisibility(true);
218
219 // Spines in Dark-Green
220 fSpineColour = new G4VisAttributes;
221 fSpineColour->SetColour(G4Colour(G4Colour(0.0, 100 / 255., 0.0)));
222 fSpineColour->SetForceSolid(true);
223 fSpineColour->SetVisibility(true);
224
225 // Whole neuron in semitransparent navy blue
226 fNeuronColour = new G4VisAttributes;
227 fNeuronColour->SetColour(G4Colour(G4Colour(0.0, 0.4, 0.8, 0.9)));
228 fNeuronColour->SetForceSolid(true);
229 fNeuronColour->SetVisibility(true);
230
231 // Placement volumes: G4examples/extended/parameterisations/gflash
232
233 // =======================================================================
234 // Structure-1: Soma
235
236 // Create Target G4Region and add logical volume
237 // Active Geant4-DNA processes in this region
238 fpRegion = new G4Region("Soma");
239 G4ProductionCuts* cuts = new G4ProductionCuts();
240 G4double defCut = 1 * nanometer;
241 cuts->SetProductionCut(defCut, "gamma");
242 cuts->SetProductionCut(defCut, "e-");
243 cuts->SetProductionCut(defCut, "e+");
244 cuts->SetProductionCut(defCut, "proton");
245 fpRegion->SetProductionCuts(cuts);
246 G4ThreeVector shift(fNeuronLoadParamz->GetshiftX(), fNeuronLoadParamz->GetshiftY(),
247 fNeuronLoadParamz->GetshiftZ());
248
249 G4int n = fNeuronLoadParamz->GetnbSomacomp();
250 if (n <= 0) {
251 G4cout << " ---- Soma not found! ---- " << G4endl;
252 }
253 else {
254 G4cout << " ---- Soma for construction: ---- " << G4endl;
255 G4cout << " Total number of compartments into Soma : " << n << G4endl;
256 fnbSomacomp = n;
257 fMassSomaTot = fNeuronLoadParamz->GetMassSomaTot();
258 fMassSomacomp.resize(n, 0.0);
259 fPosSomacomp.resize(n);
260 fsomaS.resize(n, nullptr);
261 fsomaLV.resize(n, nullptr);
262 fsomaPV.resize(n, nullptr);
263
264 for (G4int i = 0; i < n; ++i) {
265 fsomaS[i] = new G4Orb("Soma", fNeuronLoadParamz->GetRadSomacomp(i) * um);
266 // you can change parameters of Soma with a single compartment
267 G4ThreeVector pos = (fNeuronLoadParamz->GetPosSomacomp(i) - shift) * um;
268 fsomaLV[i] = new G4LogicalVolume(fsomaS[i], fpWaterMaterial, "Soma");
269 fsomaLV[i]->SetVisAttributes(fSomaColour);
270 fsomaPV[i] = new G4PVPlacement(0, // no rotation
271 pos, fsomaLV[i], "Soma", boundingLV, false, i);
272 fMassSomacomp[i] = fNeuronLoadParamz->GetMassSomacomp(i);
273 fPosSomacomp[i] = fNeuronLoadParamz->GetPosSomacomp(i);
274 fpRegion->AddRootLogicalVolume(fsomaLV[i]);
275 }
276 }
277
278 // =======================================================================
279 // Structure-2: Dendrites
280
281 n = fNeuronLoadParamz->GetnbDendritecomp();
282 if (n <= 0) {
283 G4cout << " ---- Dendrites not found! ---- " << G4endl;
284 }
285 else {
286 fnbDendritecomp = n;
287 G4cout << " ---- Dendrites for construction: ---- " << G4endl;
288 G4cout << " Total number of compartments into Dendrites: " << n << G4endl;
289
290 // Active Geant4-DNA processes in this region
291 auto region = new G4Region("Dendrites");
292 region->SetProductionCuts(cuts);
293
294 fMassDendTot = fNeuronLoadParamz->GetMassDendTot();
295 fMassDendcomp.resize(n, 0.0);
296 fDistADendSoma.resize(n, 0.0);
297 fDistBDendSoma.resize(n, 0.0);
298 fPosDendcomp.resize(n);
299 fdendriteS.resize(n, nullptr);
300 fdendriteLV.resize(n, nullptr);
301 fdendritePV.resize(n, nullptr);
302 for (G4int i = 1; i < n; ++i) {
303 fdendriteS[i] = new G4Tubs("Dendrites", 0., fNeuronLoadParamz->GetRadDendcomp(i) * um,
304 fNeuronLoadParamz->GetHeightDendcomp(i) * um / 2., 0., 2. * pi);
305 fdendriteLV[i] = new G4LogicalVolume(fdendriteS[i], fpWaterMaterial, "Dendrites");
306 fdendriteLV[i]->SetVisAttributes(fDendColour);
307
308 G4ThreeVector pos = (fNeuronLoadParamz->GetPosDendcomp(i) - shift) * um;
309 // rotation checking with function ComputeTransformation!
310 // RotationMatrix with Inverse
311 fdendritePV[i] = new G4PVPlacement(G4Transform3D(fNeuronLoadParamz->GetRotDendcomp(i), pos),
312 fdendriteLV[i], "Dendrites", boundingLV, false, i);
313 fMassDendcomp[i] = fNeuronLoadParamz->GetMassDendcomp(i);
314 fPosDendcomp[i] = fNeuronLoadParamz->GetPosDendcomp(i);
315 fDistADendSoma[i] = fNeuronLoadParamz->GetDistADendSoma(i);
316 fDistBDendSoma[i] = fNeuronLoadParamz->GetDistBDendSoma(i);
317 region->AddRootLogicalVolume(fdendriteLV[i]);
318 }
319 }
320
321 // =======================================================================
322 // Structure-3: Axon
323
324 n = fNeuronLoadParamz->GetnbAxoncomp();
325 if (n <= 0) {
326 G4cout << " ---- Axon not found! ---- " << G4endl;
327 }
328 else {
329 G4cout << " ---- Axon for construction: ---- " << G4endl;
330 G4cout << " Total number of compartments into Axon : " << n << G4endl;
331 // Active Geant4-DNA processes in this region
332 auto region = new G4Region("Axon");
333 region->SetProductionCuts(cuts);
334
335 fnbAxoncomp = n;
336 fMassAxonTot = fNeuronLoadParamz->GetMassAxonTot();
337 fMassAxoncomp.resize(n, 0.0);
338 fDistAxonsoma.resize(n, 0.0);
339 fPosAxoncomp.resize(n);
340 faxonS.resize(n, nullptr);
341 faxonLV.resize(n, nullptr);
342 faxonPV.resize(n, nullptr);
343
344 for (G4int i = 1; i < n; ++i) {
345 faxonS[i] = new G4Tubs("Axon", 0., fNeuronLoadParamz->GetRadAxoncomp(i) * um,
346 fNeuronLoadParamz->GetHeightAxoncomp(i) * um / 2., 0., 2. * pi);
347 faxonLV[i] = new G4LogicalVolume(faxonS[i], fpWaterMaterial, "Axon");
348 faxonLV[i]->SetVisAttributes(fAxonColour);
349 // RotationMatrix with Inverse
350 G4ThreeVector pos = (fNeuronLoadParamz->GetPosAxoncomp(i) - shift) * um;
351 faxonPV[i] = new G4PVPlacement(G4Transform3D(fNeuronLoadParamz->GetRotAxoncomp(i), pos),
352 faxonLV[i], "Axon", boundingLV, false, i);
353 fMassAxoncomp[i] = fNeuronLoadParamz->GetMassAxoncomp(i);
354 fPosAxoncomp[i] = fNeuronLoadParamz->GetPosAxoncomp(i);
355 fDistAxonsoma[i] = fNeuronLoadParamz->GetDistAxonsoma(i);
356 region->AddRootLogicalVolume(faxonLV[i]);
357 }
358 }
359 // =======================================================================
360 // Structure-4: Spines
361 if (fNeuronLoadParamz->GetnbSpinecomp() == 0) {
362 G4cout << " ---- Spines not found! ---- " << G4endl;
363 }
364 else {
365 G4cout << " ---- Spines for construction: ---- " << G4endl;
366 G4cout << " Total number of compartments into Spines : " << fNeuronLoadParamz->GetnbSpinecomp()
367 << G4endl;
368 }
369
370 G4cout << "\n ---- End of Neuron Construction! ---- "
371 << "\n ========================================================== \n"
372 << G4endl;
373
374 // =======================================================================
375 // Active Geant4-DNA processes in BoundingSlice with whole neuron structure
376 // fpRegion = new G4Region("BoundingSlice");
377 // fpRegion->SetProductionCuts(cuts);
378 // fpRegion->AddRootLogicalVolume(boundingLV);
379
380 return worldPV;
381 }
382