Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/nanobeam/src/TabulatedField3D.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/nanobeam/src/TabulatedField3D.cc (Version 11.3.0) and /examples/advanced/nanobeam/src/TabulatedField3D.cc (Version 10.3)


  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 // Please cite the following papers if you use     26 // Please cite the following papers if you use this software
 27 // Nucl.Instrum.Meth.B260:20-27, 2007              27 // Nucl.Instrum.Meth.B260:20-27, 2007
 28 // IEEE TNS 51, 4:1395-1401, 2004                  28 // IEEE TNS 51, 4:1395-1401, 2004
 29 //                                             << 
 30 // Based on purging magnet advanced example    << 
 31 //                                             << 
 32                                                    29 
 33 #include "TabulatedField3D.hh"                     30 #include "TabulatedField3D.hh"
 34 #include "G4SystemOfUnits.hh"                      31 #include "G4SystemOfUnits.hh"
 35 #include "G4Exp.hh"                                32 #include "G4Exp.hh"
 36                                                    33 
 37 #include "G4AutoLock.hh"                       << 
 38                                                << 
 39 namespace                                      << 
 40 {                                              << 
 41   G4Mutex myTabulatedField3DLock = G4MUTEX_INI << 
 42 }                                              << 
 43                                                << 
 44 //....oooOO0OOooo........oooOO0OOooo........oo     34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 45                                                    35 
 46 TabulatedField3D::TabulatedField3D(G4float gr1     36 TabulatedField3D::TabulatedField3D(G4float gr1, G4float gr2, G4float gr3, G4float gr4, G4int choiceModel) 
 47 {                                                  37 {    
 48                                                    38 
 49   G4cout << " ********************** " << G4en     39   G4cout << " ********************** " << G4endl;
 50   G4cout << " **** CONFIGURATION *** " << G4en     40   G4cout << " **** CONFIGURATION *** " << G4endl;
 51   G4cout << " ********************** " << G4en     41   G4cout << " ********************** " << G4endl;
 52                                                    42 
 53   G4cout<< G4endl;                                 43   G4cout<< G4endl; 
 54   G4cout << "=====> You have selected :" << G4     44   G4cout << "=====> You have selected :" << G4endl;
 55   if (choiceModel==1) G4cout<< "-> Square quad     45   if (choiceModel==1) G4cout<< "-> Square quadrupole field"<<G4endl;
 56   if (choiceModel==2) G4cout<< "-> 3D quadrupo     46   if (choiceModel==2) G4cout<< "-> 3D quadrupole field"<<G4endl;
 57   if (choiceModel==3) G4cout<< "-> Enge quadru     47   if (choiceModel==3) G4cout<< "-> Enge quadrupole field"<<G4endl;
 58   G4cout << "   G1 (T/m) = "<< gr1 << G4endl;      48   G4cout << "   G1 (T/m) = "<< gr1 << G4endl;
 59   G4cout << "   G2 (T/m) = "<< gr2 << G4endl;      49   G4cout << "   G2 (T/m) = "<< gr2 << G4endl;
 60   G4cout << "   G3 (T/m) = "<< gr3 << G4endl;      50   G4cout << "   G3 (T/m) = "<< gr3 << G4endl;
 61   G4cout << "   G4 (T/m) = "<< gr4 << G4endl;      51   G4cout << "   G4 (T/m) = "<< gr4 << G4endl;
 62                                                    52   
 63   fGradient1 = gr1;                                53   fGradient1 = gr1;
 64   fGradient2 = gr2;                                54   fGradient2 = gr2;
 65   fGradient3 = gr3;                                55   fGradient3 = gr3;
 66   fGradient4 = gr4;                                56   fGradient4 = gr4;
 67   fModel = choiceModel;                            57   fModel = choiceModel;
 68                                                    58   
 69   if (fModel==2)                                   59   if (fModel==2)
 70   {                                                60   {
 71   //                                           << 
 72   //This is a thread-local class and we have t << 
 73   //file at the same time                      << 
 74   G4AutoLock lock(&myTabulatedField3DLock);    << 
 75   //                                           << 
 76                                                << 
 77   const char * filename ="OM50.grid";              61   const char * filename ="OM50.grid";
 78                                                    62   
 79   double lenUnit= mm;                              63   double lenUnit= mm;
 80   G4cout << "\n-------------------------------     64   G4cout << "\n-----------------------------------------------------------"
 81    << "\n      3D Magnetic field from OPERA so     65    << "\n      3D Magnetic field from OPERA software "
 82    << "\n-------------------------------------     66    << "\n-----------------------------------------------------------";
 83                                                    67     
 84   G4cout << "\n ---> " "Reading the field grid <<  68   G4cout << "\n ---> " "Reading the field grid from " << filename << " ... " << endl; 
 85   std::ifstream file( filename ); // Open the  <<  69   ifstream file( filename ); // Open the file for reading.
 86                                                    70   
 87   // Read table dimensions                         71   // Read table dimensions 
 88   file >> fNx >> fNy >> fNz; // Note dodgy ord     72   file >> fNx >> fNy >> fNz; // Note dodgy order
 89                                                    73 
 90   G4cout << "  [ Number of values x,y,z: "         74   G4cout << "  [ Number of values x,y,z: " 
 91    << fNx << " " << fNy << " " << fNz << " ] "     75    << fNx << " " << fNy << " " << fNz << " ] "
 92    << G4endl;                                  <<  76    << endl;
 93                                                    77 
 94   // Set up storage space for table                78   // Set up storage space for table
 95   fXField.resize( fNx );                           79   fXField.resize( fNx );
 96   fYField.resize( fNx );                           80   fYField.resize( fNx );
 97   fZField.resize( fNx );                           81   fZField.resize( fNx );
 98   int ix, iy, iz;                                  82   int ix, iy, iz;
 99   for (ix=0; ix<fNx; ix++)                     <<  83   for (ix=0; ix<fNx; ix++) {
100   {                                            << 
101     fXField[ix].resize(fNy);                       84     fXField[ix].resize(fNy);
102     fYField[ix].resize(fNy);                       85     fYField[ix].resize(fNy);
103     fZField[ix].resize(fNy);                       86     fZField[ix].resize(fNy);
104     for (iy=0; iy<fNy; iy++)                   <<  87     for (iy=0; iy<fNy; iy++) {
105     {                                          << 
106       fXField[ix][iy].resize(fNz);                 88       fXField[ix][iy].resize(fNz);
107       fYField[ix][iy].resize(fNz);                 89       fYField[ix][iy].resize(fNz);
108       fZField[ix][iy].resize(fNz);                 90       fZField[ix][iy].resize(fNz);
109     }                                              91     }
110   }                                                92   }
111                                                    93   
112   // Read in the data                              94   // Read in the data
113   double xval,yval,zval,bx,by,bz;                  95   double xval,yval,zval,bx,by,bz;
114   double permeability; // Not used in this exa     96   double permeability; // Not used in this example.
115   for (ix=0; ix<fNx; ix++)                     <<  97   for (ix=0; ix<fNx; ix++) {
116   {                                            <<  98     for (iy=0; iy<fNy; iy++) {
117     for (iy=0; iy<fNy; iy++)                   <<  99       for (iz=0; iz<fNz; iz++) {
118     {                                          << 
119       for (iz=0; iz<fNz; iz++)                 << 
120       {                                        << 
121         file >> xval >> yval >> zval >> bx >>     100         file >> xval >> yval >> zval >> bx >> by >> bz >> permeability;
122         if ( ix==0 && iy==0 && iz==0 )         << 101         if ( ix==0 && iy==0 && iz==0 ) {
123   {                                            << 
124           fMinix = xval * lenUnit;                102           fMinix = xval * lenUnit;
125           fMiniy = yval * lenUnit;                103           fMiniy = yval * lenUnit;
126           fMiniz = zval * lenUnit;                104           fMiniz = zval * lenUnit;
127         }                                         105         }
128         fXField[ix][iy][iz] = bx ;                106         fXField[ix][iy][iz] = bx ;
129         fYField[ix][iy][iz] = by ;                107         fYField[ix][iy][iz] = by ;
130         fZField[ix][iy][iz] = bz ;                108         fZField[ix][iy][iz] = bz ;
131       }                                           109       }
132     }                                             110     }
133   }                                               111   }
134   file.close();                                   112   file.close();
135                                                   113 
136   //                                           << 
137   lock.unlock();                               << 
138   //                                           << 
139                                                << 
140   fMaxix = xval * lenUnit;                        114   fMaxix = xval * lenUnit;
141   fMaxiy = yval * lenUnit;                        115   fMaxiy = yval * lenUnit;
142   fMaxiz = zval * lenUnit;                        116   fMaxiz = zval * lenUnit;
143                                                   117 
144   G4cout << "\n ---> ... done reading " << G4e << 118   G4cout << "\n ---> ... done reading " << endl;
145                                                << 
146   // G4cout << " Read values of field from fil << 
147                                                   119 
                                                   >> 120   // G4cout << " Read values of field from file " << filename << endl; 
148   G4cout << " ---> assumed the order:  x, y, z    121   G4cout << " ---> assumed the order:  x, y, z, Bx, By, Bz "
149    << "\n ---> Min values x,y,z: "                122    << "\n ---> Min values x,y,z: " 
150    << fMinix/cm << " " << fMiniy/cm << " " <<     123    << fMinix/cm << " " << fMiniy/cm << " " << fMiniz/cm << " cm "
151    << "\n ---> Max values x,y,z: "                124    << "\n ---> Max values x,y,z: " 
152    << fMaxix/cm << " " << fMaxiy/cm << " " <<  << 125    << fMaxix/cm << " " << fMaxiy/cm << " " << fMaxiz/cm << " cm " << endl;
153                                                   126 
154   fDx = fMaxix - fMinix;                          127   fDx = fMaxix - fMinix;
155   fDy = fMaxiy - fMiniy;                          128   fDy = fMaxiy - fMiniy;
156   fDz = fMaxiz - fMiniz;                          129   fDz = fMaxiz - fMiniz;
157                                                << 
158   G4cout << "\n ---> Dif values x,y,z (range):    130   G4cout << "\n ---> Dif values x,y,z (range): " 
159    << fDx/cm << " " << fDy/cm << " " << fDz/cm    131    << fDx/cm << " " << fDy/cm << " " << fDz/cm << " cm in z "
160    << "\n------------------------------------- << 132    << "\n-----------------------------------------------------------" << endl;
161                                                   133 
162                                                   134   
163   // Table normalization                          135   // Table normalization
164                                                << 
165   for (ix=0; ix<fNx; ix++)                        136   for (ix=0; ix<fNx; ix++) 
166   {                                               137   {
167     for (iy=0; iy<fNy; iy++)                      138     for (iy=0; iy<fNy; iy++) 
168     {                                             139     {
169       for (iz=0; iz<fNz; iz++)                    140       for (iz=0; iz<fNz; iz++) 
170       {                                           141       {
171                                                   142 
172   fXField[ix][iy][iz] = (fXField[ix][iy][iz]/1    143   fXField[ix][iy][iz] = (fXField[ix][iy][iz]/197.736);
173         fYField[ix][iy][iz] = (fYField[ix][iy]    144         fYField[ix][iy][iz] = (fYField[ix][iy][iz]/197.736);
174         fZField[ix][iy][iz] = (fZField[ix][iy]    145         fZField[ix][iy][iz] = (fZField[ix][iy][iz]/197.736);
175                                                   146 
176       }                                        << 147   }
177     }                                             148     }
178   }                                               149   }
179                                                   150 
180   } // fModel==2                                  151   } // fModel==2
181                                                   152 
182 }                                                 153 }
183                                                   154  
184 //....oooOO0OOooo........oooOO0OOooo........oo    155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
185                                                   156 
186                                                   157 
187 void TabulatedField3D::GetFieldValue(const dou    158 void TabulatedField3D::GetFieldValue(const double point[4],
188               double *Bfield ) const              159               double *Bfield ) const
189 {                                                 160 { 
190   //G4cout << fGradient1 << G4endl;            << 
191   //G4cout << fGradient2 << G4endl;            << 
192   //G4cout << fGradient3 << G4endl;            << 
193   //G4cout << fGradient4 << G4endl;            << 
194   //G4cout << "---------" << G4endl;           << 
195                                                   161 
196   G4double coef, G0;                              162   G4double coef, G0;
197   G0 = 0;                                         163   G0 = 0;
198                                                   164   
199   coef=1; // for protons                          165   coef=1; // for protons
200   //coef=2; // for alphas                         166   //coef=2; // for alphas
201                                                   167 
202 //********************************************    168 //******************************************************************
203                                                   169 
204 // MAP                                            170 // MAP
205                                                << 
206 if (fModel==2)                                    171 if (fModel==2)
207 {                                                 172 {
208   Bfield[0] = 0.0;                                173   Bfield[0] = 0.0;
209   Bfield[1] = 0.0;                                174   Bfield[1] = 0.0;
210   Bfield[2] = 0.0;                                175   Bfield[2] = 0.0;
211   Bfield[3] = 0.0;                                176   Bfield[3] = 0.0;
212   Bfield[4] = 0.0;                                177   Bfield[4] = 0.0;
213   Bfield[5] = 0.0;                                178   Bfield[5] = 0.0;
214                                                   179 
215   double x = point[0];                            180   double x = point[0];
216   double y = point[1];                            181   double y = point[1];
217   double z = point[2];                            182   double z = point[2]; 
218                                                   183 
219   G4int quad;                                     184   G4int quad;
220   G4double gradient[5];                           185   G4double gradient[5];
221                                                   186 
222   gradient[0]=fGradient1*(tesla/m)/coef;          187   gradient[0]=fGradient1*(tesla/m)/coef;
223   gradient[1]=fGradient2*(tesla/m)/coef;          188   gradient[1]=fGradient2*(tesla/m)/coef;
224   gradient[2]=fGradient3*(tesla/m)/coef;          189   gradient[2]=fGradient3*(tesla/m)/coef; 
225   gradient[3]=fGradient4*(tesla/m)/coef;          190   gradient[3]=fGradient4*(tesla/m)/coef;
226   gradient[4]=-fGradient3*(tesla/m)/coef;         191   gradient[4]=-fGradient3*(tesla/m)/coef;
227                                                   192 
228   for (quad=0; quad<=4; quad++)                   193   for (quad=0; quad<=4; quad++)
229   {                                               194   {
230   if ((quad+1)==1) {z = point[2] + 3720 * mm;}    195   if ((quad+1)==1) {z = point[2] + 3720 * mm;}
231   if ((quad+1)==2) {z = point[2] + 3580 * mm;}    196   if ((quad+1)==2) {z = point[2] + 3580 * mm;}
232   if ((quad+1)==3) {z = point[2] + 330  * mm;}    197   if ((quad+1)==3) {z = point[2] + 330  * mm;}
233   if ((quad+1)==4) {z = point[2] + 190  * mm;}    198   if ((quad+1)==4) {z = point[2] + 190  * mm;}
234   if ((quad+1)==5) {z = point[2] + 50   * mm;}    199   if ((quad+1)==5) {z = point[2] + 50   * mm;}
235                                                   200 
236   // Check that the point is within the define    201   // Check that the point is within the defined region 
237                                                   202        
238   if                                              203   if 
239   (                                               204   (
240     x>=fMinix && x<=fMaxix &&                     205     x>=fMinix && x<=fMaxix &&
241     y>=fMiniy && y<=fMaxiy &&                     206     y>=fMiniy && y<=fMaxiy &&
242     z>=fMiniz && z<=fMaxiz                        207     z>=fMiniz && z<=fMaxiz 
243   )                                               208   ) 
244   {                                               209   {
245     // Position of given point within region,     210     // Position of given point within region, normalized to the range
246     // [0,1]                                      211     // [0,1]
247     double xfraction = (x - fMinix) / fDx;        212     double xfraction = (x - fMinix) / fDx;
248     double yfraction = (y - fMiniy) / fDy;        213     double yfraction = (y - fMiniy) / fDy; 
249     double zfraction = (z - fMiniz) / fDz;        214     double zfraction = (z - fMiniz) / fDz;
250                                                   215 
251     // Need addresses of these to pass to modf    216     // Need addresses of these to pass to modf below.
252     // modf uses its second argument as an OUT    217     // modf uses its second argument as an OUTPUT argument.
253     double xdindex, ydindex, zdindex;             218     double xdindex, ydindex, zdindex;
254                                                   219     
255     // Position of the point within the cuboid    220     // Position of the point within the cuboid defined by the
256     // nearest surrounding tabulated points       221     // nearest surrounding tabulated points
257     double xlocal = ( std::modf(xfraction*(fNx    222     double xlocal = ( std::modf(xfraction*(fNx-1), &xdindex));
258     double ylocal = ( std::modf(yfraction*(fNy    223     double ylocal = ( std::modf(yfraction*(fNy-1), &ydindex));
259     double zlocal = ( std::modf(zfraction*(fNz    224     double zlocal = ( std::modf(zfraction*(fNz-1), &zdindex));
260                                                   225     
261     // The indices of the nearest tabulated po    226     // The indices of the nearest tabulated point whose coordinates
262     // are all less than those of the given po    227     // are all less than those of the given point
263                                                << 228     int xindex = static_cast<int>(xdindex);
264     //int xindex = static_cast<int>(xdindex);  << 229     int yindex = static_cast<int>(ydindex);
265     //int yindex = static_cast<int>(ydindex);  << 230     int zindex = static_cast<int>(zdindex);
266     //int zindex = static_cast<int>(zdindex);  << 
267                                                << 
268     // SI 15/12/2016: modified according to bu << 
269     int xindex = static_cast<int>(std::floor(x << 
270     int yindex = static_cast<int>(std::floor(y << 
271     int zindex = static_cast<int>(std::floor(z << 
272                                                   231     
273     // Interpolated field                         232     // Interpolated field
274     Bfield[0] =                                   233     Bfield[0] =
275      (fXField[xindex  ][yindex  ][zindex  ] *     234      (fXField[xindex  ][yindex  ][zindex  ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
276       fXField[xindex  ][yindex  ][zindex+1] *     235       fXField[xindex  ][yindex  ][zindex+1] * (1-xlocal) * (1-ylocal) *    zlocal  +
277       fXField[xindex  ][yindex+1][zindex  ] *     236       fXField[xindex  ][yindex+1][zindex  ] * (1-xlocal) *    ylocal  * (1-zlocal) +
278       fXField[xindex  ][yindex+1][zindex+1] *     237       fXField[xindex  ][yindex+1][zindex+1] * (1-xlocal) *    ylocal  *    zlocal  +
279       fXField[xindex+1][yindex  ][zindex  ] *     238       fXField[xindex+1][yindex  ][zindex  ] *    xlocal  * (1-ylocal) * (1-zlocal) +
280       fXField[xindex+1][yindex  ][zindex+1] *     239       fXField[xindex+1][yindex  ][zindex+1] *    xlocal  * (1-ylocal) *    zlocal  +
281       fXField[xindex+1][yindex+1][zindex  ] *     240       fXField[xindex+1][yindex+1][zindex  ] *    xlocal  *    ylocal  * (1-zlocal) +
282       fXField[xindex+1][yindex+1][zindex+1] *     241       fXField[xindex+1][yindex+1][zindex+1] *    xlocal  *    ylocal  *    zlocal)*gradient[quad]
283       + Bfield[0];                                242       + Bfield[0];
284                                                   243       
285     Bfield[1] =                                   244     Bfield[1] =
286      (fYField[xindex  ][yindex  ][zindex  ] *     245      (fYField[xindex  ][yindex  ][zindex  ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
287       fYField[xindex  ][yindex  ][zindex+1] *     246       fYField[xindex  ][yindex  ][zindex+1] * (1-xlocal) * (1-ylocal) *    zlocal  +
288       fYField[xindex  ][yindex+1][zindex  ] *     247       fYField[xindex  ][yindex+1][zindex  ] * (1-xlocal) *    ylocal  * (1-zlocal) +
289       fYField[xindex  ][yindex+1][zindex+1] *     248       fYField[xindex  ][yindex+1][zindex+1] * (1-xlocal) *    ylocal  *    zlocal  +
290       fYField[xindex+1][yindex  ][zindex  ] *     249       fYField[xindex+1][yindex  ][zindex  ] *    xlocal  * (1-ylocal) * (1-zlocal) +
291       fYField[xindex+1][yindex  ][zindex+1] *     250       fYField[xindex+1][yindex  ][zindex+1] *    xlocal  * (1-ylocal) *    zlocal  +
292       fYField[xindex+1][yindex+1][zindex  ] *     251       fYField[xindex+1][yindex+1][zindex  ] *    xlocal  *    ylocal  * (1-zlocal) +
293       fYField[xindex+1][yindex+1][zindex+1] *     252       fYField[xindex+1][yindex+1][zindex+1] *    xlocal  *    ylocal  *    zlocal)*gradient[quad] 
294       + Bfield[1];                                253       + Bfield[1];
295                                                   254 
296     Bfield[2] =                                   255     Bfield[2] =
297      (fZField[xindex  ][yindex  ][zindex  ] *     256      (fZField[xindex  ][yindex  ][zindex  ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
298       fZField[xindex  ][yindex  ][zindex+1] *     257       fZField[xindex  ][yindex  ][zindex+1] * (1-xlocal) * (1-ylocal) *    zlocal  +
299       fZField[xindex  ][yindex+1][zindex  ] *     258       fZField[xindex  ][yindex+1][zindex  ] * (1-xlocal) *    ylocal  * (1-zlocal) +
300       fZField[xindex  ][yindex+1][zindex+1] *     259       fZField[xindex  ][yindex+1][zindex+1] * (1-xlocal) *    ylocal  *    zlocal  +
301       fZField[xindex+1][yindex  ][zindex  ] *     260       fZField[xindex+1][yindex  ][zindex  ] *    xlocal  * (1-ylocal) * (1-zlocal) +
302       fZField[xindex+1][yindex  ][zindex+1] *     261       fZField[xindex+1][yindex  ][zindex+1] *    xlocal  * (1-ylocal) *    zlocal  +
303       fZField[xindex+1][yindex+1][zindex  ] *     262       fZField[xindex+1][yindex+1][zindex  ] *    xlocal  *    ylocal  * (1-zlocal) +
304       fZField[xindex+1][yindex+1][zindex+1] *     263       fZField[xindex+1][yindex+1][zindex+1] *    xlocal  *    ylocal  *    zlocal)*gradient[quad]
305       + Bfield[2];                                264       + Bfield[2];
306                                                   265 
307      }                                            266      } 
308                                                   267 
309 } // loop on quads                                268 } // loop on quads
310                                                   269 
311 } //end  MAP                                      270 } //end  MAP
312                                                   271 
313                                                   272 
314 //********************************************    273 //******************************************************************
315 // SQUARE                                         274 // SQUARE
316                                                   275 
317 if (fModel==1)                                    276 if (fModel==1)
318 {                                                 277 {
319   Bfield[0] = 0.0;                                278   Bfield[0] = 0.0;
320   Bfield[1] = 0.0;                                279   Bfield[1] = 0.0;
321   Bfield[2] = 0.0;                                280   Bfield[2] = 0.0;
322   Bfield[3] = 0.0;                                281   Bfield[3] = 0.0;
323   Bfield[4] = 0.0;                                282   Bfield[4] = 0.0;
324   Bfield[5] = 0.0;                                283   Bfield[5] = 0.0;
325                                                   284 
326   // Field components                             285   // Field components 
327   G4double Bx = 0;                                286   G4double Bx = 0;
328   G4double By = 0;                                287   G4double By = 0;
329   G4double Bz = 0;                                288   G4double Bz = 0;
330                                                   289    
331   G4double x = point[0];                          290   G4double x = point[0];
332   G4double y = point[1];                          291   G4double y = point[1];
333   G4double z = point[2];                          292   G4double z = point[2];
334                                                   293 
335   if (z>=-3770*mm && z<=-3670*mm)  G0 = (fGrad    294   if (z>=-3770*mm && z<=-3670*mm)  G0 = (fGradient1/coef)* tesla/m;
336   if (z>=-3630*mm && z<=-3530*mm)  G0 = (fGrad    295   if (z>=-3630*mm && z<=-3530*mm)  G0 = (fGradient2/coef)* tesla/m;
337                                                   296   
338   if (z>=-380*mm  && z<=-280*mm)   G0 = (fGrad    297   if (z>=-380*mm  && z<=-280*mm)   G0 = (fGradient3/coef)* tesla/m;
339   if (z>=-240*mm  && z<=-140*mm)   G0 = (fGrad    298   if (z>=-240*mm  && z<=-140*mm)   G0 = (fGradient4/coef)* tesla/m;
340   if (z>=-100*mm  && z<=0*mm)      G0 = (-fGra    299   if (z>=-100*mm  && z<=0*mm)      G0 = (-fGradient3/coef)* tesla/m;
341                                                   300 
342   Bx = y*G0;                                      301   Bx = y*G0;
343   By = x*G0;                                      302   By = x*G0;
344   Bz = 0;                                         303   Bz = 0;
345                                                   304 
346   Bfield[0] = Bx;                                 305   Bfield[0] = Bx;
347   Bfield[1] = By;                                 306   Bfield[1] = By;
348   Bfield[2] = Bz;                                 307   Bfield[2] = Bz;
349                                                   308 
350 }                                                 309 }
351                                                   310 
352 // end SQUARE                                     311 // end SQUARE
353                                                   312 
354 //********************************************    313 //******************************************************************
355 // ENGE                                           314 // ENGE
356                                                   315 
357 if (fModel==3)                                    316 if (fModel==3)
358 {                                                 317 {
359                                                   318 
360   // X POSITION OF FIRST QUADRUPOLE               319   // X POSITION OF FIRST QUADRUPOLE
361   // G4double lineX = 0*mm;                       320   // G4double lineX = 0*mm;
362                                                   321 
363   // Z POSITION OF FIRST QUADRUPOLE               322   // Z POSITION OF FIRST QUADRUPOLE
364   G4double lineZ = -3720*mm;                      323   G4double lineZ = -3720*mm;
365                                                   324 
366   // QUADRUPOLE HALF LENGTH                       325   // QUADRUPOLE HALF LENGTH
367   // G4double quadHalfLength = 50*mm;             326   // G4double quadHalfLength = 50*mm;
368                                                   327   
369   // QUADRUPOLE CENTER COORDINATES                328   // QUADRUPOLE CENTER COORDINATES
370   G4double zoprime;                               329   G4double zoprime;
371                                                   330   
372   G4double Grad1, Grad2, Grad3, Grad4, Grad5;     331   G4double Grad1, Grad2, Grad3, Grad4, Grad5;
373   Grad1=fGradient1;                               332   Grad1=fGradient1;
374   Grad2=fGradient2;                               333   Grad2=fGradient2;
375   Grad3=fGradient3;                               334   Grad3=fGradient3;
376   Grad4=fGradient4;                               335   Grad4=fGradient4;
377   Grad5=-Grad3;                                   336   Grad5=-Grad3;  
378                                                   337 
379   Bfield[0] = 0.0;                                338   Bfield[0] = 0.0;
380   Bfield[1] = 0.0;                                339   Bfield[1] = 0.0;
381   Bfield[2] = 0.0;                                340   Bfield[2] = 0.0;
382   Bfield[3] = 0.0;                                341   Bfield[3] = 0.0;
383   Bfield[4] = 0.0;                                342   Bfield[4] = 0.0;
384   Bfield[5] = 0.0;                                343   Bfield[5] = 0.0;
385                                                   344 
386   double x = point[0];                            345   double x = point[0];
387   double y = point[1];                            346   double y = point[1];
388   double z = point[2];                            347   double z = point[2]; 
389                                                   348 
390   if ( (z>=-3900*mm && z<-3470*mm)  || (z>=-49    349   if ( (z>=-3900*mm && z<-3470*mm)  || (z>=-490*mm && z<100*mm) )
391   {                                               350   {
392   G4double Bx=0;                                  351   G4double Bx=0;
393   G4double By=0;                                  352   G4double By=0; 
394   G4double Bz=0;                                  353   G4double Bz=0;
395                                                   354   
396   // FRINGING FILED CONSTANTS                     355   // FRINGING FILED CONSTANTS
397   G4double c0[5], c1[5], c2[5], z1[5], z2[5],     356   G4double c0[5], c1[5], c2[5], z1[5], z2[5], a0[5], gradient[5];
398                                                   357   
399   // DOUBLET***************                       358   // DOUBLET***************
400                                                   359   
401   // QUADRUPOLE 1                                 360   // QUADRUPOLE 1
402   c0[0] = -10.;     // Ci are constants in Pn(    361   c0[0] = -10.;     // Ci are constants in Pn(z)=C0+C1*s+C2*s^2
403   c1[0] = 3.08874;                                362   c1[0] = 3.08874;
404   c2[0] = -0.00618654;                            363   c2[0] = -0.00618654;
405   z1[0] = 28.6834*mm;   // Fringing field lowe    364   z1[0] = 28.6834*mm;   // Fringing field lower limit
406   z2[0] = z1[0]+50*mm;    // Fringing field up    365   z2[0] = z1[0]+50*mm;    // Fringing field upper limit 
407   a0[0] = 7.5*mm;               // Bore Radius    366   a0[0] = 7.5*mm;               // Bore Radius
408   gradient[0] =Grad1*(tesla/m)/coef;              367   gradient[0] =Grad1*(tesla/m)/coef;           
409                                                   368 
410   // QUADRUPOLE 2                                 369   // QUADRUPOLE 2
411   c0[1] = -10.;     // Ci are constants in Pn(    370   c0[1] = -10.;     // Ci are constants in Pn(z)=C0+C1*s+C2*s^2
412   c1[1] = 3.08874;                                371   c1[1] = 3.08874;
413   c2[1] = -0.00618654;                            372   c2[1] = -0.00618654;
414   z1[1] = 28.6834*mm;     // Fringing field lo    373   z1[1] = 28.6834*mm;     // Fringing field lower limit
415   z2[1] = z1[1]+50*mm;    // Fringing field up    374   z2[1] = z1[1]+50*mm;    // Fringing field upper limit
416   a0[1] = 7.5*mm;               // Bore Radius    375   a0[1] = 7.5*mm;               // Bore Radius
417   gradient[1] =Grad2*(tesla/m)/coef;              376   gradient[1] =Grad2*(tesla/m)/coef;
418                                                   377     
419   // TRIPLET**********                            378   // TRIPLET**********
420                                                   379 
421   // QUADRUPOLE 3                                 380   // QUADRUPOLE 3
422   c0[2] = -10.;     // Ci are constants in Pn(    381   c0[2] = -10.;     // Ci are constants in Pn(z)=C0+C1*s+C2*s^2
423   c1[2] = 3.08874;                                382   c1[2] = 3.08874;
424   c2[2] = -0.00618654;                            383   c2[2] = -0.00618654;
425   z1[2] = 28.6834*mm;   // Fringing field lowe    384   z1[2] = 28.6834*mm;   // Fringing field lower limit
426   z2[2] = z1[2]+50*mm;    // Fringing field up    385   z2[2] = z1[2]+50*mm;    // Fringing field upper limit
427   a0[2] = 7.5*mm;               // Bore Radius    386   a0[2] = 7.5*mm;               // Bore Radius
428   gradient[2] = Grad3*(tesla/m)/coef;             387   gradient[2] = Grad3*(tesla/m)/coef;
429                                                   388 
430   // QUADRUPOLE 4                                 389   // QUADRUPOLE 4
431   c0[3] = -10.;     // Ci are constants in Pn(    390   c0[3] = -10.;     // Ci are constants in Pn(z)=C0+C1*s+C2*s^2
432   c1[3] = 3.08874;                                391   c1[3] = 3.08874;
433   c2[3] = -0.00618654;                            392   c2[3] = -0.00618654;
434   z1[3] = 28.6834*mm;   // Fringing field lowe    393   z1[3] = 28.6834*mm;   // Fringing field lower limit
435   z2[3] = z1[3]+50*mm;    // Fringing field up    394   z2[3] = z1[3]+50*mm;    // Fringing field upper limit
436   a0[3] = 7.5*mm;               // Bore Radius    395   a0[3] = 7.5*mm;               // Bore Radius
437   gradient[3] = Grad4*(tesla/m)/coef;             396   gradient[3] = Grad4*(tesla/m)/coef;
438                                                   397   
439    // QUADRUPOLE 5                                398    // QUADRUPOLE 5
440   c0[4] = -10.;     // Ci are constants in Pn(    399   c0[4] = -10.;     // Ci are constants in Pn(z)=C0+C1*s+C2*s^2
441   c1[4] = 3.08874;                                400   c1[4] = 3.08874;
442   c2[4] = -0.00618654;                            401   c2[4] = -0.00618654;
443   z1[4] = 28.6834*mm;   // Fringing field lowe    402   z1[4] = 28.6834*mm;   // Fringing field lower limit
444   z2[4] = z1[4]+50*mm;    // Fringing field up    403   z2[4] = z1[4]+50*mm;    // Fringing field upper limit
445   a0[4] = 7.5*mm;               // Bore Radius    404   a0[4] = 7.5*mm;               // Bore Radius
446   gradient[4] = Grad5*(tesla/m)/coef;             405   gradient[4] = Grad5*(tesla/m)/coef;  
447                                                   406 
448   // FIELD CREATED BY A QUADRUPOLE IN ITS LOCA    407   // FIELD CREATED BY A QUADRUPOLE IN ITS LOCAL FRAME
449   G4double Bx_local,By_local,Bz_local;            408   G4double Bx_local,By_local,Bz_local;
450   Bx_local = 0; By_local = 0; Bz_local = 0;       409   Bx_local = 0; By_local = 0; Bz_local = 0;
451                                                   410   
452   // QUADRUPOLE FRAME                             411   // QUADRUPOLE FRAME
453   G4double x_local,y_local,z_local;               412   G4double x_local,y_local,z_local;
454   x_local= 0; y_local=0; z_local=0;               413   x_local= 0; y_local=0; z_local=0;
455                                                   414 
456   G4double myVars = 0;          // For Enge fo    415   G4double myVars = 0;          // For Enge formula
457   G4double G1, G2, G3;          // For Enge fo    416   G4double G1, G2, G3;          // For Enge formula
458   G4double K1, K2, K3;    // For Enge formula     417   G4double K1, K2, K3;    // For Enge formula
459   G4double P0, P1, P2, cte; // For Enge formul    418   G4double P0, P1, P2, cte; // For Enge formula
460                                                   419 
461   K1=0;                                           420   K1=0;
462   K2=0;                                           421   K2=0;
463   K3=0;                                           422   K3=0;
464                                                   423 
465   P0=0;                                           424   P0=0;
466   P1=0;                                           425   P1=0;
467   P2=0;                                           426   P2=0;
468                                                   427 
469   G0=0;                                           428   G0=0;
470   G1=0;                                           429   G1=0;
471   G2=0;                                           430   G2=0;
472   G3=0;                                           431   G3=0;
473                                                   432 
474   cte=0;                                          433   cte=0;
475                                                   434   
476   for (G4int i=0;i<5; i++) // LOOP ON MAGNETS     435   for (G4int i=0;i<5; i++) // LOOP ON MAGNETS
477   {                                               436   {
478                                                   437  
479    if (i<2) // (if Doublet)                    << 438   if (i<2) // (if Doublet)
480    {                                           << 439   { 
481      zoprime = lineZ + i*140*mm; // centre of     440      zoprime = lineZ + i*140*mm; // centre of magnet nbr i 
482      x_local = x;                                 441      x_local = x; 
483      y_local = y;                                 442      y_local = y; 
484      z_local = (z - zoprime);                     443      z_local = (z - zoprime);
485    }                                           << 444   }
486    else    // else the current magnet is in th << 445   else    // else the current magnet is in the triplet
487    {                                           << 446   {
488      zoprime = lineZ + i*140*mm +(3150-40)*mm;    447      zoprime = lineZ + i*140*mm +(3150-40)*mm;
489                                                   448 
490      x_local = x;                                 449      x_local = x; 
491      y_local = y;                                 450      y_local = y; 
492      z_local = (z - zoprime);                     451      z_local = (z - zoprime);
493                                                   452   
494    }                                           << 453   }         
495                                                   454    
496    if ( z_local < -z2[i] || z_local > z2[i])      455    if ( z_local < -z2[i] || z_local > z2[i])  // Outside the fringing field
497    {                                              456    {
498     G0=0;                                         457     G0=0;
499     G1=0;                                         458     G1=0;
500     G2=0;                                         459     G2=0;
501     G3=0;                                         460     G3=0;
502    }                                              461    }
503                                                   462    
504    if ( (z_local>=-z1[i]) && (z_local<=z1[i])     463    if ( (z_local>=-z1[i]) && (z_local<=z1[i]) ) // inside the quadrupole but outside the fringefield
505    {                                              464    {
506     G0=gradient[i];                               465     G0=gradient[i];
507     G1=0;                                         466     G1=0;
508     G2=0;                                         467     G2=0;
509     G3=0;                                         468     G3=0;
510    }                                              469    }
511                                                   470    
512    if ( ((z_local>=-z2[i]) && (z_local<-z1[i])    471    if ( ((z_local>=-z2[i]) && (z_local<-z1[i])) ||  ((z_local>z1[i]) && (z_local<=z2[i])) ) // inside the fringefield
513    {                                              472    {
514                                                   473 
515           myVars = ( z_local - z1[i]) / a0[i];    474           myVars = ( z_local - z1[i]) / a0[i];     // se (8) p1397 TNS 51
516           if (z_local<-z1[i])  myVars = ( - z_    475           if (z_local<-z1[i])  myVars = ( - z_local - z1[i]) / a0[i];  // see (9) p1397 TNS 51
517                                                   476 
518                                                   477 
519     P0 = c0[i]+c1[i]*myVars+c2[i]*myVars*myVar    478     P0 = c0[i]+c1[i]*myVars+c2[i]*myVars*myVars;
520                                                   479 
521     P1 = c1[i]/a0[i]+2*c2[i]*(z_local-z1[i])/a    480     P1 = c1[i]/a0[i]+2*c2[i]*(z_local-z1[i])/a0[i]/a0[i]; // dP/fDz
522     if (z_local<-z1[i])  P1 = -c1[i]/a0[i]+2*c    481     if (z_local<-z1[i])  P1 = -c1[i]/a0[i]+2*c2[i]*(z_local+z1[i])/a0[i]/a0[i];
523                                                   482 
524     P2 = 2*c2[i]/a0[i]/a0[i];   // d2P/fDz2       483     P2 = 2*c2[i]/a0[i]/a0[i];   // d2P/fDz2
525                                                   484 
526     cte = 1 + G4Exp(c0[i]);   // (1+e^c0)         485     cte = 1 + G4Exp(c0[i]);   // (1+e^c0)
527                                                   486 
528     K1 = -cte*P1*G4Exp(P0)/( (1+G4Exp(P0))*(1+    487     K1 = -cte*P1*G4Exp(P0)/( (1+G4Exp(P0))*(1+G4Exp(P0)) );  // see (11) p1397 TNS 51
529                                                   488 
530     K2 = -cte*G4Exp(P0)*(         // see (12)     489     K2 = -cte*G4Exp(P0)*(         // see (12) p1397 TNS 51
531      P2/( (1+G4Exp(P0))*(1+G4Exp(P0)) )           490      P2/( (1+G4Exp(P0))*(1+G4Exp(P0)) )
532      +2*P1*K1/(1+G4Exp(P0))/cte                   491      +2*P1*K1/(1+G4Exp(P0))/cte
533      +P1*P1/(1+G4Exp(P0))/(1+G4Exp(P0))           492      +P1*P1/(1+G4Exp(P0))/(1+G4Exp(P0))
534      );                                           493      );                                                            
535                                                   494  
536     K3 = -cte*G4Exp(P0)*(       // see (13) p1    495     K3 = -cte*G4Exp(P0)*(       // see (13) p1397 TNS 51  
537      (3*P2*P1+P1*P1*P1)/(1+G4Exp(P0))/(1+G4Exp    496      (3*P2*P1+P1*P1*P1)/(1+G4Exp(P0))/(1+G4Exp(P0))
538      +4*K1*(P1*P1+P2)/(1+G4Exp(P0))/cte           497      +4*K1*(P1*P1+P2)/(1+G4Exp(P0))/cte
539      +2*P1*(K1*K1/cte/cte+K2/(1+G4Exp(P0))/cte    498      +2*P1*(K1*K1/cte/cte+K2/(1+G4Exp(P0))/cte)
540      );                                           499      );
541                                                   500     
542     G0 = gradient[i]*cte/(1+G4Exp(P0));     //    501     G0 = gradient[i]*cte/(1+G4Exp(P0));     // G = G0*K(z) see (7) p1397 TNS 51
543     G1 = gradient[i]*K1;        // dG/fDz         502     G1 = gradient[i]*K1;        // dG/fDz
544     G2 = gradient[i]*K2;        // d2G/fDz2       503     G2 = gradient[i]*K2;        // d2G/fDz2
545     G3 = gradient[i]*K3;        // d3G/fDz3       504     G3 = gradient[i]*K3;        // d3G/fDz3
546                                                   505 
547    }                                              506    }
548                                                   507     
549    Bx_local = y_local*(G0-(1./12)*(3*x_local*x    508    Bx_local = y_local*(G0-(1./12)*(3*x_local*x_local+y_local*y_local)*G2);  // see (4) p1396 TNS 51
550    By_local = x_local*(G0-(1./12)*(3*y_local*y    509    By_local = x_local*(G0-(1./12)*(3*y_local*y_local+x_local*x_local)*G2);  // see (5) p1396 TNS 51
551    Bz_local = x_local*y_local*(G1-(1./12)*(x_l    510    Bz_local = x_local*y_local*(G1-(1./12)*(x_local*x_local+y_local*y_local)*G3);  // see (6) p1396 TNS 51
552                                                   511 
553    // TOTAL MAGNETIC FIELD                        512    // TOTAL MAGNETIC FIELD
554                                                   513    
555    Bx = Bx + Bx_local ;                           514    Bx = Bx + Bx_local ;
556    By = By + By_local ;                           515    By = By + By_local ;
557    Bz = Bz + Bz_local ;                           516    Bz = Bz + Bz_local ;
558                                                   517 
559                                                   518 
560   } // LOOP ON QUADRUPOLES                        519   } // LOOP ON QUADRUPOLES 
561                                                   520   
562   Bfield[0] = Bx;                                 521   Bfield[0] = Bx;
563   Bfield[1] = By;                                 522   Bfield[1] = By;
564   Bfield[2] = Bz;                                 523   Bfield[2] = Bz;
565   }                                               524   }
566                                                   525   
567                                                   526             
568 } // end ENGE                                     527 } // end ENGE
569                                                   528 
570 }                                                 529 }
571                                                   530