lcsim-cal-calib/src/org/lcsim/cal/calib
diff -N EMClusterID.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ EMClusterID.java 5 Jun 2008 05:27:03 -0000 1.1
@@ -0,0 +1,316 @@
+package org.lcsim.cal.calib;
+import java.util.List;
+import java.util.ArrayList;
+import org.lcsim.event.EventHeader;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.event.CalorimeterHit;
+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 java.util.Calendar;
+import java.util.Date;
+import java.util.GregorianCalendar;
+import org.lcsim.conditions.ConditionsManager;
+import org.lcsim.conditions.ConditionsManager.ConditionsSetNotFoundException;
+import org.lcsim.conditions.ConditionsSet;
+import org.lcsim.event.Cluster;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.recon.emid.hmatrix.HMatrix;
+import org.lcsim.recon.emid.hmatrix.HMatrixTask;
+import org.lcsim.recon.emid.hmatrix.HMatrixBuilder;
+import org.lcsim.recon.emid.hmatrix.HMatrixConditionsConverter;
+import org.lcsim.math.chisq.ChisqProb;
+import org.lcsim.math.moments.CentralMomentsCalculator;
+import org.lcsim.recon.cluster.fixedcone.FixedConeClusterer;
+import org.lcsim.recon.cluster.fixedcone.FixedConeClusterer.FixedConeDistanceMetric;
+
+/**
+ * Reconstruction: EM Clusters
+ */
+public class EMClusterID extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+ private ITree _tree;
+ private boolean _initialized;
+ private ConditionsSet _cond;
+ private String _detName;
+ private HMatrixTask _task;
+ private int _nLayers;
+
+ private HMatrixBuilder _hmb;
+ private HMatrix _hmx;
+
+ // the number of variables in the measurement vector
+ private int _nmeas;
+
+ // the vector of measured values
+ double[] _vals;
+
+ // the mapping of physical layers to measurement space
+ double[] _layerMapping;
+
+ private DecimalFormat _myFormatter = new DecimalFormat("#.###");
+
+ private boolean _debug = true;
+
+ // where to write the HMatrix if in accumulate mode
+ private String _fileLocation = "default.hmx";
+
+ public EMClusterID()
+ {
+ this(HMatrixTask.ANALYZE);
+// this(HMatrixTask.BUILD);
+ }
+
+ public EMClusterID(HMatrixTask task)
+ {
+ _tree = aida.tree();
+ _task = task;
+
+ if(task==HMatrixTask.ANALYZE)
+ {
+ // The HMatrix could possibly change, so be sensitive to this.
+ getConditionsManager().registerConditionsConverter(new HMatrixConditionsConverter());
+ }
+
+ }
+
+ protected void process(EventHeader event)
+ {
+ super.process(event);
+
+ //FIXME: need to get the EM calorimeter names from a conditions file
+ // FIXME should get calorimeterhit collection names from a conditions file
+
+ String[] det = {"EMBarrel","EMEndcap"};
+ String[] hitsToGet = {"EcalBarrHits","EcalEndcapHits"};
+
+ if(!_initialized)
+ {
+ ConditionsManager mgr = ConditionsManager.defaultInstance();
+ try
+ {
+ _cond = mgr.getConditions("HMatrix");
+ }
+ catch(ConditionsSetNotFoundException e)
+ {
+ System.out.println("ConditionSet HMatrix not found for detector "+mgr.getDetector());
+ System.out.println("Please check that this properties file exists for this detector ");
+ }
+ // sanity check
+ String detectorNameFromFile = _cond.getString("detectorName");
+ if(!detectorNameFromFile.equals(event.getDetectorName()))
+ {
+ System.out.println("detector name from HMatrix.properties: "+detectorNameFromFile +" detector name from event "+event.getDetectorName());
+ throw new RuntimeException("detector name mismatch in HMatrix!");
+ }
+ // add the detector name to the output file
+ _fileLocation = event.getDetectorName()+"_"+_fileLocation;
+ // the vector of measurements starts as the longitudinal layers
+ _nmeas = _cond.getInt("measurementDimension");
+ // would add any additional measurements (e.g. width) here
+ _vals = new double[_nmeas];
+
+ _layerMapping = _cond.getDoubleArray(_nmeas+"layerMapping");
+// for (int i=0; i<_layerMapping.length; ++i)
+// {
+// System.out.println("_layerMapping[ "+i+" ]= "+ _layerMapping[i]);
+// }
+ // sanity check
+
+ CylindricalCalorimeter calsub = (CylindricalCalorimeter)event.getDetector().getSubdetectors().get(det[0]);
+ _nLayers = calsub.getLayering().getLayerCount();
+ if(_nLayers != _layerMapping.length)
+ {
+ System.out.println("found "+_nLayers+" layers in the "+det[0]);
+ throw new RuntimeException("layer number mismatch in EMCalorimeter!");
+ }
+
+ //FIXME key needs to be better defined
+ int key = 0;
+ if(_task==HMatrixTask.ANALYZE)
+ {
+ //FIXME need to fetch name of HMAtrix file to use from a conditions file
+ _hmx = getConditionsManager().getCachedConditions(HMatrix.class, "LongitudinalHMatrix"+_nmeas+".hmx").getCachedData();
+ }
+ else if(_task==HMatrixTask.BUILD)
+ {
+ _hmb = new HMatrixBuilder(_nmeas,key);
+ }
+ _detName = event.getDetectorName();
+
+
+ _initialized = true;
+ }
+
+ List<Cluster> photons = new ArrayList<Cluster>();
+ for(int j=0; j<det.length; ++j)
+ {
+ List<CalorimeterHit> collection = event.get(CalorimeterHit.class,hitsToGet[j]);
+
+
+ double radius = 0.5;
+ double seed = 0.;
+ double minE = 0.05;
+
+ // create the clusterer
+ FixedConeClusterer fcc = new FixedConeClusterer(radius, seed, minE, FixedConeDistanceMetric.DPHIDTHETA);
+ //cluster
+ List<Cluster> clusters = fcc.createClusters(collection);
+
+ int minHitsInCluster = 20;
+
+ // add this list of clusters to the event (for event display)
+ event.put(hitsToGet[j]+"Clusters ",clusters);
+ // Loop over all the clusters
+ int nGoodClusters = 0;
+ for (Cluster cluster : clusters)
+ {
+ int ncells = cluster.getCalorimeterHits().size();
+ double energy = cluster.getEnergy();
+ //FIXME should cut on cluster corrected energy
+ if(ncells > minHitsInCluster)
+ {
+ nGoodClusters++;
+
+ aida.cloud1D("Number of cells in cluster").fill(ncells);
+ aida.cloud2D("Number of cells in cluster vs cluster energy").fill(energy, ncells);
+
+ aida.cloud2D("Hottest cell in cluster vs cluster energy").fill(energy,hottestCellEnergy(cluster));
+// aida.cloud1D("Cluster raw energy").fill(bc.getRawEnergy() );
+ aida.cloud1D("ClusterCorrected energy").fill(energy);
+ // should be able to fetch this from the cluster...
+ double[] layerE = layerFractionalEnergies(cluster);
+ // accumulate the longitudinal energy fractions into the measurement vector...
+ for(int i=0; i<layerE.length; ++i)
+ {
+ if(_layerMapping[i] != -1)
+ {
+ _vals[(int)_layerMapping[i]] += layerE[i];
+// System.out.println("LayerE[ "+i+" ]= "+layerE[i]+" _vals[ "+(int)_layerMapping[i]+" ] "+_vals[(int)_layerMapping[i]]);
+ aida.cloud2D("Fractional Energy vs physical Layer").fill(i,layerE[i]);
+ aida.cloud2D("Fractional Energy vs mapped Layer").fill(_layerMapping[i],layerE[i]);
+
+ }
+ }
+
+ for(int i=0; i<_vals.length; ++i)
+ {
+ aida.cloud1D("Fractional Energy for mapped Layer "+i).fill(_vals[i]);
+ for(int k=i+1; k<_vals.length; ++k)
+ {
+ aida.cloud2D("Fractional Energy "+i+" vs "+k).fill(_vals[i], _vals[k]);
+ }
+ }
+
+ // have now filled the vector of measurements. need to either accumulate the HMatrix or apply it
+ if (_task==HMatrixTask.BUILD)
+ {
+ _hmb.accumulate(_vals);
+ }
+ if (_task==HMatrixTask.ANALYZE)
+ {
+ double chisq = _hmx.chisquared(_vals);
+ aida.cloud1D("Chisq").fill(chisq);
+ aida.cloud2D("Chisq vs energy").fill(energy,chisq);
+ if(chisq<500)
+ {
+ aida.cloud1D("Chisq(<500)").fill(chisq);
+ aida.cloud2D("Chisq (<500) vs energy").fill(energy,chisq);
+ }
+ aida.cloud1D("Chisq Probability").fill(ChisqProb.gammq(_nmeas,chisq));
+
+ double chisqD = _hmx.chisquaredDiagonal(_vals);
+ aida.cloud1D("ChisqD").fill(chisqD);
+ aida.cloud2D("ChisqD vs energy").fill(energy,chisqD);
+
+ double chisqDProb = ChisqProb.gammq(_nmeas,chisqD);
+ if(chisqDProb<.0000000001) chisqDProb = 0.0000000001;
+ aida.cloud1D("ChisqD Probability").fill(chisqDProb);
+ aida.cloud1D("log10 ChisqDProb").fill(Math.log10(chisqDProb));
+ }
+ 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]) );
+ //
+ // reset measurement vector...
+ //
+ for(int i=0; i<_vals.length; ++i)
+ {
+ _vals[i]=0.;
+ }
+ }// end over loop over clusters with greater tna minHits
+ }// end of loop over clusters
+ aida.cloud1D(det[j]+ "number of found clusters (above hits threshold)").fill(nGoodClusters);
+ }// end of loop over collections
+ event.put("photons",photons);
+ }
+
+ private double[] layerFractionalEnergies(Cluster clus)
+ {
+ //FIXME could reuse this array
+ double[] layerEnergies = new double[_nLayers];
+ double clusterEnergy = 0.;
+ List<CalorimeterHit> hits = clus.getCalorimeterHits();
+ for(CalorimeterHit hit : hits)
+ {
+ IDDecoder decoder = hit.getIDDecoder();
+ decoder.setID(hit.getCellID());
+ double e = hit.getCorrectedEnergy();
+ int layer = decoder.getLayer();
+// System.out.println("layer "+layer+" energy "+e);
+ clusterEnergy+=e;
+ layerEnergies[layer]+=e;
+ }
+ for(int i=0; i<_nLayers; ++i)
+ {
+ layerEnergies[i]/=clusterEnergy;
+// System.out.println("i= "+i+" layerEnergies= "+layerEnergies[i]);
+ }
+// System.out.println(clusterEnergy+" "+clus.getEnergy());
+ return layerEnergies;
+ }
+
+ private double hottestCellEnergy(Cluster clus)
+ {
+ double hottestCellEnergy = 0.;
+ List<CalorimeterHit> hits = clus.getCalorimeterHits();
+// System.out.println("New cluster with "+hits.size()+ " hits and energy "+clus.getEnergy());
+ for(CalorimeterHit hit : hits)
+ {
+ double e = hit.getCorrectedEnergy();
+ if(e>hottestCellEnergy) hottestCellEnergy=e;
+ }
+
+ return hottestCellEnergy;
+ }
+
+ protected void endOfData()
+ {
+ if (_task==HMatrixTask.BUILD)
+ {
+ _hmb.validate();
+ _hmb.write(_fileLocation,commentForHMatrix());
+ }
+ }
+
+ public void setHMatrixFileLocation(String filename)
+ {
+ _fileLocation = filename;
+ }
+
+ 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 _detName+" "+myDate+" "+System.getProperty("user.name");
+ }
+}
lcsim-cal-calib/src/org/lcsim/cal/calib
diff -N StandaloneEMClusterAnalysis.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ StandaloneEMClusterAnalysis.java 5 Jun 2008 05:27:03 -0000 1.1
@@ -0,0 +1,85 @@
+package org.lcsim.cal.calib;
+import java.io.BufferedReader;
+import org.lcsim.recon.emid.hmatrix.HMatrixTask;
+import java.io.File;
+import java.io.FileInputStream;
+import java.io.InputStreamReader;
+import java.util.ArrayList;
+import java.util.List;
+import org.lcsim.util.loop.LCIOEventSource;
+import org.lcsim.util.loop.LCSimLoop;
+
+/**
+ * A simple standalone EMClusterID example
+ */
+public class StandaloneEMClusterAnalysis
+{
+ public static void main(String[] args) throws Exception
+ {
+ if(args.length<2)
+ {
+ usage();
+ return ;
+ }
+
+ String listOfFiles = args[0];
+ List<File> filesToProcess = filesToProcess(listOfFiles);
+ LCIOEventSource src = new LCIOEventSource("EMClusterAnalysis", filesToProcess);
+
+ String task = args[1];
+ if(!(task.equals("analyze") || task.equals("build")))
+ {
+ System.out.println("'"+args[1] +"' not a recognized task. Please specify 'analyze' or 'build'");
+ return;
+ }
+ int numToProcess=-1;
+ if(args.length>2) numToProcess=Integer.parseInt(args[2]);
+
+ System.out.println("Processing "+numToProcess+" events from "+listOfFiles);
+ for(File f : filesToProcess)
+ {
+ String HMatrixName = "";
+ String[] parts = f.getName().split("_");
+ for(String s : parts)
+ {
+ if(s.startsWith("Theta")) HMatrixName+=s+"_";
+ if(s.contains("GeV")) HMatrixName+=s;
+ }
+ HMatrixName+=".hmx";
+
+ LCSimLoop loop = new LCSimLoop();
+ loop.setLCIORecordSource(f);
+ HMatrixTask taskType = HMatrixTask.ANALYZE;
+ if(task.equals("build")) taskType = HMatrixTask.BUILD;
+ EMClusterID emClusID = new EMClusterID(taskType);
+ emClusID.setHMatrixFileLocation(HMatrixName);
+ loop.add(emClusID);
+ loop.loop(numToProcess);
+ loop.dispose();
+ }
+ }
+
+ public static void usage()
+ {
+ System.out.println("This is StandaloneEMClusterAnalysis");
+ System.out.println("usage:");
+ System.out.println("java StandaloneEMClusterAnalysis listOfInputFiles build/analyze [number of events to process]");
+ }
+
+ public static List<File> filesToProcess(String listOfFiles) throws Exception
+ {
+ List<File> filesToProcess = new ArrayList<File>();
+ FileInputStream fin = new FileInputStream(listOfFiles);
+ BufferedReader br = new BufferedReader(new InputStreamReader(fin));
+ String line;
+
+ while ( (line = br.readLine()) != null)
+ {
+ File f = new File(line.trim());
+ if(!f.exists()) throw new RuntimeException("Input file "+f+ " does not exist!");
+ filesToProcess.add(f);
+ }
+
+ return filesToProcess;
+ }
+}
