7 removed files
lcsim/src/org/lcsim/contrib/EricBenavidez/EMClusterID
diff -N EMClusterID_NoE.java
--- EMClusterID_NoE.java 20 Feb 2007 20:40:21 -0000 1.2
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,494 +0,0 @@
-package org.lcsim.contrib.EricBenavidez.EMClusterID;
-
-// ************************************************************************
-// NB: THE FOLLOWING CUTS HAVE BEEN MADE
-// * no zero clusters
-// * 50 keV cell cut
-//
-// Also - Layers are remapped so that the first interacting layer
-// becomes layer 0
-//
-// ************************************************************************
-
-import java.util.*;
-import org.lcsim.geometry.util.CalorimeterIDDecoder;
-import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
-import org.lcsim.recon.cluster.fixedcone.FixedConeClusterer;
-import org.lcsim.recon.cluster.util.BasicCluster;
-import org.lcsim.util.Driver;
-import org.lcsim.util.aida.AIDA;
-import hep.aida.ITree;
-import static java.lang.Math.sqrt;
-import java.text.DecimalFormat;
-import org.lcsim.math.chisq.ChisqProb;
-import org.lcsim.event.*;
-import java.io.*;
-
-/**
- * Reconstruction: EM Clusters
- */
-public class EMClusterID_NoE extends Driver {
- private AIDA aida = AIDA.defaultInstance();
- private ITree _tree;
- private boolean _initialized = false;
- private String _detName;
- private MyCovMatrixTask _task;
- private int _nLayers;
- CalorimeterIDDecoder _decoder;
-
- private HashMap<Integer,MyCovMatrixBuilder> _cmbs;
- private HashMap<Integer,MyCovMatrix> _covs;
-
- private boolean _printChi2 = true;
-
- // the number of variables in the measurement vector
- private int _nmeas;
-
- // the index of the cluster width variables
- private int _widthIndex;
-
- private DecimalFormat _myFormatter = new DecimalFormat("#.###");
- private DecimalFormat twoDigIntFmt = new DecimalFormat("00");
-
- private FileWriter printout = null;
-
- // Path to top-level directory where the covariance matrices
- // are located. Under this directory should be a directory
- // called "CovMatrices" with the individual matrices numbered
- // based on the number of layers from the starting layer to
- // the end. (As generated by "java Standalone_NoE ...")
- //
- private String path = "/home/eric/EMClusterID";
- private String sep = System.getProperty("file.separator");
-
- public EMClusterID_NoE() {
- this(MyCovMatrixTask.ANALYZE);
- }
-
- public EMClusterID_NoE(MyCovMatrixTask task) {
- _tree = aida.tree();
- _task = task;
-
- try {
- printout = new FileWriter(path + sep + "printout");
- }
- catch (IOException e) {
- System.out.println("Failed to open file: " + e.toString());
- System.err.println("Failed to open file: " + e.toString());
- return;
- //System.exit(1);
- }
-
- try {
- _tree.mkdir("LayerClouds");
- _tree.mkdir("Chisq by starting layer");
- }
- catch (IllegalArgumentException e) {} // quietly fail if dir exists
- }
-
- protected void process(EventHeader event) {
- super.process(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();
-
- double thetaMC = 0;
- double cosThetaMC = 0;
- double phiMC = 0;
- double p2 = 0;
- double r_conv = 0;
-
- for (MCParticle mcParticle : mcParticles) {
- if (mcParticle.getGeneratorStatus() == MCParticle.FINAL_STATE) {
- double pX = mcParticle.getPX();
- double pY = mcParticle.getPY();
- double pZ = mcParticle.getPZ();
- double pT = Math.sqrt(pX*pX + pY*pY);
- p2 = pT*pT + pZ*pZ;
-
- phiMC = Math.atan2(pY, pX);
- thetaMC = Math.acos(pZ/Math.sqrt(p2));
- cosThetaMC = pZ / Math.sqrt(p2);
-
- aida.cloud1D("phiMC").fill(phiMC);
- aida.cloud1D("thetaMC").fill(thetaMC);
- aida.cloud1D("cos theta MC").fill(cosThetaMC);
- aida.cloud1D("p").fill(Math.sqrt(p2));
-
- aida.cloud1D("Radius of interaction").fill(mcParticle.getEndPoint().magnitude());
-
- // Rejects photons with radius < 1260. Here, radius is not an
- // angle but is an actual radius.
- //
- // NB: When using a different detector, verify that this
- // cut is a good one!!!
- //
- // Also records the radius of conversion from Monte Carlo
- // information.
- //
- if (mcParticle.getType().getName().equals("gamma")) {
- double r = mcParticle.getEndPoint().magnitude();
- if (r < 1260.0) {
- aida.cloud1D("Rejected Particles").fill(0);
- return;
- }
- else
- r_conv = r;
- }
- }
- }
-
- if(!_initialized) {
- CylindricalCalorimeter calsub = (CylindricalCalorimeter)event.getDetector().getSubdetectors().get("EMBarrel");
- _nLayers = calsub.getLayering().getLayerCount();
- System.out.println("found "+_nLayers+" layers in the EMBarrel");
-
- // the vector of measurements starts as the longitudinal layers
- _nmeas = _nLayers;
-
- //FIXME key needs to be better defined
- int key = 0;
- if(_task == MyCovMatrixTask.ANALYZE) {
- _covs = new HashMap<Integer,MyCovMatrix>();
-
- if (_printChi2) {
- print("layer\tf_meas\tf_exp\tchi2_contrib"+
- "\tsqrt(chi2_contrib)");
- }
- }
- else if(_task == MyCovMatrixTask.BUILD) {
- _cmbs = new HashMap<Integer,MyCovMatrixBuilder>();
- }
- _detName = event.getDetectorName();
- _initialized = true;
- }
-
- List<CalorimeterHit> collection
- = event.get(CalorimeterHit.class,"EcalBarrHits");
- _decoder = (CalorimeterIDDecoder)event.getMetaData(collection).getIDDecoder();
-
- // construct the EM Clusters
- // FIXME get these parameters from conditions file
- double radius = 0.06;
- double seed = 0.;
- double minE = 0.005;
-
- // create the clusterer
- FixedConeClusterer fcc = new FixedConeClusterer(radius, seed, minE);
- List<Cluster> clusters = fcc.createClusters(collection);
- //List<Cluster> clusters = fcc.makeClusters(collection);
-
- // add this list of clusters to the event (for event display)
- event.put("NNClusters r"+radius,clusters);
-
- // Loop over all the clusters
- int nGoodClusters = 0;
- for (Cluster clus : clusters) {
- BasicCluster cluster = (BasicCluster)clus;
- int ncells = cluster.getCalorimeterHits().size();
-
- //FIXME should cut on cluster corrected energy
- if(ncells > 5 && cluster.getRawEnergy() > 0.03) {
- nGoodClusters++;
- double energy = cluster.getRawEnergy();
-
- aida.cloud1D("Number of cells in cluster").fill(ncells);
- aida.cloud1D("energy").fill(energy);
-
- double[][] layerE_and_FE = layerEnergies(cluster);
- double[] layerE = layerE_and_FE[0];
- double[] layerFE = layerE_and_FE[1];
-
- _nmeas = _nLayers;
-
- // Remap longitudinal energy fractions to the start of the
- // shower (i.e. first nonzero layer becomes layer 0, ...,
- // and then quietly set the rest of the layers to 0)
- //
- int j = 0;
- boolean foundFirstNonzero = false;
- int firstNonzero = 0;
-
- for(int i = 0; i < _nLayers; ++i) {
- if (foundFirstNonzero) {
- layerFE[j] = layerFE[i];
- layerE[j] = layerE[i];
- j++;
-
- // if cell cut results in a 0 layer, don't count it
- // for the df in the chi2 calculation
- if (layerFE[i] <= 0.0) {
- _nmeas--;
- }
- }
- else if (layerFE[i] > 0.0) {
- foundFirstNonzero = true;
- firstNonzero = i;
- layerFE[j] = layerFE[i];
- layerE[j] = layerE[i];
- j++;
- }
- else {
- _nmeas--;
- continue;
- }
- }
-
- if (j < _nLayers) {
- for (int l = j; l < _nLayers; l++) {
- layerE[l] = 0;
- layerFE[l] = 0;
- }
- }
-
- double[] vals = new double[j];
- System.arraycopy(layerFE, 0, vals, 0, j);
-
- // have now filled the vector of measurements. need to either
- // accumulate the HMatrix or apply it
- if (_task == MyCovMatrixTask.BUILD) {
- MyCovMatrixBuilder _cmb = _cmbs.get(j);
- if (_cmb == null) {
- _cmbs.put(j, new MyCovMatrixBuilder(j,0));
- _cmb = _cmbs.get(j);
- }
- _cmb.accumulate(vals);
- }
- else if (_task == MyCovMatrixTask.ANALYZE) {
- MyCovMatrix _cov = _covs.get(j);
- if (_cov == null) {
- try {
- _cov = new MyCovMatrix(path + sep + "CovMatrices"
- + sep + j + ".cov");
- }
- /* If the covariance matrix isn't present for this dimension, continue
- * on with the next cluster. */
- catch (IOException e) {
- System.out.println("Covariance matrix for dimension " + j + " not found.");
- continue;
- }
- if (_cov.isSingular())
- continue;
-
- _covs.put(j, _cov);
- }
- else if (_cov.isSingular())
- continue;
-
- aida.cloud1D("nmeas").fill(_nmeas);
-
- double chisq = _cov.chisquared(vals);
- double chisq_prob = ChisqProb.gammq(_nmeas,chisq);
-
- double chisqD = _cov.chisquaredDiagonal(vals);
- double chisqD_prob = ChisqProb.gammq(_nmeas,chisqD);
-
- // Fill the layer clouds
- //
- for (int i = 0; i < _nLayers; i++) {
- aida.cloud2D("Fractional Energy vs Layer").fill(i,layerFE[i]);
- aida.cloud2D("Energy vs Layer").fill(i,layerE[i]);
-
- _tree.cd("LayerClouds");
- aida.cloud1D("Layer " + twoDigIntFmt.format(i) +
- " fractional energy").fill(layerFE[i]);
- aida.cloud1D("Layer " + twoDigIntFmt.format(i) +
- " energy").fill(layerE[i]);
- _tree.cd("..");
- }
-
- if (chisq_prob < 2.0E-3 || firstNonzero != 0 || chisq < 0) {
- _printChi2 = true;
- if (chisq < 0)
- System.out.println("ALERT: NEGATIVE CHI SQUARED!");
- }
-
- if (_printChi2)
- print("Event " + event.getEventNumber() + "\n");
-
- // maximum chi2 contribution
- double chi2c_max = 0;
-
- double[] avg = _cov.getAvgVector();
- double[] var = _cov.getVarVector();
-
- for (int i = 0; i < j; i++) {
-
- double sqrt_chi2c = (layerFE[i]-avg[i]) /
- Math.sqrt(var[i]);
- double chi2c = sqrt_chi2c*sqrt_chi2c;
-
- if (chi2c > chi2c_max)
- chi2c_max = chi2c;
-
- if (_printChi2)
- print(i + "\t" + layerFE[i] + "\t" + avg[i] + "\t"
- + chi2c+"\t" + sqrt_chi2c);
- }
-
- if (_printChi2) {
- print("\nTotal Chisq " + chisq);
- print("Chisq prob " + chisq_prob);
- print("Total diag Chisq " + chisqD);
- print("ChisqD prob " + chisqD_prob);
- print("First nonzero " + firstNonzero + "\n");
- }
-
- if (event.getEventNumber() >= 19)
- _printChi2 = false;
-
- double max_to_total = chi2c_max / chisq;
- double max_to_totalD = chi2c_max / chisqD;
- aida.cloud1D("Chisq max / total").fill(max_to_total);
- aida.cloud2D("Chisq max / total vs. prob").fill(chisq_prob, max_to_total);
- aida.cloud1D("ChisqD max / total").fill(max_to_totalD);
- aida.cloud2D("ChisqD max / total vs. prob").fill(chisqD_prob, max_to_totalD);
-
- aida.cloud1D("Chisq").fill(chisq);
- aida.cloud2D("Chisq vs energy").fill(energy,chisq);
- aida.cloud1D("Chisq Probability").fill(chisq_prob);
- aida.cloud1D("log(Chisq Prob)").fill(Math.log10(chisq_prob));
- aida.cloud2D("Radius vs. Chisq Prob.").fill(r_conv,chisq_prob);
-
- aida.cloud1D("ChisqD").fill(chisqD);
- aida.cloud2D("ChisqD vs energy").fill(energy,chisqD);
- aida.cloud1D("ChisqD Probability").fill(chisqD_prob);
- aida.cloud1D("log(ChisqD Prob)").fill(Math.log10(chisqD_prob));
- aida.cloud2D("Radius vs. ChisqD Prob.").fill(r_conv,chisqD_prob);
-
- double chisqND = chisq - chisqD;
- aida.cloud1D("ChisqND").fill(chisqND);
-
- _tree.cd("Chisq by starting layer");
-
- aida.cloud2D("Chisq vs. start layer").fill(firstNonzero, chisq);
- aida.cloud2D("ChisqD vs. start layer").fill(firstNonzero, chisqD);
-
- aida.cloud2D("Chisq prob vs. start layer").fill(firstNonzero, chisq_prob);
- aida.cloud2D("ChisqD prob vs. start layer").fill(firstNonzero, chisqD_prob);
- aida.cloud2D("log(Chisq prob) vs. start layer").fill(Math.log10(firstNonzero), chisq_prob);
- aida.cloud2D("log(ChisqD prob) vs. start layer").fill(Math.log10(firstNonzero), chisqD_prob);
-
- aida.cloud1D("Chisq for start layer " + twoDigIntFmt.format(firstNonzero)).fill(chisq);
- aida.cloud1D("ChisqD for start layer " + twoDigIntFmt.format(firstNonzero)).fill(chisqD);
-
- aida.cloud1D("Chisq prob for start layer " + twoDigIntFmt.format(firstNonzero)).fill(chisq_prob);
- aida.cloud1D("ChisqD prob for start layer " + twoDigIntFmt.format(firstNonzero)).fill(chisqD_prob);
- aida.cloud1D("log(Chisq prob) for start layer " + twoDigIntFmt.format(firstNonzero)).fill(Math.log10(chisq_prob));
- aida.cloud1D("log(ChisqD prob) for start layer " + twoDigIntFmt.format(firstNonzero)).fill(Math.log10(chisqD_prob));
-
- if (firstNonzero <= 7) {
- aida.cloud1D("log(Chisq prob) for start layers 0-7").fill(Math.log10(chisq_prob));
- aida.cloud1D("log(ChisqD prob) for start layers 0-7 ").fill(Math.log10(chisqD_prob));
- }
-
- _tree.cd("..");
- }
-
- double[] pos = cluster.getPosition();
- aida.cloud2D("centroid x vs y").fill(pos[0], pos[1]);
- aida.cloud1D("centroid radius").fill( sqrt(pos[0]*pos[0] + pos[1]*pos[1]) );
-
- 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.acos(z/Math.sqrt(l2));
- double cosTheta = z / Math.sqrt(l2);
-
- aida.cloud1D("phi").fill(phi);
- aida.cloud1D("theta").fill(theta);
- aida.cloud1D("cos theta").fill(cosTheta);
- aida.cloud1D("l").fill(Math.sqrt(l2));
-
- double dphi = phi - phiMC;
- if (dphi > 6.0)
- dphi -= 2.0*Math.PI;
-
- double dtheta = theta - thetaMC;
-
- aida.cloud1D("d phi").fill(dphi);
- aida.cloud1D("d theta").fill(dtheta);
- aida.cloud1D("d cos theta").fill(cosTheta - cosThetaMC);
- }
- }
- aida.cloud1D("number of found clusters (above hits threshold)").fill(nGoodClusters);
- }
-
- // Method layerEnergies: Returns a two-dimensional array indexed as
- // follows:
- // First index: 0 = layer energies, 1 = layer fractional energies
- // Second index: by layer number, starting with 0
- //
- private double[][] layerEnergies(Cluster clus) {
- double[][] retval = new double[2][_nLayers];
- double clusterEnergy = 0.;
- List<CalorimeterHit> hits = clus.getCalorimeterHits();
-
- for(CalorimeterHit hit : hits) {
- _decoder.setID(hit.getCellID());
- double e = hit.getRawEnergy();
- int layer =_decoder.getLayer();
-
- // Make a 50 keV = 5e-5 GeV energy cut
- if (e > 5.0E-5) {
- clusterEnergy+=e;
- retval[0][layer]+=e;
- retval[1][layer]+=e;
- }
- }
- for (int i = 0; i < _nLayers; i++) {
- retval[1][i] /= clusterEnergy;
- }
- return retval;
- }
-
- protected void endOfData() {
- if (_task == MyCovMatrixTask.BUILD) {
- Collection<MyCovMatrixBuilder> cmbs_coll = _cmbs.values();
-
- for (MyCovMatrixBuilder _cmb : cmbs_coll) {
- System.out.println("Validating covariance matrix of dimension "
- + _cmb.getDim());
- _cmb.validate();
- String fileLocation = path + sep + "CovMatrices" + sep +
- _cmb.getDim() + ".cov";
- System.out.println("Writing out covariance matrix to "
- + fileLocation);
- _cmb.write(fileLocation,commentForCovMatrix());
- }
- }
- try {
- printout.close();
- }
- catch (IOException e) {
- System.err.println("Failed to close file: " + e.toString());
- }
- }
-
- private String commentForCovMatrix() {
- Calendar cal = new GregorianCalendar();
- Date date = new Date();
- cal.setTime(date);
- DecimalFormat formatter = new DecimalFormat("00");
- String day = formatter.format(cal.get(Calendar.DAY_OF_MONTH));
- String month = formatter.format(cal.get(Calendar.MONTH)+1);
- String myDate = cal.get(Calendar.YEAR)+month+day;
- return _detName+" "+myDate+" "+System.getProperty("user.name");
- }
-
- private void print(String str) {
- try {
- printout.write(str+"\n");
- }
- catch (IOException e) {
- System.err.println("Failed to write to file: " + e.toString());
- }
- }
-}
lcsim/src/org/lcsim/contrib/EricBenavidez/EMClusterID
diff -N MyCovMatrix.java
--- MyCovMatrix.java 20 Feb 2007 20:40:21 -0000 1.2
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,368 +0,0 @@
-package org.lcsim.contrib.EricBenavidez.EMClusterID;
-
-import java.io.*;
-import Jama.*;
-import java.util.*;
-
-public class MyCovMatrix implements Serializable {
- private double[][] _covMat; // the covariance matrix
- private double[] _vec; // the vector of measurements
- private double[] _covMeas; // the variance on the measurements
- private int _dim; // the dimensionality
- private int _key; // the key for the matrix
- private double[] _tmp; // temporary array
- private double[] _tmp2; // temporary array
- private double[][] _invcovMat; // inverse of covariance matrix
- private boolean _isSingular;
-
- /**
- * Constructor to create a covariance matrix from raw data.
- *
- * @param dim The dimensionality of this covariance matrix.
- * @param key The key for indexing this matrix.
- * @param vec The array of average values for the measurement space.
- * @param cov The matrix for the averages.
- * @param covMat The matrix.
- */
- public MyCovMatrix(int dim, int key, double [] vec, double[] cov,
- double[][] covMat) {
- _dim = dim;
- _key = key;
- _vec = new double[_dim];
- _covMeas = new double[_dim];
- System.arraycopy(vec, 0, _vec, 0, _dim);
- System.arraycopy(cov, 0, _covMeas, 0, _dim);
- _covMat = new double[_dim][_dim];
- _isSingular = false;
-
- for(int i =0; i<_dim; ++i) {
- System.arraycopy(covMat[i], 0, _covMat[i], 0, _dim);
- }
-
- _tmp = new double[_dim];
- _tmp2 = new double[_dim];
-
- try {
- _invcovMat = invert();
- }
- catch (RuntimeException e) {
- System.out.println("Covariance matrix for dimension " + _dim + " is singular.");
- _isSingular = true;
- }
- }
-
- /**
- * Constructor to read the covariance matrix from a file.
- *
- * If the file is not present, or there is an error reading the file,
- * go ahead and throw the underlying IOException and let the caller
- * handle it.
- *
- * @param resourceFileName The name of the file where the covariance matrix is located.
- *
- */
- public MyCovMatrix(String resourceFileName) throws IOException {
- InputStream in = null;
- _isSingular = false;
-
- in = new FileInputStream(resourceFileName);
-
- BufferedReader bufferedreader
- = new BufferedReader(new InputStreamReader(in));
- readTheMatrix(bufferedreader);
- bufferedreader.close();
-
- _tmp = new double[_dim];
- _tmp2 = new double[_dim];
-
- try {
- _invcovMat = invert();
- }
- catch (RuntimeException e) {
- System.out.println("Covariance matrix for dimension " + _dim + " is singular.");
- _isSingular = true;
- }
- }
-
- public boolean isSingular() {
- return _isSingular;
- }
-
- public double[] getAvgVector() {
- double[] avgvec = new double[_dim];
- System.arraycopy(_vec, 0, avgvec, 0, _dim);
- return avgvec;
- }
-
- public double[] getVarVector() {
- double[] covvec = new double[_dim];
- System.arraycopy(_covMeas, 0, covvec, 0, _dim);
- return covvec;
- }
-
- public double[][] getCovMatrix() {
- double[][] covmat = new double[_dim][_dim];
- for (int i = 0; i < _dim; i++) {
- System.arraycopy(_covMat[i], 0, covmat[i], 0, _dim);
- }
- return covmat;
- }
-
- public double getDim() {
- return _dim;
- }
-
- /**
- * Calculates the chi-squared for the measurement compared to
- * the expected values and covariances represented by the matrix.
- *
- * @param dat The array of measured values
- * @return The value of chi-squared.
- */
- public double chisquared(double[] dat) {
- int dim = dat.length;
-
- // vector of measured-predicted values
- for(int i = 0; i < dim; ++i) {
- _tmp[i] = dat[i] - _vec[i];
- }
-
- double chisqr = 0.;
-
- // covariance matrix times difference vector
- for(int i = 0; i < dim; ++i) {
- _tmp2[i] = 0.;
- for(int j = 0; j < dim; ++j) {
- _tmp2[i] += _invcovMat[j][i] * _tmp[j];
- }
- chisqr += _tmp[i]*_tmp2[i];
- }
- return chisqr;
- }
-
- /**
- * Calculates the diagonal chi-squared for the measurement compared to
- * the expected values and covariances represented by the HMatrix,
- * neglecting correlations.
- *
- * @param dat The array of measured values
- * @return The value of the diagonal chi-squared.
- */
- public double chisquaredDiagonal(double[] dat) {
- double chisqr = 0.;
- // diagonal terms of covariance matrix times difference vector
- for(int i = 0; i < dat.length; ++i) {
- //measured-predicted values
- _tmp[i] = dat[i] - _vec[i];
- //squared, divided by the variance of the prediction
- chisqr += _tmp[i] * _tmp[i] / _covMeas[i];
- }
- return chisqr;
- }
-
- /**
- *Output stream
- *
- * @return A String representation of the object.
- */
- public String toString() {
- StringBuffer sb = new StringBuffer("CovMatrix: dimension "+_dim+" key "+_key+ "\n");
- sb.append("vector: \n");
- for(int i = 0; i<_dim; ++i) {
- sb.append(_vec[i]+" ");
- }
- sb.append("\n\nCovariance matrix: \n");
- for(int i = 0; i<_dim; ++i) {
- for(int j = 0; j<_dim; ++j) {
- sb.append(_covMat[i][j]+" ");
- }
- sb.append("\n");
- }
- return sb.toString();
- }
-
- /**
- * Writes out the covariance matrix to an ASCII file
- */
- public void write(String filename, String comment) {
- System.out.println("Writing covariance matrix to '" + filename + "'.");
- try {
- PrintWriter printwriter = new PrintWriter(new BufferedWriter(new FileWriter(filename)));
- // Print the comment
- String s2 = "# " + comment;
- printwriter.println(s2);
-
- // Now the dimension and index of the matrix
- String s3 = _dim+" "+_key;
- printwriter.println(s3);
- String s4 = "";
-
- //Now the vector of average values
- String s5 = "";
- for(int i = 0; i < _dim; i++)
- s5 = s5 + _vec[i] + " ";
- printwriter.println(s5.substring(0, s5.length() - 1));
-
- //Now the vector of variance of average values
- String s7 = "";
- for(int i = 0; i < _dim; i++)
- s7 = s7 + _covMeas[i] + " ";
- printwriter.println(s7.substring(0, s7.length() - 1));
-
- //Now the covariance matrix
- for(int i = 0; i < _dim; i++) {
- String s6 = "";
- for(int j = 0; j < _dim; j++)
- s6 = s6 + _covMat[i][j] + " ";
-
- printwriter.println(s6.substring(0, s6.length() - 1));
- }
-
- printwriter.close();
- }
- catch(IOException _ex) {
- System.err.println("MyCovMatrix::write -> Error writing to '" +
- filename + "'.");
- return;
- //System.exit(0);
- }
- System.out.println("Covariance matrix written.");
- }
-
- private void readTheMatrix(BufferedReader bufferedreader)
- throws IOException {
- boolean flag = false;
- boolean readVec = false;
- boolean readCovDiag = false;
- int i = 0;
- int ii = 0;
- int j = 0;
-
- for(String s1 = new String();
- (s1 = bufferedreader.readLine()) != null;) {
- if(s1.startsWith("#"))
- System.out.println(s1);
- else {
- if(!flag) {
- _dim = Integer.parseInt(s1.substring(0, s1.indexOf(" ")));
- _key = Integer.parseInt(s1.substring(s1.indexOf(" ") + 1, s1.length()));
-
- System.out.println("_dim = " + _dim + " _key = " + _key);
-
- _vec = new double[_dim];
- _covMeas = new double[_dim];
- _covMat = new double[_dim][_dim];
- flag = true;
- }
- else {
- if (!readVec) {
- // Read the vector of averages
- while(s1.length() > 0 && s1.indexOf(" ") != -1) {
- double d = Double.valueOf(s1.substring(0, s1.indexOf(" "))).doubleValue();
- s1 = s1.substring(s1.indexOf(" ") + 1, s1.length());
- _vec[ii] = d;
- ii++;
- }
- double d1 = Double.valueOf(s1).doubleValue();
- _vec[ii] = d1;
- readVec = true;
- ii=0;
- }
- else if(!readCovDiag) {
- // Read the vector of variance of averages
- while(s1.length() > 0 && s1.indexOf(" ") != -1) {
- double d = Double.valueOf(
- s1.substring(0,s1.indexOf(" "))).doubleValue();
- s1 = s1.substring(s1.indexOf(" ")+1, s1.length());
- _covMeas[ii] = d;
- ii++;
- }
- double d1 = Double.valueOf(s1).doubleValue();
- _covMeas[ii] = d1;
- readCovDiag = true;
- }
- else {
- // Read the covariance matrix
- while(s1.length() > 0 && s1.indexOf(" ") != -1) {
- double d = Double.valueOf(
- s1.substring(0,s1.indexOf(" "))).doubleValue();
- s1 = s1.substring(s1.indexOf(" ")+1, s1.length());
- _covMat[i][j] = d;
- j++;
- }
- double d1 = Double.valueOf(s1).doubleValue();
- _covMat[i][j] = d1;
- i++;
- j = 0;
- }
- }
- }
- }
- }
-
- // Invert the covariance matrix _covMat with the Jama standard
- // inverse() calculation.
- //
- private double[][] invert() {
- Matrix M = new Matrix(_covMat, _dim, _dim);
- Matrix Minv = M.inverse();
-
- Matrix leftEM = Minv.times(M).minus(Matrix.identity(_dim,_dim));
- Matrix rightEM = M.times(Minv).minus(Matrix.identity(_dim,_dim));
-
- double leftE = leftEM.normF();
- double rightE = rightEM.normF();
-
- System.out.println("||M^-1 * M - I||_F = " + leftE);
- System.out.println("||M * M^-1 - I||_F = " + rightE);
-
- return Minv.getArrayCopy();
- }
-
-
- // Invert the covariance matrix _covMat using its SVD M = U*S*V'.
- // M^(-1) = V*S^(-1)*U'
- //
- private double[][] invert_old() {
- Matrix M = new Matrix(_covMat, _dim, _dim);
- SingularValueDecomposition svdM = M.svd();
- Matrix Ut = svdM.getU().transpose();
- double[][] V = svdM.getV().getArrayCopy();
- int i, j;
-
- // After the loop, Sinv contains the inverse of Sigma
- double[][] Sinv = svdM.getS().getArrayCopy();
-
- for (i = 0; i < _dim; i++) {
- if (Sinv[i][i] != 0) {
- Sinv[i][i] = 1.0 / Sinv[i][i];
- }
- }
-
- // Multiply product V * Sigma^(-1) in place using that Sigma^(-1)
- // is diagonal; store in V
- //
- for (j = 0; j < _dim; j++) {
- for (i = 0; i < _dim; i++) {
- V[i][j] *= Sinv[j][j];
- }
- }
-
- Matrix VSinv = new Matrix(V, _dim, _dim);
- Matrix Minv = VSinv.times(Ut);
-
- System.out.println("Covariance matrix dimension " + _dim +
- " inversion error:");
-
- Matrix leftEM = Minv.times(M).minus(Matrix.identity(_dim,_dim));
- Matrix rightEM = M.times(Minv).minus(Matrix.identity(_dim,_dim));
-
- double leftE = leftEM.normF();
- double rightE = rightEM.normF();
-
- System.out.println("||M^-1 * M - I||_F = " + leftE);
- System.out.println("||M * M^-1 - I||_F = " + rightE);
-
- return Minv.getArrayCopy();
- }
-}
lcsim/src/org/lcsim/contrib/EricBenavidez/EMClusterID
diff -N MyCovMatrixBuilder.java
--- MyCovMatrixBuilder.java 20 Feb 2007 20:40:21 -0000 1.2
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,100 +0,0 @@
-package org.lcsim.contrib.EricBenavidez.EMClusterID;
-
-import java.io.*;
-import Jama.*;
-
-public class MyCovMatrixBuilder {
- private double[][] _cov;
- private double[] _vec;
- private double[] _vals;
- private double[][] _valsq;
- private int _dim;
- private int _key;
- private int _ndata;
- private boolean _isValid;
- private MyCovMatrix _covmat;
-
- /** Construct an <b>n</b> x <b>n</b> MyCovMatrixBuilder.
- * Leaves covariance matrix in an invalid state.
- * One must either accumulate entries to build up the
- * covariance matrix or read in the covariance matrix
- * values from a file.
- *
- * @param n The dimension of the measurement space.
- * @param key The key by which to index this matrix.
- */
- public MyCovMatrixBuilder(int n, int key) {
- _dim = n;
- _key = key;
- _cov = new double[_dim][_dim];
- _vec = new double[_dim];
- _vals = new double[_dim];
- _valsq = new double[_dim][_dim];
- _isValid = false;
- }
-
- public int getDim() {
- return _dim;
- }
-
- /**
- * Add the measurement vector to the accumulated matrix
- *
- * @param dat The array of measurements.
- */
- public void accumulate(double[] dat) {
- for(int i = 0; i < _dim; ++i) {
- for(int j = 0; j < _dim; ++j) {
- _valsq[i][j] += dat[i]*dat[j];
-
- if(i == j)
- _vals[i] += dat[i];
- }
- }
- ++_ndata;
- }
-
- /**
- * Validates the matrix by performing the averages
- */
- public void validate() {
- System.out.println("Number of data: " + _ndata);
-
- double[] _covDiag = new double[_dim];
-
- for(int i = 0; i < _dim; ++i) {
- _vec[i] = _vals[i]/_ndata;
-
- for(int j = 0; j < _dim; ++j) {
- _cov[i][j] = (_valsq[i][j]/_ndata) - ((_vals[i]/_ndata)*(_vals[j]/_ndata));
- }
- _covDiag[i] = _cov[i][i];
- }
-
- _covmat = new MyCovMatrix(_dim, _key, _vec, _covDiag, _cov);
- _isValid = true;
- }
-
- /**
- * Output Stream
- *
- * @return String representation of the object
- */
- public String toString() {
- StringBuffer sb = new StringBuffer("HMatrixBuilder: dimension "+_dim+" ");
- if(!_isValid)
- return (sb.append("INVALID").toString());
- sb.append(_cov);
- return sb.toString();
- }
-
- /**
- * Writes out the HMatrix to an ASCII file
- *
- * @param filename The file to which to write.
- * @param comment A comment for the header.
- */
- public void write(String filename, String comment) {
- _covmat.write(filename, comment);
- }
-}
lcsim/src/org/lcsim/contrib/EricBenavidez/EMClusterID
diff -N MyCovMatrixTask.java
--- MyCovMatrixTask.java 20 Feb 2007 20:40:21 -0000 1.2
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,9 +0,0 @@
-package org.lcsim.contrib.EricBenavidez.EMClusterID;
-
-/*
- * MyCovMatrixTask.java
- *
- */
-
-public enum MyCovMatrixTask
-{BUILD, ANALYZE};
lcsim/src/org/lcsim/contrib/EricBenavidez/EMClusterID
diff -N PerformanceAnalysis.java
--- PerformanceAnalysis.java 20 Feb 2007 20:40:21 -0000 1.2
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,101 +0,0 @@
-package org.lcsim.contrib.EricBenavidez.EMClusterID;
-
-import hep.aida.*;
-import org.freehep.application.*;
-import org.freehep.application.studio.*;
-
-/* Performance analysis: Plots the purity of neutrons vs. the efficiency of
- * photons for two sets of data opened at the same time.
- *
- * NB: Make sure that, before running this, that the paths
- * /gamma.aida and /neutron.aida
- * exist in the JAS3 main tree.
- *
- * This program is meant to be run under JAS3 and does not save an AIDA
- * file when its analysis is complete. A version that does so may be created
- * as time permits.
- */
-public class PerformanceAnalysis {
- public static void main(String[] args) {
- IAnalysisFactory af = IAnalysisFactory.create();
- IPlotterFactory pf = af.createPlotterFactory();
-
- ITree tree = (ITree)((Studio)Application.getApplication()).getLookup().lookup(ITree.class);
- IHistogramFactory hf = af.createHistogramFactory(tree);
-
- ICloud1D origGammaCloud = (ICloud1D)tree.find("/gamma.aida/log(Chisq Prob)");
- ICloud1D origGammaCloudD = (ICloud1D)tree.find("/gamma.aida/log(ChisqD Prob)");
- ICloud1D origNCloud = (ICloud1D)tree.find("/neutron.aida/log(ChisqD Prob)");
- ICloud1D origNCloudD = (ICloud1D)tree.find("/neutron.aida/log(Chisq Prob)");
-
- int gammaTotEnt = origGammaCloud.entries();
- int gammaTotEntD = origGammaCloudD.entries();
- int nTotEnt = origNCloud.entries();
- int nTotEntD = origNCloudD.entries();
-
- ICloud1D gammaCloud = hf.createCopy("gamma_cloud", origGammaCloud);
- ICloud1D gammaCloudD = hf.createCopy("gamma_cloudD", origGammaCloudD);
- ICloud1D nCloud = hf.createCopy("n_cloud", origNCloud);
- ICloud1D nCloudD = hf.createCopy("n_cloudD", origNCloudD);
-
- int gammaLEdge = (int)Math.floor(gammaCloud.lowerEdge());
- int gammaLEdgeD = (int)Math.floor(gammaCloudD.lowerEdge());
- int nLEdge = gammaLEdge;
- int nLEdgeD = gammaLEdgeD;
-
- gammaCloud.convert((int)Math.abs(gammaLEdge)+1, gammaLEdge-0.5, 0.5);
- gammaCloudD.convert((int)Math.abs(gammaLEdgeD)+1, gammaLEdgeD-0.5, 0.5);
- nCloud.convert((int)Math.abs(nLEdge)+1, nLEdge-0.5, 0.5);
- nCloudD.convert((int)Math.abs(nLEdgeD)+1, nLEdgeD-0.5, 0.5);
-
- IHistogram1D gammaHist = gammaCloud.histogram();
- IHistogram1D gammaHistD = gammaCloudD.histogram();
- IHistogram1D nHist = nCloud.histogram();
- IHistogram1D nHistD = nCloudD.histogram();
-
- int gammaEnt = 0;
- int gammaEntD = 0;
- int nEnt = 0;
- int nEntD = 0;
-
- ICloud2D perfPlot = hf.createCloud2D("Performance plot from chisq");
- ICloud2D perfPlotD = hf.createCloud2D("Performance plot from chisqD");
-
- for (int i = 0; i >= gammaLEdge; i--) {
- gammaEnt += gammaHist.binEntries(gammaHist.coordToIndex(i));
- gammaEntD += gammaHistD.binEntries(gammaHistD.coordToIndex(i));
- nEnt += nHist.binEntries(nHist.coordToIndex(i));
- nEntD += nHistD.binEntries(nHistD.coordToIndex(i));
-
- double gammaFrac = (double)gammaEnt / (double)gammaTotEnt;
- double gammaFracD = (double)gammaEntD / (double)gammaTotEntD;
- double nFrac = 1.0 - (double)nEnt / (double)nTotEnt;
- double nFracD = 1.0 - (double)nEntD / (double)nTotEntD;
-
- perfPlot.fill(gammaFrac, nFrac);
- perfPlotD.fill(gammaFracD, nFracD);
- }
-
- IPlotter plotter = pf.create("Performance plot from chisq");
- IPlotterStyle pStyle = plotter.currentRegion().style();
- pStyle.dataStyle().markerStyle().setColor("black");
- plotter.currentRegion().plot(perfPlot);
- plotter.show();
-
- IPlotter plotterD = pf.create("Performance plot from chisqD");
- IPlotterStyle pStyleD = plotterD.currentRegion().style();
- pStyleD.dataStyle().markerStyle().setColor("black");
- plotterD.currentRegion().plot(perfPlotD);
- plotterD.show();
-
- // Comment out the following two lines if closing off gamma.aida
- // and neutron.aida gets to be annoying.
- tree.unmount("gamma.aida");
- tree.unmount("neutron.aida");
-
- tree.rm("gamma_cloud");
- tree.rm("gamma_cloudD");
- tree.rm("n_cloud");
- tree.rm("n_cloudD");
- }
-}
lcsim/src/org/lcsim/contrib/EricBenavidez/EMClusterID
diff -N Standalone_NoE.java
--- Standalone_NoE.java 11 Sep 2007 00:21:02 -0000 1.3
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,88 +0,0 @@
-package org.lcsim.contrib.EricBenavidez.EMClusterID;
-
-import java.io.File;
-import java.io.IOException;
-import org.lcsim.util.Driver;
-import org.lcsim.util.aida.AIDA;
-import org.lcsim.util.loop.LCSimLoop;
-import org.lcsim.util.loop.LCIOEventSource;
-
-/**
- * A simple standalone EMClusterID example
- */
-public class Standalone_NoE extends Driver {
- public Standalone_NoE() {}
-
- public static void main(String[] args) throws Exception {
- MyCovMatrixTask task = MyCovMatrixTask.BUILD;
- File input;
- int numToProcess = -1;
-
- if(args.length < 1 || args.length > 3) {
- usage();
- System.exit(1);
- }
-
- // Task selection: build or analyze
- // Default behavior: build if switch not present
- // - but return if switch is invalid
- //
- if (args[0].charAt(0) == '-') {
- if (args.length < 2 || args.length > 3) {
- usage();
- System.exit(1);
- }
-
- // analyze
- if (args[0].equals("-a"))
- task = MyCovMatrixTask.ANALYZE;
- else if (args[0].equals("-b"))
- task = MyCovMatrixTask.BUILD;
- else {
- usage();
- System.exit(1);
- }
-
- input = new File(args[1]);
-
- if (args.length == 3)
- numToProcess = Integer.parseInt(args[2]);
- }
- else {
- if (args.length > 2) {
- usage();
- System.exit(1);
- }
- input = new File(args[0]);
-
- if (args.length == 2)
- numToProcess = Integer.parseInt(args[1]);
- }
-
- System.out.println("Processing " + numToProcess+" events from "+input);
- LCSimLoop loop = new LCSimLoop();
-
- LCIOEventSource evtSrc = new LCIOEventSource(input);
- loop.setRecordSource(evtSrc);
-
- System.out.println("adding the driver");
-
- loop.add(new EMClusterID_NoE(task));
-
- System.out.println("looping");
- loop.loop(numToProcess, null);
-
- loop.dispose();
- AIDA.defaultInstance().saveAs("emClusterAnalysis.aida");
- }
-
- public static void usage() {
- System.out.println("This is Standalone_NoE");
- System.out.println("Usage:");
- System.out.println(" java Standalone_NoE [-a|-b] inputFile [# events]");
- System.out.println("where");
- System.out.println(" -a specifies that task is to analyze");
- System.out.println(" -b specifies that task is to build (default)");
- System.out.println(" # events specifies the number of events to process (default -1: process all)");
- }
-}
lcsim/src/org/lcsim/contrib/EricBenavidez/EMClusterID
diff -N classpath.sh
--- classpath.sh 3 Nov 2006 21:05:34 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,33 +0,0 @@
-#!/bin/bash
-
-# LCSIM_SRC: the root of the org.lcsim source tree. Is where GeomConverter
-# and lcsim are stored.
-
-#LCSIM_SRC=/home/eric/lcsim_1.1
-LCSIM_SRC=/home/eric/lcsim_061103
-
-#for i in $HOME/.JAS3/extensions/*.jar
-#do
-# export CLASSPATH=$i:$CLASSPATH
-#done
-
-for i in $HOME/.maven/repository/*
-do
- for j in $i/jars/*.jar
- do
- export CLASSPATH=$j:$CLASSPATH
- done
-done
-
-for i in $LCSIM_SRC/GeomConverter/target/*.jar
-do
- export CLASSPATH=$i:$CLASSPATH
-done
-
-for i in $LCSIM_SRC/lcsim/target/*.jar
-do
- export CLASSPATH=$i:$CLASSPATH
-done
-
-export CLASSPATH=$PWD/Jama-1.0.2.jar:$CLASSPATH
-
CVSspam 0.2.8
