LISTSERV mailing list manager LISTSERV 16.5

Help for HPS-SVN Archives


HPS-SVN Archives

HPS-SVN Archives


HPS-SVN@LISTSERV.SLAC.STANFORD.EDU


View:

Message:

[

First

|

Previous

|

Next

|

Last

]

By Topic:

[

First

|

Previous

|

Next

|

Last

]

By Author:

[

First

|

Previous

|

Next

|

Last

]

Font:

Proportional Font

LISTSERV Archives

LISTSERV Archives

HPS-SVN Home

HPS-SVN Home

HPS-SVN  January 2016

HPS-SVN January 2016

Subject:

r4075 - in /java/trunk/users/src/main/java/org/hps/users/spaul: ./ feecc/ feecc/analysis/

From:

[log in to unmask]

Reply-To:

Notification of commits to the hps svn repository <[log in to unmask]>

Date:

Sun, 3 Jan 2016 18:18:11 -0000

Content-Type:

text/plain

Parts/Attachments:

Parts/Attachments

text/plain (1971 lines)

Author: [log in to unmask]
Date: Sun Jan  3 10:18:03 2016
New Revision: 4075

Log: (empty)

Added:
    java/trunk/users/src/main/java/org/hps/users/spaul/feecc/FEESpectrumGenerator.java
    java/trunk/users/src/main/java/org/hps/users/spaul/feecc/FormFactor.java
    java/trunk/users/src/main/java/org/hps/users/spaul/feecc/MottIntegral.java
    java/trunk/users/src/main/java/org/hps/users/spaul/feecc/MultipleScattering.java
    java/trunk/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinningXY.java
Removed:
    java/trunk/users/src/main/java/org/hps/users/spaul/feecc/analysis/
Modified:
    java/trunk/users/src/main/java/org/hps/users/spaul/StyleUtil.java
    java/trunk/users/src/main/java/org/hps/users/spaul/SumEverything.java
    java/trunk/users/src/main/java/org/hps/users/spaul/feecc/BinGenerator.java
    java/trunk/users/src/main/java/org/hps/users/spaul/feecc/CustomBinning.java
    java/trunk/users/src/main/java/org/hps/users/spaul/feecc/DisplayHistograms.java
    java/trunk/users/src/main/java/org/hps/users/spaul/feecc/EcalUtil.java
    java/trunk/users/src/main/java/org/hps/users/spaul/feecc/MakeHistograms.java
    java/trunk/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinning.java

Modified: java/trunk/users/src/main/java/org/hps/users/spaul/StyleUtil.java
 =============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/StyleUtil.java	(original)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/StyleUtil.java	Sun Jan  3 10:18:03 2016
@@ -1,10 +1,9 @@
 package org.hps.users.spaul;
 
 import java.util.Arrays;
+import java.util.Random;
 
-import hep.aida.IAnalysisFactory;
-import hep.aida.IFunction;
-import hep.aida.IPlotterRegion;
+import hep.aida.*;
 
 public class StyleUtil {
 	
@@ -27,18 +26,20 @@
 		r.style().yAxisStyle().setLabel(ly);
 		r.style().yAxisStyle().labelStyle().setFontSize(16);
 		r.style().yAxisStyle().tickLabelStyle().setFontSize(14);
-		r.style().statisticsBoxStyle().setParameter("backgroundColor", "white");
-		r.style().statisticsBoxStyle().boxStyle().backgroundStyle().setColor("White");
-		r.style().statisticsBoxStyle().boxStyle().backgroundStyle().setParameter("color", "white");
-		r.style().statisticsBoxStyle().boxStyle().backgroundStyle().setOpacity(100);
-		System.out.println(Arrays.toString(r.style().dataStyle().availableParameters()));
+		//	r.style().statisticsBoxStyle().set;
+		//System.out.println(Arrays.toString());
+		r.style().legendBoxStyle().textStyle().setFontSize(16);
+		
+		//r.style().dataStyle().showInLegendBox(false);
+		r.style().legendBoxStyle().boxStyle().foregroundStyle().setOpacity(1.0);
 		r.style().dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
 		r.style().dataStyle().fillStyle().setParameter("showZeroHeightBins", "false");
+		System.out.println(Arrays.toString(r.style().statisticsBoxStyle().boxStyle().backgroundStyle().availableParameters()));
 		//r.style().dataStyle().setParameter("showDataInStatisticsBox", "false");
 		r.style().setParameter("hist2DStyle", "colorMap");
 		//r.style().dataBoxStyle()
 	}
-	static void addHorizontalLine(IPlotterRegion r, double y, String name){
+	public static void addHorizontalLine(IPlotterRegion r, double y, String name){
 		IAnalysisFactory af = IAnalysisFactory.create();
 		IFunction f = af.createFunctionFactory(af.createTreeFactory().create()).createFunctionByName(name, "p0");
 		f.setParameter("p0", y);
@@ -53,4 +54,54 @@
 		f.setParameter("p2", p2);
 		region.plot(f);
 	}
+	public static void noFillHistogramBars(IPlotterRegion region) {
+		region.style().dataStyle().setParameter("fillHistogramBars", "false");
+		region.style().dataStyle().fillStyle().setVisible(false);
+		
+
+		region.style().dataStyle().lineStyle().setParameter("colorRotateMethod", "regionOverlayIndex");
+		region.style().dataStyle().lineStyle().setParameter("colorRotate", "black, red, green, blue");
+		region.style().dataStyle().lineStyle().setParameter("thickness", "3");
+		region.style().dataStyle().errorBarStyle().setVisible(false);
+		System.out.println(Arrays.toString(region.style().dataStyle().lineStyle().availableParameterOptions("colorRotateMethod")));
+	}
+	
+	public static void main(String arg[]){
+		IAnalysisFactory af = IAnalysisFactory.create();
+		IHistogramFactory hf = af.createHistogramFactory(af.createTreeFactory().create());
+		
+		IPlotter p = af.createPlotterFactory().create();
+		p.createRegions(1, 2);
+		IHistogram1D h1 = hf.createHistogram1D("blah", 100, -5, 5);
+		IHistogram1D h2 = hf.createHistogram1D("bleh", 100, -5, 5);
+		Random random = new Random();
+		for(int i = 0; i< 100000; i++){
+			h1.fill(random.nextGaussian());
+			h2.fill(random.nextGaussian()*2);
+		}
+		hideLegendAndStats(p.region(1));
+		noFillHistogramBars(p.region(0));
+		stylize(p.region(0), "title", "x axis label", "y axis label");
+		stylize(p.region(1), "stuff", "x axis label", "y axis label");
+		p.region(0).plot(h1);
+		p.region(0).plot(h2);
+		
+
+		IHistogram2D h3 = hf.createHistogram2D("blah", 100, -5, 5, 100, -5,5);
+		
+		for(int i = 0; i< 100000; i++){
+			h3.fill(random.nextGaussian(), random.nextGaussian());
+		}
+		
+		p.region(1).plot(h3);
+		
+		
+		
+		p.show();
+		
+	}
+	public static void hideLegendAndStats(IPlotterRegion r){
+		r.style().statisticsBoxStyle().setVisible(false);
+		r.style().legendBoxStyle().setVisible(false);
+	}
 }

Modified: java/trunk/users/src/main/java/org/hps/users/spaul/SumEverything.java
 =============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/SumEverything.java	(original)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/SumEverything.java	Sun Jan  3 10:18:03 2016
@@ -2,6 +2,7 @@
 
 import java.io.File;
 import java.io.IOException;
+import java.util.ArrayList;
 
 import hep.aida.IAnalysisFactory;
 import hep.aida.IHistogram1D;
@@ -23,20 +24,44 @@
 			oneArg(arg[0]);
 		}
 	}
-	static void oneArg(String indir) throws IllegalArgumentException, IOException{
+	static void oneArg(final String indir) throws IllegalArgumentException, IOException{
 		File outdir = new File(indir + "/sums");
 		outdir.mkdir();
-		for(File subdirf : new File(indir).listFiles()){
-			String subdir = subdirf.getAbsolutePath();
-			if(subdir.matches(".*/sums/?"))
-				continue;
-			String split[] = subdir.split("/");
-			String outfile = indir + "/sums/" + split[split.length-1] + ".aida";
-			twoArg(subdir, outfile);
+		ArrayList<Thread> threads = new ArrayList<Thread>();
+		for(final File subdirf : new File(indir).listFiles()){
+			Thread t = new Thread(){
+				public void run(){
+
+					String subdir = subdirf.getAbsolutePath();
+					if(subdir.matches(".*/sums/?"))
+						return;
+					String split[] = subdir.split("/");
+					String outfile = indir + "/sums/" + split[split.length-1] + ".aida";
+					try {
+						twoArg(subdir, outfile);
+					} catch (IllegalArgumentException e) {
+						// TODO Auto-generated catch block
+						e.printStackTrace();
+					} catch (IOException e) {
+						// TODO Auto-generated catch block
+						e.printStackTrace();
+					}
+				}
+			};
+			threads.add(t);
+			t.start();
+		}
+		for(Thread t : threads){
+			try {
+				t.join();
+			} catch (InterruptedException e) {
+				// TODO Auto-generated catch block
+				e.printStackTrace();
+			}
 		}
 		new File(indir + "/sums/total.aida").delete();
 		twoArg(outdir.getAbsolutePath(), indir + "/sums/total.aida");
-		
+
 	}
 
 	static void twoArg(String indir, String out) throws IllegalArgumentException, IOException{
@@ -45,45 +70,56 @@
 		new File(out).delete();
 		ITree tree0 = tf.create(out, "xml", false, true);
 		IHistogramFactory hf = af.createHistogramFactory(tree0);
-		
+
 
 		int j = 0;
 		for(File s : new File(indir).listFiles()){
-			ITree tree = af.createTreeFactory().create(s.getAbsolutePath(),"xml");
+			if(!s.getAbsolutePath().endsWith("aida"))
+				continue;
+			try{
+				ITree tree = af.createTreeFactory().create(s.getAbsolutePath(),"xml");
 
-			
-			if(j == 0){
-				String [] names = tree.listObjectNames();
-				tree0.mount("/tmp", tree, "/");
-				for(String name : names){
-					Object o = tree.find(name);
-					if(o instanceof IHistogram1D || o instanceof IHistogram2D)
-						tree0.cp("/tmp/" + name, name);
+
+				if(j == 0){
+					String [] names = tree.listObjectNames();
+					tree0.mount("/tmp", tree, "/");
+					for(String name : names){
+						Object o = tree.find(name);
+						if(o instanceof IHistogram1D || o instanceof IHistogram2D)
+							tree0.cp("/tmp/" + name, name);
+					}
+					tree0.unmount("/tmp");
+					tree.close();
+
 				}
-				tree0.unmount("/tmp");
+				else{
+					//tree.
+					String [] names = tree.listObjectNames();
+					tree0.mount("/tmp", tree, "/");
+					for(String name : names){
+						Object o = tree.find(name);
+						if(o instanceof IHistogram1D)
+							((IHistogram1D)tree0.find(name)).add((IHistogram1D)o);
+						if(o instanceof IHistogram2D)
+							((IHistogram2D)tree0.find(name)).add((IHistogram2D)o);
+					}
+					tree0.unmount("/tmp");
+					tree.close();
+				}
+
 				tree.close();
-				
+				j++;
+				System.out.println(j + " files have been read");
+
+
+			} catch(IllegalArgumentException e){
+				//print the filename
+				System.out.println("Exception happened at file " + s.getAbsolutePath());
+
+				e.printStackTrace();
 			}
-			else{
-				//tree.
-				String [] names = tree.listObjectNames();
-				tree0.mount("/tmp", tree, "/");
-				for(String name : names){
-					Object o = tree.find(name);
-					if(o instanceof IHistogram1D)
-						((IHistogram1D)tree0.find(name)).add((IHistogram1D)o);
-					if(o instanceof IHistogram2D)
-						((IHistogram2D)tree0.find(name)).add((IHistogram2D)o);
-				}
-				tree0.unmount("/tmp");
-				tree.close();
-			}
-
-			tree.close();
-			j++;
-			System.out.println(j + " files have been read");
+			tree0.commit();
 		}
-		tree0.commit();
 	}
 }
 

Modified: java/trunk/users/src/main/java/org/hps/users/spaul/feecc/BinGenerator.java
 =============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/feecc/BinGenerator.java	(original)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/feecc/BinGenerator.java	Sun Jan  3 10:18:03 2016
@@ -6,7 +6,7 @@
 
 public class BinGenerator {
 	public static void main(String arg[]) throws FileNotFoundException{
-		int nBins = 20;
+		int nBins = 32;
 
 		PrintStream pw = new PrintStream("generatedbins.txt");
 		pw.println(nBins);
@@ -33,17 +33,21 @@
 			bins[i] = thetaMin+dTheta*i;
 		}
 		return bins; */
-		double thetaMax = .145;
-		double thetaMin = .035;
+		double thetaMax = .200;
+		double thetaMin = .040;
 
 		double[] bins = new double[nBins +1];
-		double xMin = 1/(thetaMax*thetaMax);
+		for(int i = 0; i<nBins+1; i++){
+			bins[i] = thetaMin+i*(thetaMax-thetaMin)/nBins;
+		}
+		return bins;
+		/*double xMin = 1/(thetaMax*thetaMax);
 		double xMax = 1/(thetaMin*thetaMin);
 		for(int i = 0; i< nBins+1; i++){
 			double x = xMax - i*(xMax-xMin)/nBins;
 			bins[i] = Math.pow(x, -.5);
 		}
-		return bins;
+		return bins;*/
 	}
 	private static double[] getPhiBounds(double thetaMin, double thetaMax){
 		double phiBins[] = new double[6];
@@ -52,8 +56,14 @@
 
 		boolean prevInRange = false;
 		for(double phi = 0; phi< 3.14; phi+= dphi){
-			boolean inRange = EcalUtil.fid_ECal_spherical(thetaMin, phi) && EcalUtil.fid_ECal_spherical(thetaMax, phi)
-								&& EcalUtil.fid_ECal_spherical(thetaMin, -phi) && EcalUtil.fid_ECal_spherical(thetaMax, -phi);
+			
+			// make the angular cuts on the tracks such that the particles that go into that cut 
+			// are expected to be within 4 mm (~= 2 times the angular resolution of 1.5 mrad) of 
+			// the ecal cuts.  
+			double d = 4;
+			
+			boolean inRange = EcalUtil.fid_ECal_spherical_more_strict(thetaMin, phi, d) && EcalUtil.fid_ECal_spherical_more_strict(thetaMax, phi, d)
+								&& EcalUtil.fid_ECal_spherical_more_strict(thetaMin, -phi, d) && EcalUtil.fid_ECal_spherical_more_strict(thetaMax, -phi, d);
 			if(inRange && !prevInRange)
 				phiBins[edgeNumber++] = phi;
 			if(prevInRange && !inRange)

Modified: java/trunk/users/src/main/java/org/hps/users/spaul/feecc/CustomBinning.java
 =============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/feecc/CustomBinning.java	(original)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/feecc/CustomBinning.java	Sun Jan  3 10:18:03 2016
@@ -76,21 +76,5 @@
 		}
 		return tot;
 	}
-	/**
-	 * @param bin
-	 * @param a = 2E/M
-	 * @return the integral of 1/sin^4(th/2)*cos^2(th/2)/(1+a*sin^2(th/2)) times dPhi,
-	 * which appears in the integral of mott scattering.
-	 */
-	public double mottIntegralFactor(double a, int bin){
-		double dPhi = 0;
-		for(int i = 0; i< phiMax[bin].length; i++)
-			dPhi += 2*(phiMax[bin][0] - phiMin[bin][0]); //factor of 2 from top and bottom
-		
-		double csc2 = Math.pow(Math.sin(thetaMax[bin]/2), -2);
-		double Imax = (-csc2+(1+a)*Math.log(a+2*csc2));
-		       csc2 = Math.pow(Math.sin(thetaMin[bin]/2), -2);
-		double Imin = (-csc2+(1+a)*Math.log(a+2*csc2));
-		return 2*dPhi*(Imax-Imin);
-	}
+
 }

Modified: java/trunk/users/src/main/java/org/hps/users/spaul/feecc/DisplayHistograms.java
 =============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/feecc/DisplayHistograms.java	(original)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/feecc/DisplayHistograms.java	Sun Jan  3 10:18:03 2016
@@ -4,6 +4,7 @@
 
 public class DisplayHistograms {
 	public static void main(String arg[]) throws IllegalArgumentException, IOException{
+		//System.out.println("dognabo");
 		MakeHistograms.main(new String[]{arg[0]});
 	}
 }

Modified: java/trunk/users/src/main/java/org/hps/users/spaul/feecc/EcalUtil.java
 =============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/feecc/EcalUtil.java	(original)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/feecc/EcalUtil.java	Sun Jan  3 10:18:03 2016
@@ -12,6 +12,45 @@
 import org.lcsim.event.Cluster;
 
 public class EcalUtil {
+	
+	static int[] getCrystalIndex(Cluster c){
+		if(map == null){
+			try{
+				readMap();
+			}catch(Exception e){
+				e.printStackTrace();
+			}
+		}
+		c.getPosition();
+		int bestix = 0, bestiy = 0;
+		double bestdist = 100000;
+
+		double cx = c.getPosition()[0];
+		double cy = c.getPosition()[1];
+		for(int ix = -23; ix<= 23; ix++){
+			if(!map.containsKey(200+ix))
+				continue;
+			double x = map.get(200+ix)[0];
+			if(Math.abs(x-cx)<bestdist){
+				bestdist = Math.abs(x-cx);
+				bestix = ix;
+			}
+		}
+		
+		bestdist = 100000;
+		for(int iy = -5; iy<= 5; iy++){
+			if(!map.containsKey(100*iy + bestix))
+				continue;
+			double y = map.get(100*iy + bestix)[1];
+			if(Math.abs(y-cy)<bestdist){
+				bestdist = Math.abs(y-cy);
+				bestiy = iy;
+			}
+		}
+		return new int[]{bestix, bestiy}; 
+	}
+	
+	
 	/**
 	 * Holly's algorithm.  
 	 * @param x x position of cluster
@@ -83,6 +122,40 @@
 
 		return in_fid;
 	}
+	/**
+	 * 
+	 * @param x
+	 * @param y
+	 * @param d the additional distance from the edge of the ecal in addition
+	 *       to what is required by fid_Cal(double, double)
+	 * @return
+	 */
+	public static boolean fid_ecal_more_strict(double x, double y, double d){
+		y = Math.abs(y);
+
+		boolean in_fid = false;
+		 
+		double x_edge_low = -262.74 + d;
+		double x_edge_high = 347.7 - d;
+		double y_edge_low = 33.54 + d;
+		double y_edge_high = 75.18 - d;
+
+		double x_gap_low = -106.66 - d;
+		double x_gap_high = 42.17 + d;
+		double y_gap_high = 47.18 + d;
+
+		y = Math.abs(y);
+
+		if( x > x_edge_low && x < x_edge_high && y > y_edge_low && y < y_edge_high )
+		{
+			if( !(x > x_gap_low && x < x_gap_high && y > y_edge_low && y < y_gap_high) )
+			{
+				in_fid = true;
+			}
+		}
+
+		return in_fid;
+	}
 	static double[] toSphericalFromBeam(double pxpz, double pypz){
 		double x = pxpz, y = pypz, z = 1;
 		double beamTilt = .03057;
@@ -96,10 +169,10 @@
 		
 		return new double[]{theta, phi};
 	}
-	static Map<Integer, double[]> map = new HashMap();
+	static Map<Integer, double[]> map ;
 	static void readMap() throws FileNotFoundException{
-		Scanner s = new Scanner(new File("ecal_positions.txt"));
-
+		Scanner s = new Scanner(new File(System.getenv("HOME") + "/ecal_positions.txt"));
+		map = new HashMap();
 		
 		while(s.hasNext()){
 			int ix =s.nextInt();
@@ -154,6 +227,12 @@
 		double y = xy[1];
 		return fid_ECal(x, y);
 	}
+	public static boolean fid_ECal_spherical_more_strict(double theta, double phi, double d){
+		double[] xy = getXY(theta, phi);
+		double x = xy[0];
+		double y = xy[1];
+		return fid_ecal_more_strict(x, y, d);
+	}
 	public static void main(String arg[]){
 		double x = 0, y = 0;
 		double sp[] = getThetaPhiSpherical(x,y); 

Added: java/trunk/users/src/main/java/org/hps/users/spaul/feecc/FEESpectrumGenerator.java
 =============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/feecc/FEESpectrumGenerator.java	(added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/feecc/FEESpectrumGenerator.java	Sun Jan  3 10:18:03 2016
@@ -0,0 +1,294 @@
+package org.hps.users.spaul.feecc;
+
+import hep.aida.*;
+import hep.io.stdhep.StdhepEvent;
+import hep.io.stdhep.StdhepRecord;
+import hep.io.stdhep.StdhepWriter;
+
+import java.io.File;
+import java.io.IOException;
+import java.util.Random;
+import static java.lang.Math.*;
+
+public class FEESpectrumGenerator {
+	/*************
+	 * CONSTANTS *
+	 *************/
+	static double alpha = 1/137., hbarc = .197, q = 1.6e-19, Na= 6.022e23;
+	static double me = .000511;
+	static double amu = 931.4941;
+
+	private static double sqrt2 = Math.sqrt(2);
+
+	/***************
+	 * OTHER STUFF *
+	 ***************/
+
+	static int counter = 0;
+
+	static Random random = new Random();
+	static double thetaMax = .1;
+	static double thetaMin = .03;
+	/**
+	 * 
+	 * @return a random angle theta, with no form factor.  
+	 */
+	public static double getRandomThetaBare(){
+		double r = random.nextDouble();
+		double dTheta = (thetaMax-thetaMin)/4;
+		double theta = (thetaMax+thetaMin)/2;
+		int nIter = 11;
+		for(int i = 0; i< nIter; i++){
+			double x = MottIntegral.integral(a, theta);
+			double f = (x-xmin)/(xmax-xmin);
+			if(f<r){
+				theta += dTheta;
+			} else if(f> r){
+				theta -= dTheta;
+			}
+			dTheta /= 2;
+		}
+		return theta;
+
+
+	}
+
+	static double combine(double theta, double rms){
+		return theta + sqrt2 *rms*random.nextGaussian();
+	}
+
+
+	/**
+	 * 
+	 * @param a recoil factor = 2*E/M
+	 * @return a random angle theta, with no form factor.  
+	 */
+	public static double getRandomThetaWithFormFactor(FormFactor f){
+		while(true){
+			double theta = getRandomThetaBare();
+			double trigstuff = pow(sin(theta/2), 2);
+			double Q2 = 4*E*E*trigstuff/(1+a*trigstuff);
+			counter ++;
+			double f2 = f.getFormFactorSquared(Q2);
+			if(f2 > random.nextDouble())
+				return theta;
+		}
+	}
+
+
+
+		static double b = -1;
+	
+
+
+
+	static boolean radCorrections = false;
+	static int Nevents;
+	static double M, E, Q, aDensity;
+	static int Z;
+	static double a, xmax, xmin;
+	static boolean display;
+
+
+	private static CustomBinning cb;
+	public static void main(String arg[]) throws IOException{
+		String filename = arg[0];
+		M = Double.parseDouble(arg[1]);
+		E = Double.parseDouble(arg[2]);
+		thetaMin = Double.parseDouble(arg[3]);
+		thetaMax = Double.parseDouble(arg[4]);
+		Z = Integer.parseInt(arg[5]);
+
+		String outfile = null;
+		aDensity = Double.parseDouble(arg[6]);
+		Q = Double.parseDouble(arg[7]);
+
+		radCorrections = Boolean.parseBoolean(arg[8]);
+
+		if(arg.length>9){
+			display = true;
+			cb = new CustomBinning(new File(arg[9]));
+			outfile = arg[10];
+		}
+
+
+		a = 2*E/(M*amu);
+		xmin = MottIntegral.integral(a, thetaMin);
+		xmax = MottIntegral.integral(a, thetaMax);
+
+		double scale = Math.pow((Z*alpha*hbarc)/(2*E), 2)/100; //this should be in barns
+		double sigma = MottIntegral.mottIntegralWithFormFactor(a, scale, thetaMin, thetaMax, FormFactor.get(Z), 1000, E);
+		
+		double luminosity = Q*1e-9/q*aDensity/(M/Na)*1e-24;
+		Nevents = (int)(sigma*luminosity);
+
+
+		System.out.println(Nevents);
+
+		String title = "fee";
+		String comment = "fee";
+		StdhepWriter writer = new StdhepWriter(filename, title, comment, Nevents);
+		IAnalysisFactory af = null;
+		ITree tree = null;
+
+		if(display){
+			af = IAnalysisFactory.create();
+			tree = af.createTreeFactory().create(outfile,"xml",false,true);
+			IHistogramFactory hf = af.createHistogramFactory(tree);
+			
+			double thetaBins[] = new double[cb.nTheta+1];
+			for(int i = 0; i<cb.nTheta; i++){
+				thetaBins[i] = cb.thetaMin[i];
+			}
+
+			thetaBins[thetaBins.length-1] = cb.thetaMax[cb.nTheta-1];
+			
+			thetaHist = hf.createHistogram1D("theta", "theta", thetaBins);
+			
+			EHist = hf.createHistogram1D("energy", 50, 0, 1.3);
+			IPlotter p = af.createPlotterFactory().create();
+			p.createRegions(2,1);
+			p.region(0).plot(thetaHist);
+			p.region(1).plot(EHist);
+			p.show();
+		}
+
+		for(int i = 1; i<= Nevents; i++){
+			int nevhep = i; 
+			int nhep = 1;
+			int[] isthep = {1}; 
+			int[] idhep = {11};
+			int[] jmohep = {0,0}; 
+			int[] jdahep = {0,0};
+			double[] phep = randomMomentum();
+			double[] vhep = {0,0,0,0};
+			
+			
+			if(phep[3] > .5*E){
+				StdhepEvent event = new StdhepEvent(nevhep, nhep, isthep, idhep, jmohep, jdahep, phep, vhep);
+				if(display)
+					EHist.fill(phep[3]*(1+random.nextGaussian()*.045));
+				writer.writeRecord(event);
+			}
+			else{ // if the event has enough energy loss,
+				  // write it out as a blank event.  This way,
+				  // slic doesnt have to deal with low energy electrons
+				  // 
+				StdhepEvent event = new StdhepEvent(i, 0, 
+						new int[0],
+						new int[0],
+						new int[0],
+						new int[0],
+						new double[0],
+						new double[0]);
+				writer.writeRecord(event);
+			}
+			if(i%1000 == 0)
+				System.out.println("event " + i);
+		}
+		writer.close();
+		if(tree != null)
+			tree.commit();
+		
+		/*
+		String title = arg[0+10*set];
+		String inputFile = arg[1+10*set];
+		CustomBinning cb = new CustomBinning(new File(arg[2+10*set]));
+
+
+		double Q = Double.parseDouble(arg[3+10*set]); //in nC
+		double E = Double.parseDouble(arg[4+10*set]); //in GeV
+		double aDensity = Double.parseDouble(arg[5+10*set]); //in g/cm^2
+		double M = Double.parseDouble(arg[6+10*set]); // in amu
+		double Z = Double.parseDouble(arg[7+10*set]); 
+		double r = Double.parseDouble(arg[8+10*set]);
+		double prescale = Math.pow(2,Double.parseDouble(arg[9+10*set]));
+		 */
+		if(display){
+			ExtractFormFactors.main(new String[]{
+					"generated",
+					outfile,
+					arg[9],
+					Double.toString(Q),
+					Double.toString(E),
+					Double.toString(aDensity),
+					Double.toString(M),
+					Double.toString(Z),
+					"3", //not used
+					"0" //no prescale
+
+			});
+
+			displayResiduals();
+
+		}
+	}
+
+	private static void displayResiduals() {
+		IAnalysisFactory af = IAnalysisFactory.create();
+		IDataPointSetFactory dpsf = af.createDataPointSetFactory(af.createTreeFactory().create());
+		IDataPointSet residuals = dpsf.create("residuals", 2);
+		for(int i = 0; i<ExtractFormFactors.dpsMeasured[0].size(); i++){
+			IDataPoint gen = ExtractFormFactors.dpsMeasured[0].point(i);
+			IDataPoint the = ExtractFormFactors.dpsTheory.point(i);
+
+			IDataPoint res = residuals.addPoint();
+			res.coordinate(0).setValue(gen.coordinate(0).value());
+			res.coordinate(1).setValue(gen.coordinate(1).value()-the.coordinate(1).value());
+			double err = Math.hypot(gen.coordinate(1).errorMinus(), the.coordinate(1).errorMinus());
+			res.coordinate(1).setErrorMinus(err);
+			res.coordinate(1).setErrorPlus(err);
+		}
+		IPlotter p = af.createPlotterFactory().create();
+		p.region(0).plot(residuals);
+		p.show();
+
+	}
+
+
+
+	static IHistogram1D thetaHist = null, EHist;
+
+
+
+	private static double[] randomMomentum() {
+		double theta = getRandomThetaWithFormFactor(FormFactor.get(Z));
+		double phi = random.nextDouble()*2*PI-PI;
+		if(thetaHist != null)
+			if(cb.inRange(theta, phi))
+				thetaHist.fill(theta);
+		double trigstuff = pow(sin(theta/2.), 2);
+		double Ep = E/(1+a*trigstuff);
+		double m = .000511;
+
+
+
+		if(radCorrections){
+
+			// L.W. Mo, Y.S. Tsai, Radiative Corrections to Elastic and inelastic ep and $\mu$p scattering
+			// Stanford Linear Accelerator Center, Stanford University, Stanford, CA
+			// Reviews of Modern Physics  Volume 41 Number 1, January 1969
+
+			//otherwise apply radiative corrections
+			double Q2 = 4*E*Ep*trigstuff;
+
+			//effective number of pre and post radiation lengths caused by radiative effects
+			double t = 3/4.*(alpha/Math.PI)*(Math.log(Q2/(me*me)-1));
+
+			// rms theta for 
+			double theta_rms = MultipleScattering.getThetaRMS(t, E);
+			theta = combine(theta, theta_rms);
+			
+			if(b == -1){
+				b = MultipleScattering.b(Z);
+			}
+			double bt = b*t;
+			Ep *= MultipleScattering.getRandomEnergyFraction(bt)*MultipleScattering.getRandomEnergyFraction(bt);
+			
+		}
+		if(Double.isNaN(Ep))
+			return null;
+		double p = sqrt(Ep*Ep-m*m);
+		return new double[]{p*sin(theta)*cos(phi), p*sin(theta)*sin(phi), p*cos(theta), Ep, m};
+	}
+}

Added: java/trunk/users/src/main/java/org/hps/users/spaul/feecc/FormFactor.java
 =============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/feecc/FormFactor.java	(added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/feecc/FormFactor.java	Sun Jan  3 10:18:03 2016
@@ -0,0 +1,43 @@
+package org.hps.users.spaul.feecc;
+
+public abstract class FormFactor {
+
+	abstract double getFormFactor(double q2);
+
+	double getFormFactorSquared(double q2){
+		return Math.pow(getFormFactor(q2), 2);
+	}
+	static double hbarc = .197;
+	static FormFactor carbon = new FormFactor(){
+
+		@Override
+		public double getFormFactor(double Q2) {
+			double Z= 6;
+			double rp = 0.8786;
+			double a = 1.64;
+			double b = Math.sqrt(a*a*(1-1/12.)+rp*rp);
+			return (1-(Z-2)/(6*Z)*a*a*Q2/(hbarc*hbarc))*Math.exp(-1/4.*b*b*Q2/(hbarc*hbarc));
+		}
+
+	};
+	static FormFactor tungsten = new FormFactor(){
+
+		@Override
+		public double getFormFactor(double Q2) {
+			double Z= 74;
+			double r = 6.87;
+			double x = Math.sqrt(Q2*r*r/(hbarc*hbarc));
+			double F_calc = 3/(x*x*x)*(Math.sin(x)-x*Math.cos(x));
+			return F_calc;
+		}
+		
+	};
+	static FormFactor get(int Z){
+		if(Z == 6)
+			return carbon;
+		if(Z == 74)
+			return tungsten;
+		System.err.println("Form Factor for " + Z + " not implemented");
+		return null;
+	}
+}

Modified: java/trunk/users/src/main/java/org/hps/users/spaul/feecc/MakeHistograms.java
 =============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/feecc/MakeHistograms.java	(original)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/feecc/MakeHistograms.java	Sun Jan  3 10:18:03 2016
@@ -2,12 +2,15 @@
 
 import java.io.File;
 import java.io.IOException;
+import java.util.ArrayList;
+import java.util.Arrays;
 import java.util.Collections;
 import java.util.Comparator;
 import java.util.Date;
 import java.util.List;
 
 import org.hps.conditions.ConditionsDriver;
+import org.hps.recon.tracking.TrackType;
 import org.hps.record.triggerbank.AbstractIntData;
 import org.hps.record.triggerbank.TIData;
 import org.lcsim.event.CalorimeterHit;
@@ -15,9 +18,12 @@
 import org.lcsim.event.EventHeader;
 import org.lcsim.event.GenericObject;
 import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.Track;
+import org.lcsim.event.TrackState;
 import org.lcsim.lcio.LCIOReader;
 
 import hep.aida.IAnalysisFactory;
+import hep.aida.IDataPointSetFactory;
 import hep.aida.IHistogram1D;
 import hep.aida.IHistogram2D;
 import hep.aida.IHistogramFactory;
@@ -29,6 +35,9 @@
 import org.hps.users.spaul.StyleUtil;
 
 public class MakeHistograms {
+	
+	
+	
 	static boolean display = false;
 	static CustomBinning cb;
 	public static void main(String arg[]) throws IllegalArgumentException, IOException{
@@ -40,11 +49,11 @@
 				path = "temp.aida";
 			}
 			IAnalysisFactory af = IAnalysisFactory.create();
-			 ITreeFactory tf = af.createTreeFactory();
+			ITreeFactory tf = af.createTreeFactory();
 			ITree tree0 = tf.create(path, "xml");
 			extractHistograms(tree0);
 			setupPlotter(af);
-			
+
 		} else{
 
 			String input = arg[0];
@@ -86,104 +95,206 @@
 		}
 
 	}
-	
+
 	static IHistogram2D h1, h2, h2a, h2b, h2c;
 	static IHistogram2D h4,h4a;
-	static IHistogram1D h3, h3a;
-	static IHistogram1D h5, h6;
-	
+	static IHistogram1D h3, h3a, h3_t, h3_b;
+	static IHistogram1D h5, h5a; 
+	static IHistogram2D h6, h6a;
+	static IHistogram1D h7, h7a;
+	static IHistogram1D h8;
+	static IHistogram1D h9_t, h9_b;
+	static IHistogram1D h10_t, h10_b;
+
 	private static void extractHistograms(ITree tree0) {
+		h1 = (IHistogram2D) tree0.find("theta vs energy");
+
 		h2 = (IHistogram2D) tree0.find("theta vs phi");
 		h2a = (IHistogram2D) tree0.find("theta vs phi cut");
 		h2b = (IHistogram2D) tree0.find("theta vs phi cut alt");
 		h2c = (IHistogram2D) tree0.find("theta vs phi alt");
-		
+
 		h3 = (IHistogram1D) tree0.find("theta");
+		h3a = (IHistogram1D) tree0.find("theta isolated ");
+		h3_t = (IHistogram1D) tree0.find("theta top");
+		h3_b = (IHistogram1D) tree0.find("theta bottom");
+		
 		h4 = (IHistogram2D) tree0.find("px\\/pz vs py\\/pz");
 		h4a = (IHistogram2D) tree0.find("px\\/pz vs py\\/pz cut");
-		h5 = (IHistogram1D) tree0.find("energy");
-		
+		System.out.println(h4a.xAxis().bins());
+		h5 = (IHistogram1D) tree0.find("energy top");
+		h5 = (IHistogram1D) tree0.find("energy bottom");
+
+		h6 = (IHistogram2D) tree0.find("cluster");
+		h6a = (IHistogram2D) tree0.find("cluster matched");
+		h7 = (IHistogram1D) tree0.find("y top");
+		h7a = (IHistogram1D) tree0.find("y bottom");
+		h8 = (IHistogram1D) tree0.find("seed energy");
+		
+
+		h9_t = (IHistogram1D) tree0.find("pz top");
+		h9_b = (IHistogram1D) tree0.find("pz bottom");
+		
+
+		h10_t = (IHistogram1D) tree0.find("clustsize top");
+		h10_b = (IHistogram1D) tree0.find("clustsize bottom");
+
 	}
 	static void setupHistograms(IHistogramFactory hf){
 		//h1 = hf.createHistogram2D("px\\/pz vs py\\/pz", 160, -.16, .24, 160, -.2, .2);
-		
-		
-		
-		
+
+
+
+
 		h2 = hf.createHistogram2D("theta vs phi", 300, 0, .3, 314, -3.14, 3.14);
-		
+
 		h2a = hf.createHistogram2D("theta vs phi cut", 300, 0, .3, 314, -3.14, 3.14);
 
 		double thetaBins[] = new double[cb.nTheta+1];
 		for(int i = 0; i<cb.nTheta; i++){
 			thetaBins[i] = cb.thetaMin[i];
 		}
-		
+
 		thetaBins[thetaBins.length-1] = cb.thetaMax[cb.nTheta-1];
 
 		double phiBins[] = new double[315];
 		for(int i = 0; i<315; i++){
 			phiBins[i] = i/50.-3.14;  //every 10 mrad;
 		}
-		
+
+		double eBins[] = new double[66];
+		for(int i = 0; i<66; i++){
+			eBins[i] = i/50.;  //every 20 MeV up to 1300 MeV
+		}
+
+
+		h1 = hf.createHistogram2D("theta vs energy", "theta vs energy", thetaBins, eBins);
+
+
 		//identical to h2a, except different binning
 		h2b = hf.createHistogram2D("theta vs phi cut alt", "theta vs phi cut alt", thetaBins, phiBins);
 		h2c = hf.createHistogram2D("theta vs phi alt", "theta vs phi alt", thetaBins, phiBins);
-		
+
 		h3 = hf.createHistogram1D("theta", "theta", thetaBins);
-
-		h4 = hf.createHistogram2D("px\\/pz vs py\\/pz", 160, -.16, .24, 160, -.2, .2);
-		h4a = hf.createHistogram2D("px\\/pz vs py\\/pz cut", 160, -.16, .24, 160, -.2, .2);
-		
-		h5 = hf.createHistogram1D("energy", 75, 0, 1.5);
+		h3a = hf.createHistogram1D("theta isolated ", "theta isolated", thetaBins);
+
+		h3_t = hf.createHistogram1D("theta top", "theta top", thetaBins);
+		h3_b = hf.createHistogram1D("theta bottom", "theta bottom", thetaBins);
+
+		
+		h4 = hf.createHistogram2D("px\\/pz vs py\\/pz", 300, -.16, .24, 300, -.2, .2);
+		h4a = hf.createHistogram2D("px\\/pz vs py\\/pz cut", 300, -.16, .24, 300, -.2, .2);
+
+		h5 = hf.createHistogram1D("energy top", 75, 0, 1.5);
+		h5a = hf.createHistogram1D("energy bottom", 75, 0, 1.5);
+		
+
+		h9_t = hf.createHistogram1D("pz top", 75, 0, 1.5);
+		h9_b = hf.createHistogram1D("pz bottom", 75, 0, 1.5);
+		
+		h6 = hf.createHistogram2D("cluster", 47, -23.5, 23.5, 11, -5.5, 5.5);
+		h6a = hf.createHistogram2D("cluster matched", 47, -23.5, 23.5, 11, -5.5, 5.5);
+
+		h7 = hf.createHistogram1D("y top", 500, 0, 100);
+
+		h7a = hf.createHistogram1D("y bottom", 500, 0, 100);
+
+		h8 = hf.createHistogram1D("seed energy", 120, 0, 1.2);
+		
+		h10_t = hf.createHistogram1D("clustsize top", 10,0, 10);
+		h10_b = hf.createHistogram1D("clustsize bottom", 10,0, 10);
 	}
 	static void setupPlotter(IAnalysisFactory af){
 		IPlotterFactory pf = af.createPlotterFactory();
 		IPlotter p = pf.create();
 		p.createRegions(2,2);
 		p.region(0).plot(h2);
-
 		StyleUtil.stylize(p.region(0), "theta", "phi");
-		p.region(1).plot(h2a);
-		StyleUtil.stylize(p.region(1), "theta", "phi");
+		p.region(1).plot(h3a);
+		StyleUtil.stylize(p.region(1), "theta", "# of particles");
 		p.region(2).plot(h3);
 		StyleUtil.stylize(p.region(2), "theta", "# of particles");
 		p.region(3).plot(h5);
+		p.region(3).plot(h5a);
 		StyleUtil.stylize(p.region(3), "energy", "# of particles");
-		
+
 		p.show();
+		
 		//new window for the next plot
 		IPlotter p2 = pf.create();
 		p2.region(0).plot(h2b);
+		//IDataPointSetFactory dpsf = af.createDataPointSetFactory(af.createTreeFactory().create());
+
 		StyleUtil.stylize(p2.region(0), "theta", "phi");
-		
+
 		p2.show();
-		
+
 		//new window for the next plot
 		IPlotter p3 = pf.create();
 		p3.region(0).plot(h2c);
 		StyleUtil.stylize(p3.region(0), "theta", "phi");
-		
+
 		p3.show();
-		
+
 		//new window for the next plot
 		IPlotter p4 = pf.create();
 		p4.region(0).plot(h4);
 		StyleUtil.stylize(p4.region(0), "px/pz", "py/pz");
-		
+
 		p4.show();
-		
+
 		//new window for the next plot
 		IPlotter p5 = pf.create();
 		p5.region(0).plot(h4a);
 		StyleUtil.stylize(p5.region(0), "px/pz", "py/pz");
-		
+
 		p5.show();
+
+		IPlotter p6 = pf.create("efficiency");
+		p6.createRegions(1,2);
+		p6.region(0).plot(h6);
+		StyleUtil.stylize(p6.region(0), "ix", "iy");
+		p6.region(1).plot(h6a);
+		StyleUtil.stylize(p6.region(1), "ix", "iy");
+		p6.show();
+
+		IPlotter p7 = pf.create("theta vs energy");
+		//p6.createRegions(1,2);
+		p7.region(0).plot(h1);
+		StyleUtil.stylize(p7.region(0), "theta", "energy");
+		//		StyleUtil.stylize(p6.region(1), "ix", "iy");
+		p7.show();
+
+		IPlotter p8 = pf.create("y");
+		//p6.createRegions(1,2);
+		p8.region(0).plot(h7);
+		p8.region(0).plot(h7a);
+		StyleUtil.stylize(p8.region(0), "y", "# of particles");
+		//		StyleUtil.stylize(p6.region(1), "ix", "iy");
+		p8.show();
+		
+		IPlotter p9 = pf.create("theta: top vs. bottom");
+		//p6.createRegions(1,2);
+		p9.region(0).plot(h3_t);
+		p9.region(0).plot(h3_b);
+		StyleUtil.stylize(p9.region(0), "theta", "theta", "# of particles");
+		StyleUtil.noFillHistogramBars(p9.region(0));
+				StyleUtil.stylize(p6.region(1), "ix", "iy");
+		p9.show();
+		
+		IPlotter p10 = pf.create("seed energy");
+		//p6.createRegions(1,2);
+		p10.region(0).plot(h8);
+		StyleUtil.stylize(p10.region(0), "seed energy", "seed energy (GeV)", "# of particles");
+		//StyleUtil.noFillHistogramBars(p10.region(0));
+				StyleUtil.stylize(p6.region(1), "ix", "iy");
+		p10.show();
+
 	}
 	private static void processEvent(EventHeader event) {
-		if(event.getEventNumber() %1000 == 0)
+		if(event.getEventNumber() %10000 == 0)
 			System.out.println("event number " + event.getEventNumber());
-		
+
 		for (GenericObject gob : event.get(GenericObject.class,"TriggerBank"))
 		{
 			if (!(AbstractIntData.getTag(gob) == TIData.BANK_TAG)) continue;
@@ -194,11 +305,61 @@
 			}
 		}
 		List<ReconstructedParticle> particles = event.get(ReconstructedParticle.class, "FinalStateParticles");
-
-		for(ReconstructedParticle p : particles){
+		particles = RemoveDuplicateParticles.removeDuplicateParticles(particles);
+		outer : for(ReconstructedParticle p : particles){
+			//check if this is a duplicate particle (ie, different track same cluster)
+			/*for(ParticleInfo cand : candidates){
+				if(p.getEnergy() == cand.E && cand.isGood == true)
+					continue outer;
+			}*/
 
 			boolean isGood = addParticle(p);
-
+			if(!isGood){
+				if(p.getClusters().size()!= 0)
+					candidates.add(new ParticleInfo(0, p.getEnergy(), p.getClusters().get(0).getCalorimeterHits().get(0).getTime(), false));
+
+			}
+
+		}
+		for(ParticleInfo c : candidates){
+			if(c.isGood){
+				boolean isIsolated = true;
+				for(ParticleInfo c2 :candidates){
+					//try to remove events that have possible mollers in them. 
+					//correct for this later.  
+					if(Math.abs(c2.t - c.t) < 2.5 && c2.E < .3 && c2 != c){
+						isIsolated = false;
+						break;
+					}
+				}
+				if(isIsolated){
+					h3a.fill(c.theta);
+				}
+			}
+		}
+		candidates.clear();
+		processEventEfficiency(particles);
+	}
+
+	private static void processEventEfficiency(List<ReconstructedParticle> parts) {
+		for(ReconstructedParticle p : parts){
+			if(p.getClusters().size() == 0)
+				continue;
+			if(!(p.getEnergy() > eMin && p.getEnergy() < eMax))
+				continue;
+			Cluster c = p.getClusters().get(0);
+			if(!EcalUtil.fid_ECal(c) || c.getCalorimeterHits().size() < 3)
+				continue;
+			//good cluster?  now continue
+			int ixiy[] = EcalUtil.getCrystalIndex(c); 
+
+			h6.fill(ixiy[0], ixiy[1]);
+			if(p.getTracks().size() == 0)
+				continue;
+			Track t = p.getTracks().get(0);
+			if(t.getChi2()> maxChi2)
+				continue;
+			h6a.fill(ixiy[0], ixiy[1]);
 
 		}
 	}
@@ -209,44 +370,129 @@
 
 	static double beamTilt = .03057;
 	static double maxChi2 = 50;
+	//maximum difference between the reconstructed energy and momentum
+	static double maxdE = .3;
+
+	static double seedEnergyCut = .4;
+
+
+	static ArrayList<ParticleInfo> candidates = new ArrayList();
 	static boolean addParticle(ReconstructedParticle part){
-		
-		if(part.getTracks().size() == 0)
+
+
+
+		if(part.getCharge() != -1)
 			return false;
-		if(part.getTracks().get(0).getChi2()>maxChi2){
-			return false;
-		}
 		if(part.getClusters().size() == 0)
 			return false;
 		Cluster c = part.getClusters().get(0);
 		double time = c.getCalorimeterHits().get(0).getTime();
+
+		if(!(time>40 && time <50))
+			return false;
+		double seedEnergy = 0;
+		for(CalorimeterHit hit : c.getCalorimeterHits()){
+			if(hit.getCorrectedEnergy() > seedEnergy)
+				seedEnergy = hit.getCorrectedEnergy();
+		}
+		h8.fill(seedEnergy);
+		
+		
+		if(seedEnergy < seedEnergyCut)
+			return false;
+		
+		if(c.getPosition()[1] > 0){
+			h10_t.fill(c.getSize());
+		}
+		else{ 
+			h10_b.fill(c.getSize());
+		}
+		
+		
+		if(c.getCalorimeterHits().size() < 3)
+			return false;
+		
+
+		if(c.getEnergy() > eMin && c.getEnergy() < eMax){
+			if(c.getPosition()[1] > 0)
+				h7.fill(c.getPosition()[1]);
+			else if(c.getPosition()[1] < 0)
+				h7a.fill(-c.getPosition()[1]);
+		}
+		
 		if(EcalUtil.fid_ECal(c)){
-			if(c.getCalorimeterHits().size() < 3)
+			
+			if(part.getTracks().size() == 0)
 				return false;
-			if(time>40 && time <48)
+			Track t = part.getTracks().get(0);
+			if(t.getChi2()>maxChi2){
+				return false;
+			}
+			if(!TrackType.isGBL(t.getType()))
+				return false;
+			if(c.getPosition()[1] > 0){
 				h5.fill(c.getEnergy());
-			if(c.getEnergy() > eMin && c.getEnergy() < eMax && (time >40 && time < 48)) {
-
-				Hep3Vector p = part.getMomentum();
-				
-				double px = p.x(), pz = p.z();
-				double pxtilt = px*Math.cos(beamTilt)-pz*Math.sin(beamTilt);
-				double py = p.y();
-				double pztilt = pz*Math.cos(beamTilt)+px*Math.sin(beamTilt);
-
-				double theta = Math.atan(Math.hypot(pxtilt, py)/pztilt);
-				double phi =Math.atan2(py, pxtilt);
+			}
+			else{ 
+				h5a.fill(c.getEnergy());
+			}
+			
+			
+			Hep3Vector p = part.getMomentum();
+
+
+
+			double px = p.x(), pz = p.z();
+			double pxtilt = px*Math.cos(beamTilt)-pz*Math.sin(beamTilt);
+			double py = p.y();
+			double pztilt = pz*Math.cos(beamTilt)+px*Math.sin(beamTilt);
+
+			if(Math.abs(pztilt - c.getEnergy()) > maxdE)
+				return false;
+			if(c.getPosition()[1] > 0)
+				h9_t.fill(pztilt);
+			else
+				h9_b.fill(pztilt);
+
+			double theta = Math.atan(Math.hypot(pxtilt, py)/pztilt);
+			double phi =Math.atan2(py, pxtilt);
+			boolean inRange = cb.inRange(theta, phi);
+			if(inRange)
+				h1.fill(theta, c.getEnergy());
+
+
+
+			if(c.getEnergy() > eMin && c.getEnergy() < eMax) {
+
+
 
 				h2.fill(theta, phi);
 				h2c.fill(theta, phi);
-				
+
 				h4.fill(px/pz, py/pz);
-				
-				if(cb.inRange(theta, phi)){
+
+				if(inRange){
+
+					//System.out.println(c.getEnergy() + " " + t.getType());
+					/*for(TrackState ts : t.getTrackStates()){
+						if(ts.getLocation() == TrackState.AtIP)
+							System.out.println(Arrays.toString( 
+									ts.getReferencePoint()));
+					}*/
 					h2a.fill(theta, phi);
 					h2b.fill(theta, phi);
+					
 					h3.fill(theta);
+					if(py > 0)
+						h3_t.fill(theta);
+					else 
+						h3_b.fill(theta);
+					//if(h3_t.sumBinHeights()+h3_b.sumBinHeights() != h3.sumBinHeights())
+						//System.out.println("NABO ERROR");
+					
+					
 					h4a.fill(px/pz, py/pz);
+					candidates.add(new ParticleInfo(theta, c.getEnergy(), c.getCalorimeterHits().get(0).getTime(), true));
 				}
 
 
@@ -256,5 +502,16 @@
 		}
 		return false;
 	}
-
+	static class ParticleInfo{
+		double theta;
+		double E;
+		double t;
+		boolean isGood;
+		ParticleInfo(double theta, double E, double t, boolean isGood){
+			this.theta = theta;
+			this.E = E;
+			this.t = t;
+			this.isGood = isGood;
+		}
+	}
 }

Added: java/trunk/users/src/main/java/org/hps/users/spaul/feecc/MottIntegral.java
 =============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/feecc/MottIntegral.java	(added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/feecc/MottIntegral.java	Sun Jan  3 10:18:03 2016
@@ -0,0 +1,66 @@
+package org.hps.users.spaul.feecc;
+
+public class MottIntegral {
+	/**
+	 * @param bin
+	 * @param scale = Z^2alpha^2/(4*E^2)
+	 * @param recoil = 2E/M
+	 * @return the integral of 1/sin^4(th/2)*cos^2(th/2)/(1+a*sin^2(th/2)) times dPhi*sin(theta),
+	 * which appears in the integral of mott scattering.
+	 * 
+	 * NOTE the sin(theta) is because dOmega = dTheta*dPhi*sin(theta)
+	 */
+	static double mottIntegral(double recoil, double scale, int bin, CustomBinning cb){
+		double dPhi = 0;
+		for(int i = 0; i< cb.phiMax[bin].length; i++)
+			dPhi += 2*(cb.phiMax[bin][i] - cb.phiMin[bin][i]); //factor of 2 from top and bottom
+		
+		double Imax = integral(recoil, cb.thetaMax[bin]);
+		       
+		double Imin = integral(recoil, cb.thetaMin[bin]);
+		double dI = Imax-Imin; 
+		
+		
+		double retval = scale*dPhi*dI;
+		return retval;
+	}
+	
+	static double mottIntegral(double recoil, double scale, double thetaMin, double thetaMax){
+		double Imax = integral(recoil, thetaMax);
+	
+		double Imin = integral(recoil, thetaMin);
+		double dI = Imax-Imin; 
+		
+		double dPhi = 2*Math.PI;  //full range in phi
+		double retval = scale*dPhi*dI;
+		return retval;
+	}
+	
+	static double integral(double a, double th){
+		double sinth22 = Math.pow(Math.sin(th/2), 2);
+		return 2*(-1/sinth22+(1+a)*(Math.log(2/sinth22+2*a)));
+	}
+	
+	static double mottIntegralWithFormFactor(double recoil, double scale, double thetaMin, double thetaMax, FormFactor ff, int nsteps, double E){
+		double sum = 0;
+		double prev = -1;
+
+		double dTheta = (thetaMax-thetaMin)/nsteps;
+		for(int i = 0; i<nsteps; i++){
+			double theta = i*(thetaMax-thetaMin)/nsteps+thetaMin;
+			double I = integral(recoil, theta);
+			if(i != 0)
+			{
+				double ts = Math.pow(Math.sin(theta/2),2);
+				double q2 = 4*E*E*ts/(1+recoil*ts);
+				double f2 = ff.getFormFactorSquared(q2);
+				sum+= (I-prev)*f2;
+			}
+			prev = I;
+			
+		}
+		double dPhi = 2*Math.PI;  //full range in phi
+		double retval = scale*dPhi*sum;
+		return retval;
+	}
+}

Added: java/trunk/users/src/main/java/org/hps/users/spaul/feecc/MultipleScattering.java
 =============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/feecc/MultipleScattering.java	(added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/feecc/MultipleScattering.java	Sun Jan  3 10:18:03 2016
@@ -0,0 +1,106 @@
+package org.hps.users.spaul.feecc;
+
+import hep.aida.IAnalysisFactory;
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogramFactory;
+import hep.aida.IPlotter;
+
+import java.util.Random;
+
+public class MultipleScattering {
+	
+	
+	/* 
+	 V.L. Highland, Nucl. Instrum. Methods 129, 497 (1975); Nucl. Instrum. Methods 161, 171 (1979)
+	 
+	 G.R. Lynch and O.I. Dahl, Nucl. Instrum. Methods B58, 6 (1991).
+	 
+	 */
+	/**
+	 * @param t number of scattering lengths
+	 * @param p momentum in GeV/c
+	 * @return the rms scattering angle of an electron with beta == 1.
+	 */
+	static double getThetaRMS(double t, double p){
+		return .0136/p*Math.sqrt(t)*(1+.038*Math.log(t));
+	}
+	/**
+	 * Mo Tsai equations A.4-5
+	 * @param Z
+	 * @return
+	 */
+	static double b(int Z){
+		double xi = Math.log(1440.*Math.pow(Z, -2/3.))/Math.log(183.*Math.pow(Z,-1/3.));
+		System.out.println(xi);
+		return 4/3.*(1+1/9.*(Z+1.)/(Z+xi)/Math.log(183.*Math.pow(Z,-1/3.)));
+		
+	
+	}
+	
+	
+	public static void main(String arg[]){
+		
+		System.out.println(b(6));
+		System.out.println(b(74));
+		
+		double bt = 0.03;
+		
+		IAnalysisFactory af = IAnalysisFactory.create();
+		IHistogramFactory hf = af.createHistogramFactory(af.createTreeFactory().create());
+		IHistogram1D h = hf.createHistogram1D("x", 100, 0, 1.5);
+		IPlotter p = af.createPlotterFactory().create();
+		p.region(0).plot(h);
+		p.show();
+		for(int i = 0; i< 100000; i++){
+			h.fill(getRandomEnergyFraction(bt)*getRandomEnergyFraction(bt)*(1+r.nextGaussian()*.045)*1.056);
+		}
+	}
+	
+	static double pdf(double x, double bt){
+		return bt/(1-x)*(x+3/4.*(1-x)*(1-x))*Math.pow(Math.log(1./x), bt);
+	}
+	static Random r = new Random();
+	/**************************
+	 * CONFIGURABLE VARIABLES *
+	 **************************/
+	/**
+	 * step size in x
+	 */
+	static double dx = .001;
+	/**
+	 * starting distance from peak to begin numerical integration
+	 */
+	static double delta0 = .05;
+	
+	/**
+	 * fraction of the energy to cut off at.  (all events with less than this energy will be blanked out,
+	 * and therefore not simulated in slic.  Slic was having problems dealing with low energy events.  
+	 */
+	static double xCutoff = .5;
+	
+	static double getRandomEnergyFraction(double bt){
+		double p = r.nextDouble(); 
+		if(p< Math.pow(delta0, bt)){
+			return 1-Math.pow(p, 1/bt);
+		}
+		double sum = Math.pow(delta0,bt);
+		if(sum > p){
+			return 1;
+		}
+		double x = 1-delta0;
+		while(true){
+			double contribution = pdf(x, bt)*dx;
+			if(sum + contribution > p)
+			{
+				return x - (p-sum)/contribution*dx;
+			}
+			if(x < 0)
+			{
+				return 0;
+			}
+			x -= dx;
+			sum += contribution;
+		}
+		//return 1;
+	}
+}

Modified: java/trunk/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinning.java
 =============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinning.java	(original)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinning.java	Sun Jan  3 10:18:03 2016
@@ -2,10 +2,16 @@
 
 import java.awt.Canvas;
 import java.awt.Color;
+import java.awt.Container;
+import java.awt.Font;
 import java.awt.Graphics;
+import java.awt.Graphics2D;
+import java.awt.image.BufferedImage;
 import java.io.File;
 import java.io.FileNotFoundException;
-
+import java.io.IOException;
+
+import javax.imageio.ImageIO;
 import javax.swing.JFrame;
 
 import hep.aida.IAnalysisFactory;
@@ -27,15 +33,47 @@
 	public static void main(String arg[]) throws FileNotFoundException{
 		JFrame frame = new JFrame();
 		frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
-		frame.add(new ShowCustomBinning(new File(arg[0])));
-		frame.setSize(800, 800);
+		Canvas c = new ShowCustomBinning(new File(arg[0]));
+		String outdir = arg[1];
+		frame.add(c);
+		frame.setSize(1200, 800);
 		frame.setVisible(true);
 		
+		
+		try {
+			BufferedImage im = new BufferedImage(c.getWidth(), c.getHeight(), BufferedImage.TYPE_INT_ARGB);
+			c.paint(im.getGraphics());
+			ImageIO.write(im, "PNG", new File(outdir +"/bins.png"));
+		} catch (IOException e) {
+			// TODO Auto-generated catch block
+			e.printStackTrace();
+		}
+		
+		frame = new JFrame();
+		frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
+		c = new ShowCustomBinningXY(new File(arg[0]));
+		frame.add(c);
+		frame.setSize(1200, 615);
+		frame.setVisible(true);
+		
+		try {
+			BufferedImage im = new BufferedImage(c.getWidth(), c.getHeight(), BufferedImage.TYPE_INT_ARGB);
+			c.paint(im.getGraphics());
+			ImageIO.write(im, "PNG", new File(outdir + "/bins_xy.png"));
+		} catch (IOException e) {
+			// TODO Auto-generated catch block
+			e.printStackTrace();
+		}
 	}
 	public void paint(Graphics g){
+		g.setFont(new Font(Font.DIALOG, Font.PLAIN, 24));
+		
 		drawEcalOutline(g);
 		drawFidEcalOutline(g);
 		drawCustomBinRectangles(g);
+		g.setColor(Color.BLACK);
+		drawXAxis(g);
+		drawYAxis(g);
 	}
 	
 	void drawFidEcalOutline(Graphics g){
@@ -95,15 +133,24 @@
 	}
 	
 	void drawEcalOutline(Graphics g){
+		/*double x_edge_low = -262.74;
+		double x_edge_high = 347.7;
+		double y_edge_low = 33.54;
+		double y_edge_high = 75.18;
+
+		double x_gap_low = -106.66;
+		double x_gap_high = 42.17;
+		double y_gap_high = 47.18;*/
+		
 		g.setColor(Color.BLACK);
-		double x_edge_low = -276.50;
-		double x_edge_high = 361.55;
-		double y_edge_low = 20.17;
-		double y_edge_high = 89;
-
-		double x_gap_low = -93.30;
-		double x_gap_high = 28.93;
-		double y_gap_high = 33.12;
+		double x_edge_low = -269.56;
+		double x_edge_high = 354.52;
+		double y_edge_low = 26.72;
+		double y_edge_high = 82;
+
+		double x_gap_low = -99.84;
+		double x_gap_high = 33.35;
+		double y_gap_high = 40.36;
 		double x1,y1, x2, y2;
 		double nPoints = 500;
 		for(int i = 0; i< nPoints-1; i++){
@@ -157,6 +204,8 @@
 	}
 	Color altBin1 = new Color(0, 0, 128);
 	Color altBin2 = new Color(0,128,0);
+	Color fillBin1 = new Color(196, 196, 255);
+	Color fillBin2 = new Color(196,255,196);
 	void drawCustomBinRectangles(Graphics g){
 		for(int i = 0; i<binning.nTheta; i++){
 			g.setColor(i%2 == 0 ? altBin1 : altBin2);
@@ -167,13 +216,45 @@
 				double theta2 = binning.thetaMax[i];
 				
 				int x =getX(theta1)+1, y = getY(phi1), w =  getX(theta2)-getX(theta1), h = getY(phi2)-getY(phi1);
+
+				g.setColor(i%2 == 0 ? fillBin1 : fillBin2);
+				g.fillRect(x, y, w, h);
+				g.setColor(i%2 == 0 ? altBin1 : altBin2);
 				g.drawRect(x, y, w, h);
 				 x =getX(theta1)+1; y = getY(-phi2); w =  getX(theta2)-getX(theta1); h = getY(-phi1)-getY(-phi2);
-
-				g.drawRect(x, y, w, h);
-				
-				
+				 g.setColor(i%2 == 0 ? fillBin1 : fillBin2);
+					g.fillRect(x, y, w, h);
+					g.setColor(i%2 == 0 ? altBin1 : altBin2);
+					g.drawRect(x, y, w, h);
 			}
+			
+		}
+	}
+	void drawXAxis(Graphics g){
+		//x axis
+		g.drawString("θ", getX(.28), getY(0) - 40);
+		g.drawLine(getX(0), getY(0), getX(.28), getY(0));
+		for(int i = 0; i< 28; i++){
+			if(i%5 == 0 && i != 0){
+				g.drawString(i/100.+"", getX(i/100.)-15, getY(0)-20);
+				g.drawLine(getX(i/100.), getY(0), getX(i/100.), getY(0)-15);
+			}
+			g.drawLine(getX(i/100.), getY(0), getX(i/100.), getY(0)-5);
+		}
+	}
+	void drawYAxis(Graphics g){
+		g.drawString("φ", getX(0)+70, getY(3));
+		g.drawLine(getX(0), getY(-3), getX(0), getY(3));
+		for(int i = -30; i<= 30; i++){
+			if(i%5 == 0 && i != 0){
+				g.drawString(i/10.+"", getX(0)+20, getY(i/10.) + 5);
+
+				g.drawLine(getX(0), getY(i/10.), getX(0) + 15, getY(i/10.));
+			}
+			if(i == 0){
+				//g.drawString(i/10.+"", getX(0)+10, getY(i/10.) - 15);
+			}
+			g.drawLine(getX(0), getY(i/10.), getX(0) + 5, getY(i/10.));
 		}
 	}
 	
@@ -186,8 +267,9 @@
 		
 	}
 	int getX(double theta){
-		return (int)(this.getWidth()*theta/.3);
-	}
+		return (int)(this.getWidth()*theta/.3)+left_margin;
+	}
+	int left_margin = 20;
 	int getY(double phi){
 		return (int)(this.getHeight()*(3.2-phi)/6.4);
 	}

Added: java/trunk/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinningXY.java
 =============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinningXY.java	(added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/feecc/ShowCustomBinningXY.java	Sun Jan  3 10:18:03 2016
@@ -0,0 +1,227 @@
+package org.hps.users.spaul.feecc;
+
+import java.awt.Canvas;
+import java.awt.Color;
+import java.awt.Graphics;
+import java.awt.Polygon;
+import java.io.File;
+import java.io.FileNotFoundException;
+
+public class ShowCustomBinningXY extends ShowCustomBinning{
+	ShowCustomBinningXY(File file) throws FileNotFoundException
+	{
+		super(file);
+	}
+	@Override
+	void drawCustomBinRectangles(Graphics g){
+		for(int i = 0; i<binning.nTheta; i++){
+			g.setColor(i%2 == 0 ? altBin1 : altBin2);
+			for(int j = 0; j<binning.phiMax[i].length; j++){
+				double phi1 = binning.phiMax[i][j];
+				double phi2 = binning.phiMin[i][j];
+				double theta1 = binning.thetaMin[i];
+				double theta2 = binning.thetaMax[i];
+
+				double nK = 10;
+				for(int sign = -1; sign <=1; sign +=2){
+					for(int k = 0; k< nK; k++){
+						drawLineFromPolar(theta1, sign*(phi1+(phi2-phi1)*k/nK), theta1, sign*(phi1+(phi2-phi1)*(k+1)/nK),g);
+					}
+					for(int k = 0; k< nK; k++){
+						drawLineFromPolar(theta1 + (theta2-theta1)*k/nK, sign*phi2, theta1 + (theta2-theta1)*(k+1)/nK, sign*phi2,g);
+					}
+					for(int k = 0; k< nK; k++){
+						drawLineFromPolar(theta2, sign*(phi2+(phi1-phi2)*k/nK), theta2, sign*(phi2+(phi1-phi2)*(k+1)/nK),g);
+					}
+					for(int k = 0; k< nK; k++){
+						drawLineFromPolar(theta2 + (theta1-theta2)*k/nK, sign*phi1, theta2 + (theta1-theta2)*(k+1)/nK, sign*phi1,g);
+					}
+					
+					closePolarFigure(g, i%2 == 0 ? altBin1 : altBin2, i%2 == 0 ? fillBin1 : fillBin2);
+				}
+
+			}
+		}
+	}
+	Polygon p = null;
+	private void drawLineFromPolar(double theta1, double phi1, double theta2,
+			double phi2, Graphics g) {
+		double[] xy1 = toXY(theta1, phi1);
+		double[] xy2 = toXY(theta2, phi2);
+
+		if(p == null)
+			p = new Polygon();
+		/*g.drawLine(getX(xy1[0]), 
+				getY(xy1[1]), 
+				getX(xy2[0]),
+				getY(xy2[1]));*/
+		p.addPoint(getX(xy2[0]), getY(xy2[1]));
+	}
+	private void closePolarFigure(Graphics g, Color outlineColor, Color fillColor){
+		g.setColor(fillColor);
+		g.fillPolygon(p);
+		g.setColor(outlineColor);
+		g.drawPolygon(p);
+		p = null;
+		
+	}
+	double tilt = .03057;
+	double[] toXY(double theta, double phi){
+		double ux = Math.cos(phi)*Math.sin(theta)*Math.cos(tilt)+Math.cos(theta)*Math.sin(tilt);
+		double uy = Math.sin(phi)*Math.sin(theta);
+		double uz = Math.cos(theta)*Math.cos(tilt)-Math.cos(phi)*Math.sin(theta)*Math.sin(tilt);
+		double pxpz = ux/uz;
+		double pypz = uy/uz;
+		return new double[]{pxpz, pypz};
+	}
+	int getX(double pxpz){
+		//return (int)((pxpz -(-.16))/(.34-(-.16))*getWidth());
+		return (int)((.32-pxpz)/(.32-(-.18))*getWidth()); 
+	}
+	int getY(double pypz){
+		return getHeight()-(int)((pypz -(-.125))/(.125-(-.125))*getHeight()); 
+	}
+	void drawFidEcalOutline(Graphics g){
+		g.setColor(Color.GRAY);
+		double x_edge_low = -262.74;
+		double x_edge_high = 347.7;
+		double y_edge_low = 33.54;
+		double y_edge_high = 75.18;
+
+		double x_gap_low = -106.66;
+		double x_gap_high = 42.17;
+		double y_gap_high = 47.18;
+		double x1,y1, x2, y2;
+		double nPoints = 500;
+		for(int i = 0; i< nPoints-1; i++){
+			x1 = x_gap_high+i/nPoints*(x_edge_high-x_gap_high);
+			x2 = x_gap_high+(i+1)/nPoints*(x_edge_high-x_gap_high);
+			y1 = y_edge_low;
+			y2 = y1;
+			drawEcalFaceLine(g, x1, y1, x2, y2);
+			drawEcalFaceLine(g, x1, -y1, x2, -y2);
+			
+			
+			x1 = x_edge_low+i/nPoints*(x_gap_low-x_edge_low);
+			x2 = x_edge_low+(i+1)/nPoints*(x_gap_low-x_edge_low);
+			y1 = y2 = y_edge_low;
+			drawEcalFaceLine(g, x1, y1, x2, y2);
+			drawEcalFaceLine(g, x1, -y1, x2, -y2);
+			
+			
+			
+			x1 = x_gap_low+i/nPoints*(x_gap_high-x_gap_low);
+			x2 = x_gap_low+(i+1)/nPoints*(x_gap_high-x_gap_low);
+			y1 = y2 = y_gap_high;
+			drawEcalFaceLine(g, x1, y1, x2, y2);
+			drawEcalFaceLine(g, x1, -y1, x2, -y2);
+			
+			x1 = x_edge_low+i/nPoints*(x_edge_high-x_edge_low);
+			x2 = x_edge_low+(i+1)/nPoints*(x_edge_high-x_edge_low);
+			y1 = y2 = y_edge_high;
+			drawEcalFaceLine(g, x1, y1, x2, y2);
+			drawEcalFaceLine(g, x1, -y1, x2, -y2);
+		}
+		drawEcalFaceLine(g, x_gap_low, y_edge_low, x_gap_low, y_gap_high);
+		drawEcalFaceLine(g, x_gap_high, y_edge_low, x_gap_high, y_gap_high);
+
+		drawEcalFaceLine(g, x_edge_low, y_edge_low, x_edge_low, y_edge_high);
+		drawEcalFaceLine(g, x_edge_high, y_edge_low, x_edge_high, y_edge_high);
+		
+
+		drawEcalFaceLine(g, x_gap_low, -y_edge_low, x_gap_low, -y_gap_high);
+		drawEcalFaceLine(g, x_gap_high, -y_edge_low, x_gap_high, -y_gap_high);
+
+		drawEcalFaceLine(g, x_edge_low, -y_edge_low, x_edge_low, -y_edge_high);
+		drawEcalFaceLine(g, x_edge_high, -y_edge_low, x_edge_high, -y_edge_high);
+	
+	}
+	
+	void drawEcalOutline(Graphics g){
+		g.setColor(Color.BLACK);
+		double x_edge_low = -269.56;
+		double x_edge_high = 354.52;
+		double y_edge_low = 26.72;
+		double y_edge_high = 82;
+
+		double x_gap_low = -99.84;
+		double x_gap_high = 33.35;
+		double y_gap_high = 40.36;
+		double x1,y1, x2, y2;
+			x1 = x_gap_high;
+			x2 = x_edge_high;
+			y1 = y_edge_low;
+			y2 = y1;
+			drawEcalFaceLine(g, x1, y1, x2, y2);
+			drawEcalFaceLine(g, x1, -y1, x2, -y2);
+			
+			
+			x1 = x_edge_low;
+			x2 = x_gap_low;
+			y1 = y2 = y_edge_low;
+			drawEcalFaceLine(g, x1, y1, x2, y2);
+			drawEcalFaceLine(g, x1, -y1, x2, -y2);
+			
+			
+			
+			x1 = x_gap_low;
+			x2 = x_gap_high;
+			y1 = y2 = y_gap_high;
+			drawEcalFaceLine(g, x1, y1, x2, y2);
+			drawEcalFaceLine(g, x1, -y1, x2, -y2);
+			
+			x1 = x_edge_low;
+			x2 = x_edge_high;
+			y1 = y2 = y_edge_high;
+			drawEcalFaceLine(g, x1, y1, x2, y2);
+			drawEcalFaceLine(g, x1, -y1, x2, -y2);
+		
+		drawEcalFaceLine(g, x_gap_low, y_edge_low, x_gap_low, y_gap_high);
+		drawEcalFaceLine(g, x_gap_high, y_edge_low, x_gap_high, y_gap_high);
+
+		drawEcalFaceLine(g, x_edge_low, y_edge_low, x_edge_low, y_edge_high);
+		drawEcalFaceLine(g, x_edge_high, y_edge_low, x_edge_high, y_edge_high);
+		
+
+		drawEcalFaceLine(g, x_gap_low, -y_edge_low, x_gap_low, -y_gap_high);
+		drawEcalFaceLine(g, x_gap_high, -y_edge_low, x_gap_high, -y_gap_high);
+
+		drawEcalFaceLine(g, x_edge_low, -y_edge_low, x_edge_low, -y_edge_high);
+		drawEcalFaceLine(g, x_edge_high, -y_edge_low, x_edge_high, -y_edge_high);
+	}
+	void drawEcalFaceLine(Graphics g, double x1, double y1, double x2, double y2){
+
+		double[] h1 = EcalUtil.toHollysCoordinates(x1, y1, 1.057, 11);
+		double[] h2 =  EcalUtil.toHollysCoordinates(x2, y2, 1.057, 11);
+		g.drawLine(getX(Math.tan(h1[1])), getY(Math.tan(h1[0])), getX(Math.tan(h2[1])), getY(Math.tan(h2[0])));
+		
+	}
+	
+	void drawXAxis(Graphics g){
+		//x axis
+		g.drawString("px/pz", getX(.30)-10, getY(0) - 45);
+		g.drawLine(getX(-.15), getY(0), getX(.30), getY(0));
+		for(int i = -15; i<= 30; i++){
+			if(i%5 == 0 && i != 0){
+				g.drawString(i/100.+"", getX(i/100.)-15, getY(0)-20);
+				g.drawLine(getX(i/100.), getY(0), getX(i/100.), getY(0)-15);
+			}
+			g.drawLine(getX(i/100.), getY(0), getX(i/100.), getY(0)-5);
+		}
+	}
+	void drawYAxis(Graphics g){
+		g.drawString("py/pz", getX(0), getY(.1)-20);
+		g.drawLine(getX(0), getY(-.1), getX(0), getY(.1));
+		for(int i = -10; i<= 10; i++){
+			if(i%5 == 0 && i != 0){
+				g.drawString(i/100.+"", getX(0)+15, getY(i/100.) + 5);
+
+				g.drawLine(getX(0), getY(i/100.), getX(0) + 10, getY(i/100.));
+			}
+			if(i == 0){
+				//g.drawString(i/10.+"", getX(0)+10, getY(i/10.) - 15);
+			}
+			g.drawLine(getX(0), getY(i/100.), getX(0) + 5, getY(i/100.));
+		}
+	}
+}

Top of Message | Previous Page | Permalink

Advanced Options


Options

Log In

Log In

Get Password

Get Password


Search Archives

Search Archives


Subscribe or Unsubscribe

Subscribe or Unsubscribe


Archives

November 2017
August 2017
July 2017
January 2017
December 2016
November 2016
October 2016
September 2016
August 2016
July 2016
June 2016
May 2016
April 2016
March 2016
February 2016
January 2016
December 2015
November 2015
October 2015
September 2015
August 2015
July 2015
June 2015
May 2015
April 2015
March 2015
February 2015
January 2015
December 2014
November 2014
October 2014
September 2014
August 2014
July 2014
June 2014
May 2014
April 2014
March 2014
February 2014
January 2014
December 2013
November 2013

ATOM RSS1 RSS2



LISTSERV.SLAC.STANFORD.EDU

Secured by F-Secure Anti-Virus CataList Email List Search Powered by the LISTSERV Email List Manager

Privacy Notice, Security Notice and Terms of Use