5 modified files
lcsim/src/org/lcsim/contrib/CarstenHensel
diff -u -r1.3 -r1.4
--- EMClusterIDAnalyzer.java 24 Aug 2005 15:11:51 -0000 1.3
+++ EMClusterIDAnalyzer.java 25 Aug 2005 19:09:36 -0000 1.4
@@ -47,23 +47,27 @@
// Constructor: Initialize the class. This is the only constructor, and
// it requires no inputs.
//
+
+
+
public EMClusterIDAnalyzer() {
this.tree = aida.tree();
this.df = new DecimalFormat();
this.angleVSk = aida.cloud2D("angle vs k");
-// add sub-drivers
+ // add sub-drivers
FixedConeClustererAnalyzer fixedConeClustererAnalyzer = new FixedConeClustererAnalyzer(this.tree);
fixedConeClustererAnalyzer.setAttributes(this.SAMP_FRAC, this.debug, this.interactionRadius);
- this.add(fixedConeClustererAnalyzer);
+ //this.add(fixedConeClustererAnalyzer);
- HMatrixAnalyzer hMatrixAnalyzer = new HMatrixAnalyzer("/home/carsten/LC/Output/CH_cdcaug05_gamma_Theta90_5GeV.hmx", tree);
+ HMatrixAnalyzer hMatrixAnalyzer = new HMatrixAnalyzer("/home/carsten/LC/Output/cdcaug05_gamma_Theta90_5GeV.hmx", tree);
//HMatrixAnalyzer hMatrixAnalyzer = new HMatrixAnalyzer(tree);
this.add(hMatrixAnalyzer);
- MyHMatrixBuilder builder = new MyHMatrixBuilder();
- this.add(builder);
+ int debug = 0;
+ MyHMatrixBuilder builder = new MyHMatrixBuilder(this.debug);
+ //this.add(builder);
}
lcsim/src/org/lcsim/contrib/CarstenHensel
diff -u -r1.1 -r1.2
--- HMatrixAnalyzer.java 24 Aug 2005 15:09:26 -0000 1.1
+++ HMatrixAnalyzer.java 25 Aug 2005 19:09:37 -0000 1.2
@@ -14,6 +14,7 @@
import hep.aida.ITuple;
import hep.aida.ITupleFactory;
import hep.aida.ref.tree.TreeObjectAlreadyExistException;
+import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
@@ -55,6 +56,10 @@
private boolean attributesSet = false;
private boolean initialized = false;
private ITuple tuple;
+ private String detectorName;
+ private int debug;
+ private double samplingFraction;
+ private static final double log10inv = 1./Math.log(10.0);
/** Creates a new instance of HMatrixAnalyzer */
public HMatrixAnalyzer(ITree tree) {
@@ -73,15 +78,38 @@
}
public void initialize(EventHeader event){
+
+ this.detectorName = event.getDetectorName();
+ if (this.debug > 0)
+ System.out.println("detector is " + this.detectorName);
+
+ // how to get the sampling fraction from the conditions data base?
+ if (this.detectorName.equals("cdcaug05")) {
+ this.samplingFraction = 0.0113;
+ } else {
+ this.samplingFraction = 0.0117;
+ }
+
+ if (this.debug > 0)
+ System.out.println("Sampling Fraction set to " + this.samplingFraction);
+
+ CylindricalCalorimeter calsub = (CylindricalCalorimeter)event.getDetector().getSubdetectors().get("EMBarrel");
+ this.nLayers = calsub.getLayering().getLayerCount();
+ System.out.println("found " + this.nLayers + " layers in the EMBarrel");
+
+
+
+
this.nLayers = ((CylindricalCalorimeter)event.getDetector().getSubdetectors().get("EMBarrel")).getLayering().getLayerCount();
this.tree.cd("/");
this.tree.mkdir(this.treePath);
IAnalysisFactory af = IAnalysisFactory.create();
ITupleFactory tf = af.createTupleFactory(tree);
- String[] columnNames = { "chisq", "chisqDiag", "prob", "clusterEnergy" };
- Class[] columnClasses = { Double.TYPE, Double.TYPE, Double.TYPE, Double.TYPE };
- this.tuple = tf.create( this.treePath + "MyTuple", "myTupleLabel", columnNames, columnClasses, "");
+ String[] columnNames = { "chisq", "chisqDiag", "prob", "logProb", "probDiag", "logProbDiag", "clusterEnergy", "energyPerLayer"};
+
+ Class[] columnClasses = { Double.TYPE, Double.TYPE, Double.TYPE, Double.TYPE, Double.TYPE , Double.TYPE, Double.TYPE, Object.class};
+ this.tuple = tf.create( this.treePath + "HMatrixTuple", "ThisIsALabel", columnNames, columnClasses, "");
this.initialized = true;
}
@@ -142,20 +170,36 @@
double[] val = new double[this.nLayers + 1];
double[] energiesInLayer = this.getListOfLayerEnergies(cluster);
+ ArrayList<Double> energyPerLayer = new ArrayList<Double>();
for (int i = 0; i < energiesInLayer.length; i++) {
val[i] = energiesInLayer[i];
aida.cloud2D("Fractional Energy vs Layer").fill(i, energiesInLayer[i]);
+ energyPerLayer.add(new Double(energiesInLayer[i]));
}
-
double clusterRawEnergy = cluster.getRawEnergy();
aida.cloud1D("cluster raw energy").fill(clusterRawEnergy);
val[this.nLayers] = Math.log(clusterRawEnergy) / Math.log(10.0);;
+ val = this.getHMInputVals(cluster);
+
double chisq = this.hMatrix.chisquared(val);
double chisqDiag = this.hMatrix.chisquaredDiagonal(val);
double prob = ChisqProb.gammq(this.nLayers + 1, chisq);
+ double logProb;
+ if (prob > 0.0) {
+ logProb = Math.log10(prob);
+ } else {
+ logProb = -999.0;
+ }
double probDiag = ChisqProb.gammq(this.nLayers + 1, chisqDiag);
+ double logProbDiag;
+ if (probDiag > 0) {
+ logProbDiag = Math.log10(probDiag);
+ } else {
+ logProbDiag = -.0;
+ }
+
aida.cloud1D("chisq").fill(chisq);
aida.cloud1D("chisq diag").fill(chisqDiag);
aida.cloud1D("chisq probability").fill(prob);
@@ -166,7 +210,11 @@
this.tuple.fill(0, chisq);
this.tuple.fill(1, chisqDiag);
this.tuple.fill(2, prob);
- this.tuple.fill(3, clusterRawEnergy);
+ this.tuple.fill(3, logProb);
+ this.tuple.fill(4, probDiag);
+ this.tuple.fill(5, logProbDiag);
+ this.tuple.fill(6, clusterRawEnergy);
+ this.tuple.fill(7, (Object)energyPerLayer);
this.tuple.addRow();
}
@@ -244,6 +292,52 @@
}
}
+ private double[] getHMInputVals(BasicCluster clus) {
+ //FIXME could reuse this array
+ double[] layerEnergies = new double[this.nLayers + 1];
+ double clusterEnergy = 0.;
+ List<CalorimeterHit> hits = clus.getCalorimeterHits();
+ //System.out.println("New cluster with "+hits.size()+ " hits"); // and energy "+clus.getEnergy());
+ for(CalorimeterHit hit : hits) {
+ this.decoder.setID(hit.getCellID());
+ // GWW double e = hit.getEnergy();
+ double e = hit.getRawEnergy();
+ int layer = this.decoder.getLayer();
+
+
+ if ((this.detectorName.equals("cdcaug05")) && (layer > 20)) {
+ e = e / (this.samplingFraction * 0.5);
+ } else {
+ e = e / this.samplingFraction;
+ }
+
+ if (this.debug > 0)
+ System.out.println("layer "+layer+" energy "+e);
+
+ clusterEnergy+=e;
+ layerEnergies[layer]+=e;
+ }
+
+ if (this.debug > 0)
+ System.out.println("Cluster energy = " + clusterEnergy + "\n");
+
+ layerEnergies[this.nLayers] = Math.log(clusterEnergy)*this.log10inv;
+ for(int i=0; i < this.nLayers; ++i) {
+ layerEnergies[i]/=clusterEnergy;
+ }
+
+ //FIXME sum of cell energies does not equal cluster energy!
+ if (this.debug > 0) {
+ System.out.println(clusterEnergy+" "+clus.getEnergy());
+ for (int i = 0; i < this.nLayers + 1; i++) {
+ System.out.println(i + " " + layerEnergies[i]);
+ }
+ }
+
+ return layerEnergies;
+ }
+
+
// Method mkdirIgnoreExist: Makes a directory "dirname" in the
// ITree "tree" while discarding any exception that arises because
// the directory already exists. Also changes to directory
lcsim/src/org/lcsim/contrib/CarstenHensel
diff -u -r1.1 -r1.2
--- HMatrixAnalyzerPlotter.java 24 Aug 2005 15:10:31 -0000 1.1
+++ HMatrixAnalyzerPlotter.java 25 Aug 2005 19:09:37 -0000 1.2
@@ -24,10 +24,12 @@
private AIDA aida = AIDA.defaultInstance();
- private static final String inputFile1 = "exampleAnalysisJava.aida";
- private static final String inputFile2 = "exampleAnalysisJava2.aida";
- private static final String tupleName = "HMatrixAnalyzerMyTuple";
- private static final String[] columnNames = { "chisq", "chisqDiag", "prob", "clusterEnergy" };
+ private static final String inputFile1 = "cdcaug05_gamma_Theta90_5GeV_cdcaug05_gamma_Theta90_5GeV.hmx.aida";
+ private static final String inputFile2 = "cdcaug05_n_Theta90_20GeV_cdcaug05_gamma_Theta90_5GeV.hmx.aida";
+ private static final String tupleName = "HMatrixAnalyzerHMatrixTuple";
+ String[] columnNames = { "chisq", "chisqDiag", "prob", "logProb", "probDiag", "logProbDiag", "clusterEnergy", "energyPerLayer"};
+ private int nProbVals = 500;
+ private int maxProb = 100;
/** Creates a new instance of HMAtrixAnalyzerPlotter */
public HMatrixAnalyzerPlotter() {
@@ -46,35 +48,118 @@
IHistogramFactory hf = af.createHistogramFactory(tree);
ITree aidaMasterTree = (ITree) ((Studio) Application.getApplication()).getLookup().lookup(ITree.class);
-
-
+ aidaMasterTree.ls();
+
aidaMasterTree.cd("/" + plotter.inputFile1);
ITuple tuple1 = (ITuple)aidaMasterTree.find(plotter.tupleName);
aidaMasterTree.cd("/" + plotter.inputFile2);
ITuple tuple2 = (ITuple)aidaMasterTree.find(plotter.tupleName);
+ IHistogram1D histo1 = hf.createHistogram1D("gamma (5 GeV)", plotter.nProbVals, -1.0 * plotter.maxProb, 0.0);
+ IHistogram1D histo2 = hf.createHistogram1D("neutron (20 GeV)", plotter.nProbVals, -1.0 * plotter.maxProb, 0.0);
+ IHistogram1D histo3 = hf.createHistogram1D("log(chi^2 prob) gamma", plotter.maxProb, -1.0 * plotter.maxProb, 0.0);
+ IHistogram1D histo4 = hf.createHistogram1D("log(chu^2 prob) neutron", plotter.maxProb, -1.0 * plotter.maxProb, 0.0);
+
+ double[] histoEntries1 = new double[plotter.nProbVals];
+ double[] histoEntries2 = new double[plotter.nProbVals];
+
+ double[][] histoEntries = {histoEntries1, histoEntries2};
+
ITuple[] tupleHolder = {tuple1, tuple2};
+ IHistogram1D[] histoHolder = {histo1, histo2};
+ IHistogram1D[] histoHolder2 = {histo3, histo4};
- double[] prob = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9};
- int entries = tuple1.rows();
- tuple1.start();
for (int iTuple = 0; iTuple < tupleHolder.length; iTuple++) {
- double[] sum = plotter.integrate(tupleHolder[iTuple], 2, prob, 1.0);
- for (int i = 0; i < prob.length; i++) {
- plotter.aida.cloud2D((new Integer(iTuple)).toString()).fill(prob[i], sum[i]);
+ double[] sum;
+ sum = plotter.plotEfficiency(histoHolder[iTuple], tupleHolder[iTuple]);
+ plotter.plotLogChiSq(histoHolder2[iTuple], tupleHolder[iTuple]);
+ for (int i = 0; i < plotter.nProbVals; i++) {
+ histoEntries[iTuple][i] = sum[i];
}
}
+
+ for (int i = 0; i < histoEntries1.length; i++) {
+ System.out.println(histoEntries1[i]);
+ plotter.aida.cloud2D("eff vs eff").fill(histoEntries2[i], histoEntries1[i]);
+ }
+
+ IPlotterFactory pf = af.createPlotterFactory();
+ IPlotterStyle style = pf.createPlotterStyle();
+ style.dataStyle().setParameter("showHistogramBars","false");
+ style.dataStyle().setParameter("showDataPoints","false");
+
+
+ IPlotter p1 = af.createPlotterFactory().create("efficiency");
+
+ p1.createRegions(2,2);
+ IPlotterStyle style1 = p1.region(0).style();
+
+
+
+ style1.dataStyle().setParameter("showHistogramBars","false");
+ style1.dataStyle().setParameter("showDataPoints","false");
+ style1.setParameter("backgroundColor","white");
+ style1.setParameter("dataAreaColor","white");
+ style1.setParameter("dataAreaBorderType","");
+ style1.setParameter("showStatisticsBox","false");
+ style1.dataStyle().setParameter("showErrorBars","false");
+ style1.dataStyle().setParameter("showHistogramBars","true");
+ style1.dataStyle().setParameter("fillHistogramBars","true");
+ style1.dataStyle().setParameter("showDataPoints","false");
+ style1.dataStyle().lineStyle().setParameter("color","red");
+ style1.dataStyle().fillStyle().setParameter("color","red");
+
+ p1.show();
+
+ p1.region(0).plot(histo3, style1);
+ p1.region(1).plot(histo1, style1);
+ p1.region(2).plot(histo4, style1);
+ p1.region(3).plot(histo2, style1);
+
+
+
}
- private double[] integrate(ITuple tuple, int rowNumber, double[] probVals, double weight) {
+ private void plotLogChiSq(IHistogram1D histo,ITuple tuple) {
+ int entries = tuple.rows();
+ double val = 0;
+ tuple.start();
+ double weight = 1.0 / entries;
+ for (int i = 0; i < entries; i++) {
+ tuple.next();
+ val = tuple.getDouble(3);
+ if (val < histo.axis().lowerEdge()) {
+ histo.fill(histo.axis().lowerEdge(), weight);
+ } else {
+ histo.fill(tuple.getDouble(3), weight);
+ }
+ }
+ }
+
+
+ private double[] plotEfficiency(IHistogram1D histo,ITuple tuple) {
+
+ int entries = tuple.rows();
+ double[] sum = integrate(tuple, 3, 1.0 / entries);
+ for (int i = 0; i < this.nProbVals; i++) {
+ histo.fill( -1 * i / (this.nProbVals / this.maxProb), sum[i]);
+ }
+ return sum;
+ }
+
+
+
+
+ private double[] integrate(ITuple tuple, int rowNumber, double weight) {
int nTupleRows = tuple.rows();
- double[] integral = new double[probVals.length];
+
+ double[] integral = new double[nProbVals];
tuple.start();
for (int iRow = 0; iRow < nTupleRows; iRow++) {
tuple.next();
- for (int iProb = 0; iProb < probVals.length; iProb++) {
- if (tuple.getDouble(rowNumber) > probVals[iProb]) {
+ for (int iProb = 0; iProb < nProbVals; iProb++) {
+ if (tuple.getDouble(rowNumber) > -1.0 * iProb) {
integral[iProb] += weight;
}
}
lcsim/src/org/lcsim/contrib/CarstenHensel
diff -u -r1.3 -r1.4
--- MainLoop.java 24 Aug 2005 15:11:51 -0000 1.3
+++ MainLoop.java 25 Aug 2005 19:09:37 -0000 1.4
@@ -32,12 +32,12 @@
public static void main(String[] args) throws Exception{
LCSimLoop loop = new LCSimLoop();
- File input = new File("/home/carsten/LC/DataSamples/gamma_Theta90_5GeV_SLIC_v1r9p3_sidaug05.slcio");
+ File input = new File("/home/carsten/LC/DataSamples/gamma_Theta1-179_100MeV-10GeV-0-1000_SLIC_v1r9p3_cdcaug05.slcio");
loop.setLCIORecordSource(input);
loop.add(new EMClusterIDAnalyzer());
// File output = new File("exampleAnalysisJava.slcio");
// loop.add(new LCIODriver(output));
- loop.loop(50);
+ loop.loop(-1);
loop.dispose();
AIDA.defaultInstance().saveAs("/home/carsten/LC/Output/exampleAnalysisJava.aida");
}
lcsim/src/org/lcsim/contrib/CarstenHensel
diff -u -r1.1 -r1.2
--- MyHMatrixBuilder.java 24 Aug 2005 15:11:45 -0000 1.1
+++ MyHMatrixBuilder.java 25 Aug 2005 19:09:37 -0000 1.2
@@ -34,24 +34,65 @@
private boolean initialized = false;
private int nLayers;
private int nmeas;
+ private double cutRadius;
+ private double samplingFraction;
private int logEIndex;
private double[] vals;
private HMatrixBuilder hmb;
private String detectorName;
private CalorimeterIDDecoder decoder;
private static final double log10inv = 1./Math.log(10.0);
+ private int debug = 0;
/** Creates a new instance of MyHMatrixBuilder */
- public MyHMatrixBuilder() {
+ public MyHMatrixBuilder(int debug) {
+ this.debug = debug;
}
+ private void initialize(EventHeader event) {
+
+ this.detectorName = event.getDetectorName();
+ if (this.debug > 0)
+ System.out.println("detector is " + this.detectorName);
+
+ // how to get the sampling fraction from the conditions data base?
+ if (this.detectorName.equals("cdcaug05")) {
+ this.samplingFraction = 0.0113;
+ } else {
+ this.samplingFraction = 0.0117;
+ }
+
+ if (this.debug > 0)
+ System.out.println("Sampling Fraction set to " + this.samplingFraction);
+
+ CylindricalCalorimeter calsub = (CylindricalCalorimeter)event.getDetector().getSubdetectors().get("EMBarrel");
+ this.nLayers = calsub.getLayering().getLayerCount();
+ System.out.println("found " + this.nLayers + " layers in the EMBarrel");
+
+ // the vector of measurements starts as the longitudinal layers
+ this.nmeas = this.nLayers;
+ // add the logarithm of the energy
+ this.logEIndex = this.nmeas;
+ this.nmeas += 1;
+ // would add any additional measurements (e.g. width) here
+ this.vals = new double[this.nmeas];
+
+ //FIXME key needs to be better defined
+ int key = 0;
+
+ this.hmb = new HMatrixBuilder(this.nmeas,key);
+
+
+ this.initialized = true;
+
+
+ }
protected void process(EventHeader event) {
+ if(!this.initialized)
+ this.initialize(event);
- // Get MC Particle information - namely the energy of the final state
- // particle. It looks as if for theta=90 we can get the initial state
- // energy fairly accurately from the final state energy.
List<MCParticle> mcParticles = event.getMCParticles();
@@ -88,29 +129,8 @@
//FIXME: need to get the EM calorimeter names from a conditions file
- if(!this.initialized) {
- CylindricalCalorimeter calsub = (CylindricalCalorimeter)event.getDetector().getSubdetectors().get("EMBarrel");
- this.nLayers = calsub.getLayering().getLayerCount();
- System.out.println("found " + this.nLayers + " layers in the EMBarrel");
- // the vector of measurements starts as the longitudinal layers
- this.nmeas = this.nLayers;
- // add the logarithm of the energy
- this.logEIndex = this.nmeas;
- this.nmeas+=1;
- // would add any additional measurements (e.g. width) here
- this.vals = new double[this.nmeas];
-
- //FIXME key needs to be better defined
- int key = 0;
-
- this.hmb = new HMatrixBuilder(this.nmeas,key);
-
- this.detectorName = event.getDetectorName();
-
- this.initialized = true;
- }
-// FIXME should get calorimeterhit collection names from a conditions file
+ // FIXME should get calorimeterhit collection names from a conditions file
List<CalorimeterHit> collection = event.get(CalorimeterHit.class,"EcalBarrHits");
this.decoder = (CalorimeterIDDecoder) event.getMetaData(collection).getIDDecoder();
@@ -140,16 +160,9 @@
//GWW double energy = cluster.getEnergy();
double energy = cluster.getRawEnergy();
// should be able to fetch this from the cluster...
- double[] layerE = layerEnergies(cluster);
- // accumulate the longitudinal energy fractions...
- for(int i=0; i<layerE.length; ++i) {
- this.vals[i] = layerE[i];
- }
- // Add the correlation to the log of the cluster energy
- //We want the common logarithm (log 10) of the energy
- this.vals[this.logEIndex]=Math.log(energy)*this.log10inv;
- // have now filled the vector of measurements. need to either accumulate the HMatrix or apply it
+ this.vals = getHMInputVals(cluster);
+
this.hmb.accumulate(this.vals);
/*
StringBuffer sb = new StringBuffer();
@@ -162,33 +175,6 @@
- double[] pos = cluster.getPosition();
-
- double x = pos[0];
- double y = pos[1];
- double z = pos[2];
- double r = Math.sqrt(x*x + y*y);
- double l2 = r*r + z*z;
-
- double phi = Math.atan2(y, x);
- //double theta = Math.atan(r/z);
- //double theta = Math.atan2(r, z);
- double theta = Math.acos(z/Math.sqrt(l2));
- double cosTheta = z / Math.sqrt(l2);
-
-
-
- double dphi = phi - phiMC;
- if (dphi > 6.0)
- dphi -= 2.0*Math.PI;
-
- double dtheta = theta - thetaMC;
- //if (dtheta < -3.0)
- // dtheta += Math.PI;
-
-
-
-
}
}
@@ -202,9 +188,9 @@
}
- private double[] layerEnergies(BasicCluster clus) {
+ private double[] getHMInputVals(BasicCluster clus) {
//FIXME could reuse this array
- double[] layerEnergies = new double[this.nLayers];
+ double[] layerEnergies = new double[this.nLayers + 1];
double clusterEnergy = 0.;
List<CalorimeterHit> hits = clus.getCalorimeterHits();
//System.out.println("New cluster with "+hits.size()+ " hits"); // and energy "+clus.getEnergy());
@@ -213,18 +199,37 @@
// GWW double e = hit.getEnergy();
double e = hit.getRawEnergy();
int layer = this.decoder.getLayer();
- //System.out.println("layer "+layer+" energy "+e);
+
+
+ if ((this.detectorName.equals("cdcaug05")) && (layer > 20)) {
+ e = e / (this.samplingFraction * 0.5);
+ } else {
+ e = e / this.samplingFraction;
+ }
+
+ if (this.debug > 0)
+ System.out.println("layer "+layer+" energy "+e);
+
clusterEnergy+=e;
layerEnergies[layer]+=e;
}
- //System.out.println("Cluster energy = " + clusterEnergy + "\n");
+ if (this.debug > 0)
+ System.out.println("Cluster energy = " + clusterEnergy + "\n");
+ layerEnergies[this.nLayers] = Math.log(clusterEnergy)*this.log10inv;
for(int i=0; i < this.nLayers; ++i) {
layerEnergies[i]/=clusterEnergy;
}
+
//FIXME sum of cell energies does not equal cluster energy!
-// System.out.println(clusterEnergy+" "+clus.getEnergy());
+ if (this.debug > 0) {
+ System.out.println(clusterEnergy+" "+clus.getEnergy());
+ for (int i = 0; i < this.nLayers + 1; i++) {
+ System.out.println(i + " " + layerEnergies[i]);
+ }
+ }
+
return layerEnergies;
}
CVSspam 0.2.8
