HelicalDipoleField Class Reference

Inheritance diagram for HelicalDipoleField:

Detailed Description

HelicalDipoleField represents one placement of a helicaldipole magnet.


Public Member Functions
	HelicalDipoleField (BLCoordinateTransform &_global2local, BLCMDhelicaldipole *hd)
	constructor.
void	addFieldValue (const G4double point[4], G4double field[6]) const
	addFieldValue() adds the field for this solenoid into field[]. point[] is in global coordinates.
Private Attributes
G4double	radius
G4double	halflength
G4double	b
G4double	bprime
G4double	bpp
G4double	rr
G4double	psi0
G4double	lambda
G4double	phi0
G4double	Bsolenoid
G4double	model
G4double	ez
BLCoordinateTransform	global2local
G4RotationMatrix	rotation

Constructor & Destructor Documentation

HelicalDipoleField::HelicalDipoleField	(	BLCoordinateTransform &	_global2local,
		BLCMDhelicaldipole *	hd
	)

constructor.

00364                                                                 :
00365                                         BLElementField(), rotation()
00366 {
00367         radius = hd->radius;
00368         halflength = hd->length/2.0;
00369         b = hd->b;
00370         bprime = hd->bprime;
00371         bpp    = hd->bpp;
00372         rr     = hd->rr;
00373         model = hd->model;
00374         lambda = hd->lambda;
00375         phi0 = hd->phi0;
00376         psi0 = hd->psi0;
00377         Bsolenoid = hd->Bsolenoid;
00378         ez = hd->ez;
00379         global2local = _global2local;
00380         rotation = global2local.getRotation().inverse();
00381 
00382 //printf("HelicalDipoleField: radius=%.1f mm, halflength=%.1f mm, b=%.4f T, bprime=%.4f T/m\n\tmodel=%.1f, lambda=%.1f mm, phi0=%.1f deg, psi0=%.1f deg, Bsolenoid=%.4f T\n",radius/mm,halflength/mm,b/tesla,bprime/(tesla/meter),model,lambda/mm,phi0/deg,psi0/deg,Bsolenoid/tesla);
00383 
00384         // set global bounding box
00385         G4double local[4], global[4];
00386         local[3] = 0.0;
00387         for(int i=0; i<2; ++i) {
00388                 local[0] = (i==0 ? -1.0 : 1.0) * radius;
00389                 for(int j=0; j<2; ++j) {
00390                         local[1] = (j==0 ? -1.0 : 1.0) * radius;
00391                         for(int k=0; k<2; ++k) {
00392                                 local[2] = (k==0 ? -1.0 : 1.0) * halflength;
00393                                 global2local.getGlobal(local,global);
00394                                 setGlobalPoint(global);
00395                         }
00396                 }
00397         }
00398 }

Member Function Documentation

void HelicalDipoleField::addFieldValue	(	const G4double	point[4],
		G4double	field[6]
	)			const `[virtual]`

addFieldValue() adds the field for this solenoid into field[]. point[] is in global coordinates.

Implements BLElementField.

References b, B, bpp, bprime, Bsolenoid, cylPHI, cylRHO, CYLTOCARTESIAN(), DIPOLEFIELD(), ez, BLCoordinateTransform::getLocal(), global2local, halflength, ICOOLFIELD(), BLCoordinateTransform::isRotated(), KBB_bugz, lambda, model, MuonsInc_PITCH_CONVENTION, phi0, PI, psi0, QUADRUPOLEFIELD(), radius, rectX, rectY, rectZ, rotation, rr, and SIMPLEFIELD().

00597 {
00598   const char *CoDE="HelicalDipoleField::addFieldValue";
00599   const int X=rectX,Y=rectY,Z=rectZ;
00600   const int PHI=cylPHI,RHO=cylRHO;                      /* notation */
00601   const int node=0;
00602   int ini;
00603   G4ThreeVector B,BcylDipole,BcylQuad,BcylSext,BxyzDipole,BxyzQuad,BxyzSext;
00604   G4ThreeVector global(point[X],point[Y],point[Z]);
00605   G4ThreeVector local;
00606   G4double kz,phi,bd,rho,kH,Bprime,Bpp,pZ;
00607   G4double mul[3];
00608 
00609   global2local.getLocal(local,global);                  /*fetch local location in mm*/
00610   rho = sqrt( local[X]*local[X] + local[Y]*local[Y] );  /*mm*/
00611   
00612   if( rho>radius || fabs(local[Z]) > halflength)  return;
00613 
00614   for(ini=0;ini<=2;ini++){
00615     BcylDipole[ini]=0.0;
00616     BxyzDipole[ini]=0.0;
00617     BcylQuad[ini]=0.0;
00618     BxyzQuad[ini]=0.0;
00619     BcylSext[ini]=0.0;
00620     BxyzSext[ini]=0.0;
00621     B[ini]=0.0;
00622   }
00623   
00624   /* the Muons,Inc. preferred convention -> lambda>0 == right handed screw */
00625 
00626   kH= MuonsInc_PITCH_CONVENTION *2*PI/lambda;
00627   //  kH= sqrt(4.0*PI*PI/lambda/lambda);
00628   //  kz= kH * (local[Z]+halflength)*lambda/sqrt(lambda*lambda);
00629   kz= kH * (local[Z]+halflength);
00630   pZ = (local[Z]+halflength)*1e-3;
00631   phi= atan2( local[Y], local[X]);
00632 
00633   bd = b;   /* the radial correction for bQ, bprime, will now be done externally by the user */ 
00634 
00635   /* Simple SIN/COS Dipole Model */
00636   if( int(model)==1 )
00637     {
00638       B= SIMPLEFIELD(b,kz+phi0,Bsolenoid,model);
00639     }
00640   
00641   /* Dipole Maxwellian Model  */
00642   else if ( int(model)==2 )
00643     {   
00644       BcylDipole= DIPOLEFIELD( bd, kH, rho, phi-kz+phi0 );
00645       BxyzDipole = CYLTOCARTESIAN(BcylDipole,phi);
00646 
00647       BcylQuad= QUADRUPOLEFIELD( bprime, kH, rho, phi-kz+phi0+psi0 );
00648       BxyzQuad = CYLTOCARTESIAN(BcylQuad,phi);
00649 
00650       B[X]= BxyzDipole[X] + BxyzQuad[X];
00651       B[Y]= BxyzDipole[Y] + BxyzQuad[Y];
00652       B[Z]= BxyzDipole[Z] + BxyzQuad[Z] + Bsolenoid;
00653 
00654       if( KBB_bugz )
00655         {
00656           printf( "%s: BcylDipole=(%f,%f,%f) BcylQuad=(%f,%f,%f)\n", 
00657                   CoDE, 
00658                   BcylDipole[PHI]/tesla, BcylDipole[RHO]/tesla, BcylDipole[Z]/tesla,
00659                   BcylQuad[PHI]/tesla, BcylQuad[RHO]/tesla, BcylQuad[Z]/tesla );
00660           
00661           printf( "%s: Bxyz=(%f+%f=%f, %f+%f=%f, %f+%f+%f=%f)\n",
00662                   CoDE, 
00663                   BxyzDipole[X]/tesla, BxyzQuad[X]/tesla, B[X]/tesla,
00664                   BxyzDipole[Y]/tesla, BxyzQuad[Y]/tesla, B[Y]/tesla,
00665                   BxyzDipole[Z]/tesla, BxyzQuad[Z]/tesla, Bsolenoid, B[Z]/tesla );
00666         }
00667     }
00668 
00669   // Extended Dipole Maxwellian (ICOOL) Model
00670   else if ( int(model)==3 )
00671     {
00672       // TO DO: Extend to multiple components
00673       //        Use more realistic formula
00674       BcylDipole = bd*ICOOLFIELD(1,rho/mm,phi-kz-phi0,kH/mm,rr/mm);
00675       BxyzDipole = CYLTOCARTESIAN(BcylDipole,phi);
00676       Bprime= bprime*1e3;
00677       BcylQuad = Bprime*ICOOLFIELD(2,rho/mm,phi-kz-phi0+psi0,kH/mm,rr/mm);
00678       BxyzQuad   = CYLTOCARTESIAN(BcylQuad,phi);
00679       Bpp= bpp*1e6;
00680       BcylSext = Bpp*ICOOLFIELD(3,rho/mm,phi-kz-phi0+psi0,kH/mm,rr/mm);
00681       BxyzSext = CYLTOCARTESIAN(BcylSext,phi);
00682 
00683       B[PHI] =BcylDipole[PHI]+BcylQuad[PHI]+BcylSext[PHI];
00684       B[RHO] =BcylDipole[RHO]+BcylQuad[RHO]+BcylSext[RHO];
00685       B[Z]   =BcylDipole[Z]+BcylQuad[Z]+BcylSext[Z];
00686 
00687       B[X]= BxyzDipole[X] + BxyzQuad[X] + BxyzSext[X];  
00688       B[Y]= BxyzDipole[Y] + BxyzQuad[Y] + BxyzSext[Y];
00689       B[Z]= BxyzDipole[Z] + BxyzQuad[Z] + BxyzSext[Z] + Bsolenoid;
00690 
00691     }
00692 
00693   // Field map calculation for MANX experiment
00694   else if ( int(model)==4 )
00695     {
00696       // TO DO: Extend to multiple components
00697       //        Use more realistic formula
00698       // Add new terms to modulate bd,bq,bz (1/20/05 KY)
00699       // mul[0]=bd, mul[1]=bq, mul[2]=bz
00700 
00701       // Original values
00702       //mul[0] = 3.24178-0.369821*pZ+1.24318e-3*pZ*pZ-1.91452e-3*pZ*pZ*pZ;
00703       //mul[1] = -0.616436+6.69245e-2*pZ-5.36372e-4*pZ*pZ+3.23624e-4*pZ*pZ*pZ;
00704       //mul[2] = -11.5474+1.27729*pZ-1.03380e-2*pZ*pZ+6.48178e-3*pZ*pZ*pZ;
00705 
00706       // Values for 300 MeV/c input
00707       mul[2] = -8.53621+1.33517*pZ+3.51035e-2*pZ*pZ+6.48178e-3*pZ*pZ*pZ;
00708       mul[0] = 2.34731-0.409774*pZ-9.50945e-3*pZ*pZ-1.85969e-3*pZ*pZ*pZ;
00709       mul[1] = -0.442645+8.65446e-2*pZ+6.15396e-4*pZ*pZ+2.78646e-4*pZ*pZ*pZ;
00710 
00711       mul[0]*= bd;
00712       BcylDipole = mul[0]*ICOOLFIELD(1,rho/mm,phi-kz-phi0,kH/mm,159.155/mm);
00713       BxyzDipole = CYLTOCARTESIAN(BcylDipole,phi);
00714       mul[1]*= 1e3*bprime;
00715       BcylQuad = mul[1]*ICOOLFIELD(2,rho/mm,phi-kz-phi0+psi0,kH/mm,159.155/mm);
00716       BxyzQuad   = CYLTOCARTESIAN(BcylQuad,phi);
00717 //@@@ mul[3] is OUT OF BOUNDS! It is also UNINITIALIZED!
00718 //@@@ mul[3]*= bpp*1e6;
00719 //@@@ BcylSext = mul[3]*ICOOLFIELD(3,rho/mm,phi-kz-phi0+psi0,kH/mm,rr/mm);
00720 G4Exception("BLCMDhelicaldipole","Invalid model",FatalException, "Internal coding error");
00721 //@@@ BxyzSext = CYLTOCARTESIAN(BcylSext,phi);
00722       mul[2]*= Bsolenoid;
00723 
00724       //      printf("z=%e m1=%e m2=%e m3=%e\n",pZ,mul[0],mul[1],mul[2]);
00725 
00726       B[PHI] =BcylDipole[PHI]+BcylQuad[PHI]+BcylSext[PHI];
00727       B[RHO] =BcylDipole[RHO]+BcylQuad[RHO]+BcylSext[RHO];
00728       B[Z]   =BcylDipole[Z]+BcylQuad[Z]+BcylSext[Z];
00729 
00730       B[X]= BxyzDipole[X] + BxyzQuad[X] + BxyzSext[X];  
00731       B[Y]= BxyzDipole[Y] + BxyzQuad[Y] + BxyzSext[Y];
00732       B[Z]= BxyzDipole[Z] + BxyzQuad[Z] + BxyzSext[Z] + mul[2];
00733 
00734     }
00735   
00736   else
00737     {
00738       fprintf(stderr,"HelicalDipoleField::addFieldValue bad model#%d (1|2)\n", int(model) );
00739       fflush(stderr);
00740     }
00741   
00742   /* Rotation if applicable */
00743   
00744   if(global2local.isRotated())  B = rotation * B;
00745   
00746   field[0] += B[X];        /* update the field */
00747   field[1] += B[Y];
00748   field[2] += B[Z];
00749 
00750   field[5] += ez;
00751 
00752 }