Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-31 09:22:17 

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // Hadrontherapy advanced example for Geant4
0027 // See more at: https://twiki.cern.ch/twiki/bin/view/Geant4/AdvancedExamplesHadrontherapy
0028 
0029 #include "HadrontherapyElectricTabulatedField3D.hh"
0030 #include "G4SystemOfUnits.hh"
0031 #include "G4AutoLock.hh"
0032 
0033 namespace{  G4Mutex MyHadrontherapyLockEField=G4MUTEX_INITIALIZER;  }
0034 
0035 using namespace std;
0036 
0037 HadrontherapyElectricTabulatedField3D::HadrontherapyElectricTabulatedField3D( const char* filename, G4double exOffset, G4double eyOffset, G4double ezOffset)
0038   :feXoffset(exOffset),feYoffset(eyOffset),feZoffset(ezOffset),einvertX(false),einvertY(false),einvertZ(false)
0039 {
0040    //The format file is: X Y Z Ex Ey Ez
0041 
0042   G4double ElenUnit= cm;
0043   G4double EfieldUnit= volt/m;
0044   G4cout << "\n-----------------------------------------------------------"
0045          << "\n      Electric field"
0046          << "\n-----------------------------------------------------------";
0047 
0048   G4cout << "\n ---> " "Reading the field grid from " << filename << " ... " << G4endl;
0049   G4AutoLock lock(&MyHadrontherapyLockEField);
0050 
0051   ifstream file( filename ); // Open the file for reading.
0052 
0053   // Ignore first blank line
0054   char ebuffer[256];
0055   file.getline(ebuffer,256);
0056 
0057   // Read table dimensions
0058   file >> Enx >> Eny >> Enz; // Note dodgy order
0059 
0060   G4cout << "  [ Number of values x,y,z: "
0061      << Enx << " " << Eny << " " << Enz << " ] "
0062      << G4endl;
0063 
0064   // Set up storage space for table
0065   xEField.resize( Enx );
0066   yEField.resize( Enx );
0067   zEField.resize( Enx );
0068   G4int ix, iy, iz;
0069   for (ix=0; ix<Enx; ix++) {
0070     xEField[ix].resize(Eny);
0071     yEField[ix].resize(Eny);
0072     zEField[ix].resize(Eny);
0073   for (iy=0; iy<Eny; iy++) {
0074       xEField[ix][iy].resize(Enz);
0075       yEField[ix][iy].resize(Enz);
0076       zEField[ix][iy].resize(Enz);
0077     }
0078   }
0079 
0080   // Read in the data
0081   G4double Exval=0.;
0082   G4double Eyval=0.;
0083   G4double Ezval=0.;
0084   G4double Ex=0.;
0085   G4double Ey=0.;
0086   G4double Ez=0.;
0087        for (iz=0; iz<Enz; iz++) {
0088          for (iy=0; iy<Eny; iy++) {
0089             for (ix=0; ix<Enx; ix++) {
0090         file >> Exval >> Eyval >> Ezval >> Ex >> Ey >> Ez;
0091 
0092         if ( ix==0 && iy==0 && iz==0 ) {
0093           Eminx = Exval * ElenUnit;
0094           Eminy = Eyval * ElenUnit;
0095           Eminz = Ezval * ElenUnit;
0096         }
0097         xEField[ix][iy][iz] = Ex * EfieldUnit;
0098         yEField[ix][iy][iz] = Ey * EfieldUnit;
0099         zEField[ix][iy][iz] = Ez * EfieldUnit;
0100       }
0101     }
0102   }
0103   file.close();
0104   lock.unlock();
0105 
0106   Emaxx = Exval * ElenUnit;
0107   Emaxy = Eyval * ElenUnit;
0108   Emaxz = Ezval * ElenUnit;
0109 
0110   G4cout << "\n ---> ... done reading " << G4endl;
0111 
0112   // G4cout << " Read values of field from file " << filename << G4endl;
0113   G4cout << " ---> assumed the order:  x, y, z, Ex, Ey, Ez "
0114      << "\n ---> Min values x,y,z: "
0115      << Eminx/cm << " " << Eminy/cm << " " << Eminz/cm << " cm "
0116      << "\n ---> Max values x,y,z: "
0117      << Emaxx/cm << " " << Emaxy/cm << " " << Emaxz/cm << " cm "
0118      << "\n ---> The field will be offset in x by " << exOffset/cm << " cm "
0119          << "\n ---> The field will be offset in y by " << eyOffset/cm << " cm "
0120          << "\n ---> The field will be offset in z by " << ezOffset/cm << " cm " << G4endl;
0121 
0122   // Should really check that the limits are not the wrong way around.
0123   if (Emaxx < Eminx) {swap(Emaxx,Eminx); einvertX = true;}
0124   if (Emaxy < Eminy) {swap(Emaxy,Eminy); einvertY = true;}
0125   if (Emaxz < Eminz) {swap(Emaxz,Eminz); einvertZ = true;}
0126   G4cout << "\nAfter reordering if neccesary"
0127      << "\n ---> Min values x,y,z: "
0128      << Eminx/cm << " " << Eminy/cm << " " << Eminz/cm << " cm "
0129      << " \n ---> Max values x,y,z: "
0130      << Emaxx/cm << " " << Emaxy/cm << " " << Emaxz/cm << " cm ";
0131 
0132   dx1 = Emaxx - Eminx;
0133   dy1 = Emaxy - Eminy;
0134   dz1 = Emaxz - Eminz;
0135   G4cout << "\n ---> Dif values x,y,z (range): "
0136      << dx1/cm << " " << dy1/cm << " " << dz1/cm << " cm  "
0137      << "\n-----------------------------------------------------------" << G4endl;
0138 }
0139 
0140 void HadrontherapyElectricTabulatedField3D::GetFieldValue(const G4double Epoint[4],
0141                       G4double *Efield ) const
0142 {
0143     G4double x1 = Epoint[0] + feXoffset;
0144     G4double y1 = Epoint[1] + feYoffset;
0145     G4double z1 = Epoint[2] + feZoffset;
0146 
0147     // Position of given point within region, normalized to the range
0148     // [0,1]
0149     G4double Exfraction = (x1 - Eminx) / dx1;
0150     G4double Eyfraction = (y1 - Eminy) / dy1;
0151     G4double Ezfraction = (z1 - Eminz) / dz1;
0152 
0153     if (einvertX) { Exfraction = 1 - Exfraction;}
0154     if (einvertY) { Eyfraction = 1 - Eyfraction;}
0155     if (einvertZ) { Ezfraction = 1 - Ezfraction;}
0156 
0157     // Need addresses of these to pass to modf below.
0158     // modf uses its second argument as an OUTPUT argument.
0159     G4double exdindex, eydindex, ezdindex;
0160 
0161     // Position of the point within the cuboid defined by the
0162     // nearest surrounding tabulated points
0163     G4double exlocal = ( std::modf(Exfraction*(Enx-1), &exdindex));
0164     G4double eylocal = ( std::modf(Eyfraction*(Eny-1), &eydindex));
0165     G4double ezlocal = ( std::modf(Ezfraction*(Enz-1), &ezdindex));
0166 
0167     // The indices of the nearest tabulated point whose coordinates
0168     // are all less than those of the given point
0169     G4int exindex = static_cast<G4int>(std::floor(exdindex));
0170     G4int eyindex = static_cast<G4int>(std::floor(eydindex));
0171     G4int ezindex = static_cast<G4int>(std::floor(ezdindex));
0172 
0173     if ((exindex < 0) || (exindex >= Enx - 1) ||
0174         (eyindex < 0) || (eyindex >= Eny - 1) ||
0175         (ezindex < 0) || (ezindex >= Enz - 1))
0176     {
0177         Efield[0] = 0.0;
0178         Efield[1] = 0.0;
0179         Efield[2] = 0.0;
0180         Efield[3] = 0.0;
0181         Efield[4] = 0.0;
0182         Efield[5] = 0.0;
0183     }
0184     else
0185     {
0186 
0187 /*
0188 #ifdef DEBUG_G4intERPOLATING_FIELD
0189     G4cout << "Local x,y,z: " << exlocal << " " << eylocal << " " << ezlocal << G4endl;
0190     G4cout << "Index x,y,z: " << exindex << " " << eyindex << " " << ezindex << G4endl;
0191     G4double valx0z0, mulx0z0, valx1z0, mulx1z0;
0192     G4double valx0z1, mulx0z1, valx1z1, mulx1z1;
0193     valx0z0= table[exindex  ][0][ezindex];  mulx0z0=  (1-exlocal) * (1-ezlocal);
0194     valx1z0= table[exindex+1][0][ezindex];  mulx1z0=   exlocal    * (1-ezlocal);
0195     valx0z1= table[exindex  ][0][ezindex+1]; mulx0z1= (1-exlocal) * ezlocal;
0196     valx1z1= table[exindex+1][0][ezindex+1]; mulx1z1=  exlocal    * ezlocal;
0197 #endif
0198 */
0199         // Full 3-dimensional version
0200 
0201         Efield[0] = 0.0;
0202         Efield[1] = 0.0;
0203         Efield[2] = 0.0;
0204 
0205         Efield[3] =
0206           xEField[exindex  ][eyindex  ][ezindex  ] * (1-exlocal) * (1-eylocal) * (1-ezlocal) +
0207           xEField[exindex  ][eyindex  ][ezindex+1] * (1-exlocal) * (1-eylocal) *    ezlocal  +
0208           xEField[exindex  ][eyindex+1][ezindex  ] * (1-exlocal) *    eylocal  * (1-ezlocal) +
0209           xEField[exindex  ][eyindex+1][ezindex+1] * (1-exlocal) *    eylocal  *    ezlocal  +
0210           xEField[exindex+1][eyindex  ][ezindex  ] *    exlocal  * (1-eylocal) * (1-ezlocal) +
0211           xEField[exindex+1][eyindex  ][ezindex+1] *    exlocal  * (1-eylocal) *    ezlocal  +
0212           xEField[exindex+1][eyindex+1][ezindex  ] *    exlocal  *    eylocal  * (1-ezlocal) +
0213           xEField[exindex+1][eyindex+1][ezindex+1] *    exlocal  *    eylocal  *    ezlocal ;
0214         Efield[4] =
0215           yEField[exindex  ][eyindex  ][ezindex  ] * (1-exlocal) * (1-eylocal) * (1-ezlocal) +
0216           yEField[exindex  ][eyindex  ][ezindex+1] * (1-exlocal) * (1-eylocal) *    ezlocal  +
0217           yEField[exindex  ][eyindex+1][ezindex  ] * (1-exlocal) *    eylocal  * (1-ezlocal) +
0218           yEField[exindex  ][eyindex+1][ezindex+1] * (1-exlocal) *    eylocal  *    ezlocal  +
0219           yEField[exindex+1][eyindex  ][ezindex  ] *    exlocal  * (1-eylocal) * (1-ezlocal) +
0220           yEField[exindex+1][eyindex  ][ezindex+1] *    exlocal  * (1-eylocal) *    ezlocal  +
0221           yEField[exindex+1][eyindex+1][ezindex  ] *    exlocal  *    eylocal  * (1-ezlocal) +
0222           yEField[exindex+1][eyindex+1][ezindex+1] *    exlocal  *    eylocal  *    ezlocal ;
0223         Efield[5] =
0224           zEField[exindex  ][eyindex  ][ezindex  ] * (1-exlocal) * (1-eylocal) * (1-ezlocal) +
0225           zEField[exindex  ][eyindex  ][ezindex+1] * (1-exlocal) * (1-eylocal) *    ezlocal  +
0226           zEField[exindex  ][eyindex+1][ezindex  ] * (1-exlocal) *    eylocal  * (1-ezlocal) +
0227           zEField[exindex  ][eyindex+1][ezindex+1] * (1-exlocal) *    eylocal  *    ezlocal  +
0228           zEField[exindex+1][eyindex  ][ezindex  ] *    exlocal  * (1-eylocal) * (1-ezlocal) +
0229           zEField[exindex+1][eyindex  ][ezindex+1] *    exlocal  * (1-eylocal) *    ezlocal  +
0230           zEField[exindex+1][eyindex+1][ezindex  ] *    exlocal  *    eylocal  * (1-ezlocal) +
0231           zEField[exindex+1][eyindex+1][ezindex+1] *    exlocal  *    eylocal  *    ezlocal ;
0232   }
0233 //G4cout << "Getting electric field " << Efield[3]/(volt/m) << " " << Efield[4]/(volt/m) << " " << Efield[5]/(volt/m) << G4endl;
0234 //G4cout << "For coordinates: " << Epoint[0] << " " << Epoint[1] << " " << Epoint[2] << G4endl;
0235 
0236 /*std::ofstream WriteDataIn("ElectricFieldFC.out", std::ios::app);
0237        WriteDataIn  <<   Epoint[0]             << '\t' << "   "
0238             <<   Epoint[1]           << '\t' << "   "
0239             <<   Epoint[2]            << '\t' << "   "
0240             <<   Efield[3]/(volt/m)            << '\t' << "   "
0241             <<   Efield[4]/(volt/m)           << '\t' << "   "
0242             <<   Efield[5]/(volt/m)           << '\t' << "   "
0243             << std::endl;    */
0244 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!