Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-31 09:22:33 

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // This example is provided by the Geant4-DNA collaboration
0027 // Any report or published results obtained using the Geant4-DNA software 
0028 // shall cite the following Geant4-DNA collaboration publication:
0029 // Med. Phys. 37 (2010) 4692-4708
0030 // The Geant4-DNA web site is available at http://geant4-dna.org
0031 // 
0032 // If you use this example, please cite the following publication:
0033 // Rad. Prot. Dos. 133 (2009) 2-11
0034 //
0035 // Based on purging magnet advanced example.
0036 //
0037 
0038 #include "EMField.hh"
0039 #include "G4Exp.hh"
0040 #include "G4SystemOfUnits.hh"
0041 
0042 EMField::EMField() 
0043 {}
0044 
0045 void EMField::GetFieldValue(const double point[4], double *Bfield ) const
0046 { 
0047   // Magnetic field
0048   Bfield[0] = 0;
0049   Bfield[1] = 0;
0050   Bfield[2] = 0;
0051   
0052   // Electric field
0053   Bfield[3] = 0;
0054   Bfield[4] = 0;
0055   Bfield[5] = 0;
0056 
0057   G4double Bx = 0;
0058   G4double By = 0;
0059   G4double Bz = 0;
0060    
0061   G4double x = point[0];
0062   G4double y = point[1];
0063   G4double z = point[2];
0064 
0065 // ***********************
0066 // AIFIRA SWITCHING MAGNET
0067 // ***********************
0068   
0069   // MAGNETIC FIELD VALUE FOR 3 MeV ALPHAS
0070   G4double switchingField = 0.0589768635 * tesla ;
0071   
0072   // BEAM START
0073   G4double beamStart = -10*m;
0074 
0075   // RADIUS
0076   G4double Rp = 0.698*m;
0077 
0078   // ENTRANCE POSITION AFTER ANALYSIS MAGNET
0079   G4double zS = 975*mm;
0080   
0081   // POLE GAP
0082   G4double D = 31.8*mm;
0083   
0084   // FRINGING FIELD
0085 
0086   G4double fieldBoundary, wc0, wc1, wc2, wc3, limitMinEntrance, limitMaxEntrance, limitMinExit, limitMaxExit;
0087 
0088   limitMinEntrance = beamStart+zS-4*D;
0089   limitMaxEntrance = beamStart+zS+4*D;
0090   limitMinExit =Rp-4*D;
0091   limitMaxExit =Rp+4*D;  
0092     
0093   wc0 = 0.3835;
0094   wc1 = 2.388;
0095   wc2 = -0.8171;
0096   wc3 = 0.200;
0097 
0098   fieldBoundary=0.62;
0099 
0100   G4double ws, largeS, h, dhdlargeS, dhds, dlargeSds, dsdz, dsdx, zs0, Rs0, xcenter, zcenter;
0101   
0102 // - ENTRANCE OF SWITCHING MAGNET
0103 
0104 if ( (z >= limitMinEntrance) && (z < limitMaxEntrance) ) 
0105 {
0106   zs0 = fieldBoundary*D;
0107   ws = (-z+beamStart+zS-zs0)/D;
0108   dsdz = -1/D;
0109   dsdx = 0;
0110 
0111   largeS = wc0 + wc1*ws + wc2*ws*ws + wc3*ws*ws*ws;
0112   h = 1./(1.+G4Exp(largeS));
0113   dhdlargeS = -G4Exp(largeS)*h*h;  
0114   dlargeSds = wc1+ 2*wc2*ws + 3*wc3*ws*ws;
0115   dhds = dhdlargeS * dlargeSds;
0116       
0117   By = switchingField * h ;
0118   Bx = y*switchingField*dhds*dsdx;
0119   Bz = y*switchingField*dhds*dsdz;
0120 
0121 }
0122 
0123 // - HEART OF SWITCHING MAGNET    
0124         
0125  if ( 
0126           (z >= limitMaxEntrance)  
0127      &&   (( x*x + (z -(beamStart+zS))*(z -(beamStart+zS)) < limitMinExit*limitMinExit)) 
0128     )   
0129 {
0130    Bx=0; 
0131    By = switchingField; 
0132    Bz=0;
0133 }                               
0134     
0135 // - EXIT OF SWITCHING MAGNET
0136 
0137 if ( 
0138         (z >= limitMaxEntrance)  
0139      && (( x*x + (z -(beamStart+zS))*(z -(beamStart+zS))) >= limitMinExit*limitMinExit) 
0140      && (( x*x + (z -(beamStart+zS))*(z -(beamStart+zS))) < limitMaxExit*limitMaxExit)
0141 
0142    )    
0143 {
0144 
0145   xcenter = 0;
0146   zcenter =  beamStart+zS;
0147   
0148   Rs0 = Rp + D*fieldBoundary;
0149   ws = (std::sqrt((z-zcenter)*(z-zcenter)+(x-xcenter)*(x-xcenter)) - Rs0)/D;
0150     
0151   dsdz = (1/D)*(z-zcenter)/std::sqrt((z-zcenter)*(z-zcenter)+(x-xcenter)*(x-xcenter));
0152   dsdx = (1/D)*(x-xcenter)/std::sqrt((z-zcenter)*(z-zcenter)+(x-xcenter)*(x-xcenter));
0153 
0154   largeS = wc0 + wc1*ws + wc2*ws*ws + wc3*ws*ws*ws;
0155   h = 1./(1.+G4Exp(largeS));
0156   dhdlargeS = -G4Exp(largeS)*h*h;  
0157   dlargeSds = wc1+ 2*wc2*ws + 3*wc3*ws*ws;
0158   dhds = dhdlargeS * dlargeSds;
0159       
0160   By = switchingField * h ;
0161   Bx = y*switchingField*dhds*dsdx;
0162   Bz = y*switchingField*dhds*dsdz;
0163 
0164 }
0165 
0166 // **************************
0167 // MICROBEAM LINE QUADRUPOLES
0168 // **************************
0169  
0170   // MICROBEAM LINE ANGLE
0171   G4double lineAngle = -10*deg;
0172   
0173   // X POSITION OF FIRST QUADRUPOLE
0174   G4double lineX = -1295.59*mm;
0175 
0176   // Z POSITION OF FIRST QUADRUPOLE
0177   G4double lineZ = -1327*mm;
0178 
0179   // Adjust magnetic zone absolute position
0180   lineX = lineX + 5.24*micrometer*std::cos(-lineAngle); // 5.24 = 1.3 + 3.94 micrometer (cf. DetectorConstruction)
0181   lineZ = lineZ + 5.24*micrometer*std::sin(-lineAngle);
0182        
0183   // QUADRUPOLE HALF LENGTH
0184   G4double quadHalfLength = 75*mm;
0185   
0186   // QUADRUPOLE SPACING
0187   G4double quadSpacing = 40*mm;
0188   
0189   // QUADRUPOLE CENTER COORDINATES
0190   G4double xoprime, zoprime;
0191   
0192 if (z>=-1400*mm && z <-200*mm)
0193 {
0194   Bx=0; By=0; Bz=0;
0195   
0196   // FRINGING FILED CONSTANTS
0197   G4double c0[4], c1[4], c2[4], z1[4], z2[4], a0[4], gradient[4];
0198   
0199   // QUADRUPOLE 1
0200   c0[0] = -5.;
0201   c1[0] = 2.5;
0202   c2[0] = -0.1;
0203   z1[0] = 60*mm;
0204   z2[0] = 130*mm;
0205   a0[0] = 10*mm;
0206   gradient[0] = 3.406526 *tesla/m;
0207 
0208   // QUADRUPOLE 2
0209   c0[1] = -5.;
0210   c1[1] = 2.5;
0211   c2[1] = -0.1;
0212   z1[1] = 60*mm;
0213   z2[1] = 130*mm;
0214   a0[1] = 10*mm;
0215   gradient[1] = -8.505263 *tesla/m;
0216 
0217   // QUADRUPOLE 3
0218   c0[2] = -5.;
0219   c1[2] = 2.5;
0220   c2[2] = -0.1;
0221   z1[2] = 60*mm;
0222   z2[2] = 130*mm;
0223   a0[2] = 10*mm;
0224   gradient[2] = 8.505263 *tesla/m;
0225 
0226   // QUADRUPOLE 4
0227   c0[3] = -5.;
0228   c1[3] = 2.5;
0229   c2[3] = -0.1;
0230   z1[3] = 60*mm;
0231   z2[3] = 130*mm;
0232   a0[3] = 10*mm;
0233   gradient[3] = -3.406526*tesla/m;
0234 
0235   // FIELD CREATED BY A QUADRUPOLE IN ITS LOCAL FRAME
0236   G4double Bx_local,By_local,Bz_local;
0237   Bx_local = 0; By_local = 0; Bz_local = 0;
0238   
0239   // FIELD CREATED BY A QUADRUPOOLE IN WORLD FRAME
0240   G4double Bx_quad,By_quad,Bz_quad;
0241   Bx_quad = 0; By_quad=0; Bz_quad=0;
0242   
0243   // QUADRUPOLE FRAME
0244   G4double x_local,y_local,z_local;
0245   x_local= 0; y_local=0; z_local=0;
0246 
0247   G4double vars = 0;
0248   G4double G0, G1, G2, G3;
0249   G4double K1, K2, K3;
0250   G4double P0, P1, P2,     cte;
0251 
0252   K1=0;
0253   K2=0;
0254   K3=0;
0255   P0=0;
0256   P1=0;
0257   P2=0;
0258   G0=0;
0259   G1=0;
0260   G2=0;
0261   G3=0;
0262   cte=0;
0263 
0264   G4bool largeScattering=false;
0265   
0266   for (G4int i=0;i<4; i++) 
0267   {
0268  
0269      if (i==0) 
0270         {   xoprime = lineX + quadHalfLength*std::sin(lineAngle);
0271             zoprime = lineZ + quadHalfLength*std::cos(lineAngle);
0272 
0273             x_local = (x - xoprime) * std::cos (lineAngle) - (z - zoprime) * std::sin (lineAngle); 
0274             y_local = y; 
0275             z_local = (z - zoprime) * std::cos (lineAngle) + (x - xoprime) * std::sin (lineAngle); 
0276             if (std::sqrt(x_local*x_local+y_local*y_local)>a0[i]) largeScattering=true;
0277 
0278         }
0279          
0280      if (i==1) 
0281         {   xoprime = lineX + (3*quadHalfLength+quadSpacing)*std::sin(lineAngle);
0282             zoprime = lineZ + (3*quadHalfLength+quadSpacing)*std::cos(lineAngle);
0283 
0284             x_local = (x - xoprime) * std::cos (lineAngle) - (z - zoprime) * std::sin (lineAngle); 
0285             y_local = y; 
0286             z_local = (z - zoprime) * std::cos (lineAngle) + (x - xoprime) * std::sin (lineAngle); 
0287             if (std::sqrt(x_local*x_local+y_local*y_local)>a0[i]) largeScattering=true;
0288         }
0289 
0290      if (i==2) 
0291         {   xoprime = lineX + (5*quadHalfLength+2*quadSpacing)*std::sin(lineAngle);
0292             zoprime = lineZ + (5*quadHalfLength+2*quadSpacing)*std::cos(lineAngle);
0293 
0294             x_local = (x - xoprime) * std::cos (lineAngle) - (z - zoprime) * std::sin (lineAngle); 
0295             y_local = y; 
0296             z_local = (z - zoprime) * std::cos (lineAngle) + (x - xoprime) * std::sin (lineAngle); 
0297             if (std::sqrt(x_local*x_local+y_local*y_local)>a0[i]) largeScattering=true;
0298         }
0299      
0300      if (i==3) 
0301         {   xoprime = lineX + (7*quadHalfLength+3*quadSpacing)*std::sin(lineAngle);
0302             zoprime = lineZ + (7*quadHalfLength+3*quadSpacing)*std::cos(lineAngle);
0303 
0304             x_local = (x - xoprime) * std::cos (lineAngle) - (z - zoprime) * std::sin (lineAngle); 
0305             y_local = y; 
0306             z_local = (z - zoprime) * std::cos (lineAngle) + (x - xoprime) * std::sin (lineAngle); 
0307             if (std::sqrt(x_local*x_local+y_local*y_local)>a0[i]) largeScattering=true;
0308         }
0309 
0310      
0311      if ( z_local < -z2[i] )
0312      {
0313       G0=0;
0314       G1=0;
0315       G2=0;
0316       G3=0;
0317      }
0318      
0319      if ( z_local > z2[i] )
0320      {
0321       G0=0;
0322       G1=0;
0323       G2=0;
0324       G3=0;
0325      }
0326 
0327      if ( (z_local>=-z1[i]) & (z_local<=z1[i]) ) 
0328      {
0329       G0=gradient[i];
0330       G1=0;
0331       G2=0;
0332       G3=0;
0333      }
0334      
0335      if ( ((z_local>=-z2[i]) & (z_local<-z1[i])) ||  ((z_local>z1[i]) & (z_local<=z2[i])) ) 
0336      {
0337 
0338       vars = ( z_local - z1[i]) / a0[i] ;
0339       if (z_local<-z1[i]) vars = ( - z_local - z1[i]) / a0[i] ;
0340 
0341 
0342       P0 = c0[i]+c1[i]*vars+c2[i]*vars*vars;
0343 
0344       P1 = c1[i]/a0[i]+2*c2[i]*(z_local-z1[i])/a0[i]/a0[i];
0345       if (z_local<-z1[i])  P1 = -c1[i]/a0[i]+2*c2[i]*(z_local+z1[i])/a0[i]/a0[i];
0346 
0347       P2 = 2*c2[i]/a0[i]/a0[i];
0348 
0349       cte = 1 + G4Exp(c0[i]);
0350 
0351       K1 = -cte*P1*G4Exp(P0)/( (1+G4Exp(P0))*(1+G4Exp(P0)) );
0352 
0353       K2 = -cte*G4Exp(P0)*(
0354        P2/( (1+G4Exp(P0))*(1+G4Exp(P0)) )
0355       +2*P1*K1/(1+G4Exp(P0))/cte
0356       +P1*P1/(1+G4Exp(P0))/(1+G4Exp(P0))
0357       );
0358  
0359       K3 = -cte*G4Exp(P0)*(
0360       (3*P2*P1+P1*P1*P1)/(1+G4Exp(P0))/(1+G4Exp(P0))
0361       +4*K1*(P1*P1+P2)/(1+G4Exp(P0))/cte
0362       +2*P1*(K1*K1/cte/cte+K2/(1+G4Exp(P0))/cte)
0363        );
0364       
0365       G0 = gradient[i]*cte/(1+G4Exp(P0));
0366       G1 = gradient[i]*K1;
0367       G2 = gradient[i]*K2;
0368       G3 = gradient[i]*K3;
0369 
0370      }
0371       
0372      // PROTECTION AGAINST LARGE SCATTERING
0373 
0374      if ( largeScattering ) 
0375      {
0376       G0=0;
0377       G1=0;
0378       G2=0;
0379       G3=0;
0380      }
0381 
0382      // MAGNETIC FIELD COMPUTATION FOR EACH QUADRUPOLE
0383      
0384      Bx_local = y_local*(G0-(1./12)*(3*x_local*x_local+y_local*y_local)*G2);
0385      By_local = x_local*(G0-(1./12)*(3*y_local*y_local+x_local*x_local)*G2);
0386      Bz_local = x_local*y_local*(G1-(1./12)*(x_local*x_local+y_local*y_local)*G3);
0387 
0388      Bx_quad = Bz_local*std::sin(lineAngle)+Bx_local*std::cos(lineAngle);
0389      By_quad = By_local;
0390      Bz_quad = Bz_local*std::cos(lineAngle)-Bx_local*std::sin(lineAngle);
0391 
0392      // TOTAL MAGNETIC FIELD
0393      
0394      Bx = Bx + Bx_quad ;
0395      By = By + By_quad ;
0396      Bz = Bz + Bz_quad ;
0397 
0398   } // LOOP ON QUADRUPOLES
0399 
0400       
0401 } // END OF QUADRUPLET
0402 
0403   Bfield[0] = Bx;
0404   Bfield[1] = By;
0405   Bfield[2] = Bz;
0406 
0407 // *****************************************
0408 // ELECTRIC FIELD CREATED BY SCANNING PLATES
0409 // *****************************************
0410 
0411   Bfield[3] = 0;
0412   Bfield[4] = 0;
0413   Bfield[5] = 0;
0414 
0415   // POSITION OF EXIT OF LAST QUAD WHERE THE SCANNING PLATES START
0416 
0417   G4double electricPlateWidth1 = 5 * mm;
0418   G4double electricPlateWidth2 = 5 * mm;
0419   G4double electricPlateLength1 = 36 * mm;
0420   G4double electricPlateLength2 = 34 * mm;
0421   G4double electricPlateGap = 5 * mm;
0422   G4double electricPlateSpacing1 = 3 * mm;
0423   G4double electricPlateSpacing2 = 4 * mm;
0424 
0425   // APPLY VOLTAGE HERE IN VOLTS (no electrostatic deflection here)
0426   G4double electricPlateVoltage1 = 0 * volt;
0427   G4double electricPlateVoltage2 = 0 * volt;
0428 
0429   G4double electricFieldPlate1 = electricPlateVoltage1 / electricPlateSpacing1 ;
0430   G4double electricFieldPlate2 = electricPlateVoltage2 / electricPlateSpacing2 ;
0431 
0432   G4double  beginFirstZoneX = lineX + (8*quadHalfLength+3*quadSpacing)*std::sin(lineAngle);
0433   G4double  beginFirstZoneZ = lineZ + (8*quadHalfLength+3*quadSpacing)*std::cos(lineAngle);
0434 
0435   G4double  beginSecondZoneX = lineX + (8*quadHalfLength+3*quadSpacing+electricPlateLength1+electricPlateGap)*std::sin(lineAngle);
0436   G4double  beginSecondZoneZ = lineZ + (8*quadHalfLength+3*quadSpacing+electricPlateLength1+electricPlateGap)*std::cos(lineAngle);
0437 
0438   G4double xA, zA, xB, zB, xC, zC, xD, zD;
0439   G4double slope1, cte1, slope2, cte2, slope3, cte3, slope4, cte4;
0440  
0441   // WARNING : lineAngle < 0
0442 
0443   // FIRST PLATES
0444   
0445   xA = beginFirstZoneX + std::cos(lineAngle)*electricPlateSpacing1/2;
0446   zA = beginFirstZoneZ - std::sin(lineAngle)*electricPlateSpacing1/2;
0447 
0448   xB = xA + std::sin(lineAngle)*electricPlateLength1; 
0449   zB = zA + std::cos(lineAngle)*electricPlateLength1;
0450   
0451   xC = xB - std::cos(lineAngle)*electricPlateSpacing1;
0452   zC = zB + std::sin(lineAngle)*electricPlateSpacing1;
0453 
0454   xD = xC - std::sin(lineAngle)*electricPlateLength1; 
0455   zD = zC - std::cos(lineAngle)*electricPlateLength1;
0456   
0457   slope1 = (xB-xA)/(zB-zA);
0458   cte1 = xA - slope1 * zA;
0459   
0460   slope2 = (xC-xB)/(zC-zB);
0461   cte2 = xB - slope2 * zB;
0462   
0463   slope3 = (xD-xC)/(zD-zC);
0464   cte3 = xC - slope3 * zC;
0465   
0466   slope4 = (xA-xD)/(zA-zD);
0467   cte4 = xD - slope4 * zD;
0468   
0469    
0470   if 
0471   (
0472        x <= slope1 * z + cte1
0473     && x >= slope3 * z + cte3
0474     && x <= slope4 * z + cte4
0475     && x >= slope2 * z + cte2    
0476     && std::abs(y)<=electricPlateWidth1/2
0477   )  
0478 
0479   {
0480       Bfield[3] = electricFieldPlate1*std::cos(lineAngle);
0481       Bfield[4] = 0;
0482       Bfield[5] = -electricFieldPlate1*std::sin(lineAngle);
0483  
0484   }
0485       
0486   // SECOND PLATES
0487       
0488   xA = beginSecondZoneX + std::cos(lineAngle)*electricPlateWidth2/2;
0489   zA = beginSecondZoneZ - std::sin(lineAngle)*electricPlateWidth2/2;
0490 
0491   xB = xA + std::sin(lineAngle)*electricPlateLength2; 
0492   zB = zA + std::cos(lineAngle)*electricPlateLength2;
0493   
0494   xC = xB - std::cos(lineAngle)*electricPlateWidth2;
0495   zC = zB + std::sin(lineAngle)*electricPlateWidth2;
0496 
0497   xD = xC - std::sin(lineAngle)*electricPlateLength2; 
0498   zD = zC - std::cos(lineAngle)*electricPlateLength2;
0499   
0500   slope1 = (xB-xA)/(zB-zA);
0501   cte1 = xA - slope1 * zA;
0502   
0503   slope2 = (xC-xB)/(zC-zB);
0504   cte2 = xB - slope2 * zB;
0505   
0506   slope3 = (xD-xC)/(zD-zC);
0507   cte3 = xC - slope3 * zC;
0508   
0509   slope4 = (xA-xD)/(zA-zD);
0510   cte4 = xD - slope4 * zD;
0511 
0512   if 
0513   (     
0514        x <= slope1 * z + cte1
0515     && x >= slope3 * z + cte3
0516     && x <= slope4 * z + cte4
0517     && x >= slope2 * z + cte2    
0518     && std::abs(y)<=electricPlateSpacing2/2
0519   )
0520 
0521   {  
0522       Bfield[3] = 0;
0523       Bfield[4] = electricFieldPlate2;
0524       Bfield[5] = 0;
0525   }
0526 
0527 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!