neuron/src/DetectorConstruction.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // This example is provided by the Geant4-DNA collaboration
0027 // Any report or published results obtained using the Geant4-DNA software
0028 // shall cite the following Geant4-DNA collaboration publication:
0029 // Med. Phys. 37 (2010) 4692-4708
0030 // and papers
0031 // M. Batmunkh et al. J Radiat Res Appl Sci 8 (2015) 498-507
0032 // O. Belov et al. Physica Medica 32 (2016) 1510-1520
0033 // The Geant4-DNA web site is available at http://geant4-dna.org
0034 //
0035 // -------------------------------------------------------------------
0036 // November 2016
0037 // -------------------------------------------------------------------
0038 //
0039 /// \file DetectorConstruction.cc
0040 /// \brief Implementation of the DetectorConstruction class
0041
0042 #include "DetectorConstruction.hh"
0043
0044 #include "G4Colour.hh"
0045 #include "G4PhysicalConstants.hh"
0046 #include "G4ProductionCuts.hh"
0047 #include "G4Region.hh"
0048 #include "G4RotationMatrix.hh"
0049 #include "G4SystemOfUnits.hh"
0050 #include "G4VisAttributes.hh"
0051
0052 #include <algorithm>
0053 #include <iostream>
0054
0055 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0056
0057 DetectorConstruction::DetectorConstruction()
0058 {
0059   fNeuronLoadParamz = new NeuronLoadDataFile();
0060   DefineMaterials();
0061 }
0062
0063 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0064
0065 DetectorConstruction::~DetectorConstruction()
0066 {
0067   delete fNeuronLoadParamz;
0068 }
0069
0070 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0071
0072 void DetectorConstruction::DefineMaterials()
0073 {
0074   // Water is defined from NIST material database
0075   G4NistManager* man = G4NistManager::Instance();
0076   fpWaterMaterial = man->FindOrBuildMaterial("G4_WATER");
0077   fpWorldMaterial = man->FindOrBuildMaterial("G4_Galactic");
0078 }
0079
0080 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0081
0082 G4VPhysicalVolume* DetectorConstruction::Construct()
0083 {
0084   G4cout << " ---- Begin of Neuron Construction! ---- "
0085          << "\n"
0086          << " ==========================================================" << G4endl;
0087
0088   // ===============================================
0089   // WORLD VOLUME - filled by default material
0090   // ===============================================
0091
0092   // Dimensions of world volume are calculated as overall dimensions of neuron!
0093
0094   G4double worldSizeX;
0095   worldSizeX = fNeuronLoadParamz->GetdiagnlLength() * um;
0096
0097   if (worldSizeX <= 0.0) {
0098     worldSizeX = 10. * cm;
0099   }
0100
0101   G4double worldSizeY = worldSizeX;
0102   G4double worldSizeZ = worldSizeX;
0103   G4cout << " Side length of word volume is calculated : " << worldSizeX / um << " um" << G4endl;
0104   G4VSolid* worldS = new G4Box("World", worldSizeX / 2, worldSizeY / 2, worldSizeZ / 2);
0105
0106   G4LogicalVolume* worldLV = new G4LogicalVolume(worldS, fpWorldMaterial, "World");
0107
0108   // Visualization attributes
0109   G4VisAttributes* worldVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5, 0.1));  // Gray
0110   worldVisAtt->SetVisibility(true);
0111   worldLV->SetVisAttributes(worldVisAtt);
0112
0113   G4VPhysicalVolume* worldPV = new G4PVPlacement(0,  // no rotation
0114                                                  G4ThreeVector(),  // at (0,0,0)
0115                                                  "World",  // its name
0116                                                  worldLV,  // its logical volume
0117                                                  0,  // its mother  volume
0118                                                  false,  // no boolean operation
0119                                                  0,  // copy number
0120                                                  true);  // checking overlaps forced to YES
0121
0122   // ===============================================
0123   // HOMOGENEOUS MEDIUM - LARGE SPHERE VOLUME
0124   // ===============================================
0125
0126   // Radius of water sphere calculated as overall dimensions of neuron.
0127   fTotMassMedium = fNeuronLoadParamz->GetTotMassMedium();
0128   fTotSurfMedium = fNeuronLoadParamz->GetTotSurfMedium();
0129   G4double RadiusMedium = fNeuronLoadParamz->GetdiagnlLength() * um / 2.;
0130   G4cout << " Radius of homogeneous medium is calculated : " << RadiusMedium / um << " um"
0131          << G4endl;
0132   G4VSolid* mediumS = new G4Orb("Medium", RadiusMedium);
0133
0134   G4LogicalVolume* mediumLV = new G4LogicalVolume(mediumS, fpWaterMaterial, "Medium");
0135
0136   // Visualization attributes
0137   G4VisAttributes* mediumVisAtt = new G4VisAttributes(G4Colour(0.0, 1.0, 0.0, 0.02));  // Green
0138   // mediumVisAtt->SetForceSolid(true);
0139   // mediumVisAtt->SetForceWireframe (true);
0140   mediumVisAtt->SetForceLineSegmentsPerCircle(180);
0141   mediumVisAtt->SetVisibility(true);
0142   mediumLV->SetVisAttributes(mediumVisAtt);
0143
0144   G4VPhysicalVolume* mediumPV =
0145     new G4PVPlacement(0, G4ThreeVector(), "Medium", mediumLV, worldPV, false, 0, fCheckOverlaps);
0146
0147   // ===============================================
0148   // TARGET - BOUNDING SLICE including NEURON
0149   // ===============================================
0150
0151   // Dimensions of bounding slice volume defined as overall measure of neuron
0152
0153   G4double TargetSizeX = fNeuronLoadParamz->GetwidthB() * um;
0154   G4double TargetSizeY = fNeuronLoadParamz->GetheightB() * um;
0155   G4double TargetSizeZ = fNeuronLoadParamz->GetdepthB() * um;
0156   fTotMassSlice = fNeuronLoadParamz->GetTotMassSlice();
0157   G4cout << " Overall dimensions (um) of neuron morphology : "
0158          << "\n"
0159          << '\t' << " width = " << TargetSizeX / um << " height = " << TargetSizeY / um
0160          << " depth = " << TargetSizeZ / um << G4endl;
0161
0162   G4cout << " Volume (um3), surface (um2) and mass (ug) of Bounding Slice are"
0163          << " calculated : "
0164          << "\n"
0165          << '\t' << fNeuronLoadParamz->GetTotVolSlice() / std::pow(um, 3) << "; " << '\t'
0166          << fNeuronLoadParamz->GetTotSurfSlice() / (um * um) << "; " << '\t'
0167          << fNeuronLoadParamz->GetTotMassSlice() * 1e6 / g << "\n " << G4endl;
0168
0169   G4Box* boundingS =
0170     new G4Box("BoundingSlice", TargetSizeX / 2., TargetSizeY / 2., TargetSizeZ / 2.);
0171
0172   G4LogicalVolume* boundingLV = new G4LogicalVolume(boundingS, fpWaterMaterial, "BoundingSlice");
0173
0174   // Visualization attributes with opacity!
0175   G4VisAttributes* TargetVisAtt = new G4VisAttributes(G4Colour(1.0, 1.0, 1.0, 0.1));
0176   TargetVisAtt->SetForceSolid(true);
0177   TargetVisAtt->SetVisibility(true);
0178   boundingLV->SetVisAttributes(TargetVisAtt);
0179   new G4PVPlacement(0, G4ThreeVector(), "BoundingSlice", boundingLV, mediumPV, false, 0,
0180                     fCheckOverlaps);
0181
0182   // ===============================================
0183   // NEURON MORPHOLOGY
0184   // ===============================================
0185
0186   G4cout << " Volume (um3), surface (um2) and mass(ug) of Neuron "
0187          << "are calculated : "
0188          << "\n"
0189          << '\t' << fNeuronLoadParamz->GetTotVolNeuron() / std::pow(um, 3) << "; " << '\t'
0190          << fNeuronLoadParamz->GetTotSurfNeuron() / (um * um) << "; " << '\t'
0191          << fNeuronLoadParamz->GetTotMassNeuron() * 1e6 / g << G4endl;
0192   fTotMassNeuron = fNeuronLoadParamz->GetTotMassNeuron();
0193   G4cout << " Total number of compartments into Neuron : " << fNeuronLoadParamz->GetnbNeuroncomp()
0194          << G4endl;
0195   G4cout << " Shift values (um) for Neuron translation are calculated : "
0196          << "\n"
0197          << '\t' << " shiftX = " << fNeuronLoadParamz->GetshiftX()
0198          << " shiftY = " << fNeuronLoadParamz->GetshiftY()
0199          << " shiftZ = " << fNeuronLoadParamz->GetshiftZ() << G4endl;
0200
0201   // Soma in Violet with opacity   // 0.85,0.44,0.84
0202   fSomaColour = new G4VisAttributes;
0203   fSomaColour->SetColour(G4Colour(G4Colour(22 / 255., 200 / 255., 30 / 255.)));
0204   fSomaColour->SetForceSolid(true);  // true
0205   fSomaColour->SetVisibility(true);
0206
0207   // Dendrites in Dark-Blue
0208   fDendColour = new G4VisAttributes;
0209   fDendColour->SetColour(G4Colour(G4Colour(0.0, 0.0, 0.5)));
0210   fDendColour->SetForceSolid(true);
0211   // fDendColour->SetVisibility(true);
0212
0213   // Axon in Maroon
0214   fAxonColour = new G4VisAttributes;
0215   fAxonColour->SetColour(G4Colour(G4Colour(0.5, 0.0, 0.0)));
0216   fAxonColour->SetForceSolid(true);
0217   fAxonColour->SetVisibility(true);
0218
0219   // Spines in Dark-Green
0220   fSpineColour = new G4VisAttributes;
0221   fSpineColour->SetColour(G4Colour(G4Colour(0.0, 100 / 255., 0.0)));
0222   fSpineColour->SetForceSolid(true);
0223   fSpineColour->SetVisibility(true);
0224
0225   // Whole neuron in semitransparent navy blue
0226   fNeuronColour = new G4VisAttributes;
0227   fNeuronColour->SetColour(G4Colour(G4Colour(0.0, 0.4, 0.8, 0.9)));
0228   fNeuronColour->SetForceSolid(true);
0229   fNeuronColour->SetVisibility(true);
0230
0231   // Placement volumes: G4examples/extended/parameterisations/gflash
0232
0233   // =======================================================================
0234   // Structure-1: Soma
0235
0236   // Create Target G4Region and add logical volume
0237   // Active Geant4-DNA processes in this region
0238   fpRegion = new G4Region("Soma");
0239   G4ProductionCuts* cuts = new G4ProductionCuts();
0240   G4double defCut = 1 * nanometer;
0241   cuts->SetProductionCut(defCut, "gamma");
0242   cuts->SetProductionCut(defCut, "e-");
0243   cuts->SetProductionCut(defCut, "e+");
0244   cuts->SetProductionCut(defCut, "proton");
0245   fpRegion->SetProductionCuts(cuts);
0246   G4ThreeVector shift(fNeuronLoadParamz->GetshiftX(), fNeuronLoadParamz->GetshiftY(),
0247                       fNeuronLoadParamz->GetshiftZ());
0248
0249   G4int n = fNeuronLoadParamz->GetnbSomacomp();
0250   if (n <= 0) {
0251     G4cout << " ---- Soma not found! ---- " << G4endl;
0252   }
0253   else {
0254     G4cout << " ---- Soma for construction: ---- " << G4endl;
0255     G4cout << " Total number of compartments into Soma : " << n << G4endl;
0256     fnbSomacomp = n;
0257     fMassSomaTot = fNeuronLoadParamz->GetMassSomaTot();
0258     fMassSomacomp.resize(n, 0.0);
0259     fPosSomacomp.resize(n);
0260     fsomaS.resize(n, nullptr);
0261     fsomaLV.resize(n, nullptr);
0262     fsomaPV.resize(n, nullptr);
0263
0264     for (G4int i = 0; i < n; ++i) {
0265       fsomaS[i] = new G4Orb("Soma", fNeuronLoadParamz->GetRadSomacomp(i) * um);
0266       // you can change parameters of Soma with a single compartment
0267       G4ThreeVector pos = (fNeuronLoadParamz->GetPosSomacomp(i) - shift) * um;
0268       fsomaLV[i] = new G4LogicalVolume(fsomaS[i], fpWaterMaterial, "Soma");
0269       fsomaLV[i]->SetVisAttributes(fSomaColour);
0270       fsomaPV[i] = new G4PVPlacement(0,  // no rotation
0271                                      pos, fsomaLV[i], "Soma", boundingLV, false, i);
0272       fMassSomacomp[i] = fNeuronLoadParamz->GetMassSomacomp(i);
0273       fPosSomacomp[i] = fNeuronLoadParamz->GetPosSomacomp(i);
0274       fpRegion->AddRootLogicalVolume(fsomaLV[i]);
0275     }
0276   }
0277
0278   // =======================================================================
0279   // Structure-2: Dendrites
0280
0281   n = fNeuronLoadParamz->GetnbDendritecomp();
0282   if (n <= 0) {
0283     G4cout << " ---- Dendrites not found! ---- " << G4endl;
0284   }
0285   else {
0286     fnbDendritecomp = n;
0287     G4cout << " ---- Dendrites for construction: ---- " << G4endl;
0288     G4cout << " Total number of compartments into Dendrites: " << n << G4endl;
0289
0290     // Active Geant4-DNA processes in this region
0291     auto region = new G4Region("Dendrites");
0292     region->SetProductionCuts(cuts);
0293
0294     fMassDendTot = fNeuronLoadParamz->GetMassDendTot();
0295     fMassDendcomp.resize(n, 0.0);
0296     fDistADendSoma.resize(n, 0.0);
0297     fDistBDendSoma.resize(n, 0.0);
0298     fPosDendcomp.resize(n);
0299     fdendriteS.resize(n, nullptr);
0300     fdendriteLV.resize(n, nullptr);
0301     fdendritePV.resize(n, nullptr);
0302     for (G4int i = 1; i < n; ++i) {
0303       fdendriteS[i] = new G4Tubs("Dendrites", 0., fNeuronLoadParamz->GetRadDendcomp(i) * um,
0304                                  fNeuronLoadParamz->GetHeightDendcomp(i) * um / 2., 0., 2. * pi);
0305       fdendriteLV[i] = new G4LogicalVolume(fdendriteS[i], fpWaterMaterial, "Dendrites");
0306       fdendriteLV[i]->SetVisAttributes(fDendColour);
0307
0308       G4ThreeVector pos = (fNeuronLoadParamz->GetPosDendcomp(i) - shift) * um;
0309       // rotation checking with function ComputeTransformation!
0310       // RotationMatrix with Inverse
0311       fdendritePV[i] = new G4PVPlacement(G4Transform3D(fNeuronLoadParamz->GetRotDendcomp(i), pos),
0312                                          fdendriteLV[i], "Dendrites", boundingLV, false, i);
0313       fMassDendcomp[i] = fNeuronLoadParamz->GetMassDendcomp(i);
0314       fPosDendcomp[i] = fNeuronLoadParamz->GetPosDendcomp(i);
0315       fDistADendSoma[i] = fNeuronLoadParamz->GetDistADendSoma(i);
0316       fDistBDendSoma[i] = fNeuronLoadParamz->GetDistBDendSoma(i);
0317       region->AddRootLogicalVolume(fdendriteLV[i]);
0318     }
0319   }
0320
0321   // =======================================================================
0322   // Structure-3: Axon
0323
0324   n = fNeuronLoadParamz->GetnbAxoncomp();
0325   if (n <= 0) {
0326     G4cout << " ---- Axon not found! ---- " << G4endl;
0327   }
0328   else {
0329     G4cout << " ---- Axon for construction: ---- " << G4endl;
0330     G4cout << " Total number of compartments into Axon : " << n << G4endl;
0331     // Active Geant4-DNA processes in this region
0332     auto region = new G4Region("Axon");
0333     region->SetProductionCuts(cuts);
0334
0335     fnbAxoncomp = n;
0336     fMassAxonTot = fNeuronLoadParamz->GetMassAxonTot();
0337     fMassAxoncomp.resize(n, 0.0);
0338     fDistAxonsoma.resize(n, 0.0);
0339     fPosAxoncomp.resize(n);
0340     faxonS.resize(n, nullptr);
0341     faxonLV.resize(n, nullptr);
0342     faxonPV.resize(n, nullptr);
0343
0344     for (G4int i = 1; i < n; ++i) {
0345       faxonS[i] = new G4Tubs("Axon", 0., fNeuronLoadParamz->GetRadAxoncomp(i) * um,
0346                              fNeuronLoadParamz->GetHeightAxoncomp(i) * um / 2., 0., 2. * pi);
0347       faxonLV[i] = new G4LogicalVolume(faxonS[i], fpWaterMaterial, "Axon");
0348       faxonLV[i]->SetVisAttributes(fAxonColour);
0349       // RotationMatrix with Inverse
0350       G4ThreeVector pos = (fNeuronLoadParamz->GetPosAxoncomp(i) - shift) * um;
0351       faxonPV[i] = new G4PVPlacement(G4Transform3D(fNeuronLoadParamz->GetRotAxoncomp(i), pos),
0352                                      faxonLV[i], "Axon", boundingLV, false, i);
0353       fMassAxoncomp[i] = fNeuronLoadParamz->GetMassAxoncomp(i);
0354       fPosAxoncomp[i] = fNeuronLoadParamz->GetPosAxoncomp(i);
0355       fDistAxonsoma[i] = fNeuronLoadParamz->GetDistAxonsoma(i);
0356       region->AddRootLogicalVolume(faxonLV[i]);
0357     }
0358   }
0359   // =======================================================================
0360   // Structure-4: Spines
0361   if (fNeuronLoadParamz->GetnbSpinecomp() == 0) {
0362     G4cout << " ---- Spines not found! ---- " << G4endl;
0363   }
0364   else {
0365     G4cout << " ---- Spines for construction: ---- " << G4endl;
0366     G4cout << " Total number of compartments into Spines : " << fNeuronLoadParamz->GetnbSpinecomp()
0367            << G4endl;
0368   }
0369
0370   G4cout << "\n ---- End of Neuron Construction! ---- "
0371          << "\n ========================================================== \n"
0372          << G4endl;
0373
0374   // =======================================================================
0375   // Active Geant4-DNA processes in BoundingSlice with whole neuron structure
0376   // fpRegion = new G4Region("BoundingSlice");
0377   // fpRegion->SetProductionCuts(cuts);
0378   // fpRegion->AddRootLogicalVolume(boundingLV);
0379
0380   return worldPV;
0381 }