ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
NeuronLoadDataFile.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file NeuronLoadDataFile.cc
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 //
42 
43 #include "NeuronLoadDataFile.hh"
44 //#include "NeuronLoadMessenger.hh"
45 #include "G4VPhysicalVolume.hh"
46 #include "G4LogicalVolume.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4UnitsTable.hh"
49 #include "G4PhysicalConstants.hh"
50 #include "G4Colour.hh"
51 #include "G4VisAttributes.hh"
52 #include "G4RotationMatrix.hh"
53 #include "G4ios.hh"
54 #include <algorithm>
55 #include <fstream>
56 #include <iostream>
57 #include <limits>
58 #include <cmath>
59 #include <sstream>
60 #include <string>
61 #include <stdlib.h>
62 //define if the program is running with Geant4
63 #define GEANT4
64 #ifdef GEANT4
65 //Specific to Geant4, globals.hh is used for G4cout
66 #include "globals.hh"
67 #endif
68 #include "CommandLineParser.hh"
69 #include "G4UImanager.hh"
70 
71 using namespace std;
72 using namespace G4DNAPARSER;
73 
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
75 
76 NeuronLoadDataFile::NeuronLoadDataFile()
77 {
78  //CommandLineParser* parser = CommandLineParser::GetParser();
79  Command* commandLine(0);
80 
81  // 1. Load single neuron morphology and obtain parameters.
82  // Default SWC file name of neuron
83  fNeuronFileNameSWC = G4String("GranuleCell-Nr2.CNG.swc");
84 
85  // 2. Load neural network and obtain parameters.
86  // Default prepared data filename of neural network with single/multi-layer.
87  // Small network of 10 pyramidal neurons with single layer
88  fNeuronFileNameDATA = G4String("NeuralNETWORK.dat");
89 
90  // Load/change SWC or DAT as "CommandLineParser" class
91  if((commandLine=CommandLineParser::GetParser()->GetCommandIfActive("-swc")))
92  {
93  fNeuronFileNameSWC = G4String(commandLine->GetOption());
94  SingleNeuronSWCfile(fNeuronFileNameSWC);
95  }
96  //if (CommandLineParser::GetParser()->GetCommandIfActive("-network"))
97  // {
98  // NeuralNetworkDATAfile(fNeuronFileNameDATA);
99  // }
100  if ((commandLine =CommandLineParser::GetParser()->
101  GetCommandIfActive("-network")))
102  {
103  fNeuronFileNameDATA = G4String(commandLine->GetOption());
104  NeuralNetworkDATAfile(fNeuronFileNameDATA);
105  }
106  else
107  {
108  SingleNeuronSWCfile(fNeuronFileNameSWC);
109  }
110 
111 }
112 
113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114 
115 void NeuronLoadDataFile::SingleNeuronSWCfile (const G4String& filename)
116 {
117 
118  // -----------
119  // 12 November 2012 - code created
120  // -------------------------------------------------------------------
121  // November 2012: First model of neuron[*] adapted into Geant4 microdosimetry
122  // from Claiborne`s database[**] by M. Batmunkh.
123  // February 2013: Loading SWC file from NeuronMorpho.Org[***]
124  // suggested by L. Bayarchimeg.
125  // [*] http://lt-jds.jinr.ru/record/62124/files/lrb_e_2012.pdf
126  // [**] http://www.utsa.edu/claibornelab/
127  // [***] http://neuromorpho.org
128  // -------------------------------------------------------------------
129 
130  G4String sLine = "";
131  std::ifstream infile;
132  infile.open(filename.c_str());
133  if (!infile)
134  {
135 #ifdef GEANT4
136  G4cout<<"\n NeuronLoadDataFile::SingleNeuronSWCfile >> datafile "
137  <<filename<<" not found !!!!"<<G4endl;
138  exit(0);
139 #endif
140  }
141  else
142  {
143 #ifdef G4VERBOSE
144  G4cout<<"\n NeuronLoadDataFile::SingleNeuronSWCfile >> opening filename: "
145  << "\n" <<'\t'<<'\t'<<'\t'<<'\t'<<'\t'<<"' "<<filename
146  << " ' \n"<< G4endl;
147 #endif
148 
149  G4int nrows,nlines;
150  nrows=0; nlines=0;
151  while (getline(infile, sLine))
152  {
153  fnNn = new G4int[nrows];
154  fpNn = new G4int[nrows];
155  fnNd = new G4int[nrows];
156  fpNd = new G4int[nrows];
157  fnNa = new G4int[nrows];
158  fpNa = new G4int[nrows];
159  nrows++;
160  }
161  infile.close();
162  //G4cout << " number of tracing points: "<< nrows <<G4endl;
163  infile.open(filename.c_str());
164 
165  fnbSomacomp = 0 ; // total number of compartment into Soma
166  fnbDendritecomp = 0 ; // total number of compartment into Dendrites
167  fnbAxoncomp = 0 ; // total number of compartment into Axon
168  fnbSpinecomp = 0 ; // total number of compartment into Spines
169  G4double TotVolSoma, TotVolDend, TotVolAxon, TotVolSpine;
170  TotVolSoma=TotVolDend=TotVolAxon=TotVolSpine=0.;
171  G4double TotSurfSoma, TotSurfDend, TotSurfAxon, TotSurfSpine;
172  TotSurfSoma=TotSurfDend=TotSurfAxon=TotSurfSpine=0.;
173  G4int nNcomp; // current index of neuronal compartment
174  G4int typeNcomp; // type of neuron structures: soma, axon, dendrite, etc.
175  G4double x,y,z; // cartesian coordinates of each compartment in micrometer
176  G4double radius; // radius of each compartment in micrometer
177  G4int pNcomp; // linked compartment, indicates branch points of dendrites
178  G4double minX,minY,minZ;
179  G4double maxX,maxY,maxZ;
180  G4double maxRad = -1e+09;
181  minX=minY=minZ=1e+09;
182  maxX=maxY=maxZ=-1e+09;
183  G4double density = 1.0 * (g/cm3) ; // water medium
184  G4double Piconst = (4.0/3.0)*pi ;
185 
186  fMassSomacomp = new G4double[nrows];
187  fMassSomaTot = 0.0 ;
188  fPosSomacomp = new G4ThreeVector[nrows];
189  fRadSomacomp = new G4double[nrows];
190  G4ThreeVector * PosDendcomp= new G4ThreeVector[nrows];
191  fRadDendcomp = new G4double[nrows];
192  fHeightDendcomp= new G4double[nrows];
193  fMassDendcomp = new G4double[nrows];
194  fMassDendTot = 0.0 ;
195  fDistADendSoma = new G4double[nrows];
196  fDistBDendSoma = new G4double[nrows];
197  fPosDendcomp = new G4ThreeVector[nrows];
198  fRotDendcomp = new G4RotationMatrix[nrows];
199  G4ThreeVector * PosAxoncomp= new G4ThreeVector[nrows];
200  fRadAxoncomp = new G4double[nrows];
201  fHeightAxoncomp= new G4double[nrows];
202  fMassAxoncomp = new G4double[nrows];
203  fMassAxonTot = 0.0 ;
204  fDistAxonsoma = new G4double[nrows];
205  fPosAxoncomp = new G4ThreeVector[nrows];
206  fRotAxoncomp = new G4RotationMatrix[nrows];
207  fMassSpinecomp = new G4double[nrows];
208  fMassSpineTot = 0.0 ;
209  fPosSpinecomp = new G4ThreeVector[nrows];
210  fRadSpinecomp = new G4double[nrows];
211  G4ThreeVector * PosNeuroncomp= new G4ThreeVector[nrows];
212  fRadNeuroncomp = new G4double[nrows];
213  fHeightNeuroncomp = new G4double[nrows];
214  fDistNeuronsoma = new G4double[nrows];
215  fPosNeuroncomp = new G4ThreeVector[nrows];
216  fRotNeuroncomp = new G4RotationMatrix[nrows];
217  fPosNeuroncomp = new G4ThreeVector[nrows];
218  fRadNeuroncomp = new G4double[nrows];
219  fTypeN = new G4int[nrows];
220 
221  // to read datafile containing numbers, alphabets and symbols..,
222  while (getline(infile, sLine))
223  {
224  std::istringstream form(sLine);
225  G4String token;
226  while (getline(form, token, ':'))
227  {
228  std::istringstream found(token);
229  while (found >> nNcomp >> typeNcomp >> x >> y >> z >> radius >> pNcomp)
230  {
231  // =======================================================================
232  // to find the largest and the smallest values of compartment positions
233  // for parameters of bounding slice, sphere medium and shift of neuron.
234  if (minX > x) minX = x;
235  if (minY > y) minY = y;
236  if (minZ > z) minZ = z;
237  if (maxX < x) maxX = x;
238  if (maxY < y) maxY = y;
239  if (maxZ < z) maxZ = z;
240  // max diameter of compartments
241  if (maxRad < radius) maxRad = radius;
242 
243  // =======================================================================
244  // Soma compartments represented as Sphere or Ellipsoid solid
245  if (typeNcomp == 1)
246  {
247  // Sphere volume and surface area
248  G4double VolSomacomp = Piconst*pow(radius*um,3.) ;
249  TotVolSoma = TotVolSoma + VolSomacomp;
250  G4double SurSomacomp = 3.*Piconst*pow(radius*um,2.) ;
251  TotSurfSoma = TotSurfSoma + SurSomacomp;
252  // OR
253  // Ellipsoid volume and Approximate formula of surface area
254  //G4double VolSomacomp = Piconst*(Ra*um)*(Rb*um)*(Rc*um);
255  //G4double SurSomacomp = 3.*Piconst*pow((pow(Ra,1.6075)*pow(Rb,1.6075)+
256  //pow(Ra,1.6075)*pow(Rc,1.6075)+pow(Rb,1.6075)*pow(Rc,1.6075))/3.,0.622084);
257  fMassSomacomp[fnbSomacomp] = density*VolSomacomp;
258  fMassSomaTot = fMassSomaTot + fMassSomacomp[fnbSomacomp];
259  G4ThreeVector vSoma (x ,y ,z);
260  fPosSomacomp [fnbSomacomp] = vSoma;
261  fRadSomacomp [fnbSomacomp]= radius;
262  // no rotate
263  // OR
264  // RotationMatrix for Ellipsoid solid
265  // ....
266  fnbSomacomp++ ;
267  }
268  // =======================================================================
269  // Apical and basal dendritic compartments represented as cylinderical solid
270  if (typeNcomp == 3 || typeNcomp == 4)
271  {
272  G4ThreeVector vDend (x ,y ,z);
273  // Position and Radius of compartments
274  PosDendcomp [fnbDendritecomp] = vDend;
275  fRadDendcomp [fnbDendritecomp]= radius;
276  fnNd[fnbDendritecomp]= nNcomp-(fnbSomacomp+fnbAxoncomp)-1;
277  fpNd[fnbDendritecomp]= pNcomp-(fnbSomacomp+fnbAxoncomp)-1;
278  // To join two tracing points along the dendritic branches.
279  // To calculate length, center and rotation angles of each cylinder
280 
281  // Center-position of each cylinder
282  G4double Dendxx= PosDendcomp[fnNd[fnbDendritecomp]].x()+
283  PosDendcomp[fpNd[fnbDendritecomp]].x();
284  G4double Dendyy= PosDendcomp[fnNd[fnbDendritecomp]].y()+
285  PosDendcomp[fpNd[fnbDendritecomp]].y();
286  G4double Dendzz= PosDendcomp[fnNd[fnbDendritecomp]].z()+
287  PosDendcomp[fpNd[fnbDendritecomp]].z();
288  G4ThreeVector translmDend = G4ThreeVector(Dendxx/2. ,
289  Dendyy/2. , Dendzz/2.) ;
290  fPosDendcomp [fnbDendritecomp] = translmDend;
291  // delta of position A and position B of cylinder
292  G4double Dendx, Dendy, Dendz;
293  //primary dendritic branch should be connect with Soma
294  if (fpNd[fnbDendritecomp] == -fnbSomacomp)
295  {
296  Dendx= PosDendcomp[fnNd[fnbDendritecomp]].x()-
297  (fPosSomacomp[0].x()+fRadSomacomp[0]);
298  Dendy= PosDendcomp[fnNd[fnbDendritecomp]].y()-
299  (fPosSomacomp[0].y()+fRadSomacomp[0]);
300  Dendz= PosDendcomp[fnNd[fnbDendritecomp]].z()-
301  (fPosSomacomp[0].z()+fRadSomacomp[0]);
302  }
303  else
304  {
305  Dendx= PosDendcomp[fnNd[fnbDendritecomp]].x()-
306  PosDendcomp[fpNd[fnbDendritecomp]].x();
307  Dendy= PosDendcomp[fnNd[fnbDendritecomp]].y()-
308  PosDendcomp[fpNd[fnbDendritecomp]].y();
309  Dendz= PosDendcomp[fnNd[fnbDendritecomp]].z()-
310  PosDendcomp[fpNd[fnbDendritecomp]].z();
311  }
312  G4double lengthDendcomp = std::sqrt(Dendx*Dendx+Dendy*Dendy+Dendz*Dendz);
313  // Height of compartment
314  fHeightDendcomp [fnbDendritecomp]= lengthDendcomp;
315 
316  // Distance from Soma
317  G4double DendDisx= fPosSomacomp[0].x()-
318  fPosDendcomp [fnbDendritecomp].x();
319  G4double DendDisy= fPosSomacomp[0].y()-
320  fPosDendcomp [fnbDendritecomp].y();
321  G4double DendDisz= fPosSomacomp[0].z()-
322  fPosDendcomp [fnbDendritecomp].z();
323  if (typeNcomp == 3) fDistADendSoma[fnbDendritecomp] =
324  std::sqrt(DendDisx*DendDisx + DendDisy*DendDisy + DendDisz*DendDisz);
325  if (typeNcomp == 4) fDistBDendSoma[fnbDendritecomp] =
326  std::sqrt(DendDisx*DendDisx + DendDisy*DendDisy + DendDisz*DendDisz);
327 
328  // Cylinder volume and surface area
329  G4double VolDendcomp = pi*pow(radius*um,2)*(lengthDendcomp*um);
330  TotVolDend = TotVolDend + VolDendcomp;
331  G4double SurDendcomp = 2.*pi*radius*um*(radius+lengthDendcomp)*um;
332  TotSurfDend = TotSurfDend + SurDendcomp;
333  fMassDendcomp[fnbDendritecomp] = density*VolDendcomp;
334  fMassDendTot = fMassDendTot + fMassDendcomp[fnbDendritecomp];
335 
336  Dendx=Dendx/lengthDendcomp;
337  Dendy=Dendy/lengthDendcomp;
338  Dendz=Dendz/lengthDendcomp;
339 
340  // Euler angles of each compartment
341  G4ThreeVector directionDend = G4ThreeVector(Dendx,Dendy,Dendz);
342  G4double theta_eulerDend = directionDend.theta();
343  G4double phi_eulerDend = directionDend.phi();
344  G4double psi_eulerDend = 0;
345 
346  //Rotation Matrix, Euler constructor build inverse matrix.
347  G4RotationMatrix rotmDendInv = G4RotationMatrix(
348  phi_eulerDend+pi/2,
349  theta_eulerDend,
350  psi_eulerDend);
351  G4RotationMatrix rotmDend = rotmDendInv.inverse();
352  /*
353  // To convert from Rotation Matrix after inverse to Euler angles
354  G4double cosX = std::sqrt (rotmDend.xx()*rotmDend.xx() +
355  rotmDend.yx()*rotmDend.yx()) ;
356  G4double euX, euY, euZ;
357  if (cosX > 16*FLT_EPSILON)
358  {
359  euX = std::atan2 (rotmDend.zy(),rotmDend.zz());
360  euY = std::atan2 (-rotmDend.zx(),cosX);
361  euZ = std::atan2 (rotmDend.yx(),rotmDend.xx());
362  }
363  else
364  {
365  euX = std::atan2 (-rotmDend.yz(),rotmDend.yy());
366  euY = std::atan2 (-rotmDend.zx(),cosX);
367  euZ = 0. ;
368  }
369  G4RotationMatrix * rot = new G4RotationMatrix();
370  rot->rotateX(euX);
371  rot->rotateY(euY);
372  rot->rotateZ(euZ);
373  */
374  fRotDendcomp [fnbDendritecomp]= rotmDend ;
375 
376  fnbDendritecomp++ ;
377  }
378 
379  // =======================================================================
380  // Axon compartments represented as cylinderical solid
381  if (typeNcomp == 2)
382  {
383  G4ThreeVector vAxon (x ,y ,z);
384  // Position and Radius of compartments
385  PosAxoncomp [fnbAxoncomp] = vAxon;
386  fRadAxoncomp [fnbAxoncomp]= radius;
387  fnNa[fnbAxoncomp]= nNcomp-(fnbSomacomp+fnbDendritecomp)-1;
388  fpNa[fnbAxoncomp]= pNcomp-(fnbSomacomp+fnbDendritecomp)-1;
389  // To join two tracing points in loaded SWC data file.
390  // To calculate length, center and rotation angles of each cylinder
391 
392  // Center-position of each cylinder
393  G4double Axonxx= PosAxoncomp[fnNa[fnbAxoncomp]].x()+
394  PosAxoncomp[fpNa[fnbAxoncomp]].x();
395  G4double Axonyy= PosAxoncomp[fnNa[fnbAxoncomp]].y()+
396  PosAxoncomp[fpNa[fnbAxoncomp]].y();
397  G4double Axonzz= PosAxoncomp[fnNa[fnbAxoncomp]].z()+
398  PosAxoncomp[fpNa[fnbAxoncomp]].z();
399  G4ThreeVector translmAxon = G4ThreeVector(Axonxx/2. ,
400  Axonyy/2. , Axonzz/2.) ;
401  fPosAxoncomp [fnbAxoncomp] = translmAxon;
402  // delta of position A and position B of cylinder
403  G4double Axonx, Axony, Axonz;
404  //primary axon point should be connect with Soma
405  if (fpNa[fnbAxoncomp] == -(fnbSomacomp+fnbDendritecomp))
406  {
407  Axonx= PosAxoncomp[fnNa[fnbAxoncomp]].x()-
408  (fPosSomacomp[0].x()+fRadSomacomp[0]);
409  Axony= PosAxoncomp[fnNa[fnbAxoncomp]].y()-
410  (fPosSomacomp[0].y()+fRadSomacomp[0]);
411  Axonz= PosAxoncomp[fnNa[fnbAxoncomp]].z()-
412  (fPosSomacomp[0].z()+fRadSomacomp[0]);
413  }
414  else
415  {
416  Axonx= PosAxoncomp[fnNa[fnbAxoncomp]].x()-
417  PosAxoncomp[fpNa[fnbAxoncomp]].x();
418  Axony= PosAxoncomp[fnNa[fnbAxoncomp]].y()-
419  PosAxoncomp[fpNa[fnbAxoncomp]].y();
420  Axonz= PosAxoncomp[fnNa[fnbAxoncomp]].z()-
421  PosAxoncomp[fpNa[fnbAxoncomp]].z();
422  }
423  G4double lengthAxoncomp = std::sqrt(Axonx*Axonx+Axony*Axony+Axonz*Axonz);
424  // Height of compartment
425  fHeightAxoncomp [fnbAxoncomp]= lengthAxoncomp;
426 
427  // Distance from Soma
428  G4double AxonDisx= fPosSomacomp[0].x()-
429  fPosAxoncomp [fnbAxoncomp].x();
430  G4double AxonDisy= fPosSomacomp[0].y()-
431  fPosAxoncomp [fnbAxoncomp].y();
432  G4double AxonDisz= fPosSomacomp[0].z()-
433  fPosAxoncomp [fnbAxoncomp].z();
434  fDistAxonsoma[fnbAxoncomp] = std::sqrt(AxonDisx*AxonDisx +
435  AxonDisy*AxonDisy + AxonDisz*AxonDisz);
436 
437  // Cylinder volume and surface area
438  G4double VolAxoncomp = pi*pow(radius*um,2)*(lengthAxoncomp*um);
439  TotVolAxon = TotVolAxon + VolAxoncomp;
440  G4double SurAxoncomp = 2.*pi*radius*um*(radius+lengthAxoncomp)*um;
441  TotSurfAxon = TotSurfAxon + SurAxoncomp;
442  fMassAxoncomp[fnbAxoncomp] = density*VolAxoncomp;
443  fMassAxonTot = fMassAxonTot + fMassAxoncomp[fnbAxoncomp];
444  Axonx=Axonx/lengthAxoncomp;
445  Axony=Axony/lengthAxoncomp;
446  Axonz=Axonz/lengthAxoncomp;
447 
448  // Euler angles of each compartment
449  G4ThreeVector directionAxon = G4ThreeVector(Axonx,Axony,Axonz);
450  G4double theta_eulerAxon = directionAxon.theta();
451  G4double phi_eulerAxon = directionAxon.phi();
452  G4double psi_eulerAxon = 0;
453 
454  //Rotation Matrix, Euler constructor build inverse matrix.
455  G4RotationMatrix rotmAxonInv = G4RotationMatrix(
456  phi_eulerAxon+pi/2,
457  theta_eulerAxon,
458  psi_eulerAxon);
459  G4RotationMatrix rotmAxon = rotmAxonInv.inverse();
460  fRotAxoncomp [fnbAxoncomp]= rotmAxon ;
461 
462  fnbAxoncomp++ ;
463  }
464  // =======================================================================
465  // checking additional types
466  if (typeNcomp != 1 && typeNcomp != 2 && typeNcomp != 3 && typeNcomp != 4)
467  {
468  G4cout << " Additional types:--> "<< typeNcomp <<G4endl;
469  }
470 
471  // If tracing points including spines, user can be define spine morphology
472  // including stubby, mushroom, thin, long thin, filopodia and
473  // branched with heads and necks!
474 
475  if (typeNcomp == 5)
476  {
477  // Sphere volume and surface area
478  G4double VolSpinecomp = Piconst*pow(radius*um,3.) ;
479  TotVolSpine = TotVolSpine + VolSpinecomp;
480  G4double SurSpinecomp = 3.*Piconst*pow(radius*um,2.) ;
481  TotSurfSpine = TotSurfSpine + SurSpinecomp;
482  fMassSpinecomp[fnbSpinecomp] = density*VolSpinecomp;
483  fMassSpineTot = fMassSpineTot + fMassSpinecomp[fnbSpinecomp];
484  // OR
485  // Ellipsoid volume and Approximate formula of surface area
486  // ...
487  G4ThreeVector vSpine (x ,y ,z);
488  fPosSpinecomp [fnbSpinecomp] = vSpine;
489  fRadSpinecomp [fnbSpinecomp] = radius;
490  // no rotate
491  // OR
492  // RotationMatrix for Ellipsoid solid
493  // ....
494  fnbSpinecomp++ ;
495  }
496 
497  // =======================================================================
498  // Neuron- all compartments allocate in ThreeVector
499  fTypeN[nlines] = typeNcomp;
500  G4ThreeVector vNeuron (x ,y ,z);
501  // Position and Radius of compartments
502  PosNeuroncomp [nlines] = vNeuron;
503  fRadNeuroncomp [nlines]= radius;
504  fnNn[nlines]= nNcomp-1;
505  fpNn[nlines]= pNcomp-1;
506  // To join two tracing points in loaded SWC data file.
507  // To calculate length, center and rotation angles of each cylinder
508 
509  // Center-position of each cylinder
510  G4double Neuronxx= PosNeuroncomp[fnNn[nlines]].x()+
511  PosNeuroncomp[fpNn[nlines]].x();
512  G4double Neuronyy= PosNeuroncomp[fnNn[nlines]].y()+
513  PosNeuroncomp[fpNn[nlines]].y();
514  G4double Neuronzz= PosNeuroncomp[fnNn[nlines]].z()+
515  PosNeuroncomp[fpNn[nlines]].z();
516  G4ThreeVector translmNeuron = G4ThreeVector(Neuronxx/2. ,
517  Neuronyy/2. , Neuronzz/2.) ;
518  fPosNeuroncomp [nlines] = translmNeuron;
519  // delta of position A and position B of cylinder
520  G4double Neuronx, Neurony, Neuronz;
521  //primary point
522  if (fpNn[nlines] == -2)
523  {
524  Neuronx= PosNeuroncomp[fnNn[nlines]].x()-
525  fPosNeuroncomp[0].x();
526  Neurony= PosNeuroncomp[fnNn[nlines]].y()-
527  fPosNeuroncomp[0].y();
528  Neuronz= PosNeuroncomp[fnNn[nlines]].z()-
529  fPosNeuroncomp[0].z();
530  }
531  else
532  {
533  Neuronx= PosNeuroncomp[fnNn[nlines]].x()-
534  PosNeuroncomp[fpNn[nlines]].x();
535  Neurony= PosNeuroncomp[fnNn[nlines]].y()-
536  PosNeuroncomp[fpNn[nlines]].y();
537  Neuronz= PosNeuroncomp[fnNn[nlines]].z()-
538  PosNeuroncomp[fpNn[nlines]].z();
539  }
540  G4double lengthNeuroncomp = std::sqrt(Neuronx*Neuronx+
541  Neurony*Neurony+Neuronz*Neuronz);
542  // Height of compartment
543  fHeightNeuroncomp [nlines]= lengthNeuroncomp;
544  // Distance from Soma
545  G4double NeuronDisx= fPosNeuroncomp[0].x()-
546  fPosNeuroncomp [nlines].x();
547  G4double NeuronDisy= fPosNeuroncomp[0].y()-
548  fPosNeuroncomp [nlines].y();
549  G4double NeuronDisz= fPosNeuroncomp[0].z()-
550  fPosNeuroncomp [nlines].z();
551  fDistNeuronsoma[nlines] = std::sqrt(NeuronDisx*NeuronDisx +
552  NeuronDisy*NeuronDisy + NeuronDisz*NeuronDisz);
553  /*
554  // Cylinder volume and surface area
555  G4double VolNeuroncomp = pi*pow(radius*um,2)*(lengthNeuroncomp*um);
556  fTotVolNeuron = fTotVolNeuron + VolNeuroncomp;
557  G4double SurNeuroncomp = 2.*pi*radius*um*
558  (radius+lengthNeuroncomp)*um;
559  fTotSurfNeuron = TotSurfNeuron + SurNeuroncomp;
560  fMassNeuroncomp[nlines] = density*VolNeuroncomp;
561  MassNeuronTot = MassNeuronTot + fMassNeuroncomp[nlines]; */
562  Neuronx=Neuronx/lengthNeuroncomp;
563  Neurony=Neurony/lengthNeuroncomp;
564  Neuronz=Neuronz/lengthNeuroncomp;
565 
566  // Euler angles of each compartment
567  G4ThreeVector directionNeuron = G4ThreeVector(Neuronx,Neurony,Neuronz);
568  G4double theta_eulerNeuron = directionNeuron.theta();
569  G4double phi_eulerNeuron = directionNeuron.phi();
570  G4double psi_eulerNeuron = 0;
571 
572  //Rotation Matrix, Euler constructor build inverse matrix.
573  G4RotationMatrix rotmNeuronInv = G4RotationMatrix(
574  phi_eulerNeuron+pi/2,
575  theta_eulerNeuron,
576  psi_eulerNeuron);
577  G4RotationMatrix rotmNeuron = rotmNeuronInv.inverse();
578  fRotNeuroncomp [nlines]= rotmNeuron ;
579 
580  nlines++;
581  }
582  }
583  }
584  infile.close();
585  // =======================================================================
586 
587  fnbNeuroncomp = nlines ;
588  G4cout << " Total number of compartments into Neuron : "
589  << fnbNeuroncomp<<G4endl;
590  G4cout << "\n"<<G4endl;
591 
592  // to calculate SHIFT value for neuron translation
593  fshiftX = (minX + maxX)/2. ;
594  fshiftY = (minY + maxY)/2. ;
595  fshiftZ = (minZ + maxZ)/2. ;
596 
597  // width, height, depth of bounding slice volume
598  //maxRad = 0.0 ;
599  fwidthB = std::fabs(minX - maxX) + maxRad;
600  fheightB = std::fabs(minY - maxY) + maxRad;
601  fdepthB = std::fabs(minZ - maxZ) + maxRad;
602 
603  // diagonal length of bounding slice, that give diameter of sphere
604  // for particle direction and fluence!
605  fdiagnlLength = std::sqrt(fwidthB*fwidthB + fheightB*fheightB
606  + fdepthB*fdepthB);
607 
608  fTotVolNeuron = TotVolSoma+TotVolDend+TotVolAxon;
609  fTotSurfNeuron = TotSurfSoma+TotSurfDend+TotSurfAxon;
610  fTotMassNeuron = fMassSomaTot+fMassDendTot+fMassAxonTot;
611 
612  fTotVolSlice = fwidthB*um*fheightB*um*fdepthB*um;
613  fTotSurfSlice = 2*(fwidthB*um*fheightB*um+fheightB*um*fdepthB*um+
614  fwidthB*um*fdepthB*um);
615  fTotMassSlice = 1.0 * (g/cm3) *fTotVolSlice;
616 
617  fTotVolMedium = Piconst*pow(fdiagnlLength*um/2.,3.) ;
618  fTotSurfMedium = 3.*Piconst*pow(fdiagnlLength*um/2.,2);
619  fTotMassMedium = 1.0 * (g/cm3) *fTotVolMedium;
620 
621  // Soma in Violet with opacity
622  fSomaColour = new G4VisAttributes;
623  fSomaColour->SetColour(G4Colour(G4Colour(0.85,0.44,0.84))); // ,1.0
624  fSomaColour->SetForceSolid(true); // true
625  fSomaColour->SetVisibility(true);
626 
627  // Dendrites in Dark-Blue
628  fDendColour = new G4VisAttributes;
629  fDendColour->SetColour(G4Colour(G4Colour(0.0, 0.0, 0.5)));
630  fDendColour->SetForceSolid(true);
631  //fDendColour->SetVisibility(true);
632 
633  // Axon in Maroon
634  fAxonColour = new G4VisAttributes;
635  fAxonColour->SetColour(G4Colour(G4Colour(0.5, 0.0, 0.0)));
636  fAxonColour->SetForceSolid(true);
637  fAxonColour->SetVisibility(true);
638 
639  // Spines in Dark-Green
640  fSpineColour = new G4VisAttributes;
641  fSpineColour->SetColour(G4Colour(G4Colour(0.0 , 100/255. , 0.0)));
642  fSpineColour->SetForceSolid(true);
643  fSpineColour->SetVisibility(true);
644 
645  // Whole neuron in semitransparent navy blue
646  fNeuronColour = new G4VisAttributes;
647  fNeuronColour->SetColour(G4Colour(G4Colour(0.0,0.4,0.8,0.5)));
648  fNeuronColour->SetForceSolid(true);
649  fNeuronColour->SetVisibility(true);
650 
651  }
652 }
653 
654 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
655 
656 // Load prepared data file of neural network with single and multiple layers
657 
658 void NeuronLoadDataFile::NeuralNetworkDATAfile
659  (const G4String& filename)
660 {
661 
662  G4String sLine = "";
663  std::ifstream infile;
664  infile.open(filename.c_str());
665  if (!infile)
666  {
667 #ifdef GEANT4
668  G4cout<<" \n NeuronLoadDataFile::NeuralNetworkDATAfile >> datafile "
669  <<filename<<" not found !!!!"<<G4endl;
670  exit(0);
671 #endif
672  }
673  else
674  {
675 #ifdef G4VERBOSE
676  G4cout<< " NeuronLoadDataFile::NeuralNetworkDATAfile >> opening filename: "
677  << "\n" <<'\t'<<'\t'<<'\t'<<'\t'<<'\t'<<"' "<<filename
678  << " ' \n"<< G4endl;
679 #endif
680 
681  G4int nlines, nbSoma, nbDendrite;
682  nlines=0;
683  fnbSomacomp = 0 ; // total number of compartment into Soma
684  fnbDendritecomp = 0 ; // total number of compartment into Dendrites
685  fnbAxoncomp = 0 ; // total number of compartment into Axon
686  fnbSpinecomp = 0 ; // total number of compartment into Spines
687  G4double TotVolSoma, TotVolDend, TotVolAxon;
688  TotVolSoma=TotVolDend=TotVolAxon=0.;
689  G4double TotSurfSoma, TotSurfDend, TotSurfAxon;
690  TotSurfSoma=TotSurfDend=TotSurfAxon=0.;
691  //G4int nNmorph; // current index of neuronal morphology
692  G4int typeNcomp; // types of structure: soma, axon, apical dendrite, etc.
693  G4double x1,y1,z1,x2,y2,z2; // cartesian coordinates of each compartment
694  G4double radius; // radius of each compartment in micrometer
695  G4double height; // height of each compartment in micrometer
696  //G4double minX,minY,minZ; //minimum
697  //G4double maxX,maxY,maxZ; //maximum
698  G4double maxRad = -1e+09;
699  //minX=minY=minZ=1e+09;
700  //maxX=maxY=maxZ=-1e+09;
701  G4double density = 1.0 * (g/cm3) ; // water medium
702  G4double Piconst = (4.0/3.0)*pi ;
703 
704  while (getline(infile, sLine))
705  {
706  std::istringstream form(sLine);
707  if (nlines == 0) {
708  // to read total number of compartments
709  form >> fnbNeuroncomp >> nbSoma >> nbDendrite ;
710  fMassSomacomp = new G4double[nbSoma];
711  fMassSomaTot = 0.0 ;
712  fPosSomacomp = new G4ThreeVector[nbSoma];
713  fRadSomacomp = new G4double[nbSoma];
714  fRadDendcomp = new G4double[nbDendrite];
715  fHeightDendcomp = new G4double[nbDendrite];
716  fMassDendcomp = new G4double[nbDendrite];
717  fMassDendTot = 0.0 ;
718  fDistADendSoma = new G4double[nbDendrite];
719  fDistBDendSoma = new G4double[nbDendrite];
720  fPosDendcomp = new G4ThreeVector[nbDendrite];
721  fRotDendcomp = new G4RotationMatrix[nbDendrite];
722  }
723  // =======================================================================
724  // Soma compartments represented as Sphere or Ellipsoid solid
725  if (nlines > 0 && nlines <= nbSoma) // Total number of Soma compartments
726  {
727  form >> typeNcomp >> x1 >> y1 >> z1 >> radius ;
728  if (typeNcomp !=1) break;
729  // max diameter of compartments
730  if (maxRad < radius) maxRad = radius;
731  // Sphere volume and surface area
732  G4double VolSomacomp = Piconst*pow(radius*um,3.) ;
733  TotVolSoma = TotVolSoma + VolSomacomp;
734  G4double SurSomacomp = 3.*Piconst*pow(radius*um,2.) ;
735  TotSurfSoma = TotSurfSoma + SurSomacomp;
736  fMassSomacomp[fnbSomacomp] = density*VolSomacomp;
737  fMassSomaTot = fMassSomaTot + fMassSomacomp[fnbSomacomp];
738  // OR
739  // Ellipsoid volume and Approximate formula of surface area
740  //G4double VolSomacomp = Piconst*(Ra*um)*(Rb*um)*(Rc*um);
741  //G4double SurSomacomp = 3.*Piconst*pow((pow(Ra,1.6075)*pow(Rb,1.6075)+
742  //pow(Ra,1.6075)*pow(Rc,1.6075)+pow(Rb,1.6075)*pow(Rc,1.6075))/3.,0.622084);
743 
744  G4ThreeVector vSoma (x1 ,y1 ,z1);
745  fPosSomacomp [fnbSomacomp] = vSoma;
746  fRadSomacomp [fnbSomacomp]= radius;
747 
748  // RotationMatrix for Ellipsoid solid
749  // ....
750  fnbSomacomp++ ;
751  }
752  // =======================================================================
753  // Apical and basal dendritic compartments represented as cylinderical solid
754  if (nlines > nbSoma && nlines <= fnbNeuroncomp)
755  {
756  form >> typeNcomp >> x1 >> y1 >> z1 >> x2 >> y2 >> z2 >> radius >> height;
757  if (typeNcomp != 3 ) break; // || typeNcomp != 4
758 
759  // To calculate length, center and rotation angles of each cylinder
760  // Center-position of each cylinder
761  G4double Dendxx= x1 + x2;
762  G4double Dendyy= y1 + y2;
763  G4double Dendzz= z1 + z2;
764  G4ThreeVector translmDend = G4ThreeVector(Dendxx/2. ,
765  Dendyy/2. , Dendzz/2.) ;
766  fPosDendcomp [fnbDendritecomp] = translmDend;
767  fRadDendcomp [fnbDendritecomp]= radius;
768  G4double lengthDendcomp = height;
769  // Height of compartment
770  fHeightDendcomp [fnbDendritecomp]= lengthDendcomp;
771  // Distance from Soma
772 
773  // Cylinder volume and surface area
774  G4double VolDendcomp = pi*pow(radius*um,2)*(lengthDendcomp*um);
775  TotVolDend = TotVolDend + VolDendcomp;
776  G4double SurDendcomp = 2.*pi*radius*um*(radius+lengthDendcomp)*um;
777  TotSurfDend = TotSurfDend + SurDendcomp;
778  fMassDendcomp[fnbDendritecomp] = density*VolDendcomp;
779  fMassDendTot = fMassDendTot + fMassDendcomp[fnbDendritecomp];
780 
781  G4double Dendx= x1 - x2;
782  G4double Dendy= y1 - y2;
783  G4double Dendz= z1 - z2;
784  Dendx=Dendx/lengthDendcomp;
785  Dendy=Dendy/lengthDendcomp;
786  Dendz=Dendz/lengthDendcomp;
787 
788  // Euler angles of each compartment
789  G4ThreeVector directionDend = G4ThreeVector(Dendx,Dendy,Dendz);
790  G4double theta_eulerDend = directionDend.theta();
791  G4double phi_eulerDend = directionDend.phi();
792  G4double psi_eulerDend = 0;
793 
794  //Rotation Matrix, Euler constructor build inverse matrix.
795  G4RotationMatrix rotmDendInv = G4RotationMatrix(
796  phi_eulerDend+pi/2,
797  theta_eulerDend,
798  psi_eulerDend);
799  G4RotationMatrix rotmDend = rotmDendInv.inverse();
800 
801  fRotDendcomp [fnbDendritecomp]= rotmDend ;
802  //G4Transform3D transformDend = G4Transform3D(rotmDend,translmDend);
803  //fRotTransDendPos [fnbDendritecomp]= transformDend ;
804  fnbDendritecomp++ ;
805 
806  }
807 
808  nlines++;
809  }
810 
811  // =======================================================================
812 
813  G4cout << " Total number of compartments into Neuron : "<<
814  fnbNeuroncomp <<G4endl;
815  G4cout << "\n"<<G4endl;
816 
817  // to calculate SHIFT value for neuron translation
818  fshiftX = 0.; //(minX + maxX)/2. ;
819  fshiftY = 0.; //(minY + maxY)/2. ;
820  fshiftZ = 0.; //(minZ + maxZ)/2. ;
821 
822  // width, height, depth of bounding slice volume
823  //maxRad = 0.0 ;
824  fwidthB = 640.;
825  fheightB = 280.;
826  fdepthB = 25.;
827  // diagonal length of bounding slice, that give diameter of sphere
828  // for particle direction and fluence!
829  fdiagnlLength = std::sqrt(fwidthB*fwidthB + fheightB*fheightB
830  + fdepthB*fdepthB);
831 
832  fTotVolNeuron = TotVolSoma+TotVolDend+TotVolAxon;
833  fTotSurfNeuron = TotSurfSoma+TotSurfDend+TotSurfAxon;
834  fTotMassNeuron = fMassSomaTot+fMassDendTot+fMassAxonTot;
835 
836  fTotVolSlice = fwidthB*um*fheightB*um*fdepthB*um;
837  fTotSurfSlice = 2*(fwidthB*um*fheightB*um+fheightB*um*fdepthB*um+
838  fwidthB*um*fdepthB*um);
839  fTotMassSlice = 1.0 * (g/cm3) *fTotVolSlice;
840 
841  fTotVolMedium = Piconst*pow(fdiagnlLength*um/2.,3.) ;
842  fTotSurfMedium = 3.*Piconst*pow(fdiagnlLength*um/2.,2);
843  fTotMassMedium = 1.0 * (g/cm3) *fTotVolMedium;
844 
845  }
846  infile.close();
847 }
848 
849 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
850 
851 NeuronLoadDataFile::~NeuronLoadDataFile()
852 {
853 delete[] fMassSomacomp ;
854 delete[] fPosSomacomp ;
855 delete[] fRadSomacomp ;
856 delete[] fRadDendcomp ;
857 delete[] fHeightDendcomp;
858 delete[] fMassDendcomp ;
859 delete[] fDistADendSoma ;
860 delete[] fDistBDendSoma ;
861 delete[] fPosDendcomp ;
862 delete[] fRotDendcomp ;
863 delete[] fRadAxoncomp ;
864 delete[] fHeightAxoncomp;
865 delete[] fMassAxoncomp ;
866 delete[] fDistAxonsoma ;
867 delete[] fPosAxoncomp ;
868 delete[] fRotAxoncomp ;
869 delete[] fRadNeuroncomp ;
870 delete[] fHeightNeuroncomp;
871 delete[] fMassNeuroncomp ;
872 delete[] fDistNeuronsoma ;
873 delete[] fPosNeuroncomp ;
874 delete[] fRotNeuroncomp ;
875 }
876 
877 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
878 
879 void NeuronLoadDataFile::ComputeTransformation
880 (const G4int copyNo, G4VPhysicalVolume* physVol) const
881 {
882 /*
883 // sphere transformation
884  G4ThreeVector
885  originSoma(
886  (fPosSomacomp[copyNo].x()-fshiftX) * um,
887  (fPosSomacomp[copyNo].y()-fshiftY) * um,
888  (fPosSomacomp[copyNo].z()-fshiftZ) * um
889  );
890  physVol->SetRotation(0);
891  physVol->SetTranslation(originSoma);
892  */
893 
894 // cylinder rotation and transformation
895 
896 // to calculate Euler angles from Rotation Matrix after Inverse!
897 //
898  G4RotationMatrix rotmNeuron = G4RotationMatrix(fRotNeuroncomp[copyNo]);
899  G4double cosX = std::sqrt (rotmNeuron.xx()*rotmNeuron.xx() +
900  rotmNeuron.yx()*rotmNeuron.yx()) ;
901  G4double euX, euY, euZ;
902  if (cosX > 16*FLT_EPSILON)
903  {
904  euX = std::atan2 (rotmNeuron.zy(),rotmNeuron.zz());
905  euY = std::atan2 (-rotmNeuron.zx(),cosX);
906  euZ = std::atan2 (rotmNeuron.yx(),rotmNeuron.xx());
907  }
908  else
909  {
910  euX = std::atan2 (-rotmNeuron.yz(),rotmNeuron.yy());
911  euY = std::atan2 (-rotmNeuron.zx(),cosX);
912  euZ = 0. ;
913  }
914  G4RotationMatrix* rot = new G4RotationMatrix();
915  rot->rotateX(euX);
916  rot->rotateY(euY);
917  rot->rotateZ(euZ);
918 
919  physVol->SetRotation(rot);
920 
921 // shift of cylinder compartments
922  G4ThreeVector
923  originNeuron(
924  (fPosNeuroncomp[copyNo].x()-fshiftX) * um,
925  (fPosNeuroncomp[copyNo].y()-fshiftY) * um,
926  (fPosNeuroncomp[copyNo].z()-fshiftZ) * um
927  );
928  physVol->SetTranslation(originNeuron);
929 
930 }
931 
932 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
933 /*
934 G4VSolid* NeuronLoadDataFile::ComputeSolid(const G4int copyNo,
935  G4VPhysicalVolume* physVol )
936 {
937  G4VSolid* solid;
938  if( typeNcomp[copyNo] == 1 )
939  {
940  solid = sphereComp ;
941  }
942  else if( typeNcomp[copyNo] == 3 || typeNcomp[copyNo] == 4 )
943  {
944  solid = cylinderComp ;
945  }
946  return solid;
947 }
948 */
949 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
950 
951 void NeuronLoadDataFile::ComputeDimensions
952 (G4Tubs& fcylinderComp, const G4int copyNo, const G4VPhysicalVolume*) const
953 {
954  fcylinderComp.SetInnerRadius(0*um);
955  fcylinderComp.SetOuterRadius(fRadNeuroncomp[copyNo]*um);
956  fcylinderComp.SetZHalfLength(fHeightNeuroncomp[copyNo]*um /2.);
957  fcylinderComp.SetStartPhiAngle(0.*deg);
958  fcylinderComp.SetDeltaPhiAngle(360.*deg);
959 }
960 
961 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
962 /*
963 void NeuronLoadDataFile::ComputeDimensions
964 (G4Sphere& sphereComp, const G4int copyNo, const G4VPhysicalVolume*) const
965 {
966  fsphereComp.SetInnerRadius(0);
967  fsphereComp.SetOuterRadius(fRadSomacomp[copyNo] * um);
968  fsphereComp.SetStartPhiAngle(0.*deg);
969  fsphereComp.SetDeltaPhiAngle(360.*deg);
970  fsphereComp.SetStartThetaAngle(0.*deg);
971  fsphereComp.SetDeltaThetaAngle(180.*deg);
972 }
973 */
974 /*
975 #if 1
976 (G4Sphere& somaS, const G4int copyNo, const G4VPhysicalVolume* physVol) const
977 #else
978 (G4Tubs& dendritesS, const G4int copyNo, const G4VPhysicalVolume* ) const
979 #endif
980 {
981  G4LogicalVolume* LogicalVolume = physVol->GetLogicalVolume();
982 
983  G4PhysicalVolume* somaPV = LogicalVolume->GetDaughter(0);
984  G4LogicalVolume* somaLV = somaPV->GetLogicalVolume();
985  G4Sphere* somaS = (G4Sphere*)somaLV->GetSolid();
986 
987  G4PhysicalVolume* dendritesPV = LogicalVolume->GetDaughter(0);
988  G4LogicalVolume* dendritesLV = dendritesPV->GetLogicalVolume();
989  G4Tubs* dendritesS = (G4Tubs*)dendritesLV->GetSolid();
990 }
991 */
992 
993 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......