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