| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // This example is provided by the Geant4-DNA 26 // This example is provided by the Geant4-DNA collaboration
27 // Any report or published results obtained us 27 // Any report or published results obtained using the Geant4-DNA software
28 // shall cite the following Geant4-DNA collabo 28 // shall cite the following Geant4-DNA collaboration publication:
29 // Med. Phys. 37 (2010) 4692-4708 29 // Med. Phys. 37 (2010) 4692-4708
30 // and papers 30 // and papers
31 // M. Batmunkh et al. J Radiat Res Appl Sci 8 31 // M. Batmunkh et al. J Radiat Res Appl Sci 8 (2015) 498-507
32 // O. Belov et al. Physica Medica 32 (2016) 15 << 32 // O. Belov et al. Physica Medica 32 (2016) 1510-1520
33 // The Geant4-DNA web site is available at htt 33 // The Geant4-DNA web site is available at http://geant4-dna.org
34 // << 34 //
35 // ------------------------------------------- 35 // -------------------------------------------------------------------
36 // November 2016 36 // November 2016
37 // ------------------------------------------- 37 // -------------------------------------------------------------------
38 // 38 //
39 // << 39 //
40 /// \file NeuronLoadDataFile.cc 40 /// \file NeuronLoadDataFile.cc
41 /// \brief Implementation of the NeuronLoadDat 41 /// \brief Implementation of the NeuronLoadDataFile class
42 42
43 #include "NeuronLoadDataFile.hh" 43 #include "NeuronLoadDataFile.hh"
44 << 44 //#include "NeuronLoadMessenger.hh"
45 #include "CommandLineParser.hh" << 45 #include "G4VPhysicalVolume.hh"
46 <<
47 #include "G4Colour.hh" <<
48 #include "G4LogicalVolume.hh" 46 #include "G4LogicalVolume.hh"
49 #include "G4PhysicalConstants.hh" <<
50 #include "G4RotationMatrix.hh" <<
51 #include "G4SystemOfUnits.hh" 47 #include "G4SystemOfUnits.hh"
52 #include "G4UImanager.hh" <<
53 #include "G4UnitsTable.hh" 48 #include "G4UnitsTable.hh"
54 #include "G4VPhysicalVolume.hh" << 49 #include "G4PhysicalConstants.hh"
>> 50 #include "G4Colour.hh"
55 #include "G4VisAttributes.hh" 51 #include "G4VisAttributes.hh"
>> 52 #include "G4RotationMatrix.hh"
56 #include "G4ios.hh" 53 #include "G4ios.hh"
57 #include "globals.hh" << 54 #include <algorithm>
58 <<
59 #include <algorithm> <<
60 #include <fstream> 55 #include <fstream>
61 #include <iostream> 56 #include <iostream>
62 #include <limits> 57 #include <limits>
>> 58 #include <cmath>
63 #include <sstream> 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"
64 70
65 using namespace std; 71 using namespace std;
66 using namespace G4DNAPARSER; 72 using namespace G4DNAPARSER;
67 73
68 //....oooOO0OOooo........oooOO0OOooo........oo 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69 << 75
70 NeuronLoadDataFile::NeuronLoadDataFile() 76 NeuronLoadDataFile::NeuronLoadDataFile()
71 { 77 {
72 Command* commandLine(nullptr); << 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 }
73 110
74 // 1. Load single neuron morphology and obta <<
75 // Default SWC file name of neuron <<
76 fNeuronFileNameSWC = "GranuleCell-Nr2.CNG.sw <<
77 <<
78 // 2. Load neural network and obtain paramet <<
79 // Default prepared data filename of neural <<
80 // Small network of 10 pyramidal neurons wit <<
81 fNeuronFileNameDATA = "NeuralNETWORK.dat"; <<
82 <<
83 // Load/change SWC or DAT as "CommandLinePar <<
84 if ((commandLine = CommandLineParser::GetPar <<
85 fNeuronFileNameSWC = G4String(commandLine- <<
86 } <<
87 if ((commandLine = CommandLineParser::GetPar <<
88 auto ss = G4String(commandLine->GetOption( <<
89 if ("" != ss) { <<
90 fNeuronFileNameDATA = ss; <<
91 } <<
92 NeuralNetworkDATAfile(fNeuronFileNameDATA) <<
93 } <<
94 else { <<
95 SingleNeuronSWCfile(fNeuronFileNameSWC); <<
96 } <<
97 } 111 }
98 112
99 //....oooOO0OOooo........oooOO0OOooo........oo 113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100 114
101 void NeuronLoadDataFile::SingleNeuronSWCfile(c << 115 void NeuronLoadDataFile::SingleNeuronSWCfile (const G4String& filename)
102 { 116 {
>> 117
103 // ----------- 118 // -----------
104 // 12 November 2012 - code created 119 // 12 November 2012 - code created
105 // ----------------------------------------- 120 // -------------------------------------------------------------------
106 // November 2012: First model of neuron[*] a << 121 // November 2012: First model of neuron[*] adapted into Geant4 microdosimetry
107 // from Claiborne`s database[ 122 // from Claiborne`s database[**] by M. Batmunkh.
108 // February 2013: Loading SWC file from Neur << 123 // February 2013: Loading SWC file from NeuronMorpho.Org[***]
109 // suggested by L. Bayarchime 124 // suggested by L. Bayarchimeg.
110 // [*] http://lt-jds.jinr.ru/record/62124/fi 125 // [*] http://lt-jds.jinr.ru/record/62124/files/lrb_e_2012.pdf
111 // [**] http://www.utsa.edu/claibornelab/ 126 // [**] http://www.utsa.edu/claibornelab/
112 // [***] http://neuromorpho.org 127 // [***] http://neuromorpho.org
113 // ----------------------------------------- 128 // -------------------------------------------------------------------
114 129
115 G4String sLine = ""; 130 G4String sLine = "";
116 std::ifstream infile; 131 std::ifstream infile;
117 infile.open(filename.c_str()); 132 infile.open(filename.c_str());
118 if (!infile.is_open()) { << 133 if (!infile)
119 G4ExceptionDescription ed; << 134 {
120 ed << "Datafile " << filename << " is not << 135 #ifdef GEANT4
121 G4Exception("NeuronLoadDataFile::SingleNeu << 136 G4cout<<"\n NeuronLoadDataFile::SingleNeuronSWCfile >> datafile "
122 "Check file path"); << 137 <<filename<<" not found !!!!"<<G4endl;
123 return; << 138 exit(0);
>> 139 #endif
124 } 140 }
125 G4int nrows, nlines; << 141 else
126 nrows = 0; << 142 {
127 nlines = 0; << 143 #ifdef G4VERBOSE
128 for (;;) { << 144 G4cout<<"\n NeuronLoadDataFile::SingleNeuronSWCfile >> opening filename: "
129 getline(infile, sLine); << 145 << "\n" <<'\t'<<'\t'<<'\t'<<'\t'<<'\t'<<"' "<<filename
130 if (infile.eof()) { << 146 << " ' \n"<< G4endl;
131 break; << 147 #endif
132 } << 148
133 ++nrows; << 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++ ;
134 } 377 }
135 infile.close(); << 378
136 << 379 // =======================================================================
137 G4cout << "NeuronLoadDataFile::SingleNeuronS << 380 // Axon compartments represented as cylinderical solid
138 << filename << G4endl; << 381 if (typeNcomp == 2)
139 << 382 {
140 infile.open(filename.c_str()); << 383 G4ThreeVector vAxon (x ,y ,z);
141 << 384 // Position and Radius of compartments
142 G4double TotVolSoma, TotVolDend, TotVolAxon, << 385 PosAxoncomp [fnbAxoncomp] = vAxon;
143 TotVolSoma = TotVolDend = TotVolAxon = TotVo << 386 fRadAxoncomp [fnbAxoncomp]= radius;
144 G4double TotSurfSoma, TotSurfDend, TotSurfAx << 387 fnNa[fnbAxoncomp]= nNcomp-(fnbSomacomp+fnbDendritecomp)-1;
145 TotSurfSoma = TotSurfDend = TotSurfAxon = To << 388 fpNa[fnbAxoncomp]= pNcomp-(fnbSomacomp+fnbDendritecomp)-1;
146 G4int nNcomp; // current index of neuronal << 389 // To join two tracing points in loaded SWC data file.
147 G4int typeNcomp; // type of neuron structur << 390 // To calculate length, center and rotation angles of each cylinder
148 G4double x, y, z; // cartesian coordinates << 391
149 G4double radius; // radius of each compartm << 392 // Center-position of each cylinder
150 G4int pNcomp; // linked compartment, indica << 393 G4double Axonxx= PosAxoncomp[fnNa[fnbAxoncomp]].x()+
151 G4double minX, minY, minZ; << 394 PosAxoncomp[fpNa[fnbAxoncomp]].x();
152 G4double maxX, maxY, maxZ; << 395 G4double Axonyy= PosAxoncomp[fnNa[fnbAxoncomp]].y()+
153 G4double maxRad = -1e+09; << 396 PosAxoncomp[fpNa[fnbAxoncomp]].y();
154 minX = minY = minZ = 1e+09; << 397 G4double Axonzz= PosAxoncomp[fnNa[fnbAxoncomp]].z()+
155 maxX = maxY = maxZ = -1e+09; << 398 PosAxoncomp[fpNa[fnbAxoncomp]].z();
156 G4double density = 1.0 * (g / cm3); // wate << 399 G4ThreeVector translmAxon = G4ThreeVector(Axonxx/2. ,
157 G4double Piconst = (4.0 / 3.0) * pi; << 400 Axonyy/2. , Axonzz/2.) ;
158 << 401 fPosAxoncomp [fnbAxoncomp] = translmAxon;
159 fMassSomacomp.resize(nrows, 0); << 402 // delta of position A and position B of cylinder
160 fPosSomacomp.resize(nrows); << 403 G4double Axonx, Axony, Axonz;
161 fRadSomacomp.resize(nrows, 0); << 404 //primary axon point should be connect with Soma
162 std::vector<G4ThreeVector> PosDendcomp(nrows << 405 if (fpNa[fnbAxoncomp] == -(fnbSomacomp+fnbDendritecomp))
163 fRadDendcomp.resize(nrows, 0); << 406 {
164 fHeightDendcomp.resize(nrows, 0); << 407 Axonx= PosAxoncomp[fnNa[fnbAxoncomp]].x()-
165 fMassDendcomp.resize(nrows, 0); << 408 (fPosSomacomp[0].x()+fRadSomacomp[0]);
166 fDistADendSoma.resize(nrows, 0); << 409 Axony= PosAxoncomp[fnNa[fnbAxoncomp]].y()-
167 fDistBDendSoma.resize(nrows, 0); << 410 (fPosSomacomp[0].y()+fRadSomacomp[0]);
168 fPosDendcomp.resize(nrows); << 411 Axonz= PosAxoncomp[fnNa[fnbAxoncomp]].z()-
169 fRotDendcomp.resize(nrows); << 412 (fPosSomacomp[0].z()+fRadSomacomp[0]);
170 std::vector<G4ThreeVector> PosAxoncomp(nrows << 413 }
171 fRadAxoncomp.resize(nrows, 0); << 414 else
172 fHeightAxoncomp.resize(nrows, 0); << 415 {
173 fMassAxoncomp.resize(nrows, 0); << 416 Axonx= PosAxoncomp[fnNa[fnbAxoncomp]].x()-
174 fDistAxonsoma.resize(nrows, 0); << 417 PosAxoncomp[fpNa[fnbAxoncomp]].x();
175 fPosAxoncomp.resize(nrows); << 418 Axony= PosAxoncomp[fnNa[fnbAxoncomp]].y()-
176 fRotAxoncomp.resize(nrows); << 419 PosAxoncomp[fpNa[fnbAxoncomp]].y();
177 fMassSpinecomp.resize(nrows, 0); << 420 Axonz= PosAxoncomp[fnNa[fnbAxoncomp]].z()-
178 fPosSpinecomp.resize(nrows); << 421 PosAxoncomp[fpNa[fnbAxoncomp]].z();
179 fRadSpinecomp.resize(nrows, 0); << 422 }
180 fRadNeuroncomp.resize(nrows, 0); << 423 G4double lengthAxoncomp = std::sqrt(Axonx*Axonx+Axony*Axony+Axonz*Axonz);
181 fHeightNeuroncomp.resize(nrows, 0); << 424 // Height of compartment
182 fDistNeuronsoma.resize(nrows, 0); << 425 fHeightAxoncomp [fnbAxoncomp]= lengthAxoncomp;
183 fPosNeuroncomp.resize(nrows); << 426
184 fRotNeuroncomp.resize(nrows); << 427 // Distance from Soma
185 fPosNeuroncomp.resize(nrows); << 428 G4double AxonDisx= fPosSomacomp[0].x()-
186 fRadNeuroncomp.resize(nrows, 0); << 429 fPosAxoncomp [fnbAxoncomp].x();
187 fTypeN.resize(nrows, 0); << 430 G4double AxonDisy= fPosSomacomp[0].y()-
188 G4ThreeVector base; << 431 fPosAxoncomp [fnbAxoncomp].y();
189 << 432 G4double AxonDisz= fPosSomacomp[0].z()-
190 // to read datafile containing numbers, alph << 433 fPosAxoncomp [fnbAxoncomp].z();
191 for (;;) { << 434 fDistAxonsoma[fnbAxoncomp] = std::sqrt(AxonDisx*AxonDisx +
192 getline(infile, sLine); << 435 AxonDisy*AxonDisy + AxonDisz*AxonDisz);
193 if (infile.eof()) { << 436
194 break; << 437 // Cylinder volume and surface area
195 } << 438 G4double VolAxoncomp = pi*pow(radius*um,2)*(lengthAxoncomp*um);
196 if ("#" == sLine.substr(0, 1)) { << 439 TotVolAxon = TotVolAxon + VolAxoncomp;
197 continue; << 440 G4double SurAxoncomp = 2.*pi*radius*um*(radius+lengthAxoncomp)*um;
198 }; << 441 TotSurfAxon = TotSurfAxon + SurAxoncomp;
199 << 442 fMassAxoncomp[fnbAxoncomp] = density*VolAxoncomp;
200 std::istringstream form(sLine); << 443 fMassAxonTot = fMassAxonTot + fMassAxoncomp[fnbAxoncomp];
201 form >> nNcomp >> typeNcomp >> x >> y >> z << 444 Axonx=Axonx/lengthAxoncomp;
202 /* << 445 Axony=Axony/lengthAxoncomp;
203 G4cout << "NeuronLoadDataFile::SingleNeuronS << 446 Axonz=Axonz/lengthAxoncomp;
204 << typeNcomp << " nNcomp=" << nNcomp << 447
205 << fnbSomacomp << " N2=" << fnbAxonco << 448 // Euler angles of each compartment
206 */ << 449 G4ThreeVector directionAxon = G4ThreeVector(Axonx,Axony,Axonz);
207 // ======================================= << 450 G4double theta_eulerAxon = directionAxon.theta();
208 // to find the largest and the smallest va << 451 G4double phi_eulerAxon = directionAxon.phi();
209 // for parameters of bounding slice, spher << 452 G4double psi_eulerAxon = 0;
210 if (minX > x) minX = x; << 453
211 if (minY > y) minY = y; << 454 //Rotation Matrix, Euler constructor build inverse matrix.
212 if (minZ > z) minZ = z; << 455 G4RotationMatrix rotmAxonInv = G4RotationMatrix(
213 if (maxX < x) maxX = x; << 456 phi_eulerAxon+pi/2,
214 if (maxY < y) maxY = y; << 457 theta_eulerAxon,
215 if (maxZ < z) maxZ = z; << 458 psi_eulerAxon);
216 // max diameter of compartments << 459 G4RotationMatrix rotmAxon = rotmAxonInv.inverse();
217 if (maxRad < radius) maxRad = radius; << 460 fRotAxoncomp [fnbAxoncomp]= rotmAxon ;
218 << 461
219 // ======================================= << 462 fnbAxoncomp++ ;
220 // Soma compartments represented as Sphere << 463 }
221 if (typeNcomp == 1) { << 464 // =======================================================================
222 // Sphere volume and surface area << 465 // checking additional types
223 G4double VolSomacomp = Piconst * pow(rad << 466 if (typeNcomp != 1 && typeNcomp != 2 && typeNcomp != 3 && typeNcomp != 4)
224 TotVolSoma = TotVolSoma + VolSomacomp; << 467 {
225 G4double SurSomacomp = 3. * Piconst * po << 468 G4cout << " Additional types:--> "<< typeNcomp <<G4endl;
226 TotSurfSoma = TotSurfSoma + SurSomacomp; <<
227 fMassSomacomp[fnbSomacomp] = density * V <<
228 fMassSomaTot = fMassSomaTot + fMassSomac <<
229 G4ThreeVector vSoma(x, y, z); <<
230 fPosSomacomp[fnbSomacomp] = vSoma; <<
231 fRadSomacomp[fnbSomacomp] = radius; <<
232 if (0 == fnbSomacomp) { <<
233 base = G4ThreeVector(fRadSomacomp[0], <<
234 } <<
235 ++fnbSomacomp; <<
236 } <<
237 <<
238 // ======================================= <<
239 // Apical and basal dendritic compartments <<
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 <<
246 // To calculate length, center and rotat <<
247 fPosDendcomp[fnbDendritecomp] = vDend; <<
248 // delta of position A and position B of <<
249 G4ThreeVector dend; <<
250 // primary dendritic branch should be co <<
251 if (0 == fnbDendritecomp) { <<
252 dend = PosDendcomp[fnbDendritecomp] - <<
253 } <<
254 else { <<
255 dend = PosDendcomp[fnbDendritecomp] - <<
256 } <<
257 // Height of compartment <<
258 G4double lengthDendcomp = dend.mag(); <<
259 fHeightDendcomp[fnbDendritecomp] = lengt <<
260 <<
261 // Distance from Soma <<
262 G4ThreeVector dendDis = fPosSomacomp[0] <<
263 if (typeNcomp == 3) fDistADendSoma[fnbDe <<
264 if (typeNcomp == 4) fDistBDendSoma[fnbDe <<
265 <<
266 // Cylinder volume and surface area <<
267 G4double VolDendcomp = pi * pow(radius * <<
268 TotVolDend = TotVolDend + VolDendcomp; <<
269 G4double SurDendcomp = 2. * pi * radius <<
270 TotSurfDend = TotSurfDend + SurDendcomp; <<
271 fMassDendcomp[fnbDendritecomp] = density <<
272 fMassDendTot = fMassDendTot + fMassDendc <<
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 bu <<
282 G4RotationMatrix rotmDendInv = <<
283 G4RotationMatrix(phi_eulerDend + pi / <<
284 fRotDendcomp[fnbDendritecomp] = rotmDend <<
285 ++fnbDendritecomp; <<
286 } <<
287 <<
288 // ======================================= <<
289 // Axon compartments represented as cylind <<
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 <<
296 // To calculate length, center and rotat <<
297 <<
298 // delta of position A and position B of <<
299 G4ThreeVector Axon; <<
300 // primary axon point should be connect <<
301 if (0 == fnbAxoncomp) { <<
302 Axon = PosAxoncomp[fnbAxoncomp] - fPos <<
303 } <<
304 else { <<
305 Axon = PosAxoncomp[fnbAxoncomp] - PosA <<
306 } <<
307 G4double lengthAxoncomp = Axon.mag(); <<
308 // Height of compartment <<
309 fHeightAxoncomp[fnbAxoncomp] = lengthAxo <<
310 <<
311 // Distance from Soma <<
312 G4ThreeVector AxonDis = fPosSomacomp[0] <<
313 fDistAxonsoma[fnbAxoncomp] = AxonDis.mag <<
314 <<
315 // Cylinder volume and surface area <<
316 G4double VolAxoncomp = pi * pow(radius * <<
317 TotVolAxon = TotVolAxon + VolAxoncomp; <<
318 G4double SurAxoncomp = 2. * pi * radius <<
319 TotSurfAxon = TotSurfAxon + SurAxoncomp; <<
320 fMassAxoncomp[fnbAxoncomp] = density * V <<
321 fMassAxonTot += fMassAxoncomp[fnbAxoncom <<
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 bu <<
330 G4RotationMatrix rotmAxonInv = <<
331 G4RotationMatrix(phi_eulerAxon + pi / <<
332 G4RotationMatrix rotmAxon = rotmAxonInv. <<
333 fRotAxoncomp[fnbAxoncomp] = rotmAxon; <<
334 ++fnbAxoncomp; <<
335 } <<
336 // ======================================= <<
337 // checking additional types <<
338 if (typeNcomp != 1 && typeNcomp != 2 && ty <<
339 G4cout << " Additional types:--> " << t <<
340 } <<
341 <<
342 // If tracing points including spines, use <<
343 // including stubby, mushroom, thin, long <<
344 // branched with heads and necks! <<
345 <<
346 if (typeNcomp == 5) { <<
347 // Sphere volume and surface area <<
348 G4double VolSpinecomp = Piconst * pow(ra <<
349 TotVolSpine = TotVolSpine + VolSpinecomp <<
350 G4double SurSpinecomp = 3. * Piconst * p <<
351 TotSurfSpine = TotSurfSpine + SurSpineco <<
352 fMassSpinecomp[fnbSpinecomp] = density * <<
353 fMassSpineTot = fMassSpineTot + fMassSpi <<
354 // OR <<
355 // Ellipsoid volume and Approximate for <<
356 // ... <<
357 G4ThreeVector vSpine(x, y, z); <<
358 fPosSpinecomp[fnbSpinecomp] = vSpine; <<
359 fRadSpinecomp[fnbSpinecomp] = radius; <<
360 ++fnbSpinecomp; <<
361 } <<
362 ++nlines; <<
363 } 469 }
364 infile.close(); << 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
365 // ========================================= 497 // =======================================================================
366 << 498 // Neuron- all compartments allocate in ThreeVector
367 fnbNeuroncomp = nlines; << 499 fTypeN[nlines] = typeNcomp;
368 G4cout << " Total number of compartments int << 500 G4ThreeVector vNeuron (x ,y ,z);
369 G4cout << G4endl; << 501 // Position and Radius of compartments
370 << 502 PosNeuroncomp [nlines] = vNeuron;
371 // to calculate SHIFT value for neuron trans << 503 fRadNeuroncomp [nlines]= radius;
372 fshiftX = (minX + maxX) / 2.; << 504 fnNn[nlines]= nNcomp-1;
373 fshiftY = (minY + maxY) / 2.; << 505 fpNn[nlines]= pNcomp-1;
374 fshiftZ = (minZ + maxZ) / 2.; << 506 // To join two tracing points in loaded SWC data file.
375 << 507 // To calculate length, center and rotation angles of each cylinder
376 // width, height, depth of bounding slice vo << 508
377 fwidthB = std::fabs(minX - maxX) + maxRad; << 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;
378 fheightB = std::fabs(minY - maxY) + maxRad; 600 fheightB = std::fabs(minY - maxY) + maxRad;
379 fdepthB = std::fabs(minZ - maxZ) + maxRad; << 601 fdepthB = std::fabs(minZ - maxZ) + maxRad;
380 602
381 // diagonal length of bounding slice, that g << 603 // diagonal length of bounding slice, that give diameter of sphere
382 // for particle direction and fluence! << 604 // for particle direction and fluence!
383 fdiagnlLength = std::sqrt(fwidthB * fwidthB << 605 fdiagnlLength = std::sqrt(fwidthB*fwidthB + fheightB*fheightB
384 << 606 + fdepthB*fdepthB);
385 fTotVolNeuron = TotVolSoma + TotVolDend + To << 607
386 fTotSurfNeuron = TotSurfSoma + TotSurfDend + << 608 fTotVolNeuron = TotVolSoma+TotVolDend+TotVolAxon;
387 fTotMassNeuron = fMassSomaTot + fMassDendTot << 609 fTotSurfNeuron = TotSurfSoma+TotSurfDend+TotSurfAxon;
388 << 610 fTotMassNeuron = fMassSomaTot+fMassDendTot+fMassAxonTot;
389 fTotVolSlice = fwidthB * um * fheightB * um << 611
390 fTotSurfSlice = << 612 fTotVolSlice = fwidthB*um*fheightB*um*fdepthB*um;
391 2 * (fwidthB * um * fheightB * um + fheigh << 613 fTotSurfSlice = 2*(fwidthB*um*fheightB*um+fheightB*um*fdepthB*um+
392 fTotMassSlice = 1.0 * (g / cm3) * fTotVolSli << 614 fwidthB*um*fdepthB*um);
393 << 615 fTotMassSlice = 1.0 * (g/cm3) *fTotVolSlice;
394 fTotVolMedium = Piconst * pow(fdiagnlLength << 616
395 fTotSurfMedium = 3. * Piconst * pow(fdiagnlL << 617 fTotVolMedium = Piconst*pow(fdiagnlLength*um/2.,3.) ;
396 fTotMassMedium = 1.0 * (g / cm3) * fTotVolMe << 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 }
397 } 652 }
398 653
399 //....oooOO0OOooo........oooOO0OOooo........oo 654 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
400 655
401 // Load prepared data file of neural network w 656 // Load prepared data file of neural network with single and multiple layers
402 void NeuronLoadDataFile::NeuralNetworkDATAfile << 657
403 { << 658 void NeuronLoadDataFile::NeuralNetworkDATAfile
>> 659 (const G4String& filename)
>> 660 {
>> 661
404 G4String sLine = ""; 662 G4String sLine = "";
405 std::ifstream infile; 663 std::ifstream infile;
406 infile.open(filename.c_str()); 664 infile.open(filename.c_str());
407 if (!infile.is_open()) { << 665 if (!infile)
408 G4ExceptionDescription ed; << 666 {
409 ed << "Datafile " << filename << " is not << 667 #ifdef GEANT4
410 G4Exception("NeuronLoadDataFile::NeuralNet << 668 G4cout<<" \n NeuronLoadDataFile::NeuralNetworkDATAfile >> datafile "
411 "Check file path"); << 669 <<filename<<" not found !!!!"<<G4endl;
412 return; << 670 exit(0);
>> 671 #endif
413 } 672 }
414 G4cout << "NeuronLoadDataFile::NeuralNetwork << 673 else
415 << 674 {
416 G4int nlines, nbSoma, nbDendrite; << 675 #ifdef G4VERBOSE
417 nlines = 0; << 676 G4cout<< " NeuronLoadDataFile::NeuralNetworkDATAfile >> opening filename: "
418 fnbSomacomp = 0; // total number of compart << 677 << "\n" <<'\t'<<'\t'<<'\t'<<'\t'<<'\t'<<"' "<<filename
419 fnbDendritecomp = 0; // total number of com << 678 << " ' \n"<< G4endl;
420 fnbAxoncomp = 0; // total number of compart << 679 #endif
421 fnbSpinecomp = 0; // total number of compar << 680
422 G4double TotVolSoma, TotVolDend, TotVolAxon; << 681 G4int nlines, nbSoma, nbDendrite;
423 TotVolSoma = TotVolDend = TotVolAxon = 0.; << 682 nlines=0;
424 G4double TotSurfSoma, TotSurfDend, TotSurfAx << 683 fnbSomacomp = 0 ; // total number of compartment into Soma
425 TotSurfSoma = TotSurfDend = TotSurfAxon = 0. << 684 fnbDendritecomp = 0 ; // total number of compartment into Dendrites
426 G4int typeNcomp; // types of structure: som << 685 fnbAxoncomp = 0 ; // total number of compartment into Axon
427 G4double x1, y1, z1, x2, y2, z2; // cartesi << 686 fnbSpinecomp = 0 ; // total number of compartment into Spines
428 G4double radius; // radius of each compartm << 687 G4double TotVolSoma, TotVolDend, TotVolAxon;
429 G4double height; // height of each compartm << 688 TotVolSoma=TotVolDend=TotVolAxon=0.;
430 G4double maxRad = -1e+09; << 689 G4double TotSurfSoma, TotSurfDend, TotSurfAxon;
431 G4double density = 1.0 * (g / cm3); // wate << 690 TotSurfSoma=TotSurfDend=TotSurfAxon=0.;
432 G4double Piconst = (4.0 / 3.0) * pi; << 691 //G4int nNmorph; // current index of neuronal morphology
433 << 692 G4int typeNcomp; // types of structure: soma, axon, apical dendrite, etc.
434 for (;;) { << 693 G4double x1,y1,z1,x2,y2,z2; // cartesian coordinates of each compartment
435 getline(infile, sLine); << 694 G4double radius; // radius of each compartment in micrometer
436 if (infile.eof()) { << 695 G4double height; // height of each compartment in micrometer
437 break; << 696 //G4double minX,minY,minZ; //minimum
438 } << 697 //G4double maxX,maxY,maxZ; //maximum
439 std::istringstream form(sLine); << 698 G4double maxRad = -1e+09;
440 if (nlines == 0) { << 699 //minX=minY=minZ=1e+09;
441 // to read total number of compartments << 700 //maxX=maxY=maxZ=-1e+09;
442 form >> fnbNeuroncomp >> nbSoma >> nbDen << 701 G4double density = 1.0 * (g/cm3) ; // water medium
443 fMassSomacomp.resize(nbSoma, 0); << 702 G4double Piconst = (4.0/3.0)*pi ;
444 fPosSomacomp.resize(nbSoma); << 703
445 fRadSomacomp.resize(nbSoma, 0); << 704 while (getline(infile, sLine))
446 fRadDendcomp.resize(nbDendrite, 0); << 705 {
447 fHeightDendcomp.resize(nbDendrite, 0); << 706 std::istringstream form(sLine);
448 fMassDendcomp.resize(nbDendrite, 0); << 707 if (nlines == 0) {
449 fDistADendSoma.resize(nbDendrite, 0); << 708 // to read total number of compartments
450 fDistBDendSoma.resize(nbDendrite, 0); << 709 form >> fnbNeuroncomp >> nbSoma >> nbDendrite ;
451 fPosDendcomp.resize(nbDendrite); << 710 fMassSomacomp = new G4double[nbSoma];
452 fRotDendcomp.resize(nbDendrite); << 711 fMassSomaTot = 0.0 ;
453 } << 712 fPosSomacomp = new G4ThreeVector[nbSoma];
454 // ======================================= << 713 fRadSomacomp = new G4double[nbSoma];
455 // Soma compartments represented as Sphere << 714 fRadDendcomp = new G4double[nbDendrite];
456 if (nlines > 0 && nlines <= nbSoma) { << 715 fHeightDendcomp = new G4double[nbDendrite];
457 form >> typeNcomp >> x1 >> y1 >> z1 >> r << 716 fMassDendcomp = new G4double[nbDendrite];
458 if (typeNcomp != 1) break; << 717 fMassDendTot = 0.0 ;
459 // max diameter of compartments << 718 fDistADendSoma = new G4double[nbDendrite];
460 if (maxRad < radius) maxRad = radius; << 719 fDistBDendSoma = new G4double[nbDendrite];
461 // Sphere volume and surface area << 720 fPosDendcomp = new G4ThreeVector[nbDendrite];
462 G4double VolSomacomp = Piconst * pow(rad << 721 fRotDendcomp = new G4RotationMatrix[nbDendrite];
463 TotVolSoma = TotVolSoma + VolSomacomp; << 722 }
464 G4double SurSomacomp = 3. * Piconst * po << 723 // =======================================================================
465 TotSurfSoma = TotSurfSoma + SurSomacomp; << 724 // Soma compartments represented as Sphere or Ellipsoid solid
466 fMassSomacomp[fnbSomacomp] = density * V << 725 if (nlines > 0 && nlines <= nbSoma) // Total number of Soma compartments
467 fMassSomaTot = fMassSomaTot + fMassSomac << 726 {
468 << 727 form >> typeNcomp >> x1 >> y1 >> z1 >> radius ;
469 G4ThreeVector vSoma(x1, y1, z1); << 728 if (typeNcomp !=1) break;
470 fPosSomacomp[fnbSomacomp] = vSoma; << 729 // max diameter of compartments
471 fRadSomacomp[fnbSomacomp] = radius; << 730 if (maxRad < radius) maxRad = radius;
472 ++fnbSomacomp; << 731 // Sphere volume and surface area
473 } << 732 G4double VolSomacomp = Piconst*pow(radius*um,3.) ;
474 // ======================================= << 733 TotVolSoma = TotVolSoma + VolSomacomp;
475 // Apical and basal dendritic compartments << 734 G4double SurSomacomp = 3.*Piconst*pow(radius*um,2.) ;
476 if (nlines > nbSoma && nlines <= fnbNeuron << 735 TotSurfSoma = TotSurfSoma + SurSomacomp;
477 form >> typeNcomp >> x1 >> y1 >> z1 >> x << 736 fMassSomacomp[fnbSomacomp] = density*VolSomacomp;
478 if (typeNcomp != 3) break; // || typeNc << 737 fMassSomaTot = fMassSomaTot + fMassSomacomp[fnbSomacomp];
479 << 738 // OR
480 // To calculate length, center and rotat << 739 // Ellipsoid volume and Approximate formula of surface area
481 // Center-position of each cylinder << 740 //G4double VolSomacomp = Piconst*(Ra*um)*(Rb*um)*(Rc*um);
482 G4double Dendxx = x1 + x2; << 741 //G4double SurSomacomp = 3.*Piconst*pow((pow(Ra,1.6075)*pow(Rb,1.6075)+
483 G4double Dendyy = y1 + y2; << 742 //pow(Ra,1.6075)*pow(Rc,1.6075)+pow(Rb,1.6075)*pow(Rc,1.6075))/3.,0.622084);
484 G4double Dendzz = z1 + z2; << 743
485 G4ThreeVector translmDend = G4ThreeVecto << 744 G4ThreeVector vSoma (x1 ,y1 ,z1);
486 fPosDendcomp[fnbDendritecomp] = translmD << 745 fPosSomacomp [fnbSomacomp] = vSoma;
487 fRadDendcomp[fnbDendritecomp] = radius; << 746 fRadSomacomp [fnbSomacomp]= radius;
488 G4double lengthDendcomp = height; << 747
489 // Height of compartment << 748 // RotationMatrix for Ellipsoid solid
490 fHeightDendcomp[fnbDendritecomp] = lengt << 749 // ....
491 // Distance from Soma << 750 fnbSomacomp++ ;
492 <<
493 // Cylinder volume and surface area <<
494 G4double VolDendcomp = pi * pow(radius * <<
495 TotVolDend = TotVolDend + VolDendcomp; <<
496 G4double SurDendcomp = 2. * pi * radius <<
497 TotSurfDend = TotSurfDend + SurDendcomp; <<
498 fMassDendcomp[fnbDendritecomp] = density <<
499 fMassDendTot = fMassDendTot + fMassDendc <<
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 = G4ThreeVec <<
510 G4double theta_eulerDend = directionDend <<
511 G4double phi_eulerDend = directionDend.p <<
512 G4double psi_eulerDend = 0; <<
513 <<
514 // Rotation Matrix, Euler constructor bu <<
515 G4RotationMatrix rotmDendInv = <<
516 G4RotationMatrix(phi_eulerDend + pi / <<
517 G4RotationMatrix rotmDend = rotmDendInv. <<
518 <<
519 fRotDendcomp[fnbDendritecomp] = rotmDend <<
520 ++fnbDendritecomp; <<
521 } <<
522 ++nlines; <<
523 } 751 }
524 <<
525 // ========================================= 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 }
526 848
527 G4cout << " Total number of compartments int << 849 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
528 <<
529 // to calculate SHIFT value for neuron trans <<
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 vo <<
535 fwidthB = 640.; <<
536 fheightB = 280.; <<
537 fdepthB = 25.; <<
538 // diagonal length of bounding slice, that g <<
539 // for particle direction and fluence! <<
540 fdiagnlLength = std::sqrt(fwidthB * fwidthB <<
541 <<
542 fTotVolNeuron = TotVolSoma + TotVolDend + To <<
543 fTotSurfNeuron = TotSurfSoma + TotSurfDend + <<
544 fTotMassNeuron = fMassSomaTot + fMassDendTot <<
545 <<
546 fTotVolSlice = fwidthB * um * fheightB * um <<
547 fTotSurfSlice = <<
548 2 * (fwidthB * um * fheightB * um + fheigh <<
549 fTotMassSlice = 1.0 * (g / cm3) * fTotVolSli <<
550 <<
551 fTotVolMedium = Piconst * pow(fdiagnlLength <<
552 fTotSurfMedium = 3. * Piconst * pow(fdiagnlL <<
553 fTotMassMedium = 1.0 * (g / cm3) * fTotVolMe <<
554 850
555 infile.close(); << 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 ;
556 } 875 }
557 876
558 //....oooOO0OOooo........oooOO0OOooo........oo 877 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
559 878
560 void NeuronLoadDataFile::ComputeTransformation << 879 void NeuronLoadDataFile::ComputeTransformation
>> 880 (const G4int copyNo, G4VPhysicalVolume* physVol) const
561 { 881 {
562 // to calculate Euler angles from Rotation M << 882 /*
563 // << 883 // sphere transformation
564 G4RotationMatrix rotmNeuron = G4RotationMatr << 884 G4ThreeVector
565 G4double cosX = std::sqrt(rotmNeuron.xx() * << 885 originSoma(
566 G4double euX, euY, euZ; << 886 (fPosSomacomp[copyNo].x()-fshiftX) * um,
567 if (cosX > 16 * FLT_EPSILON) { << 887 (fPosSomacomp[copyNo].y()-fshiftY) * um,
568 euX = std::atan2(rotmNeuron.zy(), rotmNeur << 888 (fPosSomacomp[copyNo].z()-fshiftZ) * um
569 euY = std::atan2(-rotmNeuron.zx(), cosX); << 889 );
570 euZ = std::atan2(rotmNeuron.yx(), rotmNeur << 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());
571 } 907 }
572 else { << 908 else
573 euX = std::atan2(-rotmNeuron.yz(), rotmNeu << 909 {
574 euY = std::atan2(-rotmNeuron.zx(), cosX); << 910 euX = std::atan2 (-rotmNeuron.yz(),rotmNeuron.yy());
575 euZ = 0.; << 911 euY = std::atan2 (-rotmNeuron.zx(),cosX);
>> 912 euZ = 0. ;
576 } 913 }
577 G4RotationMatrix* rot = new G4RotationMatrix << 914 G4RotationMatrix* rot = new G4RotationMatrix();
578 rot->rotateX(euX); << 915 rot->rotateX(euX);
579 rot->rotateY(euY); << 916 rot->rotateY(euY);
580 rot->rotateZ(euZ); << 917 rot->rotateZ(euZ);
581 << 918
582 physVol->SetRotation(rot); << 919 physVol->SetRotation(rot);
583 << 920
584 // shift of cylinder compartments << 921 // shift of cylinder compartments
585 G4ThreeVector originNeuron((fPosNeuroncomp[c << 922 G4ThreeVector
586 (fPosNeuroncomp[c << 923 originNeuron(
587 (fPosNeuroncomp[c << 924 (fPosNeuroncomp[copyNo].x()-fshiftX) * um,
588 physVol->SetTranslation(originNeuron); << 925 (fPosNeuroncomp[copyNo].y()-fshiftY) * um,
589 } << 926 (fPosNeuroncomp[copyNo].z()-fshiftZ) * um
>> 927 );
>> 928 physVol->SetTranslation(originNeuron);
>> 929
>> 930 }
590 931
591 //....oooOO0OOooo........oooOO0OOooo........oo 932 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
592 << 933 /*
593 void NeuronLoadDataFile::ComputeDimensions(G4T << 934 G4VSolid* NeuronLoadDataFile::ComputeSolid(const G4int copyNo,
594 con << 935 G4VPhysicalVolume* physVol )
595 { 936 {
596 fcylinderComp.SetInnerRadius(0 * um); << 937 G4VSolid* solid;
597 fcylinderComp.SetOuterRadius(fRadNeuroncomp[ << 938 if( typeNcomp[copyNo] == 1 )
598 fcylinderComp.SetZHalfLength(fHeightNeuronco << 939 {
599 fcylinderComp.SetStartPhiAngle(0. * deg); << 940 solid = sphereComp ;
600 fcylinderComp.SetDeltaPhiAngle(360. * deg); << 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();
601 } 990 }
>> 991 */
602 992
603 //....oooOO0OOooo........oooOO0OOooo........oo 993 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
604 994