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