lcsim/src/org/lcsim/contrib/CarstenHensel
diff -N MyHMatrixBuilder.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MyHMatrixBuilder.java 24 Aug 2005 15:11:45 -0000 1.1
@@ -0,0 +1,244 @@
+/*
+ * 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.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.geometry.subdetector.CalorimeterIDDecoder;
+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 int logEIndex;
+ private double[] vals;
+ private HMatrixBuilder hmb;
+ private String detectorName;
+ private CalorimeterIDDecoder decoder;
+ private static final double log10inv = 1./Math.log(10.0);
+
+ /** Creates a new instance of MyHMatrixBuilder */
+ public MyHMatrixBuilder() {
+ }
+
+
+ protected void process(EventHeader 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;
+
+ 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
+ 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
+ List<CalorimeterHit> collection = event.get(CalorimeterHit.class,"EcalBarrHits");
+ this.decoder = (CalorimeterIDDecoder) 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<BasicCluster> clusters = fcc.makeClusters(collection, this.decoder);
+
+ // 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 (BasicCluster 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 = 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.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);
+ */
+
+
+
+ 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;
+
+
+
+
+ }
+
+ }
+ }
+
+ 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[] layerEnergies(BasicCluster clus) {
+ //FIXME could reuse this array
+ double[] layerEnergies = new double[this.nLayers];
+ 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();
+ //System.out.println("layer "+layer+" energy "+e);
+ clusterEnergy+=e;
+ layerEnergies[layer]+=e;
+ }
+
+ //System.out.println("Cluster energy = " + clusterEnergy + "\n");
+
+ 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());
+ 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