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