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