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