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