Print

Print


Commit in lcsim-cal-calib on MAIN
build.sh+17added 1.1
maven.xml+8added 1.1
pom.xml+83added 1.1
project.properties+22added 1.1
project.xml+187added 1.1
src/org/lcsim/cal/calib/Calibrate.java+98added 1.1
                       /SamplingFractionAnalysisDriver.java+298added 1.1
+713
7 added files
First implementation of calorimeter calibration, with sampling fractions

lcsim-cal-calib
build.sh added at 1.1
diff -N build.sh
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ build.sh	19 May 2008 21:31:58 -0000	1.1
@@ -0,0 +1,17 @@
+#!/bin/sh
+
+checkcyg=$(uname | grep CYGWIN)
+if [ -n "$checkcyg" ]
+then 
+  install_dir=$(cygpath -w $(pwd))
+else
+  install_dir=$(pwd)
+fi
+
+echo "install_dir=$install_dir"
+
+maven -Drun.install=${install_dir} \
+      -Dmaven.test.skip=true clean \
+      jas:install jar:install run:install
+
+#./test-single.sh org.lcsim.slic.diagnostics.SlicDiagnosticsDriverTest

lcsim-cal-calib
maven.xml added at 1.1
diff -N maven.xml
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ maven.xml	19 May 2008 21:31:58 -0000	1.1
@@ -0,0 +1,8 @@
+<?xml version="1.0" encoding="UTF-8"?>
+
+<project default="jar">
+<postGoal name="xdoc:register-reports">
+  <attainGoal name="maven-checkstyle-plugin:deregister"/>  
+  <attainGoal name="maven-license-plugin:deregister"/>
+</postGoal>
+</project>

lcsim-cal-calib
pom.xml added at 1.1
diff -N pom.xml
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ pom.xml	19 May 2008 21:31:58 -0000	1.1
@@ -0,0 +1,83 @@
+<?xml version="1.0" encoding="UTF-8"?><project>
+  <modelVersion>4.0.0</modelVersion>
+  <groupId>org.lcsim</groupId>
+  <artifactId>lcsim-cal-calib</artifactId>
+  <name>lcsim Calorimeter Calibration</name>
+  <version>0.1-SNAPSHOT</version>
+  <description>An org.lcsim based calorimeter calibration package.</description>
+  <build>
+    <sourceDirectory>src</sourceDirectory>
+    <defaultGoal>install</defaultGoal>
+    <testSourceDirectory>test</testSourceDirectory>
+    <testResources>
+      <testResource>
+        <directory>test</directory>
+      </testResource>
+    </testResources>
+    <pluginManagement>
+      <plugins>
+        <plugin>
+          <artifactId>maven-compiler-plugin</artifactId>
+          <configuration>
+            <excludes>
+              <exclude>org/lcsim/contrib/**</exclude>
+            </excludes>
+            <source>1.5</source>
+            <target>1.5</target>
+          </configuration>
+        </plugin>
+      </plugins>
+    </pluginManagement>
+  </build>
+  <repositories>
+    <repository>
+      <id>freehep-maven</id>
+      <name>Maven FreeHEP</name>
+      <url>http://java.freehep.org/maven2</url>
+    </repository>
+  </repositories>
+  <dependencies>
+    <dependency>
+      <groupId>org.freehep</groupId>
+      <artifactId>freehep-jaida</artifactId>
+      <version>3.3.0-5</version> 
+    </dependency>
+    <dependency>
+      <groupId>org.lcsim</groupId>
+      <artifactId>lcsim</artifactId>
+      <version>1.2-SNAPSHOT</version>
+    </dependency>
+    <dependency>
+      <groupId>commons-cli</groupId>
+      <artifactId>commons-cli</artifactId>
+      <version>1.0</version>
+      <type>jar</type>
+    </dependency>
+    <dependency>
+      <groupId>commons-io</groupId>
+      <artifactId>commons-io</artifactId>
+      <version>1.1</version>
+      <type>jar</type>
+    </dependency>
+  </dependencies>
+  <reporting>
+    <plugins>
+      <plugin>
+        <artifactId>maven-javadoc-plugin</artifactId>
+        <configuration>
+          <source>1.5</source>
+          <links>
+            <link>http://java.sun.com/j2se/1.5.0/docs/api/</link>
+          </links>
+        </configuration>
+      </plugin>
+    </plugins>
+  </reporting>
+  <distributionManagement>
+    <repository>
+      <id>freehep-maven-deploy</id>
+      <name>FreeHEP Maven Repository</name>
+      <url>scpexe://svn.freehep.org/nfs/slac/g/jas/maven2</url>
+    </repository>
+  </distributionManagement>
+</project>
\ No newline at end of file

lcsim-cal-calib
project.properties added at 1.1
diff -N project.properties
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ project.properties	19 May 2008 21:31:58 -0000	1.1
@@ -0,0 +1,22 @@
+# Main class
+maven.jar.mainclass=org.lcsim.cal.calib.Calibrate
+
+# remove repository
+maven.repo.remote = http://mirrors.ibiblio.org/pub/mirrors/maven/,http://java.freehep.org/maven,http://www.lcsim.org/maven
+
+# JavaDoc links
+maven.javadoc.links=http://java.sun.com/j2se/1.5.0/docs/api/,http://java.freehep.org/lib/freehep/api/
+
+# skip tests?
+maven.test.skip = false
+
+maven.compile.source=1.5
+maven.compile.target=1.5
+maven.javadoc.source=1.5
+
+# increase max heap size for run script
+run.args=-Xmx1024m
+
+# increase max heap size for tests
+#maven.junit.jvmargs=-Xmx1024m
+maven.junit.jvmargs=-Xmx1536m

lcsim-cal-calib
project.xml added at 1.1
diff -N project.xml
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ project.xml	19 May 2008 21:31:58 -0000	1.1
@@ -0,0 +1,187 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project>
+  <pomVersion>3</pomVersion>
+  <artifactId>lcsim-cal-calib</artifactId>
+  <groupId>lcsim</groupId>
+  <currentVersion>0.1</currentVersion>
+  <organization>
+    <name>SLAC</name>
+    <url>http://www.slac.stanford.edu</url>
+  </organization>
+  <description>An org.lcsim based calorimeter calibration package.</description>
+  <shortDescription>Calorimeter Calibration Package</shortDescription>
+  <url>http://www.lcsim.org</url>
+  <issueTrackingUrl>http://jira.slac.stanford.edu/browse/SLIC</issueTrackingUrl>
+  <repository>
+    <connection>scm:cvs:pserver:[log in to unmask]:/cvs/lcd:SlicDiagnostics</connection>
+  </repository>
+  <name>lcsim-cal-calib</name>
+  <package>org.lcsim</package>
+  <inceptionYear>2008</inceptionYear>
+
+  <dependencies>
+    <dependency>
+      <groupId>lcsim</groupId>
+      <artifactId>lcsim</artifactId>
+      <version>1.1</version>
+      <url>http://www.lcsim.org</url>
+    </dependency>
+    <dependency>
+      <groupId>maven</groupId>
+      <artifactId>maven-scm-plugin</artifactId>
+      <version>1.5</version>
+      <url>http://www.ibiblio.org/maven</url>
+      <type>plugin</type>
+    </dependency>
+    <dependency>
+      <groupId>freehep</groupId>
+      <artifactId>freehep-jas-plugin</artifactId>
+      <version>1.0.2</version>
+      <url>http://java.freehep.org/freehep-jas-plugin</url>
+      <type>plugin</type>
+    </dependency>
+    <dependency>
+      <groupId>jdom</groupId>
+      <artifactId>jdom</artifactId>
+      <version>1.0</version>
+      <url>http://www.jdom.org/</url>
+    </dependency>
+    <dependency>
+      <groupId>jel</groupId>
+      <artifactId>jel</artifactId>
+      <version>0.9.10</version>
+      <url>http://galaxy.fzu.cz/JEL/</url>
+    </dependency>
+    <dependency>
+      <groupId>aida</groupId>
+      <artifactId>aida</artifactId>
+      <version>3.3</version>
+      <url>http://aida.freehep.org</url>
+    </dependency>
+    <dependency>
+      <groupId>aida</groupId>
+      <artifactId>aida-dev</artifactId>
+      <version>3.3</version>
+      <url>http://aida.freehep.org</url>
+    </dependency>
+    <dependency>
+      <groupId>net.java.dev</groupId>
+      <artifactId>truezip</artifactId>
+      <version>6.6</version>
+	<url>https://truezip.dev.java.net</url>
+    </dependency>
+    <dependency>
+      <groupId>freehep</groupId>
+      <artifactId>freehep-hep</artifactId>
+      <version>04112007</version>
+      <url>http://java.freehep.org</url>
+    </dependency>
+    <dependency>
+      <groupId>freehep</groupId>
+      <artifactId>freehep-base</artifactId>
+      <version>26092007</version>
+      <url>http://java.freehep.org</url>
+    </dependency>
+    <dependency>
+      <groupId>org.freehep</groupId>
+      <artifactId>freehep-physics</artifactId>
+      <version>2.1</version>
+      <url>http://java.freehep.org</url>
+      <properties>
+        <jas3.bundle>true</jas3.bundle>
+      </properties>
+    </dependency>
+    <dependency>
+      <groupId>org.freehep</groupId>
+      <artifactId>freehep-sio</artifactId>
+      <version>2.0</version>
+      <url>http://java.freehep.org</url>
+      <properties>
+        <jas3.bundle>true</jas3.bundle>
+      </properties>
+    </dependency>
+    <dependency>
+      <groupId>freehep</groupId>
+      <artifactId>freehep-heprep</artifactId>
+      <version>02032005</version>
+      <url>http://java.freehep.org</url>
+    </dependency>
+    <dependency>
+      <groupId>freehep</groupId>
+      <artifactId>freehep-jheprep</artifactId>
+      <version>02032005</version>
+      <url>http://java.freehep.org</url>
+    </dependency>
+    <dependency>
+      <groupId>freehep</groupId>
+      <artifactId>freehep-run-plugin</artifactId>
+      <version>1.1.1</version>
+      <url>http://java.freehep.org/maven/freehep/plugins</url>
+      <type>plugin</type>
+    </dependency>
+    <dependency>
+      <groupId>lcsim</groupId>
+      <artifactId>GeomConverter</artifactId>
+      <version>1.1</version>
+      <url>http://www.lcsim.org</url>
+      <properties>
+	<jas3.bundle>true</jas3.bundle>
+      </properties>
+    </dependency>
+    <dependency>
+      <groupId>netbeans</groupId>
+      <artifactId>openide-lookup</artifactId>
+      <version>1.9-patched-1.0</version>
+      <url>http://aidatld.freehep.org</url>
+    </dependency>
+    <dependency>
+      <groupId>jas</groupId>
+      <artifactId>jas3</artifactId>
+      <version>26092007</version>
+      <url>http://jas.freehep.org</url>
+    </dependency>
+    <dependency>
+      <groupId>commons-math</groupId>
+      <artifactId>commons-math</artifactId>
+      <version>1.0</version>
+      <type>jar</type>
+    </dependency>
+    <dependency>
+      <groupId>commons-cli</groupId>
+      <artifactId>commons-cli</artifactId>
+      <version>1.1</version>
+      <type>jar</type>
+    </dependency>
+    <dependency>
+      <groupId>commons-io</groupId>
+      <artifactId>commons-io</artifactId>
+      <version>1.1</version>
+      <type>jar</type>
+    </dependency>
+  </dependencies>
+  <build>
+    <sourceDirectory>src</sourceDirectory>
+    <sourceModifications>
+      <sourceModification>
+	<className>fakeClass</className>
+	<excludes>
+	  <exclude>org/lcsim/plugin/web/**</exclude>
+	</excludes>
+      </sourceModification>
+    </sourceModifications>
+    <unitTestSourceDirectory>test</unitTestSourceDirectory>
+    <resources>
+      <resource>
+	<directory>src</directory>
+	<includes>
+	  <include>org/lcsim/detector/**</include>
+	  <include>PLUGIN-inf/plugins.xml</include>
+	  <include>org/lcsim/**/*.menus</include>
+	  <!-- Temporary for DigiSim -->
+	  <include>org/lcsim/**/*.steer</include>
+	  <include>org/lcsim/plugin/web/**</include>
+	</includes>
+      </resource>
+    </resources>
+  </build>
+</project>
\ No newline at end of file

lcsim-cal-calib/src/org/lcsim/cal/calib
Calibrate.java added at 1.1
diff -N Calibrate.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ Calibrate.java	19 May 2008 21:31:58 -0000	1.1
@@ -0,0 +1,98 @@
+/*
+ * Calibrate.java
+ *
+ * Created on May 19, 2008, 11:50 AM
+ *
+ * $Id: Calibrate.java,v 1.1 2008/05/19 21:31:58 ngraf Exp $
+ */
+
+package org.lcsim.cal.calib;
+
+import java.io.BufferedReader;
+import java.io.File;
+import java.io.FileInputStream;
+import java.io.InputStreamReader;
+import java.util.ArrayList;
+import java.util.List;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.util.loop.LCIOEventSource;
+import org.lcsim.util.loop.LCSimLoop;
+
+/**
+ *
+ * @author Norman Graf
+ */
+public class Calibrate
+{
+    
+    /** Creates a new instance of Calibrate */
+    public Calibrate()
+    {
+    }
+    
+    public static void main(String[] args) throws Exception
+    {
+/*        
+        if(args.length<1)
+        {
+            usage();
+            return ;
+        }
+
+        String listOfFiles = args[0];
+ */
+        String listOfFiles = "C:/lcddata/gamma_Theta90_sid01_calibrate.txt";
+        List<File> filesToProcess = filesToProcess(listOfFiles);
+        LCIOEventSource src = new LCIOEventSource("SamplingFractionAnalysis", filesToProcess);
+        
+        int numToProcess=1;
+        if(args.length>1) numToProcess=Integer.parseInt(args[1]);
+        
+        System.out.println("Processing "+numToProcess+" events from "+listOfFiles);
+        LCSimLoop loop = new LCSimLoop();
+        loop.setLCIORecordSource(src);
+        System.out.println("adding the driver");
+        loop.add(new SamplingFractionAnalysisDriver());
+        System.out.println("looping");
+        try
+        {
+            loop.loop(numToProcess);
+        }
+        catch(Exception e)
+        {
+            e.printStackTrace();
+        }
+        System.out.println("done looping");
+        loop.dispose();
+        
+        int truncate = listOfFiles.lastIndexOf(".");
+        String listOfFilesName = listOfFiles.substring(0,truncate);
+        System.out.println(listOfFilesName);
+        AIDA.defaultInstance().saveAs("SamplingFractionAnalysis_"+listOfFilesName+".aida");
+    }
+    
+    public static void usage()
+    {
+        System.out.println("This is Calibrate");
+        System.out.println("usage:");
+        System.out.println("java Calibrate listOfInputFiles [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;
+    }
+    
+}

lcsim-cal-calib/src/org/lcsim/cal/calib
SamplingFractionAnalysisDriver.java added at 1.1
diff -N SamplingFractionAnalysisDriver.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ SamplingFractionAnalysisDriver.java	19 May 2008 21:31:58 -0000	1.1
@@ -0,0 +1,298 @@
+/*
+ * SamplingFractionAnalysisDriver.java
+ *
+ * Created on May 19, 2008, 11:54 AM
+ *
+ * $Id: SamplingFractionAnalysisDriver.java,v 1.1 2008/05/19 21:31:58 ngraf Exp $
+ */
+
+package org.lcsim.cal.calib;
+
+import Jama.Matrix;
+import hep.aida.ITree;
+import hep.physics.vec.Hep3Vector;
+import java.util.ArrayList;
+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.geometry.subdetector.CylindricalEndcapCalorimeter;
+import org.lcsim.recon.cluster.fixedcone.FixedConeClusterer;
+import org.lcsim.recon.cluster.fixedcone.FixedConeClusterer.FixedConeDistanceMetric;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+import static java.lang.Math.abs;
+import static java.lang.Math.sqrt;
+
+/**
+ *
+ * @author Norman Graf
+ */
+public class SamplingFractionAnalysisDriver extends Driver
+{
+    
+    // we will accumulate the raw energy values in three depths:
+    // 1. Layers (0)1 through (20)21 of the EM calorimeter (note that if layer 0 is massless, SF==1.)
+    // 2. last ten layers of the EM calorimeter
+    // 3. the hadron calorimeter
+    //
+    private double[][] _acc = new double[3][3];
+    private double[] _vec = new double[3];
+    
+    // let's use a clusterer to remove effects of calorimeter cells hit far, far away.
+    // use the only cross-detector clusterer we have:
+    private FixedConeClusterer _fcc;
+    
+    private AIDA aida = AIDA.defaultInstance();
+    private ITree _tree;
+    
+    private boolean _initialized;
+    
+    /** Creates a new instance of SamplingFractionAnalysisDriver */
+    public SamplingFractionAnalysisDriver()
+    {
+        _tree = aida.tree();
+    }
+    
+    
+    protected void process(EventHeader event)
+    {
+        // TODO make these values runtime definable
+        String[] det = {"EMBarrel","EMEndcap"};
+        String[] collNames = {"EcalBarrHits", "EcalEndcapHits", "HcalBarrHits", "HcalEndcapHits"};
+        double[] mipHistMaxBinVal = {.0005, .0005, .005, .005};
+        double timeCut = 100.; // cut on energy depositions coming more than 100 ns late
+        double ECalMipCut = .0001/3.; // determined from silicon using muons at normal incidence
+        double HCalMipCut = .0008/3.; // determined from scintillator using muons at normal incidence
+        double[] mipCut = {ECalMipCut, ECalMipCut, HCalMipCut, HCalMipCut};
+        
+        // TODO fix this dependence on EM calorimeter geometry
+        boolean skipFirstLayer = false;
+        int firstEmStartLayer = 0;
+        int secondEmStartLayer = 19;
+        
+        double emCalInnerRadius = 0.;
+        double emCalInnerZ = 0.;
+        if(!_initialized)
+        {
+            double radius = .5;
+            double seed = 0.;//.1;
+            double minE = .05; //.25;
+            _fcc = new FixedConeClusterer(radius, seed, minE, FixedConeDistanceMetric.DPHIDTHETA);
+            _initialized = true;
+            // detector geometries here...
+            // barrel
+            CylindricalCalorimeter calsubBarrel = (CylindricalCalorimeter)event.getDetector().getSubdetectors().get(det[0]);
+            // TODO remove this hardcoded dependence on the first layer
+            if(calsubBarrel.getLayering().getLayer(0).getSlices().get(0).isSensitive())
+            {
+                skipFirstLayer = true;
+                firstEmStartLayer += 1;
+                secondEmStartLayer += 1;
+            }
+//            Layering layering = calsubBarrel.getLayering();
+//            for(int i=0; i<layering.size(); ++i)
+//            {
+//                Layer l = layering.getLayer(i);
+//                System.out.println("layering "+i);
+//                List<LayerSlice> slices = l.getSlices();
+//                for(int j=0; j<slices.size(); ++j)
+//                {
+//                    LayerSlice slice = slices.get(j);
+//                    System.out.println("Layer "+i+" slice "+j+" is "+ slice.getMaterial().getName() +" and "+(slice.isSensitive() ? " is sensitive" : ""));
+//                }
+//            }
+            emCalInnerRadius = calsubBarrel.getInnerRadius();
+            //endcap
+            CylindricalEndcapCalorimeter calsubEndcap = (CylindricalEndcapCalorimeter)event.getDetector().getSubdetectors().get(det[1]);
+            emCalInnerZ = abs(calsubEndcap.getZMin());
+            if(skipFirstLayer) System.out.println("processing "+event.getDetectorName()+" with an em calorimeter with a massless first gap");
+            System.out.println("Calorimeter bounds: r= "+emCalInnerRadius+ " z= "+emCalInnerZ);
+            System.out.println("initialized...");
+        }
+        // organize the histogram tree by species and energy
+        List<MCParticle> mcparts = event.getMCParticles();
+        MCParticle mcpart = mcparts.get(mcparts.size()-1);
+        String particleType = mcpart.getType().getName();
+        double mcEnergy = mcpart.getEnergy();
+        long mcIntegerEnergy = Math.round(mcEnergy);
+        boolean meV = false;
+        if(mcEnergy<.99)
+        {
+            mcIntegerEnergy = Math.round(mcEnergy*1000);
+            meV = true;
+        }
+        
+        _tree.mkdirs(particleType);
+        _tree.cd(particleType);
+        _tree.mkdirs(mcIntegerEnergy+(meV ? "_MeV": "_GeV"));
+        _tree.cd(mcIntegerEnergy+(meV ? "_MeV": "_GeV"));
+        
+        // this analysis is intended for single particle calorimeter response.
+        // let's make sure that the primary particle did not interact in the
+        // tracker...
+        Hep3Vector endpoint = mcpart.getEndPoint();
+        // this is just crap. Why not use SpacePoint?
+        double radius = sqrt(endpoint.x()*endpoint.x()+endpoint.y()*endpoint.y());
+        double z = endpoint.z();
+//        System.out.println("Input MCParticle endpoint: r="+radius+" z= "+z);
+        
+        boolean doit = true;
+        if(radius<emCalInnerRadius && abs(z) < emCalInnerZ) doit = false;
+        if(doit)
+        {
+            // now let's check the em calorimeters...
+            // get all of the calorimeter hits...
+            List<CalorimeterHit> allHits = new ArrayList<CalorimeterHit>();
+            // and the list after cuts.
+            List<CalorimeterHit> hitsToCluster = new ArrayList<CalorimeterHit>();
+            int i = 0;
+            for(String name : collNames)
+            {
+//                System.out.println("fetching "+name+" from the event");
+                List<CalorimeterHit> hits = event.get(CalorimeterHit.class, name);
+//                System.out.println(name+ " has "+hits.size()+" hits");
+                // let's look at the hits and see if we need to cut on energy or time...
+                for(CalorimeterHit hit: hits)
+                {
+                    aida.histogram1D(name+" raw calorimeter cell energy",100, 0., mipHistMaxBinVal[i]).fill(hit.getRawEnergy());
+                    aida.histogram1D(name+" raw calorimeter cell energy full range",100, 0., 0.2).fill(hit.getRawEnergy());
+//                    aida.cloud1D(name+" raw calorimeter cell energies").fill(hit.getRawEnergy());
+                    aida.histogram1D(name+" calorimeter cell time",100,0., 200.).fill(hit.getTime());
+                    if(hit.getTime()<timeCut)
+                    {
+                        if(hit.getRawEnergy()>mipCut[i])
+                        {
+                            hitsToCluster.add(hit);
+                        }
+                    }
+                }
+                allHits.addAll(hits);
+                ++i;
+            }
+//            System.out.println("ready to cluster "+hitsToCluster.size()+ " hits");
+            
+            // cluster the hits
+            List<Cluster> clusters = _fcc.createClusters(hitsToCluster);
+//            System.out.println("found "+clusters.size()+" clusters");
+            aida.histogram1D("number of found clusters", 10, -0.5, 9.5).fill(clusters.size());
+            for(Cluster c : clusters)
+            {
+//                System.out.println(c);
+                aida.cloud1D("cluster energy").fill(c.getEnergy());
+            }
+            
+            // proceed only if we found a single cluster above threshold
+            if(clusters.size()==1)
+            {
+                Cluster c = clusters.get(0);
+                
+                aida.cloud1D("Single cluster events cluster energy").fill(c.getEnergy());
+                aida.cloud1D("Single cluster events cluster number of cells").fill(c.getCalorimeterHits().size());
+                
+                double clusterEnergy = c.getEnergy();
+                double mcMass = mcpart.getType().getMass();
+                // subtract the mass to get kinetic energy...
+                double expectedEnergy = sqrt(mcEnergy*mcEnergy - mcMass*mcMass);
+//                System.out.println(mcpart.getType().getName()+" "+expectedEnergy);
+                aida.cloud1D("measured - predicted energy").fill(clusterEnergy - expectedEnergy);
+                
+                // let's now break down the cluster by component.
+                // this analysis uses:
+                // 1.) first 20 EM layers
+                // 2.) next 10 EM layers
+                // 3.) Had layers
+                List<CalorimeterHit> hits = c.getCalorimeterHits();
+                double[] vals = new double[3];
+                double clusterRawEnergy = 0.;
+                for(CalorimeterHit hit : hits)
+                {
+                    long id = hit.getCellID();
+                    IDDecoder decoder = hit.getIDDecoder();
+                    decoder.setID(id);
+                    int layer = decoder.getLayer();
+                    String detectorName = decoder.getSubdetector().getName();
+                    int type = 0;
+                    if(detectorName.startsWith("EM"))
+                    {
+                        if(layer>=firstEmStartLayer && layer<secondEmStartLayer)
+                        {
+                            type = 0;
+                        }
+                        else
+                        {
+                            type = 1;
+                        }
+                    }
+                    if(detectorName.startsWith("HAD"))
+                    {
+                        type = 2;
+                    }
+                    clusterRawEnergy += hit.getRawEnergy();
+                    vals[type]+=hit.getRawEnergy();
+                } // end of loop over hits in cluster
+                // set up linear least squares:
+                // expectedEnergy = a*E1 + b*E2 +c*E3
+                for(int j=0; j<3; ++j)
+                {
+                    _vec[j] += expectedEnergy*vals[j];
+                    for(int k=0; k<3; ++k)
+                    {
+                        _acc[j][k] += vals[j]*vals[k];
+                    }
+                }
+            } // end of single cluster cut
+            
+            event.put("All Calorimeter Hits",allHits);
+            event.put("Hits To Cluster",hitsToCluster);
+            event.put("Found Clusters",clusters);
+        }// end of check on decays outside tracker volume
+        _tree.cd("/");
+    }
+    
+    protected void endOfData()
+    {
+        System.out.println("done! endOfData.");
+        // calculate the sampling fractions...
+        Matrix A = new Matrix(_acc, 3, 3);
+        A.print(6,4);
+        Matrix b = new Matrix(3,1);
+        for(int i=0; i<3; ++i)
+        {
+            b.set(i,0,_vec[i]);
+        }
+        b.print(6,4);
+        try
+        {
+            Matrix x = A.solve(b);
+            x.print(6, 4);
+        }
+        catch(Exception e)
+        {
+            e.printStackTrace();
+            System.out.println("try reducing dimensionality...");
+            Matrix Ap = new Matrix(_acc, 2, 2);
+            Ap.print(6,4);
+            Matrix bp = new Matrix(2,1);
+            for(int i=0; i<2; ++i)
+            {
+                bp.set(i,0,_vec[i]);
+            }
+            bp.print(6,4);
+            try
+            {
+                Matrix x = Ap.solve(bp);
+                x.print(6, 4);
+            }
+            catch(Exception e2)
+            {
+                e2.printStackTrace();
+            }
+        }
+        
+    }
+}
CVSspam 0.2.8