GeomConverter/src/org/lcsim/geometry/segmentation
diff -N EcalBarrelCartesianGridXY.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ EcalBarrelCartesianGridXY.java 18 Nov 2009 12:48:18 -0000 1.1
@@ -0,0 +1,425 @@
+package org.lcsim.geometry.segmentation;
+
+import java.util.ArrayList;
+import java.util.List;
+
+import org.jdom.DataConversionException;
+import org.jdom.Element;
+import org.lcsim.detector.IDetectorElement;
+import org.lcsim.detector.IDetectorElementContainer;
+import org.lcsim.detector.identifier.IExpandedIdentifier;
+import org.lcsim.detector.identifier.ExpandedIdentifier;
+import org.lcsim.detector.identifier.IIdentifier;
+import org.lcsim.detector.identifier.IIdentifierHelper;
+import org.lcsim.detector.identifier.Identifier;
+import org.lcsim.detector.solids.Box;
+import org.lcsim.geometry.util.IDDescriptor;
+
+/**
+ * XY Cartesian grid segmentation for overlapping staves.
+ *
+ * @author cassell
+ */
+public class EcalBarrelCartesianGridXY extends CartesianGridXY
+{
+ private double gridSizeX = 0;
+ private double gridSizeY = 0;
+
+ private int xIndex = -1;
+ private int yIndex = -1;
+ private int moduleIndex = -1;
+ private int sliceIndex = -1;
+ private int nmodules = 12;
+
+ private int[] validXplusG;
+ private int[] validXminusG;
+ private int[] validXplusP;
+ private int[] validXminusP;
+ private int validZplus;
+ private int validZminus;
+
+ private double xc0 = 0.;
+ private double[] yc;
+
+ private double sinth = 0.;
+ private double costh = 0.;
+ private double tanth = 0.;
+ private double cotth = 0.;
+ private double secth = 0.;
+ private double cscth = 0.;
+
+ private int nlayers = 0;
+ private int[] activeSlice;
+
+ public EcalBarrelCartesianGridXY(Element node) throws DataConversionException
+ {
+ super(node);
+
+ if (node.getAttribute("gridSizeX") != null)
+ {
+ gridSizeX = node.getAttribute("gridSizeX").getDoubleValue();
+ }
+ else
+ {
+ throw new RuntimeException("Missing gridSizeX parameter.");
+ }
+
+ if (node.getAttribute("gridSizeY") != null)
+ {
+ gridSizeY = node.getAttribute("gridSizeY").getDoubleValue();
+ }
+ else
+ {
+ throw new RuntimeException("Missing gridSizeY parameter.");
+ }
+ }
+
+ public long[] getNeighbourIDs(int layerRange, int xRange, int yRange)
+ {
+ if (this.getDecoder().getID() == 0)
+ throw new RuntimeException("No current ID is set.");
+
+ // Get the IdentifierHelper.
+ IIdentifierHelper helper = detector.getDetectorElement().getIdentifierHelper();
+
+ // Current packed id.
+ IIdentifier currId = new Identifier(this.getDecoder().getID());
+
+ // Create an ExpandedIdentifier for the current id.
+ IExpandedIdentifier thisId = helper.unpack(currId);
+// Check if we need initialization
+ if(validXplusP == null)
+ {
+ initializeMappings(this.getNumberOfLayers());
+ }
+ // Set values for current id.
+ int currLayer = helper.getValue(currId, layerIndex);
+ int currX = helper.getValue(currId, xIndex);
+ int currY = helper.getValue(currId, yIndex);
+ int currModule = helper.getValue(currId, moduleIndex);
+ // Create a new ExpandedIdenfier from the base Id.
+ ExpandedIdentifier neighborId = new ExpandedIdentifier(thisId);
+ // Find the neighbors within the current module
+ List<Long> nl = getNeighbourIDs(layerRange, xRange, yRange, helper, neighborId, currLayer, currX, currY);
+ // Check for border problems
+ boolean bp1 = false;
+ boolean bp2 = false;
+ // If the layer range includes layer 0 or negative, and the
+ // X bin exceeds the maximum X bin in a "pseudo" geometry
+ // then there is possible neighbors in module +1
+ if( ((currLayer - layerRange) < 1) && ((currX + xRange) > validXplusP[0]) )bp1 = true;
+ // If the X bin can be less than the minimum valid X bin for
+ // the minimum layer then there is possible neighbors in module
+ // -1
+ else if( (currX - xRange) <= validXminusG[Math.min(currLayer+layerRange,nlayers-1)])bp2 = true;
+ // otherwise no border problem and just return the neighbors in
+ // the current module.
+ else
+ {
+ long result[] = new long[nl.size()];
+ int i = 0;
+ for (Long id : nl)
+ {
+ result[i] = id;
+ i++;
+ }
+ return result;
+ }
+ if(bp1)
+ {
+ // Layer range includes l<=0, xbin large enough to overlap
+ // with module m+1. Neighbors in module m already added.
+ // Now find the neighbors in module m+1. The neighbor
+ // volume is a rectangle, so find cells is in module m+1
+ // cotained by the rectangle.
+ //
+ // Convert the corners of the rectangle layer, xbin to
+ // actual coordinates.
+ double lsep = yc[1] - yc[0];
+ // The maximum y coordinate is the layer 0 y coordinate +
+ // half the layer separation.
+ double ycmax = yc[0] + lsep/2.;
+ // The minimum y coordinate is the layer 0 y coordinate -
+ // half the layer separation - the layer separation*# negative
+ // layers.
+ double ycmin = yc[0] - lsep/2. +lsep*(currLayer - layerRange);
+ // The maximum/minimum X coordinate is the maximum/minimum
+ // Xbin+-.5 times the grid size, + x0
+ double xcmax = xc0 + gridSizeX*(currX + xRange + .5);
+ double xcmin = xc0 + gridSizeX*(currX - xRange - .5);
+ // rotate the coordinate system to module +1, and find the layer
+ // range in that module.
+ double minyp = ycmin*costh + xcmin*sinth;
+ double maxyp = ycmax*costh + xcmax*sinth;
+ int minl = 0;
+ int maxl = 0;
+ for(int il=0;il<nlayers;il++)
+ {
+ if(yc[il] < minyp)minl++;
+ if(yc[il] > maxyp)break;
+ maxl++;
+ }
+ // Set the module
+ int neighborModule = (currModule +1)%nmodules;
+ neighborId.setValue(moduleIndex, neighborModule);
+ // Loop over the layer range
+ for(int neighborLayer = minl;neighborLayer<maxl;neighborLayer++)
+ {
+ neighborId.setValue(layerIndex, neighborLayer);
+ // Ugly, need to set the slice bin
+ int iSlice = 2;
+ if(neighborLayer == 0)iSlice = 0;
+ neighborId.setValue(sliceIndex, iSlice);
+ // In this module's coordinate system, calculate the
+ // intersection of y=yc(layer) with the box to find the
+ // X coordinate range inside the box.
+ double xpmax = Math.min(xcmax*secth - yc[neighborLayer]*tanth,
+ yc[neighborLayer]*cotth - ycmin*cscth);
+ double xpmin = Math.max(xcmin*secth - yc[neighborLayer]*tanth,
+ yc[neighborLayer]*cotth - ycmax*cscth);
+ // Convert X coordinate to X bin.
+ int xbmin = (int) ((xpmin - xc0 + gridSizeX*.5)/gridSizeX);
+ int xbmax = (int) ((xpmax - xc0 - gridSizeX*.5)/gridSizeX);
+ // Loop over X bins
+ for(int neighborX=xbmin;neighborX <= xbmax;neighborX++)
+ {
+ neighborId.setValue(xIndex, neighborX);
+ // Loop over the Z range
+ for (int iy=-yRange; iy<=yRange; iy++)
+ {
+ // Compute y value.
+ int neighborY = currY + iy;
+ neighborId.setValue(yIndex, neighborY);
+ // Add the neighbor
+ nl.add(helper.pack(neighborId).getValue());
+ }
+ }
+ }
+ long result[] = new long[nl.size()];
+ int i = 0;
+ for (Long id : nl)
+ {
+ result[i] = id;
+ i++;
+ }
+ return result;
+ }
+ // To get here we have overlap with module m-1. Add the
+ // neighbors in module m-1.
+ int neighborModule = (nmodules+currModule -1)%nmodules;
+ neighborId.setValue(moduleIndex, neighborModule);
+ double ycmax = 0.;
+ double ycmin = 0.;
+ // The layer separation depends on where we are
+ if(currLayer - layerRange < 1)
+ {
+ double lsep = yc[1] - yc[0];
+ ycmin = yc[0] + (currLayer - layerRange)*lsep - lsep/2.;
+ }
+ else
+ {
+ double lsep = yc[currLayer - layerRange] - yc[currLayer - layerRange - 1];
+ ycmin = yc[currLayer - layerRange] - lsep/2.;
+ }
+ if(currLayer + layerRange > nlayers - 2)
+ {
+ double lsep = yc[nlayers - 1] - yc[nlayers - 2];
+ ycmax = yc[nlayers - 1] + lsep/2.;
+ }
+ else
+ {
+ double lsep = yc[currLayer + layerRange + 1] - yc[currLayer + layerRange];
+ ycmax = yc[currLayer + layerRange] + lsep/2.;
+ }
+ // Find the coordinates of the box in module m-1
+ double xcmax = xc0 + gridSizeX*(currX + xRange + .5);
+ double xcmin = xc0 + gridSizeX*(currX - xRange - .5);
+ double minyp = ycmin*costh - xcmax*sinth;
+ double maxyp = ycmax*costh - xcmin*sinth;
+ int minl = 0;
+ int maxl = 0;
+ for(int il=0;il<nlayers;il++)
+ {
+ if(yc[il] < minyp)minl++;
+ if(yc[il] > maxyp)break;
+ maxl++;
+ }
+ // Loop over layer range
+ for(int neighborLayer = minl;neighborLayer<maxl;neighborLayer++)
+ {
+ neighborId.setValue(layerIndex, neighborLayer);
+ int iSlice = 2;
+ if(neighborLayer == 0)iSlice = 0;
+ neighborId.setValue(sliceIndex, iSlice);
+ // Find intersection of y=yc[layer] with box to determine
+ // X range
+ double xpmax = Math.min(xcmax*secth + yc[neighborLayer]*tanth,
+ -yc[neighborLayer]*cotth + ycmax*cscth);
+ double xpmin = Math.max(xcmin*secth + yc[neighborLayer]*tanth,
+ -yc[neighborLayer]*cotth + ycmin*cscth);
+ int xbmin = (int) ((xpmin - xc0 + gridSizeX*.5)/gridSizeX);
+ int xbmax = (int) ((xpmax - xc0 - gridSizeX*.5)/gridSizeX);
+ xbmax = Math.min(xbmax,validXplusG[neighborLayer]);
+ // Loop and add neighbors
+ for(int neighborX=xbmin;neighborX <= xbmax;neighborX++)
+ {
+ neighborId.setValue(xIndex, neighborX);
+ for(int neighborY=Math.max(validZminus,currY-yRange);
+ neighborY<=Math.min(validZplus,currY+yRange);
+ neighborY++)
+ {
+ neighborId.setValue(yIndex, neighborY);
+ nl.add(helper.pack(neighborId).getValue());
+ }
+ }
+ }
+ long result[] = new long[nl.size()];
+ int i = 0;
+ for (Long id : nl)
+ {
+ result[i] = id;
+ i++;
+ }
+ return result;
+ }
+ protected List<Long> getNeighbourIDs(int layerRange, int xRange, int yRange, IIdentifierHelper helper, ExpandedIdentifier neighborId, int currLayer, int currX, int currY)
+ {
+ // Find neighbors within the currect module
+ List<Long> rl = new ArrayList<Long>();
+ for(int neighborLayer=Math.max(0,currLayer-layerRange);
+ neighborLayer<=Math.min(nlayers-1,currLayer+layerRange);
+ neighborLayer++)
+ {
+ neighborId.setValue(layerIndex, neighborLayer);
+ int iSlice = 2;
+ if(neighborLayer == 0)iSlice = 0;
+ neighborId.setValue(sliceIndex, iSlice);
+ for(int neighborX=Math.max(validXminusG[neighborLayer],currX-xRange);
+ neighborX<=Math.min(validXplusG[neighborLayer],currX+xRange);
+ neighborX++)
+ {
+ neighborId.setValue(xIndex, neighborX);
+ for(int neighborY=Math.max(validZminus,currY-yRange);
+ neighborY<=Math.min(validZplus,currY+yRange);
+ neighborY++)
+ {
+ neighborId.setValue(yIndex, neighborY);
+ rl.add(helper.pack(neighborId).getValue());
+ }
+ }
+ }
+ return rl;
+ }
+ public void setupGeomFields(IDDescriptor id)
+ {
+ super.setupGeomFields(id);
+ xIndex = id.indexOf("x");
+ yIndex = id.indexOf("y");
+ moduleIndex = id.indexOf("module");
+ layerIndex = id.indexOf("layer");
+ sliceIndex = id.indexOf("slice");
+ }
+
+ protected void initializeMappings(int nl)
+ {
+ // Initialize all the needed parameters for the neighbor
+ // finding across borders
+ // Get number of layers.
+ nlayers = nl;
+ // The valid X bins vary by layer
+ validXplusP = new int[this.getNumberOfLayers()];
+ validXminusP = new int[this.getNumberOfLayers()];
+ validXplusG = new int[this.getNumberOfLayers()];
+ validXminusG = new int[this.getNumberOfLayers()];
+ activeSlice = new int[this.getNumberOfLayers()];
+ // Calculate the trig functions needed for the rotations
+ sinth = Math.sin(2.*Math.PI/nmodules);
+ costh = Math.cos(2.*Math.PI/nmodules);
+ tanth = sinth/costh;
+ cotth = costh/sinth;
+ secth = 1./costh;
+ cscth = 1./sinth;
+// Get the helper and identifiers needed to manipulate the id
+ IIdentifierHelper helper = detector.getDetectorElement().getIdentifierHelper();
+ IIdentifier currId = new Identifier(this.getDecoder().getID());
+ IExpandedIdentifier thisId = helper.unpack(currId);
+ ExpandedIdentifier pseudoId = new ExpandedIdentifier(thisId);
+// Set the module and xbin to 0 to find pseudo-lengths and offsets
+ pseudoId.setValue(moduleIndex,0);
+ pseudoId.setValue(xIndex,0);
+ pseudoId.setValue(layerIndex,0);
+ pseudoId.setValue(sliceIndex,0);
+ // Save the current ID
+ long save = this.getDecoder().getID();
+ this.getDecoder().setID(helper.pack(pseudoId).getValue());
+ this.computeGlobalPosition();
+ // xc0 is the global x coordinate or the center of the module
+ xc0 = this.globalPosition[0] - gridSizeX/2.;
+ // yc is the global y coordinate of each layer
+ yc = new double[nlayers];
+ // xhlp is the half length of each layer
+ double[] xhlp = new double[nlayers];
+ for(int i=0;i<nlayers;i++)
+ {
+// Find the half lengths per layer and convert to x index.
+ pseudoId.setValue(layerIndex,i);
+ if(i == 1)pseudoId.setValue(sliceIndex,2);
+ this.getDecoder().setID(helper.pack(pseudoId).getValue());
+ this.computeGlobalPosition();
+ yc[i] = this.globalPosition[1];
+ // Find the box containing this cell
+ IIdentifier geomId = makeGeometryIdentifier(helper.pack(pseudoId).getValue());
+ IDetectorElementContainer deSrch = getSubdetector().getDetectorElement().findDetectorElement(geomId);
+ IDetectorElement de = deSrch.get(0);
+ if (de.getGeometry().getLogicalVolume().getSolid() instanceof Box)
+ {
+ Box sensorBox = (Box)de.getGeometry().getLogicalVolume().getSolid();
+ // Assume the Z extent is the same for all layers
+ if(i == 0)
+ {
+ double yhl = sensorBox.getYHalfLength();
+ // Convert half length to valid index
+ int nvalidy = (int) (yhl/gridSizeY) + 1;
+ validZplus = nvalidy - 1;
+ validZminus = -nvalidy;
+ }
+ // Convert X half length to valid X bins
+ double xhl = sensorBox.getXHalfLength();
+ int nvalidx = (int)(xhl/gridSizeX) + 1;
+ validXplusG[i] = nvalidx -1;
+ validXminusG[i] = -nvalidx;
+ // Compute valid X indices contained in a pseudo
+ // geometry of a regular N-gon
+ xhlp[i] = yc[i]*Math.tan(Math.PI/nmodules);
+ validXplusP[i] = (int)((xhlp[i] - xc0)/gridSizeX);
+ validXminusP[i] = -(int)((xhlp[i] + xc0)/gridSizeX) - 1;
+ }
+ else
+ {
+ throw new RuntimeException("Don't know how to bounds check solid " + de.getGeometry().getLogicalVolume().getSolid().getName() + ".");
+ }
+ }
+// Put back the cellID
+ this.getDecoder().setID(save);
+ this.computeGlobalPosition();
+ }
+ public int getPseudoLayer()
+ {
+ if(validXplusP == null)
+ {
+ initializeMappings(this.getNumberOfLayers());
+ }
+// Get the helper and identifiers needed to manipulate the id
+ IIdentifierHelper helper = detector.getDetectorElement().getIdentifierHelper();
+ IIdentifier currId = new Identifier(this.getDecoder().getID());
+ IExpandedIdentifier thisId = helper.unpack(currId);
+ int xbin = thisId.getValue(xIndex);
+ int layer = thisId.getValue(layerIndex);
+ if(xbin > validXplusP[layer])
+ {
+ // TODO: - calculate pseudoLayer
+ return layer;
+ }
+ return layer;
+ }
+}
GeomConverter/src/org/lcsim/geometry/segmentation
diff -N RegularNgonCartesianGridXY.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ RegularNgonCartesianGridXY.java 18 Nov 2009 12:48:18 -0000 1.1
@@ -0,0 +1,223 @@
+package org.lcsim.geometry.segmentation;
+
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+
+import java.util.ArrayList;
+import java.util.List;
+
+import org.jdom.DataConversionException;
+import org.jdom.Element;
+import org.lcsim.detector.IDetectorElement;
+import org.lcsim.detector.IDetectorElementContainer;
+import org.lcsim.detector.identifier.IExpandedIdentifier;
+import org.lcsim.detector.identifier.ExpandedIdentifier;
+import org.lcsim.detector.identifier.IIdentifier;
+import org.lcsim.detector.identifier.IIdentifierHelper;
+import org.lcsim.detector.identifier.Identifier;
+import org.lcsim.detector.solids.Box;
+import org.lcsim.detector.solids.Inside;
+import org.lcsim.geometry.util.IDDescriptor;
+
+/**
+ * Extend XY Cartesian grid segmentation for neighboring across
+ * module boundaries.
+ *
+ * @author cassell
+ */
+public class RegularNgonCartesianGridXY extends CartesianGridXY
+{
+ private double gridSizeX = 0;
+ private double gridSizeY = 0;
+ private int nmodules = 12;
+
+ private int xIndex = -1;
+ private int yIndex = -1;
+ private int moduleIndex = -1;
+
+ private int[] nvalidx;
+ private boolean[] borderCellIsDuplicate;
+
+ public RegularNgonCartesianGridXY(Element node) throws DataConversionException
+ {
+ super(node);
+
+ if (node.getAttribute("gridSizeX") != null)
+ {
+ gridSizeX = node.getAttribute("gridSizeX").getDoubleValue();
+ }
+ else
+ {
+ throw new RuntimeException("Missing gridSizeX parameter.");
+ }
+
+ if (node.getAttribute("gridSizeY") != null)
+ {
+ gridSizeY = node.getAttribute("gridSizeY").getDoubleValue();
+ }
+ else
+ {
+ throw new RuntimeException("Missing gridSizeY parameter.");
+ }
+ }
+
+ public long[] getNeighbourIDs(int layerRange, int xRange, int yRange)
+ {
+ if (this.getDecoder().getID() == 0)
+ throw new RuntimeException("No current ID is set.");
+
+ // Get the IdentifierHelper.
+ IIdentifierHelper helper = detector.getDetectorElement().getIdentifierHelper();
+
+ // Get number of layers.
+ int nlayers = this.getNumberOfLayers();
+
+ // Current packed id.
+ IIdentifier currId = new Identifier(this.getDecoder().getID());
+ // Create an ExpandedIdentifier for the current id.
+ IExpandedIdentifier thisId = helper.unpack(currId);
+ // Initialize the border arrays the first time in
+ if (nvalidx == null)
+ {
+ initializeBorders(thisId);
+ }
+ // Set values for current id.
+ int currLayer = helper.getValue(currId, layerIndex);
+ int currX = helper.getValue(currId, xIndex);
+ int currY = helper.getValue(currId, yIndex);
+ int currModule = helper.getValue(currId, moduleIndex);
+ // Create return array.
+ List<Long> neighbors = new ArrayList<Long>();
+ // Create a new ExpandedIdenfier from the base Id.
+ ExpandedIdentifier neighborId = new ExpandedIdentifier(thisId);
+ ExpandedIdentifier dupId = new ExpandedIdentifier(thisId);
+
+ // Loop over layer range.
+ for (int ilayer=-layerRange; ilayer<=layerRange; ilayer++)
+ {
+ // Compute layer value.
+ int neighborLayer = currLayer + ilayer;
+
+ if (neighborLayer >= 0 && neighborLayer < nlayers)
+ {
+ neighborId.setValue(layerIndex, neighborLayer);
+ // Loop over X range.
+ for (int ix=-xRange; ix<=xRange; ix++)
+ {
+ // Need to check border
+ int neighborX = currX + ix;
+ int neighborModule = currModule;
+ int dupX = neighborX;
+ int dupModule = neighborModule;
+ boolean dup = false;
+ if(neighborX < -nvalidx[neighborLayer]+1)
+ {
+ if(neighborX == -nvalidx[neighborLayer])
+ {
+ if(borderCellIsDuplicate[neighborLayer])
+ {
+ dup = true;
+ dupX = -neighborX -1;
+ dupModule = (nmodules + neighborModule -1)%nmodules;
+ }
+ }
+ else
+ {
+ neighborModule = (nmodules + neighborModule -1)%nmodules;
+ neighborX = 2*nvalidx[neighborLayer] + neighborX;
+ if(borderCellIsDuplicate[neighborLayer])neighborX--;
+ }
+ }
+ else if(neighborX > nvalidx[neighborLayer]-2)
+ {
+ if(neighborX == nvalidx[neighborLayer]-1)
+ {
+ if(borderCellIsDuplicate[neighborLayer])
+ {
+ dup = true;
+ dupX = -neighborX -1;
+ dupModule = (neighborModule +1)%nmodules;
+ }
+ }
+ else
+ {
+ neighborModule = (neighborModule +1)%nmodules;
+ neighborX = -2*nvalidx[neighborLayer] + neighborX;
+ if(borderCellIsDuplicate[neighborLayer])neighborX++;
+ }
+ }
+ neighborId.setValue(xIndex, neighborX);
+ // Loop over Y range.
+ for (int iy=-yRange; iy<=yRange; iy++)
+ {
+ // Compute y value.
+ int neighborY = currY + iy;
+ // Set the neighbor fields.
+ neighborId.setValue(yIndex, neighborY);
+ if(neighborModule != currModule)neighborId.setValue(moduleIndex, neighborModule);
+ // Add the neighbor id to the return array.
+ neighbors.add(helper.pack(neighborId).getValue());
+ // If we have a duplicate cell, add it
+ if(dup)
+ {
+ dupId.setValue(yIndex, neighborY);
+ dupId.setValue(xIndex, dupX);
+ dupId.setValue(layerIndex, neighborLayer);
+ dupId.setValue(moduleIndex, dupModule);
+ neighbors.add(helper.pack(dupId).getValue());
+ }
+ }
+ }
+ }
+ }
+
+ long result[] = new long[neighbors.size()];
+ int i = 0;
+ for (Long id : neighbors)
+ {
+ result[i] = id;
+ i++;
+ }
+ return result;
+ }
+ protected void initializeBorders(IExpandedIdentifier thisId)
+ {
+ // nvalidx is the number of valid xbins in the elements half
+ // length per layer.
+ nvalidx = new int[this.getNumberOfLayers()];
+ // borderCellIsDuplicate - true if we treat the partial cell\
+ // on each side of the boundary as a single cell
+ borderCellIsDuplicate = new boolean[this.getNumberOfLayers()];
+ IIdentifierHelper helper = detector.getDetectorElement().getIdentifierHelper();
+ ExpandedIdentifier testId = new ExpandedIdentifier(thisId);
+ for(int layer=0;layer<this.getNumberOfLayers();layer++)
+ {
+ testId.setValue(layerIndex, layer);
+ IIdentifier geomId = makeGeometryIdentifier(helper.pack(testId).getValue());
+ IDetectorElementContainer deSrch = getSubdetector().getDetectorElement().findDetectorElement(geomId);
+ IDetectorElement de = deSrch.get(0);
+ if (de.getGeometry().getLogicalVolume().getSolid() instanceof Box)
+ {
+ Box sensorBox = (Box)de.getGeometry().getLogicalVolume().getSolid();
+ double xhl = sensorBox.getXHalfLength();
+ nvalidx[layer] = (int)(xhl/gridSizeX) + 1;
+ double bcw = (nvalidx[layer]*gridSizeX - xhl)/gridSizeX;
+ borderCellIsDuplicate[layer] = false;
+ if(bcw < .5)borderCellIsDuplicate[layer] = true;
+ }
+ else
+ {
+ throw new RuntimeException("Don't know how to bounds check solid " + de.getGeometry().getLogicalVolume().getSolid().getName() + ".");
+ }
+ }
+ }
+ public void setupGeomFields(IDDescriptor id)
+ {
+ super.setupGeomFields(id);
+ xIndex = id.indexOf("x");
+ yIndex = id.indexOf("y");
+ layerIndex = id.indexOf("layer");
+ moduleIndex = id.indexOf("module");
+ }
+
+}
\ No newline at end of file
GeomConverter/src/org/lcsim/geometry/segmentation
diff -N RegularNgonCartesianGridXZ.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ RegularNgonCartesianGridXZ.java 18 Nov 2009 12:48:18 -0000 1.1
@@ -0,0 +1,243 @@
+package org.lcsim.geometry.segmentation;
+
+import java.util.ArrayList;
+import java.util.List;
+
+import org.jdom.DataConversionException;
+import org.jdom.Element;
+import org.lcsim.detector.IDetectorElement;
+import org.lcsim.detector.IDetectorElementContainer;
+import org.lcsim.detector.identifier.IExpandedIdentifier;
+import org.lcsim.detector.identifier.ExpandedIdentifier;
+import org.lcsim.detector.identifier.IIdentifier;
+import org.lcsim.detector.identifier.IIdentifierHelper;
+import org.lcsim.detector.identifier.Identifier;
+import org.lcsim.detector.solids.Trd;
+import org.lcsim.geometry.util.IDDescriptor;
+
+/**
+ * XZ Cartesian grid segmentation.
+ *
+ * @author cassell
+ */
+public class RegularNgonCartesianGridXZ extends CartesianGridXZ
+{
+ private double gridSizeX = 0;
+ private double gridSizeZ = 0;
+ private int nmodules = 12;
+
+ private int xIndex = -1;
+ private int zIndex = -1;
+ private int moduleIndex = -1;
+ private int barrelIndex = -1;
+
+ private int nzbins = 0;
+ private int zbinmax = 0;
+ private int zbinmin = 0;
+ private int[] nvalidx;
+ private boolean[] borderCellIsDuplicate;
+
+ private double xhl0 = 0.;
+ private double xslope = 0.;
+
+ private int ecsign = -1;
+ public RegularNgonCartesianGridXZ(Element node) throws DataConversionException
+ {
+ super(node);
+
+ if (node.getAttribute("gridSizeX") != null)
+ {
+ gridSizeX = node.getAttribute("gridSizeX").getDoubleValue();
+ }
+ else
+ {
+ throw new RuntimeException("Missing gridSizeX parameter.");
+ }
+
+ if (node.getAttribute("gridSizeZ") != null)
+ {
+ gridSizeZ = node.getAttribute("gridSizeZ").getDoubleValue();
+ }
+ else
+ {
+ throw new RuntimeException("Missing gridSizeZ parameter.");
+ }
+ }
+
+ public long[] getNeighbourIDs(int layerRange, int xRange, int zRange)
+ {
+ if (this.getDecoder().getID() == 0)
+ throw new RuntimeException("No current ID is set.");
+
+ // Get the IdentifierHelper.
+ IIdentifierHelper helper = detector.getDetectorElement().getIdentifierHelper();
+
+ // Get number of layers.
+ int nlayers = this.getNumberOfLayers();
+
+ // Current packed id.
+ IIdentifier currId = new Identifier(this.getDecoder().getID());
+ // Create an ExpandedIdentifier for the current id.
+ IExpandedIdentifier thisId = helper.unpack(currId);
+ // Initialize the border arrays the first time in
+ if (nvalidx == null)
+ {
+ initializeBorders(thisId);
+ }
+ // Set values for current id.
+ int currLayer = helper.getValue(currId, layerIndex);
+ int currX = helper.getValue(currId, xIndex);
+ int currZ = helper.getValue(currId, zIndex);
+ int currModule = helper.getValue(currId, moduleIndex);
+ // Create return array.
+ List<Long> neighbors = new ArrayList<Long>();
+ // Create a new ExpandedIdenfier from the base Id.
+ ExpandedIdentifier neighborId = new ExpandedIdentifier(thisId);
+ ExpandedIdentifier dupId = new ExpandedIdentifier(thisId);
+
+ // Loop over Z range.
+ for (int iz=-zRange; iz<=zRange; iz++)
+ {
+ // Compute y value.
+ int neighborZ = currZ + iz;
+ if(neighborZ >= zbinmin && neighborZ <= zbinmax)
+ {
+ neighborId.setValue(zIndex, neighborZ);
+ int Zindex = neighborZ + nzbins/2;
+ // Loop over X range.
+ for (int ix=-xRange; ix<=xRange; ix++)
+ {
+ // Need to check border
+ int neighborX = currX + ix;
+ int neighborModule = currModule;
+ int dupX = neighborX;
+ int dupModule = neighborModule;
+ boolean dup = false;
+ if(neighborX < -nvalidx[Zindex]+1)
+ {
+ if(neighborX == -nvalidx[Zindex])
+ {
+ if(borderCellIsDuplicate[Zindex])
+ {
+ dup = true;
+ dupX = -neighborX -1;
+ dupModule = (nmodules + neighborModule -ecsign)%nmodules;
+ }
+ }
+ else
+ {
+ neighborModule = (nmodules + neighborModule -ecsign)%nmodules;
+ neighborX = 2*nvalidx[Zindex] + neighborX;
+ if(borderCellIsDuplicate[Zindex])neighborX--;
+ }
+ }
+ else if(neighborX > nvalidx[Zindex]-2)
+ {
+ if(neighborX == nvalidx[Zindex]-1)
+ {
+ if(borderCellIsDuplicate[Zindex])
+ {
+ dup = true;
+ dupX = -neighborX -1;
+ dupModule = (nmodules + neighborModule +ecsign)%nmodules;
+ }
+ }
+ else
+ {
+ neighborModule = (nmodules + neighborModule +ecsign)%nmodules;
+ neighborX = -2*nvalidx[Zindex] + neighborX;
+ if(borderCellIsDuplicate[Zindex])neighborX++;
+ }
+ }
+ neighborId.setValue(xIndex, neighborX);
+ // Loop over layer range.
+ for (int ilayer=-layerRange; ilayer<=layerRange; ilayer++)
+ {
+ // Compute layer value.
+ int neighborLayer = currLayer + ilayer;
+
+ if (neighborLayer >= 0 && neighborLayer < nlayers)
+ {
+
+ // Set the neighbor fields.
+ neighborId.setValue(layerIndex, neighborLayer);
+ if(neighborModule != currModule)neighborId.setValue(moduleIndex, neighborModule);
+ // Add the neighbor id to the return array.
+ neighbors.add(helper.pack(neighborId).getValue());
+ // If we have a duplicate cell, add it
+ if(dup)
+ {
+ dupId.setValue(zIndex, neighborZ);
+ dupId.setValue(xIndex, dupX);
+ dupId.setValue(layerIndex, neighborLayer);
+ dupId.setValue(moduleIndex, dupModule);
+ neighbors.add(helper.pack(dupId).getValue());
+ }
+ }
+ }
+ }
+ }
+ }
+
+ long result[] = new long[neighbors.size()];
+ int i = 0;
+ for (Long id : neighbors)
+ {
+ result[i] = id;
+ i++;
+ }
+ return result;
+ }
+ protected void initializeBorders(IExpandedIdentifier thisId)
+ {
+ IIdentifierHelper helper = detector.getDetectorElement().getIdentifierHelper();
+ IIdentifier id = new Identifier(this.getDecoder().getID());
+ ExpandedIdentifier testId = new ExpandedIdentifier(thisId);
+ int barrelVal = helper.getValue(id, barrelIndex);
+ ecsign = 1;
+ if(barrelVal == 2)ecsign = -1;
+ testId.setValue(xIndex, 0);
+ IIdentifier geomId = makeGeometryIdentifier(helper.pack(testId).getValue());
+ IDetectorElementContainer deSrch = getSubdetector().getDetectorElement().findDetectorElement(geomId);
+ IDetectorElement de = deSrch.get(0);
+ if (de.getGeometry().getLogicalVolume().getSolid() instanceof Trd)
+ {
+ Trd sensorTrd = (Trd)de.getGeometry().getLogicalVolume().getSolid();
+ double zhl = sensorTrd.getZHalfLength();
+ double xhl1 = sensorTrd.getXHalfLength1();
+ double xhl2 = sensorTrd.getXHalfLength2();
+ xhl0 = (xhl1+xhl2)/2.;
+ xslope = (xhl2 - xhl0)*gridSizeZ/zhl;
+ nzbins = 2*((int)(zhl/gridSizeZ) + 1);
+ zbinmax = nzbins/2 -1;
+ zbinmin = -nzbins/2;
+ nvalidx = new int[nzbins];
+ borderCellIsDuplicate = new boolean[nzbins];
+ for(int zi=0;zi<nzbins;zi++)
+ {
+ int zbin = zi - nzbins/2;
+ double xhl = xhl0 + xslope*(zbin + .5);
+ nvalidx[zi] = (int)(xhl/gridSizeX) + 1;
+ double bcw = (nvalidx[zi]*gridSizeX - xhl)/gridSizeX;
+ borderCellIsDuplicate[zi] = false;
+ if(bcw < .5)borderCellIsDuplicate[zi] = true;
+ }
+ }
+ else
+ {
+ throw new RuntimeException("Don't know how to bounds check solid " + de.getGeometry().getLogicalVolume().getSolid().getName() + ".");
+ }
+ }
+
+ protected void setupGeomFields(IDDescriptor id)
+ {
+ super.setupGeomFields(id);
+ xIndex = id.indexOf("x");
+ zIndex = id.indexOf("z");
+ moduleIndex = id.indexOf("module");
+ layerIndex = id.indexOf("layer");
+ barrelIndex = id.indexOf("barrel");
+
+ }
+
+}
\ No newline at end of file
