Print

Print


Commit in mcd-analysis/src/main/java/org/lcsim/mcd/analysis on MAIN
DRCorrection/alexDRCorr/proton_resolution_steering.xml+33added 1.1
                       /CorrectionFunctionSet.java+171added 1.1
                       /correctionFunctions.java+65added 1.1
                       /electronCorrection.java+326added 1.1
                       /ElectronCalibrationDriver.java+64added 1.1
                       /DRResolutionDriver.java+65added 1.1
                       /DualCorrection.java+353added 1.1
                       /Resolution.java+359added 1.1
                       /proton_calibration_steering.xml+33added 1.1
                       /cutConfiguration.java+83added 1.1
                       /electron_calibration_steering.xml+32added 1.1
                       /DRCalibrationDriver.java+65added 1.1
DRCorrection/DualCorrection.java+11-21.2 -> 1.3
            /ElectronCorrection.java+1-11.2 -> 1.3
ILCstdhepdriver.java+118-91.1 -> 1.2
+1779-12
12 added + 3 modified, total 15 files
Added new package, alexDRCorr, to implement certain changes to the analysis. C/S correction functions are now polynomials and stored as files and the 'DetectorConfiguration' is replaced by a simpler, more appropriate interface. I also use an S-only correction to compare against C/S.

mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
proton_resolution_steering.xml added at 1.1
diff -N proton_resolution_steering.xml
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ proton_resolution_steering.xml	18 Apr 2014 20:52:20 -0000	1.1
@@ -0,0 +1,33 @@
+<lcsim xmlns:xs="http://www.w3.org/2001/XMLSchema-instance"
+       xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/lcsim/1.0/lcsim.xsd">
+    
+    <classpath>
+        <jar>/home/aconway/NetBeansProjects/mcd-analysis/target/mcd-analysis-1.0-SNAPSHOT.jar</jar>
+    </classpath>
+
+    <inputFiles>
+        <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-1gev.slcio</file>
+        <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-2gev.slcio</file>
+        <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-5gev.slcio</file>
+        <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-10gev.slcio</file>
+        <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-25gev.slcio</file>
+        <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-50gev.slcio</file>
+    </inputFiles>
+    <control>
+        <numberOfEvents>-1</numberOfEvents>
+    </control>
+
+    <execute>
+        <driver name="DRResolutionDriver"/>
+    </execute>
+    <drivers>
+        <driver name="DRResolutionDriver" 
+                type="org.lcsim.mcd.analysis.DRCorrection.alexDRCorr.DRResolutionDriver">
+            <myP_name>proton</myP_name>
+            <myS_thresh>0.0</myS_thresh>
+            <myC_thresh>0.0</myC_thresh>
+            <myT_thresh>3.0</myT_thresh>
+            <myTag>default</myTag>
+        </driver>
+    </drivers>
+</lcsim>

mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
CorrectionFunctionSet.java added at 1.1
diff -N CorrectionFunctionSet.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ CorrectionFunctionSet.java	18 Apr 2014 20:52:20 -0000	1.1
@@ -0,0 +1,171 @@
+package org.lcsim.mcd.analysis.DRCorrection.alexDRCorr;
+
+import hep.aida.IFunction;
+import hep.aida.IFunctionFactory;
+import java.io.BufferedWriter;
+import java.io.File;
+import java.io.FileWriter;
+import java.io.IOException;
+import java.nio.charset.StandardCharsets;
+import java.nio.file.FileSystems;
+import java.nio.file.Files;
+import java.nio.file.Path;
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Scanner;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import org.lcsim.util.aida.AIDA;
+/**
+ * Read and write polynomial correction functions to a file,
+ * provide correction function set (HashMap<String,IFunction>) to drivers
+ * using a cutConfiguration object.
+ * 
+ * Lines in the files have the following format:
+ * ,<fn name>,<parm name>,<parm value>,<parm name>,...,\n
+ * @author aconway
+ */
+public class CorrectionFunctionSet {
+    
+    private AIDA aida;
+    private IFunctionFactory functionFactory;
+    private HashMap<String,IFunction> map;
+    private Path path;
+    
+    public CorrectionFunctionSet(cutConfiguration cutco) {
+        aida = AIDA.defaultInstance();
+        functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+        
+        String fname = cutco.toString() + ".txt";
+        path = FileSystems.getDefault().getPath("functions",fname);
+        String anyname = cutco.toStringAnyP() + ".txt";
+        Path anypath = FileSystems.getDefault().getPath("functions", anyname);
+        if (Files.isReadable(path)) {
+            map = loadFunctionsFromFile(path);
+        } else if (Files.isReadable(anypath)) {
+            map = loadFunctionsFromFile(anypath);
+            this.saveFunctionsToFile();
+        } else {
+            System.out.println("!!! Warning: Could not load functions from " + path.toString());
+            System.out.println("!!! Creating empty HashMap");
+            map = new HashMap<String,IFunction>();
+            this.saveFunctionsToFile();
+        }
+    }
+    
+    public CorrectionFunctionSet(cutConfiguration cutco, HashMap<String,IFunction> inmap) {
+        aida = AIDA.defaultInstance();
+        functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+        String fname = cutco.toString() + ".txt";
+        path = FileSystems.getDefault().getPath("functions",fname);
+        map = inmap;
+        this.saveFunctionsToFile();
+    }
+    
+    public CorrectionFunctionSet(cutConfiguration cutco, boolean DefaultEC) {
+        aida = AIDA.defaultInstance();
+        functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+        String fname = cutco.toString() + ".txt";
+        path = FileSystems.getDefault().getPath("functions",fname);
+        map = new HashMap<String,IFunction>();
+        if (DefaultEC) {
+            IFunction SC = functionFactory.createFunctionByName("SC", "p1");
+            IFunction CC = functionFactory.createFunctionByName("CC", "p1");
+            double[] defSCparms = {0.0,1.0};
+            double[] defCCparms = {3.6605e-3,1.87977e-5};
+            SC.setParameters(defSCparms);
+            CC.setParameters(defCCparms);
+            map.put("SC", SC);
+            map.put("CC", CC);
+        }
+    }
+    
+    public void addFunction(String name, IFunction func) {
+        this.map.put(name, func);
+        this.saveFunctionsToFile();
+    }
+    
+    public IFunction getFunction(String name) {
+        return this.map.get(name);
+    }
+    
+    private HashMap loadFunctionsFromFile(Path path) {
+        HashMap rv = new HashMap<String,IFunction>();
+        try {
+            List<String> lines = Files.readAllLines(path, StandardCharsets.UTF_8);
+            for (String line : lines) {
+                Scanner sc = new Scanner(line);
+                sc.useDelimiter(",");
+                String type = null;
+                String name = sc.next();
+                List<Double> parmlist = new ArrayList<Double>();
+                List<String> namelist = new ArrayList<String>();
+                while (sc.hasNext()){
+                    if (sc.hasNextDouble()) {
+                        parmlist.add(sc.nextDouble());
+                    } else {
+                        namelist.add(sc.next());
+                    }
+                }
+                if (name.matches("CS\\d.*")) {
+                    type = "p4";
+                } else if (name.matches("[CS]C")) {
+                    type = "p1";
+                } else if (name.matches("S\\d.*")) {
+                    type = "G";
+                } else {
+                    type = "p" + String.valueOf(parmlist.size()-1);
+                }
+                assert (parmlist.size() == namelist.size());
+                IFunction func = this.functionFactory.createFunctionByName(name, type);
+                for (int i=0; i<parmlist.size(); i++) {
+                    func.setParameter(namelist.get(i), parmlist.get(i));
+                }
+                rv.put(name, func);
+            }
+        } catch (IOException ex) {
+            System.out.println("!!! ERROR: Could not read file at " + path.toString());
+            Logger.getLogger(CorrectionFunctionSet.class.getName()).log(Level.SEVERE, null, ex);
+            System.exit(0);
+        }
+        return rv;
+    }
+    
+    private void saveFunctionsToFile() {
+        System.out.println("Saving functions to file " + this.path.toString());
+        try {
+            File ofile = this.path.toFile();
+            ofile.setWritable(true);
+            FileWriter fstream = new FileWriter(ofile, false);
+            BufferedWriter out = new BufferedWriter(fstream);
+            boolean first = true;
+            try {
+                for (String key : this.map.keySet()) {
+                    System.out.println("\nWriting function " + key);
+                    if (first) {
+                        first = false;
+                    } else {
+                        out.newLine();
+                    }
+                    IFunction func = this.map.get(key);
+                    String[] pnames = func.parameterNames();
+                    double[] params = func.parameters();
+                    out.write("," + key + ",");
+                    for (int i=0; i<params.length; i++) {
+                        System.out.print(pnames[i] + ", " + String.valueOf(params[i]) + ", ");
+                        out.write(pnames[i] + "," + String.valueOf(params[i]) + ",");
+                    }
+                }
+                out.close();
+                System.out.println("\nSuccessfully saved corrections to " +this.path.toString());
+            } catch (IOException ex) {
+                System.out.println("!!! ERROR: Error writing to file " + this.path.toString());
+                Logger.getLogger(CorrectionFunctionSet.class.getName()).log(Level.SEVERE, null, ex);
+            }
+        } catch (IOException ex) {
+            System.out.println("!!! ERROR: Could not open file for writing " + this.path.toString());
+            Logger.getLogger(CorrectionFunctionSet.class.getName()).log(Level.SEVERE, null, ex);
+        }
+    }
+}

mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
correctionFunctions.java added at 1.1
diff -N correctionFunctions.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ correctionFunctions.java	18 Apr 2014 20:52:20 -0000	1.1
@@ -0,0 +1,65 @@
+package org.lcsim.mcd.analysis.DRCorrection.alexDRCorr;
+
+import hep.aida.IFunction;
+import hep.aida.IFunctionFactory;
+import java.util.ArrayList;
+import java.util.Dictionary;
+import java.util.HashMap;
+import java.util.Map;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author aconway
+ */
+public class correctionFunctions {
+    
+    private AIDA aida;
+    private IFunctionFactory functionFactory;
+    private static Map<cutConfiguration, Map<String,IFunction>> map;
+    
+    public correctionFunctions() {
+        aida = AIDA.defaultInstance();
+        functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+        map = new HashMap<cutConfiguration, Map<String,IFunction>>();
+        Map submap = new HashMap<String,IFunction>();
+        
+        cutConfiguration cutcon = new cutConfiguration(0.0,0.0,0.0,"any","default");
+        double[] defSCparms = {0.0,1.0};
+        double[] defCCparms = {3.6605e-3,1.87977e-5};
+        IFunction SC = functionFactory.createFunctionByName("SC", "p1");
+        IFunction CC = functionFactory.createFunctionByName("CC", "p1");
+        SC.setParameters(defSCparms);
+        CC.setParameters(defCCparms);
+        submap.put("SC", SC);
+        submap.put("CC", CC);
+        double[] def1gparms = {0.0,0.0,0.0,0.0,0.0};
+        double[] def2gparms = {0.0,0.0,0.0,0.0,0.0};
+        double[] def5gparms = {0.0,0.0,0.0,0.0,0.0};
+        double[] def10gparms = {0.0,0.0,0.0,0.0,0.0};
+        double[] def25gparms = {0.0,0.0,0.0,0.0,0.0};
+        double[] def50gparms = {0.0,0.0,0.0,0.0,0.0};
+        double[] def100gparms = {0.0,0.0,0.0,0.0,0.0};
+        IFunction cs1gev = functionFactory.createFunctionByName("1gev", "p4");
+        IFunction cs2gev = functionFactory.createFunctionByName("2gev", "p4");
+        IFunction cs5gev = functionFactory.createFunctionByName("5gev", "p4");
+        IFunction cs10gev = functionFactory.createFunctionByName("10gev", "p4");
+        IFunction cs25gev = functionFactory.createFunctionByName("25gev", "p4");
+        IFunction cs50gev = functionFactory.createFunctionByName("50gev", "p4");
+        IFunction cs100gev = functionFactory.createFunctionByName("100gev", "p4");
+        cs1gev.setParameters(def1gparms);
+        cs2gev.setParameters(def2gparms);
+        cs5gev.setParameters(def5gparms);
+        cs10gev.setParameters(def10gparms);
+        cs25gev.setParameters(def25gparms);
+        cs50gev.setParameters(def50gparms);
+        cs100gev.setParameters(def100gparms);
+        submap.put("cs1gev", cs1gev);
+        submap.put("cs2gev", cs2gev);
+        submap.put("cs5gev", cs5gev);
+        submap.put("cs10gev", cs10gev);
+        submap.put("cs25gev", cs25gev);
+        submap.put("cs50gev", cs50gev);
+        submap.put("cs100gev", cs100gev);
+    }
+}

mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
electronCorrection.java added at 1.1
diff -N electronCorrection.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ electronCorrection.java	18 Apr 2014 20:52:20 -0000	1.1
@@ -0,0 +1,326 @@
+package org.lcsim.mcd.analysis.DRCorrection.alexDRCorr;
+
+import hep.aida.ICloud1D;
+import hep.aida.IDataPoint;
+import hep.aida.IDataPointSet;
+import hep.aida.IDataPointSetFactory;
+import hep.aida.IFitFactory;
+import hep.aida.IFitResult;
+import hep.aida.IFitter;
+import hep.aida.IFunction;
+import hep.aida.IFunctionFactory;
+import hep.aida.IHistogram1D;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import java.io.BufferedWriter;
+import java.io.FileWriter;
+import java.io.IOException;
+import java.util.List;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author aconway
+ */
+public class electronCorrection extends Driver {
+    String AIDAFile;
+    AIDA aida;
+    IFunctionFactory functionFactory;
+    IFitFactory fitFactory;
+    IDataPointSetFactory dpsFactory;
+    
+    double E_in;
+    double E_kin;
+    String Particlename;
+    Integer Ein;
+    Integer Ein_prev;
+    boolean firstEvent;
+    boolean first;
+    
+    ICloud1D Edep;
+    ICloud1D Eceren;
+    IFunction gauss;
+    IDataPointSet dps_emean = null;
+    IDataPointSet dps_cerenemean = null;
+    IDataPointSet dps_cerene = null;
+    IFitter jminuit;
+    IFitResult jminuitResult;
+    int point;
+    int point_cerenemean;
+    int point_cerene;
+    
+    cutConfiguration cutCon;
+    CorrectionFunctionSet corfus;
+    
+    public electronCorrection() {        
+        Ein_prev = 0;
+        firstEvent = true;
+        first = true;
+        
+        aida = AIDA.defaultInstance();
+        fitFactory = aida.analysisFactory().createFitFactory();
+        functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+        dpsFactory = aida.analysisFactory().createDataPointSetFactory(aida.tree());
+        gauss = functionFactory.createFunctionByName("gauss", "G");
+        dps_emean = dpsFactory.create("dps_emean", "electron response mean", 2);
+        dps_cerenemean = dpsFactory.create("dps_cerenemean", "electron response mean", 2);
+        dps_cerene = dpsFactory.create("dps_cerene_Chi2", "electron response mean", 2);
+        point = 0;
+        point_cerenemean = 0;
+    }
+    
+    @Override
+    protected void startOfData() {
+    }
+    
+    @Override
+    protected void process(EventHeader event) {
+        E_in = 0.0;
+        E_kin = 0.0;
+        double Chrg = 0.0;
+        double P_in = 0.0;
+        Hep3Vector origin = new BasicHep3Vector(0,0,0);
+        String DirName = null;
+        List<MCParticle> particles = event.get(MCParticle.class, event.MC_PARTICLES);
+        for (MCParticle particle : particles) {
+            if (particle.getParents().isEmpty()) {
+                E_in = E_in + particle.getEnergy();
+                E_kin = E_kin + particle.getEnergy() - particle.getMass();
+                P_in = particle.getMomentum().magnitude();
+                Particlename = particle.getType().getName();
+                origin = particle.getOrigin();
+                Chrg = particle.getCharge();
+            }
+            break;
+        }
+        if (E_kin < 0.999) {
+            Ein = (int) Math.floor(1000.0*E_kin + 0.5d);
+        } else {
+            Ein = (int) Math.floor(E_kin + 0.5d);
+        }
+        String E_str = Ein.toString();
+        DirName = Particlename.concat(E_str);
+        if (!Ein.equals(Ein_prev)) {
+            if (first) {
+                first = false;
+            } else {
+                System.out.println("electronCorrection:First Event");
+                convertandfit();
+            }
+            Ein_prev = Ein;
+            System.out.println("First Event:");
+            System.out.println("E_in:  " + E_in);
+            System.out.println("E_kin:  " + E_kin);
+            System.out.println("Name:  " + Particlename);
+            aida.tree().cd("/");
+            aida.tree().mkdir(DirName);
+            aida.tree().cd(DirName);
+            Edep = aida.histogramFactory().createCloud1D("Edep", "Energy Cloud", 100000, "autoConvert = false");
+            Eceren = aida.histogramFactory().createCloud1D("Eceren", "Cerenkov Cloud", 100000, "autoConvert = false");
+            System.out.println("DirName:  " + DirName);
+        }
+        double dif_curv = 0.0;
+        if (Chrg != 0.0) {
+            double rad_curv = P_in * 667.128;
+            double len_curv = rad_curv * 2.0 * Math.asin(0.667 / rad_curv);
+            dif_curv = len_curv - 1.334;
+        }
+        List<List<SimCalorimeterHit>> simCalorimeterHitCollections = event.get(SimCalorimeterHit.class);
+        double sumEEdep = 0.0;
+        double sumECeren = 0.0;
+        for (List<SimCalorimeterHit> simCalorimeterHits : simCalorimeterHitCollections) {
+            String CollectionName = event.getMetaData(simCalorimeterHits).getName();
+            if (CollectionName.contains("Edep")) {
+                for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+                    if (cutCon.T_thresh > 0.0 || cutCon.S_thresh > 0.0) {
+                        Hep3Vector diffvec = VecOp.sub(calorimeterHit.getPositionVec(), origin);
+                        int ncontribs = calorimeterHit.getMCParticleCount();
+                        for (int i=0; i<ncontribs; i++) {
+                            double tdep = calorimeterHit.getContributedTime(i);
+                            double tofl = tdep - (diffvec.magnitude()+dif_curv)/299.79;
+                            double edep = calorimeterHit.getContributedEnergy(i);
+                            if ( (tofl < cutCon.T_thresh) && (edep > cutCon.S_thresh/1000.0) ) {
+                                sumEEdep += edep;
+                            }
+                        }
+                    }
+                    else {
+                        sumEEdep = sumEEdep + calorimeterHit.getRawEnergy();
+                    }
+                }
+            } // end if Edep
+            if (CollectionName.contains("Opti")) {
+                for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+                    if (cutCon.T_thresh > 0.0 || cutCon.C_thresh > 0.0) {
+                        Hep3Vector diffvec = VecOp.sub(calorimeterHit.getPositionVec(), origin);
+                        int ncontribs = calorimeterHit.getMCParticleCount();
+                        for (int i=0; i<ncontribs; i++) {
+                            double tdep = calorimeterHit.getContributedTime(i);
+                            double tofl = tdep - (diffvec.magnitude()+dif_curv)/299.79;
+                            double edep = calorimeterHit.getContributedEnergy(i);
+                            if ( (tofl < cutCon.T_thresh) && (edep > cutCon.C_thresh) ) {
+                                sumECeren += edep;
+                            }
+                        }
+                    }
+                    else {
+                        sumECeren = sumECeren + calorimeterHit.getRawEnergy();
+                    }
+                }
+            }
+        }
+        Edep.fill(sumEEdep);
+        Eceren.fill(sumECeren);
+    }
+    
+    @Override
+    protected void endOfData() {
+        System.out.println("electronCorrection:End of Data:");
+        convertandfit();
+        aida.tree().cd("/");
+        IFunction line = functionFactory.createFunctionByName("line", "p1");
+        IFitter linefit = fitFactory.createFitter("Chi2", "jminuit");
+        IFitResult resultline = linefit.fit(dps_emean, line);
+        functionFactory.cloneFunction("e mean fitted line ", resultline.fittedFunction());
+        corfus.addFunction("SC",resultline.fittedFunction());
+        
+        double[] fresult = resultline.fittedParameters();
+        System.out.println("Ionization scale results:");
+        System.out.println("Chi2 = " + resultline.quality());
+        System.out.println(fresult[0] + " , " + fresult[1]);
+        //
+        // now deal with the cerenkov stuff:
+        //
+        String fitter = "Chi2";
+        System.out.println("Fitter:  " + fitter);
+        resultline = linefit.fit(dps_cerene, line);
+        String fname = "cerenkov  " + fitter + " fitted line";
+        System.out.println(fname);
+        functionFactory.cloneFunction(fname, resultline.fittedFunction());
+        corfus.addFunction("CC",resultline.fittedFunction());
+        
+        System.out.println(fitter + ":  " + resultline.quality());
+        fresult = resultline.fittedParameters();
+        System.out.println("Cerenkov scale results:");
+        System.out.println("Fitter:  " + fitter);
+        System.out.println("Chi2 = " + resultline.quality());
+        System.out.println(fresult[0] + " , " + fresult[1]);
+        try {
+            aida.saveAs(AIDAFile);
+        } catch (IOException ex) {
+            Logger.getLogger(electronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+        }
+    }
+    
+    protected void convertandfit() {
+        System.out.println("electronCorrection:convertandfit");
+        IHistogram1D Edep_conv;
+        if (Edep.isConverted()) {
+            Edep_conv = Edep.histogram();
+        } else {
+            System.out.println("Converting EDep");
+            double meanc = Edep.mean();
+            double rmsc = Edep.rms();
+            double nsigmas = 3.;
+            int nbins = 100;
+
+            double minx = meanc - nsigmas * rmsc;
+            double maxx = meanc + nsigmas * rmsc;
+            Edep.setConversionParameters(nbins, minx, maxx);
+            Edep.convertToHistogram();
+            Edep_conv = Edep.histogram();
+        }
+
+        dps_emean.addPoint();
+        IDataPoint dp = dps_emean.point(point);
+        dp.coordinate(1).setValue(Ein_prev);
+        dp.coordinate(1).setErrorPlus(Ein_prev * 0.001);
+        dp.coordinate(1).setErrorMinus(Ein_prev * 0.001);
+        dp.coordinate(0).setValue(Edep_conv.mean());
+        dp.coordinate(0).setErrorPlus(Edep_conv.rms());
+        dp.coordinate(0).setErrorMinus(Edep_conv.rms());
+        point++;
+
+        //
+        // now deal with cerenkov response:
+        //
+        IHistogram1D Eceren_conv;
+
+        if (Eceren.isConverted()) {
+            Eceren_conv = Eceren.histogram();
+        } else {
+            System.out.println("Converting Eceren");
+            double meanc = Eceren.mean();
+            double rmsc = Eceren.rms();
+            double nsigmas = 5.;
+            int nbins = 100;
+            double minx = meanc - nsigmas * rmsc;
+            double maxx = meanc + nsigmas * rmsc;
+            Eceren.setConversionParameters(nbins, minx, maxx);
+            Eceren.convertToHistogram();
+            Eceren_conv = Eceren.histogram();
+        }
+
+        dps_cerenemean.addPoint();
+        IDataPoint dp_cerenemean = dps_cerenemean.point(point_cerenemean);
+        dp_cerenemean.coordinate(1).setValue(Ein_prev);
+        dp_cerenemean.coordinate(1).setErrorPlus(Ein_prev * 0.001);
+        dp_cerenemean.coordinate(1).setErrorMinus(Ein_prev * 0.001);
+        dp_cerenemean.coordinate(0).setValue(Eceren_conv.mean());
+        dp_cerenemean.coordinate(0).setErrorPlus(Eceren_conv.rms());
+        dp_cerenemean.coordinate(0).setErrorMinus(Eceren_conv.rms());
+        point_cerenemean++;
+
+        String fitter = "Chi2";
+        System.out.println("Fitter:  " + fitter);
+        gauss.setParameter("amplitude", Eceren_conv.maxBinHeight());
+        gauss.setParameter("mean", Eceren_conv.mean());
+        gauss.setParameter("sigma", Eceren_conv.rms());
+        jminuit = fitFactory.createFitter(fitter, "jminuit");
+        jminuitResult = jminuit.fit(Eceren_conv, gauss);
+        System.out.println("jminuit " + fitter + ":  " + jminuitResult.quality());
+        double[] result = jminuitResult.fittedParameters();
+        double[] errors = jminuitResult.errors();
+        String functionname = "ceren fitted gauss  " + fitter;
+        System.out.println(functionname);
+        functionFactory.cloneFunction(functionname, jminuitResult.fittedFunction());
+        System.out.println(result[0] + "," + result[1] + "," + result[2]);
+        IDataPoint dp_cerene = dps_cerene.addPoint();
+        dp_cerene.coordinate(0).setValue(result[1]);
+        dp_cerene.coordinate(0).setErrorPlus(Math.abs(result[2]));
+        dp_cerene.coordinate(0).setErrorMinus(Math.abs(result[2]));
+        dp_cerene.coordinate(1).setValue(Ein_prev);
+        dp_cerene.coordinate(1).setErrorPlus(Ein_prev * 0.001);
+        dp_cerene.coordinate(1).setErrorMinus(Ein_prev * 0.001);
+        point_cerene++;
+
+        try {
+            aida.saveAs(AIDAFile);
+        } catch (IOException ex) {
+            Logger.getLogger(electronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+        }
+    }
+    
+    public void setCorrFu(CorrectionFunctionSet corrs) {
+        System.out.println("electronCorrection:setCorrFu");
+        this.corfus = corrs;
+    }
+    
+    public void setCutConfiguration(cutConfiguration cutco) {
+        System.out.println("electronCorrection:setCutConfiguration");
+        this.cutCon = cutco;
+    }
+    
+    public void setMyAIDAFilename(String AIDAFile) {
+        System.out.println("electronCorrection:setMyAIDAFilename");
+        this.AIDAFile = AIDAFile;
+        System.out.println(AIDAFile);
+    }
+}

mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
ElectronCalibrationDriver.java added at 1.1
diff -N ElectronCalibrationDriver.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ ElectronCalibrationDriver.java	18 Apr 2014 20:52:20 -0000	1.1
@@ -0,0 +1,64 @@
+package org.lcsim.mcd.analysis.DRCorrection.alexDRCorr;
+
+import org.lcsim.util.Driver;
+
+/**
+ *
+ * @author aconway
+ */
+public class ElectronCalibrationDriver extends Driver {
+    
+    private String P_name = "any";
+    private String tag = "default";
+    private Double C_thresh = 0.0;
+    private Double S_thresh = 0.0;
+    private Double T_thresh = 0.0;
+    private String AIDAFile = null;
+    
+    electronCorrection elec;
+    
+    public ElectronCalibrationDriver() {
+        elec = new electronCorrection();
+        add(elec);
+    }
+    
+    @Override
+    public void startOfData() {
+        System.out.println("ElectronCalibrationDriver:startOfData");
+        cutConfiguration cutco;
+        if (tag.equals("default")) {
+            cutco = new cutConfiguration(C_thresh, S_thresh, T_thresh, "any");
+        } else {
+            cutco = new cutConfiguration(C_thresh, S_thresh, T_thresh, "any", tag);
+        }
+        AIDAFile = "EC-" + cutco.toString() + ".aida";
+        elec.setMyAIDAFilename(AIDAFile);
+        elec.setCutConfiguration(cutco);
+        elec.setCorrFu(new CorrectionFunctionSet(cutco, false));
+    }
+    
+    public void setMyP_name(String pname) {
+        System.out.println("ElectronCalibrationDriver:setMyP_name = " + pname);
+        this.P_name = pname;
+    }
+    
+    public void setMyTag(String tag) {
+        System.out.println("ElectronCalibrationDriver:setMyTag = " + tag);
+        this.tag = tag;
+    }
+    
+    public void setMyC_thresh(double cthresh) {
+        System.out.println("ElectronCalibrationDriver:setMyC_thresh = " + cthresh);
+        this.C_thresh = cthresh;
+    }
+    
+    public void setMyS_thresh(double sthresh) {
+        System.out.println("ElectronCalibrationDriver:setMyS_thresh = " + sthresh);
+        this.S_thresh = sthresh;
+    }
+    
+    public void setMyT_thresh(double tthresh) {
+        System.out.println("ElectronCalibrationDriver:setMyT_thresh = " + tthresh);
+        this.T_thresh = tthresh;
+    }
+}

mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
DRResolutionDriver.java added at 1.1
diff -N DRResolutionDriver.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ DRResolutionDriver.java	18 Apr 2014 20:52:20 -0000	1.1
@@ -0,0 +1,65 @@
+package org.lcsim.mcd.analysis.DRCorrection.alexDRCorr;
+
+import org.lcsim.util.Driver;
+
+/**
+ *
+ * @author aconway
+ */
+public class DRResolutionDriver extends Driver {
+    
+    private String P_name = "any";
+    private String tag = "default";
+    private Double C_thresh = 0.0;
+    private Double S_thresh = 0.0;
+    private Double T_thresh = 0.0;
+    private String AIDAFile = null;
+    
+    Resolution resol;
+    
+    public DRResolutionDriver() {
+        resol = new Resolution();
+        add(resol);
+    }
+    
+    @Override
+    public void startOfData() {
+        System.out.println("DRResolutionDriver:startOfData");
+        cutConfiguration cutco;
+        if (tag.equals("default")) {
+            cutco = new cutConfiguration(C_thresh, S_thresh, T_thresh, P_name);
+        } else {
+            cutco = new cutConfiguration(C_thresh, S_thresh, T_thresh, P_name, tag);
+        }
+        AIDAFile = "res-" + cutco.toString() + ".aida";
+        resol.setMyAIDAFilename(AIDAFile);
+        resol.setCutConfiguration(cutco);
+        resol.setCorrFu(new CorrectionFunctionSet(cutco));
+        resol.startOfData();
+    }
+    
+    public void setMyP_name(String pname) {
+        System.out.println("DRResolutionDriver:setMyP_name = " + pname);
+        this.P_name = pname;
+    }
+    
+    public void setMyTag(String tag) {
+        System.out.println("DRResolutionDriver:setMyTag = " + tag);
+        this.tag = tag;
+    }
+    
+    public void setMyC_thresh(double cthresh) {
+        System.out.println("DRResolutionDriver:setMyC_thresh = " + cthresh);
+        this.C_thresh = cthresh;
+    }
+    
+    public void setMyS_thresh(double sthresh) {
+        System.out.println("DRResolutionDriver:setMyS_thresh = " + sthresh);
+        this.S_thresh = sthresh;
+    }
+    
+    public void setMyT_thresh(double tthresh) {
+        System.out.println("DRResolutionDriver:setMyT_thresh = " + tthresh);
+        this.T_thresh = tthresh;
+    }
+}

mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
DualCorrection.java added at 1.1
diff -N DualCorrection.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ DualCorrection.java	18 Apr 2014 20:52:20 -0000	1.1
@@ -0,0 +1,353 @@
+package org.lcsim.mcd.analysis.DRCorrection.alexDRCorr;
+
+import hep.aida.ICloud1D;
+import hep.aida.ICloud2D;
+import hep.aida.IDataPoint;
+import hep.aida.IDataPointSet;
+import hep.aida.IDataPointSetFactory;
+import hep.aida.IFitFactory;
+import hep.aida.IFitResult;
+import hep.aida.IFitter;
+import hep.aida.IFunction;
+import hep.aida.IFunctionFactory;
+import hep.aida.IHistogram1D;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import java.io.IOException;
+import java.util.List;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author aconway
+ */
+public class DualCorrection extends Driver {
+    
+    String AIDAFile;
+    AIDA aida = AIDA.defaultInstance();
+    IFunctionFactory functionFactory;
+    IFitFactory fitFactory;
+    IDataPointSetFactory dpsFactory;
+    
+    ICloud1D Edep;
+    ICloud1D Eceren;
+    ICloud1D Edep_frac;
+    ICloud1D Edep_cor;
+    ICloud1D Eceren_cor;
+    ICloud2D c_efrac_ratio;
+    ICloud2D c_Ceren_vs_Edep;
+    ICloud1D[] slice;
+    IHistogram1D[] conv_slice;
+    IDataPointSet dps_ratio;
+    IFunction gauss;
+    IFunction poly;
+    IFitter jminuit;
+    IFitResult jminuitResult;
+    double binwidth = 0.05;
+    int maxbin = (int) Math.rint(1.5 / binwidth);
+    
+    double E_in;
+    double E_kin;
+    String Particlename;
+    Integer Ein;
+    Integer Ein_prev;
+    boolean firstEvent;
+    boolean first;
+    
+    cutConfiguration cutCon;
+    CorrectionFunctionSet corfus;
+    
+    public DualCorrection() {
+        Ein_prev = 0;
+        firstEvent = true;
+        first = true;
+        
+        fitFactory = aida.analysisFactory().createFitFactory();
+        functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+        dpsFactory = aida.analysisFactory().createDataPointSetFactory(aida.tree());
+        
+        gauss = functionFactory.createFunctionByName("gauss", "G");
+        jminuit = fitFactory.createFitter("Chi2", "jminuit");
+        poly = functionFactory.createFunctionByName("poly", "p4");
+        poly.setParameter("p0", 0.63);
+        poly.setParameter("p1", 0.54);
+        poly.setParameter("p2", -0.58);
+        poly.setParameter("p3", 0.5);
+        poly.setParameter("p4", 0.0);
+    }
+    
+    // initialize variables, etc.
+    @Override
+    protected void startOfData() {
+        
+    }
+    
+    // Sum up Edep and Eceren, run fits
+    @Override
+    protected void process(EventHeader event) {
+        E_in = 0.0;
+        E_kin = 0.0;
+        double Chrg = 0.0;
+        double P_in = 0.0;
+        Hep3Vector origin = new BasicHep3Vector(0,0,0);
+        String DirName = null;
+        List<MCParticle> particles = event.get(MCParticle.class, event.MC_PARTICLES);
+        for (MCParticle particle : particles) {
+            if (particle.getParents().isEmpty()) {
+                E_in = E_in + particle.getEnergy();
+                E_kin = E_kin + particle.getEnergy() - particle.getMass();
+                P_in = particle.getMomentum().magnitude();
+                Particlename = particle.getType().getName();
+                origin = particle.getOrigin();
+                Chrg = particle.getCharge();
+            }
+            break;
+        }
+        if (E_kin < 0.999) {
+            Ein = (int) Math.floor(1000.0*E_kin + 0.5d);
+        } else {
+            Ein = (int) Math.floor(E_kin + 0.5d);
+        }
+        String E_str = Ein.toString();
+        DirName = Particlename.concat(E_str);
+        if (!Ein.equals(Ein_prev)) {
+            if (first) {
+                first = false;
+            } else {
+                System.out.println("Process event: ");
+                convertandfit(slice, conv_slice);
+            }
+            Ein_prev = Ein;
+            System.out.println("First Event:");
+            System.out.println("E_in:  " + E_in);
+            System.out.println("E_kin:  " + E_kin);
+            System.out.println("Name:  " + Particlename);
+            aida.tree().cd("/");
+            aida.tree().mkdir(DirName);
+            aida.tree().cd(DirName);
+            Edep = aida.histogramFactory().createCloud1D("Edep", "Energy Cloud", 100000, "autoConvert = false");
+            Edep_frac = aida.histogramFactory().createCloud1D("Edep frac", "Energy frac Cloud", 100000, "autoConvert = false");
+            Eceren = aida.histogramFactory().createCloud1D("Eceren", "Cerenkov Cloud", 100000, "autoConvert = false");
+            Edep_cor = aida.histogramFactory().createCloud1D("Edep_cor", "corrected Energy Cloud", 100000, "autoConvert = false");
+            Eceren_cor = aida.histogramFactory().createCloud1D("Eceren_cor", "corrected Cerenkov Cloud", 100000, "autoConvert = false");
+            c_Ceren_vs_Edep = aida.histogramFactory().createCloud2D("c_Ceren_vs_Edep", "ceren vs Edep", 100000, "autoConvert = false");
+            slice = new ICloud1D[maxbin + 1];
+            conv_slice = new IHistogram1D[maxbin + 1];
+            for (int ibin = 0; ibin < maxbin + 1; ibin++) {
+                String cloudy = "ratio_" + ibin;
+                slice[ibin] = aida.histogramFactory().createCloud1D(cloudy, cloudy, 500000, "autoConvert = false");
+            }
+            System.out.println("DirName:  " + DirName);
+        }
+        double dif_curv = 0.0;
+        if (Chrg != 0.0) {
+            double rad_curv = P_in * 667.128;
+            double len_curv = rad_curv * 2.0 * Math.asin(0.667 / rad_curv);
+            dif_curv = len_curv - 1.334;
+        }
+        List<List<SimCalorimeterHit>> simCalorimeterHitCollections = event.get(SimCalorimeterHit.class);
+        double sumEEdep = 0.0;
+        double sumECeren = 0.0;
+        for (List<SimCalorimeterHit> simCalorimeterHits : simCalorimeterHitCollections) {
+            String CollectionName = event.getMetaData(simCalorimeterHits).getName();
+            if (CollectionName.contains("Edep")) {
+                for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+                    if (cutCon.T_thresh > 0.0 || cutCon.S_thresh > 0.0) {
+                        Hep3Vector diffvec = VecOp.sub(calorimeterHit.getPositionVec(), origin);
+                        int ncontribs = calorimeterHit.getMCParticleCount();
+                        for (int i=0; i<ncontribs; i++) {
+                            double tdep = calorimeterHit.getContributedTime(i);
+                            double tofl = tdep - (diffvec.magnitude()+dif_curv)/299.79;
+                            double edep = calorimeterHit.getContributedEnergy(i);
+                            if ( (tofl < cutCon.T_thresh) && (edep > cutCon.S_thresh/1000.0) ) {
+                                sumEEdep += edep;
+                            }
+                        }
+                    }
+                    else {
+                        sumEEdep = sumEEdep + calorimeterHit.getRawEnergy();
+                    }
+                }
+            } // end if Edep
+            if (CollectionName.contains("Opti")) {
+                for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+                    if (cutCon.T_thresh > 0.0 || cutCon.C_thresh > 0.0) {
+                        Hep3Vector diffvec = VecOp.sub(calorimeterHit.getPositionVec(), origin);
+                        int ncontribs = calorimeterHit.getMCParticleCount();
+                        for (int i=0; i<ncontribs; i++) {
+                            double tdep = calorimeterHit.getContributedTime(i);
+                            double tofl = tdep - (diffvec.magnitude()+dif_curv)/299.79;
+                            double edep = calorimeterHit.getContributedEnergy(i);
+                            if ( (tofl < cutCon.T_thresh) && (edep > cutCon.C_thresh) ) {
+                                sumECeren += edep;
+                            }
+                        }
+                    }
+                    else {
+                        sumECeren = sumECeren + calorimeterHit.getRawEnergy();
+                    }
+                }
+            }
+        }
+        Edep.fill(sumEEdep);
+        Eceren.fill(sumECeren);
+        Edep_frac.fill(sumEEdep/E_kin);
+        
+        double[] xval = {0.0};
+        xval[0] = sumEEdep;
+        double sumEEdep_cor = corfus.getFunction("SC").value(xval);
+        xval[0] = sumECeren;
+        double sumECeren_cor = corfus.getFunction("CC").value(xval);
+        Edep_cor.fill(sumEEdep_cor);
+        Eceren_cor.fill(sumECeren_cor);
+        double ratio = sumECeren_cor / sumEEdep_cor;
+        double fraction = sumEEdep_cor / E_kin;
+        aida.cloud1D("frac").fill(fraction);
+        aida.cloud1D("c_CerenEdep_ratio").fill(ratio);
+        c_Ceren_vs_Edep.fill(sumECeren_cor, sumEEdep_cor);
+        aida.cloud2D("c_efrac_ratio").fill(ratio, fraction);
+        
+        int bin = (int) Math.rint(ratio / binwidth);
+        if (bin < 0) {
+            bin = 0;
+        }
+        if (bin > maxbin - 1) {
+            bin = maxbin;
+        }
+        slice[bin].fill(fraction);
+    }
+    
+    // Run final fits
+    @Override
+    protected void endOfData() {
+        System.out.println("DualCorrection: endOfData");
+        convertandfit(slice, conv_slice);
+    }
+    
+    protected void convertandfit(ICloud1D[] slices, IHistogram1D[] conv_slices) {
+        System.out.println("DualCorrection: convert and fit:");
+
+        String fitter = "Chi2";
+        String dpsname = "dps_ratio_" + fitter;
+        dps_ratio = dpsFactory.create(dpsname, "ratio", 2);
+
+        for (int i = 0; i < slices.length; i++) {
+            if (slices[i].isConverted()) {
+                System.out.println("convert and fit already converted");
+                conv_slices[i] = slices[i].histogram();
+            } else {
+                double meanc = slices[i].mean();
+                double rmsc = slices[i].rms();
+                double nsigmas = 3.;
+                int nbins = 100;
+                double minx = meanc - nsigmas * rmsc;
+                double maxx = meanc + nsigmas * rmsc;
+                slices[i].setConversionParameters(nbins, minx, maxx);
+                slices[i].convertToHistogram();
+                conv_slices[i] = slices[i].histogram();
+            }
+            int entries = conv_slices[i].entries();
+            if (entries > 300) {
+                System.out.println("Fitter:  " + fitter);
+                double meanguess = 0;
+                for (int j = 0; j < conv_slices[i].axis().bins(); j++) {
+                    if (conv_slices[i].binHeight(j) == conv_slices[i].maxBinHeight()) {
+                        meanguess = conv_slices[i].axis().binCenter(j);
+                    }
+                }
+                gauss.setParameter("amplitude", conv_slices[i].maxBinHeight());
+                gauss.setParameter("mean", meanguess);
+                gauss.setParameter("sigma", conv_slices[i].rms() / 2.0);
+                jminuit = fitFactory.createFitter(fitter, "jminuit");
+                try {
+                    jminuitResult = jminuit.fit(conv_slices[i], gauss);
+                } catch (RuntimeException ex) {
+                    System.out.println("Runtime Exception caught in fit!");
+                    System.out.println(ex.toString());
+                } catch (Exception ex) {
+                    System.out.println("Exception caught in fit!");
+                    System.out.println(ex.toString());
+                }
+                System.out.println("jminuit " + fitter + ":  " + jminuitResult.quality());
+                double[] result = jminuitResult.fittedParameters();
+                String functionname = "ratio fitted gauss  slices: " + i + "  " + fitter;
+                System.out.println(functionname);
+                functionFactory.cloneFunction(functionname, jminuitResult.fittedFunction());
+                System.out.println(result[0] + "," + result[1] + "," + result[2]);
+                IDataPoint dp_cerene = dps_ratio.addPoint();
+                double ratio = (double) (i * binwidth) + 0.5 * binwidth;
+                dp_cerene.coordinate(0).setValue(ratio);
+                dp_cerene.coordinate(0).setErrorPlus(binwidth * 0.5);
+                dp_cerene.coordinate(0).setErrorMinus(binwidth * 0.5);
+                dp_cerene.coordinate(1).setValue(result[1]);
+                dp_cerene.coordinate(1).setErrorPlus(Math.abs(result[2]));
+                dp_cerene.coordinate(1).setErrorMinus(Math.abs(result[2]));
+            }
+        }
+        //Now fit the data point sets: 
+        System.out.println("Fitter:  " + fitter);
+        jminuitResult = jminuit.fit(dps_ratio, poly);
+        String fname = "correction  " + fitter + " fitted poly";
+        System.out.println(fname);
+        functionFactory.cloneFunction(fname, jminuitResult.fittedFunction());
+        System.out.println(fitter + ":  " + jminuitResult.quality());
+        corfus.addFunction("CS" + Ein_prev.toString() + "gev", jminuitResult.fittedFunction());
+        
+        // Now get Edep-only correction
+        double meanc = Edep_frac.mean();
+        double rmsc = Edep_frac.rms();
+        double nsigmas = 3.;
+        int nbins = 100;
+        double minx = meanc - nsigmas * rmsc;
+        double maxx = meanc + nsigmas * rmsc;
+        Edep_frac.setConversionParameters(nbins, minx, maxx);
+        Edep_frac.convertToHistogram();
+        IHistogram1D Edep_conv = Edep_frac.histogram();
+        gauss.setParameter("amplitude", Edep_conv.maxBinHeight());
+        gauss.setParameter("mean", Edep_conv.mean());
+        gauss.setParameter("sigma", Edep_conv.rms());
+        jminuit = fitFactory.createFitter(fitter, "jminuit");
+        try {
+            jminuitResult = jminuit.fit(Edep_conv, gauss);
+        } catch (RuntimeException ex) {
+            System.out.println("Runtime Exception caught in fit!");
+            System.out.println(ex.toString());
+        } catch (Exception ex) {
+            System.out.println("Exception caught in fit!");
+            System.out.println(ex.toString());
+        }
+        fname = "S" + Ein_prev.toString() + "gev";
+        functionFactory.cloneFunction("fname", jminuitResult.fittedFunction());
+        corfus.addFunction(fname, jminuitResult.fittedFunction());
+
+        try {
+            aida.saveAs(AIDAFile);
+            System.out.println("AIDAFile Saved.");
+        } catch (IOException ex) {
+            Logger.getLogger(org.lcsim.mcd.analysis.DRCorrection.DualCorrection.class.getName()).log(Level.SEVERE, null, ex);
+        }
+    }
+    
+    public void setCorrFu(CorrectionFunctionSet corrs) {
+        System.out.println("DualCorrection:setCorrFu");
+        this.corfus = corrs;
+    }
+    
+    public void setCutConfiguration(cutConfiguration cutco) {
+        System.out.println("DualCorrection:setCutConfiguration");
+        this.cutCon = cutco;
+    }
+    
+    public void setMyAIDAFilename(String AIDAFile) {
+        System.out.println("DualCorrection:setMyAIDAFilename");
+        this.AIDAFile = AIDAFile;
+        System.out.println(AIDAFile);
+    }
+}

mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
Resolution.java added at 1.1
diff -N Resolution.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ Resolution.java	18 Apr 2014 20:52:20 -0000	1.1
@@ -0,0 +1,359 @@
+package org.lcsim.mcd.analysis.DRCorrection.alexDRCorr;
+
+/**
+ *
+ * @author wenzel
+ * @date
+ */
+import hep.aida.ICloud1D;
+import hep.aida.IDataPoint;
+import hep.aida.IDataPointSet;
+import hep.aida.IDataPointSetFactory;
+import hep.aida.IFitFactory;
+import hep.aida.IFitResult;
+import hep.aida.IFitter;
+import hep.aida.IFunction;
+import hep.aida.IFunctionFactory;
+import hep.aida.IHistogram1D;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import java.io.IOException;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import org.lcsim.util.aida.AIDA;
+import java.util.List;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.util.Driver;
+
+public class Resolution extends Driver {
+
+    private AIDA aida = AIDA.defaultInstance();
+    String AIDAFile;
+    String file_name;
+    IFunctionFactory functionFactory;
+    IFitFactory fitFactory;
+    IDataPointSetFactory dpsFactory;
+    IFunction gauss;
+    static boolean first;
+    static boolean firstEvent;
+    static double E_in;
+    Integer Ein;
+    Integer Ein_prev;
+    double E_kin;
+    String Particlename;
+    ICloud1D Edep;
+    ICloud1D ECeren;
+    ICloud1D Edep_cor;
+    ICloud1D Edep_dcor;
+    ICloud1D Edep_scor;
+    ICloud1D Eceren_cor;
+    double nsigmas=5.0;
+    int nbins=100;
+    IDataPointSet[] dps_Correctede = null;
+    IDataPointSet[] dps_DRResolution = null;
+    IDataPointSet[] dps_SRResolution = null;
+    IDataPointSet dps_ECorrFrac = null;
+    IDataPointSet dps_CCorrFrac = null;
+    IDataPointSet[] dps_DCorrFrac = null;
+    int point_Correctede[];
+    double[] result;
+    double errors[];
+//    String[] Fitters = {"Chi2", "leastsquares"};
+    String[] Fitters = {"Chi2"};
+    // available are "Chi2", "leastsquares", "bml", "cleverchi2"
+    // but "bml", "cleverchi2" seem to produce rubish
+
+    IFitter jminuit;
+    IFitResult jminuitResult;
+    IFunction E_cor;
+    IFunction C_cor;
+    IFunction D_cor;
+    double xval[] = {10.};
+    
+    
+    cutConfiguration cutCon;
+    CorrectionFunctionSet corfus;
+
+    @Override
+    protected void startOfData() {
+        System.out.println("Start of Data:");
+        dps_Correctede = new IDataPointSet[Fitters.length];
+        dps_DRResolution = new IDataPointSet[Fitters.length];
+        dps_SRResolution = new IDataPointSet[Fitters.length];
+        dps_DCorrFrac = new IDataPointSet[Fitters.length];
+        point_Correctede = new int[Fitters.length];
+        fitFactory = aida.analysisFactory().createFitFactory();
+        functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+        dpsFactory = aida.analysisFactory().createDataPointSetFactory(aida.tree());
+        dps_ECorrFrac = dpsFactory.create("dps_ECorrFrac_", 2);
+        dps_CCorrFrac = dpsFactory.create("dps_CCorrFrac_", 2);
+        for (int i = 0; i < Fitters.length; i++) {
+            String dpsname = "dps_Correctede_" + Fitters[i];
+            point_Correctede[i] = 0;
+            dps_Correctede[i] = dpsFactory.create(dpsname, "electron response mean", 2);
+            dps_DRResolution[i] = dpsFactory.create("dps_DRResolution_"+Fitters[i], 2);
+            dps_SRResolution[i] = dpsFactory.create("dps_SRResolution_"+Fitters[i], 2);
+            dps_DCorrFrac[i] = dpsFactory.create("dps_DCorrFrac_"+Fitters[i], 2);
+        }
+
+        Ein_prev = 0;
+        firstEvent = true;
+        first = true;
+        gauss = functionFactory.createFunctionByName("gauss", "G");
+    }
+
+    @Override
+    protected void process(EventHeader event) {
+        E_in = 0.0;
+        E_kin = 0.0;
+        double Chrg = 0.0;
+        double P_in = 0.0;
+        Hep3Vector origin = new BasicHep3Vector(0,0,0);
+        String DirName = null;
+        List<MCParticle> particles = event.get(MCParticle.class, event.MC_PARTICLES);
+        for (MCParticle particle : particles) {
+            if (particle.getParents().isEmpty()) {
+                E_in = E_in + particle.getEnergy();
+                E_kin = E_kin + particle.getEnergy() - particle.getMass();
+                P_in = particle.getMomentum().magnitude();
+                Particlename = particle.getType().getName();
+                origin = particle.getOrigin();
+                Chrg = particle.getCharge();
+            }
+            break;
+        }
+        if (E_kin < 0.999) {
+            Ein = (int) Math.floor(1000.0*E_kin + 0.5d);
+        } else {
+            Ein = (int) Math.floor(E_kin + 0.5d);
+        }
+        String E_str = Ein.toString();
+        DirName = Particlename.concat(E_str);
+        if (!Ein.equals(Ein_prev)) {
+            if (first) {
+                first = false;
+            } else {
+                System.out.println("E_in:  " + E_in);
+                System.out.println("Ein_prev:  " + Ein_prev);
+                convertandfit();
+            }
+            Ein_prev = Ein;
+            System.out.println("First Event:");
+            System.out.println("E_in:  " + E_in);
+            System.out.println("E_kin:  " + E_kin);
+            System.out.println("Name:  " + Particlename);
+            aida.tree().cd("/");
+            aida.tree().mkdir(DirName);
+            aida.tree().cd(DirName);
+            Edep = aida.histogramFactory().createCloud1D("Edep", "uncorrected Energy Cloud", 100000, "autoConvert = false");
+            ECeren = aida.histogramFactory().createCloud1D("ECeren", "uncorrected Cerenkov Cloud", 100000, "autoConvert = false");
+            Edep = aida.histogramFactory().createCloud1D("Edep", "Energy Cloud", 100000, "autoConvert = false");
+            Edep_cor = aida.histogramFactory().createCloud1D("Edep_cor", "corrected Energy Cloud", 100000, "autoConvert = false");
+            Edep_dcor = aida.histogramFactory().createCloud1D("Edep_dcor", "dual readout corrected Energy Cloud", 100000, "autoConvert = false");
+            Edep_scor = aida.histogramFactory().createCloud1D("Edep_scor", "single readout corrected Energy Cloud", 100000, "autoConvert = false");
+            Eceren_cor = aida.histogramFactory().createCloud1D("Eceren_cor", "corrected Cerenkov Cloud", 100000, "autoConvert = false");
+            System.out.println("DirName:  " + DirName);
+        }
+
+        double dif_curv = 0.0;
+        if (Chrg != 0.0) {
+            double rad_curv = P_in * 667.128;
+            double len_curv = rad_curv * 2.0 * Math.asin(0.667 / rad_curv);
+            dif_curv = len_curv - 1.334;
+        }
+        List<List<SimCalorimeterHit>> simCalorimeterHitCollections = event.get(SimCalorimeterHit.class);
+        double sumEEdep = 0.0;
+        double sumECeren = 0.0;
+        for (List<SimCalorimeterHit> simCalorimeterHits : simCalorimeterHitCollections) {
+            String CollectionName = event.getMetaData(simCalorimeterHits).getName();
+            if (CollectionName.contains("Edep")) {
+                for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+                    if (cutCon.T_thresh > 0.0 || cutCon.S_thresh > 0.0) {
+                        Hep3Vector diffvec = VecOp.sub(calorimeterHit.getPositionVec(), origin);
+                        int ncontribs = calorimeterHit.getMCParticleCount();
+                        for (int i=0; i<ncontribs; i++) {
+                            double tdep = calorimeterHit.getContributedTime(i);
+                            double tofl = tdep - (diffvec.magnitude()+dif_curv)/299.79;
+                            double edep = calorimeterHit.getContributedEnergy(i);
+                            if ( (tofl < cutCon.T_thresh) && (edep > cutCon.S_thresh/1000.0) ) {
+                                sumEEdep += edep;
+                            }
+                        }
+                    }
+                    else {
+                        sumEEdep = sumEEdep + calorimeterHit.getRawEnergy();
+                    }
+                }
+            } // end if Edep
+            if (CollectionName.contains("Opti")) {
+                for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+                    if (cutCon.T_thresh > 0.0 || cutCon.C_thresh > 0.0) {
+                        Hep3Vector diffvec = VecOp.sub(calorimeterHit.getPositionVec(), origin);
+                        int ncontribs = calorimeterHit.getMCParticleCount();
+                        for (int i=0; i<ncontribs; i++) {
+                            double tdep = calorimeterHit.getContributedTime(i);
+                            double tofl = tdep - (diffvec.magnitude()+dif_curv)/299.79;
+                            double edep = calorimeterHit.getContributedEnergy(i);
+                            if ( (tofl < cutCon.T_thresh) && (edep > cutCon.C_thresh) ) {
+                                sumECeren += edep;
+                            }
+                        }
+                    }
+                    else {
+                        sumECeren = sumECeren + calorimeterHit.getRawEnergy();
+                    }
+                }
+            }
+        }
+        Edep.fill(sumEEdep);
+        ECeren.fill(sumECeren);
+        
+        xval[0] = sumEEdep;
+        double sumEEdep_cor = corfus.getFunction("SC").value(xval);
+        xval[0] = sumECeren;
+        double sumECeren_cor = corfus.getFunction("CC").value(xval);
+        Edep_cor.fill(sumEEdep_cor);
+        Eceren_cor.fill(sumECeren_cor);
+        
+        double ratio = sumECeren_cor / sumEEdep_cor;
+        double fraction = sumEEdep_cor / E_in;
+        xval[0] = ratio;
+        String cfname = "CS" + Ein.toString() + "gev";
+        double cfac = corfus.getFunction(cfname).value(xval);
+        Edep_dcor.fill(sumEEdep_cor / cfac);
+        String sfname = "S" + Ein.toString() + "gev";
+        double sfac = corfus.getFunction(sfname).parameter("mean");
+        Edep_scor.fill(sumEEdep_cor / sfac);
+    }
+
+    @Override
+    protected void endOfData() {
+        System.out.println("End of Data:");
+        convertandfit();
+        try {
+            aida.saveAs(AIDAFile);
+        } catch (IOException ex) {
+            Logger.getLogger(Resolution.class.getName()).log(Level.SEVERE, null, ex);
+        }
+    }
+
+    @Override
+    protected void resume() {
+        System.out.println("resume:");
+        firstEvent = true;
+        aida.cloud1D("c_Edep_energy").reset();
+    }
+
+    protected void convertandfit() {
+        System.out.println("convert and fit:");
+        IHistogram1D Edep_dcor_conv;
+        if (Edep_dcor.isConverted()) {
+            Edep_dcor_conv = Edep_dcor.histogram();
+        } else {
+            double meanc = Edep_dcor.mean();
+            double rmsc = Edep_dcor.rms();
+            double minx = meanc - nsigmas * rmsc;
+            double maxx = meanc + nsigmas * rmsc;
+            Edep_dcor.setConversionParameters(nbins, minx, maxx);
+            Edep_dcor.convertToHistogram();
+            Edep_dcor_conv = Edep_dcor.histogram();
+        }
+        gauss.setParameter("amplitude", Edep_dcor_conv.maxBinHeight());
+        gauss.setParameter("mean", Edep_dcor_conv.mean());
+        gauss.setParameter("sigma", Edep_dcor_conv.rms());
+        System.out.print(Particlename + "  " + Ein_prev + " GeV\n");
+        System.out.print(Edep_dcor_conv.maxBinHeight() + "," + Edep_dcor_conv.mean() + "," + Edep_dcor_conv.rms() + "\n");
+
+        for (int i = 0; i < Fitters.length; i++) {
+            System.out.println("Fitter:  " + Fitters[i]);
+            gauss.setParameter("amplitude", Edep_dcor_conv.maxBinHeight());
+            gauss.setParameter("mean", Edep_dcor_conv.mean());
+            gauss.setParameter("sigma", Edep_dcor_conv.rms());
+            jminuit = fitFactory.createFitter(Fitters[i], "jminuit");
+            jminuitResult = jminuit.fit(Edep_dcor_conv, gauss);
+            System.out.println("jminuit " + Fitters[i] + ":  " + jminuitResult.quality());
+            result = jminuitResult.fittedParameters();
+            errors = jminuitResult.errors();
+            double[] errsplus = jminuitResult.errorsPlus();
+            double[] errsminu = jminuitResult.errorsMinus();
+            String functionname = "Corrected fitted gauss  " + Fitters[i];
+            System.out.println(functionname);
+            functionFactory.cloneFunction(functionname, jminuitResult.fittedFunction());
+            IDataPoint dp_Correctede = dps_Correctede[i].addPoint();
+            dp_Correctede.coordinate(0).setValue(result[1]);
+            dp_Correctede.coordinate(0).setErrorPlus(Math.abs(result[2]));
+            dp_Correctede.coordinate(0).setErrorMinus(Math.abs(result[2]));
+            dp_Correctede.coordinate(1).setValue(Ein_prev);
+            dp_Correctede.coordinate(1).setErrorPlus(Ein_prev * 0.001);
+            dp_Correctede.coordinate(1).setErrorMinus(Ein_prev * 0.001);
+            point_Correctede[i]++;
+            IDataPoint dp_Resolution = dps_DRResolution[i].addPoint();
+            double sigma = Math.abs(result[2]);
+            dp_Resolution.coordinate(0).setValue(1.0/Math.sqrt(Ein_prev));
+            dp_Resolution.coordinate(1).setValue(Math.abs(sigma)/Ein_prev);
+            dp_Resolution.coordinate(1).setErrorMinus(errsplus[2]/Ein_prev);
+            dp_Resolution.coordinate(1).setErrorPlus(errsminu[2]/Ein_prev);
+            IDataPoint dp_DCorrFrac = dps_DCorrFrac[i].addPoint();
+            double mean = result[1];
+            dp_DCorrFrac.coordinate(0).setValue(Ein_prev);
+            dp_DCorrFrac.coordinate(1).setValue(mean/Ein_prev);
+            dp_DCorrFrac.coordinate(1).setErrorMinus(sigma);
+            dp_DCorrFrac.coordinate(1).setErrorPlus(sigma);
+            
+            IHistogram1D Edep_scor_conv;
+            if (Edep_scor.isConverted()) {
+                Edep_scor_conv = Edep_scor.histogram();
+            } else {
+                double meanc = Edep_scor.mean();
+                double rmsc = Edep_scor.rms();
+                double minx = meanc - nsigmas * rmsc;
+                double maxx = meanc + nsigmas * rmsc;
+                Edep_scor.setConversionParameters(nbins, minx, maxx);
+                Edep_scor.convertToHistogram();
+                Edep_scor_conv = Edep_scor.histogram();
+            }
+            gauss.setParameter("amplitude", Edep_scor_conv.maxBinHeight());
+            gauss.setParameter("mean", Edep_scor_conv.mean());
+            gauss.setParameter("sigma", Edep_scor_conv.rms());
+            System.out.print(Particlename + "  " + Ein_prev + " GeV\n");
+            System.out.print(Edep_scor_conv.maxBinHeight() + "," + Edep_scor_conv.mean() + "," + Edep_scor_conv.rms() + "\n");
+            jminuit = fitFactory.createFitter("Chi2", "jminuit");
+            jminuitResult = jminuit.fit(Edep_dcor_conv, gauss);
+            System.out.println("jminuit " + Fitters[i] + ":  " + jminuitResult.quality());
+            result = jminuitResult.fittedParameters();
+            errors = jminuitResult.errors();
+            IDataPoint dp_SRResolution = dps_SRResolution[i].addPoint();
+            sigma = Math.abs(result[2]);
+            dp_Resolution.coordinate(0).setValue(1.0/Math.sqrt(Ein_prev));
+            dp_Resolution.coordinate(1).setValue(Math.abs(sigma)/Ein_prev);
+            dp_Resolution.coordinate(1).setErrorMinus(errsplus[2]/Ein_prev);
+            dp_Resolution.coordinate(1).setErrorPlus(errsminu[2]/Ein_prev);
+        }
+        
+
+        try {
+            aida.saveAs(AIDAFile);
+        } catch (IOException ex) {
+            Logger.getLogger(Resolution.class.getName()).log(Level.SEVERE, null, ex);
+        }
+    }
+
+    public void setCorrFu(CorrectionFunctionSet corrs) {
+        System.out.println("DualCorrection:setCorrFu");
+        this.corfus = corrs;
+    }
+    
+    public void setCutConfiguration(cutConfiguration cutco) {
+        System.out.println("DualCorrection:setCutConfiguration");
+        this.cutCon = cutco;
+    }
+    
+    public void setMyAIDAFilename(String AIDAFile) {
+        System.out.println("DualCorrection:setMyAIDAFilename");
+        this.AIDAFile = AIDAFile;
+        System.out.println(AIDAFile);
+    }
+}

mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
proton_calibration_steering.xml added at 1.1
diff -N proton_calibration_steering.xml
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ proton_calibration_steering.xml	18 Apr 2014 20:52:20 -0000	1.1
@@ -0,0 +1,33 @@
+<lcsim xmlns:xs="http://www.w3.org/2001/XMLSchema-instance"
+       xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/lcsim/1.0/lcsim.xsd">
+    
+    <classpath>
+        <jar>/home/aconway/NetBeansProjects/mcd-analysis/target/mcd-analysis-1.0-SNAPSHOT.jar</jar>
+    </classpath>
+
+    <inputFiles>
+        <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-1gev.slcio</file>
+        <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-2gev.slcio</file>
+        <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-5gev.slcio</file>
+        <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-10gev.slcio</file>
+        <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-25gev.slcio</file>
+        <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-50gev.slcio</file>
+    </inputFiles>
+    <control>
+        <numberOfEvents>-1</numberOfEvents>
+    </control>
+
+    <execute>
+        <driver name="DRCalibrationDriver"/>
+    </execute>
+    <drivers>
+        <driver name="DRCalibrationDriver" 
+                type="org.lcsim.mcd.analysis.DRCorrection.alexDRCorr.DRCalibrationDriver">
+            <myP_name>proton</myP_name>
+            <myS_thresh>0.0</myS_thresh>
+            <myC_thresh>0.0</myC_thresh>
+            <myT_thresh>3.0</myT_thresh>
+            <myTag>default</myTag>
+        </driver>
+    </drivers>
+</lcsim>

mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
cutConfiguration.java added at 1.1
diff -N cutConfiguration.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ cutConfiguration.java	18 Apr 2014 20:52:20 -0000	1.1
@@ -0,0 +1,83 @@
+/*
+ * To change this license header, choose License Headers in Project Properties.
+ * To change this template file, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.mcd.analysis.DRCorrection.alexDRCorr;
+
+/**
+ *
+ * @author aconway
+ */
+public class cutConfiguration {
+    public Double C_thresh;
+    public Double S_thresh;
+    public Double T_thresh;
+    public String P_name;
+    public String tag;
+    
+    public cutConfiguration(Double Ct, Double St, Double Tt, String name) {
+        C_thresh = Ct;
+        S_thresh = St;
+        T_thresh = Tt;
+        P_name = name;
+        tag = "default";
+    }
+    
+    public cutConfiguration(Double Ct, Double St, Double Tt, String name, String t) {
+        C_thresh = Ct;
+        S_thresh = St;
+        T_thresh = Tt;
+        P_name = name;
+        tag = t;
+    }
+    
+    @Override
+    public boolean equals(Object o) {
+        if ((o instanceof cutConfiguration) &&
+                (((cutConfiguration) o).P_name.equals(this.P_name)) &&
+                (((cutConfiguration) o).tag.equals(this.tag)) &&
+                (((cutConfiguration) o).C_thresh.compareTo(this.C_thresh) == 0) &&
+                (((cutConfiguration) o).S_thresh.compareTo(this.S_thresh) == 0) &&
+                (((cutConfiguration) o).T_thresh.compareTo(this.T_thresh) == 0))
+        {
+            return true;
+        } else {
+            return false;
+        }
+    }
+    
+    @Override
+    public int hashCode() {
+        int hash = 7;
+        hash = 53 * hash + (this.P_name != null ? this.P_name.hashCode() : 0);
+        hash = 53 * hash + (this.C_thresh != null ? this.C_thresh.hashCode() : 0);
+        hash = 53 * hash + (this.S_thresh != null ? this.S_thresh.hashCode() : 0);
+        hash = 53 * hash + (this.T_thresh != null ? this.T_thresh.hashCode() : 0);
+        return hash;
+    }
+    
+    @Override
+    public String toString() {
+        String rv = P_name + "-" +
+                String.valueOf(this.T_thresh) + "ns-" +
+                String.valueOf(this.S_thresh) + "mev-" +
+                String.valueOf(this.C_thresh);
+        if (!this.tag.equals("default")) {
+            rv  = rv + "-" + this.tag;
+        }
+        return rv;
+    }
+    
+    public String toStringAnyP() {
+        String rv = "any" + "-" +
+                String.valueOf(this.T_thresh) + "ns-" +
+                String.valueOf(this.S_thresh) + "mev-" +
+                String.valueOf(this.C_thresh);
+        if (!this.tag.equals("default")) {
+            rv  = rv + "-" + this.tag;
+        }
+        return rv;
+    }
+}

mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
electron_calibration_steering.xml added at 1.1
diff -N electron_calibration_steering.xml
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ electron_calibration_steering.xml	18 Apr 2014 20:52:20 -0000	1.1
@@ -0,0 +1,32 @@
+<lcsim xmlns:xs="http://www.w3.org/2001/XMLSchema-instance"
+       xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/lcsim/1.0/lcsim.xsd">
+    
+    <classpath>
+        <jar>/home/aconway/NetBeansProjects/mcd-analysis/target/mcd-analysis-1.0-SNAPSHOT.jar</jar>
+    </classpath>
+    <inputFiles>
+        <file>/home/aconway/fermi/eventdata/slcio/singleparticle/electron/electron-90deg-1gev.slcio</file>
+        <file>/home/aconway/fermi/eventdata/slcio/singleparticle/electron/electron-90deg-2gev.slcio</file>
+        <file>/home/aconway/fermi/eventdata/slcio/singleparticle/electron/electron-90deg-5gev.slcio</file>
+        <file>/home/aconway/fermi/eventdata/slcio/singleparticle/electron/electron-90deg-10gev.slcio</file>
+        <file>/home/aconway/fermi/eventdata/slcio/singleparticle/electron/electron-90deg-25gev.slcio</file>
+        <file>/home/aconway/fermi/eventdata/slcio/singleparticle/electron/electron-90deg-50gev.slcio</file>
+        <file>/home/aconway/fermi/eventdata/slcio/singleparticle/electron/electron-90deg-100gev.slcio</file>
+    </inputFiles>
+    <control>
+        <numberOfEvents>-1</numberOfEvents>
+    </control>
+    <execute>
+        <driver name="ElectronCalibrationDriver"/>
+    </execute>
+    <drivers>
+        <driver name="ElectronCalibrationDriver"
+                type="org.lcsim.mcd.analysis.DRCorrection.alexDRCorr.ElectronCalibrationDriver">
+            <myP_name>any</myP_name>
+            <myS_thresh>0.0</myS_thresh>
+            <myC_thresh>0.0</myC_thresh>
+            <myT_thresh>0.0</myT_thresh>
+            <myTag>default</myTag>
+        </driver>
+    </drivers>
+</lcsim>

mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
DRCalibrationDriver.java added at 1.1
diff -N DRCalibrationDriver.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ DRCalibrationDriver.java	18 Apr 2014 20:52:20 -0000	1.1
@@ -0,0 +1,65 @@
+package org.lcsim.mcd.analysis.DRCorrection.alexDRCorr;
+
+import org.lcsim.util.Driver;
+
+/**
+ *
+ * @author aconway
+ */
+public class DRCalibrationDriver extends Driver {
+    
+    private String P_name = "any";
+    private String tag = "default";
+    private Double C_thresh = 0.0;
+    private Double S_thresh = 0.0;
+    private Double T_thresh = 0.0;
+    private String AIDAFile = null;
+    
+    DualCorrection dual;
+    
+    public DRCalibrationDriver() {
+        dual = new DualCorrection();
+        add(dual);
+    }
+    
+    @Override
+    public void startOfData() {
+        System.out.println("DRCalibrationDriver:startOfData");
+        cutConfiguration cutco;
+        if (tag.equals("default")) {
+            cutco = new cutConfiguration(C_thresh, S_thresh, T_thresh, P_name);
+        } else {
+            cutco = new cutConfiguration(C_thresh, S_thresh, T_thresh, P_name, tag);
+        }
+        AIDAFile = "DR-" + cutco.toString() + ".aida";
+        dual.setMyAIDAFilename(AIDAFile);
+        dual.setCutConfiguration(cutco);
+        dual.setCorrFu(new CorrectionFunctionSet(cutco, true));
+        dual.startOfData();
+    }
+    
+    public void setMyP_name(String pname) {
+        System.out.println("ElectronCalibrationDriver:setMyP_name = " + pname);
+        this.P_name = pname;
+    }
+    
+    public void setMyTag(String tag) {
+        System.out.println("ElectronCalibrationDriver:setMyTag = " + tag);
+        this.tag = tag;
+    }
+    
+    public void setMyC_thresh(double cthresh) {
+        System.out.println("ElectronCalibrationDriver:setMyC_thresh = " + cthresh);
+        this.C_thresh = cthresh;
+    }
+    
+    public void setMyS_thresh(double sthresh) {
+        System.out.println("ElectronCalibrationDriver:setMyS_thresh = " + sthresh);
+        this.S_thresh = sthresh;
+    }
+    
+    public void setMyT_thresh(double tthresh) {
+        System.out.println("ElectronCalibrationDriver:setMyT_thresh = " + tthresh);
+        this.T_thresh = tthresh;
+    }
+}

mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection
DualCorrection.java 1.2 -> 1.3
diff -u -r1.2 -r1.3
--- DualCorrection.java	7 Apr 2014 15:13:48 -0000	1.2
+++ DualCorrection.java	18 Apr 2014 20:52:20 -0000	1.3
@@ -434,15 +434,24 @@
             if (entries > 300) {
                 for (int ii = 0; ii < Fitters.length; ii++) {
                     System.out.println("Fitter:  " + Fitters[ii]);
+                    double meanguess = 0;
+                    for (int j=0; j<conv_slices[i].axis().bins(); j++) {
+                        if (conv_slices[i].binHeight(j) == conv_slices[i].maxBinHeight()) {
+                            meanguess = conv_slices[i].axis().binCenter(j);
+                        }
+                    }
                     gauss.setParameter("amplitude", conv_slices[i].maxBinHeight());
-                    gauss.setParameter("mean", conv_slices[i].mean());
-                    gauss.setParameter("sigma", conv_slices[i].rms());
+                    gauss.setParameter("mean", meanguess);
+                    gauss.setParameter("sigma", conv_slices[i].rms()/2.0);
                     jminuit = fitFactory.createFitter(Fitters[ii], "jminuit");
                     try {
                         jminuitResult = jminuit.fit(conv_slices[i], gauss);
                     } catch (RuntimeException ex) {
                         System.out.println("Runtime Exception caught in fit!");
                         System.out.println(ex.toString());
+                    } catch (Exception ex) {
+                        System.out.println("Exception caught in f!");
+                        System.out.println(ex.toString());
                     }
                     System.out.println("jminuit " + Fitters[ii] + ":  " + jminuitResult.quality());
                     result = jminuitResult.fittedParameters();

mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection
ElectronCorrection.java 1.2 -> 1.3
diff -u -r1.2 -r1.3
--- ElectronCorrection.java	7 Apr 2014 15:13:48 -0000	1.2
+++ ElectronCorrection.java	18 Apr 2014 20:52:20 -0000	1.3
@@ -391,4 +391,4 @@
         this.file_name = file_name;
         System.out.println(file_name);
     }
-}
+}
\ No newline at end of file

mcd-analysis/src/main/java/org/lcsim/mcd/analysis
ILCstdhepdriver.java 1.1 -> 1.2
diff -u -r1.1 -r1.2
--- ILCstdhepdriver.java	7 Apr 2014 15:13:49 -0000	1.1
+++ ILCstdhepdriver.java	18 Apr 2014 20:52:20 -0000	1.2
@@ -1,17 +1,15 @@
-/*
- * To change this license header, choose License Headers in Project Properties.
- * To change this template file, choose Tools | Templates
- * and open the template in the editor.
- */
-
 //package org.lcsim.mcd.analysis;
+import hep.aida.IAxis;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
 import java.util.List;
 import java.util.ArrayList;
 import org.lcsim.event.EventHeader;
 import org.lcsim.event.MCParticle;
-import org.lcsim.event.base.BaseMCParticle;
 import org.lcsim.util.Driver;
 import org.lcsim.util.aida.AIDA;
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogram2D;
 
 /**
  *
@@ -21,6 +19,26 @@
 public class ILCstdhepdriver extends Driver {
     private AIDA aida = AIDA.defaultInstance();
     
+    private List<String> hadrons = new ArrayList<String>();
+    
+    IHistogram1D allhadsen = aida.histogram1D("sum energy wtd by frac energy",
+            25,0,50);
+    IHistogram2D allhadsfrac = aida.histogram2D("energy vs frac energy",
+            20,0,50,
+            25,0.0,0.5);
+    
+    public ILCstdhepdriver() {
+        hadrons.add("p");
+        hadrons.add("p_bar");
+        hadrons.add("pi+");
+        hadrons.add("pi-");
+        hadrons.add("n");
+        hadrons.add("n_bar");
+        hadrons.add("K+");
+        hadrons.add("K-");
+        hadrons.add("K0_L");
+    }
+    
     protected void process(EventHeader event) {
         List<MCParticle> particles = event.get(MCParticle.class,"MCParticle");
         List<MCParticle> children = new ArrayList<MCParticle>();
@@ -40,14 +58,105 @@
         List<MCParticle> fsChildren = getAllFSChildren(higgs);
 
         String pfix = "/"+decayname+"/";
+        double sumEShower = 0.0;
+        double sumEBarrel = 0.0;
+        double sumEHadron = 0.0;
+        for (MCParticle p : fsChildren) {
+            sumEShower += p.getEnergy();
+            double logen = Math.log10(p.getEnergy());
+            double logp = Math.log10(p.getMomentum().magnitude());
+            if (hitsBarrel(p)) {
+                sumEBarrel += p.getEnergy();
+                String name = p.getType().getName();
+                if (hadrons.contains(name)) {
+                    sumEHadron += p.getEnergy();
+                }
+                aida.cloud1D(pfix+"hits/"+p.getType().getName()+"/energy").fill(p.getEnergy());
+                aida.cloud1D(pfix+"hits/"+p.getType().getName()+"/momentum").fill(p.getMomentum().magnitude());
+                aida.cloud1D(pfix+"hits/"+p.getType().getName()+"/log energy").fill(logen);
+                aida.cloud1D(pfix+"hits/"+p.getType().getName()+"/log momentum").fill(logp);
+            } else {
+                aida.cloud1D(pfix+"misses/"+p.getType().getName()+"/energy").fill(p.getEnergy());
+                aida.cloud1D(pfix+"misses/"+p.getType().getName()+"/momentum").fill(p.getMomentum().magnitude());
+                aida.cloud1D(pfix+"misses/"+p.getType().getName()+"/log energy").fill(logen);
+                aida.cloud1D(pfix+"misses/"+p.getType().getName()+"/log momentum").fill(logp);
+            }
+        }
+        aida.cloud1D(pfix+"tot shower energy").fill(sumEShower);
+        aida.cloud1D(pfix+"frac barrel energy").fill(sumEBarrel/sumEShower);
+        List<String> pnames = new ArrayList<String>();
         for (MCParticle p : fsChildren) {
             String particlename = p.getType().getName();
             aida.cloud1D(pfix+particlename+"/momentum").fill(p.getMomentum().magnitude());
             aida.cloud1D(pfix+particlename+"/energy").fill(p.getEnergy());
             aida.cloud1D(pfix+particlename+"/log momentum").fill(Math.log10(p.getMomentum().magnitude()));
             aida.cloud1D(pfix+particlename+"/log energy").fill(Math.log10(p.getEnergy()));
-            // TODO: look at relative composition of shower particles and energies. Is most of the energy in an event in low energy hadrons?
-            
+            if (hitsBarrel(p)) {
+                aida.cloud1D(pfix+particlename+"/log energy wtd by frac energy").fill(Math.log10(p.getEnergy()), p.getEnergy()/sumEBarrel);
+                aida.cloud1D(pfix+particlename+"/sum log energy wtd by frac energy").fill(Math.log10(p.getEnergy()), p.getEnergy()/sumEBarrel);
+                aida.cloud1D(pfix+"all/sum log energy wtd by frac energy").fill(Math.log10(p.getEnergy()), p.getEnergy()/sumEBarrel);
+                aida.cloud1D(pfix+particlename+"/energy wtd by frac energy").fill(p.getEnergy(), p.getEnergy()/sumEBarrel);
+                aida.cloud1D(pfix+particlename+"/sum energy wtd by frac energy").fill(p.getEnergy(), p.getEnergy()/sumEBarrel);
+                aida.cloud1D(pfix+"all/sum energy wtd by frac energy").fill(p.getEnergy(), p.getEnergy()/sumEBarrel);
+                if (!pnames.contains(particlename)) {
+                    pnames.add(particlename);
+                }
+                if (hadrons.contains(particlename)) {
+                    aida.cloud1D(pfix+"allhads/sum energy wtd by frac energy").fill(p.getEnergy(), p.getEnergy()/sumEHadron);
+                    allhadsen.fill(p.getEnergy(),p.getEnergy()/sumEHadron);
+                }
+            }
+        }
+        for (String pname : pnames) {
+            String opname = pfix+pname+"/energy wtd by frac energy";
+            if (!aida.cloud1D(opname).isConverted()) {
+                aida.cloud1D(opname).convertToHistogram();
+            }
+            IHistogram1D ophist = aida.cloud1D(opname).histogram();
+            IAxis opax = ophist.axis();
+            for (int i=0; i<opax.bins(); i++) {
+                double frac = ophist.binHeight(i);
+                double en = ophist.binMean(i);
+                if (frac > 0.005) {
+                    aida.cloud2D(pfix+pname+"/log energy vs frac energy").fill(en,frac);
+                }
+            }
+            aida.cloud1D(opname).reset();
+        }
+        String opname = pfix+"allhads/sum energy wtd by frac energy";
+        if (!aida.cloud1D(opname).isConverted()) {
+            aida.cloud1D(opname).convertToHistogram();
+        }
+        IHistogram1D ophist = aida.cloud1D(opname).histogram();
+        IAxis opax = ophist.axis();
+        for (int i=0; i<opax.bins(); i++) {
+            double frac = ophist.binHeight(i);
+            double en = ophist.axis().binCenter(i);
+            if (frac > 0.00) {
+                aida.cloud2D(pfix+"allhads/log energy vs frac energy").fill(en,frac);
+            }
+        }
+        aida.cloud1D(opname).reset();
+        for (int i=0; i<allhadsen.axis().bins(); i++) {
+            double frac = allhadsen.binHeight(i);
+            double en = allhadsen.axis().binCenter(i);
+            if (frac > 0.00005) {
+                aida.cloud2D("cloud energy vs frac energy").fill(en, frac);
+                allhadsfrac.fill(en,frac);
+            }
+        }
+        allhadsen.reset();
+    }
+    
+    public boolean hitsBarrel (MCParticle mcp) {
+        if (mcp.getCharge() == 0.0) {
+            return true;
+        } else {
+            Hep3Vector origin = new BasicHep3Vector(mcp.getOriginX(),mcp.getOriginY(),0.0);
+            Hep3Vector pt = new BasicHep3Vector(mcp.getMomentum().x(),mcp.getMomentum().y(),0.0);
+            // diameter of curvature in 5T field in mm
+            double diam = pt.magnitude()*1334.256;        
+            return diam > 1350.0 - origin.magnitude();
         }
     }
     
CVSspam 0.2.12


Use REPLY-ALL to reply to list

To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1