lcsim-contrib/src/main/java/org/lcsim/contrib/HansWenzel/DualCorrection
diff -N Correction.java
--- Correction.java 27 Sep 2009 23:55:08 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,309 +0,0 @@
-package org.lcsim.contrib.HansWenzel.DualCorrection;
-
-/**
- *
- * @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 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 Correction 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;
- ICloud1D Edep_cor;
- ICloud1D Edep_dcor;
- ICloud1D Edep_dcor_new;
- ICloud1D Eceren_cor;
- double nsigmas;
- int nbins;
- IDataPointSet[] dps_Correctede = null;
- int point_Correctede[];
- int point;
- int point_Correctedemean;
- 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.};
-
- @Override
- protected void startOfData() {
- System.out.println("Start of Data:");
- dps_Correctede = new IDataPointSet[Fitters.length];
- point_Correctede = 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_Correctede_" + Fitters[i];
- point_Correctede[i] = 0;
- dps_Correctede[i] = dpsf.create(dpsname, "electron response mean", 2);
- }
- E_cor = Electron_Correction();
- C_cor = Cerenkov_Correction();
- D_cor = Dual_Correction();
- D_cor_new = Dual_Correction_new();
- Ein_prev = 0;
- firstEvent = true;
- first = true;
- // Create a two dimensional IDataPointSet.
- point = 0;
- point_Correctedemean = 0;
- gauss = functionFactory.createFunctionByName("gauss", "G");
- fstream = null;
- try {
- fstream = new FileWriter(file_name);
- } catch (IOException ex) {
- Logger.getLogger(Correction.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", "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_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");
- 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 Corrected
- } // 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;
- 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);
-
-
-
- }
-
- @Override
- protected void endOfData() {
- System.out.println("End of Data:");
- convertandfit();
-
- try {
- aida.saveAs(AIDAFile);
- } catch (IOException ex) {
- Logger.getLogger(Correction.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_new_conv;
-
- if (Edep_dcor_new.isConverted()) {
- Edep_dcor_new_conv = Edep_dcor_new.histogram();
- } else {
- System.out.println("Converting Edep_dcor_new");
-
- double meanc = Edep_dcor_new.mean();
- double rmsc = Edep_dcor_new.rms();
- nsigmas = 5.;
- nbins = 100;
- double minx = meanc - nsigmas * rmsc;
- double maxx = meanc + nsigmas * rmsc;
- Edep_dcor_new.setConversionParameters(nbins, minx, maxx);
- Edep_dcor_new.convertToHistogram();
- Edep_dcor_new_conv = Edep_dcor_new.histogram();
- }
- gauss.setParameter("amplitude", Edep_dcor_new_conv.maxBinHeight());
- gauss.setParameter("mean", Edep_dcor_new_conv.mean());
- gauss.setParameter("sigma", Edep_dcor_new_conv.rms());
- try {
- out.write(Particlename + " " + Ein_prev + " GeV\n");
- out.write(Edep_dcor_new_conv.maxBinHeight() + "," + Edep_dcor_new_conv.mean() + "," + Edep_dcor_new_conv.rms() + "\n");
- } catch (IOException ex) {
- Logger.getLogger(Correction.class.getName()).log(Level.SEVERE, null, ex);
- }
-
- for (int i = 0; i < Fitters.length; i++) {
- System.out.println("Fitter: " + Fitters[i]);
- gauss.setParameter("amplitude", Edep_dcor_new_conv.maxBinHeight());
- gauss.setParameter("mean", Edep_dcor_new_conv.mean());
- gauss.setParameter("sigma", Edep_dcor_new_conv.rms());
- jminuit = fitFactory.createFitter(Fitters[i], "jminuit");
- jminuitResult = jminuit.fit(Edep_dcor_new_conv, gauss);
- System.out.println("jminuit " + Fitters[i] + ": " + jminuitResult.quality());
- result = jminuitResult.fittedParameters();
- errors = jminuitResult.errors();
- String functionname = "Corrected 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(Correction.class.getName()).log(Level.SEVERE, null, ex);
- }
- 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]++;
- }
-
- try {
- aida.saveAs(AIDAFile);
- } catch (IOException ex) {
- Logger.getLogger(Correction.class.getName()).log(Level.SEVERE, null, ex);
- }
- }
-
- public void setMyFilename(String filename) {
- System.out.println("setMyFilename");
- file_name = filename;
- System.out.println(file_name);
- }
- public void setMyVariable(String jpgfile) {
- System.out.println("setMyVariable");
- AIDAFile = jpgfile;
- System.out.println(AIDAFile);
- }
- 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.65191);
- ECor.setParameter("p1", 0.37036);
- ECor.setParameter("p2", -0.69285);
- ECor.setParameter("p3", 0.62872);
- return ECor;
- }
-}