lcsim-contrib/src/main/java/org/lcsim/contrib/HansWenzel/DualCorrection
diff -N DualCorrection.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ DualCorrection.java 27 Sep 2009 18:22:03 -0000 1.2
@@ -0,0 +1,420 @@
+package org.lcsim.contrib.HansWenzel.DualCorrection;
+
+/**
+ *
+ * @author wenzel
+ * @date
+ */
+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.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 Edep_dcor_new;
+ 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_cerene;
+ int point_cerene[];
+ double[] result;
+ double errors[];
+ String[] Fitters = {"Chi2", "leastsquares", "bml", "cleverchi2"};
+ IFitter jminuit;
+ IFitResult jminuitResult;
+ IFunction E_cor;
+ IFunction C_cor;
+ IFunction D_cor;
+ IFunction D_cor_new;
+ 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();
+ dps_cerene = new IDataPointSet[Fitters.length];
+ point_cerene = new int[Fitters.length];
+ dpsf = aida.analysisFactory().createDataPointSetFactory(aida.tree());
+ functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+ gauss = functionFactory.createFunctionByName("gauss", "G");
+ E_cor = Electron_Correction();
+ C_cor = Cerenkov_Correction();
+ D_cor = Dual_Correction();
+ D_cor_new = Dual_Correction_new();
+ 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 {
+ System.out.println("Process event: ");
+ fitprofile(aida.profile1D("ratio"));
+ 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");
+ 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");
+ Edep_dcor_new = aida.histogramFactory().createCloud1D("Edep_dcor_new", "new 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);
+ cfac = D_cor_new.value(xval);
+ Edep_dcor_new.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(slice, conv_slice);
+ aida.tree().cd("/");
+ fitprofile(prof_combined);
+ convertandfit(slice_comb, conv_slice_comb);
+ 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(ICloud1D[] slices, IHistogram1D[] conv_slices) {
+ 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 < Fitters.length; i++) {
+ String dpsname = "dps_ratio_" + Fitters[i];
+ point_cerene[i] = 0;
+ dps_cerene[i] = dpsf.create(dpsname, "ratio", 2);
+ }
+ for (int i = 0; i < slices.length - 3; i++) {
+ if (slices[i].isConverted()) {
+ conv_slices[i] = slices[i].histogram();
+ } else {
+ System.out.println("Converting EDep");
+ double meanc = slices[i].mean();
+ double rmsc = slices[i].rms();
+ nsigmas = 3.;
+ 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 > 100) {
+ for (int ii = 0; ii < Fitters.length; ii++) {
+ System.out.println("Fitter: " + Fitters[ii]);
+ gauss.setParameter("amplitude", conv_slices[i].maxBinHeight());
+ gauss.setParameter("mean", conv_slices[i].mean());
+ gauss.setParameter("sigma", conv_slices[i].rms());
+ jminuit = fitFactory.createFitter(Fitters[ii], "jminuit");
+ jminuitResult = jminuit.fit(conv_slices[i], gauss);
+ System.out.println("jminuit " + Fitters[ii] + ": " + jminuitResult.quality());
+ result = jminuitResult.fittedParameters();
+ errors = jminuitResult.errors();
+ String functionname = "ratio fitted gauss slices: " + i + " " + Fitters[ii];
+ 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[ii].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]));
+ point_cerene[ii]++;
+ }
+ }
+ }
+ //Now fit the data point sets:
+ for (int ii = 0; ii < Fitters.length; ii++) {
+ 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);
+ System.out.println("Fitter: " + Fitters[ii]);
+ jminuitResult = jminuit.fit(dps_cerene[ii], poly);
+ String fname = "correction " + Fitters[ii] + " fitted poly";
+ System.out.println(fname);
+ functionFactory.cloneFunction(fname, jminuitResult.fittedFunction());
+ System.out.println(Fitters[ii] + ": " + jminuitResult.quality());
+ }
+ 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;
+ }
+
+ IFunction Dual_Correction_new() {
+ IFunction ECor = functionFactory.createFunctionByName("DualCor_new", "p3");
+ ECor.setParameter("p0", 0.62198);
+ ECor.setParameter("p1", 0.45651);
+ ECor.setParameter("p2", -0.64669);
+ ECor.setParameter("p3", 0.53322);
+ return ECor;
+ }
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/HansWenzel/DualCorrection
diff -N ElectronCorrection.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ElectronCorrection.java 27 Sep 2009 18:22:03 -0000 1.2
@@ -0,0 +1,335 @@
+package org.lcsim.contrib.HansWenzel.DualCorrection;
+
+/**
+ *
+ * @author wenzel
+ */
+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 java.io.BufferedWriter;
+import java.io.FileWriter;
+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 ElectronCorrection extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+ String AIDAFile;
+ String file_name;
+ FileWriter fstream;
+ BufferedWriter out;
+ IFunctionFactory functionFactory;
+ IFitFactory fitFactory;
+ IDataPointSetFactory dpsf;
+ IFunction gauss;
+ boolean first;
+ boolean firstEvent;
+ double E_in;
+ Integer Ein;
+ Integer Ein_prev;
+ double E_kin;
+ String Particlename;
+ ICloud1D Edep;
+ ICloud1D Eceren;
+ double nsigmas;
+ int nbins;
+ IDataPointSet dps_emean = null;
+ IDataPointSet dps_cerenemean = null;
+ IDataPointSet[] dps_cerene = null;
+ int point_cerene[];
+ int point;
+ int point_cerenemean;
+ double[] result;
+ double errors[];
+ String[] Fitters = {"Chi2", "leastsquares", "bml", "cleverchi2"};
+ IFitter jminuit;
+ IFitResult jminuitResult;
+ @Override
+ protected void startOfData() {
+ System.out.println("Start of Data:");
+ dps_cerene = new IDataPointSet[Fitters.length];
+ point_cerene = new int[Fitters.length];
+ fitFactory = aida.analysisFactory().createFitFactory();
+ functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+ dpsf = aida.analysisFactory().createDataPointSetFactory(aida.tree());
+ for (int i = 0; i < Fitters.length; i++) {
+ String dpsname = "dps_cerene_" + Fitters[i];
+ point_cerene[i] = 0;
+ dps_cerene[i] = dpsf.create(dpsname, "electron response mean", 2);
+ }
+ Ein_prev = 0;
+ firstEvent = true;
+ first = true;
+ // Create a two dimensional IDataPointSet.
+ dps_emean = dpsf.create("dps_emean", "electron response mean", 2);
+ dps_cerenemean = dpsf.create("dps_cerenemean", "electron response mean", 2);
+ point = 0;
+ point_cerenemean = 0;
+ gauss = functionFactory.createFunctionByName("gauss", "G");
+ fstream = null;
+ try {
+ fstream = new FileWriter(file_name);
+ } catch (IOException ex) {
+ Logger.getLogger(ElectronCorrection.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 {
+ 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);
+ }
+ 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.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);
+ Eceren.fill(sumECeren);
+ }
+
+ @Override
+ protected void endOfData() {
+ System.out.println("End of Data:");
+ convertandfit();
+ aida.tree().cd("/");
+ try {
+ out.close();
+ } catch (IOException ex) {
+ Logger.getLogger(ElectronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ 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());
+ System.out.println("Chi2=" + resultline.quality());
+ //
+ // now deal with the cerenkov stuff:
+ //
+ for (int i = 0; i < Fitters.length; i++) {
+ System.out.println("Fitter: " + Fitters[i]);
+ resultline = linefit.fit(dps_cerene[i], line);
+ String fname = "cerenkov " + Fitters[i] + " fitted line";
+ System.out.println(fname);
+ functionFactory.cloneFunction(fname, resultline.fittedFunction());
+ System.out.println(Fitters[i] + ": " + resultline.quality());
+ }
+ try {
+ aida.saveAs(AIDAFile);
+ } catch (IOException ex) {
+ Logger.getLogger(ElectronCorrection.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_conv;
+ if (Edep.isConverted()) {
+ Edep_conv = Edep.histogram();
+ } else {
+ System.out.println("Converting EDep");
+ double meanc = Edep.mean();
+ double rmsc = Edep.rms();
+ nsigmas = 3.;
+ nbins = 100;
+ rmsc = meanc * 0.001;
+ double minx = meanc - nsigmas * rmsc;
+ double maxx = meanc + nsigmas * rmsc;
+ Edep.setConversionParameters(nbins, minx, maxx);
+ Edep.convertToHistogram();
+ Edep_conv = Edep.histogram();
+ }
+ gauss.setParameter("amplitude", Edep_conv.maxBinHeight());
+ gauss.setParameter("mean", Edep_conv.mean());
+ gauss.setParameter("sigma", Edep_conv.rms());
+ try {
+ out.write(Particlename + " " + Ein_prev + " GeV\n");
+ out.write(Edep_conv.maxBinHeight() + "," + Edep_conv.mean() + "," + Edep_conv.rms() + "\n");
+ } catch (IOException ex) {
+ Logger.getLogger(ElectronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ 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++;
+ // chi 2 fit:
+ System.out.println("chi2 fit:");
+ jminuit = fitFactory.createFitter("Chi2", "jminuit");
+ jminuitResult = jminuit.fit(Edep_conv, gauss);
+ System.out.println(Edep_conv.maxBinHeight() + " , " + Edep_conv.mean() + " , " + Edep_conv.rms());
+ System.out.println("jminuit Chi2=" + jminuitResult.quality());
+ result = jminuitResult.fittedParameters();
+ functionFactory.cloneFunction("fitted gauss (jminuitchi2)", 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(ElectronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+
+ // least squares fit:
+ System.out.println("least squares fit:");
+ gauss.setParameter("amplitude", Edep_conv.maxBinHeight());
+ gauss.setParameter("mean", Edep_conv.mean());
+ gauss.setParameter("sigma", Edep_conv.rms());
+ IFitter jminuitls = fitFactory.createFitter("leastsquares", "jminuit");
+ IFitResult jminuitResultls = jminuitls.fit(Edep_conv, gauss);
+ System.out.println("jminuit ls=" + jminuitResultls.quality());
+ double[] resultls = jminuitResultls.fittedParameters();
+ functionFactory.cloneFunction("fitted gauss (jminuitls)", jminuitResultls.fittedFunction());
+ System.out.println(resultls[0] + "," + resultls[1] + "," + resultls[2]);
+ try {
+ out.write(resultls[0] + "," + resultls[1] + "," + resultls[2] + "\n");
+ } catch (IOException ex) {
+ Logger.getLogger(ElectronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ //
+ // 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();
+ nsigmas = 5.;
+ nbins = 100;
+ double minx = meanc - nsigmas * rmsc;
+ double maxx = meanc + nsigmas * rmsc;
+ Eceren.setConversionParameters(nbins, minx, maxx);
+ Eceren.convertToHistogram();
+ Eceren_conv = Eceren.histogram();
+ }
+ gauss.setParameter("amplitude", Eceren_conv.maxBinHeight());
+ gauss.setParameter("mean", Eceren_conv.mean());
+ gauss.setParameter("sigma", Eceren_conv.rms());
+ try {
+ out.write(Particlename + " " + Ein_prev + " GeV\n");
+ out.write(Eceren_conv.maxBinHeight() + "," + Eceren_conv.mean() + "," + Eceren_conv.rms() + "\n");
+ } catch (IOException ex) {
+ Logger.getLogger(ElectronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ 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++;
+ 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(ElectronCorrection.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(ElectronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+
+ public void setMyFilename(String filename) {
+ System.out.println("setMyFilename");
+ file_name = filename;
+ System.out.println(file_name);
+ }
+}