| 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 // This example is provided by the Geant4-DNA 26 // This example is provided by the Geant4-DNA collaboration
27 // Any report or published results obtained us 27 // Any report or published results obtained using the Geant4-DNA software
28 // shall cite the following Geant4-DNA collabo 28 // shall cite the following Geant4-DNA collaboration publication:
29 // Med. Phys. 37 (2010) 4692-4708 29 // Med. Phys. 37 (2010) 4692-4708
30 // The Geant4-DNA web site is available at htt 30 // The Geant4-DNA web site is available at http://geant4-dna.org
31 // 31 //
32 // If you use this example, please cite the fo 32 // If you use this example, please cite the following publication:
33 // Rad. Prot. Dos. 133 (2009) 2-11 33 // Rad. Prot. Dos. 133 (2009) 2-11
34 // 34 //
35 // Based on purging magnet advanced example. 35 // Based on purging magnet advanced example.
36 // 36 //
37 37
38 #include "EMField.hh" 38 #include "EMField.hh"
39 #include "G4Exp.hh" 39 #include "G4Exp.hh"
40 #include "G4SystemOfUnits.hh" 40 #include "G4SystemOfUnits.hh"
41 41
>> 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 43
>> 44
42 EMField::EMField() 45 EMField::EMField()
43 {} 46 {}
44 47
>> 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 49
45 void EMField::GetFieldValue(const double point 50 void EMField::GetFieldValue(const double point[4], double *Bfield ) const
46 { 51 {
47 // Magnetic field 52 // Magnetic field
48 Bfield[0] = 0; 53 Bfield[0] = 0;
49 Bfield[1] = 0; 54 Bfield[1] = 0;
50 Bfield[2] = 0; 55 Bfield[2] = 0;
51 56
52 // Electric field 57 // Electric field
53 Bfield[3] = 0; 58 Bfield[3] = 0;
54 Bfield[4] = 0; 59 Bfield[4] = 0;
55 Bfield[5] = 0; 60 Bfield[5] = 0;
56 61
57 G4double Bx = 0; 62 G4double Bx = 0;
58 G4double By = 0; 63 G4double By = 0;
59 G4double Bz = 0; 64 G4double Bz = 0;
60 65
61 G4double x = point[0]; 66 G4double x = point[0];
62 G4double y = point[1]; 67 G4double y = point[1];
63 G4double z = point[2]; 68 G4double z = point[2];
64 69
65 // *********************** 70 // ***********************
66 // AIFIRA SWITCHING MAGNET 71 // AIFIRA SWITCHING MAGNET
67 // *********************** 72 // ***********************
68 73
69 // MAGNETIC FIELD VALUE FOR 3 MeV ALPHAS 74 // MAGNETIC FIELD VALUE FOR 3 MeV ALPHAS
70 G4double switchingField = 0.0589768635 * tes << 75 // G4double switchingField = 0.0589768635 * tesla ;
>> 76 G4double switchingField = 0.0590201 * tesla ;
71 77
72 // BEAM START 78 // BEAM START
73 G4double beamStart = -10*m; 79 G4double beamStart = -10*m;
74 80
75 // RADIUS 81 // RADIUS
76 G4double Rp = 0.698*m; 82 G4double Rp = 0.698*m;
77 83
78 // ENTRANCE POSITION AFTER ANALYSIS MAGNET 84 // ENTRANCE POSITION AFTER ANALYSIS MAGNET
79 G4double zS = 975*mm; 85 G4double zS = 975*mm;
80 86
81 // POLE GAP 87 // POLE GAP
82 G4double D = 31.8*mm; 88 G4double D = 31.8*mm;
83 89
84 // FRINGING FIELD 90 // FRINGING FIELD
85 91
86 G4double fieldBoundary, wc0, wc1, wc2, wc3, 92 G4double fieldBoundary, wc0, wc1, wc2, wc3, limitMinEntrance, limitMaxEntrance, limitMinExit, limitMaxExit;
87 93
88 limitMinEntrance = beamStart+zS-4*D; 94 limitMinEntrance = beamStart+zS-4*D;
89 limitMaxEntrance = beamStart+zS+4*D; 95 limitMaxEntrance = beamStart+zS+4*D;
90 limitMinExit =Rp-4*D; 96 limitMinExit =Rp-4*D;
91 limitMaxExit =Rp+4*D; 97 limitMaxExit =Rp+4*D;
92 98
93 wc0 = 0.3835; 99 wc0 = 0.3835;
94 wc1 = 2.388; 100 wc1 = 2.388;
95 wc2 = -0.8171; 101 wc2 = -0.8171;
96 wc3 = 0.200; 102 wc3 = 0.200;
97 103
98 fieldBoundary=0.62; 104 fieldBoundary=0.62;
99 105
100 G4double ws, largeS, h, dhdlargeS, dhds, dla 106 G4double ws, largeS, h, dhdlargeS, dhds, dlargeSds, dsdz, dsdx, zs0, Rs0, xcenter, zcenter;
101 107
102 // - ENTRANCE OF SWITCHING MAGNET 108 // - ENTRANCE OF SWITCHING MAGNET
103 109
104 if ( (z >= limitMinEntrance) && (z < limitMaxE 110 if ( (z >= limitMinEntrance) && (z < limitMaxEntrance) )
105 { 111 {
106 zs0 = fieldBoundary*D; 112 zs0 = fieldBoundary*D;
107 ws = (-z+beamStart+zS-zs0)/D; 113 ws = (-z+beamStart+zS-zs0)/D;
108 dsdz = -1/D; 114 dsdz = -1/D;
109 dsdx = 0; 115 dsdx = 0;
110 116
111 largeS = wc0 + wc1*ws + wc2*ws*ws + wc3*ws*w 117 largeS = wc0 + wc1*ws + wc2*ws*ws + wc3*ws*ws*ws;
112 h = 1./(1.+G4Exp(largeS)); 118 h = 1./(1.+G4Exp(largeS));
113 dhdlargeS = -G4Exp(largeS)*h*h; 119 dhdlargeS = -G4Exp(largeS)*h*h;
114 dlargeSds = wc1+ 2*wc2*ws + 3*wc3*ws*ws; 120 dlargeSds = wc1+ 2*wc2*ws + 3*wc3*ws*ws;
115 dhds = dhdlargeS * dlargeSds; 121 dhds = dhdlargeS * dlargeSds;
116 122
117 By = switchingField * h ; 123 By = switchingField * h ;
118 Bx = y*switchingField*dhds*dsdx; 124 Bx = y*switchingField*dhds*dsdx;
119 Bz = y*switchingField*dhds*dsdz; 125 Bz = y*switchingField*dhds*dsdz;
120 126
121 } 127 }
122 128
123 // - HEART OF SWITCHING MAGNET 129 // - HEART OF SWITCHING MAGNET
124 130
125 if ( 131 if (
126 (z >= limitMaxEntrance) 132 (z >= limitMaxEntrance)
127 && (( x*x + (z -(beamStart+zS))*(z -(be 133 && (( x*x + (z -(beamStart+zS))*(z -(beamStart+zS)) < limitMinExit*limitMinExit))
128 ) 134 )
129 { 135 {
130 Bx=0; 136 Bx=0;
131 By = switchingField; 137 By = switchingField;
132 Bz=0; 138 Bz=0;
133 } 139 }
134 140
135 // - EXIT OF SWITCHING MAGNET 141 // - EXIT OF SWITCHING MAGNET
136 142
137 if ( 143 if (
138 (z >= limitMaxEntrance) 144 (z >= limitMaxEntrance)
139 && (( x*x + (z -(beamStart+zS))*(z -(beam 145 && (( x*x + (z -(beamStart+zS))*(z -(beamStart+zS))) >= limitMinExit*limitMinExit)
140 && (( x*x + (z -(beamStart+zS))*(z -(beam 146 && (( x*x + (z -(beamStart+zS))*(z -(beamStart+zS))) < limitMaxExit*limitMaxExit)
141 147
142 ) 148 )
143 { 149 {
144 150
145 xcenter = 0; 151 xcenter = 0;
146 zcenter = beamStart+zS; 152 zcenter = beamStart+zS;
147 153
148 Rs0 = Rp + D*fieldBoundary; 154 Rs0 = Rp + D*fieldBoundary;
149 ws = (std::sqrt((z-zcenter)*(z-zcenter)+(x-x 155 ws = (std::sqrt((z-zcenter)*(z-zcenter)+(x-xcenter)*(x-xcenter)) - Rs0)/D;
150 156
151 dsdz = (1/D)*(z-zcenter)/std::sqrt((z-zcente 157 dsdz = (1/D)*(z-zcenter)/std::sqrt((z-zcenter)*(z-zcenter)+(x-xcenter)*(x-xcenter));
152 dsdx = (1/D)*(x-xcenter)/std::sqrt((z-zcente 158 dsdx = (1/D)*(x-xcenter)/std::sqrt((z-zcenter)*(z-zcenter)+(x-xcenter)*(x-xcenter));
153 159
154 largeS = wc0 + wc1*ws + wc2*ws*ws + wc3*ws*w 160 largeS = wc0 + wc1*ws + wc2*ws*ws + wc3*ws*ws*ws;
155 h = 1./(1.+G4Exp(largeS)); 161 h = 1./(1.+G4Exp(largeS));
156 dhdlargeS = -G4Exp(largeS)*h*h; 162 dhdlargeS = -G4Exp(largeS)*h*h;
157 dlargeSds = wc1+ 2*wc2*ws + 3*wc3*ws*ws; 163 dlargeSds = wc1+ 2*wc2*ws + 3*wc3*ws*ws;
158 dhds = dhdlargeS * dlargeSds; 164 dhds = dhdlargeS * dlargeSds;
159 165
160 By = switchingField * h ; 166 By = switchingField * h ;
161 Bx = y*switchingField*dhds*dsdx; 167 Bx = y*switchingField*dhds*dsdx;
162 Bz = y*switchingField*dhds*dsdz; 168 Bz = y*switchingField*dhds*dsdz;
163 169
164 } 170 }
165 171
166 // ************************** 172 // **************************
167 // MICROBEAM LINE QUADRUPOLES 173 // MICROBEAM LINE QUADRUPOLES
168 // ************************** 174 // **************************
169 175
170 // MICROBEAM LINE ANGLE 176 // MICROBEAM LINE ANGLE
171 G4double lineAngle = -10*deg; 177 G4double lineAngle = -10*deg;
172 178
173 // X POSITION OF FIRST QUADRUPOLE 179 // X POSITION OF FIRST QUADRUPOLE
174 G4double lineX = -1295.59*mm; 180 G4double lineX = -1295.59*mm;
175 181
176 // Z POSITION OF FIRST QUADRUPOLE 182 // Z POSITION OF FIRST QUADRUPOLE
177 G4double lineZ = -1327*mm; 183 G4double lineZ = -1327*mm;
178 184
179 // Adjust magnetic zone absolute position 185 // Adjust magnetic zone absolute position
180 lineX = lineX + 5.24*micrometer*std::cos(-li 186 lineX = lineX + 5.24*micrometer*std::cos(-lineAngle); // 5.24 = 1.3 + 3.94 micrometer (cf. DetectorConstruction)
181 lineZ = lineZ + 5.24*micrometer*std::sin(-li 187 lineZ = lineZ + 5.24*micrometer*std::sin(-lineAngle);
182 188
183 // QUADRUPOLE HALF LENGTH 189 // QUADRUPOLE HALF LENGTH
184 G4double quadHalfLength = 75*mm; 190 G4double quadHalfLength = 75*mm;
185 191
186 // QUADRUPOLE SPACING 192 // QUADRUPOLE SPACING
187 G4double quadSpacing = 40*mm; 193 G4double quadSpacing = 40*mm;
188 194
189 // QUADRUPOLE CENTER COORDINATES 195 // QUADRUPOLE CENTER COORDINATES
190 G4double xoprime, zoprime; 196 G4double xoprime, zoprime;
191 197
192 if (z>=-1400*mm && z <-200*mm) 198 if (z>=-1400*mm && z <-200*mm)
193 { 199 {
194 Bx=0; By=0; Bz=0; 200 Bx=0; By=0; Bz=0;
195 201
196 // FRINGING FILED CONSTANTS 202 // FRINGING FILED CONSTANTS
197 G4double c0[4], c1[4], c2[4], z1[4], z2[4], 203 G4double c0[4], c1[4], c2[4], z1[4], z2[4], a0[4], gradient[4];
198 204
199 // QUADRUPOLE 1 205 // QUADRUPOLE 1
200 c0[0] = -5.; 206 c0[0] = -5.;
201 c1[0] = 2.5; 207 c1[0] = 2.5;
202 c2[0] = -0.1; 208 c2[0] = -0.1;
203 z1[0] = 60*mm; 209 z1[0] = 60*mm;
204 z2[0] = 130*mm; 210 z2[0] = 130*mm;
205 a0[0] = 10*mm; 211 a0[0] = 10*mm;
206 gradient[0] = 3.406526 *tesla/m; 212 gradient[0] = 3.406526 *tesla/m;
207 213
208 // QUADRUPOLE 2 214 // QUADRUPOLE 2
209 c0[1] = -5.; 215 c0[1] = -5.;
210 c1[1] = 2.5; 216 c1[1] = 2.5;
211 c2[1] = -0.1; 217 c2[1] = -0.1;
212 z1[1] = 60*mm; 218 z1[1] = 60*mm;
213 z2[1] = 130*mm; 219 z2[1] = 130*mm;
214 a0[1] = 10*mm; 220 a0[1] = 10*mm;
215 gradient[1] = -8.505263 *tesla/m; 221 gradient[1] = -8.505263 *tesla/m;
216 222
217 // QUADRUPOLE 3 223 // QUADRUPOLE 3
218 c0[2] = -5.; 224 c0[2] = -5.;
219 c1[2] = 2.5; 225 c1[2] = 2.5;
220 c2[2] = -0.1; 226 c2[2] = -0.1;
221 z1[2] = 60*mm; 227 z1[2] = 60*mm;
222 z2[2] = 130*mm; 228 z2[2] = 130*mm;
223 a0[2] = 10*mm; 229 a0[2] = 10*mm;
224 gradient[2] = 8.505263 *tesla/m; 230 gradient[2] = 8.505263 *tesla/m;
225 231
226 // QUADRUPOLE 4 232 // QUADRUPOLE 4
227 c0[3] = -5.; 233 c0[3] = -5.;
228 c1[3] = 2.5; 234 c1[3] = 2.5;
229 c2[3] = -0.1; 235 c2[3] = -0.1;
230 z1[3] = 60*mm; 236 z1[3] = 60*mm;
231 z2[3] = 130*mm; 237 z2[3] = 130*mm;
232 a0[3] = 10*mm; 238 a0[3] = 10*mm;
233 gradient[3] = -3.406526*tesla/m; 239 gradient[3] = -3.406526*tesla/m;
234 240
235 // FIELD CREATED BY A QUADRUPOLE IN ITS LOCA 241 // FIELD CREATED BY A QUADRUPOLE IN ITS LOCAL FRAME
236 G4double Bx_local,By_local,Bz_local; 242 G4double Bx_local,By_local,Bz_local;
237 Bx_local = 0; By_local = 0; Bz_local = 0; 243 Bx_local = 0; By_local = 0; Bz_local = 0;
238 244
239 // FIELD CREATED BY A QUADRUPOOLE IN WORLD F 245 // FIELD CREATED BY A QUADRUPOOLE IN WORLD FRAME
240 G4double Bx_quad,By_quad,Bz_quad; 246 G4double Bx_quad,By_quad,Bz_quad;
241 Bx_quad = 0; By_quad=0; Bz_quad=0; 247 Bx_quad = 0; By_quad=0; Bz_quad=0;
242 248
243 // QUADRUPOLE FRAME 249 // QUADRUPOLE FRAME
244 G4double x_local,y_local,z_local; 250 G4double x_local,y_local,z_local;
245 x_local= 0; y_local=0; z_local=0; 251 x_local= 0; y_local=0; z_local=0;
246 252
247 G4double vars = 0; 253 G4double vars = 0;
248 G4double G0, G1, G2, G3; 254 G4double G0, G1, G2, G3;
249 G4double K1, K2, K3; 255 G4double K1, K2, K3;
250 G4double P0, P1, P2, cte; 256 G4double P0, P1, P2, cte;
251 257
252 K1=0; 258 K1=0;
253 K2=0; 259 K2=0;
254 K3=0; 260 K3=0;
255 P0=0; 261 P0=0;
256 P1=0; 262 P1=0;
257 P2=0; 263 P2=0;
258 G0=0; 264 G0=0;
259 G1=0; 265 G1=0;
260 G2=0; 266 G2=0;
261 G3=0; 267 G3=0;
262 cte=0; 268 cte=0;
263 269
264 G4bool largeScattering=false; 270 G4bool largeScattering=false;
265 271
266 for (G4int i=0;i<4; i++) 272 for (G4int i=0;i<4; i++)
267 { 273 {
268 274
269 if (i==0) 275 if (i==0)
270 { xoprime = lineX + quadHalfLength*std::si 276 { xoprime = lineX + quadHalfLength*std::sin(lineAngle);
271 zoprime = lineZ + quadHalfLength*std::co 277 zoprime = lineZ + quadHalfLength*std::cos(lineAngle);
272 278
273 x_local = (x - xoprime) * std::cos (line 279 x_local = (x - xoprime) * std::cos (lineAngle) - (z - zoprime) * std::sin (lineAngle);
274 y_local = y; 280 y_local = y;
275 z_local = (z - zoprime) * std::cos (line 281 z_local = (z - zoprime) * std::cos (lineAngle) + (x - xoprime) * std::sin (lineAngle);
276 if (std::sqrt(x_local*x_local+y_local*y_ 282 if (std::sqrt(x_local*x_local+y_local*y_local)>a0[i]) largeScattering=true;
277 283
278 } 284 }
279 285
280 if (i==1) 286 if (i==1)
281 { xoprime = lineX + (3*quadHalfLength+quad 287 { xoprime = lineX + (3*quadHalfLength+quadSpacing)*std::sin(lineAngle);
282 zoprime = lineZ + (3*quadHalfLength+quad 288 zoprime = lineZ + (3*quadHalfLength+quadSpacing)*std::cos(lineAngle);
283 289
284 x_local = (x - xoprime) * std::cos (line 290 x_local = (x - xoprime) * std::cos (lineAngle) - (z - zoprime) * std::sin (lineAngle);
285 y_local = y; 291 y_local = y;
286 z_local = (z - zoprime) * std::cos (line 292 z_local = (z - zoprime) * std::cos (lineAngle) + (x - xoprime) * std::sin (lineAngle);
287 if (std::sqrt(x_local*x_local+y_local*y_ 293 if (std::sqrt(x_local*x_local+y_local*y_local)>a0[i]) largeScattering=true;
288 } 294 }
289 295
290 if (i==2) 296 if (i==2)
291 { xoprime = lineX + (5*quadHalfLength+2*qu 297 { xoprime = lineX + (5*quadHalfLength+2*quadSpacing)*std::sin(lineAngle);
292 zoprime = lineZ + (5*quadHalfLength+2*qu 298 zoprime = lineZ + (5*quadHalfLength+2*quadSpacing)*std::cos(lineAngle);
293 299
294 x_local = (x - xoprime) * std::cos (line 300 x_local = (x - xoprime) * std::cos (lineAngle) - (z - zoprime) * std::sin (lineAngle);
295 y_local = y; 301 y_local = y;
296 z_local = (z - zoprime) * std::cos (line 302 z_local = (z - zoprime) * std::cos (lineAngle) + (x - xoprime) * std::sin (lineAngle);
297 if (std::sqrt(x_local*x_local+y_local*y_ 303 if (std::sqrt(x_local*x_local+y_local*y_local)>a0[i]) largeScattering=true;
298 } 304 }
299 305
300 if (i==3) 306 if (i==3)
301 { xoprime = lineX + (7*quadHalfLength+3*qu 307 { xoprime = lineX + (7*quadHalfLength+3*quadSpacing)*std::sin(lineAngle);
302 zoprime = lineZ + (7*quadHalfLength+3*qu 308 zoprime = lineZ + (7*quadHalfLength+3*quadSpacing)*std::cos(lineAngle);
303 309
304 x_local = (x - xoprime) * std::cos (line 310 x_local = (x - xoprime) * std::cos (lineAngle) - (z - zoprime) * std::sin (lineAngle);
305 y_local = y; 311 y_local = y;
306 z_local = (z - zoprime) * std::cos (line 312 z_local = (z - zoprime) * std::cos (lineAngle) + (x - xoprime) * std::sin (lineAngle);
307 if (std::sqrt(x_local*x_local+y_local*y_ 313 if (std::sqrt(x_local*x_local+y_local*y_local)>a0[i]) largeScattering=true;
308 } 314 }
309 315
310 316
311 if ( z_local < -z2[i] ) 317 if ( z_local < -z2[i] )
312 { 318 {
313 G0=0; 319 G0=0;
314 G1=0; 320 G1=0;
315 G2=0; 321 G2=0;
316 G3=0; 322 G3=0;
317 } 323 }
318 324
319 if ( z_local > z2[i] ) 325 if ( z_local > z2[i] )
320 { 326 {
321 G0=0; 327 G0=0;
322 G1=0; 328 G1=0;
323 G2=0; 329 G2=0;
324 G3=0; 330 G3=0;
325 } 331 }
326 332
327 if ( (z_local>=-z1[i]) & (z_local<=z1[i]) ) 333 if ( (z_local>=-z1[i]) & (z_local<=z1[i]) )
328 { 334 {
329 G0=gradient[i]; 335 G0=gradient[i];
330 G1=0; 336 G1=0;
331 G2=0; 337 G2=0;
332 G3=0; 338 G3=0;
333 } 339 }
334 340
335 if ( ((z_local>=-z2[i]) & (z_local<-z1[i])) 341 if ( ((z_local>=-z2[i]) & (z_local<-z1[i])) || ((z_local>z1[i]) & (z_local<=z2[i])) )
336 { 342 {
337 343
338 vars = ( z_local - z1[i]) / a0[i] ; 344 vars = ( z_local - z1[i]) / a0[i] ;
339 if (z_local<-z1[i]) vars = ( - z_local - 345 if (z_local<-z1[i]) vars = ( - z_local - z1[i]) / a0[i] ;
340 346
341 347
342 P0 = c0[i]+c1[i]*vars+c2[i]*vars*vars; 348 P0 = c0[i]+c1[i]*vars+c2[i]*vars*vars;
343 349
344 P1 = c1[i]/a0[i]+2*c2[i]*(z_local-z1[i])/a 350 P1 = c1[i]/a0[i]+2*c2[i]*(z_local-z1[i])/a0[i]/a0[i];
345 if (z_local<-z1[i]) P1 = -c1[i]/a0[i]+2*c 351 if (z_local<-z1[i]) P1 = -c1[i]/a0[i]+2*c2[i]*(z_local+z1[i])/a0[i]/a0[i];
346 352
347 P2 = 2*c2[i]/a0[i]/a0[i]; 353 P2 = 2*c2[i]/a0[i]/a0[i];
348 354
349 cte = 1 + G4Exp(c0[i]); 355 cte = 1 + G4Exp(c0[i]);
350 356
351 K1 = -cte*P1*G4Exp(P0)/( (1+G4Exp(P0))*(1+ 357 K1 = -cte*P1*G4Exp(P0)/( (1+G4Exp(P0))*(1+G4Exp(P0)) );
352 358
353 K2 = -cte*G4Exp(P0)*( 359 K2 = -cte*G4Exp(P0)*(
354 P2/( (1+G4Exp(P0))*(1+G4Exp(P0)) ) 360 P2/( (1+G4Exp(P0))*(1+G4Exp(P0)) )
355 +2*P1*K1/(1+G4Exp(P0))/cte 361 +2*P1*K1/(1+G4Exp(P0))/cte
356 +P1*P1/(1+G4Exp(P0))/(1+G4Exp(P0)) 362 +P1*P1/(1+G4Exp(P0))/(1+G4Exp(P0))
357 ); 363 );
358 364
359 K3 = -cte*G4Exp(P0)*( 365 K3 = -cte*G4Exp(P0)*(
360 (3*P2*P1+P1*P1*P1)/(1+G4Exp(P0))/(1+G4Exp( 366 (3*P2*P1+P1*P1*P1)/(1+G4Exp(P0))/(1+G4Exp(P0))
361 +4*K1*(P1*P1+P2)/(1+G4Exp(P0))/cte 367 +4*K1*(P1*P1+P2)/(1+G4Exp(P0))/cte
362 +2*P1*(K1*K1/cte/cte+K2/(1+G4Exp(P0))/cte) 368 +2*P1*(K1*K1/cte/cte+K2/(1+G4Exp(P0))/cte)
363 ); 369 );
364 370
365 G0 = gradient[i]*cte/(1+G4Exp(P0)); 371 G0 = gradient[i]*cte/(1+G4Exp(P0));
366 G1 = gradient[i]*K1; 372 G1 = gradient[i]*K1;
367 G2 = gradient[i]*K2; 373 G2 = gradient[i]*K2;
368 G3 = gradient[i]*K3; 374 G3 = gradient[i]*K3;
369 375
370 } 376 }
371 377
372 // PROTECTION AGAINST LARGE SCATTERING 378 // PROTECTION AGAINST LARGE SCATTERING
373 379
374 if ( largeScattering ) 380 if ( largeScattering )
375 { 381 {
376 G0=0; 382 G0=0;
377 G1=0; 383 G1=0;
378 G2=0; 384 G2=0;
379 G3=0; 385 G3=0;
380 } 386 }
381 387
382 // MAGNETIC FIELD COMPUTATION FOR EACH QUAD 388 // MAGNETIC FIELD COMPUTATION FOR EACH QUADRUPOLE
383 389
384 Bx_local = y_local*(G0-(1./12)*(3*x_local*x 390 Bx_local = y_local*(G0-(1./12)*(3*x_local*x_local+y_local*y_local)*G2);
385 By_local = x_local*(G0-(1./12)*(3*y_local*y 391 By_local = x_local*(G0-(1./12)*(3*y_local*y_local+x_local*x_local)*G2);
386 Bz_local = x_local*y_local*(G1-(1./12)*(x_l 392 Bz_local = x_local*y_local*(G1-(1./12)*(x_local*x_local+y_local*y_local)*G3);
387 393
388 Bx_quad = Bz_local*std::sin(lineAngle)+Bx_l 394 Bx_quad = Bz_local*std::sin(lineAngle)+Bx_local*std::cos(lineAngle);
389 By_quad = By_local; 395 By_quad = By_local;
390 Bz_quad = Bz_local*std::cos(lineAngle)-Bx_l 396 Bz_quad = Bz_local*std::cos(lineAngle)-Bx_local*std::sin(lineAngle);
391 397
392 // TOTAL MAGNETIC FIELD 398 // TOTAL MAGNETIC FIELD
393 399
394 Bx = Bx + Bx_quad ; 400 Bx = Bx + Bx_quad ;
395 By = By + By_quad ; 401 By = By + By_quad ;
396 Bz = Bz + Bz_quad ; 402 Bz = Bz + Bz_quad ;
397 403
398 } // LOOP ON QUADRUPOLES 404 } // LOOP ON QUADRUPOLES
399 405
400 406
401 } // END OF QUADRUPLET 407 } // END OF QUADRUPLET
402 408
403 Bfield[0] = Bx; 409 Bfield[0] = Bx;
404 Bfield[1] = By; 410 Bfield[1] = By;
405 Bfield[2] = Bz; 411 Bfield[2] = Bz;
406 412
407 // ***************************************** 413 // *****************************************
408 // ELECTRIC FIELD CREATED BY SCANNING PLATES 414 // ELECTRIC FIELD CREATED BY SCANNING PLATES
409 // ***************************************** 415 // *****************************************
410 416
411 Bfield[3] = 0; 417 Bfield[3] = 0;
412 Bfield[4] = 0; 418 Bfield[4] = 0;
413 Bfield[5] = 0; 419 Bfield[5] = 0;
414 420
415 // POSITION OF EXIT OF LAST QUAD WHERE THE S 421 // POSITION OF EXIT OF LAST QUAD WHERE THE SCANNING PLATES START
416 422
417 G4double electricPlateWidth1 = 5 * mm; 423 G4double electricPlateWidth1 = 5 * mm;
418 G4double electricPlateWidth2 = 5 * mm; 424 G4double electricPlateWidth2 = 5 * mm;
419 G4double electricPlateLength1 = 36 * mm; 425 G4double electricPlateLength1 = 36 * mm;
420 G4double electricPlateLength2 = 34 * mm; 426 G4double electricPlateLength2 = 34 * mm;
421 G4double electricPlateGap = 5 * mm; 427 G4double electricPlateGap = 5 * mm;
422 G4double electricPlateSpacing1 = 3 * mm; 428 G4double electricPlateSpacing1 = 3 * mm;
423 G4double electricPlateSpacing2 = 4 * mm; 429 G4double electricPlateSpacing2 = 4 * mm;
424 430
425 // APPLY VOLTAGE HERE IN VOLTS (no electrost 431 // APPLY VOLTAGE HERE IN VOLTS (no electrostatic deflection here)
426 G4double electricPlateVoltage1 = 0 * volt; 432 G4double electricPlateVoltage1 = 0 * volt;
427 G4double electricPlateVoltage2 = 0 * volt; 433 G4double electricPlateVoltage2 = 0 * volt;
428 434
429 G4double electricFieldPlate1 = electricPlate 435 G4double electricFieldPlate1 = electricPlateVoltage1 / electricPlateSpacing1 ;
430 G4double electricFieldPlate2 = electricPlate 436 G4double electricFieldPlate2 = electricPlateVoltage2 / electricPlateSpacing2 ;
431 437
432 G4double beginFirstZoneX = lineX + (8*quadH 438 G4double beginFirstZoneX = lineX + (8*quadHalfLength+3*quadSpacing)*std::sin(lineAngle);
433 G4double beginFirstZoneZ = lineZ + (8*quadH 439 G4double beginFirstZoneZ = lineZ + (8*quadHalfLength+3*quadSpacing)*std::cos(lineAngle);
434 440
435 G4double beginSecondZoneX = lineX + (8*quad 441 G4double beginSecondZoneX = lineX + (8*quadHalfLength+3*quadSpacing+electricPlateLength1+electricPlateGap)*std::sin(lineAngle);
436 G4double beginSecondZoneZ = lineZ + (8*quad 442 G4double beginSecondZoneZ = lineZ + (8*quadHalfLength+3*quadSpacing+electricPlateLength1+electricPlateGap)*std::cos(lineAngle);
437 443
438 G4double xA, zA, xB, zB, xC, zC, xD, zD; 444 G4double xA, zA, xB, zB, xC, zC, xD, zD;
439 G4double slope1, cte1, slope2, cte2, slope3, 445 G4double slope1, cte1, slope2, cte2, slope3, cte3, slope4, cte4;
440 446
441 // WARNING : lineAngle < 0 447 // WARNING : lineAngle < 0
442 448
443 // FIRST PLATES 449 // FIRST PLATES
444 450
445 xA = beginFirstZoneX + std::cos(lineAngle)*e 451 xA = beginFirstZoneX + std::cos(lineAngle)*electricPlateSpacing1/2;
446 zA = beginFirstZoneZ - std::sin(lineAngle)*e 452 zA = beginFirstZoneZ - std::sin(lineAngle)*electricPlateSpacing1/2;
447 453
448 xB = xA + std::sin(lineAngle)*electricPlateL 454 xB = xA + std::sin(lineAngle)*electricPlateLength1;
449 zB = zA + std::cos(lineAngle)*electricPlateL 455 zB = zA + std::cos(lineAngle)*electricPlateLength1;
450 456
451 xC = xB - std::cos(lineAngle)*electricPlateS 457 xC = xB - std::cos(lineAngle)*electricPlateSpacing1;
452 zC = zB + std::sin(lineAngle)*electricPlateS 458 zC = zB + std::sin(lineAngle)*electricPlateSpacing1;
453 459
454 xD = xC - std::sin(lineAngle)*electricPlateL 460 xD = xC - std::sin(lineAngle)*electricPlateLength1;
455 zD = zC - std::cos(lineAngle)*electricPlateL 461 zD = zC - std::cos(lineAngle)*electricPlateLength1;
456 462
457 slope1 = (xB-xA)/(zB-zA); 463 slope1 = (xB-xA)/(zB-zA);
458 cte1 = xA - slope1 * zA; 464 cte1 = xA - slope1 * zA;
459 465
460 slope2 = (xC-xB)/(zC-zB); 466 slope2 = (xC-xB)/(zC-zB);
461 cte2 = xB - slope2 * zB; 467 cte2 = xB - slope2 * zB;
462 468
463 slope3 = (xD-xC)/(zD-zC); 469 slope3 = (xD-xC)/(zD-zC);
464 cte3 = xC - slope3 * zC; 470 cte3 = xC - slope3 * zC;
465 471
466 slope4 = (xA-xD)/(zA-zD); 472 slope4 = (xA-xD)/(zA-zD);
467 cte4 = xD - slope4 * zD; 473 cte4 = xD - slope4 * zD;
468 474
469 475
470 if 476 if
471 ( 477 (
472 x <= slope1 * z + cte1 478 x <= slope1 * z + cte1
473 && x >= slope3 * z + cte3 479 && x >= slope3 * z + cte3
474 && x <= slope4 * z + cte4 480 && x <= slope4 * z + cte4
475 && x >= slope2 * z + cte2 481 && x >= slope2 * z + cte2
476 && std::abs(y)<=electricPlateWidth1/2 482 && std::abs(y)<=electricPlateWidth1/2
477 ) 483 )
478 484
479 { 485 {
480 Bfield[3] = electricFieldPlate1*std::cos 486 Bfield[3] = electricFieldPlate1*std::cos(lineAngle);
481 Bfield[4] = 0; 487 Bfield[4] = 0;
482 Bfield[5] = -electricFieldPlate1*std::si 488 Bfield[5] = -electricFieldPlate1*std::sin(lineAngle);
483 489
484 } 490 }
485 491
486 // SECOND PLATES 492 // SECOND PLATES
487 493
488 xA = beginSecondZoneX + std::cos(lineAngle)* 494 xA = beginSecondZoneX + std::cos(lineAngle)*electricPlateWidth2/2;
489 zA = beginSecondZoneZ - std::sin(lineAngle)* 495 zA = beginSecondZoneZ - std::sin(lineAngle)*electricPlateWidth2/2;
490 496
491 xB = xA + std::sin(lineAngle)*electricPlateL 497 xB = xA + std::sin(lineAngle)*electricPlateLength2;
492 zB = zA + std::cos(lineAngle)*electricPlateL 498 zB = zA + std::cos(lineAngle)*electricPlateLength2;
493 499
494 xC = xB - std::cos(lineAngle)*electricPlateW 500 xC = xB - std::cos(lineAngle)*electricPlateWidth2;
495 zC = zB + std::sin(lineAngle)*electricPlateW 501 zC = zB + std::sin(lineAngle)*electricPlateWidth2;
496 502
497 xD = xC - std::sin(lineAngle)*electricPlateL 503 xD = xC - std::sin(lineAngle)*electricPlateLength2;
498 zD = zC - std::cos(lineAngle)*electricPlateL 504 zD = zC - std::cos(lineAngle)*electricPlateLength2;
499 505
500 slope1 = (xB-xA)/(zB-zA); 506 slope1 = (xB-xA)/(zB-zA);
501 cte1 = xA - slope1 * zA; 507 cte1 = xA - slope1 * zA;
502 508
503 slope2 = (xC-xB)/(zC-zB); 509 slope2 = (xC-xB)/(zC-zB);
504 cte2 = xB - slope2 * zB; 510 cte2 = xB - slope2 * zB;
505 511
506 slope3 = (xD-xC)/(zD-zC); 512 slope3 = (xD-xC)/(zD-zC);
507 cte3 = xC - slope3 * zC; 513 cte3 = xC - slope3 * zC;
508 514
509 slope4 = (xA-xD)/(zA-zD); 515 slope4 = (xA-xD)/(zA-zD);
510 cte4 = xD - slope4 * zD; 516 cte4 = xD - slope4 * zD;
511 517
512 if 518 if
513 ( 519 (
514 x <= slope1 * z + cte1 520 x <= slope1 * z + cte1
515 && x >= slope3 * z + cte3 521 && x >= slope3 * z + cte3
516 && x <= slope4 * z + cte4 522 && x <= slope4 * z + cte4
517 && x >= slope2 * z + cte2 523 && x >= slope2 * z + cte2
518 && std::abs(y)<=electricPlateSpacing2/2 524 && std::abs(y)<=electricPlateSpacing2/2
519 ) 525 )
520 526
521 { 527 {
522 Bfield[3] = 0; 528 Bfield[3] = 0;
523 Bfield[4] = electricFieldPlate2; 529 Bfield[4] = electricFieldPlate2;
524 Bfield[5] = 0; 530 Bfield[5] = 0;
525 } 531 }
>> 532
>> 533 //
526 534
527 } 535 }
528 536