lcsim-contrib/src/main/java/org/lcsim/contrib/HansWenzel/DualCorrection
diff -u -r1.1.2.1 -r1.1.2.2
--- DualCorrection.java 25 Sep 2009 15:07:05 -0000 1.1.2.1
+++ DualCorrection.java 25 Sep 2009 19:33:38 -0000 1.1.2.2
@@ -7,6 +7,7 @@
*/
import hep.aida.ICloud1D;
import hep.aida.ICloud2D;
+import hep.aida.IDataPoint;
import hep.aida.IDataPointSet;
import hep.aida.IDataPointSetFactory;
import hep.aida.IFitFactory;
@@ -52,6 +53,7 @@
ICloud1D Eceren;
ICloud1D Edep_cor;
ICloud1D Edep_dcor;
+ ICloud1D Edep_dcor_new;
ICloud1D Eceren_cor;
ICloud2D c_efrac_ratio;
ICloud1D[] slice_comb;
@@ -62,20 +64,8 @@
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;
+ IDataPointSet[] dps_cerene;
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"};
@@ -84,6 +74,7 @@
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);
@@ -92,11 +83,20 @@
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());
+ for (int i = 0; i < Fitters.length; i++) {
+ String dpsname = "dps_cerene_" + Fitters[i];
+ point_cerene[i] = 0;
+ dps_cerene[i] = dpsf.create(dpsname, "ratio", 2);
+ }
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");
@@ -157,7 +157,13 @@
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");
+ for (int i = 0; i < Fitters.length; i++) {
+ String dpsname = "dps_cerene_" + Fitters[i];
+ point_cerene[i] = 0;
+ dps_cerene[i] = dpsf.create(dpsname, "ratio", 2);
+ }
slice = new ICloud1D[maxbin + 2];
conv_slice = new IHistogram1D[maxbin + 2];
for (int ibin = 0; ibin < maxbin + 2; ibin++) {
@@ -203,9 +209,10 @@
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;
@@ -224,6 +231,7 @@
convertandfit();
aida.tree().cd("/");
fitprofile(prof_combined);
+ convertandfit();
try {
out.close();
} catch (IOException ex) {
@@ -254,10 +262,10 @@
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));
+ /* 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();
@@ -283,7 +291,7 @@
} catch (IOException ex) {
Logger.getLogger(DualCorrection.class.getName()).log(Level.SEVERE, null, ex);
}
- for (int i = 0; i < slice.length; i++) {
+ for (int i = 0; i < slice.length - 3; i++) {
if (slice[i].isConverted()) {
conv_slice[i] = slice[i].histogram();
} else {
@@ -298,42 +306,54 @@
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);
+ int entries = conv_slice[i].entries();
+ if (entries > 100) {
+ for (int ii = 0; ii < Fitters.length; ii++) {
+ System.out.println("Fitter: " + Fitters[ii]);
+ gauss.setParameter("amplitude", conv_slice[i].maxBinHeight());
+ gauss.setParameter("mean", conv_slice[i].mean());
+ gauss.setParameter("sigma", conv_slice[i].rms());
+ jminuit = fitFactory.createFitter(Fitters[ii], "jminuit");
+ jminuitResult = jminuit.fit(conv_slice[i], gauss);
+ System.out.println("jminuit " + Fitters[ii] + ": " + jminuitResult.quality());
+ result = jminuitResult.fittedParameters();
+ errors = jminuitResult.errors();
+ String functionname = "ratio fitted gauss slice: " + 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]++;
+ }
+ }
}
- 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]++;
+ //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) {
@@ -392,4 +412,13 @@
ECor.setParameter("p3", -1.0873650332172629);
return ECor;
}
-}
\ No newline at end of file
+
+ 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;
+ }
+}