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