GeomConverter/src/org/lcsim/geometry/segmentation

diff -u -r1.6 -r1.7
--- NonprojectiveCylinder.java	28 Jun 2005 23:15:05 -0000	1.6
+++ NonprojectiveCylinder.java	5 Jul 2005 06:48:57 -0000	1.7
@@ -19,6 +19,8 @@

 import org.lcsim.geometry.util.IDDescriptor;
 import org.lcsim.geometry.util.IDEncoder;
 import org.lcsim.geometry.subdetector.CylindricalCalorimeter;

+import org.lcsim.geometry.layer.Layer;
+import org.lcsim.geometry.layer.LayerStack;

 
 /**
  * @author jeremym

@@ -74,14 +76,6 @@

 	return gridSizeZ;
     }
 

-    public double getZMin() {
-	return ((CylindricalCalorimeter)detector).getZMin();
-    }
-
-    public double getZMax() {
-	return ((CylindricalCalorimeter)detector).getZMax();
-    }
-

     public double getPhi()
     {
         return (((double)getValue(phiIndex)) + 0.5) * computeDeltaPhiForLayer();

@@ -152,44 +146,200 @@

     public long[] getNeighbourIDs(int layerRange, int zRange, int phiRange)
     {
       encoder.setValues(values);

+//       int klay = decoder.getValue("layer");
+//       int kz   = decoder.getValue("z");
+//       int kphi = decoder.getValue("phi");
+//       System.out.println("NeighborID: ref="+klay+" "+kz+" "+kphi+" (hex "
+// 			 +Long.toHexString(encoder.getID())+")");

 
       int nMax = (2*layerRange + 1)*(2*zRange + 1)*(2*phiRange + 1) - 1;
       long[] result = new long[nMax];
       

+      // theta and phi are used to find central neighbors in other layers
+      double theta = getTheta();
+      double phi = getPhi();
+

       int zBins = (int)Math.floor( (getZMax()-getZMin()) / getGridSizeZ() );

-      System.out.println("In z: zmin="+getZMin()+", zmax="+getZMax()
-			 +", gridZ= "+getGridSizeZ()+", zBins="+zBins);

 
       int size = 0;

-      for (int i=-layerRange;i<=layerRange;i++) {
-	int l = values[layerIndex] + i;

+      for (int i=-layerRange; i<=layerRange; ++i) {
+	int ilay = values[layerIndex] + i;
+
+	if (ilay<0 || ilay>=detector.getLayerCount()) continue;
+	encoder.setValue(layerIndex,ilay);
+
+	double dphi = this.computeDeltaPhiForLayer(ilay);
+	int phiBins =  (int)Math.floor(2*Math.PI / dphi);

-	if (l<0 || l>=detector.getLayerCount()) continue;
-	encoder.setValue(layerIndex,l);

+	double cylR = getRadiusSensitiveMid( ilay );
+	double x = cylR * Math.cos(phi);
+	double y = cylR * Math.sin(phi);
+	double z = cylR / Math.tan(theta);
+// 	System.out.println("layer "+ilay+", cylR="+cylR
+// 			   +", theta="+theta+", phi="+phi
+// 			   +", rvec=( "+x+", "+y+", "+z+")");
+
+	long id = this.findCellContainingXYZ(x,y,z);
+	if(id==0) continue;
+
+	// save indices in a new array, as values[] keeps the original ref
+	decoder.setID(id);
+	int[] tmpval = new int[ values.length ];
+        decoder.getValues(tmpval);

-	double dphi = this.computeDeltaPhiForLayer(i);
-	int phiBins =  (int)Math.floor(2*Math.PI / dphi) ;
-	for (int j=-zRange;j<=zRange;j++) {
-          int t = values[zIndex] + j;

+	for (int j=-zRange; j<=zRange; ++j) {
+          int iz = tmpval[zIndex] + j;

-	  if (t<0 || t>=zBins) continue;
-	  encoder.setValue(zIndex,t);

+	  if (iz<-zBins/2 || iz>=zBins/2) continue;
+	  encoder.setValue(zIndex,iz);

-	  for (int k=-phiRange;k<=phiRange;k++) {

+	  for (int k=-phiRange; k<=phiRange; ++k) {
+	    // skip reference cell

 	    if (i==0 && j==0 && k==0) continue;
 

-	    int p = values[phiIndex] + k;
-	    if (p<0 || p>=phiBins) continue;

+	    int iphi = tmpval[phiIndex] + k;
+	    if(iphi<0) iphi += phiBins;
+	    if(iphi>=phiBins) iphi -= phiBins;
+	    if (iphi<0 || iphi>=phiBins) continue; 

-	    result[size++] = encoder.setValue(phiIndex,p);

+// 	    System.out.println("Adding neighbor: "+ilay+" "+iz+" "+iphi);
+	    result[size++] = encoder.setValue(phiIndex,iphi);

 	  }
 	}
       }
       if(size < result.length) {

+	int klay = decoder.getValue("layer");
+	int kz   = decoder.getValue("z");
+	int kphi = decoder.getValue("phi");
+// 	System.out.println("Neighbors="+size+"/"+result.length
+// 			   +", ref="+klay+" "+kz+" "+kphi);
+

 	long[] temp = new long[size];
 	System.arraycopy(result,0,temp,0,size);
 	result = temp;
       }
       return result;
     }

+
+    // Not for public use, but this is needed to calculate positions
+    // and number of cells, etc.
+    private double getZMin() {
+	return ((CylindricalCalorimeter)detector).getZMin();
+    }
+
+    private double getZMax() {
+	return ((CylindricalCalorimeter)detector).getZMax();
+    }
+
+    private double getRMin() {
+	return ((CylindricalCalorimeter)detector).getInnerRadius();
+    }
+
+    private double getRMax() {
+	return ((CylindricalCalorimeter)detector).getOuterRadius();
+    }
+
+    /**
+     * Return the cell which contains a given point (x,y,z), or null.
+     * Side effect: encoder will be reset to the ID returned by this method.
+     *
+     * @param x,y,z cartesian coordinates of the point
+     * @return ID of cell containing the point (maybe either in
+     * absorber or live material)
+     */
+    public long findCellContainingXYZ(double x, double y, double z) {
+	// validate point
+	if(z<getZMin()) {
+	    System.out.println("find: fail z<zmin: z="+z);
+	    return 0;
+	}
+	if(z>getZMax()) {
+	    System.out.println("find: fail z>zmax: z="+z);
+	    return 0;
+	}
+	double rho = Math.sqrt(x*x+y*y);
+	if(rho<getRMin()) {
+	    System.out.println("find: fail rho<Rmin: rho="+rho);
+	    return 0;
+	}
+	if(rho>getRMax()) {
+	    System.out.println("find: fail rho>Rmax: rho="+rho);
+	    return 0;
+	}
+
+	// ok, point is valid, so a valid ID should be returned
+	int ilay = getLayerBin(rho);
+	int iz = getZBin(z);
+
+	double phi = Math.atan2( y, x );
+	if(phi<0) phi += 2*Math.PI;
+	int iphi = getPhiBin(ilay,phi);
+
+	encoder.setValue("layer",ilay);
+	encoder.setValue(zIndex,iz);
+	encoder.setValue(phiIndex,iphi);
+
+	return encoder.getID();
+    }
+
+    public int getLayerBin(double r) {
+	// In order to be general, we should not assume that all
+	// layers have the same thickness.  Therefore, one has to
+	// guess the starting layer (based on average thickness), and
+	// then navigate trhough layers until one finds the right one
+	double depth = r - getRMin();
+	double mean_t = (getRMax()-getRMin()) / detector.getLayerCount();
+
+	int ilay = (int)Math.floor( depth / mean_t );
+	LayerStack stack = detector.getLayers();
+	Layer layer = stack.getLayer(ilay);
+	double depHi = stack.getThicknessToLayerBack(ilay);
+	double depLo = depHi - layer.getTotalThickness();
+	for(;;) {
+	  if( depth>depLo && depth<=depHi ) return ilay;
+	  if( depth<=depLo ) {
+	      --ilay;
+	      depHi = depLo;
+	      layer = stack.getLayer(ilay);
+	      depLo -= layer.getTotalThickness();
+	  }
+	  if( depth>depHi ) {
+	      ++ilay;
+	      depLo = depHi;
+	      layer = stack.getLayer(ilay);
+	      depHi += layer.getTotalThickness();
+	  }
+	}
+    }
+
+    public int getPhiBin(int ilay, double phi) {
+      // phic is phi at center of cells with iphi=0
+      double deltaPhi = this.computeDeltaPhiForLayer(ilay);
+      double phic = deltaPhi / 2;
+      int iphi = (int)Math.floor((phi-phic)/deltaPhi + 0.5);
+      return iphi;
+    }
+
+    public int getZBin(double z) {
+      // zc is z at center of cells with iz=0
+      int numz = (int)Math.floor((getZMax()-getZMin()) / gridSizeZ);
+      double zc = 0;  // for odd numz
+      if( numz%2 == 0 ) zc = gridSizeZ / 2;
+      int iz = (int)Math.floor((z-zc)/gridSizeZ + 0.5);
+      return iz;
+    }
+
+    /**
+     * Returns cylindrical radius to center of sensitive slice of any layer
+     * @param layer layer index
+     */
+    private double getRadiusSensitiveMid(int ilay) {
+	LayerStack stack = detector.getLayers();
+	Layer layer = stack.getLayer(ilay);
+
+	double preLayers = 0;
+	if(ilay>0) preLayers = stack.getSectionThickness(0,ilay-1);
+
+	return this.getRMin() + preLayers + layer.getThicknessToSensitiveMid();
+    }

 }