Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/medical/dna/neuron/src/NeuronLoadDataFile.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 //
 39 //
 40 /// \file NeuronLoadDataFile.cc
 41 /// \brief Implementation of the NeuronLoadDataFile class
 42 
 43 #include "NeuronLoadDataFile.hh"
 44 
 45 #include "CommandLineParser.hh"
 46 
 47 #include "G4Colour.hh"
 48 #include "G4LogicalVolume.hh"
 49 #include "G4PhysicalConstants.hh"
 50 #include "G4RotationMatrix.hh"
 51 #include "G4SystemOfUnits.hh"
 52 #include "G4UImanager.hh"
 53 #include "G4UnitsTable.hh"
 54 #include "G4VPhysicalVolume.hh"
 55 #include "G4VisAttributes.hh"
 56 #include "G4ios.hh"
 57 #include "globals.hh"
 58 
 59 #include <algorithm>
 60 #include <fstream>
 61 #include <iostream>
 62 #include <limits>
 63 #include <sstream>
 64 
 65 using namespace std;
 66 using namespace G4DNAPARSER;
 67 
 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 69 
 70 NeuronLoadDataFile::NeuronLoadDataFile()
 71 {
 72   Command* commandLine(nullptr);
 73 
 74   // 1. Load single neuron morphology and obtain parameters.
 75   // Default SWC file name of neuron
 76   fNeuronFileNameSWC = "GranuleCell-Nr2.CNG.swc";
 77 
 78   // 2. Load neural network and obtain parameters.
 79   // Default prepared data filename of neural network with single/multi-layer.
 80   // Small network of 10 pyramidal neurons with single layer
 81   fNeuronFileNameDATA = "NeuralNETWORK.dat";
 82 
 83   // Load/change SWC or DAT as "CommandLineParser" class
 84   if ((commandLine = CommandLineParser::GetParser()->GetCommandIfActive("-swc"))) {
 85     fNeuronFileNameSWC = G4String(commandLine->GetOption());
 86   }
 87   if ((commandLine = CommandLineParser::GetParser()->GetCommandIfActive("-network"))) {
 88     auto ss = G4String(commandLine->GetOption());
 89     if ("" != ss) {
 90       fNeuronFileNameDATA = ss;
 91     }
 92     NeuralNetworkDATAfile(fNeuronFileNameDATA);
 93   }
 94   else {
 95     SingleNeuronSWCfile(fNeuronFileNameSWC);
 96   }
 97 }
 98 
 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100 
101 void NeuronLoadDataFile::SingleNeuronSWCfile(const G4String& filename)
102 {
103   // -----------
104   // 12 November 2012 - code created
105   // -------------------------------------------------------------------
106   // November 2012: First model of neuron[*] adapted into Geant4 microdosimetry
107   //                from Claiborne`s database[**] by M. Batmunkh.
108   // February 2013: Loading SWC file from NeuronMorpho.Org[***]
109   //                suggested by L. Bayarchimeg.
110   // [*] http://lt-jds.jinr.ru/record/62124/files/lrb_e_2012.pdf
111   // [**] http://www.utsa.edu/claibornelab/
112   // [***] http://neuromorpho.org
113   // -------------------------------------------------------------------
114 
115   G4String sLine = "";
116   std::ifstream infile;
117   infile.open(filename.c_str());
118   if (!infile.is_open()) {
119     G4ExceptionDescription ed;
120     ed << "Datafile " << filename << " is not opened!";
121     G4Exception("NeuronLoadDataFile::SingleNeuronSWCfile()", "dna014", FatalException, ed,
122                 "Check file path");
123     return;
124   }
125   G4int nrows, nlines;
126   nrows = 0;
127   nlines = 0;
128   for (;;) {
129     getline(infile, sLine);
130     if (infile.eof()) {
131       break;
132     }
133     ++nrows;
134   }
135   infile.close();
136 
137   G4cout << "NeuronLoadDataFile::SingleNeuronSWCfile: number of strings=" << nrows << " in file "
138          << filename << G4endl;
139 
140   infile.open(filename.c_str());
141 
142   G4double TotVolSoma, TotVolDend, TotVolAxon, TotVolSpine;
143   TotVolSoma = TotVolDend = TotVolAxon = TotVolSpine = 0.;
144   G4double TotSurfSoma, TotSurfDend, TotSurfAxon, TotSurfSpine;
145   TotSurfSoma = TotSurfDend = TotSurfAxon = TotSurfSpine = 0.;
146   G4int nNcomp;  // current index of neuronal compartment
147   G4int typeNcomp;  // type of neuron structures: soma, axon, dendrite, etc.
148   G4double x, y, z;  // cartesian coordinates of each compartment in micrometer
149   G4double radius;  // radius of each compartment in micrometer
150   G4int pNcomp;  // linked compartment, indicates branch points of dendrites
151   G4double minX, minY, minZ;
152   G4double maxX, maxY, maxZ;
153   G4double maxRad = -1e+09;
154   minX = minY = minZ = 1e+09;
155   maxX = maxY = maxZ = -1e+09;
156   G4double density = 1.0 * (g / cm3);  // water medium
157   G4double Piconst = (4.0 / 3.0) * pi;
158 
159   fMassSomacomp.resize(nrows, 0);
160   fPosSomacomp.resize(nrows);
161   fRadSomacomp.resize(nrows, 0);
162   std::vector<G4ThreeVector> PosDendcomp(nrows);
163   fRadDendcomp.resize(nrows, 0);
164   fHeightDendcomp.resize(nrows, 0);
165   fMassDendcomp.resize(nrows, 0);
166   fDistADendSoma.resize(nrows, 0);
167   fDistBDendSoma.resize(nrows, 0);
168   fPosDendcomp.resize(nrows);
169   fRotDendcomp.resize(nrows);
170   std::vector<G4ThreeVector> PosAxoncomp(nrows);
171   fRadAxoncomp.resize(nrows, 0);
172   fHeightAxoncomp.resize(nrows, 0);
173   fMassAxoncomp.resize(nrows, 0);
174   fDistAxonsoma.resize(nrows, 0);
175   fPosAxoncomp.resize(nrows);
176   fRotAxoncomp.resize(nrows);
177   fMassSpinecomp.resize(nrows, 0);
178   fPosSpinecomp.resize(nrows);
179   fRadSpinecomp.resize(nrows, 0);
180   fRadNeuroncomp.resize(nrows, 0);
181   fHeightNeuroncomp.resize(nrows, 0);
182   fDistNeuronsoma.resize(nrows, 0);
183   fPosNeuroncomp.resize(nrows);
184   fRotNeuroncomp.resize(nrows);
185   fPosNeuroncomp.resize(nrows);
186   fRadNeuroncomp.resize(nrows, 0);
187   fTypeN.resize(nrows, 0);
188   G4ThreeVector base;
189 
190   // to read datafile containing numbers, alphabets and symbols..,
191   for (;;) {
192     getline(infile, sLine);
193     if (infile.eof()) {
194       break;
195     }
196     if ("#" == sLine.substr(0, 1)) {
197       continue;
198     };
199 
200     std::istringstream form(sLine);
201     form >> nNcomp >> typeNcomp >> x >> y >> z >> radius >> pNcomp;
202     /*
203   G4cout << "NeuronLoadDataFile::SingleNeuronSWCfile: typeNcomp="
204          << typeNcomp << " nNcomp=" << nNcomp << " pNcomp=" << pNcomp << " N1="
205          << fnbSomacomp << " N2=" << fnbAxoncomp << " N3=" << fnbDendritecomp << G4endl;
206     */
207     // =======================================================================
208     // to find the largest and the smallest values of compartment positions
209     // for parameters of bounding slice, sphere medium and shift of neuron.
210     if (minX > x) minX = x;
211     if (minY > y) minY = y;
212     if (minZ > z) minZ = z;
213     if (maxX < x) maxX = x;
214     if (maxY < y) maxY = y;
215     if (maxZ < z) maxZ = z;
216     // max diameter of compartments
217     if (maxRad < radius) maxRad = radius;
218 
219     // =======================================================================
220     // Soma compartments represented as Sphere or Ellipsoid solid
221     if (typeNcomp == 1) {
222       //  Sphere volume and surface area
223       G4double VolSomacomp = Piconst * pow(radius * um, 3);
224       TotVolSoma = TotVolSoma + VolSomacomp;
225       G4double SurSomacomp = 3. * Piconst * pow(radius * um, 2);
226       TotSurfSoma = TotSurfSoma + SurSomacomp;
227       fMassSomacomp[fnbSomacomp] = density * VolSomacomp;
228       fMassSomaTot = fMassSomaTot + fMassSomacomp[fnbSomacomp];
229       G4ThreeVector vSoma(x, y, z);
230       fPosSomacomp[fnbSomacomp] = vSoma;
231       fRadSomacomp[fnbSomacomp] = radius;
232       if (0 == fnbSomacomp) {
233         base = G4ThreeVector(fRadSomacomp[0], fRadSomacomp[0], fRadSomacomp[0]);
234       }
235       ++fnbSomacomp;
236     }
237 
238     // =======================================================================
239     // Apical and basal dendritic compartments represented as cylinderical solid
240     if (typeNcomp == 3 || typeNcomp == 4) {
241       G4ThreeVector vDend(x, y, z);
242       // Position and Radius of compartments
243       PosDendcomp[fnbDendritecomp] = vDend;
244       fRadDendcomp[fnbDendritecomp] = radius;
245       // To join two tracing points along the dendritic branches.
246       // To calculate length, center and rotation angles of each cylinder
247       fPosDendcomp[fnbDendritecomp] = vDend;  // translmDend;
248       // delta of position A and position B of cylinder
249       G4ThreeVector dend;
250       // primary dendritic branch should be connect with Soma
251       if (0 == fnbDendritecomp) {
252         dend = PosDendcomp[fnbDendritecomp] - fPosSomacomp[0] - base;
253       }
254       else {
255         dend = PosDendcomp[fnbDendritecomp] - PosDendcomp[fnbDendritecomp - 1];
256       }
257       // Height of compartment
258       G4double lengthDendcomp = dend.mag();
259       fHeightDendcomp[fnbDendritecomp] = lengthDendcomp;
260 
261       // Distance from Soma
262       G4ThreeVector dendDis = fPosSomacomp[0] - fPosDendcomp[fnbDendritecomp];
263       if (typeNcomp == 3) fDistADendSoma[fnbDendritecomp] = dendDis.mag();
264       if (typeNcomp == 4) fDistBDendSoma[fnbDendritecomp] = dendDis.mag();
265 
266       //  Cylinder volume and surface area
267       G4double VolDendcomp = pi * pow(radius * um, 2) * (lengthDendcomp * um);
268       TotVolDend = TotVolDend + VolDendcomp;
269       G4double SurDendcomp = 2. * pi * radius * um * (radius + lengthDendcomp) * um;
270       TotSurfDend = TotSurfDend + SurDendcomp;
271       fMassDendcomp[fnbDendritecomp] = density * VolDendcomp;
272       fMassDendTot = fMassDendTot + fMassDendcomp[fnbDendritecomp];
273 
274       dend = dend.unit();
275 
276       // Euler angles of each compartment
277       G4double theta_eulerDend = dend.theta();
278       G4double phi_eulerDend = dend.phi();
279       G4double psi_eulerDend = 0;
280 
281       // Rotation Matrix, Euler constructor build inverse matrix.
282       G4RotationMatrix rotmDendInv =
283         G4RotationMatrix(phi_eulerDend + pi / 2, theta_eulerDend, psi_eulerDend);
284       fRotDendcomp[fnbDendritecomp] = rotmDendInv.inverse();
285       ++fnbDendritecomp;
286     }
287 
288     // =======================================================================
289     // Axon compartments represented as cylinderical solid
290     if (typeNcomp == 2 || typeNcomp == 7) {
291       G4ThreeVector vAxon(x, y, z);
292       // Position and Radius of compartments
293       PosAxoncomp[fnbAxoncomp] = vAxon;
294       fRadAxoncomp[fnbAxoncomp] = radius;
295       // To join two tracing points in loaded SWC data file.
296       // To calculate length, center and rotation angles of each cylinder
297 
298       // delta of position A and position B of cylinder
299       G4ThreeVector Axon;
300       // primary axon point should be connect with Soma
301       if (0 == fnbAxoncomp) {
302         Axon = PosAxoncomp[fnbAxoncomp] - fPosSomacomp[0] - base;
303       }
304       else {
305         Axon = PosAxoncomp[fnbAxoncomp] - PosAxoncomp[fnbAxoncomp - 1];
306       }
307       G4double lengthAxoncomp = Axon.mag();
308       // Height of compartment
309       fHeightAxoncomp[fnbAxoncomp] = lengthAxoncomp;
310 
311       // Distance from Soma
312       G4ThreeVector AxonDis = fPosSomacomp[0] - fPosAxoncomp[fnbAxoncomp];
313       fDistAxonsoma[fnbAxoncomp] = AxonDis.mag();
314 
315       //  Cylinder volume and surface area
316       G4double VolAxoncomp = pi * pow(radius * um, 2) * (lengthAxoncomp * um);
317       TotVolAxon = TotVolAxon + VolAxoncomp;
318       G4double SurAxoncomp = 2. * pi * radius * um * (radius + lengthAxoncomp) * um;
319       TotSurfAxon = TotSurfAxon + SurAxoncomp;
320       fMassAxoncomp[fnbAxoncomp] = density * VolAxoncomp;
321       fMassAxonTot += fMassAxoncomp[fnbAxoncomp];
322       Axon = Axon.unit();
323 
324       // Euler angles of each compartment
325       G4double theta_eulerAxon = Axon.theta();
326       G4double phi_eulerAxon = Axon.phi();
327       G4double psi_eulerAxon = 0;
328 
329       // Rotation Matrix, Euler constructor build inverse matrix.
330       G4RotationMatrix rotmAxonInv =
331         G4RotationMatrix(phi_eulerAxon + pi / 2, theta_eulerAxon, psi_eulerAxon);
332       G4RotationMatrix rotmAxon = rotmAxonInv.inverse();
333       fRotAxoncomp[fnbAxoncomp] = rotmAxon;
334       ++fnbAxoncomp;
335     }
336     // =======================================================================
337     // checking additional types
338     if (typeNcomp != 1 && typeNcomp != 2 && typeNcomp != 3 && typeNcomp != 4) {
339       G4cout << " Additional types:-->  " << typeNcomp << G4endl;
340     }
341 
342     // If tracing points including spines, user can be define spine morphology
343     // including stubby, mushroom, thin, long thin, filopodia and
344     // branched with heads and necks!
345 
346     if (typeNcomp == 5) {
347       //  Sphere volume and surface area
348       G4double VolSpinecomp = Piconst * pow(radius * um, 3.);
349       TotVolSpine = TotVolSpine + VolSpinecomp;
350       G4double SurSpinecomp = 3. * Piconst * pow(radius * um, 2.);
351       TotSurfSpine = TotSurfSpine + SurSpinecomp;
352       fMassSpinecomp[fnbSpinecomp] = density * VolSpinecomp;
353       fMassSpineTot = fMassSpineTot + fMassSpinecomp[fnbSpinecomp];
354       // OR
355       //  Ellipsoid volume and Approximate formula of surface area
356       // ...
357       G4ThreeVector vSpine(x, y, z);
358       fPosSpinecomp[fnbSpinecomp] = vSpine;
359       fRadSpinecomp[fnbSpinecomp] = radius;
360       ++fnbSpinecomp;
361     }
362     ++nlines;
363   }
364   infile.close();
365   // =======================================================================
366 
367   fnbNeuroncomp = nlines;
368   G4cout << " Total number of compartments into Neuron : " << fnbNeuroncomp << G4endl;
369   G4cout << G4endl;
370 
371   // to calculate SHIFT value for neuron translation
372   fshiftX = (minX + maxX) / 2.;
373   fshiftY = (minY + maxY) / 2.;
374   fshiftZ = (minZ + maxZ) / 2.;
375 
376   // width, height, depth of bounding slice volume
377   fwidthB = std::fabs(minX - maxX) + maxRad;
378   fheightB = std::fabs(minY - maxY) + maxRad;
379   fdepthB = std::fabs(minZ - maxZ) + maxRad;
380 
381   // diagonal length of bounding slice, that give diameter of sphere
382   // for particle direction and fluence!
383   fdiagnlLength = std::sqrt(fwidthB * fwidthB + fheightB * fheightB + fdepthB * fdepthB);
384 
385   fTotVolNeuron = TotVolSoma + TotVolDend + TotVolAxon;
386   fTotSurfNeuron = TotSurfSoma + TotSurfDend + TotSurfAxon;
387   fTotMassNeuron = fMassSomaTot + fMassDendTot + fMassAxonTot;
388 
389   fTotVolSlice = fwidthB * um * fheightB * um * fdepthB * um;
390   fTotSurfSlice =
391     2 * (fwidthB * um * fheightB * um + fheightB * um * fdepthB * um + fwidthB * um * fdepthB * um);
392   fTotMassSlice = 1.0 * (g / cm3) * fTotVolSlice;
393 
394   fTotVolMedium = Piconst * pow(fdiagnlLength * um / 2., 3.);
395   fTotSurfMedium = 3. * Piconst * pow(fdiagnlLength * um / 2., 2);
396   fTotMassMedium = 1.0 * (g / cm3) * fTotVolMedium;
397 }
398 
399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
400 
401 // Load prepared data file of neural network with single and multiple layers
402 void NeuronLoadDataFile::NeuralNetworkDATAfile(const G4String& filename)
403 {
404   G4String sLine = "";
405   std::ifstream infile;
406   infile.open(filename.c_str());
407   if (!infile.is_open()) {
408     G4ExceptionDescription ed;
409     ed << "Datafile " << filename << " is not opened!";
410     G4Exception("NeuronLoadDataFile::NeuralNetworkDATAfile()", "dna014", FatalException, ed,
411                 "Check file path");
412     return;
413   }
414   G4cout << "NeuronLoadDataFile::NeuralNetworkDATAfile: opened " << filename << G4endl;
415 
416   G4int nlines, nbSoma, nbDendrite;
417   nlines = 0;
418   fnbSomacomp = 0;  // total number of compartment into Soma
419   fnbDendritecomp = 0;  // total number of compartment into Dendrites
420   fnbAxoncomp = 0;  // total number of compartment into Axon
421   fnbSpinecomp = 0;  // total number of compartment into Spines
422   G4double TotVolSoma, TotVolDend, TotVolAxon;
423   TotVolSoma = TotVolDend = TotVolAxon = 0.;
424   G4double TotSurfSoma, TotSurfDend, TotSurfAxon;
425   TotSurfSoma = TotSurfDend = TotSurfAxon = 0.;
426   G4int typeNcomp;  // types of structure: soma, axon, apical dendrite, etc.
427   G4double x1, y1, z1, x2, y2, z2;  // cartesian coordinates of each compartment
428   G4double radius;  // radius of each compartment in micrometer
429   G4double height;  // height of each compartment in micrometer
430   G4double maxRad = -1e+09;
431   G4double density = 1.0 * (g / cm3);  // water medium
432   G4double Piconst = (4.0 / 3.0) * pi;
433 
434   for (;;) {
435     getline(infile, sLine);
436     if (infile.eof()) {
437       break;
438     }
439     std::istringstream form(sLine);
440     if (nlines == 0) {
441       // to read total number of compartments
442       form >> fnbNeuroncomp >> nbSoma >> nbDendrite;
443       fMassSomacomp.resize(nbSoma, 0);
444       fPosSomacomp.resize(nbSoma);
445       fRadSomacomp.resize(nbSoma, 0);
446       fRadDendcomp.resize(nbDendrite, 0);
447       fHeightDendcomp.resize(nbDendrite, 0);
448       fMassDendcomp.resize(nbDendrite, 0);
449       fDistADendSoma.resize(nbDendrite, 0);
450       fDistBDendSoma.resize(nbDendrite, 0);
451       fPosDendcomp.resize(nbDendrite);
452       fRotDendcomp.resize(nbDendrite);
453     }
454     // =======================================================================
455     // Soma compartments represented as Sphere or Ellipsoid solid
456     if (nlines > 0 && nlines <= nbSoma) {
457       form >> typeNcomp >> x1 >> y1 >> z1 >> radius;
458       if (typeNcomp != 1) break;
459       // max diameter of compartments
460       if (maxRad < radius) maxRad = radius;
461       //  Sphere volume and surface area
462       G4double VolSomacomp = Piconst * pow(radius * um, 3.);
463       TotVolSoma = TotVolSoma + VolSomacomp;
464       G4double SurSomacomp = 3. * Piconst * pow(radius * um, 2.);
465       TotSurfSoma = TotSurfSoma + SurSomacomp;
466       fMassSomacomp[fnbSomacomp] = density * VolSomacomp;
467       fMassSomaTot = fMassSomaTot + fMassSomacomp[fnbSomacomp];
468 
469       G4ThreeVector vSoma(x1, y1, z1);
470       fPosSomacomp[fnbSomacomp] = vSoma;
471       fRadSomacomp[fnbSomacomp] = radius;
472       ++fnbSomacomp;
473     }
474     // =======================================================================
475     // Apical and basal dendritic compartments represented as cylinderical solid
476     if (nlines > nbSoma && nlines <= fnbNeuroncomp) {
477       form >> typeNcomp >> x1 >> y1 >> z1 >> x2 >> y2 >> z2 >> radius >> height;
478       if (typeNcomp != 3) break;  // || typeNcomp != 4
479 
480       // To calculate length, center and rotation angles of each cylinder
481       // Center-position of each cylinder
482       G4double Dendxx = x1 + x2;
483       G4double Dendyy = y1 + y2;
484       G4double Dendzz = z1 + z2;
485       G4ThreeVector translmDend = G4ThreeVector(Dendxx / 2., Dendyy / 2., Dendzz / 2.);
486       fPosDendcomp[fnbDendritecomp] = translmDend;
487       fRadDendcomp[fnbDendritecomp] = radius;
488       G4double lengthDendcomp = height;
489       // Height of compartment
490       fHeightDendcomp[fnbDendritecomp] = lengthDendcomp;
491       // Distance from Soma
492 
493       //  Cylinder volume and surface area
494       G4double VolDendcomp = pi * pow(radius * um, 2) * (lengthDendcomp * um);
495       TotVolDend = TotVolDend + VolDendcomp;
496       G4double SurDendcomp = 2. * pi * radius * um * (radius + lengthDendcomp) * um;
497       TotSurfDend = TotSurfDend + SurDendcomp;
498       fMassDendcomp[fnbDendritecomp] = density * VolDendcomp;
499       fMassDendTot = fMassDendTot + fMassDendcomp[fnbDendritecomp];
500 
501       G4double Dendx = x1 - x2;
502       G4double Dendy = y1 - y2;
503       G4double Dendz = z1 - z2;
504       Dendx = Dendx / lengthDendcomp;
505       Dendy = Dendy / lengthDendcomp;
506       Dendz = Dendz / lengthDendcomp;
507 
508       // Euler angles of each compartment
509       G4ThreeVector directionDend = G4ThreeVector(Dendx, Dendy, Dendz);
510       G4double theta_eulerDend = directionDend.theta();
511       G4double phi_eulerDend = directionDend.phi();
512       G4double psi_eulerDend = 0;
513 
514       // Rotation Matrix, Euler constructor build inverse matrix.
515       G4RotationMatrix rotmDendInv =
516         G4RotationMatrix(phi_eulerDend + pi / 2, theta_eulerDend, psi_eulerDend);
517       G4RotationMatrix rotmDend = rotmDendInv.inverse();
518 
519       fRotDendcomp[fnbDendritecomp] = rotmDend;
520       ++fnbDendritecomp;
521     }
522     ++nlines;
523   }
524 
525   // =======================================================================
526 
527   G4cout << " Total number of compartments into Neuron : " << fnbNeuroncomp << G4endl;
528 
529   // to calculate SHIFT value for neuron translation
530   fshiftX = 0.;  //(minX + maxX)/2. ;
531   fshiftY = 0.;  //(minY + maxY)/2. ;
532   fshiftZ = 0.;  //(minZ + maxZ)/2. ;
533 
534   // width, height, depth of bounding slice volume
535   fwidthB = 640.;
536   fheightB = 280.;
537   fdepthB = 25.;
538   // diagonal length of bounding slice, that give diameter of sphere
539   // for particle direction and fluence!
540   fdiagnlLength = std::sqrt(fwidthB * fwidthB + fheightB * fheightB + fdepthB * fdepthB);
541 
542   fTotVolNeuron = TotVolSoma + TotVolDend + TotVolAxon;
543   fTotSurfNeuron = TotSurfSoma + TotSurfDend + TotSurfAxon;
544   fTotMassNeuron = fMassSomaTot + fMassDendTot + fMassAxonTot;
545 
546   fTotVolSlice = fwidthB * um * fheightB * um * fdepthB * um;
547   fTotSurfSlice =
548     2 * (fwidthB * um * fheightB * um + fheightB * um * fdepthB * um + fwidthB * um * fdepthB * um);
549   fTotMassSlice = 1.0 * (g / cm3) * fTotVolSlice;
550 
551   fTotVolMedium = Piconst * pow(fdiagnlLength * um / 2., 3.);
552   fTotSurfMedium = 3. * Piconst * pow(fdiagnlLength * um / 2., 2);
553   fTotMassMedium = 1.0 * (g / cm3) * fTotVolMedium;
554 
555   infile.close();
556 }
557 
558 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
559 
560 void NeuronLoadDataFile::ComputeTransformation(const G4int copyNo, G4VPhysicalVolume* physVol) const
561 {
562   // to calculate Euler angles from Rotation Matrix after Inverse!
563   //
564   G4RotationMatrix rotmNeuron = G4RotationMatrix(fRotNeuroncomp[copyNo]);
565   G4double cosX = std::sqrt(rotmNeuron.xx() * rotmNeuron.xx() + rotmNeuron.yx() * rotmNeuron.yx());
566   G4double euX, euY, euZ;
567   if (cosX > 16 * FLT_EPSILON) {
568     euX = std::atan2(rotmNeuron.zy(), rotmNeuron.zz());
569     euY = std::atan2(-rotmNeuron.zx(), cosX);
570     euZ = std::atan2(rotmNeuron.yx(), rotmNeuron.xx());
571   }
572   else {
573     euX = std::atan2(-rotmNeuron.yz(), rotmNeuron.yy());
574     euY = std::atan2(-rotmNeuron.zx(), cosX);
575     euZ = 0.;
576   }
577   G4RotationMatrix* rot = new G4RotationMatrix();
578   rot->rotateX(euX);
579   rot->rotateY(euY);
580   rot->rotateZ(euZ);
581 
582   physVol->SetRotation(rot);
583 
584   // shift of cylinder compartments
585   G4ThreeVector originNeuron((fPosNeuroncomp[copyNo].x() - fshiftX) * um,
586                              (fPosNeuroncomp[copyNo].y() - fshiftY) * um,
587                              (fPosNeuroncomp[copyNo].z() - fshiftZ) * um);
588   physVol->SetTranslation(originNeuron);
589 }
590 
591 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
592 
593 void NeuronLoadDataFile::ComputeDimensions(G4Tubs& fcylinderComp, const G4int copyNo,
594                                            const G4VPhysicalVolume*) const
595 {
596   fcylinderComp.SetInnerRadius(0 * um);
597   fcylinderComp.SetOuterRadius(fRadNeuroncomp[copyNo] * um);
598   fcylinderComp.SetZHalfLength(fHeightNeuroncomp[copyNo] * um / 2.);
599   fcylinderComp.SetStartPhiAngle(0. * deg);
600   fcylinderComp.SetDeltaPhiAngle(360. * deg);
601 }
602 
603 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
604