 |
 |
| Commit in GeomConverter on MAIN | |||
| testResources/org/lcsim/geometry/field/Cartesian3DMagneticFieldMapTest.xml | +35 | added 1.1 | |
| test/org/lcsim/geometry/field/Cartesian3DMagneticFieldMapTest.java | +73 | added 1.1 | |
| src/org/lcsim/geometry/field/Cartesian3DMagneticFieldMap.java | +271 | added 1.1 | |
| +379 | |||
Implementation of a Cartesian 3D Magnetic Field Map with field values Bx, By, Bz defined on a regular grid of x,y,z using bilinear interpolation for intermediate points.
diff -N Cartesian3DMagneticFieldMapTest.xml --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ Cartesian3DMagneticFieldMapTest.xml 3 Sep 2013 22:04:01 -0000 1.1 @@ -0,0 +1,35 @@
+<lccdd xmlns:compact="http://www.lcsim.org/schemas/compact/1.0" + xmlns:xs="http://www.w3.org/2001/XMLSchema-instance" + xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd"> + + <info name="FieldMap3DTest"/> + + <define> + <constant name="cm" value="10"/> + + <constant name="world_side" value="30000" /> + <constant name="world_x" value="world_side" /> + <constant name="world_y" value="world_side" /> + <constant name="world_z" value="world_side" /> + + <constant name="tracking_region_radius" value="1100.0*cm"/> + <constant name="tracking_region_zmax" value="150.0*cm"/> + </define> + + <materials> + </materials> + + <detectors> + </detectors> + + <readouts> + </readouts> + + <fields> + <field type="Cartesian3DMagneticFieldMap" name="Cartesian3DMagneticFieldMapTest" + filename="./testResources/org/lcsim/geometry/field/HPS_b18d36_unfolded.dat" + offsetX="1.0" + offsetY="2.0" + offsetZ="3.0"/> + </fields> +</lccdd>
diff -N Cartesian3DMagneticFieldMapTest.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ Cartesian3DMagneticFieldMapTest.java 3 Sep 2013 22:04:01 -0000 1.1 @@ -0,0 +1,73 @@
+package org.lcsim.geometry.field;
+
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import java.io.InputStream;
+import java.util.Map;
+import static junit.framework.Assert.assertEquals;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.FieldMap;
+import org.lcsim.geometry.GeometryReader;
+import org.lcsim.geometry.compact.Field;
+
+/**
+ *
+ * @author Norman A Graf
+ *
+ * @version $Id:
+ */
+public class Cartesian3DMagneticFieldMapTest extends FieldTest
+{
+
+ public Cartesian3DMagneticFieldMapTest(String name)
+ {
+ super(name);
+ }
+
+ public void testRead() throws Exception
+ {
+ InputStream in = this.getClass().getResourceAsStream("/org/lcsim/geometry/field/Cartesian3DMagneticFieldMapTest.xml");
+ GeometryReader reader = new GeometryReader();
+ Detector det = reader.read(in);
+ FieldMap map = det.getFieldMap();
+
+ double[] off = {1., 2., 3.};
+ double[] fp = {-25.0, -8.9, -150.};
+ double[] B = new double[3];
+
+ // field at first map position should be zero since we offset the field
+ testFieldAt(map, fp[0], fp[1], fp[2], 0, 0, 0);
+ // field at first map position with offset included should equal first field values
+ testFieldAt(map, fp[0] + off[0], fp[1] + off[1], fp[2] + off[2], 0, -0.0019, 0);
+ // field at the origin
+ testFieldAt(map, 0.0 + off[0], -8.9 + off[1], 0.0 + off[2], 0., -0.5006, 0.);
+ // this field map is invariant in y, test this...
+ testFieldAt(map, 0.0 + off[0], 0. + off[1], 0.0 + off[2], 0., -0.5006, 0.);
+ testFieldAt(map, 0.0 + off[0], 8.9 + off[1], 0.0 + off[2], 0., -0.5006, 0.);
+
+ //TODO check interpolation more rigorously
+ //check all variations of accessor methods (why do we have SO many?!
+ double[] p = {fp[0] + off[0], fp[1] + off[1], fp[2] + off[2]};
+ Hep3Vector h3vP = new BasicHep3Vector(fp[0] + off[0], fp[1] + off[1], fp[2] + off[2]);
+
+ map.getField(p, B);
+ Hep3Vector h3vB = map.getField(h3vP);
+ assertEquals(B[0], h3vB.x());
+ assertEquals(B[1], h3vB.y());
+ assertEquals(B[2], h3vB.z());
+
+ double[] B2 = map.getField(p);
+ for (int i = 0; i < 3; ++i) {
+ assertEquals(B[i], B2[i]);
+ }
+
+ // test specifics of Cartesian3DMagneticFieldMap
+ Map<String, Field> fields = det.getFields();
+ assertEquals(fields.size(), 1);
+ Cartesian3DMagneticFieldMap cmap = (Cartesian3DMagneticFieldMap) fields.get("Cartesian3DMagneticFieldMapTest");
+ double[] offsets = cmap.globalOffset();
+ assertEquals(offsets[0], 1.);
+ assertEquals(offsets[1], 2.);
+ assertEquals(offsets[2], 3.);
+ }
+}
diff -N Cartesian3DMagneticFieldMap.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ Cartesian3DMagneticFieldMap.java 3 Sep 2013 22:04:01 -0000 1.1 @@ -0,0 +1,271 @@
+package org.lcsim.geometry.field;
+
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import java.io.BufferedReader;
+import java.io.FileInputStream;
+import java.io.IOException;
+import java.io.InputStream;
+import java.io.InputStreamReader;
+import static java.lang.Math.sqrt;
+import java.util.StringTokenizer;
+import org.jdom.Element;
+import org.jdom.JDOMException;
+
+/**
+ *
+ * @author Norman A Graf
+ *
+ * @version $Id:
+ */
+public class Cartesian3DMagneticFieldMap extends AbstractFieldMap
+{
+
+ private double[][][] _xField;
+ private double[][][] _yField;
+ private double[][][] _zField;
+ // The dimensions of the table
+ private int _nx, _ny, _nz;
+ // The physical limits of the defined region
+ private double _minx, _maxx, _miny, _maxy, _minz, _maxz;
+ // The physical extent of the defined region
+ private double _dx, _dy, _dz;
+ // Offsets if field map is not in global coordinates
+ private double _xOffset;
+ private double _yOffset;
+ private double _zOffset;
+ // maximum field strength
+ private double _bMax;
+ private double[] _Bfield = new double[3];
+ String _filename;
+
+ public Cartesian3DMagneticFieldMap(Element node) throws JDOMException
+ {
+ super(node);
+ _xOffset = node.getAttribute("offsetX").getDoubleValue();
+ _yOffset = node.getAttribute("offsetY").getDoubleValue();
+ _zOffset = node.getAttribute("offsetZ").getDoubleValue();
+
+ setup(node.getAttributeValue("filename"));
+ }
+
+ private void setup(String filename)
+ {
+ System.out.println("-----------------------------------------------------------");
+ System.out.println("Cartesian3DMagneticFieldMap ");
+ System.out.println("-----------------------------------------------------------");
+
+ System.out.println("Reading the field grid from " + filename + " ... ");
+
+ InputStream fis;
+ BufferedReader br;
+ String line;
+
+ try {
+ fis = new FileInputStream(filename);
+ br = new BufferedReader(new InputStreamReader(fis));
+ // ignore the first blank line
+ line = br.readLine();
+ // next line has table dimensions
+ line = br.readLine();
+ // read in the table dimensions of the file
+ StringTokenizer st = new StringTokenizer(line, " ");
+ _nx = Integer.parseInt(st.nextToken());
+ _ny = Integer.parseInt(st.nextToken());
+ _nz = Integer.parseInt(st.nextToken());
+
+
+ // Set up storage space for table
+ _xField = new double[_nx + 1][_ny + 1][_nz + 1];
+ _yField = new double[_nx + 1][_ny + 1][_nz + 1];
+ _zField = new double[_nx + 1][_ny + 1][_nz + 1];
+
+ // Ignore other header information
+ // The first line whose second character is '0' is considered to
+ // be the last line of the header.
+ do {
+ line = br.readLine();
+ st = new StringTokenizer(line, " ");
+ } while (!st.nextToken().trim().equals("0"));
+
+ // now ready to read in the values in the table
+ // format is:
+ // x y z Bx By Bz
+ //
+ int ix, iy, iz;
+ double xval = 0.;
+ double yval = 0.;
+ double zval = 0.;
+ double bx, by, bz;
+ for (ix = 0; ix < _nx; ix++) {
+ for (iy = 0; iy < _ny; iy++) {
+ for (iz = 0; iz < _nz; iz++) {
+ line = br.readLine();
+ st = new StringTokenizer(line, " ");
+ xval = Double.parseDouble(st.nextToken());
+ yval = Double.parseDouble(st.nextToken());
+ zval = Double.parseDouble(st.nextToken());
+ bx = Double.parseDouble(st.nextToken());
+ by = Double.parseDouble(st.nextToken());
+ bz = Double.parseDouble(st.nextToken());
+ if (ix == 0 && iy == 0 && iz == 0) {
+ _minx = xval;
+ _miny = yval;
+ _minz = zval;
+ }
+ _xField[ix][iy][iz] = bx;
+ _yField[ix][iy][iz] = by;
+ _zField[ix][iy][iz] = bz;
+ double b = bx * bx + by * by + bz * bz;
+ if (b > _bMax) {
+ _bMax = b;
+ }
+ }
+ }
+ }
+ _bMax = sqrt(_bMax);
+
+ _maxx = xval;
+ _maxy = yval;
+ _maxz = zval;
+
+ System.out.println("\n ---> ... done reading ");
+ System.out.println(" ---> assumed the order: x, y, z, Bx, By, Bz "
+ + "\n ---> Min values x,y,z: "
+ + _minx + " " + _miny + " " + _minz + " cm "
+ + "\n ---> Max values x,y,z: "
+ + _maxx + " " + _maxy + " " + _maxz + " cm "
+ + "\n Maximum Field strength: " + _bMax + " "
+ + "\n ---> The field will be offset by " + _xOffset + " " + _yOffset + " " + _zOffset + " cm ");
+
+ _dx = _maxx - _minx;
+ _dy = _maxy - _miny;
+ _dz = _maxz - _minz;
+ System.out.println("\n ---> Range of values x,y,z: "
+ + _dx + " " + _dy + " " + _dz + " cm in z "
+ + "\n-----------------------------------------------------------");
+
+ br.close();
+ } catch (IOException x) {
+ throw new RuntimeException(x);
+ }
+
+ }
+
+ @Override
+ public void getField(double[] position, double[] b)
+ {
+ getField(position[0], position[1], position[2]);
+ System.arraycopy(_Bfield, 0, b, 0, 3);
+ }
+
+ @Override
+ public Hep3Vector getField(Hep3Vector position)
+ {
+ getField(position.x(), position.y(), position.z());
+ return new BasicHep3Vector(_Bfield[0], _Bfield[1], _Bfield[2]);
+ }
+
+ @Override
+ public double[] getField(double[] position)
+ {
+ getField(position[0], position[1], position[2]);
+ double[] field = {_Bfield[0], _Bfield[1], _Bfield[2]};
+ return field;
+ }
+
+ @Override
+ void getField(double x, double y, double z, BasicHep3Vector field)
+ {
+ getField(x, y, z);
+ field.setV(_Bfield[0], _Bfield[1], _Bfield[2]);
+ }
+
+ public double[] globalOffset()
+ {
+ return new double[]{_xOffset, _yOffset, _zOffset};
+ }
+
+ private void getField(double x, double y, double z)
+ {
+ // allow for offsets
+ x -= _xOffset;
+ y -= _yOffset;
+ z -= _zOffset;
+ // Check that the point is within the defined region
+ if (x >= _minx && x <= _maxx
+ && y >= _miny && y <= _maxy
+ && z >= _minz && z <= _maxz) {
+
+ // Position of given point within region, normalized to the range
+ // [0,1]
+ double xfraction = (x - _minx) / _dx;
+ double yfraction = (y - _miny) / _dy;
+ double zfraction = (z - _minz) / _dz;
+
+ //double xdindex, ydindex, zdindex;
+
+ // Position of the point within the cuboid defined by the
+ // nearest surrounding tabulated points
+ double[] xmodf = modf(xfraction * (_nx - 1));
+ double[] ymodf = modf(yfraction * (_ny - 1));
+ double[] zmodf = modf(zfraction * (_nz - 1));
+
+ // The indices of the nearest tabulated point whose coordinates
+ // are all less than those of the given point
+
+ int xindex = (int) xmodf[0];
+ int yindex = (int) ymodf[0];
+ int zindex = (int) zmodf[0];
+ double xlocal = xmodf[1];
+ double ylocal = ymodf[1];
+ double zlocal = zmodf[1];
+ // bilinear interpolation
+ _Bfield[0] =
+ _xField[xindex][yindex][zindex] * (1 - xlocal) * (1 - ylocal) * (1 - zlocal)
+ + _xField[xindex][yindex][zindex + 1] * (1 - xlocal) * (1 - ylocal) * zlocal
+ + _xField[xindex][yindex + 1][zindex] * (1 - xlocal) * ylocal * (1 - zlocal)
+ + _xField[xindex][yindex + 1][zindex + 1] * (1 - xlocal) * ylocal * zlocal
+ + _xField[xindex + 1][yindex][zindex] * xlocal * (1 - ylocal) * (1 - zlocal)
+ + _xField[xindex + 1][yindex][zindex + 1] * xlocal * (1 - ylocal) * zlocal
+ + _xField[xindex + 1][yindex + 1][zindex] * xlocal * ylocal * (1 - zlocal)
+ + _xField[xindex + 1][yindex + 1][zindex + 1] * xlocal * ylocal * zlocal;
+ _Bfield[1] =
+ _yField[xindex][yindex][zindex] * (1 - xlocal) * (1 - ylocal) * (1 - zlocal)
+ + _yField[xindex][yindex][zindex + 1] * (1 - xlocal) * (1 - ylocal) * zlocal
+ + _yField[xindex][yindex + 1][zindex] * (1 - xlocal) * ylocal * (1 - zlocal)
+ + _yField[xindex][yindex + 1][zindex + 1] * (1 - xlocal) * ylocal * zlocal
+ + _yField[xindex + 1][yindex][zindex] * xlocal * (1 - ylocal) * (1 - zlocal)
+ + _yField[xindex + 1][yindex][zindex + 1] * xlocal * (1 - ylocal) * zlocal
+ + _yField[xindex + 1][yindex + 1][zindex] * xlocal * ylocal * (1 - zlocal)
+ + _yField[xindex + 1][yindex + 1][zindex + 1] * xlocal * ylocal * zlocal;
+ _Bfield[2] =
+ _zField[xindex][yindex][zindex] * (1 - xlocal) * (1 - ylocal) * (1 - zlocal)
+ + _zField[xindex][yindex][zindex + 1] * (1 - xlocal) * (1 - ylocal) * zlocal
+ + _zField[xindex][yindex + 1][zindex] * (1 - xlocal) * ylocal * (1 - zlocal)
+ + _zField[xindex][yindex + 1][zindex + 1] * (1 - xlocal) * ylocal * zlocal
+ + _zField[xindex + 1][yindex][zindex] * xlocal * (1 - ylocal) * (1 - zlocal)
+ + _zField[xindex + 1][yindex][zindex + 1] * xlocal * (1 - ylocal) * zlocal
+ + _zField[xindex + 1][yindex + 1][zindex] * xlocal * ylocal * (1 - zlocal)
+ + _zField[xindex + 1][yindex + 1][zindex + 1] * xlocal * ylocal * zlocal;
+
+ } else {
+ _Bfield[0] = 0.0;
+ _Bfield[1] = 0.0;
+ _Bfield[2] = 0.0;
+ }
+ }
+
+ //TODO pass double[] as argument to minimize internal array creation
+ private double[] modf(double fullDouble)
+ {
+ int intVal = (int) fullDouble;
+ double remainder = fullDouble - intVal;
+
+ double[] retVal = new double[2];
+ retVal[0] = intVal;
+ retVal[1] = remainder;
+
+ return retVal;
+ }
+}
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1