| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // Please cite the following papers if you use << 26 // -------------------------------------------------------------------
27 // Nucl.Instrum.Meth.B260:20-27, 2007 << 27 // $Id: TabulatedField3D.cc,v 1.4 2009/03/02 16:47:52 gunter Exp $
28 // IEEE TNS 51, 4:1395-1401, 2004 << 28 // -------------------------------------------------------------------
29 // <<
30 // Based on purging magnet advanced example <<
31 // <<
32 29
33 #include "TabulatedField3D.hh" 30 #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 31
44 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45 33
46 TabulatedField3D::TabulatedField3D(G4float gr1 34 TabulatedField3D::TabulatedField3D(G4float gr1, G4float gr2, G4float gr3, G4float gr4, G4int choiceModel)
47 { 35 {
48 36
49 G4cout << " ********************** " << G4en 37 G4cout << " ********************** " << G4endl;
50 G4cout << " **** CONFIGURATION *** " << G4en 38 G4cout << " **** CONFIGURATION *** " << G4endl;
51 G4cout << " ********************** " << G4en 39 G4cout << " ********************** " << G4endl;
52 40
53 G4cout<< G4endl; 41 G4cout<< G4endl;
54 G4cout << "=====> You have selected :" << G4 42 G4cout << "=====> You have selected :" << G4endl;
55 if (choiceModel==1) G4cout<< "-> Square quad 43 if (choiceModel==1) G4cout<< "-> Square quadrupole field"<<G4endl;
56 if (choiceModel==2) G4cout<< "-> 3D quadrupo 44 if (choiceModel==2) G4cout<< "-> 3D quadrupole field"<<G4endl;
57 if (choiceModel==3) G4cout<< "-> Enge quadru 45 if (choiceModel==3) G4cout<< "-> Enge quadrupole field"<<G4endl;
58 G4cout << " G1 (T/m) = "<< gr1 << G4endl; 46 G4cout << " G1 (T/m) = "<< gr1 << G4endl;
59 G4cout << " G2 (T/m) = "<< gr2 << G4endl; 47 G4cout << " G2 (T/m) = "<< gr2 << G4endl;
60 G4cout << " G3 (T/m) = "<< gr3 << G4endl; 48 G4cout << " G3 (T/m) = "<< gr3 << G4endl;
61 G4cout << " G4 (T/m) = "<< gr4 << G4endl; 49 G4cout << " G4 (T/m) = "<< gr4 << G4endl;
62 50
63 fGradient1 = gr1; << 51 gradient1 = gr1;
64 fGradient2 = gr2; << 52 gradient2 = gr2;
65 fGradient3 = gr3; << 53 gradient3 = gr3;
66 fGradient4 = gr4; << 54 gradient4 = gr4;
67 fModel = choiceModel; << 55 model = choiceModel;
68 56
69 if (fModel==2) << 57 if (model==2)
70 { 58 {
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"; 59 const char * filename ="OM50.grid";
78 60
79 double lenUnit= mm; 61 double lenUnit= mm;
80 G4cout << "\n------------------------------- 62 G4cout << "\n-----------------------------------------------------------"
81 << "\n 3D Magnetic field from OPERA so 63 << "\n 3D Magnetic field from OPERA software "
82 << "\n------------------------------------- 64 << "\n-----------------------------------------------------------";
83 65
84 G4cout << "\n ---> " "Reading the field grid << 66 G4cout << "\n ---> " "Reading the field grid from " << filename << " ... " << endl;
85 std::ifstream file( filename ); // Open the << 67 ifstream file( filename ); // Open the file for reading.
86 68
87 // Read table dimensions 69 // Read table dimensions
88 file >> fNx >> fNy >> fNz; // Note dodgy ord << 70 file >> nx >> ny >> nz; // Note dodgy order
89 71
90 G4cout << " [ Number of values x,y,z: " 72 G4cout << " [ Number of values x,y,z: "
91 << fNx << " " << fNy << " " << fNz << " ] " << 73 << nx << " " << ny << " " << nz << " ] "
92 << G4endl; << 74 << endl;
93 75
94 // Set up storage space for table 76 // Set up storage space for table
95 fXField.resize( fNx ); << 77 xField.resize( nx );
96 fYField.resize( fNx ); << 78 yField.resize( nx );
97 fZField.resize( fNx ); << 79 zField.resize( nx );
98 int ix, iy, iz; 80 int ix, iy, iz;
99 for (ix=0; ix<fNx; ix++) << 81 for (ix=0; ix<nx; ix++) {
100 { << 82 xField[ix].resize(ny);
101 fXField[ix].resize(fNy); << 83 yField[ix].resize(ny);
102 fYField[ix].resize(fNy); << 84 zField[ix].resize(ny);
103 fZField[ix].resize(fNy); << 85 for (iy=0; iy<ny; iy++) {
104 for (iy=0; iy<fNy; iy++) << 86 xField[ix][iy].resize(nz);
105 { << 87 yField[ix][iy].resize(nz);
106 fXField[ix][iy].resize(fNz); << 88 zField[ix][iy].resize(nz);
107 fYField[ix][iy].resize(fNz); <<
108 fZField[ix][iy].resize(fNz); <<
109 } 89 }
110 } 90 }
111 91
112 // Read in the data 92 // Read in the data
113 double xval,yval,zval,bx,by,bz; 93 double xval,yval,zval,bx,by,bz;
114 double permeability; // Not used in this exa 94 double permeability; // Not used in this example.
115 for (ix=0; ix<fNx; ix++) << 95 for (ix=0; ix<nx; ix++) {
116 { << 96 for (iy=0; iy<ny; iy++) {
117 for (iy=0; iy<fNy; iy++) << 97 for (iz=0; iz<nz; iz++) {
118 { <<
119 for (iz=0; iz<fNz; iz++) <<
120 { <<
121 file >> xval >> yval >> zval >> bx >> 98 file >> xval >> yval >> zval >> bx >> by >> bz >> permeability;
122 if ( ix==0 && iy==0 && iz==0 ) << 99 if ( ix==0 && iy==0 && iz==0 ) {
123 { << 100 minx = xval * lenUnit;
124 fMinix = xval * lenUnit; << 101 miny = yval * lenUnit;
125 fMiniy = yval * lenUnit; << 102 minz = zval * lenUnit;
126 fMiniz = zval * lenUnit; <<
127 } 103 }
128 fXField[ix][iy][iz] = bx ; << 104 xField[ix][iy][iz] = bx ;
129 fYField[ix][iy][iz] = by ; << 105 yField[ix][iy][iz] = by ;
130 fZField[ix][iy][iz] = bz ; << 106 zField[ix][iy][iz] = bz ;
131 } 107 }
132 } 108 }
133 } 109 }
134 file.close(); 110 file.close();
135 111
136 // << 112 maxx = xval * lenUnit;
137 lock.unlock(); << 113 maxy = yval * lenUnit;
138 // << 114 maxz = zval * lenUnit;
139 <<
140 fMaxix = xval * lenUnit; <<
141 fMaxiy = yval * lenUnit; <<
142 fMaxiz = zval * lenUnit; <<
143 115
144 G4cout << "\n ---> ... done reading " << G4e << 116 G4cout << "\n ---> ... done reading " << endl;
145 <<
146 // G4cout << " Read values of field from fil <<
147 117
>> 118 // G4cout << " Read values of field from file " << filename << endl;
148 G4cout << " ---> assumed the order: x, y, z 119 G4cout << " ---> assumed the order: x, y, z, Bx, By, Bz "
149 << "\n ---> Min values x,y,z: " 120 << "\n ---> Min values x,y,z: "
150 << fMinix/cm << " " << fMiniy/cm << " " << << 121 << minx/cm << " " << miny/cm << " " << minz/cm << " cm "
151 << "\n ---> Max values x,y,z: " 122 << "\n ---> Max values x,y,z: "
152 << fMaxix/cm << " " << fMaxiy/cm << " " << << 123 << maxx/cm << " " << maxy/cm << " " << maxz/cm << " cm " << endl;
153 <<
154 fDx = fMaxix - fMinix; <<
155 fDy = fMaxiy - fMiniy; <<
156 fDz = fMaxiz - fMiniz; <<
157 124
>> 125 dx = maxx - minx;
>> 126 dy = maxy - miny;
>> 127 dz = maxz - minz;
158 G4cout << "\n ---> Dif values x,y,z (range): 128 G4cout << "\n ---> Dif values x,y,z (range): "
159 << fDx/cm << " " << fDy/cm << " " << fDz/cm << 129 << dx/cm << " " << dy/cm << " " << dz/cm << " cm in z "
160 << "\n------------------------------------- << 130 << "\n-----------------------------------------------------------" << endl;
161 131
162 132
163 // Table normalization 133 // Table normalization
164 << 134 for (ix=0; ix<nx; ix++)
165 for (ix=0; ix<fNx; ix++) <<
166 { 135 {
167 for (iy=0; iy<fNy; iy++) << 136 for (iy=0; iy<ny; iy++)
168 { 137 {
169 for (iz=0; iz<fNz; iz++) << 138 for (iz=0; iz<nz; iz++)
170 { 139 {
171 140
172 fXField[ix][iy][iz] = (fXField[ix][iy][iz]/1 << 141 xField[ix][iy][iz] = (xField[ix][iy][iz]/197.736);
173 fYField[ix][iy][iz] = (fYField[ix][iy] << 142 yField[ix][iy][iz] = (yField[ix][iy][iz]/197.736);
174 fZField[ix][iy][iz] = (fZField[ix][iy] << 143 zField[ix][iy][iz] = (zField[ix][iy][iz]/197.736);
175 144
176 } << 145 }
177 } 146 }
178 } 147 }
179 148
180 } // fModel==2 << 149 } // model==2
181 150
182 } 151 }
183 152
184 //....oooOO0OOooo........oooOO0OOooo........oo 153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
185 154
186 155
187 void TabulatedField3D::GetFieldValue(const dou 156 void TabulatedField3D::GetFieldValue(const double point[4],
188 double *Bfield ) const 157 double *Bfield ) const
189 { 158 {
190 //G4cout << fGradient1 << G4endl; <<
191 //G4cout << fGradient2 << G4endl; <<
192 //G4cout << fGradient3 << G4endl; <<
193 //G4cout << fGradient4 << G4endl; <<
194 //G4cout << "---------" << G4endl; <<
195 159
196 G4double coef, G0; 160 G4double coef, G0;
197 G0 = 0; 161 G0 = 0;
198 162
199 coef=1; // for protons << 163 coef=1; //protons
200 //coef=2; // for alphas << 164 //coef=2; // alphas
201 165
202 //******************************************** 166 //******************************************************************
203 167
204 // MAP 168 // MAP
205 << 169 if (model==2)
206 if (fModel==2) <<
207 { 170 {
208 Bfield[0] = 0.0; 171 Bfield[0] = 0.0;
209 Bfield[1] = 0.0; 172 Bfield[1] = 0.0;
210 Bfield[2] = 0.0; 173 Bfield[2] = 0.0;
211 Bfield[3] = 0.0; 174 Bfield[3] = 0.0;
212 Bfield[4] = 0.0; 175 Bfield[4] = 0.0;
213 Bfield[5] = 0.0; 176 Bfield[5] = 0.0;
214 177
215 double x = point[0]; 178 double x = point[0];
216 double y = point[1]; 179 double y = point[1];
217 double z = point[2]; 180 double z = point[2];
218 181
219 G4int quad; 182 G4int quad;
220 G4double gradient[5]; 183 G4double gradient[5];
221 184
222 gradient[0]=fGradient1*(tesla/m)/coef; << 185 gradient[0]=gradient1*(tesla/m)/coef;
223 gradient[1]=fGradient2*(tesla/m)/coef; << 186 gradient[1]=gradient2*(tesla/m)/coef;
224 gradient[2]=fGradient3*(tesla/m)/coef; << 187 gradient[2]=gradient3*(tesla/m)/coef;
225 gradient[3]=fGradient4*(tesla/m)/coef; << 188 gradient[3]=gradient4*(tesla/m)/coef;
226 gradient[4]=-fGradient3*(tesla/m)/coef; << 189 gradient[4]=-gradient3*(tesla/m)/coef;
227 190
228 for (quad=0; quad<=4; quad++) 191 for (quad=0; quad<=4; quad++)
229 { 192 {
230 if ((quad+1)==1) {z = point[2] + 3720 * mm;} 193 if ((quad+1)==1) {z = point[2] + 3720 * mm;}
231 if ((quad+1)==2) {z = point[2] + 3580 * mm;} 194 if ((quad+1)==2) {z = point[2] + 3580 * mm;}
232 if ((quad+1)==3) {z = point[2] + 330 * mm;} 195 if ((quad+1)==3) {z = point[2] + 330 * mm;}
233 if ((quad+1)==4) {z = point[2] + 190 * mm;} 196 if ((quad+1)==4) {z = point[2] + 190 * mm;}
234 if ((quad+1)==5) {z = point[2] + 50 * mm;} 197 if ((quad+1)==5) {z = point[2] + 50 * mm;}
235 198
236 // Check that the point is within the define 199 // Check that the point is within the defined region
237 200
238 if 201 if
239 ( 202 (
240 x>=fMinix && x<=fMaxix && << 203 x>=minx && x<=maxx &&
241 y>=fMiniy && y<=fMaxiy && << 204 y>=miny && y<=maxy &&
242 z>=fMiniz && z<=fMaxiz << 205 z>=minz && z<=maxz
243 ) 206 )
244 { 207 {
245 // Position of given point within region, 208 // Position of given point within region, normalized to the range
246 // [0,1] 209 // [0,1]
247 double xfraction = (x - fMinix) / fDx; << 210 double xfraction = (x - minx) / dx;
248 double yfraction = (y - fMiniy) / fDy; << 211 double yfraction = (y - miny) / dy;
249 double zfraction = (z - fMiniz) / fDz; << 212 double zfraction = (z - minz) / dz;
250 213
251 // Need addresses of these to pass to modf 214 // Need addresses of these to pass to modf below.
252 // modf uses its second argument as an OUT 215 // modf uses its second argument as an OUTPUT argument.
253 double xdindex, ydindex, zdindex; 216 double xdindex, ydindex, zdindex;
254 217
255 // Position of the point within the cuboid 218 // Position of the point within the cuboid defined by the
256 // nearest surrounding tabulated points 219 // nearest surrounding tabulated points
257 double xlocal = ( std::modf(xfraction*(fNx << 220 double xlocal = ( std::modf(xfraction*(nx-1), &xdindex));
258 double ylocal = ( std::modf(yfraction*(fNy << 221 double ylocal = ( std::modf(yfraction*(ny-1), &ydindex));
259 double zlocal = ( std::modf(zfraction*(fNz << 222 double zlocal = ( std::modf(zfraction*(nz-1), &zdindex));
260 223
261 // The indices of the nearest tabulated po 224 // The indices of the nearest tabulated point whose coordinates
262 // are all less than those of the given po 225 // are all less than those of the given point
263 << 226 int xindex = static_cast<int>(xdindex);
264 //int xindex = static_cast<int>(xdindex); << 227 int yindex = static_cast<int>(ydindex);
265 //int yindex = static_cast<int>(ydindex); << 228 int zindex = static_cast<int>(zdindex);
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 229
273 // Interpolated field 230 // Interpolated field
274 Bfield[0] = 231 Bfield[0] =
275 (fXField[xindex ][yindex ][zindex ] * << 232 (xField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
276 fXField[xindex ][yindex ][zindex+1] * << 233 xField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
277 fXField[xindex ][yindex+1][zindex ] * << 234 xField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
278 fXField[xindex ][yindex+1][zindex+1] * << 235 xField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
279 fXField[xindex+1][yindex ][zindex ] * << 236 xField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
280 fXField[xindex+1][yindex ][zindex+1] * << 237 xField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
281 fXField[xindex+1][yindex+1][zindex ] * << 238 xField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
282 fXField[xindex+1][yindex+1][zindex+1] * << 239 xField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal)*gradient[quad]
283 + Bfield[0]; 240 + Bfield[0];
284 241
285 Bfield[1] = 242 Bfield[1] =
286 (fYField[xindex ][yindex ][zindex ] * << 243 (yField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
287 fYField[xindex ][yindex ][zindex+1] * << 244 yField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
288 fYField[xindex ][yindex+1][zindex ] * << 245 yField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
289 fYField[xindex ][yindex+1][zindex+1] * << 246 yField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
290 fYField[xindex+1][yindex ][zindex ] * << 247 yField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
291 fYField[xindex+1][yindex ][zindex+1] * << 248 yField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
292 fYField[xindex+1][yindex+1][zindex ] * << 249 yField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
293 fYField[xindex+1][yindex+1][zindex+1] * << 250 yField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal)*gradient[quad]
294 + Bfield[1]; 251 + Bfield[1];
295 252
296 Bfield[2] = 253 Bfield[2] =
297 (fZField[xindex ][yindex ][zindex ] * << 254 (zField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
298 fZField[xindex ][yindex ][zindex+1] * << 255 zField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
299 fZField[xindex ][yindex+1][zindex ] * << 256 zField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
300 fZField[xindex ][yindex+1][zindex+1] * << 257 zField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
301 fZField[xindex+1][yindex ][zindex ] * << 258 zField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
302 fZField[xindex+1][yindex ][zindex+1] * << 259 zField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
303 fZField[xindex+1][yindex+1][zindex ] * << 260 zField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
304 fZField[xindex+1][yindex+1][zindex+1] * << 261 zField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal)*gradient[quad]
305 + Bfield[2]; 262 + Bfield[2];
306 263
307 } 264 }
308 265
309 } // loop on quads 266 } // loop on quads
310 267
311 } //end MAP 268 } //end MAP
312 269
313 270
314 //******************************************** 271 //******************************************************************
315 // SQUARE 272 // SQUARE
316 273
317 if (fModel==1) << 274 if (model==1)
318 { 275 {
319 Bfield[0] = 0.0; 276 Bfield[0] = 0.0;
320 Bfield[1] = 0.0; 277 Bfield[1] = 0.0;
321 Bfield[2] = 0.0; 278 Bfield[2] = 0.0;
322 Bfield[3] = 0.0; 279 Bfield[3] = 0.0;
323 Bfield[4] = 0.0; 280 Bfield[4] = 0.0;
324 Bfield[5] = 0.0; 281 Bfield[5] = 0.0;
325 282
326 // Field components 283 // Field components
327 G4double Bx = 0; 284 G4double Bx = 0;
328 G4double By = 0; 285 G4double By = 0;
329 G4double Bz = 0; 286 G4double Bz = 0;
330 287
331 G4double x = point[0]; 288 G4double x = point[0];
332 G4double y = point[1]; 289 G4double y = point[1];
333 G4double z = point[2]; 290 G4double z = point[2];
334 291
335 if (z>=-3770*mm && z<=-3670*mm) G0 = (fGrad << 292 if (z>=-3770*mm && z<=-3670*mm) G0 = (gradient1/coef)* tesla/m;
336 if (z>=-3630*mm && z<=-3530*mm) G0 = (fGrad << 293 if (z>=-3630*mm && z<=-3530*mm) G0 = (gradient2/coef)* tesla/m;
337 294
338 if (z>=-380*mm && z<=-280*mm) G0 = (fGrad << 295 if (z>=-380*mm && z<=-280*mm) G0 = (gradient3/coef)* tesla/m;
339 if (z>=-240*mm && z<=-140*mm) G0 = (fGrad << 296 if (z>=-240*mm && z<=-140*mm) G0 = (gradient4/coef)* tesla/m;
340 if (z>=-100*mm && z<=0*mm) G0 = (-fGra << 297 if (z>=-100*mm && z<=0*mm) G0 = (-gradient3/coef)* tesla/m;
341 298
342 Bx = y*G0; 299 Bx = y*G0;
343 By = x*G0; 300 By = x*G0;
344 Bz = 0; 301 Bz = 0;
345 302
346 Bfield[0] = Bx; 303 Bfield[0] = Bx;
347 Bfield[1] = By; 304 Bfield[1] = By;
348 Bfield[2] = Bz; 305 Bfield[2] = Bz;
349 306
350 } 307 }
351 308
352 // end SQUARE 309 // end SQUARE
353 310
354 //******************************************** 311 //******************************************************************
355 // ENGE 312 // ENGE
356 313
357 if (fModel==3) << 314 if (model==3)
358 { 315 {
359 316
360 // X POSITION OF FIRST QUADRUPOLE 317 // X POSITION OF FIRST QUADRUPOLE
361 // G4double lineX = 0*mm; << 318 //G4double lineX = 0*mm;
362 319
363 // Z POSITION OF FIRST QUADRUPOLE 320 // Z POSITION OF FIRST QUADRUPOLE
364 G4double lineZ = -3720*mm; 321 G4double lineZ = -3720*mm;
365 322
366 // QUADRUPOLE HALF LENGTH 323 // QUADRUPOLE HALF LENGTH
367 // G4double quadHalfLength = 50*mm; << 324 //G4double quadHalfLength = 50*mm;
368 325
369 // QUADRUPOLE CENTER COORDINATES 326 // QUADRUPOLE CENTER COORDINATES
370 G4double zoprime; 327 G4double zoprime;
371 328
372 G4double Grad1, Grad2, Grad3, Grad4, Grad5; 329 G4double Grad1, Grad2, Grad3, Grad4, Grad5;
373 Grad1=fGradient1; << 330 Grad1=gradient1;
374 Grad2=fGradient2; << 331 Grad2=gradient2;
375 Grad3=fGradient3; << 332 Grad3=gradient3;
376 Grad4=fGradient4; << 333 Grad4=gradient4;
377 Grad5=-Grad3; 334 Grad5=-Grad3;
378 335
379 Bfield[0] = 0.0; 336 Bfield[0] = 0.0;
380 Bfield[1] = 0.0; 337 Bfield[1] = 0.0;
381 Bfield[2] = 0.0; 338 Bfield[2] = 0.0;
382 Bfield[3] = 0.0; 339 Bfield[3] = 0.0;
383 Bfield[4] = 0.0; 340 Bfield[4] = 0.0;
384 Bfield[5] = 0.0; 341 Bfield[5] = 0.0;
385 342
386 double x = point[0]; 343 double x = point[0];
387 double y = point[1]; 344 double y = point[1];
388 double z = point[2]; 345 double z = point[2];
389 346
390 if ( (z>=-3900*mm && z<-3470*mm) || (z>=-49 347 if ( (z>=-3900*mm && z<-3470*mm) || (z>=-490*mm && z<100*mm) )
391 { 348 {
392 G4double Bx=0; 349 G4double Bx=0;
393 G4double By=0; 350 G4double By=0;
394 G4double Bz=0; 351 G4double Bz=0;
395 352
396 // FRINGING FILED CONSTANTS 353 // FRINGING FILED CONSTANTS
397 G4double c0[5], c1[5], c2[5], z1[5], z2[5], 354 G4double c0[5], c1[5], c2[5], z1[5], z2[5], a0[5], gradient[5];
398 355
399 // DOUBLET*************** 356 // DOUBLET***************
400 <<
401 // QUADRUPOLE 1 357 // QUADRUPOLE 1
402 c0[0] = -10.; // Ci are constants in Pn( << 358 c0[0] = -10.; // Ci are constants in the; Pn(z)=C0+C1*s+C2*s^2
403 c1[0] = 3.08874; 359 c1[0] = 3.08874;
404 c2[0] = -0.00618654; 360 c2[0] = -0.00618654;
405 z1[0] = 28.6834*mm; // Fringing field lowe << 361 z1[0] = 28.6834*mm; //Fringing field lower limit
406 z2[0] = z1[0]+50*mm; // Fringing field up << 362 z2[0] = z1[0]+50*mm; //Fringing field upper limit
407 a0[0] = 7.5*mm; // Bore Radius << 363 a0[0] = 7.5*mm; //Bore Radius
408 gradient[0] =Grad1*(tesla/m)/coef; 364 gradient[0] =Grad1*(tesla/m)/coef;
409 365
410 // QUADRUPOLE 2 366 // QUADRUPOLE 2
411 c0[1] = -10.; // Ci are constants in Pn( << 367 c0[1] = -10.; // Ci are constants in the; Pn(z)=C0+C1*s+C2*s^2
412 c1[1] = 3.08874; 368 c1[1] = 3.08874;
413 c2[1] = -0.00618654; 369 c2[1] = -0.00618654;
414 z1[1] = 28.6834*mm; // Fringing field lo << 370 z1[1] = 28.6834*mm; //Fringing field lower limit
415 z2[1] = z1[1]+50*mm; // Fringing field up << 371 z2[1] = z1[1]+50*mm; //Fringing field upper limit
416 a0[1] = 7.5*mm; // Bore Radius << 372 a0[1] = 7.5*mm; //Bore Radius
417 gradient[1] =Grad2*(tesla/m)/coef; 373 gradient[1] =Grad2*(tesla/m)/coef;
418 374
419 // TRIPLET********** 375 // TRIPLET**********
420 <<
421 // QUADRUPOLE 3 376 // QUADRUPOLE 3
422 c0[2] = -10.; // Ci are constants in Pn( << 377 c0[2] = -10.; // Ci are constants in the; Pn(z)=C0+C1*s+C2*s^2
423 c1[2] = 3.08874; 378 c1[2] = 3.08874;
424 c2[2] = -0.00618654; 379 c2[2] = -0.00618654;
425 z1[2] = 28.6834*mm; // Fringing field lowe << 380 z1[2] = 28.6834*mm; //Fringing field lower limit
426 z2[2] = z1[2]+50*mm; // Fringing field up << 381 z2[2] = z1[2]+50*mm; //Fringing field upper limit
427 a0[2] = 7.5*mm; // Bore Radius << 382 a0[2] = 7.5*mm; //Bore Radius
428 gradient[2] = Grad3*(tesla/m)/coef; 383 gradient[2] = Grad3*(tesla/m)/coef;
429 384
430 // QUADRUPOLE 4 385 // QUADRUPOLE 4
431 c0[3] = -10.; // Ci are constants in Pn( << 386 c0[3] = -10.; // Ci are constants in the; Pn(z)=C0+C1*s+C2*s^2
432 c1[3] = 3.08874; 387 c1[3] = 3.08874;
433 c2[3] = -0.00618654; 388 c2[3] = -0.00618654;
434 z1[3] = 28.6834*mm; // Fringing field lowe << 389 z1[3] = 28.6834*mm; //Fringing field lower limit
435 z2[3] = z1[3]+50*mm; // Fringing field up << 390 z2[3] = z1[3]+50*mm; //Fringing field upper limit
436 a0[3] = 7.5*mm; // Bore Radius << 391 a0[3] = 7.5*mm; //Bore Radius
437 gradient[3] = Grad4*(tesla/m)/coef; 392 gradient[3] = Grad4*(tesla/m)/coef;
438 393
439 // QUADRUPOLE 5 394 // QUADRUPOLE 5
440 c0[4] = -10.; // Ci are constants in Pn( << 395 c0[4] = -10.; // Ci are constants in the; Pn(z)=C0+C1*s+C2*s^2
441 c1[4] = 3.08874; 396 c1[4] = 3.08874;
442 c2[4] = -0.00618654; 397 c2[4] = -0.00618654;
443 z1[4] = 28.6834*mm; // Fringing field lowe << 398 z1[4] = 28.6834*mm; //Fringing field lower limit
444 z2[4] = z1[4]+50*mm; // Fringing field up << 399 z2[4] = z1[4]+50*mm; //Fringing field upper limit
445 a0[4] = 7.5*mm; // Bore Radius << 400 a0[4] = 7.5*mm; //Bore Radius
446 gradient[4] = Grad5*(tesla/m)/coef; 401 gradient[4] = Grad5*(tesla/m)/coef;
447 402
448 // FIELD CREATED BY A QUADRUPOLE IN ITS LOCA 403 // FIELD CREATED BY A QUADRUPOLE IN ITS LOCAL FRAME
449 G4double Bx_local,By_local,Bz_local; 404 G4double Bx_local,By_local,Bz_local;
450 Bx_local = 0; By_local = 0; Bz_local = 0; 405 Bx_local = 0; By_local = 0; Bz_local = 0;
451 406
>> 407 // FIELD CREATED BY A QUADRUPOOLE IN WORLD FRAME
>> 408 G4double Bx_quad,By_quad,Bz_quad;
>> 409 Bx_quad = 0; By_quad=0; Bz_quad=0;
>> 410
452 // QUADRUPOLE FRAME 411 // QUADRUPOLE FRAME
453 G4double x_local,y_local,z_local; 412 G4double x_local,y_local,z_local;
454 x_local= 0; y_local=0; z_local=0; 413 x_local= 0; y_local=0; z_local=0;
455 414
456 G4double myVars = 0; // For Enge fo << 415 G4double s = 0; // For Enges formula
457 G4double G1, G2, G3; // For Enge fo << 416 G4double G0, G1, G2, G3; // For Enges formula
458 G4double K1, K2, K3; // For Enge formula << 417 G4double K0, K1, K2, K3; // For Enges formula
459 G4double P0, P1, P2, cte; // For Enge formul << 418 G4double P0, P1, P2, P3, cte; // For Enges formula
460 419
>> 420 K0=0;
461 K1=0; 421 K1=0;
462 K2=0; 422 K2=0;
463 K3=0; 423 K3=0;
464 <<
465 P0=0; 424 P0=0;
466 P1=0; 425 P1=0;
467 P2=0; 426 P2=0;
468 << 427 P3=0;
469 G0=0; 428 G0=0;
470 G1=0; 429 G1=0;
471 G2=0; 430 G2=0;
472 G3=0; 431 G3=0;
473 <<
474 cte=0; 432 cte=0;
475 433
476 for (G4int i=0;i<5; i++) // LOOP ON MAGNETS 434 for (G4int i=0;i<5; i++) // LOOP ON MAGNETS
477 { 435 {
478 436
479 if (i<2) // (if Doublet) << 437 if (i<2) // (if Doublet)
480 { << 438 {
481 zoprime = lineZ + i*140*mm; // centre of << 439 zoprime = lineZ + i*140*mm; //centre of magnet nbr i
482 x_local = x; << 440 x_local = x;
483 y_local = y; << 441 y_local = y;
484 z_local = (z - zoprime); << 442 z_local = (z - zoprime);
485 } << 443 }
486 else // else the current magnet is in th << 444 else // else the current magnet is in the triplet
487 { << 445 {
488 zoprime = lineZ + i*140*mm +(3150-40)*mm; << 446 zoprime = lineZ + i*140*mm +(3150-40)*mm;
489 << 447
490 x_local = x; << 448 x_local = x;
491 y_local = y; << 449 y_local = y;
492 z_local = (z - zoprime); << 450 z_local = (z - zoprime);
493 451
494 } << 452 }
>> 453
495 454
496 if ( z_local < -z2[i] || z_local > z2[i]) 455 if ( z_local < -z2[i] || z_local > z2[i]) // Outside the fringing field
497 { 456 {
498 G0=0; 457 G0=0;
499 G1=0; 458 G1=0;
500 G2=0; 459 G2=0;
501 G3=0; 460 G3=0;
502 } 461 }
503 462
504 if ( (z_local>=-z1[i]) && (z_local<=z1[i]) << 463
>> 464 if ( (z_local>=-z1[i]) & (z_local<=z1[i]) ) // inside the quadrupole but outside the fringefield
505 { 465 {
506 G0=gradient[i]; 466 G0=gradient[i];
507 G1=0; 467 G1=0;
508 G2=0; 468 G2=0;
509 G3=0; 469 G3=0;
510 } 470 }
511 471
512 if ( ((z_local>=-z2[i]) && (z_local<-z1[i]) << 472 if ( ((z_local>=-z2[i]) & (z_local<-z1[i])) || ((z_local>z1[i]) & (z_local<=z2[i])) ) // Inside the fringefield
513 { 473 {
514 474
515 myVars = ( z_local - z1[i]) / a0[i]; << 475 s = ( z_local - z1[i]) / a0[i]; // se (8)
516 if (z_local<-z1[i]) myVars = ( - z_ << 476 if (z_local<-z1[i]) s = ( - z_local - z1[i]) / a0[i]; // se (9) p1397 Incerti et.al.
>> 477
517 478
>> 479 P0 = c0[i]+c1[i]*s+c2[i]*s*s;
518 480
519 P0 = c0[i]+c1[i]*myVars+c2[i]*myVars*myVar << 481 P1 = c1[i]/a0[i]+2*c2[i]*(z_local-z1[i])/a0[i]/a0[i]; //dP/dz
>> 482 if (z_local<-z1[i]) P1 = -c1[i]/a0[i]+2*c2[i]*(z_local+z1[i])/a0[i]/a0[i]; // --"--
520 483
521 P1 = c1[i]/a0[i]+2*c2[i]*(z_local-z1[i])/a << 484 P2 = 2*c2[i]/a0[i]/a0[i]; // d2P/dz2
522 if (z_local<-z1[i]) P1 = -c1[i]/a0[i]+2*c <<
523 485
524 P2 = 2*c2[i]/a0[i]/a0[i]; // d2P/fDz2 << 486 P3 = 0; // d3P/dw3 ??
525 487
526 cte = 1 + G4Exp(c0[i]); // (1+e^c0) << 488 cte = 1 + std::exp(c0[i]); // (1+e^c0)
527 489
528 K1 = -cte*P1*G4Exp(P0)/( (1+G4Exp(P0))*(1+ << 490 K1 = -cte*P1*std::exp(P0)/( (1+std::exp(P0))*(1+std::exp(P0)) ); // se (11) p1397 Incerti et.al.
529 491
530 K2 = -cte*G4Exp(P0)*( // see (12) << 492 K2 = -cte*std::exp(P0)*( // se (12) p1397 Incerti et.al.
531 P2/( (1+G4Exp(P0))*(1+G4Exp(P0)) ) << 493 P2/( (1+std::exp(P0))*(1+std::exp(P0)) )
532 +2*P1*K1/(1+G4Exp(P0))/cte << 494 +2*P1*K1/(1+std::exp(P0))/cte
533 +P1*P1/(1+G4Exp(P0))/(1+G4Exp(P0)) << 495 +P1*P1/(1+std::exp(P0))/(1+std::exp(P0))
534 ); << 496 );
535 497
536 K3 = -cte*G4Exp(P0)*( // see (13) p1 << 498 K3 = -cte*std::exp(P0)*( // se (13) p1397 Incerti et.al
537 (3*P2*P1+P1*P1*P1)/(1+G4Exp(P0))/(1+G4Exp << 499 (3*P2*P1+P1*P1*P1)/(1+std::exp(P0))/(1+std::exp(P0))
538 +4*K1*(P1*P1+P2)/(1+G4Exp(P0))/cte << 500 +4*K1*(P1*P1+P2)/(1+std::exp(P0))/cte
539 +2*P1*(K1*K1/cte/cte+K2/(1+G4Exp(P0))/cte << 501 +2*P1*(K1*K1/cte/cte+K2/(1+std::exp(P0))/cte)
540 ); << 502 );
541 503
542 G0 = gradient[i]*cte/(1+G4Exp(P0)); // << 504 G0 = gradient[i]*cte/(1+std::exp(P0)); // G = G0*K(z) , se (7) p1397 Incerti et.al
543 G1 = gradient[i]*K1; // dG/fDz << 505 G1 = gradient[i]*K1; // dG/dz
544 G2 = gradient[i]*K2; // d2G/fDz2 << 506 G2 = gradient[i]*K2; // d2G/dz2
545 G3 = gradient[i]*K3; // d3G/fDz3 << 507 G3 = gradient[i]*K3; // d3G/dz3
546 508
547 } 509 }
548 510
549 Bx_local = y_local*(G0-(1./12)*(3*x_local*x << 511 Bx_local = y_local*(G0-(1./12)*(3*x_local*x_local+y_local*y_local)*G2); // se (4) p1396 Incerti et.al
550 By_local = x_local*(G0-(1./12)*(3*y_local*y << 512 By_local = x_local*(G0-(1./12)*(3*y_local*y_local+x_local*x_local)*G2); // se (5) p1396 Incerti et.al
551 Bz_local = x_local*y_local*(G1-(1./12)*(x_l << 513 Bz_local = x_local*y_local*(G1-(1./12)*(x_local*x_local+y_local*y_local)*G3); // se (6) p1396 Incerti et.al
552 514
553 // TOTAL MAGNETIC FIELD 515 // TOTAL MAGNETIC FIELD
554 516
555 Bx = Bx + Bx_local ; 517 Bx = Bx + Bx_local ;
556 By = By + By_local ; 518 By = By + By_local ;
557 Bz = Bz + Bz_local ; 519 Bz = Bz + Bz_local ;
558 520
559 521
560 } // LOOP ON QUADRUPOLES 522 } // LOOP ON QUADRUPOLES
561 523
562 Bfield[0] = Bx; 524 Bfield[0] = Bx;
563 Bfield[1] = By; 525 Bfield[1] = By;
564 Bfield[2] = Bz; 526 Bfield[2] = Bz;
565 } 527 }
566 528
567 529
568 } // end ENGE 530 } // end ENGE
569 531
570 } 532 }
571 533