| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // This example is provided by the Geant4-DNA 26 // This example is provided by the Geant4-DNA collaboration
27 // Any report or published results obtained us 27 // Any report or published results obtained using the Geant4-DNA software
28 // shall cite the following Geant4-DNA collabo 28 // shall cite the following Geant4-DNA collaboration publication:
29 // Med. Phys. 37 (2010) 4692-4708 29 // Med. Phys. 37 (2010) 4692-4708
30 // and papers 30 // and papers
31 // M. Batmunkh et al. J Radiat Res Appl Sci 8 31 // M. Batmunkh et al. J Radiat Res Appl Sci 8 (2015) 498-507
32 // O. Belov et al. Physica Medica 32 (2016) 15 << 32 // O. Belov et al. Physica Medica 32 (2016) 1510-1520
33 // The Geant4-DNA web site is available at htt 33 // The Geant4-DNA web site is available at http://geant4-dna.org
34 // << 34 //
35 // ------------------------------------------- 35 // -------------------------------------------------------------------
36 // November 2016 36 // November 2016
37 // ------------------------------------------- 37 // -------------------------------------------------------------------
38 // 38 //
39 /// \file DetectorConstruction.cc << 39 // $ID$
>> 40 /// \file DetectorConstruction.cc
40 /// \brief Implementation of the DetectorConst 41 /// \brief Implementation of the DetectorConstruction class
41 42
42 #include "DetectorConstruction.hh" 43 #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" 44 #include "G4SystemOfUnits.hh"
>> 45 #include "G4Region.hh"
>> 46 #include "G4ProductionCuts.hh"
>> 47 #include "G4PhysicalConstants.hh"
>> 48 #include "G4Colour.hh"
50 #include "G4VisAttributes.hh" 49 #include "G4VisAttributes.hh"
51 << 50 #include "G4RotationMatrix.hh"
52 #include <algorithm> << 51 #include <algorithm>
53 #include <iostream> 52 #include <iostream>
54 53
55 //....oooOO0OOooo........oooOO0OOooo........oo 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 55
57 DetectorConstruction::DetectorConstruction() 56 DetectorConstruction::DetectorConstruction()
>> 57 :G4VUserDetectorConstruction(),
>> 58 fpDefaultMaterial(0), fpWaterMaterial(0),
>> 59 fpRegion(0), fCheckOverlaps(false)
58 { 60 {
59 fNeuronLoadParamz = new NeuronLoadDataFile() << 61
60 DefineMaterials(); << 62 }
61 } <<
62 63
63 //....oooOO0OOooo........oooOO0OOooo........oo 64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64 65
65 DetectorConstruction::~DetectorConstruction() 66 DetectorConstruction::~DetectorConstruction()
66 { 67 {
67 delete fNeuronLoadParamz; << 68
68 } 69 }
69 70
70 //....oooOO0OOooo........oooOO0OOooo........oo 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71 72
72 void DetectorConstruction::DefineMaterials() << 73 G4VPhysicalVolume* DetectorConstruction::Construct()
73 { 74 {
>> 75
>> 76 // load Neuron positions and obtain parameters!
>> 77 fNeuronLoadParamz = new NeuronLoadDataFile() ;
>> 78
>> 79 DefineMaterials();
>> 80 return ConstructDetector();
>> 81
>> 82 }
>> 83
>> 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 85
>> 86 void DetectorConstruction::DefineMaterials()
>> 87 {
74 // Water is defined from NIST material datab 88 // Water is defined from NIST material database
75 G4NistManager* man = G4NistManager::Instance << 89 G4NistManager * man = G4NistManager::Instance();
76 fpWaterMaterial = man->FindOrBuildMaterial(" << 90 G4Material * H2O = man->FindOrBuildMaterial("G4_WATER");
77 fpWorldMaterial = man->FindOrBuildMaterial(" << 91 // Default materials in setup.
>> 92 fpWaterMaterial = H2O;
>> 93 G4Material * Vacuum = man->FindOrBuildMaterial("G4_Galactic");
>> 94 //G4Material * Air = man->FindOrBuildMaterial("G4_AIR");
>> 95 fpDefaultMaterial = Vacuum;
78 } 96 }
79 97
80 //....oooOO0OOooo........oooOO0OOooo........oo 98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
81 99
82 G4VPhysicalVolume* DetectorConstruction::Const << 100 G4VPhysicalVolume* DetectorConstruction::ConstructDetector()
83 { 101 {
84 G4cout << " ---- Begin of Neuron Constructio << 102 G4cout <<" ---- Begin of Neuron Construction! ---- "
85 << "\n" << 103 <<"\n"<<" =========================================================="<<G4endl;
86 << " ================================ <<
87 104
88 // ========================================= << 105 // ===============================================
89 // WORLD VOLUME - filled by default material 106 // WORLD VOLUME - filled by default material
90 // ========================================= 107 // ===============================================
91 108
92 // Dimensions of world volume are calculated << 109 // Dimensions of world volume are calculated as overall dimensions of neuron!
93 << 110
94 G4double worldSizeX; << 111 G4double worldSizeX;
95 worldSizeX = fNeuronLoadParamz->GetdiagnlLen << 112 worldSizeX = 1.0*fNeuronLoadParamz->GetdiagnlLength()*um;
96 << 113
97 if (worldSizeX <= 0.0) { << 114 if (!worldSizeX)
98 worldSizeX = 10. * cm; << 115 {
99 } << 116 worldSizeX = 10.*cm;
100 << 117 }
101 G4double worldSizeY = worldSizeX; << 118
102 G4double worldSizeZ = worldSizeX; << 119 G4double worldSizeY = worldSizeX;
103 G4cout << " Side length of word volume is ca << 120 G4double worldSizeZ = worldSizeX;
104 G4VSolid* worldS = new G4Box("World", worldS << 121 G4cout << " Side length of word volume is calculated : "
105 << 122 << worldSizeX/um <<" um"<< G4endl;
106 G4LogicalVolume* worldLV = new G4LogicalVolu << 123 G4VSolid* worldS = new G4Box("World", //its name
>> 124 worldSizeX/2, worldSizeY/2, worldSizeZ/2); //its size
>> 125
>> 126 G4LogicalVolume* worldLV = new G4LogicalVolume(worldS, //its solid
>> 127 fpDefaultMaterial, //its material
>> 128 "World"); //its name
107 129
108 // Visualization attributes 130 // Visualization attributes
109 G4VisAttributes* worldVisAtt = new G4VisAttr << 131 G4VisAttributes* worldVisAtt =
>> 132 new G4VisAttributes(G4Colour(0.5,0.5,0.5,0.1)); //Gray
>> 133 //worldVisAtt->SetForceSolid(true);
110 worldVisAtt->SetVisibility(true); 134 worldVisAtt->SetVisibility(true);
111 worldLV->SetVisAttributes(worldVisAtt); 135 worldLV->SetVisAttributes(worldVisAtt);
112 136
113 G4VPhysicalVolume* worldPV = new G4PVPlaceme << 137
114 << 138 G4VPhysicalVolume* worldPV = new G4PVPlacement(0, //no rotation
115 << 139 G4ThreeVector(), //at (0,0,0)
116 << 140 "World", //its name
117 << 141 worldLV, //its logical volume
118 << 142 0, //its mother volume
119 << 143 false, //no boolean operation
120 << 144 0, //copy number
>> 145 true); // checking overlaps forced to YES
121 146
122 // ========================================= 147 // ===============================================
123 // HOMOGENEOUS MEDIUM - LARGE SPHERE VOLUME 148 // HOMOGENEOUS MEDIUM - LARGE SPHERE VOLUME
124 // ========================================= 149 // ===============================================
125 150
126 // Radius of water sphere calculated as over 151 // Radius of water sphere calculated as overall dimensions of neuron.
127 fTotMassMedium = fNeuronLoadParamz->GetTotMa << 152 fTotMassMedium = fNeuronLoadParamz->GetTotMassMedium() ;
128 fTotSurfMedium = fNeuronLoadParamz->GetTotSu << 153 fTotSurfMedium = fNeuronLoadParamz->GetTotSurfMedium() ;
129 G4double RadiusMedium = fNeuronLoadParamz->G << 154 G4double RadiusMedium = fNeuronLoadParamz->GetdiagnlLength()*um / 2.;
130 G4cout << " Radius of homogeneous medium is << 155 G4cout << " Radius of homogeneous medium is calculated : "
131 << G4endl; << 156 << RadiusMedium/um <<" um"<< G4endl;
132 G4VSolid* mediumS = new G4Orb("Medium", Radi << 157 G4VSolid* mediumS = new G4Orb("Medium", RadiusMedium);
133 << 158 //G4VSolid* mediumS = new G4Box("Medium", RadiusMedium*um,
134 G4LogicalVolume* mediumLV = new G4LogicalVol << 159 // RadiusMedium*um,RadiusMedium*um);
135 << 160
>> 161 G4LogicalVolume* mediumLV =
>> 162 new G4LogicalVolume(mediumS, fpWaterMaterial, "Medium");
>> 163
136 // Visualization attributes 164 // Visualization attributes
137 G4VisAttributes* mediumVisAtt = new G4VisAtt << 165 G4VisAttributes* mediumVisAtt =
138 // mediumVisAtt->SetForceSolid(true); << 166 new G4VisAttributes(G4Colour(0.0,1.0,0.0,0.02)); //Green
139 // mediumVisAtt->SetForceWireframe (true); << 167 //mediumVisAtt->SetForceSolid(true);
>> 168 //mediumVisAtt->SetForceWireframe (true);
140 mediumVisAtt->SetForceLineSegmentsPerCircle( 169 mediumVisAtt->SetForceLineSegmentsPerCircle(180);
141 mediumVisAtt->SetVisibility(true); 170 mediumVisAtt->SetVisibility(true);
142 mediumLV->SetVisAttributes(mediumVisAtt); 171 mediumLV->SetVisAttributes(mediumVisAtt);
143 << 172
144 G4VPhysicalVolume* mediumPV = << 173 G4VPhysicalVolume* mediumPV = new G4PVPlacement(0,
145 new G4PVPlacement(0, G4ThreeVector(), "Med << 174 G4ThreeVector(),
>> 175 "Medium",
>> 176 mediumLV,
>> 177 worldPV,
>> 178 false,
>> 179 0,
>> 180 fCheckOverlaps);
146 181
147 // ========================================= 182 // ===============================================
148 // TARGET - BOUNDING SLICE including NEURON << 183 // TARGET - BOUNDING SLICE including NEURON
149 // ========================================= 184 // ===============================================
150 185
151 // Dimensions of bounding slice volume defin 186 // Dimensions of bounding slice volume defined as overall measure of neuron
152 << 187
153 G4double TargetSizeX = fNeuronLoadParamz->Ge << 188 G4double TargetSizeX = fNeuronLoadParamz->GetwidthB()*um;
154 G4double TargetSizeY = fNeuronLoadParamz->Ge << 189 G4double TargetSizeY = fNeuronLoadParamz->GetheightB()*um;
155 G4double TargetSizeZ = fNeuronLoadParamz->Ge << 190 G4double TargetSizeZ = fNeuronLoadParamz->GetdepthB()*um;
156 fTotMassSlice = fNeuronLoadParamz->GetTotMas << 191 fTotMassSlice = fNeuronLoadParamz->GetTotMassSlice() ;
157 G4cout << " Overall dimensions (um) of neuro << 192 G4cout << " Overall dimensions (um) of neuron morphology : " << "\n"
158 << "\n" << 193 << '\t'<< " width = " <<TargetSizeX/um<< " height = " <<
159 << '\t' << " width = " << TargetSizeX << 194 TargetSizeY/um
160 << " depth = " << TargetSizeZ / um << << 195 << " depth = " <<TargetSizeZ/um<<G4endl;
161 << 196
162 G4cout << " Volume (um3), surface (um2) and 197 G4cout << " Volume (um3), surface (um2) and mass (ug) of Bounding Slice are"
163 << " calculated : " << 198 << " calculated : " << "\n"
164 << "\n" << 199 << '\t'<<fNeuronLoadParamz->GetTotVolSlice()/std::pow(um,3)<<"; "<<'\t'
165 << '\t' << fNeuronLoadParamz->GetTotV << 200 <<fNeuronLoadParamz->GetTotSurfSlice()/(um*um)
166 << fNeuronLoadParamz->GetTotSurfSlice << 201 <<"; "<<'\t'<<fNeuronLoadParamz->GetTotMassSlice()*1e6/g<< "\n "<<G4endl;
167 << fNeuronLoadParamz->GetTotMassSlice << 202
168 << 203 G4Box* boundingS = new G4Box("BoundingSlice",
169 G4Box* boundingS = << 204 TargetSizeX/2.,TargetSizeY/2.,TargetSizeZ/2.);
170 new G4Box("BoundingSlice", TargetSizeX / 2 <<
171 205
172 G4LogicalVolume* boundingLV = new G4LogicalV << 206 G4LogicalVolume* boundingLV =
>> 207 new G4LogicalVolume(boundingS,fpWaterMaterial, "BoundingSlice");
173 208
174 // Visualization attributes with opacity! 209 // Visualization attributes with opacity!
175 G4VisAttributes* TargetVisAtt = new G4VisAtt << 210 G4VisAttributes* TargetVisAtt =
>> 211 new G4VisAttributes(G4Colour(1.0,1.0,1.0,0.1));
176 TargetVisAtt->SetForceSolid(true); 212 TargetVisAtt->SetForceSolid(true);
177 TargetVisAtt->SetVisibility(true); 213 TargetVisAtt->SetVisibility(true);
178 boundingLV->SetVisAttributes(TargetVisAtt); << 214 boundingLV->SetVisAttributes(TargetVisAtt);
179 new G4PVPlacement(0, G4ThreeVector(), "Bound << 215 //G4VPhysicalVolume* boundingPV =
180 fCheckOverlaps); << 216 new G4PVPlacement(0,
>> 217 G4ThreeVector(),
>> 218 "BoundingSlice",
>> 219 boundingLV,
>> 220 mediumPV,
>> 221 false,
>> 222 0,
>> 223 fCheckOverlaps);
181 224
182 // ========================================= 225 // ===============================================
183 // NEURON MORPHOLOGY 226 // NEURON MORPHOLOGY
184 // ========================================= 227 // ===============================================
185 228
186 G4cout << " Volume (um3), surface (um2) and << 229 G4cout << " Volume (um3), surface (um2) and mass(ug) of Neuron "
187 << "are calculated : " << 230 << "are calculated : "<< "\n"
188 << "\n" << 231 << '\t'<<fNeuronLoadParamz->GetTotVolNeuron()/std::pow(um,3)<<"; "<<'\t'
189 << '\t' << fNeuronLoadParamz->GetTotV << 232 <<fNeuronLoadParamz->GetTotSurfNeuron()/(um*um)
190 << fNeuronLoadParamz->GetTotSurfNeuro << 233 <<"; "<<'\t'<<fNeuronLoadParamz->GetTotMassNeuron()*1e6/g<<G4endl;
191 << fNeuronLoadParamz->GetTotMassNeuro << 234 fTotMassNeuron = fNeuronLoadParamz->GetTotMassNeuron() ;
192 fTotMassNeuron = fNeuronLoadParamz->GetTotMa << 235 G4cout << " Total number of compartments into Neuron : "
193 G4cout << " Total number of compartments int << 236 <<fNeuronLoadParamz->GetnbNeuroncomp()<<G4endl;
194 << G4endl; <<
195 G4cout << " Shift values (um) for Neuron tra 237 G4cout << " Shift values (um) for Neuron translation are calculated : "
196 << "\n" << 238 << "\n"
197 << '\t' << " shiftX = " << fNeuronLoa << 239 << '\t' << " shiftX = " <<fNeuronLoadParamz->GetshiftX()<< " shiftY = "
198 << " shiftY = " << fNeuronLoadParamz- << 240 << fNeuronLoadParamz->GetshiftY()
199 << " shiftZ = " << fNeuronLoadParamz- << 241 << " shiftZ = " <<fNeuronLoadParamz->GetshiftZ()<< "\n"<< G4endl;
200 << 242
201 // Soma in Violet with opacity // 0.85,0.4 << 243 // Soma in Violet with opacity // 0.85,0.44,0.84
202 fSomaColour = new G4VisAttributes; << 244 fSomaColour = new G4VisAttributes;
203 fSomaColour->SetColour(G4Colour(G4Colour(22 << 245 fSomaColour->SetColour(G4Colour(G4Colour(22/255. , 200/255. , 30/255.)));
204 fSomaColour->SetForceSolid(true); // true << 246 fSomaColour->SetForceSolid(true); // true
205 fSomaColour->SetVisibility(true); << 247 fSomaColour->SetVisibility(true);
206 << 248
207 // Dendrites in Dark-Blue << 249 // Dendrites in Dark-Blue
208 fDendColour = new G4VisAttributes; << 250 fDendColour = new G4VisAttributes;
209 fDendColour->SetColour(G4Colour(G4Colour(0.0 << 251 fDendColour->SetColour(G4Colour(G4Colour(0.0, 0.0, 0.5)));
210 fDendColour->SetForceSolid(true); << 252 fDendColour->SetForceSolid(true);
211 // fDendColour->SetVisibility(true); << 253 //fDendColour->SetVisibility(true);
212 << 254
213 // Axon in Maroon << 255 // Axon in Maroon
214 fAxonColour = new G4VisAttributes; << 256 fAxonColour = new G4VisAttributes;
215 fAxonColour->SetColour(G4Colour(G4Colour(0.5 << 257 fAxonColour->SetColour(G4Colour(G4Colour(0.5, 0.0, 0.0)));
216 fAxonColour->SetForceSolid(true); << 258 fAxonColour->SetForceSolid(true);
217 fAxonColour->SetVisibility(true); << 259 fAxonColour->SetVisibility(true);
218 << 260
219 // Spines in Dark-Green << 261 // Spines in Dark-Green
220 fSpineColour = new G4VisAttributes; << 262 fSpineColour = new G4VisAttributes;
221 fSpineColour->SetColour(G4Colour(G4Colour(0. << 263 fSpineColour->SetColour(G4Colour(G4Colour(0.0 , 100/255. , 0.0)));
222 fSpineColour->SetForceSolid(true); << 264 fSpineColour->SetForceSolid(true);
223 fSpineColour->SetVisibility(true); << 265 fSpineColour->SetVisibility(true);
224 << 266
225 // Whole neuron in semitransparent navy blue << 267 // Whole neuron in semitransparent navy blue
226 fNeuronColour = new G4VisAttributes; << 268 fNeuronColour = new G4VisAttributes;
227 fNeuronColour->SetColour(G4Colour(G4Colour(0 << 269 fNeuronColour->SetColour(G4Colour(G4Colour(0.0,0.4,0.8,0.9)));
228 fNeuronColour->SetForceSolid(true); << 270 fNeuronColour->SetForceSolid(true);
229 fNeuronColour->SetVisibility(true); << 271 fNeuronColour->SetVisibility(true);
230 << 272
231 // Placement volumes: G4examples/extended/pa << 273 // Placement volumes: G4examples/extended/parameterisations/gflash
232 274
233 // ========================================= 275 // =======================================================================
234 // Structure-1: Soma 276 // Structure-1: Soma
235 277
236 // Create Target G4Region and add logical vo 278 // Create Target G4Region and add logical volume
237 // Active Geant4-DNA processes in this regio << 279 // Active Geant4-DNA processes in this region
238 fpRegion = new G4Region("Soma"); << 280 fpRegion = new G4Region("Soma");
239 G4ProductionCuts* cuts = new G4ProductionCut 281 G4ProductionCuts* cuts = new G4ProductionCuts();
240 G4double defCut = 1 * nanometer; << 282 G4double defCut = 1*nanometer;
241 cuts->SetProductionCut(defCut, "gamma"); << 283 cuts->SetProductionCut(defCut,"gamma");
242 cuts->SetProductionCut(defCut, "e-"); << 284 cuts->SetProductionCut(defCut,"e-");
243 cuts->SetProductionCut(defCut, "e+"); << 285 cuts->SetProductionCut(defCut,"e+");
244 cuts->SetProductionCut(defCut, "proton"); << 286 cuts->SetProductionCut(defCut,"proton");
>> 287
>> 288 fnbSomacomp = fNeuronLoadParamz->GetnbSomacomp() ;
>> 289 fMassSomaTot = fNeuronLoadParamz->GetMassSomaTot() ;
>> 290 fMassSomacomp = new G4double[fnbSomacomp];
>> 291 fPosSomacomp = new G4ThreeVector[fnbSomacomp];
>> 292 if (fNeuronLoadParamz->GetnbSomacomp()==0)
>> 293 {
>> 294 G4cout <<" ---- Soma not found! ---- "<< G4endl;
>> 295 }
>> 296 else
>> 297 {
>> 298 G4cout <<" ---- Soma for construction: ---- "<< G4endl;
>> 299 G4cout << " Total number of compartments into Soma : "
>> 300 <<fNeuronLoadParamz->GetnbSomacomp()<<G4endl;
>> 301 for (G4int i=0; i<fNeuronLoadParamz->GetnbSomacomp() ; i++)
>> 302 {
>> 303 fsomaS [i] = new G4Orb("Soma", fNeuronLoadParamz->GetRadSomacomp(i)* um);
>> 304 //fsomaS [i] = new G4Ellipsoid("Soma", Ra *um, Rb *um, Rc *um, 0,0) ;
>> 305 // you can change parameters of Soma with a single compartment
>> 306 fsomaLV[i] = new G4LogicalVolume(fsomaS[i], fpWaterMaterial, "Soma");
>> 307 fsomaLV[i] ->SetVisAttributes(fSomaColour);
>> 308 fsomaPV[i] = new G4PVPlacement(
>> 309 0, // no rotation
>> 310 G4ThreeVector(
>> 311 (fNeuronLoadParamz->GetPosSomacomp(i).x()-fNeuronLoadParamz->GetshiftX())*um,
>> 312 (fNeuronLoadParamz->GetPosSomacomp(i).y()-fNeuronLoadParamz->GetshiftY())*um,
>> 313 (fNeuronLoadParamz->GetPosSomacomp(i).z()-fNeuronLoadParamz->GetshiftZ())*um),
>> 314 fsomaLV[i],
>> 315 "Soma",
>> 316 boundingLV,
>> 317 false,
>> 318 i);
>> 319 fMassSomacomp[i] = fNeuronLoadParamz->GetMassSomacomp(i) ;
>> 320 fPosSomacomp[i] = fNeuronLoadParamz->GetPosSomacomp(i) ;
245 fpRegion->SetProductionCuts(cuts); 321 fpRegion->SetProductionCuts(cuts);
246 G4ThreeVector shift(fNeuronLoadParamz->Getsh << 322 fpRegion->AddRootLogicalVolume(fsomaLV[i]);
247 fNeuronLoadParamz->Getsh <<
248 323
249 G4int n = fNeuronLoadParamz->GetnbSomacomp() << 324 }
250 if (n <= 0) { << 325 }
251 G4cout << " ---- Soma not found! ---- " << <<
252 } <<
253 else { <<
254 G4cout << " ---- Soma for construction: -- <<
255 G4cout << " Total number of compartments i <<
256 fnbSomacomp = n; <<
257 fMassSomaTot = fNeuronLoadParamz->GetMassS <<
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", fNeuronLoa <<
266 // you can change parameters of Soma wit <<
267 G4ThreeVector pos = (fNeuronLoadParamz-> <<
268 fsomaLV[i] = new G4LogicalVolume(fsomaS[ <<
269 fsomaLV[i]->SetVisAttributes(fSomaColour <<
270 fsomaPV[i] = new G4PVPlacement(0, // no <<
271 pos, fsom <<
272 fMassSomacomp[i] = fNeuronLoadParamz->Ge <<
273 fPosSomacomp[i] = fNeuronLoadParamz->Get <<
274 fpRegion->AddRootLogicalVolume(fsomaLV[i <<
275 } <<
276 } <<
277 326
278 // ========================================= << 327 // =======================================================================
279 // Structure-2: Dendrites 328 // Structure-2: Dendrites
280 329
281 n = fNeuronLoadParamz->GetnbDendritecomp(); << 330 // Active Geant4-DNA processes in this region
282 if (n <= 0) { << 331 fpRegion = new G4Region("Dendrites");
283 G4cout << " ---- Dendrites not found! ---- << 332
>> 333 fnbDendritecomp = fNeuronLoadParamz->GetnbDendritecomp() ;
>> 334 fMassDendTot = fNeuronLoadParamz->GetMassDendTot() ;
>> 335 fMassDendcomp = new G4double[fnbDendritecomp];
>> 336 fDistADendSoma = new G4double[fnbDendritecomp];
>> 337 fDistBDendSoma = new G4double[fnbDendritecomp];
>> 338 fPosDendcomp = new G4ThreeVector[fnbDendritecomp];
>> 339 if (fNeuronLoadParamz->GetnbDendritecomp()==0)
>> 340 {
>> 341 G4cout <<" ---- Dendrites not found! ---- "<< G4endl;
284 } 342 }
285 else { << 343 else
286 fnbDendritecomp = n; << 344 {
287 G4cout << " ---- Dendrites for constructio << 345 G4cout <<" ---- Dendrites for construction: ---- "<< G4endl;
288 G4cout << " Total number of compartments i << 346 G4cout << " Total number of compartments into Dendrites : "
289 << 347 <<fNeuronLoadParamz->GetnbDendritecomp()<<G4endl;
290 // Active Geant4-DNA processes in this reg << 348 for (G4int i=1; i<fNeuronLoadParamz->GetnbDendritecomp() ; i++)
291 auto region = new G4Region("Dendrites"); << 349 {
292 region->SetProductionCuts(cuts); << 350 fdendriteS [i] = new G4Tubs( "Dendrites",
293 << 351 0., fNeuronLoadParamz->GetRadDendcomp(i)*um,
294 fMassDendTot = fNeuronLoadParamz->GetMassD << 352 fNeuronLoadParamz->GetHeightDendcomp(i)*um/2.,
295 fMassDendcomp.resize(n, 0.0); << 353 0., 2.*pi );
296 fDistADendSoma.resize(n, 0.0); << 354 fdendriteLV[i] = new G4LogicalVolume(fdendriteS[i],fpWaterMaterial,
297 fDistBDendSoma.resize(n, 0.0); << 355 "Dendrites");
298 fPosDendcomp.resize(n); << 356 fdendriteLV[i] ->SetVisAttributes(fDendColour);
299 fdendriteS.resize(n, nullptr); << 357
300 fdendriteLV.resize(n, nullptr); << 358 fdendritePV[i] = new G4PVPlacement(
301 fdendritePV.resize(n, nullptr); << 359 // rot, // rotation checking with function ComputeTransformation!
302 for (G4int i = 1; i < n; ++i) { << 360 G4Transform3D( fNeuronLoadParamz->
303 fdendriteS[i] = new G4Tubs("Dendrites", << 361 GetRotDendcomp(i), // RotationMatrix with Inverse
304 fNeuronLoadPa << 362 G4ThreeVector(
305 fdendriteLV[i] = new G4LogicalVolume(fde << 363 (fNeuronLoadParamz->GetPosDendcomp(i).x()-fNeuronLoadParamz->GetshiftX())*um,
306 fdendriteLV[i]->SetVisAttributes(fDendCo << 364 (fNeuronLoadParamz->GetPosDendcomp(i).y()-fNeuronLoadParamz->GetshiftY())*um,
307 << 365 (fNeuronLoadParamz->GetPosDendcomp(i).z()-fNeuronLoadParamz->GetshiftZ())*um)),
308 G4ThreeVector pos = (fNeuronLoadParamz-> << 366 fdendriteLV[i],
309 // rotation checking with function Compu << 367 "Dendrites",
310 // RotationMatrix with Inverse << 368 boundingLV,
311 fdendritePV[i] = new G4PVPlacement(G4Tra << 369 false,
312 fdend << 370 i);
313 fMassDendcomp[i] = fNeuronLoadParamz->Ge << 371 fMassDendcomp[i] = fNeuronLoadParamz->GetMassDendcomp(i) ;
314 fPosDendcomp[i] = fNeuronLoadParamz->Get << 372 fPosDendcomp[i] = fNeuronLoadParamz->GetPosDendcomp(i) ;
315 fDistADendSoma[i] = fNeuronLoadParamz->G << 373 fDistADendSoma[i] = fNeuronLoadParamz->GetDistADendSoma(i) ;
316 fDistBDendSoma[i] = fNeuronLoadParamz->G << 374 fDistBDendSoma[i] = fNeuronLoadParamz->GetDistBDendSoma(i) ;
317 region->AddRootLogicalVolume(fdendriteLV << 375 fpRegion->SetProductionCuts(cuts);
>> 376 fpRegion->AddRootLogicalVolume(fdendriteLV[i]);
>> 377
318 } 378 }
319 } 379 }
320 << 380
321 // ========================================= 381 // =======================================================================
322 // Structure-3: Axon << 382 // Structure-3: Axon
>> 383
>> 384 // Active Geant4-DNA processes in this region
>> 385 fpRegion = new G4Region("Axon");
323 386
324 n = fNeuronLoadParamz->GetnbAxoncomp(); << 387 fnbAxoncomp = fNeuronLoadParamz->GetnbAxoncomp() ;
325 if (n <= 0) { << 388 fMassAxonTot = fNeuronLoadParamz->GetMassAxonTot() ;
326 G4cout << " ---- Axon not found! ---- " << << 389 fMassAxoncomp = new G4double[fnbAxoncomp];
>> 390 fDistAxonsoma = new G4double[fnbAxoncomp];
>> 391 fPosAxoncomp = new G4ThreeVector[fnbAxoncomp];
>> 392 if (fNeuronLoadParamz->GetnbAxoncomp()==0)
>> 393 {
>> 394 G4cout <<" ---- Axon not found! ---- "<< G4endl;
327 } 395 }
328 else { << 396 else
329 G4cout << " ---- Axon for construction: -- << 397 {
330 G4cout << " Total number of compartments i << 398 G4cout <<" ---- Axon for construction: ---- "<< G4endl;
331 // Active Geant4-DNA processes in this reg << 399 G4cout << " Total number of compartments into Axon : "
332 auto region = new G4Region("Axon"); << 400 << fNeuronLoadParamz->GetnbAxoncomp() <<G4endl;
333 region->SetProductionCuts(cuts); << 401 for (G4int i=1; i< fNeuronLoadParamz->GetnbAxoncomp() ; i++)
334 << 402 {
335 fnbAxoncomp = n; << 403 faxonS [i] = new G4Tubs( "Axon",
336 fMassAxonTot = fNeuronLoadParamz->GetMassA << 404 0., fNeuronLoadParamz->GetRadAxoncomp(i)*um,
337 fMassAxoncomp.resize(n, 0.0); << 405 fNeuronLoadParamz->GetHeightAxoncomp(i)*um/2.,
338 fDistAxonsoma.resize(n, 0.0); << 406 0., 2.*pi );
339 fPosAxoncomp.resize(n); << 407 faxonLV[i] = new G4LogicalVolume(faxonS[i], fpWaterMaterial, "Axon");
340 faxonS.resize(n, nullptr); << 408 faxonLV[i] ->SetVisAttributes(fAxonColour);
341 faxonLV.resize(n, nullptr); << 409 faxonPV[i] = new G4PVPlacement(G4Transform3D(
342 faxonPV.resize(n, nullptr); << 410 fNeuronLoadParamz->GetRotAxoncomp(i), // RotationMatrix with Inverse
343 << 411 G4ThreeVector(
344 for (G4int i = 1; i < n; ++i) { << 412 (fNeuronLoadParamz->GetPosAxoncomp(i).x()-fNeuronLoadParamz->GetshiftX())*um,
345 faxonS[i] = new G4Tubs("Axon", 0., fNeur << 413 (fNeuronLoadParamz->GetPosAxoncomp(i).y()-fNeuronLoadParamz->GetshiftY())*um,
346 fNeuronLoadParamz << 414 (fNeuronLoadParamz->GetPosAxoncomp(i).z()-fNeuronLoadParamz->GetshiftZ())*um)),
347 faxonLV[i] = new G4LogicalVolume(faxonS[ << 415 faxonLV[i],
348 faxonLV[i]->SetVisAttributes(fAxonColour << 416 "Axon",
349 // RotationMatrix with Inverse << 417 boundingLV,
350 G4ThreeVector pos = (fNeuronLoadParamz-> << 418 false,
351 faxonPV[i] = new G4PVPlacement(G4Transfo << 419 i);
352 faxonLV[i << 420 fMassAxoncomp[i] = fNeuronLoadParamz->GetMassAxoncomp(i) ;
353 fMassAxoncomp[i] = fNeuronLoadParamz->Ge << 421 fPosAxoncomp[i] = fNeuronLoadParamz->GetPosAxoncomp(i) ;
354 fPosAxoncomp[i] = fNeuronLoadParamz->Get << 422 fDistAxonsoma[i] = fNeuronLoadParamz->GetDistAxonsoma(i) ;
355 fDistAxonsoma[i] = fNeuronLoadParamz->Ge << 423 fpRegion->SetProductionCuts(cuts);
356 region->AddRootLogicalVolume(faxonLV[i]) << 424 fpRegion->AddRootLogicalVolume(faxonLV[i]);
>> 425
357 } 426 }
358 } << 427 }
359 // ========================================= 428 // =======================================================================
360 // Structure-4: Spines << 429 // Structure-4: Spines
361 if (fNeuronLoadParamz->GetnbSpinecomp() == 0 << 430 if (fNeuronLoadParamz->GetnbSpinecomp()==0)
362 G4cout << " ---- Spines not found! ---- " << 431 {
>> 432 G4cout <<" ---- Spines not found! ---- "<< G4endl;
363 } 433 }
364 else { << 434 else
365 G4cout << " ---- Spines for construction: << 435 {
366 G4cout << " Total number of compartments i << 436 G4cout <<" ---- Spines for construction: ---- "<< G4endl;
367 << G4endl; << 437 G4cout << " Total number of compartments into Spines : "
>> 438 << fNeuronLoadParamz->GetnbSpinecomp() <<G4endl;
368 } 439 }
369 << 440
370 G4cout << "\n ---- End of Neuron Constructio << 441 G4cout <<"\n ---- End of Neuron Construction! ---- "
371 << "\n ============================== << 442 << "\n ========================================================== \n"
372 << G4endl; << 443 << G4endl;
373 << 444
374 // ========================================= 445 // =======================================================================
375 // Active Geant4-DNA processes in BoundingSl 446 // Active Geant4-DNA processes in BoundingSlice with whole neuron structure
376 // fpRegion = new G4Region("BoundingSlice"); << 447 //fpRegion = new G4Region("BoundingSlice");
377 // fpRegion->SetProductionCuts(cuts); << 448 //fpRegion->SetProductionCuts(cuts);
378 // fpRegion->AddRootLogicalVolume(boundingLV << 449 // fpRegion->AddRootLogicalVolume(boundingLV);
379 450
380 return worldPV; 451 return worldPV;
381 } 452 }
382 453