Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/hadrontherapy/src/HadrontherapyMagneticField3D.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/advanced/hadrontherapy/src/HadrontherapyMagneticField3D.cc (Version 11.3.0) and /examples/advanced/hadrontherapy/src/HadrontherapyMagneticField3D.cc (Version 10.5.p1)


  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 // Hadrontherapy advanced example for Geant4       26 // Hadrontherapy advanced example for Geant4
 27 // See more at: https://twiki.cern.ch/twiki/bi     27 // See more at: https://twiki.cern.ch/twiki/bin/view/Geant4/AdvancedExamplesHadrontherapy
 28                                                    28 
 29 #include "HadrontherapyMagneticField3D.hh"         29 #include "HadrontherapyMagneticField3D.hh"
 30 #include "G4SystemOfUnits.hh"                      30 #include "G4SystemOfUnits.hh"
 31 #include "G4AutoLock.hh"                           31 #include "G4AutoLock.hh"
 32                                                    32 
 33 namespace{  G4Mutex MyHadrontherapyLock=G4MUTE     33 namespace{  G4Mutex MyHadrontherapyLock=G4MUTEX_INITIALIZER;  }
 34                                                    34 
 35 using namespace std;                           << 
 36                                                << 
 37 HadrontherapyMagneticField3D::HadrontherapyMag     35 HadrontherapyMagneticField3D::HadrontherapyMagneticField3D( const char* filename, double xOffset )
 38   :fXoffset(xOffset),invertX(false),invertY(fa     36   :fXoffset(xOffset),invertX(false),invertY(false),invertZ(false)
 39 {                                                  37 {
 40    //The format file is: X Y Z Ex Ey Ez            38    //The format file is: X Y Z Ex Ey Ez
 41                                                    39 
 42   double lenUnit= meter;                           40   double lenUnit= meter;
 43   double fieldUnit= tesla;                         41   double fieldUnit= tesla;
 44   G4cout << "\n-------------------------------     42   G4cout << "\n-----------------------------------------------------------"
 45    << "\n      Magnetic field"                     43    << "\n      Magnetic field"
 46    << "\n-------------------------------------     44    << "\n-----------------------------------------------------------";
 47                                                    45 
 48                                                    46 
 49   G4cout << "\n ---> " "Reading the field grid <<  47   G4cout << "\n ---> " "Reading the field grid from " << filename << " ... " << endl;
 50   G4AutoLock lock(&MyHadrontherapyLock);           48   G4AutoLock lock(&MyHadrontherapyLock);
 51                                                    49 
 52   ifstream file( filename ); // Open the file      50   ifstream file( filename ); // Open the file for reading.
 53                                                    51 
 54   // Ignore first blank line                       52   // Ignore first blank line
 55   char buffer[256];                                53   char buffer[256];
 56   file.getline(buffer,256);                        54   file.getline(buffer,256);
 57                                                    55 
 58   // Read table dimensions                         56   // Read table dimensions
 59   file >> nx >> ny >> nz; // Note dodgy order      57   file >> nx >> ny >> nz; // Note dodgy order
 60                                                    58 
 61   G4cout << "  [ Number of values x,y,z: "         59   G4cout << "  [ Number of values x,y,z: "
 62    << nx << " " << ny << " " << nz << " ] "        60    << nx << " " << ny << " " << nz << " ] "
 63    << G4endl;                                  <<  61    << endl;
 64                                                    62 
 65   // Set up storage space for table                63   // Set up storage space for table
 66   xField.resize( nx );                             64   xField.resize( nx );
 67   yField.resize( nx );                             65   yField.resize( nx );
 68   zField.resize( nx );                             66   zField.resize( nx );
 69   int ix, iy, iz;                                  67   int ix, iy, iz;
 70   for (ix=0; ix<nx; ix++) {                        68   for (ix=0; ix<nx; ix++) {
 71     xField[ix].resize(ny);                         69     xField[ix].resize(ny);
 72     yField[ix].resize(ny);                         70     yField[ix].resize(ny);
 73     zField[ix].resize(ny);                         71     zField[ix].resize(ny);
 74     for (iy=0; iy<ny; iy++) {                      72     for (iy=0; iy<ny; iy++) {
 75       xField[ix][iy].resize(nz);                   73       xField[ix][iy].resize(nz);
 76       yField[ix][iy].resize(nz);                   74       yField[ix][iy].resize(nz);
 77       zField[ix][iy].resize(nz);                   75       zField[ix][iy].resize(nz);
 78     }                                              76     }
 79   }                                                77   }
 80                                                    78 
 81   // Read in the data                              79   // Read in the data
 82   G4double xval=0.;                                80   G4double xval=0.;
 83   G4double yval=0.;                                81   G4double yval=0.;
 84   G4double zval=0.;                                82   G4double zval=0.;
 85   G4double bx=0.;                                  83   G4double bx=0.;
 86   G4double by=0.;                                  84   G4double by=0.;
 87   G4double bz=0.;                                  85   G4double bz=0.;
 88   for (ix=0; ix<nx; ix++) {                        86   for (ix=0; ix<nx; ix++) {
 89     for (iy=0; iy<ny; iy++) {                      87     for (iy=0; iy<ny; iy++) {
 90       for (iz=0; iz<nz; iz++) {                    88       for (iz=0; iz<nz; iz++) {
 91         file >> xval >> yval >> zval >> bx >>      89         file >> xval >> yval >> zval >> bx >> by >> bz ;
 92         if ( ix==0 && iy==0 && iz==0 ) {           90         if ( ix==0 && iy==0 && iz==0 ) {
 93           minx = xval * lenUnit;                   91           minx = xval * lenUnit;
 94           miny = yval * lenUnit;                   92           miny = yval * lenUnit;
 95           minz = zval * lenUnit;                   93           minz = zval * lenUnit;
 96         }                                          94         }
 97         xField[ix][iy][iz] = bx * fieldUnit;       95         xField[ix][iy][iz] = bx * fieldUnit;
 98         yField[ix][iy][iz] = by * fieldUnit;       96         yField[ix][iy][iz] = by * fieldUnit;
 99         zField[ix][iy][iz] = bz * fieldUnit;       97         zField[ix][iy][iz] = bz * fieldUnit;
100       }                                            98       }
101     }                                              99     }
102   }                                               100   }
103   file.close();                                   101   file.close();
104                                                   102 
105   lock.unlock();                                  103   lock.unlock();
106                                                   104 
107   maxx = xval * lenUnit;                          105   maxx = xval * lenUnit;
108   maxy = yval * lenUnit;                          106   maxy = yval * lenUnit;
109   maxz = zval * lenUnit;                          107   maxz = zval * lenUnit;
110                                                   108 
111   G4cout << "\n ---> ... done reading " << G4e << 109   G4cout << "\n ---> ... done reading " << endl;
112                                                   110 
113   // G4cout << " Read values of field from fil << 111   // G4cout << " Read values of field from file " << filename << endl;
114   G4cout << " ---> assumed the order:  x, y, z    112   G4cout << " ---> assumed the order:  x, y, z, Bx, By, Bz "
115    << "\n ---> Min values x,y,z: "                113    << "\n ---> Min values x,y,z: "
116    << minx/cm << " " << miny/cm << " " << minz    114    << minx/cm << " " << miny/cm << " " << minz/cm << " cm "
117    << "\n ---> Max values x,y,z: "                115    << "\n ---> Max values x,y,z: "
118    << maxx/cm << " " << maxy/cm << " " << maxz    116    << maxx/cm << " " << maxy/cm << " " << maxz/cm << " cm "
119    << "\n ---> The field will be offset by " < << 117    << "\n ---> The field will be offset by " << xOffset/cm << " cm " << endl;
120                                                   118 
121   // Should really check that the limits are n    119   // Should really check that the limits are not the wrong way around.
122   if (maxx < minx) {swap(maxx,minx); invertX =    120   if (maxx < minx) {swap(maxx,minx); invertX = true;}
123   if (maxy < miny) {swap(maxy,miny); invertY =    121   if (maxy < miny) {swap(maxy,miny); invertY = true;}
124   if (maxz < minz) {swap(maxz,minz); invertZ =    122   if (maxz < minz) {swap(maxz,minz); invertZ = true;}
125   G4cout << "\nAfter reordering if neccesary"     123   G4cout << "\nAfter reordering if neccesary"
126    << "\n ---> Min values x,y,z: "                124    << "\n ---> Min values x,y,z: "
127    << minx/cm << " " << miny/cm << " " << minz    125    << minx/cm << " " << miny/cm << " " << minz/cm << " cm "
128    << " \n ---> Max values x,y,z: "               126    << " \n ---> Max values x,y,z: "
129    << maxx/cm << " " << maxy/cm << " " << maxz    127    << maxx/cm << " " << maxy/cm << " " << maxz/cm << " cm ";
130                                                   128 
131   dx = maxx - minx;                               129   dx = maxx - minx;
132   dy = maxy - miny;                               130   dy = maxy - miny;
133   dz = maxz - minz;                               131   dz = maxz - minz;
134   G4cout << "\n ---> Dif values x,y,z (range):    132   G4cout << "\n ---> Dif values x,y,z (range): "
135    << dx/cm << " " << dy/cm << " " << dz/cm <<    133    << dx/cm << " " << dy/cm << " " << dz/cm << " cm in z "
136    << "\n------------------------------------- << 134    << "\n-----------------------------------------------------------" << endl;
137 }                                                 135 }
138                                                   136 
139 void HadrontherapyMagneticField3D::GetFieldVal    137 void HadrontherapyMagneticField3D::GetFieldValue(const double point[4],
140               double *Bfield ) const              138               double *Bfield ) const
141 {                                                 139 {
142     double x = point[0]+ fXoffset;                140     double x = point[0]+ fXoffset;
143     double y = point[1];                          141     double y = point[1];
144     double z = point[2];                          142     double z = point[2];
145                                                   143 
146     // Position of given point within region,     144     // Position of given point within region, normalized to the range
147     // [0,1]                                      145     // [0,1]
148     double xfraction = (x - minx) / dx;           146     double xfraction = (x - minx) / dx;
149     double yfraction = (y - miny) / dy;           147     double yfraction = (y - miny) / dy;
150     double zfraction = (z - minz) / dz;           148     double zfraction = (z - minz) / dz;
151                                                   149 
152     if (invertX) { xfraction = 1 - xfraction;}    150     if (invertX) { xfraction = 1 - xfraction;}
153     if (invertY) { yfraction = 1 - yfraction;}    151     if (invertY) { yfraction = 1 - yfraction;}
154     if (invertZ) { zfraction = 1 - zfraction;}    152     if (invertZ) { zfraction = 1 - zfraction;}
155                                                   153 
156     // Need addresses of these to pass to modf    154     // Need addresses of these to pass to modf below.
157     // modf uses its second argument as an OUT    155     // modf uses its second argument as an OUTPUT argument.
158     double xdindex, ydindex, zdindex;             156     double xdindex, ydindex, zdindex;
159                                                   157 
160     // Position of the point within the cuboid    158     // Position of the point within the cuboid defined by the
161     // nearest surrounding tabulated points       159     // nearest surrounding tabulated points
162     double xlocal = ( std::modf(xfraction*(nx-    160     double xlocal = ( std::modf(xfraction*(nx-1), &xdindex));
163     double ylocal = ( std::modf(yfraction*(ny-    161     double ylocal = ( std::modf(yfraction*(ny-1), &ydindex));
164     double zlocal = ( std::modf(zfraction*(nz-    162     double zlocal = ( std::modf(zfraction*(nz-1), &zdindex));
165                                                   163 
166     // The indices of the nearest tabulated po    164     // The indices of the nearest tabulated point whose coordinates
167     // are all less than those of the given po    165     // are all less than those of the given point
168     int xindex = static_cast<int>(std::floor(x    166     int xindex = static_cast<int>(std::floor(xdindex));
169     int yindex = static_cast<int>(std::floor(y    167     int yindex = static_cast<int>(std::floor(ydindex));
170     int zindex = static_cast<int>(std::floor(z    168     int zindex = static_cast<int>(std::floor(zdindex));
171                                                   169 
172       // Check that the point is within the de    170       // Check that the point is within the defined region
173     if ((xindex < 0) || (xindex >= nx - 1) ||     171     if ((xindex < 0) || (xindex >= nx - 1) ||
174         (yindex < 0) || (yindex >= ny - 1) ||     172         (yindex < 0) || (yindex >= ny - 1) ||
175         (zindex < 0) || (zindex >= nz - 1))       173         (zindex < 0) || (zindex >= nz - 1))
176     {                                             174     {
177         Bfield[0] = 0.0;                          175         Bfield[0] = 0.0;
178         Bfield[1] = 0.0;                          176         Bfield[1] = 0.0;
179         Bfield[2] = 0.0;                          177         Bfield[2] = 0.0;
180     }                                             178     }
181     else                                          179     else
182     {                                             180     {
183                                                   181 
184 #ifdef DEBUG_INTERPOLATING_FIELD                  182 #ifdef DEBUG_INTERPOLATING_FIELD
185         G4cout << "Local x,y,z: " << xlocal << << 183         G4cout << "Local x,y,z: " << xlocal << " " << ylocal << " " << zlocal << endl;
186         G4cout << "Index x,y,z: " << xindex << << 184         G4cout << "Index x,y,z: " << xindex << " " << yindex << " " << zindex << endl;
187         double valx0z0, mulx0z0, valx1z0, mulx    185         double valx0z0, mulx0z0, valx1z0, mulx1z0;
188         double valx0z1, mulx0z1, valx1z1, mulx    186         double valx0z1, mulx0z1, valx1z1, mulx1z1;
189         valx0z0= table[xindex  ][0][zindex];      187         valx0z0= table[xindex  ][0][zindex];  mulx0z0=  (1-xlocal) * (1-zlocal);
190         valx1z0= table[xindex+1][0][zindex];      188         valx1z0= table[xindex+1][0][zindex];  mulx1z0=   xlocal    * (1-zlocal);
191         valx0z1= table[xindex  ][0][zindex+1];    189         valx0z1= table[xindex  ][0][zindex+1]; mulx0z1= (1-xlocal) * zlocal;
192         valx1z1= table[xindex+1][0][zindex+1];    190         valx1z1= table[xindex+1][0][zindex+1]; mulx1z1=  xlocal    * zlocal;
193 #endif                                            191 #endif
194                                                   192 
195         // Full 3-dimensional version             193         // Full 3-dimensional version
196         Bfield[0] =                               194         Bfield[0] =
197           xField[xindex  ][yindex  ][zindex  ]    195           xField[xindex  ][yindex  ][zindex  ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
198           xField[xindex  ][yindex  ][zindex+1]    196           xField[xindex  ][yindex  ][zindex+1] * (1-xlocal) * (1-ylocal) *    zlocal  +
199           xField[xindex  ][yindex+1][zindex  ]    197           xField[xindex  ][yindex+1][zindex  ] * (1-xlocal) *    ylocal  * (1-zlocal) +
200           xField[xindex  ][yindex+1][zindex+1]    198           xField[xindex  ][yindex+1][zindex+1] * (1-xlocal) *    ylocal  *    zlocal  +
201           xField[xindex+1][yindex  ][zindex  ]    199           xField[xindex+1][yindex  ][zindex  ] *    xlocal  * (1-ylocal) * (1-zlocal) +
202           xField[xindex+1][yindex  ][zindex+1]    200           xField[xindex+1][yindex  ][zindex+1] *    xlocal  * (1-ylocal) *    zlocal  +
203           xField[xindex+1][yindex+1][zindex  ]    201           xField[xindex+1][yindex+1][zindex  ] *    xlocal  *    ylocal  * (1-zlocal) +
204           xField[xindex+1][yindex+1][zindex+1]    202           xField[xindex+1][yindex+1][zindex+1] *    xlocal  *    ylocal  *    zlocal ;
205                                                   203 
206         Bfield[1] =                               204         Bfield[1] =
207           yField[xindex  ][yindex  ][zindex  ]    205           yField[xindex  ][yindex  ][zindex  ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
208           yField[xindex  ][yindex  ][zindex+1]    206           yField[xindex  ][yindex  ][zindex+1] * (1-xlocal) * (1-ylocal) *    zlocal  +
209           yField[xindex  ][yindex+1][zindex  ]    207           yField[xindex  ][yindex+1][zindex  ] * (1-xlocal) *    ylocal  * (1-zlocal) +
210           yField[xindex  ][yindex+1][zindex+1]    208           yField[xindex  ][yindex+1][zindex+1] * (1-xlocal) *    ylocal  *    zlocal  +
211           yField[xindex+1][yindex  ][zindex  ]    209           yField[xindex+1][yindex  ][zindex  ] *    xlocal  * (1-ylocal) * (1-zlocal) +
212           yField[xindex+1][yindex  ][zindex+1]    210           yField[xindex+1][yindex  ][zindex+1] *    xlocal  * (1-ylocal) *    zlocal  +
213           yField[xindex+1][yindex+1][zindex  ]    211           yField[xindex+1][yindex+1][zindex  ] *    xlocal  *    ylocal  * (1-zlocal) +
214           yField[xindex+1][yindex+1][zindex+1]    212           yField[xindex+1][yindex+1][zindex+1] *    xlocal  *    ylocal  *    zlocal ;
215                                                   213 
216         Bfield[2] =                               214         Bfield[2] =
217           zField[xindex  ][yindex  ][zindex  ]    215           zField[xindex  ][yindex  ][zindex  ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
218           zField[xindex  ][yindex  ][zindex+1]    216           zField[xindex  ][yindex  ][zindex+1] * (1-xlocal) * (1-ylocal) *    zlocal  +
219           zField[xindex  ][yindex+1][zindex  ]    217           zField[xindex  ][yindex+1][zindex  ] * (1-xlocal) *    ylocal  * (1-zlocal) +
220           zField[xindex  ][yindex+1][zindex+1]    218           zField[xindex  ][yindex+1][zindex+1] * (1-xlocal) *    ylocal  *    zlocal  +
221           zField[xindex+1][yindex  ][zindex  ]    219           zField[xindex+1][yindex  ][zindex  ] *    xlocal  * (1-ylocal) * (1-zlocal) +
222           zField[xindex+1][yindex  ][zindex+1]    220           zField[xindex+1][yindex  ][zindex+1] *    xlocal  * (1-ylocal) *    zlocal  +
223           zField[xindex+1][yindex+1][zindex  ]    221           zField[xindex+1][yindex+1][zindex  ] *    xlocal  *    ylocal  * (1-zlocal) +
224           zField[xindex+1][yindex+1][zindex+1]    222           zField[xindex+1][yindex+1][zindex+1] *    xlocal  *    ylocal  *    zlocal ;
225     }                                             223     }
226                                                   224 
227 //In order to obtain the output file with the     225 //In order to obtain the output file with the magnetic components read from a particle passing in the magnetic field
228 /*  std::ofstream MagneticField("MagneticField    226 /*  std::ofstream MagneticField("MagneticField.out", std::ios::app);
229      MagneticField<<   Bfield[0] << '\t' << "     227      MagneticField<<   Bfield[0] << '\t' << "   "
230       <<   Bfield[1] << '\t' << "    "            228       <<   Bfield[1] << '\t' << "    "
231       <<   Bfield[2] << '\t' << "   "             229       <<   Bfield[2] << '\t' << "   "
232       <<   point[0] << '\t' << "   "              230       <<   point[0] << '\t' << "   "
233       <<   point[1] << '\t' << "    "             231       <<   point[1] << '\t' << "    "
234       <<   point[2] << '\t' << "   "              232       <<   point[2] << '\t' << "   "
235       << std::endl;*/                          << 233       << G4endl;*/
236                                                   234 
237 }                                                 235 }
238                                                   236