Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/nanobeam/src/TabulatedField3D.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 
 46 TabulatedField3D::TabulatedField3D(G4float gr1, G4float gr2, G4float gr3, G4float gr4, G4int choiceModel) 
 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 
187 void TabulatedField3D::GetFieldValue(const double point[4],
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 }
571