lcsim-contrib/src/main/java/org/lcsim/contrib/HansWenzel/DualCorrection
diff -u -r1.1 -r1.2
--- DRFunctionFactory.java 16 Dec 2009 17:20:26 -0000 1.1
+++ DRFunctionFactory.java 21 Dec 2009 22:45:04 -0000 1.2
@@ -44,12 +44,21 @@
// the parameters of the correction functions from a file.
//
DetectorConfiguration dualcor = new DetectorConfiguration("ccal02", "BGO", 7.13, 1.65, "Digis", 0.02, 0.02, "LCPhys");
- IFunction Corfu = functionFactory.createFunctionByName("dc_ccal02_digis_LCPhys", "p3");
- Corfu.setParameter("p0", 0.52415);
- Corfu.setParameter("p1", 0.37964);
- Corfu.setParameter("p2", -0.60463);
- Corfu.setParameter("p3", 0.66544);
+ //IFunction Corfu = functionFactory.createFunctionByName("dc_ccal02_digis_LCPhys", "p3");
+ //Corfu.setParameter("p0", 0.52415);
+ //Corfu.setParameter("p1", 0.37964);
+ //Corfu.setParameter("p2", -0.60463);
+ //Corfu.setParameter("p3", 0.66544);
+ //al.add(Corfu);
+ IFunction Corfu = functionFactory.createFunctionByName("dc_ccal02_digis_LCPhys", "p4");
+ Corfu.setParameter("p0", 0.51442);
+ Corfu.setParameter("p1", 0.45316);
+ Corfu.setParameter("p2", -0.60437);
+ Corfu.setParameter("p3", 0.53277);
+ Corfu.setParameter("p4", 0.077606);
al.add(Corfu);
+
+
IFunction Cerfu = functionFactory.createFunctionByName("cc_ccal02_digis_LCPhys", "p1");
Cerfu.setParameter("p0", 8.5701e-3);
Cerfu.setParameter("p1", 7677.7);
@@ -135,6 +144,7 @@
}
}
// get dual readout correction function
+
public IFunction getdcFunction(DetectorConfiguration dc) {
IFunction refunc = null;
if (mp.containsKey(dc)) {
@@ -144,6 +154,7 @@
return refunc;
}
// get cerenkov scale
+
public IFunction getccFunction(DetectorConfiguration dc) {
IFunction refunc = null;
if (mp.containsKey(dc)) {
@@ -153,6 +164,7 @@
return refunc;
}
// get electron scale:
+
public IFunction getecFunction(DetectorConfiguration dc) {
IFunction refunc = null;
if (mp.containsKey(dc)) {
@@ -162,6 +174,7 @@
return refunc;
}
/// check if there is an entry for a given DetectorConfiguration
+
boolean checkdc(DetectorConfiguration dc) {
if (mp.containsKey(dc)) {
return true;
lcsim-contrib/src/main/java/org/lcsim/contrib/HansWenzel/DualCorrection
diff -u -r1.9 -r1.10
--- DualCorrection.java 17 Dec 2009 21:02:12 -0000 1.9
+++ DualCorrection.java 21 Dec 2009 22:45:04 -0000 1.10
@@ -67,18 +67,18 @@
IProfile1D prof_combined;
double nsigmas;
int nbins;
- IDataPointSet[] dps_cerene;
- int point_cerene[];
+ IDataPointSet[] dps_ratio;
+ int point_ratio[];
double[] result;
double errors[];
-
String[] Fitters = {"Chi2", "leastsquares"};
- // Available are: "Chi2", "leastsquares", "bml", "cleverchi2";
+ // Available are: "Chi2", "leastsquares", "bml", "cleverchi2";
// but bml and cleverchi2 results don't make too much sense
IFitter jminuit;
IFitResult jminuitResult;
IFunction E_cor; // electron Ionisation correction function
IFunction C_cor; // electron cerenkov correction function
+ IFunction poly;
double xval[] = {10.};
double binwidth = 0.04;
int maxbin = (int) Math.rint(1.0 / binwidth);
@@ -98,11 +98,18 @@
protected void startOfData() {
System.out.println("DualCorrection:Start of Data");
fitFactory = aida.analysisFactory().createFitFactory();
- dps_cerene = new IDataPointSet[Fitters.length];
- point_cerene = new int[Fitters.length];
+ dps_ratio = new IDataPointSet[Fitters.length];
+ point_ratio = new int[Fitters.length];
dpsFactory = aida.analysisFactory().createDataPointSetFactory(aida.tree());
functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
gauss = functionFactory.createFunctionByName("gauss", "G");
+ poly = functionFactory.createFunctionByName("poly", "p4");
+ 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);
+ poly.setParameter("p4", 0.0);
DRFunctionFactory mapoff = DRFunctionFactory.getInstance();
DetectorConfiguration dc = new DetectorConfiguration(Detector, Material,
Density, rindex,
@@ -114,13 +121,12 @@
} else {
System.exit(0);
}
-// aida.profile1D("ratio", 20, 0.4, 1.);
dlength = aida.histogramFactory().createCloud1D("dlength", "decay length combined ", 500000, "autoConvert = false");
- prof_combined = aida.histogramFactory().createProfile1D("ratio_comb", "ratio combined", 20, 0.4, 1.);
+ prof_combined = aida.histogramFactory().createProfile1D("ratio_comb", "ratio combined", 30, 0.2, 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++) {
+ slice_comb = new ICloud1D[maxbin + 1];
+ conv_slice_comb = new IHistogram1D[maxbin + 1];
+ for (int ibin = 0; ibin < maxbin + 1; ibin++) {
String cloudy = "ratio_" + ibin;
slice_comb[ibin] = aida.histogramFactory().createCloud1D(cloudy, cloudy, 500000, "autoConvert = false");
}
@@ -142,7 +148,6 @@
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();
@@ -155,7 +160,6 @@
if (decaylength > 500.) {
return;
}
-
if (particle.getProductionTime() == 0.0) {
Particlename = particle.getType().getName();
}
@@ -165,12 +169,12 @@
Ein = (int) Math.floor(E_kin + 0.5d);
String E_str = Ein.toString();
DirName = Particlename.concat(E_str);
- if (Ein != Ein_prev||Particlename.equals(Part_prev)!=true) {
+ if (Ein != Ein_prev || Particlename.equals(Part_prev) != true) {
if (first) {
first = false;
} else {
System.out.println("Process event: ");
- fitprofile(aida.profile1D("ratio"));
+ fitprofile(aida.profile1D("ratio"), "dodo");
convertandfit(slice, conv_slice);
}
Ein_prev = Ein;
@@ -187,9 +191,9 @@
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 + 2];
- conv_slice = new IHistogram1D[maxbin + 2];
- for (int ibin = 0; ibin < maxbin + 2; ibin++) {
+ 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");
}
@@ -199,7 +203,7 @@
List<List<SimCalorimeterHit>> simCalorimeterHitCollections = event.get(SimCalorimeterHit.class);
double sumEEdep = 0.0;
double sumECeren = 0.0;
- aida.profile1D("ratio", 20, 0.4, 1.);
+ aida.profile1D("ratio", 30, 0.2, 1.);
for (List<SimCalorimeterHit> simCalorimeterHits : simCalorimeterHitCollections) {
String ColName = event.getMetaData(simCalorimeterHits).getName();
if (ColName.startsWith("Edep_") && ColName.endsWith("DigiHits")) {
@@ -236,8 +240,8 @@
if (bin < 0) {
bin = 0;
}
- if (bin > maxbin) {
- bin = maxbin + 1;
+ if (bin > maxbin - 1) {
+ bin = maxbin;
}
slice_comb[bin].fill(fraction);
slice[bin].fill(fraction);
@@ -245,11 +249,10 @@
@Override
protected void endOfData() {
- System.out.println("End of Data:");
-// fitprofile(aida.profile1D("ratio"));
+ System.out.println("DualCorrection: endOfData");
convertandfit(slice, conv_slice);
aida.tree().cd("/");
- fitprofile(prof_combined);
+ fitprofile(prof_combined, "profcombined");
convertandfit(slice_comb, conv_slice_comb);
try {
out.close();
@@ -270,14 +273,10 @@
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);
+ protected void fitprofile(IProfile1D profile, String Functionname) {
+
+ System.out.println("DualCorrection: fitprofile");
+
int profbins = profile.axis().bins();
double profxmin = profile.axis().lowerEdge();
double profxmax = profile.axis().upperEdge();
@@ -285,7 +284,7 @@
jminuitResult = jminuit.fit(profile, poly);
System.out.println("jminuit Chi2=" + jminuitResult.quality());
result = jminuitResult.fittedParameters();
- functionFactory.cloneFunction("fitted poly (jminuit)", jminuitResult.fittedFunction());
+ functionFactory.cloneFunction(Functionname, jminuitResult.fittedFunction());
System.out.println("correction function: " + result[0] + "," + result[1] + "," + result[2] + "," + result[3]);
try {
out.write(Particlename + " " + Ein_prev + " GeV\n");
@@ -293,15 +292,10 @@
} 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:");
+ System.out.println("DualCorrection: convert and fit:");
try {
out.write(Particlename + " " + Ein_prev + " GeV\n");
} catch (IOException ex) {
@@ -309,10 +303,10 @@
}
for (int i = 0; i < Fitters.length; i++) {
String dpsname = "dps_ratio_" + Fitters[i];
- point_cerene[i] = 0;
- dps_cerene[i] = dpsFactory.create(dpsname, "ratio", 2);
+ point_ratio[i] = 0;
+ dps_ratio[i] = dpsFactory.create(dpsname, "ratio", 2);
}
- for (int i = 0; i < slices.length - 3; i++) {
+ 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();
@@ -348,7 +342,7 @@
} catch (IOException ex) {
Logger.getLogger(DualCorrection.class.getName()).log(Level.SEVERE, null, ex);
}
- IDataPoint dp_cerene = dps_cerene[ii].addPoint();
+ IDataPoint dp_cerene = dps_ratio[ii].addPoint();
double ratio = (double) (i * binwidth) + 0.5 * binwidth;
dp_cerene.coordinate(0).setValue(ratio);
dp_cerene.coordinate(0).setErrorPlus(binwidth * 0.5);
@@ -356,30 +350,19 @@
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]++;
+ point_ratio[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);
+ jminuitResult = jminuit.fit(dps_ratio[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 setMyAIDAFilename(String jpgfile) {