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