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.));
+ }
+ }
+}
