Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 // -------------------------------------------------------------------
27 // $Id$
28 // -------------------------------------------------------------------
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  gradient1 = gr1;
53  gradient2 = gr2;
54  gradient3 = gr3;
55  gradient4 = gr4;
56  model = choiceModel;
57 
58  if (model==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 >> nx >> ny >> nz; // Note dodgy order
72 
73  G4cout << " [ Number of values x,y,z: "
74  << nx << " " << ny << " " << nz << " ] "
75  << endl;
76 
77  // Set up storage space for table
78  xField.resize( nx );
79  yField.resize( nx );
80  zField.resize( nx );
81  int ix, iy, iz;
82  for (ix=0; ix<nx; ix++) {
83  xField[ix].resize(ny);
84  yField[ix].resize(ny);
85  zField[ix].resize(ny);
86  for (iy=0; iy<ny; iy++) {
87  xField[ix][iy].resize(nz);
88  yField[ix][iy].resize(nz);
89  zField[ix][iy].resize(nz);
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<nx; ix++) {
97  for (iy=0; iy<ny; iy++) {
98  for (iz=0; iz<nz; iz++) {
99  file >> xval >> yval >> zval >> bx >> by >> bz >> permeability;
100  if ( ix==0 && iy==0 && iz==0 ) {
101  minx = xval * lenUnit;
102  miny = yval * lenUnit;
103  minz = zval * lenUnit;
104  }
105  xField[ix][iy][iz] = bx ;
106  yField[ix][iy][iz] = by ;
107  zField[ix][iy][iz] = bz ;
108  }
109  }
110  }
111  file.close();
112 
113  maxx = xval * lenUnit;
114  maxy = yval * lenUnit;
115  maxz = 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  << minx/cm << " " << miny/cm << " " << minz/cm << " cm "
123  << "\n ---> Max values x,y,z: "
124  << maxx/cm << " " << maxy/cm << " " << maxz/cm << " cm " << endl;
125 
126  dx = maxx - minx;
127  dy = maxy - miny;
128  dz = maxz - minz;
129  G4cout << "\n ---> Dif values x,y,z (range): "
130  << dx/cm << " " << dy/cm << " " << dz/cm << " cm in z "
131  << "\n-----------------------------------------------------------" << endl;
132 
133 
134  // Table normalization
135  for (ix=0; ix<nx; ix++)
136  {
137  for (iy=0; iy<ny; iy++)
138  {
139  for (iz=0; iz<nz; iz++)
140  {
141 
142  xField[ix][iy][iz] = (xField[ix][iy][iz]/197.736);
143  yField[ix][iy][iz] = (yField[ix][iy][iz]/197.736);
144  zField[ix][iy][iz] = (zField[ix][iy][iz]/197.736);
145 
146  }
147  }
148  }
149 
150  } // model==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; //protons
165  //coef=2; // alphas
166 
167 //******************************************************************
168 
169 // MAP
170 if (model==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]=gradient1*(tesla/m)/coef;
187  gradient[1]=gradient2*(tesla/m)/coef;
188  gradient[2]=gradient3*(tesla/m)/coef;
189  gradient[3]=gradient4*(tesla/m)/coef;
190  gradient[4]=-gradient3*(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>=minx && x<=maxx &&
205  y>=miny && y<=maxy &&
206  z>=minz && z<=maxz
207  )
208  {
209  // Position of given point within region, normalized to the range
210  // [0,1]
211  double xfraction = (x - minx) / dx;
212  double yfraction = (y - miny) / dy;
213  double zfraction = (z - minz) / dz;
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*(nx-1), &xdindex));
222  double ylocal = ( std::modf(yfraction*(ny-1), &ydindex));
223  double zlocal = ( std::modf(zfraction*(nz-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  (xField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
234  xField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
235  xField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
236  xField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
237  xField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
238  xField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
239  xField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
240  xField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal)*gradient[quad]
241  + Bfield[0];
242 
243  Bfield[1] =
244  (yField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
245  yField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
246  yField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
247  yField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
248  yField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
249  yField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
250  yField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
251  yField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal)*gradient[quad]
252  + Bfield[1];
253 
254  Bfield[2] =
255  (zField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
256  zField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
257  zField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
258  zField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
259  zField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
260  zField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
261  zField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
262  zField[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 (model==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 = (gradient1/coef)* tesla/m;
294  if (z>=-3630*mm && z<=-3530*mm) G0 = (gradient2/coef)* tesla/m;
295 
296  if (z>=-380*mm && z<=-280*mm) G0 = (gradient3/coef)* tesla/m;
297  if (z>=-240*mm && z<=-140*mm) G0 = (gradient4/coef)* tesla/m;
298  if (z>=-100*mm && z<=0*mm) G0 = (-gradient3/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 (model==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=gradient1;
332  Grad2=gradient2;
333  Grad3=gradient3;
334  Grad4=gradient4;
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  // QUADRUPOLE 1
359  c0[0] = -10.; // Ci are constants in the; Pn(z)=C0+C1*s+C2*s^2
360  c1[0] = 3.08874;
361  c2[0] = -0.00618654;
362  z1[0] = 28.6834*mm; //Fringing field lower limit
363  z2[0] = z1[0]+50*mm; //Fringing field upper limit
364  a0[0] = 7.5*mm; //Bore Radius
365  gradient[0] =Grad1*(tesla/m)/coef;
366 
367  // QUADRUPOLE 2
368  c0[1] = -10.; // Ci are constants in the; Pn(z)=C0+C1*s+C2*s^2
369  c1[1] = 3.08874;
370  c2[1] = -0.00618654;
371  z1[1] = 28.6834*mm; //Fringing field lower limit
372  z2[1] = z1[1]+50*mm; //Fringing field upper limit
373  a0[1] = 7.5*mm; //Bore Radius
374  gradient[1] =Grad2*(tesla/m)/coef;
375 
376  // TRIPLET**********
377  // QUADRUPOLE 3
378  c0[2] = -10.; // Ci are constants in the; Pn(z)=C0+C1*s+C2*s^2
379  c1[2] = 3.08874;
380  c2[2] = -0.00618654;
381  z1[2] = 28.6834*mm; //Fringing field lower limit
382  z2[2] = z1[2]+50*mm; //Fringing field upper limit
383  a0[2] = 7.5*mm; //Bore Radius
384  gradient[2] = Grad3*(tesla/m)/coef;
385 
386  // QUADRUPOLE 4
387  c0[3] = -10.; // Ci are constants in the; Pn(z)=C0+C1*s+C2*s^2
388  c1[3] = 3.08874;
389  c2[3] = -0.00618654;
390  z1[3] = 28.6834*mm; //Fringing field lower limit
391  z2[3] = z1[3]+50*mm; //Fringing field upper limit
392  a0[3] = 7.5*mm; //Bore Radius
393  gradient[3] = Grad4*(tesla/m)/coef;
394 
395  // QUADRUPOLE 5
396  c0[4] = -10.; // Ci are constants in the; Pn(z)=C0+C1*s+C2*s^2
397  c1[4] = 3.08874;
398  c2[4] = -0.00618654;
399  z1[4] = 28.6834*mm; //Fringing field lower limit
400  z2[4] = z1[4]+50*mm; //Fringing field upper limit
401  a0[4] = 7.5*mm; //Bore Radius
402  gradient[4] = Grad5*(tesla/m)/coef;
403 
404  // FIELD CREATED BY A QUADRUPOLE IN ITS LOCAL FRAME
405  G4double Bx_local,By_local,Bz_local;
406  Bx_local = 0; By_local = 0; Bz_local = 0;
407 
408  // FIELD CREATED BY A QUADRUPOOLE IN WORLD FRAME
409  //unused G4double Bx_quad,By_quad,Bz_quad;
410  //unsued Bx_quad = 0; By_quad=0; Bz_quad=0;
411 
412  // QUADRUPOLE FRAME
413  G4double x_local,y_local,z_local;
414  x_local= 0; y_local=0; z_local=0;
415 
416  //G4double vars = 0; // For Enges formula
417  G4double G1, G2, G3; // For Enges formula
418  //G4double K0, K1, K2, K3; // For Enges formula
419  //G4double P0, P1, P2, P3, cte; // For Enges formula
420  G4double K1, K2, K3; // For Enges formula
421  G4double P0, P1, P2, cte; // For Enges formula
422 
423  //K0=0;
424  K1=0;
425  K2=0;
426  K3=0;
427  P0=0;
428  P1=0;
429  P2=0;
430  //P3=0;
431  G0=0;
432  G1=0;
433  G2=0;
434  G3=0;
435  cte=0;
436 
437  for (G4int i=0;i<5; i++) // LOOP ON MAGNETS
438  {
439 
440  if (i<2) // (if Doublet)
441  {
442  zoprime = lineZ + i*140*mm; //centre of magnet nbr i
443  x_local = x;
444  y_local = y;
445  z_local = (z - zoprime);
446  }
447  else // else the current magnet is in the triplet
448  {
449  zoprime = lineZ + i*140*mm +(3150-40)*mm;
450 
451  x_local = x;
452  y_local = y;
453  z_local = (z - zoprime);
454 
455  }
456 
457 
458  if ( z_local < -z2[i] || z_local > z2[i]) // Outside the fringing field
459  {
460  G0=0;
461  G1=0;
462  G2=0;
463  G3=0;
464  }
465 
466 
467  if ( (z_local>=-z1[i]) & (z_local<=z1[i]) ) // inside the quadrupole but outside the fringefield
468  {
469  G0=gradient[i];
470  G1=0;
471  G2=0;
472  G3=0;
473  }
474 
475  if ( ((z_local>=-z2[i]) & (z_local<-z1[i])) || ((z_local>z1[i]) & (z_local<=z2[i])) ) // Inside the fringefield
476  {
477 
478  //vars = ( z_local - z1[i]) / a0[i]; // se (8)
479  //if (z_local<-z1[i]) vars = ( - z_local - z1[i]) / a0[i]; // se (9) p1397 Incerti et.al.
480 
481 
482  P0 = c0[i]+c1[i]*s+c2[i]*s*s;
483 
484  P1 = c1[i]/a0[i]+2*c2[i]*(z_local-z1[i])/a0[i]/a0[i]; //dP/dz
485  if (z_local<-z1[i]) P1 = -c1[i]/a0[i]+2*c2[i]*(z_local+z1[i])/a0[i]/a0[i]; // --"--
486 
487  P2 = 2*c2[i]/a0[i]/a0[i]; // d2P/dz2
488 
489  //P3 = 0; // d3P/dw3 ??
490 
491  cte = 1 + std::exp(c0[i]); // (1+e^c0)
492 
493  K1 = -cte*P1*std::exp(P0)/( (1+std::exp(P0))*(1+std::exp(P0)) ); // se (11) p1397 Incerti et.al.
494 
495  K2 = -cte*std::exp(P0)*( // se (12) p1397 Incerti et.al.
496  P2/( (1+std::exp(P0))*(1+std::exp(P0)) )
497  +2*P1*K1/(1+std::exp(P0))/cte
498  +P1*P1/(1+std::exp(P0))/(1+std::exp(P0))
499  );
500 
501  K3 = -cte*std::exp(P0)*( // se (13) p1397 Incerti et.al
502  (3*P2*P1+P1*P1*P1)/(1+std::exp(P0))/(1+std::exp(P0))
503  +4*K1*(P1*P1+P2)/(1+std::exp(P0))/cte
504  +2*P1*(K1*K1/cte/cte+K2/(1+std::exp(P0))/cte)
505  );
506 
507  G0 = gradient[i]*cte/(1+std::exp(P0)); // G = G0*K(z) , se (7) p1397 Incerti et.al
508  G1 = gradient[i]*K1; // dG/dz
509  G2 = gradient[i]*K2; // d2G/dz2
510  G3 = gradient[i]*K3; // d3G/dz3
511 
512  }
513 
514  Bx_local = y_local*(G0-(1./12)*(3*x_local*x_local+y_local*y_local)*G2); // se (4) p1396 Incerti et.al
515  By_local = x_local*(G0-(1./12)*(3*y_local*y_local+x_local*x_local)*G2); // se (5) p1396 Incerti et.al
516  Bz_local = x_local*y_local*(G1-(1./12)*(x_local*x_local+y_local*y_local)*G3); // se (6) p1396 Incerti et.al
517 
518  // TOTAL MAGNETIC FIELD
519 
520  Bx = Bx + Bx_local ;
521  By = By + By_local ;
522  Bz = Bz + Bz_local ;
523 
524 
525  } // LOOP ON QUADRUPOLES
526 
527  Bfield[0] = Bx;
528  Bfield[1] = By;
529  Bfield[2] = Bz;
530  }
531 
532 
533 } // end ENGE
534 
535 }
Generated on Sat May 25 2013 14:32:14 for Geant4 by   doxygen 1.8.4