ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DetectorConstruction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DetectorConstruction.cc
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 //
41 
42 #include "DetectorConstruction.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4Region.hh"
45 #include "G4ProductionCuts.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4Colour.hh"
48 #include "G4VisAttributes.hh"
49 #include "G4RotationMatrix.hh"
50 #include <algorithm>
51 #include <iostream>
52 
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54 
55 DetectorConstruction::DetectorConstruction()
56 :G4VUserDetectorConstruction(),
57  fpDefaultMaterial(0), fpWaterMaterial(0),
58  fpRegion(0), fCheckOverlaps(false)
59 {
60 
61 }
62 
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64 
65 DetectorConstruction::~DetectorConstruction()
66 {
67 
68 }
69 
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71 
72 G4VPhysicalVolume* DetectorConstruction::Construct()
73 {
74 
75  // load Neuron positions and obtain parameters!
76  fNeuronLoadParamz = new NeuronLoadDataFile() ;
77 
78  DefineMaterials();
79  return ConstructDetector();
80 
81 }
82 
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
84 
85 void DetectorConstruction::DefineMaterials()
86 {
87  // Water is defined from NIST material database
88  G4NistManager * man = G4NistManager::Instance();
89  G4Material * H2O = man->FindOrBuildMaterial("G4_WATER");
90  // Default materials in setup.
91  fpWaterMaterial = H2O;
92  G4Material * Vacuum = man->FindOrBuildMaterial("G4_Galactic");
93  //G4Material * Air = man->FindOrBuildMaterial("G4_AIR");
94  fpDefaultMaterial = Vacuum;
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
98 
99 G4VPhysicalVolume* DetectorConstruction::ConstructDetector()
100 {
101  G4cout <<" ---- Begin of Neuron Construction! ---- "
102  <<"\n"<<" =========================================================="<<G4endl;
103 
104  // ===============================================
105  // WORLD VOLUME - filled by default material
106  // ===============================================
107 
108  // Dimensions of world volume are calculated as overall dimensions of neuron!
109 
110  G4double worldSizeX;
111  worldSizeX = 1.0*fNeuronLoadParamz->GetdiagnlLength()*um;
112 
113  if (!worldSizeX)
114  {
115  worldSizeX = 10.*cm;
116  }
117 
118  G4double worldSizeY = worldSizeX;
119  G4double worldSizeZ = worldSizeX;
120  G4cout << " Side length of word volume is calculated : "
121  << worldSizeX/um <<" um"<< G4endl;
122  G4VSolid* worldS = new G4Box("World", //its name
123  worldSizeX/2, worldSizeY/2, worldSizeZ/2); //its size
124 
125  G4LogicalVolume* worldLV = new G4LogicalVolume(worldS, //its solid
126  fpDefaultMaterial, //its material
127  "World"); //its name
128 
129  // Visualization attributes
130  G4VisAttributes* worldVisAtt =
131  new G4VisAttributes(G4Colour(0.5,0.5,0.5,0.1)); //Gray
132  //worldVisAtt->SetForceSolid(true);
133  worldVisAtt->SetVisibility(true);
134  worldLV->SetVisAttributes(worldVisAtt);
135 
136 
137  G4VPhysicalVolume* worldPV = new G4PVPlacement(0, //no rotation
138  G4ThreeVector(), //at (0,0,0)
139  "World", //its name
140  worldLV, //its logical volume
141  0, //its mother volume
142  false, //no boolean operation
143  0, //copy number
144  true); // checking overlaps forced to YES
145 
146  // ===============================================
147  // HOMOGENEOUS MEDIUM - LARGE SPHERE VOLUME
148  // ===============================================
149 
150  // Radius of water sphere calculated as overall dimensions of neuron.
151  fTotMassMedium = fNeuronLoadParamz->GetTotMassMedium() ;
152  fTotSurfMedium = fNeuronLoadParamz->GetTotSurfMedium() ;
153  G4double RadiusMedium = fNeuronLoadParamz->GetdiagnlLength()*um / 2.;
154  G4cout << " Radius of homogeneous medium is calculated : "
155  << RadiusMedium/um <<" um"<< G4endl;
156  G4VSolid* mediumS = new G4Orb("Medium", RadiusMedium);
157  //G4VSolid* mediumS = new G4Box("Medium", RadiusMedium*um,
158  // RadiusMedium*um,RadiusMedium*um);
159 
160  G4LogicalVolume* mediumLV =
161  new G4LogicalVolume(mediumS, fpWaterMaterial, "Medium");
162 
163  // Visualization attributes
164  G4VisAttributes* mediumVisAtt =
165  new G4VisAttributes(G4Colour(0.0,1.0,0.0,0.02)); //Green
166  //mediumVisAtt->SetForceSolid(true);
167  //mediumVisAtt->SetForceWireframe (true);
168  mediumVisAtt->SetForceLineSegmentsPerCircle(180);
169  mediumVisAtt->SetVisibility(true);
170  mediumLV->SetVisAttributes(mediumVisAtt);
171 
172  G4VPhysicalVolume* mediumPV = new G4PVPlacement(0,
173  G4ThreeVector(),
174  "Medium",
175  mediumLV,
176  worldPV,
177  false,
178  0,
179  fCheckOverlaps);
180 
181  // ===============================================
182  // TARGET - BOUNDING SLICE including NEURON
183  // ===============================================
184 
185  // Dimensions of bounding slice volume defined as overall measure of neuron
186 
187  G4double TargetSizeX = fNeuronLoadParamz->GetwidthB()*um;
188  G4double TargetSizeY = fNeuronLoadParamz->GetheightB()*um;
189  G4double TargetSizeZ = fNeuronLoadParamz->GetdepthB()*um;
190  fTotMassSlice = fNeuronLoadParamz->GetTotMassSlice() ;
191  G4cout << " Overall dimensions (um) of neuron morphology : " << "\n"
192  << '\t'<< " width = " <<TargetSizeX/um<< " height = " <<
193  TargetSizeY/um
194  << " depth = " <<TargetSizeZ/um<<G4endl;
195 
196  G4cout << " Volume (um3), surface (um2) and mass (ug) of Bounding Slice are"
197  << " calculated : " << "\n"
198  << '\t'<<fNeuronLoadParamz->GetTotVolSlice()/std::pow(um,3)<<"; "<<'\t'
199  <<fNeuronLoadParamz->GetTotSurfSlice()/(um*um)
200  <<"; "<<'\t'<<fNeuronLoadParamz->GetTotMassSlice()*1e6/g<< "\n "<<G4endl;
201 
202  G4Box* boundingS = new G4Box("BoundingSlice",
203  TargetSizeX/2.,TargetSizeY/2.,TargetSizeZ/2.);
204 
205  G4LogicalVolume* boundingLV =
206  new G4LogicalVolume(boundingS,fpWaterMaterial, "BoundingSlice");
207 
208  // Visualization attributes with opacity!
209  G4VisAttributes* TargetVisAtt =
210  new G4VisAttributes(G4Colour(1.0,1.0,1.0,0.1));
211  TargetVisAtt->SetForceSolid(true);
212  TargetVisAtt->SetVisibility(true);
213  boundingLV->SetVisAttributes(TargetVisAtt);
214  //G4VPhysicalVolume* boundingPV =
215  new G4PVPlacement(0,
216  G4ThreeVector(),
217  "BoundingSlice",
218  boundingLV,
219  mediumPV,
220  false,
221  0,
222  fCheckOverlaps);
223 
224  // ===============================================
225  // NEURON MORPHOLOGY
226  // ===============================================
227 
228  G4cout << " Volume (um3), surface (um2) and mass(ug) of Neuron "
229  << "are calculated : "<< "\n"
230  << '\t'<<fNeuronLoadParamz->GetTotVolNeuron()/std::pow(um,3)<<"; "<<'\t'
231  <<fNeuronLoadParamz->GetTotSurfNeuron()/(um*um)
232  <<"; "<<'\t'<<fNeuronLoadParamz->GetTotMassNeuron()*1e6/g<<G4endl;
233  fTotMassNeuron = fNeuronLoadParamz->GetTotMassNeuron() ;
234  G4cout << " Total number of compartments into Neuron : "
235  <<fNeuronLoadParamz->GetnbNeuroncomp()<<G4endl;
236  G4cout << " Shift values (um) for Neuron translation are calculated : "
237  << "\n"
238  << '\t' << " shiftX = " <<fNeuronLoadParamz->GetshiftX()<< " shiftY = "
239  << fNeuronLoadParamz->GetshiftY()
240  << " shiftZ = " <<fNeuronLoadParamz->GetshiftZ()<< "\n"<< G4endl;
241 
242  // Soma in Violet with opacity // 0.85,0.44,0.84
243  fSomaColour = new G4VisAttributes;
244  fSomaColour->SetColour(G4Colour(G4Colour(22/255. , 200/255. , 30/255.)));
245  fSomaColour->SetForceSolid(true); // true
246  fSomaColour->SetVisibility(true);
247 
248  // Dendrites in Dark-Blue
249  fDendColour = new G4VisAttributes;
250  fDendColour->SetColour(G4Colour(G4Colour(0.0, 0.0, 0.5)));
251  fDendColour->SetForceSolid(true);
252  //fDendColour->SetVisibility(true);
253 
254  // Axon in Maroon
255  fAxonColour = new G4VisAttributes;
256  fAxonColour->SetColour(G4Colour(G4Colour(0.5, 0.0, 0.0)));
257  fAxonColour->SetForceSolid(true);
258  fAxonColour->SetVisibility(true);
259 
260  // Spines in Dark-Green
261  fSpineColour = new G4VisAttributes;
262  fSpineColour->SetColour(G4Colour(G4Colour(0.0 , 100/255. , 0.0)));
263  fSpineColour->SetForceSolid(true);
264  fSpineColour->SetVisibility(true);
265 
266  // Whole neuron in semitransparent navy blue
267  fNeuronColour = new G4VisAttributes;
268  fNeuronColour->SetColour(G4Colour(G4Colour(0.0,0.4,0.8,0.9)));
269  fNeuronColour->SetForceSolid(true);
270  fNeuronColour->SetVisibility(true);
271 
272  // Placement volumes: G4examples/extended/parameterisations/gflash
273 
274  // =======================================================================
275  // Structure-1: Soma
276 
277  // Create Target G4Region and add logical volume
278  // Active Geant4-DNA processes in this region
279  fpRegion = new G4Region("Soma");
280  G4ProductionCuts* cuts = new G4ProductionCuts();
281  G4double defCut = 1*nanometer;
282  cuts->SetProductionCut(defCut,"gamma");
283  cuts->SetProductionCut(defCut,"e-");
284  cuts->SetProductionCut(defCut,"e+");
285  cuts->SetProductionCut(defCut,"proton");
286 
287  fnbSomacomp = fNeuronLoadParamz->GetnbSomacomp() ;
288  fMassSomaTot = fNeuronLoadParamz->GetMassSomaTot() ;
289  fMassSomacomp = new G4double[fnbSomacomp];
290  fPosSomacomp = new G4ThreeVector[fnbSomacomp];
291  if (fNeuronLoadParamz->GetnbSomacomp()==0)
292  {
293  G4cout <<" ---- Soma not found! ---- "<< G4endl;
294  }
295  else
296  {
297  G4cout <<" ---- Soma for construction: ---- "<< G4endl;
298  G4cout << " Total number of compartments into Soma : "
299  <<fNeuronLoadParamz->GetnbSomacomp()<<G4endl;
300  for (G4int i=0; i<fNeuronLoadParamz->GetnbSomacomp() ; i++)
301  {
302  fsomaS [i] = new G4Orb("Soma", fNeuronLoadParamz->GetRadSomacomp(i)* um);
303  //fsomaS [i] = new G4Ellipsoid("Soma", Ra *um, Rb *um, Rc *um, 0,0) ;
304  // you can change parameters of Soma with a single compartment
305  fsomaLV[i] = new G4LogicalVolume(fsomaS[i], fpWaterMaterial, "Soma");
306  fsomaLV[i] ->SetVisAttributes(fSomaColour);
307  fsomaPV[i] = new G4PVPlacement(
308  0, // no rotation
309  G4ThreeVector(
310  (fNeuronLoadParamz->GetPosSomacomp(i).x()-fNeuronLoadParamz->GetshiftX())*um,
311  (fNeuronLoadParamz->GetPosSomacomp(i).y()-fNeuronLoadParamz->GetshiftY())*um,
312  (fNeuronLoadParamz->GetPosSomacomp(i).z()-fNeuronLoadParamz->GetshiftZ())*um),
313  fsomaLV[i],
314  "Soma",
315  boundingLV,
316  false,
317  i);
318  fMassSomacomp[i] = fNeuronLoadParamz->GetMassSomacomp(i) ;
319  fPosSomacomp[i] = fNeuronLoadParamz->GetPosSomacomp(i) ;
320  fpRegion->SetProductionCuts(cuts);
321  fpRegion->AddRootLogicalVolume(fsomaLV[i]);
322 
323  }
324  }
325 
326  // =======================================================================
327  // Structure-2: Dendrites
328 
329  // Active Geant4-DNA processes in this region
330  fpRegion = new G4Region("Dendrites");
331 
332  fnbDendritecomp = fNeuronLoadParamz->GetnbDendritecomp() ;
333  fMassDendTot = fNeuronLoadParamz->GetMassDendTot() ;
334  fMassDendcomp = new G4double[fnbDendritecomp];
335  fDistADendSoma = new G4double[fnbDendritecomp];
336  fDistBDendSoma = new G4double[fnbDendritecomp];
337  fPosDendcomp = new G4ThreeVector[fnbDendritecomp];
338  if (fNeuronLoadParamz->GetnbDendritecomp()==0)
339  {
340  G4cout <<" ---- Dendrites not found! ---- "<< G4endl;
341  }
342  else
343  {
344  G4cout <<" ---- Dendrites for construction: ---- "<< G4endl;
345  G4cout << " Total number of compartments into Dendrites : "
346  <<fNeuronLoadParamz->GetnbDendritecomp()<<G4endl;
347  for (G4int i=1; i<fNeuronLoadParamz->GetnbDendritecomp() ; i++)
348  {
349  fdendriteS [i] = new G4Tubs( "Dendrites",
350  0., fNeuronLoadParamz->GetRadDendcomp(i)*um,
351  fNeuronLoadParamz->GetHeightDendcomp(i)*um/2.,
352  0., 2.*pi );
353  fdendriteLV[i] = new G4LogicalVolume(fdendriteS[i],fpWaterMaterial,
354  "Dendrites");
355  fdendriteLV[i] ->SetVisAttributes(fDendColour);
356 
357  fdendritePV[i] = new G4PVPlacement(
358  // rot, // rotation checking with function ComputeTransformation!
359  G4Transform3D( fNeuronLoadParamz->
360  GetRotDendcomp(i), // RotationMatrix with Inverse
361  G4ThreeVector(
362  (fNeuronLoadParamz->GetPosDendcomp(i).x()-fNeuronLoadParamz->GetshiftX())*um,
363  (fNeuronLoadParamz->GetPosDendcomp(i).y()-fNeuronLoadParamz->GetshiftY())*um,
364  (fNeuronLoadParamz->GetPosDendcomp(i).z()-fNeuronLoadParamz->GetshiftZ())*um)),
365  fdendriteLV[i],
366  "Dendrites",
367  boundingLV,
368  false,
369  i);
370  fMassDendcomp[i] = fNeuronLoadParamz->GetMassDendcomp(i) ;
371  fPosDendcomp[i] = fNeuronLoadParamz->GetPosDendcomp(i) ;
372  fDistADendSoma[i] = fNeuronLoadParamz->GetDistADendSoma(i) ;
373  fDistBDendSoma[i] = fNeuronLoadParamz->GetDistBDendSoma(i) ;
374  fpRegion->SetProductionCuts(cuts);
375  fpRegion->AddRootLogicalVolume(fdendriteLV[i]);
376 
377  }
378  }
379 
380  // =======================================================================
381  // Structure-3: Axon
382 
383  // Active Geant4-DNA processes in this region
384  fpRegion = new G4Region("Axon");
385 
386  fnbAxoncomp = fNeuronLoadParamz->GetnbAxoncomp() ;
387  fMassAxonTot = fNeuronLoadParamz->GetMassAxonTot() ;
388  fMassAxoncomp = new G4double[fnbAxoncomp];
389  fDistAxonsoma = new G4double[fnbAxoncomp];
390  fPosAxoncomp = new G4ThreeVector[fnbAxoncomp];
391  if (fNeuronLoadParamz->GetnbAxoncomp()==0)
392  {
393  G4cout <<" ---- Axon not found! ---- "<< G4endl;
394  }
395  else
396  {
397  G4cout <<" ---- Axon for construction: ---- "<< G4endl;
398  G4cout << " Total number of compartments into Axon : "
399  << fNeuronLoadParamz->GetnbAxoncomp() <<G4endl;
400  for (G4int i=1; i< fNeuronLoadParamz->GetnbAxoncomp() ; i++)
401  {
402  faxonS [i] = new G4Tubs( "Axon",
403  0., fNeuronLoadParamz->GetRadAxoncomp(i)*um,
404  fNeuronLoadParamz->GetHeightAxoncomp(i)*um/2.,
405  0., 2.*pi );
406  faxonLV[i] = new G4LogicalVolume(faxonS[i], fpWaterMaterial, "Axon");
407  faxonLV[i] ->SetVisAttributes(fAxonColour);
408  faxonPV[i] = new G4PVPlacement(G4Transform3D(
409  fNeuronLoadParamz->GetRotAxoncomp(i), // RotationMatrix with Inverse
410  G4ThreeVector(
411  (fNeuronLoadParamz->GetPosAxoncomp(i).x()-fNeuronLoadParamz->GetshiftX())*um,
412  (fNeuronLoadParamz->GetPosAxoncomp(i).y()-fNeuronLoadParamz->GetshiftY())*um,
413  (fNeuronLoadParamz->GetPosAxoncomp(i).z()-fNeuronLoadParamz->GetshiftZ())*um)),
414  faxonLV[i],
415  "Axon",
416  boundingLV,
417  false,
418  i);
419  fMassAxoncomp[i] = fNeuronLoadParamz->GetMassAxoncomp(i) ;
420  fPosAxoncomp[i] = fNeuronLoadParamz->GetPosAxoncomp(i) ;
421  fDistAxonsoma[i] = fNeuronLoadParamz->GetDistAxonsoma(i) ;
422  fpRegion->SetProductionCuts(cuts);
423  fpRegion->AddRootLogicalVolume(faxonLV[i]);
424 
425  }
426  }
427  // =======================================================================
428  // Structure-4: Spines
429  if (fNeuronLoadParamz->GetnbSpinecomp()==0)
430  {
431  G4cout <<" ---- Spines not found! ---- "<< G4endl;
432  }
433  else
434  {
435  G4cout <<" ---- Spines for construction: ---- "<< G4endl;
436  G4cout << " Total number of compartments into Spines : "
437  << fNeuronLoadParamz->GetnbSpinecomp() <<G4endl;
438  }
439 
440  G4cout <<"\n ---- End of Neuron Construction! ---- "
441  << "\n ========================================================== \n"
442  << G4endl;
443 
444  // =======================================================================
445  // Active Geant4-DNA processes in BoundingSlice with whole neuron structure
446  //fpRegion = new G4Region("BoundingSlice");
447  //fpRegion->SetProductionCuts(cuts);
448  // fpRegion->AddRootLogicalVolume(boundingLV);
449 
450  return worldPV;
451 }