diff -N DualCorrection.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ DualCorrection.java	25 Sep 2009 15:07:05 -0000	1.1.2.1
@@ -0,0 +1,395 @@

+package org.lcsim.contrib.HansWenzel.DualCorrection;
+
+/**
+ *
+ * @author wenzel
+ * @date
+ */
+import hep.aida.ICloud1D;
+import hep.aida.ICloud2D;
+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.aida.IProfile1D;
+import java.io.BufferedWriter;
+import java.io.FileWriter;
+import java.io.IOException;
+import java.text.DecimalFormat;
+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 DualCorrection extends Driver {
+
+    private AIDA aida = AIDA.defaultInstance();
+    String AIDAFile;
+    String file_name;
+    String Physics_List;
+    double Gun_radius;
+    FileWriter fstream;
+    BufferedWriter out;
+    IFunctionFactory functionFactory;
+    IFitFactory fitFactory;
+    IDataPointSetFactory dpsf;
+    boolean first;
+    boolean firstEvent;
+    double E_in;
+    Integer Ein;
+    Integer Ein_prev;
+    double E_kin;
+    String Particlename;
+    ICloud1D Edep;
+    ICloud1D Eceren;
+    ICloud1D Edep_cor;
+    ICloud1D Edep_dcor;
+    ICloud1D Eceren_cor;
+    ICloud2D c_efrac_ratio;
+    ICloud1D[] slice_comb;
+    ICloud1D[] slice;
+    IHistogram1D[] conv_slice_comb;
+    IHistogram1D[] conv_slice;
+    IFunction gauss;
+    IProfile1D prof_combined;
+    double nsigmas;
+    int nbins;
+    IDataPointSet dps_emean = null;
+    IDataPointSet dps_cerenemean = null;
+    IDataPointSet dps_cerenechi2 = null;
+    IDataPointSet dps_cerenels = null;
+    IDataPointSet dps_cerenecc = null;
+    IDataPointSet[] dps_cerene = null;
+    int point_cerene[];
+    IDataPointSet dps_cereneml = null;
+    int point;
+    int point_cerenemean;
+    int point_cerenechi2;
+    int point_cerenels;
+    int point_cerenecc;
+    int point_cereneml;
+    double[] result;
+    double errors[];
+    String[] Fitters = {"Chi2", "leastsquares", "bml", "cleverchi2"};
+    IFitter jminuit;
+    IFitResult jminuitResult;
+    IFunction E_cor;
+    IFunction C_cor;
+    IFunction D_cor;
+    double xval[] = {10.};
+    double binwidth = 0.04;
+    int maxbin = (int) Math.rint(1.0 / binwidth);
+
+    @Override
+    protected void startOfData() {
+        System.out.println("Start of Data:");
+        fitFactory = aida.analysisFactory().createFitFactory();
+        functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+        gauss = functionFactory.createFunctionByName("gauss", "G");
+        E_cor = Electron_Correction();
+        C_cor = Cerenkov_Correction();
+        D_cor = Dual_Correction();
+        aida.profile1D("ratio", 20, 0.4, 1.);
+        prof_combined = aida.histogramFactory().createProfile1D("ratio", "ratio combined", 20, 0.4, 1.);
+        c_efrac_ratio = aida.histogramFactory().createCloud2D("c_efrac_ratio", "c_efrac_ratio combined ", 500000, "autoConvert = false");
+        slice_comb = new ICloud1D[maxbin + 2];
+        conv_slice_comb = new IHistogram1D[maxbin + 2];
+        for (int ibin = 0; ibin < maxbin + 2; ibin++) {
+            String cloudy = "ratio_" + ibin;
+            slice_comb[ibin] = aida.histogramFactory().createCloud1D(cloudy, cloudy, 500000, "autoConvert = false");
+        }
+        Ein_prev = 0;
+        firstEvent = true;
+        first = true;
+        fstream = null;
+        try {
+            fstream = new FileWriter(file_name);
+        } catch (IOException ex) {
+            Logger.getLogger(DualCorrection.class.getName()).log(Level.SEVERE, null, ex);
+        }
+        out = new BufferedWriter(fstream);
+    }
+
+    @Override
+    protected void process(EventHeader event) {
+        E_in = 0.0;
+        E_kin = 0.0;
+        String DirName = null;
+        List<MCParticle> particles = event.get(MCParticle.class, event.MC_PARTICLES);
+
+        for (MCParticle particle : particles) {
+            if (particle.getGeneratorStatus() == 0) {
+                E_in = E_in + particle.getEnergy();
+                E_kin = E_kin + particle.getEnergy() - particle.getMass();
+                if (particle.getProductionTime() == 0.0) {
+                    Particlename = particle.getType().getName();
+                }
+            }
+            break;
+        }
+        Ein = (int) Math.floor(E_kin + 0.5d);
+        String E_str = Ein.toString();
+        DirName = Particlename.concat(E_str);
+        if (Ein != Ein_prev) {
+            if (first) {
+                first = false;
+            } else {
+                fitprofile(aida.profile1D("ratio"));
+                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");
+            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");
+            Eceren_cor = aida.histogramFactory().createCloud1D("Eceren_cor", "corrected Cerenkov Cloud", 100000, "autoConvert = false");
+            slice = new ICloud1D[maxbin + 2];
+            conv_slice = new IHistogram1D[maxbin + 2];
+            for (int ibin = 0; ibin < maxbin + 2; ibin++) {
+                String cloudy = "ratio_" + ibin;
+                slice[ibin] = aida.histogramFactory().createCloud1D(cloudy, cloudy, 500000, "autoConvert = false");
+            }
+            System.out.println("DirName:  " + DirName);
+        }
+        List<List<SimCalorimeterHit>> simCalorimeterHitCollections = event.get(SimCalorimeterHit.class);
+        double sumEEdep = 0.0;
+        double sumECeren = 0.0;
+
+        aida.profile1D("ratio", 20, 0.4, 1.);
+        for (List<SimCalorimeterHit> simCalorimeterHits : simCalorimeterHitCollections) {
+            String CollectionName = event.getMetaData(simCalorimeterHits).getName();
+            if (CollectionName.substring(0, 5).equals("Edep_")) {
+                for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+                    sumEEdep = sumEEdep + calorimeterHit.getRawEnergy();
+                }
+            } // end if Edep
+            if (CollectionName.substring(0, 6).equals("Ceren_")) {
+                for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+                    sumECeren = sumECeren + calorimeterHit.getRawEnergy();
+                }
+            } // end if Ceren
+        }     // end loop over calorimeter hit collections
+
+        Edep.fill(sumEEdep);
+        xval[0] = sumEEdep;
+        double sumEEdep_cor = E_cor.value(xval);
+        Edep_cor.fill(sumEEdep_cor);
+        Eceren.fill(sumECeren);
+        xval[0] = sumECeren;
+        double sumECeren_cor = C_cor.value(xval);
+        Eceren_cor.fill(sumECeren_cor);
+        double ratio = sumECeren_cor / sumEEdep_cor;
+        double fraction = sumEEdep_cor / E_in;
+        aida.cloud1D("frac").fill(fraction);
+        aida.cloud1D("c_CerenEdep_ratio").fill(ratio);
+        aida.cloud2D("c_Ceren_vs_Edep").fill(sumECeren_cor, sumEEdep_cor);
+        aida.cloud2D("c_efrac_ratio").fill(ratio, fraction);
+        aida.profile1D("ratio").fill(ratio, fraction);
+        xval[0] = ratio;
+        double cfac = D_cor.value(xval);
+        Edep_dcor.fill(sumEEdep_cor / cfac);
+        prof_combined.fill(ratio, fraction);
+        c_efrac_ratio.fill(ratio, fraction);
+
+        int bin = (int) Math.rint(ratio / binwidth);
+        if (bin < 0) {
+            bin = 0;
+        }
+        if (bin > maxbin) {
+            bin = maxbin + 1;
+        }
+        slice_comb[bin].fill(fraction);
+        slice[bin].fill(fraction);
+    }
+
+    @Override
+    protected void endOfData() {
+        System.out.println("End of Data:");
+        fitprofile(aida.profile1D("ratio"));
+        convertandfit();
+        aida.tree().cd("/");
+        fitprofile(prof_combined);
+        try {
+            out.close();
+        } catch (IOException ex) {
+            Logger.getLogger(DualCorrection.class.getName()).log(Level.SEVERE, null, ex);
+        }
+        try {
+            aida.saveAs(AIDAFile);
+        } catch (IOException ex) {
+            Logger.getLogger(DualCorrection.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 fitprofile(IProfile1D profile) {
+        System.out.println("fit profile:");
+        IFunction poly = functionFactory.createFunctionByName("poly", "p3");
+        jminuit = fitFactory.createFitter("Chi2", "jminuit");
+        poly.setParameter("p0", 0.63);
+        poly.setParameter("p1", 0.54);
+        poly.setParameter("p2", -0.58);
+        poly.setParameter("p3", 0.5);
+        int profbins = profile.axis().bins();
+        double profxmin = profile.axis().lowerEdge();
+        double profxmax = profile.axis().upperEdge();
+/*        for (int i = 0; i < profbins; i++) {
+            System.out.println(profile.binEntries(i));
+        }
+*/
+        jminuitResult = jminuit.fit(profile, poly);
+        System.out.println("jminuit Chi2=" + jminuitResult.quality());
+        result = jminuitResult.fittedParameters();
+        functionFactory.cloneFunction("fitted poly (jminuit)", jminuitResult.fittedFunction());
+        System.out.println("correction function:  " + result[0] + "," + result[1] + "," + result[2] + "," + result[3]);
+        try {
+            out.write(Particlename + "  " + Ein_prev + " GeV\n");
+            out.write(result[0] + "," + result[1] + "," + result[2] + "," + result[3]);
+        } catch (IOException ex) {
+            Logger.getLogger(DualCorrection.class.getName()).log(Level.SEVERE, null, ex);
+        }
+        try {
+            aida.saveAs(AIDAFile);
+        } catch (IOException ex) {
+            Logger.getLogger(DualCorrection.class.getName()).log(Level.SEVERE, null, ex);
+        }
+    }
+
+    protected void convertandfit() {
+        System.out.println("convert and fit:");
+        try {
+            out.write(Particlename + "  " + Ein_prev + " GeV\n");
+        } catch (IOException ex) {
+            Logger.getLogger(DualCorrection.class.getName()).log(Level.SEVERE, null, ex);
+        }
+        for (int i = 0; i < slice.length; i++) {
+            if (slice[i].isConverted()) {
+                conv_slice[i] = slice[i].histogram();
+            } else {
+                System.out.println("Converting EDep");
+                double meanc = slice[i].mean();
+                double rmsc = slice[i].rms();
+                nsigmas = 3.;
+                nbins = 100;
+                double minx = meanc - nsigmas * rmsc;
+                double maxx = meanc + nsigmas * rmsc;
+                slice[i].setConversionParameters(nbins, minx, maxx);
+                slice[i].convertToHistogram();
+                conv_slice[i] = slice[i].histogram();
+            }
+
+            gauss.setParameter("amplitude", conv_slice[i].maxBinHeight());
+            gauss.setParameter("mean", conv_slice[i].mean());
+            gauss.setParameter("sigma", conv_slice[i].rms());
+            System.out.println(conv_slice[i].maxBinHeight() + " , " + conv_slice[i].mean() + " , " + conv_slice[i].rms());
+        }
+        /*
+        for (int i = 0; i < Fitters.length; i++) {
+        System.out.println("Fitter:  " + Fitters[i]);
+        gauss.setParameter("amplitude", Eceren_conv.maxBinHeight());
+        gauss.setParameter("mean", Eceren_conv.mean());
+        gauss.setParameter("sigma", Eceren_conv.rms());
+        jminuit = fitFactory.createFitter(Fitters[i], "jminuit");
+        jminuitResult = jminuit.fit(Eceren_conv, gauss);
+        System.out.println("jminuit " + Fitters[i] + ":  " + jminuitResult.quality());
+        result = jminuitResult.fittedParameters();
+        errors = jminuitResult.errors();
+        String functionname = "ceren fitted gauss  " + Fitters[i];
+        System.out.println(functionname);
+        functionFactory.cloneFunction(functionname, jminuitResult.fittedFunction());
+        System.out.println(result[0] + "," + result[1] + "," + result[2]);
+        try {
+        out.write(result[0] + "," + result[1] + "," + result[2] + "\n");
+        } catch (IOException ex) {
+        Logger.getLogger(DualCorrection.class.getName()).log(Level.SEVERE, null, ex);
+        }
+        IDataPoint dp_cerene = dps_cerene[i].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[i]++;
+        }
+         */
+        try {
+            aida.saveAs(AIDAFile);
+        } catch (IOException ex) {
+            Logger.getLogger(DualCorrection.class.getName()).log(Level.SEVERE, null, ex);
+        }
+    }
+
+    public void setMyVariable(String jpgfile) {
+        System.out.println("setMyVariable");
+        AIDAFile = jpgfile;
+        System.out.println(AIDAFile);
+    }
+
+    public void setMyFilename(String filename) {
+        System.out.println("setMyFilename");
+        file_name = filename;
+        System.out.println(file_name);
+    }
+
+    public void setMyPhysicsList(String list) {
+        System.out.println("setMyPhysicsList");
+        Physics_List = list;
+        System.out.println(list);
+    }
+
+    public void setMyradius(double radius) {
+        System.out.println("setMyPhysicsList");
+        Gun_radius = radius;
+        System.out.println(Gun_radius);
+    }
+
+    double roundTwoDecimals(double d) {
+        DecimalFormat twoDForm = new DecimalFormat("#.##");
+        return Double.valueOf(twoDForm.format(d));
+    }
+
+    IFunction Electron_Correction() {
+        IFunction ECor = functionFactory.createFunctionByName("ECor", "p1");
+        ECor.setParameter("p0", -1.9834e-3);
+        ECor.setParameter("p1", 1.0019);
+        return ECor;
+    }
+
+    IFunction Cerenkov_Correction() {
+        IFunction ECor = functionFactory.createFunctionByName("ECerenCor", "p1");
+        ECor.setParameter("p0", -2.4961e-4);
+        ECor.setParameter("p1", 7655.0);
+        return ECor;
+    }
+
+    IFunction Dual_Correction() {
+        IFunction ECor = functionFactory.createFunctionByName("DualCor", "p3");
+        ECor.setParameter("p0", 1.0900819852754733);
+        ECor.setParameter("p1", -1.72608413712137);
+        ECor.setParameter("p2", 2.602954159670887);
+        ECor.setParameter("p3", -1.0873650332172629);
+        return ECor;
+    }
+}

\ No newline at end of file