MyBLFieldMap Class Reference

#include <MyBLFieldMap.hh>

Public Member Functions
	MyBLFieldMap ()
	default constructor.
virtual	~MyBLFieldMap ()
	destructor.
bool	readFile (G4String filename)
bool	writeFile (G4String filename, G4String comment="")
void	getFieldValue (const G4double local[4], G4double field[6], G4double current=1.0, G4double gradient=1.0)
void	getBoundingPoint (int i, G4double point[4])
	getBoundingPoint() returns the i-th bounding point of the map.
bool	hasB ()
	hasB() returns true if this map has a nonzero B field.
bool	hasE ()
	hasE() returns true if this map has a nonzero E field.
bool	createGridMap (G4double X0, G4double Y0, G4double Z0, G4double dX, G4double dY, G4double dZ, int nX, int nY, int nZ, class G4ElectroMagneticField *field)
bool	createCylinderMap (G4double Z0, G4double dR, G4double dZ, int nR, int nZ, class G4ElectroMagneticField *field)
bool	createTimeDependence (int n, G4double t[], G4double b[], G4double e[]=0, G4double period=-1.0)
bool	getTimeFactor (G4double t, G4double b, G4double e)
Private Attributes
G4int	maxline
G4double	current
G4double	gradient
G4double	normB
G4double	normE
class FieldMapImpl *	impl
class TimeImpl *	time
Friends
class	FieldMapPlacement

Detailed Description

class MyBLFieldMap implements a general field map, both B and E.

Initially this class simply reads an input file to define the map. Eventually it will also generate a map automatically from a BLCoil, and from an arbitrary collection of magnets (to improve tracking efficiency for a large number of overlapping solenoids).

Field components are interpolated in the map; the values used are: Bx = (interpolated value)*normB*(element_current/current_param) Ex = (interpolated value)*normE*(element_gradient/gradient_param) [other components are similar] where element_current and element_gradient are from the element definition command in the input file, and normB, normE, curren_param, and gradient_param come from the input file. The presence of both norm and current/gradient in the input file is to accommodate diverse sources of input files and units.

Outside the map the fields are zero, so if the map is truncated you should ensure that no particles are tracked outside the map in that region. This happens, for instance, for a beam solenoid truncated in radius to the inside of the coil -- particles which enter the coil will see a zero field, so the coil should be set to kill them; particles outside the map along z but inside the bore will also see a zero field, but that is usually OK as long as the map extends far enough along Z to include the nonzero field region.

Input File format: Blank lines, and lines beginning with # or * are comments. Lines beginning with * are printed to stdout. Units are mm for coordinates, Tesla for B, and MegaVolts/meter for E; use normB and normE if the data points use different units.

The input file starts with a set of commands to define the parameters of the map, followed by blocks of lines containing the values of the field components. The field component names depend on the type of map (grid: Bx,By,Bz,Ex,Ey,Ez; cylinder: Br,Bz,Er,Ez). Each command has a specific list of arguments to define parameters of the map.

BEWARE: the parsing is not exhaustive. For instance, invalid arguments are silently ignored (which means you must verify the spelling and capitalization of argument names). Correct inputs will yield correct results, but invalid inputs may not be detected and may yield seemingly-correct but unintended results.

The first command is usually a param command, which has the following arguments: maxline The maximum number of characters per line (default=1023) current The current corresponding to this map (default=1.0) gradient The gradient corresponding to the map (default=1.0) normE A normalization factor for E components (default=1.0) normB A normalization factor for B components (default=1.0)

Two types of maps are implemented: grid and cylinder.

grid maps are a 3-D grid, with each block of data being a single X-Y plane; within a block the lines are Y and the columns of each line are values along X. The grid command has the following arguments: X0 The X value for the first value in each line Y0 The Y value for the first line in each block Z0 The Z value for the first block of each field component nX The number of columns per line nY The number of lines per block nZ The number of blocks per field component dX The X increment between values in each line dY The Y increment between lines dZ The Z increment between blocks tolerance The tolerance for pointwise data (default=0.01 mm) After the grid command, the following optional commands can be given: extendX flip=... extendY flip=... extendZ flip=... These commands permit a half-map to be extended to the full map around X=0, Y=0, or Z=0 respectively. The optional flip argument is a comma-separated list of field components whose signs will be inverted for negative values of the coordinate. For example, "extendZ flip=Bx,Ex" means the map from Z=0 to Z=(nZ-1)*dZ is extended symmetrically around Z=0 to negative Z values, flipping the signs of Bx and Ex when Z<0. This could be followed by "extendX flip=Bx,Ex", and the field flips will be the products of both commands.

cylinder maps are a 2-D map with rotational symmetry around the Z axis. Each field component has a single block with lines being Z and the columns being R. The cylinder command has the following arguments: Z0 The Z value for the first line in each block nR The number of columns per line nZ The number of lines per block dR The R increment between colums dZ The Z increment between lines tolerance The tolerance for pointwise data (default=0.01 mm) After the cylinder command, the following optional commands can be given: extendZ flip=... This command behaves the same as for the grid map.

After the commands, each block consists of a line containing the name of the field component, followed by the lines of the block. The values within a line can be separated by whitespace or a ',' followed by optional whitespace. Field components that are not given are set to 0.0 everywhere. Missing values will be considered to be 0.0. For grid maps the first block is for Z=Z0, and successive blocks increment Z by dZ; the first line in a block is for Y=Y0 and the first column in each line is for X=X0. For cylinder maps, the first line in each block is for Z=Z0 and the first column in each line is for R=0.

Instead of the blocked input format, a pointwise data format can be used. This is introduced by a line containing the command "data", followed by the individual points of the map, one per line. for a grid field, each line contains values for X,Y,Z,Bx,By,Bz,Ex,Ey,Ez separated by either a comma and optional whitespece or by whitespace. for a cylinder field each line contains values for R,Z,Br,Bz,Er,Ez. The order of the points does not matter; omitted grid points will be 0.0, and for duplicates the last entry wins. If there is no E field, the Ex,Ey,Ez or Er,Ez entries should be omitted on every line. NOTE: every line's X,Y,Z or R,Z must be on a grid point as specified by the arguments to the grid or cylinder commands, to within the tolerance specified; if not, an error message is printed and the input line is ignored.

For time=dependent fields, the "time" command is used: time [period=12] period, if given, is in nanoseconds, and causes the interval [0,period) to be extended forever (before and after the values given). Because of the interpolation used, at least two points beyond the interval boundaries should be provided; there need not be a point at either boundary (but usually there are). Following the time command are lines containing 2 or 3 doubles: t B E where t is the time (nanoseconds), and B and E are factors for the fields. If E is omitted, the value for B is used. These values will be interpolated in time with a cubic spline that can handle either uniform or non-uniform spacing of points along t. The time command can come either before or after the cylinder or grid commands, but not within either of their sequences. Note a cubic spline is used to interpolate between points, and that can cause over/under-shoot near an abrupt change. Combined with period= this gives an excellent representation of sinewave/cosinewave.

Note that time dependence can currently only be specified via the time command in an input file (i.e. not programmable method exists).

Example block input file: * this is an example MyBLFieldMap input file, suitable for a solenoid # grid interval is 1 cm. # The region of validity is -390<=Z<=390 and 0<=R<=90 param normB=1.0 current=1.0 cylinder Z0=0.0 nR=10 nZ=40 dR=10.0 dZ=10.0 extendZ flip=Br Bz ... 40 lines of 10 values, Z=0 thru Z=390 Br ... 40 lines of 10 values, Z=0 thru Z=390 --EOF--

Example pointwise input file: * this is an example MyBLFieldMap input file, suitable for a solenoid # grid interval is 1 cm. # The region of validity is -390<=Z<=390 and 0<=R<=90 param normB=1.0 current=1.0 cylinder Z0=0.0 nR=10 nZ=40 dR=10.0 dZ=10.0 extendZ flip=Br data ... 400 lines of 4 values, giving R,Z,Br,Bz --EOF--

Definition at line 196 of file MyBLFieldMap.hh.

Constructor & Destructor Documentation

MyBLFieldMap::MyBLFieldMap ( )

default constructor.

Definition at line 267 of file MyBLFieldMap.cc.

References current, gradient, impl, maxline, normB, normE, and time.

00267                            {
00268     maxline = 1024;
00269     current = 1.0;
00270     gradient = 1.0;
00271     normB = 1.0;
00272     normE = 1.0;
00273     impl = 0;
00274     time = 0;
00275 }

MyBLFieldMap::~MyBLFieldMap ( ) [virtual]

destructor.

Definition at line 278 of file MyBLFieldMap.cc.

References impl, and time.

00278                             {
00279     if(impl) delete impl;
00280     impl = 0;
00281     if(time) delete time;
00282     time = 0;
00283 }

Member Function Documentation

bool MyBLFieldMap::createCylinderMap	(	G4double	Z0,
		G4double	dR,
		G4double	dZ,
		int	nR,
		int	nZ,
		class G4ElectroMagneticField *	field
	)

createCylinderMap() will create a 2-d map from a pre-defined cylindrically-symmetric field. NOTE: cylindrical symmetry is ASSUMED; the field in the Y=0,X>=0 half-plane is used (Bx->Br...). returns false if error.

Definition at line 481 of file MyBLFieldMap.cc.

References current, d2string(), gradient, i2string(), impl, maxline, normB, normE, and CylinderImpl::setField().

00483 00484 00485 00486     } 00487 00488 00489 00490 00491 00492 00493 00494 00495 00496 00497 00498 00499 00500 00501 00502 00503 00504 00505 00506 00507 00508 00509 00510 00511 00512 00513 00514 00515 00516 00517         } 00518     } 00519 00520 00521 }

class="fragment">00482 { if(impl) { delete impl; impl = 0; maxline = 128; current = 1.0; gradient = 1.0; normB = 1.0; normE = 1.0; MyBLArgumentVector argv; MyBLArgumentMap args; args["Z0"] = d2string(Z0); args["dR"] = d2string(dR); args["dZ"] = d2string(dZ); args["nR"] = i2string(nR); args["nZ"] = i2string(nZ); CylinderImpl *cyl = new CylinderImpl(argv,args); impl = cyl; // use the Y=0 plane for the R,Z plane bool retval = true; G4double pos[4], field[6]; pos[3] = 0.0; for(int i=0; i<nR; ++i) { pos[0] = 0.0 + i*dR; pos[1] = 0.0; for(int k=0; k<nZ; ++k) { pos[2] = Z0 + k*dZ; emField->GetFieldValue(pos,field); if(!cyl->setField(pos[0],pos[2],field[0],field[2], field[3],field[5],0)) retval = false; return false;

bool MyBLFieldMap::createGridMap	(	G4double	X0,
		G4double	Y0,
		G4double	Z0,
		G4double	dX,
		G4double	dY,
		G4double	dZ,
		int	nX,
		int	nY,
		int	nZ,
		class G4ElectroMagneticField *	field
	)

createGridMap() will create a 3-d map from a pre-defined field. returns false if error.

Definition at line 431 of file MyBLFieldMap.cc.

References current, d2string(), gradient, i2string(), impl, maxline, normB, normE, and GridImpl::setField().

00433                                  {
00434     if(impl) {
00435         delete impl;
00436         impl = 0;
00437     }
00438     maxline = 128;
00439     current = 1.0;
00440     gradient = 1.0;
00441     normB = 1.0;
00442     normE = 1.0;
00443 
00444     MyBLArgumentVector argv;
00445     MyBLArgumentMap args;
00446     args["X0"] = d2string(X0);
00447     args["Y0"] = d2string(Y0);
00448     args["Z0"] = d2string(Z0);
00449     args["dX"] = d2string(dX);
00450     args["dY"] = d2string(dY);
00451     args["dZ"] = d2string(dZ);
00452     args["nX"] = i2string(nX);
00453     args["nY"] = i2string(nY);
00454     args["nZ"] = i2string(nZ);
00455 
00456     GridImpl *grid = new GridImpl(argv,args);
00457     impl = grid;
00458 
00459     bool retval = true;
00460     G4double pos[4], field[6];
00461     pos[3] = 0.0;
00462     for(int i=0; i<nX; ++i) {
00463         pos[0] = X0 + i*dX;
00464         for(int j=0; j<nY; ++j) {
00465             pos[1] = Y0 + j*dY;
00466             for(int k=0; k<nZ; ++k) {
00467                 pos[2] = Z0 + k*dZ;
00468                 emField->GetFieldValue(pos,field);
00469                 if(!grid->setField(pos[0],pos[1],pos[2],
00470                     field[0],field[1],field[2],
00471                     field[3],field[4],field[5],0))
00472                     retval = false;
00473             }
00474         }
00475     }
00476 
00477     return retval;
00478 }

bool MyBLFieldMap::createTimeDependence	(	int	n,
		G4double	t[],
		G4double	b[],
		G4double	e[] = `0`,
		G4double	period = `-1.0`
	)

createTimeDependence() will apply the time dependence given. n is the # elements in the arrays; t[] is the time value for each point, b[] is the factor for B, and e[] is the factor for E. returns false if error.

Definition at line 552 of file MyBLFieldMap.cc.

References TimeImpl::setPeriod(), and time.

00553                                {
00554     if(time) return false;
00555     time = new TimeImpl(n,t,b,e);
00556     if(period > 0.0)
00557         time->setPeriod(period);
00558     return true;
00559 }

void MyBLFieldMap::getBoundingPoint	(	int	i,
		G4double	point[4]
	)

getBoundingPoint() returns the i-th bounding point of the map.

Definition at line 364 of file MyBLFieldMap.cc.

References FieldMapImpl::getBoundingPoint(), and impl.

00364                                                             {
00365     impl->getBoundingPoint(i,point);
00366 }

void MyBLFieldMap::getFieldValue	(	const G4double	local[4],
		G4double	field[6],
		G4double	current = `1.0`,
		G4double	gradient = `1.0`
	)

getFieldValue() gets the map's field to the field[] array. local[] must be LOCAL coordinates; local[3] (time) is ignored unless there is a time dependence specified. (elements using MyBLFieldMap perform the global->local conversion)

Definition at line 286 of file MyBLFieldMap.cc.

References current, TimeImpl::factorB(), TimeImpl::factorE(), FieldMapImpl::getFieldValue(), gradient, impl, normB, normE, and time.

Referenced by MyFieldMap::addFieldValue().

00287                                        {
00288     if(!impl)
00289         throw "MyBLFieldMap::getFieldValue called, no implementation";
00290 
00291     G4double thisField[6];
00292     impl->getFieldValue(local,thisField);
00293     G4double timeB=1.0, timeE=1.0;
00294     if(time) {
00295         timeB = time->factorB(local[3]);
00296         timeE = time->factorE(local[3]);
00297     }
00298     field[0] = thisField[0] * normB * timeB * _current/current;
00299     field[1] = thisField[1] * normB * timeB * _current/current;
00300     field[2] = thisField[2] * normB * timeB * _current/current;
00301     field[3] = thisField[3] * normE * timeE * _gradient/gradient;
00302     field[4] = thisField[4] * normE * timeE * _gradient/gradient;
00303     field[5] = thisField[5] * normE * timeE * _gradient/gradient;
00304 }

bool MyBLFieldMap::getTimeFactor	(	G4double	t,
		G4double *	b,
		G4double *	e
	)

getTimeFactor() returns the time factors for B and E at time t. returns false if error.

Definition at line 562 of file MyBLFieldMap.cc.

References TimeImpl::factorB(), TimeImpl::factorE(), and time.

00562                                                                      {
00563     if(b) *b = (time ? time->factorB(t) : 1.0);
00564     if(e) *e = (time ? time->factorE(t) : 1.0);
00565     return true;
00566 }

bool MyBLFieldMap::hasB ( )

hasB() returns true if this map has a nonzero B field.

Definition at line 369 of file MyBLFieldMap.cc.

References FieldMapImpl::hasB(), and impl.

Referenced by MyFieldMap::addFieldValue().

00369                         {
00370     return impl->hasB();
00371 }

bool MyBLFieldMap::hasE ( )

hasE() returns true if this map has a nonzero E field.

Definition at line 374 of file MyBLFieldMap.cc.

References FieldMapImpl::hasE(), and impl.

Referenced by MyFieldMap::addFieldValue().

00374                         {
00375     return impl->hasE();
00376 }

bool MyBLFieldMap::readFile ( G4String filename )

readFile() reads a file to initialize the map. Returns true if OK, false on error.

Definition at line 307 of file MyBLFieldMap.cc.

References argDouble(), argInt(), BLFuncs, InputFile::close(), current, InputFile::filename(), InputFile::getline(), InputFile::good(), gradient, FieldMapImpl::handleCommand(), impl, InputFile::linenumber(), maxline, normB, normE, MyBLFuncs::parseArgs(), TimeImpl::readTime(), InputFile::setMaxline(), and time.

Referenced by MyFieldMap::MyFieldMap().

00307                                              {
00308     InputFile in(filename,maxline);
00309     if(!in.good()) {
00310         throw (G4String("MyBLFieldMap::readFile: cannot open file ") + in.filename());
00311 
00312     }
00313     printf("MyBLFieldMap: reading file '%s'\n",in.filename());
00314 
00315     bool retval = true;
00316 
00317     char *line;
00318     while((line=in.getline()) != 0) {
00319         MyBLArgumentVector argv;
00320         MyBLArgumentMap namedArgs;
00321         if(BLFuncs.parseArgs(line,argv,namedArgs) < 0)
00322             goto invalid;
00323         if(argv[0] == "") continue;
00324         if(argv[0] == "param") {
00325             argInt(maxline,"maxline",namedArgs);
00326             argDouble(current,"current",namedArgs);
00327             argDouble(gradient,"gradient",namedArgs);
00328             argDouble(normB,"normB",namedArgs);
00329             argDouble(normE,"normE",namedArgs);
00330             in.setMaxline(maxline);
00331         }
00332         else if(argv[0] == "grid") {
00333             if(impl) goto invalid;
00334             impl = new GridImpl(argv,namedArgs);
00335         }
00336         else if(argv[0] == "cylinder") {
00337             if(impl) goto invalid;
00338             impl = new CylinderImpl(argv,namedArgs);
00339         }
00340         else if(argv[0] == "time") {
00341             if(time) goto invalid;
00342             time = TimeImpl::readTime(in,argv,namedArgs);
00343             if(!time) goto invalid;
00344         }
00345         else if(impl) {
00346             if(!impl->handleCommand(in,argv,namedArgs))
00347                 goto invalid;
00348         }
00349         else {
00350             invalid:        G4cerr << "MyBLFieldMap file "  << in.filename()
00351                 << " line " << in.linenumber()
00352                 << " invalid command " << argv[0].c_str()
00353                 << G4endl;
00354             retval = false;
00355             break;
00356         }
00357     }
00358     in.close();
00359 
00360     return retval;
00361 }

bool MyBLFieldMap::writeFile	(	G4String	filename,
		G4String	comment = `""`
	)

writeFile() writes the map to a file. Returns true if OK, false on error.

Definition at line 524 of file MyBLFieldMap.cc.

References current, gradient, impl, maxline, normB, normE, and FieldMapImpl::writeFile().

00524                                                                 {
00525     if(!impl) return false;
00526 
00527     FILE *f = fopen(filename.c_str(),"r");
00528     if(f) {
00529         fclose(f);
00530         G4Exception("MyBLFieldMap","Output File Exists",FatalException,
00531             filename.c_str());
00532     }
00533 
00534     f = fopen(filename,"w");
00535     if(!f) {
00536         fprintf(stderr,"MyBLFieldMap::writeFile CANNOT WRITE file '%s'\n",
00537             filename.c_str());
00538         return false;
00539     }
00540 
00541     fprintf(f,"# %s\n",comment.c_str());
00542     fprintf(f,"param maxline=%d current=%g gradient=%g normB=%g normE=%g\n",
00543         maxline,current,gradient,normB,normE);
00544 
00545     bool retval = impl->writeFile(f);
00546 
00547     fclose(f);
00548     return retval;
00549 }