lcsim/src/org/lcsim/contrib/CarstenHensel
diff -N HMatrixAnalyzer.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HMatrixAnalyzer.java 24 Aug 2005 15:09:26 -0000 1.1
@@ -0,0 +1,264 @@
+/*
+ * 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.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;
+
+ /** 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.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, "");
+
+ 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<BasicCluster> clusters = fcClusterer.makeClusters(ecalBarrHits, decoder);
+
+ // 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 (BasicCluster cluster : clusters) {
+
+
+ if(cluster.getCalorimeterHits().size() > 5) { // number of cells
+
+ double[] val = new double[this.nLayers + 1];
+ double[] energiesInLayer = this.getListOfLayerEnergies(cluster);
+ for (int i = 0; i < energiesInLayer.length; i++) {
+ val[i] = energiesInLayer[i];
+ aida.cloud2D("Fractional Energy vs Layer").fill(i, energiesInLayer[i]);
+ }
+
+ double clusterRawEnergy = cluster.getRawEnergy();
+ aida.cloud1D("cluster raw energy").fill(clusterRawEnergy);
+
+ val[this.nLayers] = Math.log(clusterRawEnergy) / Math.log(10.0);;
+
+ double chisq = this.hMatrix.chisquared(val);
+ double chisqDiag = this.hMatrix.chisquaredDiagonal(val);
+ double prob = ChisqProb.gammq(this.nLayers + 1, chisq);
+ double probDiag = ChisqProb.gammq(this.nLayers + 1, chisqDiag);
+ 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, clusterRawEnergy);
+ 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(BasicCluster 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);
+ }
+ }
+
+ // 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);
+ }
+
+
+
+}