Geant4 Cross Reference

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

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

Diff markup

Differences between /examples/extended/medical/dna/neuron/src/NeuronLoadDataFile.cc (Version 11.3.0) and /examples/extended/medical/dna/neuron/src/NeuronLoadDataFile.cc (Version 10.7)


  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