| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Please cite the following papers if you use
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_INI
42 }
43
44 //....oooOO0OOooo........oooOO0OOooo........oo
45
46 TabulatedField3D::TabulatedField3D(G4float gr1
47 {
48
49 G4cout << " ********************** " << G4en
50 G4cout << " **** CONFIGURATION *** " << G4en
51 G4cout << " ********************** " << G4en
52
53 G4cout<< G4endl;
54 G4cout << "=====> You have selected :" << G4
55 if (choiceModel==1) G4cout<< "-> Square quad
56 if (choiceModel==2) G4cout<< "-> 3D quadrupo
57 if (choiceModel==3) G4cout<< "-> Enge quadru
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 t
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 so
82 << "\n-------------------------------------
83
84 G4cout << "\n ---> " "Reading the field grid
85 std::ifstream file( filename ); // Open the
86
87 // Read table dimensions
88 file >> fNx >> fNy >> fNz; // Note dodgy ord
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 exa
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 >>
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 " << G4e
145
146 // G4cout << " Read values of field from fil
147
148 G4cout << " ---> assumed the order: x, y, z
149 << "\n ---> Min values x,y,z: "
150 << fMinix/cm << " " << fMiniy/cm << " " <<
151 << "\n ---> Max values x,y,z: "
152 << fMaxix/cm << " " << fMaxiy/cm << " " <<
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
160 << "\n-------------------------------------
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]/1
173 fYField[ix][iy][iz] = (fYField[ix][iy]
174 fZField[ix][iy][iz] = (fZField[ix][iy]
175
176 }
177 }
178 }
179
180 } // fModel==2
181
182 }
183
184 //....oooOO0OOooo........oooOO0OOooo........oo
185
186
187 void TabulatedField3D::GetFieldValue(const dou
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 define
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,
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
252 // modf uses its second argument as an OUT
253 double xdindex, ydindex, zdindex;
254
255 // Position of the point within the cuboid
256 // nearest surrounding tabulated points
257 double xlocal = ( std::modf(xfraction*(fNx
258 double ylocal = ( std::modf(yfraction*(fNy
259 double zlocal = ( std::modf(zfraction*(fNz
260
261 // The indices of the nearest tabulated po
262 // are all less than those of the given po
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 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
273 // Interpolated field
274 Bfield[0] =
275 (fXField[xindex ][yindex ][zindex ] *
276 fXField[xindex ][yindex ][zindex+1] *
277 fXField[xindex ][yindex+1][zindex ] *
278 fXField[xindex ][yindex+1][zindex+1] *
279 fXField[xindex+1][yindex ][zindex ] *
280 fXField[xindex+1][yindex ][zindex+1] *
281 fXField[xindex+1][yindex+1][zindex ] *
282 fXField[xindex+1][yindex+1][zindex+1] *
283 + Bfield[0];
284
285 Bfield[1] =
286 (fYField[xindex ][yindex ][zindex ] *
287 fYField[xindex ][yindex ][zindex+1] *
288 fYField[xindex ][yindex+1][zindex ] *
289 fYField[xindex ][yindex+1][zindex+1] *
290 fYField[xindex+1][yindex ][zindex ] *
291 fYField[xindex+1][yindex ][zindex+1] *
292 fYField[xindex+1][yindex+1][zindex ] *
293 fYField[xindex+1][yindex+1][zindex+1] *
294 + Bfield[1];
295
296 Bfield[2] =
297 (fZField[xindex ][yindex ][zindex ] *
298 fZField[xindex ][yindex ][zindex+1] *
299 fZField[xindex ][yindex+1][zindex ] *
300 fZField[xindex ][yindex+1][zindex+1] *
301 fZField[xindex+1][yindex ][zindex ] *
302 fZField[xindex+1][yindex ][zindex+1] *
303 fZField[xindex+1][yindex+1][zindex ] *
304 fZField[xindex+1][yindex+1][zindex+1] *
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 = (fGrad
336 if (z>=-3630*mm && z<=-3530*mm) G0 = (fGrad
337
338 if (z>=-380*mm && z<=-280*mm) G0 = (fGrad
339 if (z>=-240*mm && z<=-140*mm) G0 = (fGrad
340 if (z>=-100*mm && z<=0*mm) G0 = (-fGra
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>=-49
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],
398
399 // DOUBLET***************
400
401 // QUADRUPOLE 1
402 c0[0] = -10.; // Ci are constants in Pn(
403 c1[0] = 3.08874;
404 c2[0] = -0.00618654;
405 z1[0] = 28.6834*mm; // Fringing field lowe
406 z2[0] = z1[0]+50*mm; // Fringing field up
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(
412 c1[1] = 3.08874;
413 c2[1] = -0.00618654;
414 z1[1] = 28.6834*mm; // Fringing field lo
415 z2[1] = z1[1]+50*mm; // Fringing field up
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(
423 c1[2] = 3.08874;
424 c2[2] = -0.00618654;
425 z1[2] = 28.6834*mm; // Fringing field lowe
426 z2[2] = z1[2]+50*mm; // Fringing field up
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(
432 c1[3] = 3.08874;
433 c2[3] = -0.00618654;
434 z1[3] = 28.6834*mm; // Fringing field lowe
435 z2[3] = z1[3]+50*mm; // Fringing field up
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(
441 c1[4] = 3.08874;
442 c2[4] = -0.00618654;
443 z1[4] = 28.6834*mm; // Fringing field lowe
444 z2[4] = z1[4]+50*mm; // Fringing field up
445 a0[4] = 7.5*mm; // Bore Radius
446 gradient[4] = Grad5*(tesla/m)/coef;
447
448 // FIELD CREATED BY A QUADRUPOLE IN ITS LOCA
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 fo
457 G4double G1, G2, G3; // For Enge fo
458 G4double K1, K2, K3; // For Enge formula
459 G4double P0, P1, P2, cte; // For Enge formul
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
482 x_local = x;
483 y_local = y;
484 z_local = (z - zoprime);
485 }
486 else // else the current magnet is in th
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])
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])
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])
513 {
514
515 myVars = ( z_local - z1[i]) / a0[i];
516 if (z_local<-z1[i]) myVars = ( - z_
517
518
519 P0 = c0[i]+c1[i]*myVars+c2[i]*myVars*myVar
520
521 P1 = c1[i]/a0[i]+2*c2[i]*(z_local-z1[i])/a
522 if (z_local<-z1[i]) P1 = -c1[i]/a0[i]+2*c
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+
529
530 K2 = -cte*G4Exp(P0)*( // see (12)
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) p1
537 (3*P2*P1+P1*P1*P1)/(1+G4Exp(P0))/(1+G4Exp
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)); //
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
550 By_local = x_local*(G0-(1./12)*(3*y_local*y
551 Bz_local = x_local*y_local*(G1-(1./12)*(x_l
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 }
571