neuron/src/NeuronLoadDataFile.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 //
0040 /// \file NeuronLoadDataFile.cc
0041 /// \brief Implementation of the NeuronLoadDataFile class
0042
0043 #include "NeuronLoadDataFile.hh"
0044
0045 #include "CommandLineParser.hh"
0046
0047 #include "G4Colour.hh"
0048 #include "G4LogicalVolume.hh"
0049 #include "G4PhysicalConstants.hh"
0050 #include "G4RotationMatrix.hh"
0051 #include "G4SystemOfUnits.hh"
0052 #include "G4UImanager.hh"
0053 #include "G4UnitsTable.hh"
0054 #include "G4VPhysicalVolume.hh"
0055 #include "G4VisAttributes.hh"
0056 #include "G4ios.hh"
0057 #include "globals.hh"
0058
0059 #include <algorithm>
0060 #include <fstream>
0061 #include <iostream>
0062 #include <limits>
0063 #include <sstream>
0064
0065 using namespace std;
0066 using namespace G4DNAPARSER;
0067
0068 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0069
0070 NeuronLoadDataFile::NeuronLoadDataFile()
0071 {
0072   Command* commandLine(nullptr);
0073
0074   // 1. Load single neuron morphology and obtain parameters.
0075   // Default SWC file name of neuron
0076   fNeuronFileNameSWC = "GranuleCell-Nr2.CNG.swc";
0077
0078   // 2. Load neural network and obtain parameters.
0079   // Default prepared data filename of neural network with single/multi-layer.
0080   // Small network of 10 pyramidal neurons with single layer
0081   fNeuronFileNameDATA = "NeuralNETWORK.dat";
0082
0083   // Load/change SWC or DAT as "CommandLineParser" class
0084   if ((commandLine = CommandLineParser::GetParser()->GetCommandIfActive("-swc"))) {
0085     fNeuronFileNameSWC = G4String(commandLine->GetOption());
0086   }
0087   if ((commandLine = CommandLineParser::GetParser()->GetCommandIfActive("-network"))) {
0088     auto ss = G4String(commandLine->GetOption());
0089     if ("" != ss) {
0090       fNeuronFileNameDATA = ss;
0091     }
0092     NeuralNetworkDATAfile(fNeuronFileNameDATA);
0093   }
0094   else {
0095     SingleNeuronSWCfile(fNeuronFileNameSWC);
0096   }
0097 }
0098
0099 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0100
0101 void NeuronLoadDataFile::SingleNeuronSWCfile(const G4String& filename)
0102 {
0103   // -----------
0104   // 12 November 2012 - code created
0105   // -------------------------------------------------------------------
0106   // November 2012: First model of neuron[*] adapted into Geant4 microdosimetry
0107   //                from Claiborne`s database[**] by M. Batmunkh.
0108   // February 2013: Loading SWC file from NeuronMorpho.Org[***]
0109   //                suggested by L. Bayarchimeg.
0110   // [*] http://lt-jds.jinr.ru/record/62124/files/lrb_e_2012.pdf
0111   // [**] http://www.utsa.edu/claibornelab/
0112   // [***] http://neuromorpho.org
0113   // -------------------------------------------------------------------
0114
0115   G4String sLine = "";
0116   std::ifstream infile;
0117   infile.open(filename.c_str());
0118   if (!infile.is_open()) {
0119     G4ExceptionDescription ed;
0120     ed << "Datafile " << filename << " is not opened!";
0121     G4Exception("NeuronLoadDataFile::SingleNeuronSWCfile()", "dna014", FatalException, ed,
0122                 "Check file path");
0123     return;
0124   }
0125   G4int nrows, nlines;
0126   nrows = 0;
0127   nlines = 0;
0128   for (;;) {
0129     getline(infile, sLine);
0130     if (infile.eof()) {
0131       break;
0132     }
0133     ++nrows;
0134   }
0135   infile.close();
0136
0137   G4cout << "NeuronLoadDataFile::SingleNeuronSWCfile: number of strings=" << nrows << " in file "
0138          << filename << G4endl;
0139
0140   infile.open(filename.c_str());
0141
0142   G4double TotVolSoma, TotVolDend, TotVolAxon, TotVolSpine;
0143   TotVolSoma = TotVolDend = TotVolAxon = TotVolSpine = 0.;
0144   G4double TotSurfSoma, TotSurfDend, TotSurfAxon, TotSurfSpine;
0145   TotSurfSoma = TotSurfDend = TotSurfAxon = TotSurfSpine = 0.;
0146   G4int nNcomp;  // current index of neuronal compartment
0147   G4int typeNcomp;  // type of neuron structures: soma, axon, dendrite, etc.
0148   G4double x, y, z;  // cartesian coordinates of each compartment in micrometer
0149   G4double radius;  // radius of each compartment in micrometer
0150   G4int pNcomp;  // linked compartment, indicates branch points of dendrites
0151   G4double minX, minY, minZ;
0152   G4double maxX, maxY, maxZ;
0153   G4double maxRad = -1e+09;
0154   minX = minY = minZ = 1e+09;
0155   maxX = maxY = maxZ = -1e+09;
0156   G4double density = 1.0 * (g / cm3);  // water medium
0157   G4double Piconst = (4.0 / 3.0) * pi;
0158
0159   fMassSomacomp.resize(nrows, 0);
0160   fPosSomacomp.resize(nrows);
0161   fRadSomacomp.resize(nrows, 0);
0162   std::vector<G4ThreeVector> PosDendcomp(nrows);
0163   fRadDendcomp.resize(nrows, 0);
0164   fHeightDendcomp.resize(nrows, 0);
0165   fMassDendcomp.resize(nrows, 0);
0166   fDistADendSoma.resize(nrows, 0);
0167   fDistBDendSoma.resize(nrows, 0);
0168   fPosDendcomp.resize(nrows);
0169   fRotDendcomp.resize(nrows);
0170   std::vector<G4ThreeVector> PosAxoncomp(nrows);
0171   fRadAxoncomp.resize(nrows, 0);
0172   fHeightAxoncomp.resize(nrows, 0);
0173   fMassAxoncomp.resize(nrows, 0);
0174   fDistAxonsoma.resize(nrows, 0);
0175   fPosAxoncomp.resize(nrows);
0176   fRotAxoncomp.resize(nrows);
0177   fMassSpinecomp.resize(nrows, 0);
0178   fPosSpinecomp.resize(nrows);
0179   fRadSpinecomp.resize(nrows, 0);
0180   fRadNeuroncomp.resize(nrows, 0);
0181   fHeightNeuroncomp.resize(nrows, 0);
0182   fDistNeuronsoma.resize(nrows, 0);
0183   fPosNeuroncomp.resize(nrows);
0184   fRotNeuroncomp.resize(nrows);
0185   fPosNeuroncomp.resize(nrows);
0186   fRadNeuroncomp.resize(nrows, 0);
0187   fTypeN.resize(nrows, 0);
0188   G4ThreeVector base;
0189
0190   // to read datafile containing numbers, alphabets and symbols..,
0191   for (;;) {
0192     getline(infile, sLine);
0193     if (infile.eof()) {
0194       break;
0195     }
0196     if ("#" == sLine.substr(0, 1)) {
0197       continue;
0198     };
0199
0200     std::istringstream form(sLine);
0201     form >> nNcomp >> typeNcomp >> x >> y >> z >> radius >> pNcomp;
0202     /*
0203   G4cout << "NeuronLoadDataFile::SingleNeuronSWCfile: typeNcomp="
0204          << typeNcomp << " nNcomp=" << nNcomp << " pNcomp=" << pNcomp << " N1="
0205          << fnbSomacomp << " N2=" << fnbAxoncomp << " N3=" << fnbDendritecomp << G4endl;
0206     */
0207     // =======================================================================
0208     // to find the largest and the smallest values of compartment positions
0209     // for parameters of bounding slice, sphere medium and shift of neuron.
0210     if (minX > x) minX = x;
0211     if (minY > y) minY = y;
0212     if (minZ > z) minZ = z;
0213     if (maxX < x) maxX = x;
0214     if (maxY < y) maxY = y;
0215     if (maxZ < z) maxZ = z;
0216     // max diameter of compartments
0217     if (maxRad < radius) maxRad = radius;
0218
0219     // =======================================================================
0220     // Soma compartments represented as Sphere or Ellipsoid solid
0221     if (typeNcomp == 1) {
0222       //  Sphere volume and surface area
0223       G4double VolSomacomp = Piconst * pow(radius * um, 3);
0224       TotVolSoma = TotVolSoma + VolSomacomp;
0225       G4double SurSomacomp = 3. * Piconst * pow(radius * um, 2);
0226       TotSurfSoma = TotSurfSoma + SurSomacomp;
0227       fMassSomacomp[fnbSomacomp] = density * VolSomacomp;
0228       fMassSomaTot = fMassSomaTot + fMassSomacomp[fnbSomacomp];
0229       G4ThreeVector vSoma(x, y, z);
0230       fPosSomacomp[fnbSomacomp] = vSoma;
0231       fRadSomacomp[fnbSomacomp] = radius;
0232       if (0 == fnbSomacomp) {
0233         base = G4ThreeVector(fRadSomacomp[0], fRadSomacomp[0], fRadSomacomp[0]);
0234       }
0235       ++fnbSomacomp;
0236     }
0237
0238     // =======================================================================
0239     // Apical and basal dendritic compartments represented as cylinderical solid
0240     if (typeNcomp == 3 || typeNcomp == 4) {
0241       G4ThreeVector vDend(x, y, z);
0242       // Position and Radius of compartments
0243       PosDendcomp[fnbDendritecomp] = vDend;
0244       fRadDendcomp[fnbDendritecomp] = radius;
0245       // To join two tracing points along the dendritic branches.
0246       // To calculate length, center and rotation angles of each cylinder
0247       fPosDendcomp[fnbDendritecomp] = vDend;  // translmDend;
0248       // delta of position A and position B of cylinder
0249       G4ThreeVector dend;
0250       // primary dendritic branch should be connect with Soma
0251       if (0 == fnbDendritecomp) {
0252         dend = PosDendcomp[fnbDendritecomp] - fPosSomacomp[0] - base;
0253       }
0254       else {
0255         dend = PosDendcomp[fnbDendritecomp] - PosDendcomp[fnbDendritecomp - 1];
0256       }
0257       // Height of compartment
0258       G4double lengthDendcomp = dend.mag();
0259       fHeightDendcomp[fnbDendritecomp] = lengthDendcomp;
0260
0261       // Distance from Soma
0262       G4ThreeVector dendDis = fPosSomacomp[0] - fPosDendcomp[fnbDendritecomp];
0263       if (typeNcomp == 3) fDistADendSoma[fnbDendritecomp] = dendDis.mag();
0264       if (typeNcomp == 4) fDistBDendSoma[fnbDendritecomp] = dendDis.mag();
0265
0266       //  Cylinder volume and surface area
0267       G4double VolDendcomp = pi * pow(radius * um, 2) * (lengthDendcomp * um);
0268       TotVolDend = TotVolDend + VolDendcomp;
0269       G4double SurDendcomp = 2. * pi * radius * um * (radius + lengthDendcomp) * um;
0270       TotSurfDend = TotSurfDend + SurDendcomp;
0271       fMassDendcomp[fnbDendritecomp] = density * VolDendcomp;
0272       fMassDendTot = fMassDendTot + fMassDendcomp[fnbDendritecomp];
0273
0274       dend = dend.unit();
0275
0276       // Euler angles of each compartment
0277       G4double theta_eulerDend = dend.theta();
0278       G4double phi_eulerDend = dend.phi();
0279       G4double psi_eulerDend = 0;
0280
0281       // Rotation Matrix, Euler constructor build inverse matrix.
0282       G4RotationMatrix rotmDendInv =
0283         G4RotationMatrix(phi_eulerDend + pi / 2, theta_eulerDend, psi_eulerDend);
0284       fRotDendcomp[fnbDendritecomp] = rotmDendInv.inverse();
0285       ++fnbDendritecomp;
0286     }
0287
0288     // =======================================================================
0289     // Axon compartments represented as cylinderical solid
0290     if (typeNcomp == 2 || typeNcomp == 7) {
0291       G4ThreeVector vAxon(x, y, z);
0292       // Position and Radius of compartments
0293       PosAxoncomp[fnbAxoncomp] = vAxon;
0294       fRadAxoncomp[fnbAxoncomp] = radius;
0295       // To join two tracing points in loaded SWC data file.
0296       // To calculate length, center and rotation angles of each cylinder
0297
0298       // delta of position A and position B of cylinder
0299       G4ThreeVector Axon;
0300       // primary axon point should be connect with Soma
0301       if (0 == fnbAxoncomp) {
0302         Axon = PosAxoncomp[fnbAxoncomp] - fPosSomacomp[0] - base;
0303       }
0304       else {
0305         Axon = PosAxoncomp[fnbAxoncomp] - PosAxoncomp[fnbAxoncomp - 1];
0306       }
0307       G4double lengthAxoncomp = Axon.mag();
0308       // Height of compartment
0309       fHeightAxoncomp[fnbAxoncomp] = lengthAxoncomp;
0310
0311       // Distance from Soma
0312       G4ThreeVector AxonDis = fPosSomacomp[0] - fPosAxoncomp[fnbAxoncomp];
0313       fDistAxonsoma[fnbAxoncomp] = AxonDis.mag();
0314
0315       //  Cylinder volume and surface area
0316       G4double VolAxoncomp = pi * pow(radius * um, 2) * (lengthAxoncomp * um);
0317       TotVolAxon = TotVolAxon + VolAxoncomp;
0318       G4double SurAxoncomp = 2. * pi * radius * um * (radius + lengthAxoncomp) * um;
0319       TotSurfAxon = TotSurfAxon + SurAxoncomp;
0320       fMassAxoncomp[fnbAxoncomp] = density * VolAxoncomp;
0321       fMassAxonTot += fMassAxoncomp[fnbAxoncomp];
0322       Axon = Axon.unit();
0323
0324       // Euler angles of each compartment
0325       G4double theta_eulerAxon = Axon.theta();
0326       G4double phi_eulerAxon = Axon.phi();
0327       G4double psi_eulerAxon = 0;
0328
0329       // Rotation Matrix, Euler constructor build inverse matrix.
0330       G4RotationMatrix rotmAxonInv =
0331         G4RotationMatrix(phi_eulerAxon + pi / 2, theta_eulerAxon, psi_eulerAxon);
0332       G4RotationMatrix rotmAxon = rotmAxonInv.inverse();
0333       fRotAxoncomp[fnbAxoncomp] = rotmAxon;
0334       ++fnbAxoncomp;
0335     }
0336     // =======================================================================
0337     // checking additional types
0338     if (typeNcomp != 1 && typeNcomp != 2 && typeNcomp != 3 && typeNcomp != 4) {
0339       G4cout << " Additional types:-->  " << typeNcomp << G4endl;
0340     }
0341
0342     // If tracing points including spines, user can be define spine morphology
0343     // including stubby, mushroom, thin, long thin, filopodia and
0344     // branched with heads and necks!
0345
0346     if (typeNcomp == 5) {
0347       //  Sphere volume and surface area
0348       G4double VolSpinecomp = Piconst * pow(radius * um, 3.);
0349       TotVolSpine = TotVolSpine + VolSpinecomp;
0350       G4double SurSpinecomp = 3. * Piconst * pow(radius * um, 2.);
0351       TotSurfSpine = TotSurfSpine + SurSpinecomp;
0352       fMassSpinecomp[fnbSpinecomp] = density * VolSpinecomp;
0353       fMassSpineTot = fMassSpineTot + fMassSpinecomp[fnbSpinecomp];
0354       // OR
0355       //  Ellipsoid volume and Approximate formula of surface area
0356       // ...
0357       G4ThreeVector vSpine(x, y, z);
0358       fPosSpinecomp[fnbSpinecomp] = vSpine;
0359       fRadSpinecomp[fnbSpinecomp] = radius;
0360       ++fnbSpinecomp;
0361     }
0362     ++nlines;
0363   }
0364   infile.close();
0365   // =======================================================================
0366
0367   fnbNeuroncomp = nlines;
0368   G4cout << " Total number of compartments into Neuron : " << fnbNeuroncomp << G4endl;
0369   G4cout << G4endl;
0370
0371   // to calculate SHIFT value for neuron translation
0372   fshiftX = (minX + maxX) / 2.;
0373   fshiftY = (minY + maxY) / 2.;
0374   fshiftZ = (minZ + maxZ) / 2.;
0375
0376   // width, height, depth of bounding slice volume
0377   fwidthB = std::fabs(minX - maxX) + maxRad;
0378   fheightB = std::fabs(minY - maxY) + maxRad;
0379   fdepthB = std::fabs(minZ - maxZ) + maxRad;
0380
0381   // diagonal length of bounding slice, that give diameter of sphere
0382   // for particle direction and fluence!
0383   fdiagnlLength = std::sqrt(fwidthB * fwidthB + fheightB * fheightB + fdepthB * fdepthB);
0384
0385   fTotVolNeuron = TotVolSoma + TotVolDend + TotVolAxon;
0386   fTotSurfNeuron = TotSurfSoma + TotSurfDend + TotSurfAxon;
0387   fTotMassNeuron = fMassSomaTot + fMassDendTot + fMassAxonTot;
0388
0389   fTotVolSlice = fwidthB * um * fheightB * um * fdepthB * um;
0390   fTotSurfSlice =
0391     2 * (fwidthB * um * fheightB * um + fheightB * um * fdepthB * um + fwidthB * um * fdepthB * um);
0392   fTotMassSlice = 1.0 * (g / cm3) * fTotVolSlice;
0393
0394   fTotVolMedium = Piconst * pow(fdiagnlLength * um / 2., 3.);
0395   fTotSurfMedium = 3. * Piconst * pow(fdiagnlLength * um / 2., 2);
0396   fTotMassMedium = 1.0 * (g / cm3) * fTotVolMedium;
0397 }
0398
0399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0400
0401 // Load prepared data file of neural network with single and multiple layers
0402 void NeuronLoadDataFile::NeuralNetworkDATAfile(const G4String& filename)
0403 {
0404   G4String sLine = "";
0405   std::ifstream infile;
0406   infile.open(filename.c_str());
0407   if (!infile.is_open()) {
0408     G4ExceptionDescription ed;
0409     ed << "Datafile " << filename << " is not opened!";
0410     G4Exception("NeuronLoadDataFile::NeuralNetworkDATAfile()", "dna014", FatalException, ed,
0411                 "Check file path");
0412     return;
0413   }
0414   G4cout << "NeuronLoadDataFile::NeuralNetworkDATAfile: opened " << filename << G4endl;
0415
0416   G4int nlines, nbSoma, nbDendrite;
0417   nlines = 0;
0418   fnbSomacomp = 0;  // total number of compartment into Soma
0419   fnbDendritecomp = 0;  // total number of compartment into Dendrites
0420   fnbAxoncomp = 0;  // total number of compartment into Axon
0421   fnbSpinecomp = 0;  // total number of compartment into Spines
0422   G4double TotVolSoma, TotVolDend, TotVolAxon;
0423   TotVolSoma = TotVolDend = TotVolAxon = 0.;
0424   G4double TotSurfSoma, TotSurfDend, TotSurfAxon;
0425   TotSurfSoma = TotSurfDend = TotSurfAxon = 0.;
0426   G4int typeNcomp;  // types of structure: soma, axon, apical dendrite, etc.
0427   G4double x1, y1, z1, x2, y2, z2;  // cartesian coordinates of each compartment
0428   G4double radius;  // radius of each compartment in micrometer
0429   G4double height;  // height of each compartment in micrometer
0430   G4double maxRad = -1e+09;
0431   G4double density = 1.0 * (g / cm3);  // water medium
0432   G4double Piconst = (4.0 / 3.0) * pi;
0433
0434   for (;;) {
0435     getline(infile, sLine);
0436     if (infile.eof()) {
0437       break;
0438     }
0439     std::istringstream form(sLine);
0440     if (nlines == 0) {
0441       // to read total number of compartments
0442       form >> fnbNeuroncomp >> nbSoma >> nbDendrite;
0443       fMassSomacomp.resize(nbSoma, 0);
0444       fPosSomacomp.resize(nbSoma);
0445       fRadSomacomp.resize(nbSoma, 0);
0446       fRadDendcomp.resize(nbDendrite, 0);
0447       fHeightDendcomp.resize(nbDendrite, 0);
0448       fMassDendcomp.resize(nbDendrite, 0);
0449       fDistADendSoma.resize(nbDendrite, 0);
0450       fDistBDendSoma.resize(nbDendrite, 0);
0451       fPosDendcomp.resize(nbDendrite);
0452       fRotDendcomp.resize(nbDendrite);
0453     }
0454     // =======================================================================
0455     // Soma compartments represented as Sphere or Ellipsoid solid
0456     if (nlines > 0 && nlines <= nbSoma) {
0457       form >> typeNcomp >> x1 >> y1 >> z1 >> radius;
0458       if (typeNcomp != 1) break;
0459       // max diameter of compartments
0460       if (maxRad < radius) maxRad = radius;
0461       //  Sphere volume and surface area
0462       G4double VolSomacomp = Piconst * pow(radius * um, 3.);
0463       TotVolSoma = TotVolSoma + VolSomacomp;
0464       G4double SurSomacomp = 3. * Piconst * pow(radius * um, 2.);
0465       TotSurfSoma = TotSurfSoma + SurSomacomp;
0466       fMassSomacomp[fnbSomacomp] = density * VolSomacomp;
0467       fMassSomaTot = fMassSomaTot + fMassSomacomp[fnbSomacomp];
0468
0469       G4ThreeVector vSoma(x1, y1, z1);
0470       fPosSomacomp[fnbSomacomp] = vSoma;
0471       fRadSomacomp[fnbSomacomp] = radius;
0472       ++fnbSomacomp;
0473     }
0474     // =======================================================================
0475     // Apical and basal dendritic compartments represented as cylinderical solid
0476     if (nlines > nbSoma && nlines <= fnbNeuroncomp) {
0477       form >> typeNcomp >> x1 >> y1 >> z1 >> x2 >> y2 >> z2 >> radius >> height;
0478       if (typeNcomp != 3) break;  // || typeNcomp != 4
0479
0480       // To calculate length, center and rotation angles of each cylinder
0481       // Center-position of each cylinder
0482       G4double Dendxx = x1 + x2;
0483       G4double Dendyy = y1 + y2;
0484       G4double Dendzz = z1 + z2;
0485       G4ThreeVector translmDend = G4ThreeVector(Dendxx / 2., Dendyy / 2., Dendzz / 2.);
0486       fPosDendcomp[fnbDendritecomp] = translmDend;
0487       fRadDendcomp[fnbDendritecomp] = radius;
0488       G4double lengthDendcomp = height;
0489       // Height of compartment
0490       fHeightDendcomp[fnbDendritecomp] = lengthDendcomp;
0491       // Distance from Soma
0492
0493       //  Cylinder volume and surface area
0494       G4double VolDendcomp = pi * pow(radius * um, 2) * (lengthDendcomp * um);
0495       TotVolDend = TotVolDend + VolDendcomp;
0496       G4double SurDendcomp = 2. * pi * radius * um * (radius + lengthDendcomp) * um;
0497       TotSurfDend = TotSurfDend + SurDendcomp;
0498       fMassDendcomp[fnbDendritecomp] = density * VolDendcomp;
0499       fMassDendTot = fMassDendTot + fMassDendcomp[fnbDendritecomp];
0500
0501       G4double Dendx = x1 - x2;
0502       G4double Dendy = y1 - y2;
0503       G4double Dendz = z1 - z2;
0504       Dendx = Dendx / lengthDendcomp;
0505       Dendy = Dendy / lengthDendcomp;
0506       Dendz = Dendz / lengthDendcomp;
0507
0508       // Euler angles of each compartment
0509       G4ThreeVector directionDend = G4ThreeVector(Dendx, Dendy, Dendz);
0510       G4double theta_eulerDend = directionDend.theta();
0511       G4double phi_eulerDend = directionDend.phi();
0512       G4double psi_eulerDend = 0;
0513
0514       // Rotation Matrix, Euler constructor build inverse matrix.
0515       G4RotationMatrix rotmDendInv =
0516         G4RotationMatrix(phi_eulerDend + pi / 2, theta_eulerDend, psi_eulerDend);
0517       G4RotationMatrix rotmDend = rotmDendInv.inverse();
0518
0519       fRotDendcomp[fnbDendritecomp] = rotmDend;
0520       ++fnbDendritecomp;
0521     }
0522     ++nlines;
0523   }
0524
0525   // =======================================================================
0526
0527   G4cout << " Total number of compartments into Neuron : " << fnbNeuroncomp << G4endl;
0528
0529   // to calculate SHIFT value for neuron translation
0530   fshiftX = 0.;  //(minX + maxX)/2. ;
0531   fshiftY = 0.;  //(minY + maxY)/2. ;
0532   fshiftZ = 0.;  //(minZ + maxZ)/2. ;
0533
0534   // width, height, depth of bounding slice volume
0535   fwidthB = 640.;
0536   fheightB = 280.;
0537   fdepthB = 25.;
0538   // diagonal length of bounding slice, that give diameter of sphere
0539   // for particle direction and fluence!
0540   fdiagnlLength = std::sqrt(fwidthB * fwidthB + fheightB * fheightB + fdepthB * fdepthB);
0541
0542   fTotVolNeuron = TotVolSoma + TotVolDend + TotVolAxon;
0543   fTotSurfNeuron = TotSurfSoma + TotSurfDend + TotSurfAxon;
0544   fTotMassNeuron = fMassSomaTot + fMassDendTot + fMassAxonTot;
0545
0546   fTotVolSlice = fwidthB * um * fheightB * um * fdepthB * um;
0547   fTotSurfSlice =
0548     2 * (fwidthB * um * fheightB * um + fheightB * um * fdepthB * um + fwidthB * um * fdepthB * um);
0549   fTotMassSlice = 1.0 * (g / cm3) * fTotVolSlice;
0550
0551   fTotVolMedium = Piconst * pow(fdiagnlLength * um / 2., 3.);
0552   fTotSurfMedium = 3. * Piconst * pow(fdiagnlLength * um / 2., 2);
0553   fTotMassMedium = 1.0 * (g / cm3) * fTotVolMedium;
0554
0555   infile.close();
0556 }
0557
0558 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0559
0560 void NeuronLoadDataFile::ComputeTransformation(const G4int copyNo, G4VPhysicalVolume* physVol) const
0561 {
0562   // to calculate Euler angles from Rotation Matrix after Inverse!
0563   //
0564   G4RotationMatrix rotmNeuron = G4RotationMatrix(fRotNeuroncomp[copyNo]);
0565   G4double cosX = std::sqrt(rotmNeuron.xx() * rotmNeuron.xx() + rotmNeuron.yx() * rotmNeuron.yx());
0566   G4double euX, euY, euZ;
0567   if (cosX > 16 * FLT_EPSILON) {
0568     euX = std::atan2(rotmNeuron.zy(), rotmNeuron.zz());
0569     euY = std::atan2(-rotmNeuron.zx(), cosX);
0570     euZ = std::atan2(rotmNeuron.yx(), rotmNeuron.xx());
0571   }
0572   else {
0573     euX = std::atan2(-rotmNeuron.yz(), rotmNeuron.yy());
0574     euY = std::atan2(-rotmNeuron.zx(), cosX);
0575     euZ = 0.;
0576   }
0577   G4RotationMatrix* rot = new G4RotationMatrix();
0578   rot->rotateX(euX);
0579   rot->rotateY(euY);
0580   rot->rotateZ(euZ);
0581
0582   physVol->SetRotation(rot);
0583
0584   // shift of cylinder compartments
0585   G4ThreeVector originNeuron((fPosNeuroncomp[copyNo].x() - fshiftX) * um,
0586                              (fPosNeuroncomp[copyNo].y() - fshiftY) * um,
0587                              (fPosNeuroncomp[copyNo].z() - fshiftZ) * um);
0588   physVol->SetTranslation(originNeuron);
0589 }
0590
0591 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0592
0593 void NeuronLoadDataFile::ComputeDimensions(G4Tubs& fcylinderComp, const G4int copyNo,
0594                                            const G4VPhysicalVolume*) const
0595 {
0596   fcylinderComp.SetInnerRadius(0 * um);
0597   fcylinderComp.SetOuterRadius(fRadNeuroncomp[copyNo] * um);
0598   fcylinderComp.SetZHalfLength(fHeightNeuroncomp[copyNo] * um / 2.);
0599   fcylinderComp.SetStartPhiAngle(0. * deg);
0600   fcylinderComp.SetDeltaPhiAngle(360. * deg);
0601 }
0602
0603 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......