6 removed files
lcsim/src/org/lcsim/contrib/CarstenHensel
diff -N EMClusterIDAnalyzer.java
--- EMClusterIDAnalyzer.java 25 Aug 2005 19:09:36 -0000 1.4
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,100 +0,0 @@
-/*
- * EMClusterIDAnalyzer.java
- *
- * Created on August 21, 2005, 4:16 PM
- *
- * To change this template, choose Tools | Options and locate the template under
- * the Source Creation and Management node. Right-click the template and choose
- * Open. You can then make changes to the template in the Source Editor.
- */
-
-package org.lcsim.contrib.CarstenHensel;
-import org.lcsim.event.*;
-import org.lcsim.util.Driver;
-import org.lcsim.util.aida.AIDA;
-import hep.aida.*;
-import hep.aida.ICloud2D;
-import java.text.*;
-
-
-/**
- *
- *
- *
- *
- */
-public class EMClusterIDAnalyzer extends Driver {
- // This is convenient for creating and filling histograms "on-the-fly".
- private AIDA aida = AIDA.defaultInstance();
- private IAnalysisFactory af;
- private IHistogramFactory hf;
- private ITree tree;
- private DecimalFormat df;
- private static final double SAMP_FRAC = 0.012;
-
- // debugging flag:
- // 0 no print out
- private static final int debug = 0;
- private static final int histoLevel = 0;
- private static final double interactionRadius = 1260.0;
- private static int eventCounter = 0;
- private static int passedEvents = 0;
-
- //carsten's histograms:
- private ICloud2D angleVSk;
- private IHistogram2D angleVSkHist;
-
- // 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
- FixedConeClustererAnalyzer fixedConeClustererAnalyzer = new FixedConeClustererAnalyzer(this.tree);
- fixedConeClustererAnalyzer.setAttributes(this.SAMP_FRAC, this.debug, this.interactionRadius);
- //this.add(fixedConeClustererAnalyzer);
-
- HMatrixAnalyzer hMatrixAnalyzer = new HMatrixAnalyzer("/home/carsten/LC/Output/cdcaug05_gamma_Theta90_5GeV.hmx", tree);
- //HMatrixAnalyzer hMatrixAnalyzer = new HMatrixAnalyzer(tree);
- this.add(hMatrixAnalyzer);
-
- int debug = 0;
- MyHMatrixBuilder builder = new MyHMatrixBuilder(this.debug);
- //this.add(builder);
- }
-
-
-
-
-
-
-
-
- protected void startOfData(){
- super.startOfData();
- }
- public void process(EventHeader event) {
- super.process(event);
- this.eventCounter++;
- if ((this.eventCounter % 50) == 0) {
- System.out.println("processing event " + this.eventCounter);
- }
- }
-
-
- protected void endOfData() {
- super.endOfData();
- System.out.println("Events analyzed " + this.eventCounter);
- System.out.println("Events passed pre-selection " + this.passedEvents);
-
- }
-}
-
-
lcsim/src/org/lcsim/contrib/CarstenHensel
diff -N FixedConeClustererAnalyzer.java
--- FixedConeClustererAnalyzer.java 14 Feb 2007 00:16:04 -0000 1.4
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,165 +0,0 @@
-/*
- * FixedConeClustererAnalyzer.java
- *
- * Created on August 21, 2005, 8:08 PM
- *
- *
- */
-
-package org.lcsim.contrib.CarstenHensel;
-import hep.aida.ITree;
-import hep.aida.ref.tree.TreeObjectAlreadyExistException;
-
-import java.text.DecimalFormat;
-import java.text.FieldPosition;
-import java.util.HashMap;
-import java.util.List;
-import java.util.Map;
-
-import org.lcsim.event.CalorimeterHit;
-import org.lcsim.event.Cluster;
-import org.lcsim.event.EventHeader;
-import org.lcsim.event.MCParticle;
-import org.lcsim.geometry.IDDecoder;
-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;
-
-/**
- *
- * @author carsten
- */
-public class FixedConeClustererAnalyzer extends Driver{
- private AIDA aida = AIDA.defaultInstance();
- private ITree tree;
- private static final String treePath = "FixedConeClustererAnalyzer";
- private DecimalFormat df;
- private double SAMP_FRAC;
- private double radiusCut;
- private int debug;
- private int eventCounter = 0;
- private boolean initialized = false;
- private FixedConeClusterer fcClusterer;
- //private CalorimeterIDDecoder decoder;
- private IDDecoder decoder;
-
- /** Creates a new instance of FixedConeClustererAnalyzer */
- public FixedConeClustererAnalyzer() {
- }
-
- public FixedConeClustererAnalyzer(ITree tree) {
- this();
- this.tree = tree;
- }
-
- public void setAttributes(double samplingFraction, int debug, double radiusCut){
- this.SAMP_FRAC = samplingFraction;
- this.debug = debug;
- this.radiusCut = radiusCut;
- }
-
- public boolean eventPassed(List<MCParticle> mcParticles) {
- boolean passed = false;
- for (MCParticle mcParticle : mcParticles) {
- if (mcParticle.getGeneratorStatus() == MCParticle.FINAL_STATE) {
-
- /*
- * Rejects photons with radius < 1260. Here, radius is not an
- * angle but is an actual radius.
- *
- * 'pre-selection'
- */
-
- if (mcParticle.getType().getName().equals("gamma")
- && mcParticle.getEndPoint().magnitude() < this.radiusCut) {
- passed = false;
- } else {
- passed = true;
- }
- }
- }
- return passed;
- }
-
- public void initialize() {
- this.tree.cd("/");
- this.mkdirIgnoreExist(this.tree, this.treePath);
- df = new DecimalFormat();
- df.setMaximumFractionDigits(2);
- df.setMinimumFractionDigits(2);
- this.initialized = true;
- }
-
- public void process(EventHeader event) {
-
- if (!this.initialized)
- this.initialize();
-
- List<MCParticle> mcParticles = event.getMCParticles();
- if (!this.eventPassed(mcParticles))
- return;
-
- this.tree.cd("/" + this.treePath);
-
- List<CalorimeterHit> ecalBarrHits = event.get(CalorimeterHit.class,"EcalBarrHits");
- List<CalorimeterHit> endcapHits = event.get(CalorimeterHit.class,"EcalEndcapHits");
-
- double eMax = 0.0;
- for (CalorimeterHit hit : endcapHits)
- eMax += hit.getRawEnergy();
-
- for (CalorimeterHit hit : ecalBarrHits)
- eMax += hit.getRawEnergy();
-
-
-
- this.decoder = event.getMetaData(ecalBarrHits).getIDDecoder();
-
- // create a map of cells keyed on their index
- Map<Long, CalorimeterHit> hitMap = new HashMap<Long, CalorimeterHit>();
-
- double initSeed = 0.0;
- double initEmin = 0.005;
- for (double initRadius = 0.01; initRadius < 0.1; initRadius += 0.01) {
- // set the tree path
- StringBuffer buf = new StringBuffer();
- buf.delete(0, buf.length());
- buf.insert(0, "r_");
- df.format(initRadius, buf, new FieldPosition(df.INTEGER_FIELD));
- this.mkdirIgnoreExist(this.tree, buf.toString());
-
- this.fcClusterer = new FixedConeClusterer(initRadius, initSeed, initEmin);
-
- List<Cluster> clusters = fcClusterer.createClusters(ecalBarrHits);
- for (Cluster cluster : clusters) {
- if(cluster.getCalorimeterHits().size() > 5) { // number of cells
- aida.cloud1D("cluster energy over total energy").fill(((BasicCluster)cluster).getRawEnergy() / eMax);
- }
- }
- this.tree.cd("/" + this.treePath);
- }
- }
-
-
-
-
- protected void endOfData() {
- System.out.println(this.eventCounter + " events analyzed by FCCA");
- }
- // 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
- // "dirname" if it is found. If not, a runtime error results.
- //
- public void mkdirIgnoreExist(ITree tree, String dirname) {
- try {
- tree.mkdir(dirname);
- } catch (IllegalArgumentException e) {
- if (! (e instanceof TreeObjectAlreadyExistException))
- throw e;
- }
- tree.cd(dirname);
- }
-
-}
lcsim/src/org/lcsim/contrib/CarstenHensel
diff -N HMatrixAnalyzer.java
--- HMatrixAnalyzer.java 14 Feb 2007 00:16:04 -0000 1.3
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,358 +0,0 @@
-/*
- * HMatrixAnalyzer.java
- *
- * Created on August 22, 2005, 9:53 AM
- *
- * To change this template, choose Tools | Options and locate the template under
- * the Source Creation and Management node. Right-click the template and choose
- * Open. You can then make changes to the template in the Source Editor.
- */
-
-package org.lcsim.contrib.CarstenHensel;
-import hep.aida.IAnalysisFactory;
-import hep.aida.ITree;
-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;
-import org.lcsim.conditions.ConditionsManager;
-import org.lcsim.conditions.ConditionsManager.ConditionsNotFoundException;
-import org.lcsim.event.CalorimeterHit;
-import org.lcsim.event.Cluster;
-import org.lcsim.event.EventHeader;
-import org.lcsim.event.MCParticle;
-import org.lcsim.geometry.CalorimeterIDDecoder;
-import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
-import org.lcsim.math.chisq.ChisqProb;
-import org.lcsim.recon.cluster.fixedcone.FixedConeClusterer;
-import org.lcsim.recon.cluster.util.BasicCluster;
-import org.lcsim.recon.emid.hmatrix.HMatrix;
-import org.lcsim.recon.emid.hmatrix.HMatrixConditionsConverter;
-import org.lcsim.util.Driver;
-import org.lcsim.util.aida.AIDA;
-
-
-/**
- *
- * @author carsten
- */
-public class HMatrixAnalyzer extends Driver {
-
- private AIDA aida = AIDA.defaultInstance();
- private ITree tree;
- private static final String treePath = "HMatrixAnalyzer";
- private HMatrix hMatrix;
- private String hMatrixInputFileName = "";
- private int eventCounter;
- private ConditionsManager mgr = ConditionsManager.defaultInstance();
- private FixedConeClusterer fcClusterer;
- private double SAMP_FRAC = 0.012;
- private double radiusCut;
- private int nLayers;
- private CalorimeterIDDecoder decoder;
- 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) {
- this.tree = tree;
- this.readHMatrix();
- this.setAttributes();
- }
-
-
-
- public HMatrixAnalyzer(String inputFileName, ITree tree) {
- this.tree = tree;
- this.hMatrixInputFileName = inputFileName;
- this.readHMatrix();
- this.setAttributes();
- }
-
- 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", "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;
- }
-
- private void setAttributes(){
-
- double initRadius = 0.06;
- double initSeed = 0.0;
- double initEmin = 0.005;
- this.fcClusterer = new FixedConeClusterer(initRadius, initSeed, initEmin);
- this.radiusCut = 1260.0;
- }
-
- public void process(EventHeader event) {
- // set remaining attribitues which can't be set in initialize method...
- if (!this.initialized)
- this.initialize(event);
-
- List<MCParticle> mcParticles = event.getMCParticles();
- // check event 'pre-selection'
- if (!this.eventPassed(mcParticles)) {
- return;
- }
-
- this.tree.cd("/" + this.treePath);
-
- this.eventCounter++;
- List<CalorimeterHit> ecalBarrHits = event.get(CalorimeterHit.class,"EcalBarrHits");
- List<CalorimeterHit> endcapHits = event.get(CalorimeterHit.class,"EcalEndcapHits");
- this.decoder = (CalorimeterIDDecoder) event.getMetaData(ecalBarrHits).getIDDecoder();
-
- // create a map of cells keyed on their index
- Map<Long, CalorimeterHit> hitMap = new HashMap<Long, CalorimeterHit>();
-
- double sumBarrHits = 0.;
- double sumEndcapHits = 0.; // NOTE: Only nonzero for neutrons.
-
-
-
-
- List<Cluster> clusters = fcClusterer.createClusters(ecalBarrHits);
-
- // Get all the clusters and cluster sets to prepare to loop over them.
- List<List<Cluster>> clusterSets = event.get(Cluster.class);
-
-
- // Accumulators for a raw (uncorrected for sampling fraction)
- // momentum 4-vector sum, a raw energy sum, and a corrected (for
- // sampling fraction) energy sum
-
-
-
- // Loop through all clusters
- for (Cluster cluster : clusters) {
-
-
- if(cluster.getCalorimeterHits().size() > 5) { // number of cells
-
- 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 = ((BasicCluster)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);
- aida.cloud1D("chisq diag probability").fill(probDiag);
- aida.cloud2D("chisq probability vs chisq diag probability").fill(prob, probDiag);
- aida.cloud2D("chisq vs energy").fill(clusterRawEnergy,chisq);
- aida.cloud2D("prob vs energy").fill(clusterRawEnergy,chisq);
- this.tuple.fill(0, chisq);
- this.tuple.fill(1, chisqDiag);
- this.tuple.fill(2, prob);
- 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();
- }
-
-
-
- }
-
- }
-
- protected void endOfData() {
- System.out.println(this.eventCounter + " events analyzed by HMatrixAnalyzer");
- }
-
- private boolean eventPassed(List<MCParticle> mcParticles) {
- boolean passed = false;
- for (MCParticle mcParticle : mcParticles) {
- if (mcParticle.getGeneratorStatus() == MCParticle.FINAL_STATE) {
-
- /*
- * Rejects photons with radius < 1260. Here, radius is not an
- * angle but is an actual radius.
- *
- * 'pre-selection'
- */
-
- if (mcParticle.getType().getName().equals("gamma")
- && mcParticle.getEndPoint().magnitude() < this.radiusCut) {
- passed = false;
- } else {
- passed = true;
- }
- }
- }
- return passed;
- }
-
-
-
- private double[] getListOfLayerEnergies(Cluster cluster) {
- //System.out.println("Layers " + this.nLayers + " " + this.hMatrix.toString());
- double[] energyPerLayer = new double[this.nLayers];
- double clusterEnergy = 0.0;
- List<CalorimeterHit> hits = cluster.getCalorimeterHits();
- for(CalorimeterHit hit : hits){
- this.decoder.setID(hit.getCellID());
- double e = hit.getRawEnergy();
- int layer = this.decoder.getLayer();
- clusterEnergy+=e;
- energyPerLayer[layer]+=e;
- }
- for (int i = 0; i < energyPerLayer.length; i++) {
- energyPerLayer[i] /= clusterEnergy;
- }
- return energyPerLayer;
-
- }
-
-
-
-
- private void readHMatrix(){
- if (this.hMatrixInputFileName == "") {
- //read default HMatrix
- System.out.println("read default HMatrix");
- try {
- this.mgr.setDetector("cdcaug05", 0);
- } catch (ConditionsNotFoundException e) {
- System.out.println("detector not found");
- throw new RuntimeException("Please see http://confluence.slac.stanford.edu/display/ilc/Conditions+database", e);
- }
- this.mgr.registerConditionsConverter(new HMatrixConditionsConverter());
- this.hMatrix = this.mgr.getCachedConditions(HMatrix.class, "LongitudinalHMatrix.hmx").getCachedData();
- } else {
- this.hMatrix = HMatrix.read(this.hMatrixInputFileName);
- }
- }
-
- private double[] getHMInputVals(Cluster 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
- // "dirname" if it is found. If not, a runtime error results.
- //
- public void mkdirIgnoreExist(ITree tree, String dirname) {
- try {
- tree.mkdir(dirname);
- } catch (IllegalArgumentException e) {
- if (! (e instanceof TreeObjectAlreadyExistException))
- throw e;
- }
- tree.cd(dirname);
- }
-
-
-
-}
lcsim/src/org/lcsim/contrib/CarstenHensel
diff -N HMatrixAnalyzerPlotter.java
--- HMatrixAnalyzerPlotter.java 20 Feb 2007 20:41:11 -0000 1.3
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,173 +0,0 @@
-/*
- * HMAtrixAnalyzerPlotter.java
- *
- * Created on August 23, 2005, 2:44 PM
- *
- * To change this template, choose Tools | Options and locate the template under
- * the Source Creation and Management node. Right-click the template and choose
- * Open. You can then make changes to the template in the Source Editor.
- */
-
-package org.lcsim.contrib.CarstenHensel;
-
-import hep.aida.*;
-import org.freehep.application.Application;
-import org.freehep.application.studio.Studio;
-import org.lcsim.util.aida.AIDA;
-
-
-/**
- *
- * @author carsten
- */
-public class HMatrixAnalyzerPlotter {
-
-
- private AIDA aida = AIDA.defaultInstance();
- 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() {
- }
-
-
-
-
-
-
- public static void main(String[] argv) {
-
- HMatrixAnalyzerPlotter plotter = new HMatrixAnalyzerPlotter();
- IAnalysisFactory af = IAnalysisFactory.create();
- ITree tree = af.createTreeFactory().create();
- 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};
-
-
- for (int iTuple = 0; iTuple < tupleHolder.length; iTuple++) {
- 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 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[nProbVals];
- tuple.start();
- for (int iRow = 0; iRow < nTupleRows; iRow++) {
- tuple.next();
- for (int iProb = 0; iProb < nProbVals; iProb++) {
- if (tuple.getDouble(rowNumber) > -1.0 * iProb) {
- integral[iProb] += weight;
- }
- }
- }
- return integral;
- }
-
-}
-
-
-
lcsim/src/org/lcsim/contrib/CarstenHensel
diff -N MainLoop.java
--- MainLoop.java 11 Sep 2007 00:21:03 -0000 1.6
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,43 +0,0 @@
-/*
- * MainLoop.java
- *
- * Created on August 21, 2005, 3:59 PM
- *
- * Java wrapper to enable running outside of JAS3
- * 16-JUL-2005 Jan Strube
- * from a response to the JAS mailing list by Tony Johnson
- */
-
-package org.lcsim.contrib.CarstenHensel;
-
-/**
- *
- * @author carsten
- */
-
-import java.io.File;
-
-import org.lcsim.util.Driver;
-import org.lcsim.util.aida.AIDA;
-import org.lcsim.util.loop.LCSimLoop;
-
-public class MainLoop extends Driver{
-
- /** Creates a new instance of MainLoop */
- public MainLoop() {
- }
-
-
- public static void main(String[] args) throws Exception{
- LCSimLoop loop = new LCSimLoop();
- 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(-1, null);
- loop.dispose();
- AIDA.defaultInstance().saveAs("/home/carsten/LC/Output/exampleAnalysisJava.aida");
- }
-}
-
lcsim/src/org/lcsim/contrib/CarstenHensel
diff -N MyHMatrixBuilder.java
--- MyHMatrixBuilder.java 14 Feb 2007 00:16:04 -0000 1.3
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,251 +0,0 @@
-/*
- * MyHMatrixBuilder.java
- *
- * Created on August 22, 2005, 6:11 PM
- *
- * To change this template, choose Tools | Options and locate the template under
- * the Source Creation and Management node. Right-click the template and choose
- * Open. You can then make changes to the template in the Source Editor.
- */
-
-package org.lcsim.contrib.CarstenHensel;
-
-import java.text.DecimalFormat;
-import java.util.Calendar;
-import java.util.Date;
-import java.util.GregorianCalendar;
-import java.util.List;
-
-import org.lcsim.event.CalorimeterHit;
-import org.lcsim.event.Cluster;
-import org.lcsim.event.EventHeader;
-import org.lcsim.event.MCParticle;
-import org.lcsim.geometry.IDDecoder;
-import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
-import org.lcsim.recon.cluster.fixedcone.FixedConeClusterer;
-import org.lcsim.recon.cluster.util.BasicCluster;
-import org.lcsim.recon.emid.hmatrix.HMatrixBuilder;
-import org.lcsim.util.Driver;
-
-/**
- *
- * @author carsten
- */
-public class MyHMatrixBuilder extends Driver{
-
- 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 IDDecoder decoder;
- private static final double log10inv = 1./Math.log(10.0);
- private int debug = 0;
-
- /** Creates a new instance of 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);
-
-
- List<MCParticle> mcParticles = event.getMCParticles();
-
- double thetaMC = 0;
- double cosThetaMC = 0;
- double phiMC = 0;
- double p2 = 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.atan(pT/pZ);
- //thetaMC = Math.atan2(pT, pZ);
- thetaMC = Math.acos(pZ/Math.sqrt(p2));
- cosThetaMC = pZ / Math.sqrt(p2);
-
-
-
- // Rejects photons with radius < 1260. Here, radius is not an
- // angle but is an actual radius.
- //
- if (mcParticle.getType().getName().equals("gamma")
- && mcParticle.getEndPoint().magnitude() < 1260.0) {
- return;
- }
- }
- }
-
-
- //FIXME: need to get the EM calorimeter 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 = event.getMetaData(collection).getIDDecoder();
-
-
-
- // construct the EM Clusters
- // in the future will fetch from the event...
- //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);
- //cluster
- List<Cluster> clusters = fcc.createClusters(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 cluster : clusters) {
- int ncells = cluster.getCalorimeterHits().size();
- //FIXME should cut on cluster corrected energy
- if(ncells>5) {
- nGoodClusters++;
- //GWW double energy = cluster.getEnergy();
- double energy = ((BasicCluster)cluster).getRawEnergy();
- // should be able to fetch this from the cluster...
- this.vals = getHMInputVals(cluster);
-
-
- this.hmb.accumulate(this.vals);
-/*
- StringBuffer sb = new StringBuffer();
- for(int i=0; i<_vals.length; ++i)
- {
- sb.append(_myFormatter.format(_vals[i])+" ");
- }
- System.out.println(sb);
- */
-
-
-
- }
-
- }
- }
-
- protected void endOfData() {
- this.hmb.validate();
- String fileLocation = System.getProperty("HMATRIX_FILE", "CH.hmx");
- System.out.println("Writing out HMatrix to "+fileLocation);
- this.hmb.write("/home/carsten/LC/Output/CH.hmx", this.commentForHMatrix());
-
- }
-
- private double[] getHMInputVals(Cluster 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;
- }
-
-
-
- private String commentForHMatrix() {
- 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 this.detectorName+" "+myDate+" "+System.getProperty("user.name");
- }
-
-}
\ No newline at end of file
CVSspam 0.2.8
