12 added + 3 modified, total 15 files
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
diff -N proton_resolution_steering.xml
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ proton_resolution_steering.xml 18 Apr 2014 20:52:20 -0000 1.1
@@ -0,0 +1,33 @@
+<lcsim xmlns:xs="http://www.w3.org/2001/XMLSchema-instance"
+ xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/lcsim/1.0/lcsim.xsd">
+
+ <classpath>
+ <jar>/home/aconway/NetBeansProjects/mcd-analysis/target/mcd-analysis-1.0-SNAPSHOT.jar</jar>
+ </classpath>
+
+ <inputFiles>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-1gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-2gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-5gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-10gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-25gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-50gev.slcio</file>
+ </inputFiles>
+ <control>
+ <numberOfEvents>-1</numberOfEvents>
+ </control>
+
+ <execute>
+ <driver name="DRResolutionDriver"/>
+ </execute>
+ <drivers>
+ <driver name="DRResolutionDriver"
+ type="org.lcsim.mcd.analysis.DRCorrection.alexDRCorr.DRResolutionDriver">
+ <myP_name>proton</myP_name>
+ <myS_thresh>0.0</myS_thresh>
+ <myC_thresh>0.0</myC_thresh>
+ <myT_thresh>3.0</myT_thresh>
+ <myTag>default</myTag>
+ </driver>
+ </drivers>
+</lcsim>
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
diff -N CorrectionFunctionSet.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ CorrectionFunctionSet.java 18 Apr 2014 20:52:20 -0000 1.1
@@ -0,0 +1,171 @@
+package org.lcsim.mcd.analysis.DRCorrection.alexDRCorr;
+
+import hep.aida.IFunction;
+import hep.aida.IFunctionFactory;
+import java.io.BufferedWriter;
+import java.io.File;
+import java.io.FileWriter;
+import java.io.IOException;
+import java.nio.charset.StandardCharsets;
+import java.nio.file.FileSystems;
+import java.nio.file.Files;
+import java.nio.file.Path;
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Scanner;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import org.lcsim.util.aida.AIDA;
+/**
+ * Read and write polynomial correction functions to a file,
+ * provide correction function set (HashMap<String,IFunction>) to drivers
+ * using a cutConfiguration object.
+ *
+ * Lines in the files have the following format:
+ * ,<fn name>,<parm name>,<parm value>,<parm name>,...,\n
+ * @author aconway
+ */
+public class CorrectionFunctionSet {
+
+ private AIDA aida;
+ private IFunctionFactory functionFactory;
+ private HashMap<String,IFunction> map;
+ private Path path;
+
+ public CorrectionFunctionSet(cutConfiguration cutco) {
+ aida = AIDA.defaultInstance();
+ functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+
+ String fname = cutco.toString() + ".txt";
+ path = FileSystems.getDefault().getPath("functions",fname);
+ String anyname = cutco.toStringAnyP() + ".txt";
+ Path anypath = FileSystems.getDefault().getPath("functions", anyname);
+ if (Files.isReadable(path)) {
+ map = loadFunctionsFromFile(path);
+ } else if (Files.isReadable(anypath)) {
+ map = loadFunctionsFromFile(anypath);
+ this.saveFunctionsToFile();
+ } else {
+ System.out.println("!!! Warning: Could not load functions from " + path.toString());
+ System.out.println("!!! Creating empty HashMap");
+ map = new HashMap<String,IFunction>();
+ this.saveFunctionsToFile();
+ }
+ }
+
+ public CorrectionFunctionSet(cutConfiguration cutco, HashMap<String,IFunction> inmap) {
+ aida = AIDA.defaultInstance();
+ functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+ String fname = cutco.toString() + ".txt";
+ path = FileSystems.getDefault().getPath("functions",fname);
+ map = inmap;
+ this.saveFunctionsToFile();
+ }
+
+ public CorrectionFunctionSet(cutConfiguration cutco, boolean DefaultEC) {
+ aida = AIDA.defaultInstance();
+ functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+ String fname = cutco.toString() + ".txt";
+ path = FileSystems.getDefault().getPath("functions",fname);
+ map = new HashMap<String,IFunction>();
+ if (DefaultEC) {
+ IFunction SC = functionFactory.createFunctionByName("SC", "p1");
+ IFunction CC = functionFactory.createFunctionByName("CC", "p1");
+ double[] defSCparms = {0.0,1.0};
+ double[] defCCparms = {3.6605e-3,1.87977e-5};
+ SC.setParameters(defSCparms);
+ CC.setParameters(defCCparms);
+ map.put("SC", SC);
+ map.put("CC", CC);
+ }
+ }
+
+ public void addFunction(String name, IFunction func) {
+ this.map.put(name, func);
+ this.saveFunctionsToFile();
+ }
+
+ public IFunction getFunction(String name) {
+ return this.map.get(name);
+ }
+
+ private HashMap loadFunctionsFromFile(Path path) {
+ HashMap rv = new HashMap<String,IFunction>();
+ try {
+ List<String> lines = Files.readAllLines(path, StandardCharsets.UTF_8);
+ for (String line : lines) {
+ Scanner sc = new Scanner(line);
+ sc.useDelimiter(",");
+ String type = null;
+ String name = sc.next();
+ List<Double> parmlist = new ArrayList<Double>();
+ List<String> namelist = new ArrayList<String>();
+ while (sc.hasNext()){
+ if (sc.hasNextDouble()) {
+ parmlist.add(sc.nextDouble());
+ } else {
+ namelist.add(sc.next());
+ }
+ }
+ if (name.matches("CS\\d.*")) {
+ type = "p4";
+ } else if (name.matches("[CS]C")) {
+ type = "p1";
+ } else if (name.matches("S\\d.*")) {
+ type = "G";
+ } else {
+ type = "p" + String.valueOf(parmlist.size()-1);
+ }
+ assert (parmlist.size() == namelist.size());
+ IFunction func = this.functionFactory.createFunctionByName(name, type);
+ for (int i=0; i<parmlist.size(); i++) {
+ func.setParameter(namelist.get(i), parmlist.get(i));
+ }
+ rv.put(name, func);
+ }
+ } catch (IOException ex) {
+ System.out.println("!!! ERROR: Could not read file at " + path.toString());
+ Logger.getLogger(CorrectionFunctionSet.class.getName()).log(Level.SEVERE, null, ex);
+ System.exit(0);
+ }
+ return rv;
+ }
+
+ private void saveFunctionsToFile() {
+ System.out.println("Saving functions to file " + this.path.toString());
+ try {
+ File ofile = this.path.toFile();
+ ofile.setWritable(true);
+ FileWriter fstream = new FileWriter(ofile, false);
+ BufferedWriter out = new BufferedWriter(fstream);
+ boolean first = true;
+ try {
+ for (String key : this.map.keySet()) {
+ System.out.println("\nWriting function " + key);
+ if (first) {
+ first = false;
+ } else {
+ out.newLine();
+ }
+ IFunction func = this.map.get(key);
+ String[] pnames = func.parameterNames();
+ double[] params = func.parameters();
+ out.write("," + key + ",");
+ for (int i=0; i<params.length; i++) {
+ System.out.print(pnames[i] + ", " + String.valueOf(params[i]) + ", ");
+ out.write(pnames[i] + "," + String.valueOf(params[i]) + ",");
+ }
+ }
+ out.close();
+ System.out.println("\nSuccessfully saved corrections to " +this.path.toString());
+ } catch (IOException ex) {
+ System.out.println("!!! ERROR: Error writing to file " + this.path.toString());
+ Logger.getLogger(CorrectionFunctionSet.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ } catch (IOException ex) {
+ System.out.println("!!! ERROR: Could not open file for writing " + this.path.toString());
+ Logger.getLogger(CorrectionFunctionSet.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
diff -N correctionFunctions.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ correctionFunctions.java 18 Apr 2014 20:52:20 -0000 1.1
@@ -0,0 +1,65 @@
+package org.lcsim.mcd.analysis.DRCorrection.alexDRCorr;
+
+import hep.aida.IFunction;
+import hep.aida.IFunctionFactory;
+import java.util.ArrayList;
+import java.util.Dictionary;
+import java.util.HashMap;
+import java.util.Map;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author aconway
+ */
+public class correctionFunctions {
+
+ private AIDA aida;
+ private IFunctionFactory functionFactory;
+ private static Map<cutConfiguration, Map<String,IFunction>> map;
+
+ public correctionFunctions() {
+ aida = AIDA.defaultInstance();
+ functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+ map = new HashMap<cutConfiguration, Map<String,IFunction>>();
+ Map submap = new HashMap<String,IFunction>();
+
+ cutConfiguration cutcon = new cutConfiguration(0.0,0.0,0.0,"any","default");
+ double[] defSCparms = {0.0,1.0};
+ double[] defCCparms = {3.6605e-3,1.87977e-5};
+ IFunction SC = functionFactory.createFunctionByName("SC", "p1");
+ IFunction CC = functionFactory.createFunctionByName("CC", "p1");
+ SC.setParameters(defSCparms);
+ CC.setParameters(defCCparms);
+ submap.put("SC", SC);
+ submap.put("CC", CC);
+ double[] def1gparms = {0.0,0.0,0.0,0.0,0.0};
+ double[] def2gparms = {0.0,0.0,0.0,0.0,0.0};
+ double[] def5gparms = {0.0,0.0,0.0,0.0,0.0};
+ double[] def10gparms = {0.0,0.0,0.0,0.0,0.0};
+ double[] def25gparms = {0.0,0.0,0.0,0.0,0.0};
+ double[] def50gparms = {0.0,0.0,0.0,0.0,0.0};
+ double[] def100gparms = {0.0,0.0,0.0,0.0,0.0};
+ IFunction cs1gev = functionFactory.createFunctionByName("1gev", "p4");
+ IFunction cs2gev = functionFactory.createFunctionByName("2gev", "p4");
+ IFunction cs5gev = functionFactory.createFunctionByName("5gev", "p4");
+ IFunction cs10gev = functionFactory.createFunctionByName("10gev", "p4");
+ IFunction cs25gev = functionFactory.createFunctionByName("25gev", "p4");
+ IFunction cs50gev = functionFactory.createFunctionByName("50gev", "p4");
+ IFunction cs100gev = functionFactory.createFunctionByName("100gev", "p4");
+ cs1gev.setParameters(def1gparms);
+ cs2gev.setParameters(def2gparms);
+ cs5gev.setParameters(def5gparms);
+ cs10gev.setParameters(def10gparms);
+ cs25gev.setParameters(def25gparms);
+ cs50gev.setParameters(def50gparms);
+ cs100gev.setParameters(def100gparms);
+ submap.put("cs1gev", cs1gev);
+ submap.put("cs2gev", cs2gev);
+ submap.put("cs5gev", cs5gev);
+ submap.put("cs10gev", cs10gev);
+ submap.put("cs25gev", cs25gev);
+ submap.put("cs50gev", cs50gev);
+ submap.put("cs100gev", cs100gev);
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
diff -N electronCorrection.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ electronCorrection.java 18 Apr 2014 20:52:20 -0000 1.1
@@ -0,0 +1,326 @@
+package org.lcsim.mcd.analysis.DRCorrection.alexDRCorr;
+
+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 hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import java.io.BufferedWriter;
+import java.io.FileWriter;
+import java.io.IOException;
+import java.util.List;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author aconway
+ */
+public class electronCorrection extends Driver {
+ String AIDAFile;
+ AIDA aida;
+ IFunctionFactory functionFactory;
+ IFitFactory fitFactory;
+ IDataPointSetFactory dpsFactory;
+
+ double E_in;
+ double E_kin;
+ String Particlename;
+ Integer Ein;
+ Integer Ein_prev;
+ boolean firstEvent;
+ boolean first;
+
+ ICloud1D Edep;
+ ICloud1D Eceren;
+ IFunction gauss;
+ IDataPointSet dps_emean = null;
+ IDataPointSet dps_cerenemean = null;
+ IDataPointSet dps_cerene = null;
+ IFitter jminuit;
+ IFitResult jminuitResult;
+ int point;
+ int point_cerenemean;
+ int point_cerene;
+
+ cutConfiguration cutCon;
+ CorrectionFunctionSet corfus;
+
+ public electronCorrection() {
+ Ein_prev = 0;
+ firstEvent = true;
+ first = true;
+
+ aida = AIDA.defaultInstance();
+ fitFactory = aida.analysisFactory().createFitFactory();
+ functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+ dpsFactory = aida.analysisFactory().createDataPointSetFactory(aida.tree());
+ gauss = functionFactory.createFunctionByName("gauss", "G");
+ dps_emean = dpsFactory.create("dps_emean", "electron response mean", 2);
+ dps_cerenemean = dpsFactory.create("dps_cerenemean", "electron response mean", 2);
+ dps_cerene = dpsFactory.create("dps_cerene_Chi2", "electron response mean", 2);
+ point = 0;
+ point_cerenemean = 0;
+ }
+
+ @Override
+ protected void startOfData() {
+ }
+
+ @Override
+ protected void process(EventHeader event) {
+ E_in = 0.0;
+ E_kin = 0.0;
+ double Chrg = 0.0;
+ double P_in = 0.0;
+ Hep3Vector origin = new BasicHep3Vector(0,0,0);
+ String DirName = null;
+ List<MCParticle> particles = event.get(MCParticle.class, event.MC_PARTICLES);
+ for (MCParticle particle : particles) {
+ if (particle.getParents().isEmpty()) {
+ E_in = E_in + particle.getEnergy();
+ E_kin = E_kin + particle.getEnergy() - particle.getMass();
+ P_in = particle.getMomentum().magnitude();
+ Particlename = particle.getType().getName();
+ origin = particle.getOrigin();
+ Chrg = particle.getCharge();
+ }
+ break;
+ }
+ if (E_kin < 0.999) {
+ Ein = (int) Math.floor(1000.0*E_kin + 0.5d);
+ } else {
+ Ein = (int) Math.floor(E_kin + 0.5d);
+ }
+ String E_str = Ein.toString();
+ DirName = Particlename.concat(E_str);
+ if (!Ein.equals(Ein_prev)) {
+ if (first) {
+ first = false;
+ } else {
+ System.out.println("electronCorrection:First Event");
+ 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);
+ }
+ double dif_curv = 0.0;
+ if (Chrg != 0.0) {
+ double rad_curv = P_in * 667.128;
+ double len_curv = rad_curv * 2.0 * Math.asin(0.667 / rad_curv);
+ dif_curv = len_curv - 1.334;
+ }
+ 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.contains("Edep")) {
+ for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+ if (cutCon.T_thresh > 0.0 || cutCon.S_thresh > 0.0) {
+ Hep3Vector diffvec = VecOp.sub(calorimeterHit.getPositionVec(), origin);
+ int ncontribs = calorimeterHit.getMCParticleCount();
+ for (int i=0; i<ncontribs; i++) {
+ double tdep = calorimeterHit.getContributedTime(i);
+ double tofl = tdep - (diffvec.magnitude()+dif_curv)/299.79;
+ double edep = calorimeterHit.getContributedEnergy(i);
+ if ( (tofl < cutCon.T_thresh) && (edep > cutCon.S_thresh/1000.0) ) {
+ sumEEdep += edep;
+ }
+ }
+ }
+ else {
+ sumEEdep = sumEEdep + calorimeterHit.getRawEnergy();
+ }
+ }
+ } // end if Edep
+ if (CollectionName.contains("Opti")) {
+ for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+ if (cutCon.T_thresh > 0.0 || cutCon.C_thresh > 0.0) {
+ Hep3Vector diffvec = VecOp.sub(calorimeterHit.getPositionVec(), origin);
+ int ncontribs = calorimeterHit.getMCParticleCount();
+ for (int i=0; i<ncontribs; i++) {
+ double tdep = calorimeterHit.getContributedTime(i);
+ double tofl = tdep - (diffvec.magnitude()+dif_curv)/299.79;
+ double edep = calorimeterHit.getContributedEnergy(i);
+ if ( (tofl < cutCon.T_thresh) && (edep > cutCon.C_thresh) ) {
+ sumECeren += edep;
+ }
+ }
+ }
+ else {
+ sumECeren = sumECeren + calorimeterHit.getRawEnergy();
+ }
+ }
+ }
+ }
+ Edep.fill(sumEEdep);
+ Eceren.fill(sumECeren);
+ }
+
+ @Override
+ protected void endOfData() {
+ System.out.println("electronCorrection:End of Data:");
+ convertandfit();
+ aida.tree().cd("/");
+ 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());
+ corfus.addFunction("SC",resultline.fittedFunction());
+
+ double[] fresult = resultline.fittedParameters();
+ System.out.println("Ionization scale results:");
+ System.out.println("Chi2 = " + resultline.quality());
+ System.out.println(fresult[0] + " , " + fresult[1]);
+ //
+ // now deal with the cerenkov stuff:
+ //
+ String fitter = "Chi2";
+ System.out.println("Fitter: " + fitter);
+ resultline = linefit.fit(dps_cerene, line);
+ String fname = "cerenkov " + fitter + " fitted line";
+ System.out.println(fname);
+ functionFactory.cloneFunction(fname, resultline.fittedFunction());
+ corfus.addFunction("CC",resultline.fittedFunction());
+
+ System.out.println(fitter + ": " + resultline.quality());
+ fresult = resultline.fittedParameters();
+ System.out.println("Cerenkov scale results:");
+ System.out.println("Fitter: " + fitter);
+ System.out.println("Chi2 = " + resultline.quality());
+ System.out.println(fresult[0] + " , " + fresult[1]);
+ try {
+ aida.saveAs(AIDAFile);
+ } catch (IOException ex) {
+ Logger.getLogger(electronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+
+ protected void convertandfit() {
+ System.out.println("electronCorrection:convertandfit");
+ IHistogram1D Edep_conv;
+ if (Edep.isConverted()) {
+ Edep_conv = Edep.histogram();
+ } else {
+ System.out.println("Converting EDep");
+ double meanc = Edep.mean();
+ double rmsc = Edep.rms();
+ double nsigmas = 3.;
+ int nbins = 100;
+
+ double minx = meanc - nsigmas * rmsc;
+ double maxx = meanc + nsigmas * rmsc;
+ Edep.setConversionParameters(nbins, minx, maxx);
+ Edep.convertToHistogram();
+ Edep_conv = Edep.histogram();
+ }
+
+ 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++;
+
+ //
+ // 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();
+ double nsigmas = 5.;
+ int nbins = 100;
+ double minx = meanc - nsigmas * rmsc;
+ double maxx = meanc + nsigmas * rmsc;
+ Eceren.setConversionParameters(nbins, minx, maxx);
+ Eceren.convertToHistogram();
+ Eceren_conv = Eceren.histogram();
+ }
+
+ 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++;
+
+ String fitter = "Chi2";
+ System.out.println("Fitter: " + fitter);
+ gauss.setParameter("amplitude", Eceren_conv.maxBinHeight());
+ gauss.setParameter("mean", Eceren_conv.mean());
+ gauss.setParameter("sigma", Eceren_conv.rms());
+ jminuit = fitFactory.createFitter(fitter, "jminuit");
+ jminuitResult = jminuit.fit(Eceren_conv, gauss);
+ System.out.println("jminuit " + fitter + ": " + jminuitResult.quality());
+ double[] result = jminuitResult.fittedParameters();
+ double[] errors = jminuitResult.errors();
+ String functionname = "ceren fitted gauss " + fitter;
+ System.out.println(functionname);
+ functionFactory.cloneFunction(functionname, jminuitResult.fittedFunction());
+ System.out.println(result[0] + "," + result[1] + "," + result[2]);
+ IDataPoint dp_cerene = dps_cerene.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++;
+
+ try {
+ aida.saveAs(AIDAFile);
+ } catch (IOException ex) {
+ Logger.getLogger(electronCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+
+ public void setCorrFu(CorrectionFunctionSet corrs) {
+ System.out.println("electronCorrection:setCorrFu");
+ this.corfus = corrs;
+ }
+
+ public void setCutConfiguration(cutConfiguration cutco) {
+ System.out.println("electronCorrection:setCutConfiguration");
+ this.cutCon = cutco;
+ }
+
+ public void setMyAIDAFilename(String AIDAFile) {
+ System.out.println("electronCorrection:setMyAIDAFilename");
+ this.AIDAFile = AIDAFile;
+ System.out.println(AIDAFile);
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
diff -N ElectronCalibrationDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ElectronCalibrationDriver.java 18 Apr 2014 20:52:20 -0000 1.1
@@ -0,0 +1,64 @@
+package org.lcsim.mcd.analysis.DRCorrection.alexDRCorr;
+
+import org.lcsim.util.Driver;
+
+/**
+ *
+ * @author aconway
+ */
+public class ElectronCalibrationDriver extends Driver {
+
+ private String P_name = "any";
+ private String tag = "default";
+ private Double C_thresh = 0.0;
+ private Double S_thresh = 0.0;
+ private Double T_thresh = 0.0;
+ private String AIDAFile = null;
+
+ electronCorrection elec;
+
+ public ElectronCalibrationDriver() {
+ elec = new electronCorrection();
+ add(elec);
+ }
+
+ @Override
+ public void startOfData() {
+ System.out.println("ElectronCalibrationDriver:startOfData");
+ cutConfiguration cutco;
+ if (tag.equals("default")) {
+ cutco = new cutConfiguration(C_thresh, S_thresh, T_thresh, "any");
+ } else {
+ cutco = new cutConfiguration(C_thresh, S_thresh, T_thresh, "any", tag);
+ }
+ AIDAFile = "EC-" + cutco.toString() + ".aida";
+ elec.setMyAIDAFilename(AIDAFile);
+ elec.setCutConfiguration(cutco);
+ elec.setCorrFu(new CorrectionFunctionSet(cutco, false));
+ }
+
+ public void setMyP_name(String pname) {
+ System.out.println("ElectronCalibrationDriver:setMyP_name = " + pname);
+ this.P_name = pname;
+ }
+
+ public void setMyTag(String tag) {
+ System.out.println("ElectronCalibrationDriver:setMyTag = " + tag);
+ this.tag = tag;
+ }
+
+ public void setMyC_thresh(double cthresh) {
+ System.out.println("ElectronCalibrationDriver:setMyC_thresh = " + cthresh);
+ this.C_thresh = cthresh;
+ }
+
+ public void setMyS_thresh(double sthresh) {
+ System.out.println("ElectronCalibrationDriver:setMyS_thresh = " + sthresh);
+ this.S_thresh = sthresh;
+ }
+
+ public void setMyT_thresh(double tthresh) {
+ System.out.println("ElectronCalibrationDriver:setMyT_thresh = " + tthresh);
+ this.T_thresh = tthresh;
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
diff -N DRResolutionDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ DRResolutionDriver.java 18 Apr 2014 20:52:20 -0000 1.1
@@ -0,0 +1,65 @@
+package org.lcsim.mcd.analysis.DRCorrection.alexDRCorr;
+
+import org.lcsim.util.Driver;
+
+/**
+ *
+ * @author aconway
+ */
+public class DRResolutionDriver extends Driver {
+
+ private String P_name = "any";
+ private String tag = "default";
+ private Double C_thresh = 0.0;
+ private Double S_thresh = 0.0;
+ private Double T_thresh = 0.0;
+ private String AIDAFile = null;
+
+ Resolution resol;
+
+ public DRResolutionDriver() {
+ resol = new Resolution();
+ add(resol);
+ }
+
+ @Override
+ public void startOfData() {
+ System.out.println("DRResolutionDriver:startOfData");
+ cutConfiguration cutco;
+ if (tag.equals("default")) {
+ cutco = new cutConfiguration(C_thresh, S_thresh, T_thresh, P_name);
+ } else {
+ cutco = new cutConfiguration(C_thresh, S_thresh, T_thresh, P_name, tag);
+ }
+ AIDAFile = "res-" + cutco.toString() + ".aida";
+ resol.setMyAIDAFilename(AIDAFile);
+ resol.setCutConfiguration(cutco);
+ resol.setCorrFu(new CorrectionFunctionSet(cutco));
+ resol.startOfData();
+ }
+
+ public void setMyP_name(String pname) {
+ System.out.println("DRResolutionDriver:setMyP_name = " + pname);
+ this.P_name = pname;
+ }
+
+ public void setMyTag(String tag) {
+ System.out.println("DRResolutionDriver:setMyTag = " + tag);
+ this.tag = tag;
+ }
+
+ public void setMyC_thresh(double cthresh) {
+ System.out.println("DRResolutionDriver:setMyC_thresh = " + cthresh);
+ this.C_thresh = cthresh;
+ }
+
+ public void setMyS_thresh(double sthresh) {
+ System.out.println("DRResolutionDriver:setMyS_thresh = " + sthresh);
+ this.S_thresh = sthresh;
+ }
+
+ public void setMyT_thresh(double tthresh) {
+ System.out.println("DRResolutionDriver:setMyT_thresh = " + tthresh);
+ this.T_thresh = tthresh;
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
diff -N DualCorrection.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ DualCorrection.java 18 Apr 2014 20:52:20 -0000 1.1
@@ -0,0 +1,353 @@
+package org.lcsim.mcd.analysis.DRCorrection.alexDRCorr;
+
+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.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import java.io.IOException;
+import java.util.List;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author aconway
+ */
+public class DualCorrection extends Driver {
+
+ String AIDAFile;
+ AIDA aida = AIDA.defaultInstance();
+ IFunctionFactory functionFactory;
+ IFitFactory fitFactory;
+ IDataPointSetFactory dpsFactory;
+
+ ICloud1D Edep;
+ ICloud1D Eceren;
+ ICloud1D Edep_frac;
+ ICloud1D Edep_cor;
+ ICloud1D Eceren_cor;
+ ICloud2D c_efrac_ratio;
+ ICloud2D c_Ceren_vs_Edep;
+ ICloud1D[] slice;
+ IHistogram1D[] conv_slice;
+ IDataPointSet dps_ratio;
+ IFunction gauss;
+ IFunction poly;
+ IFitter jminuit;
+ IFitResult jminuitResult;
+ double binwidth = 0.05;
+ int maxbin = (int) Math.rint(1.5 / binwidth);
+
+ double E_in;
+ double E_kin;
+ String Particlename;
+ Integer Ein;
+ Integer Ein_prev;
+ boolean firstEvent;
+ boolean first;
+
+ cutConfiguration cutCon;
+ CorrectionFunctionSet corfus;
+
+ public DualCorrection() {
+ Ein_prev = 0;
+ firstEvent = true;
+ first = true;
+
+ fitFactory = aida.analysisFactory().createFitFactory();
+ functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+ dpsFactory = aida.analysisFactory().createDataPointSetFactory(aida.tree());
+
+ gauss = functionFactory.createFunctionByName("gauss", "G");
+ jminuit = fitFactory.createFitter("Chi2", "jminuit");
+ poly = functionFactory.createFunctionByName("poly", "p4");
+ 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);
+ }
+
+ // initialize variables, etc.
+ @Override
+ protected void startOfData() {
+
+ }
+
+ // Sum up Edep and Eceren, run fits
+ @Override
+ protected void process(EventHeader event) {
+ E_in = 0.0;
+ E_kin = 0.0;
+ double Chrg = 0.0;
+ double P_in = 0.0;
+ Hep3Vector origin = new BasicHep3Vector(0,0,0);
+ String DirName = null;
+ List<MCParticle> particles = event.get(MCParticle.class, event.MC_PARTICLES);
+ for (MCParticle particle : particles) {
+ if (particle.getParents().isEmpty()) {
+ E_in = E_in + particle.getEnergy();
+ E_kin = E_kin + particle.getEnergy() - particle.getMass();
+ P_in = particle.getMomentum().magnitude();
+ Particlename = particle.getType().getName();
+ origin = particle.getOrigin();
+ Chrg = particle.getCharge();
+ }
+ break;
+ }
+ if (E_kin < 0.999) {
+ Ein = (int) Math.floor(1000.0*E_kin + 0.5d);
+ } else {
+ Ein = (int) Math.floor(E_kin + 0.5d);
+ }
+ String E_str = Ein.toString();
+ DirName = Particlename.concat(E_str);
+ if (!Ein.equals(Ein_prev)) {
+ if (first) {
+ first = false;
+ } else {
+ System.out.println("Process event: ");
+ 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");
+ Edep_frac = aida.histogramFactory().createCloud1D("Edep frac", "Energy frac 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");
+ 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 + 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");
+ }
+ System.out.println("DirName: " + DirName);
+ }
+ double dif_curv = 0.0;
+ if (Chrg != 0.0) {
+ double rad_curv = P_in * 667.128;
+ double len_curv = rad_curv * 2.0 * Math.asin(0.667 / rad_curv);
+ dif_curv = len_curv - 1.334;
+ }
+ 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.contains("Edep")) {
+ for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+ if (cutCon.T_thresh > 0.0 || cutCon.S_thresh > 0.0) {
+ Hep3Vector diffvec = VecOp.sub(calorimeterHit.getPositionVec(), origin);
+ int ncontribs = calorimeterHit.getMCParticleCount();
+ for (int i=0; i<ncontribs; i++) {
+ double tdep = calorimeterHit.getContributedTime(i);
+ double tofl = tdep - (diffvec.magnitude()+dif_curv)/299.79;
+ double edep = calorimeterHit.getContributedEnergy(i);
+ if ( (tofl < cutCon.T_thresh) && (edep > cutCon.S_thresh/1000.0) ) {
+ sumEEdep += edep;
+ }
+ }
+ }
+ else {
+ sumEEdep = sumEEdep + calorimeterHit.getRawEnergy();
+ }
+ }
+ } // end if Edep
+ if (CollectionName.contains("Opti")) {
+ for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+ if (cutCon.T_thresh > 0.0 || cutCon.C_thresh > 0.0) {
+ Hep3Vector diffvec = VecOp.sub(calorimeterHit.getPositionVec(), origin);
+ int ncontribs = calorimeterHit.getMCParticleCount();
+ for (int i=0; i<ncontribs; i++) {
+ double tdep = calorimeterHit.getContributedTime(i);
+ double tofl = tdep - (diffvec.magnitude()+dif_curv)/299.79;
+ double edep = calorimeterHit.getContributedEnergy(i);
+ if ( (tofl < cutCon.T_thresh) && (edep > cutCon.C_thresh) ) {
+ sumECeren += edep;
+ }
+ }
+ }
+ else {
+ sumECeren = sumECeren + calorimeterHit.getRawEnergy();
+ }
+ }
+ }
+ }
+ Edep.fill(sumEEdep);
+ Eceren.fill(sumECeren);
+ Edep_frac.fill(sumEEdep/E_kin);
+
+ double[] xval = {0.0};
+ xval[0] = sumEEdep;
+ double sumEEdep_cor = corfus.getFunction("SC").value(xval);
+ xval[0] = sumECeren;
+ double sumECeren_cor = corfus.getFunction("CC").value(xval);
+ Edep_cor.fill(sumEEdep_cor);
+ Eceren_cor.fill(sumECeren_cor);
+ double ratio = sumECeren_cor / sumEEdep_cor;
+ double fraction = sumEEdep_cor / E_kin;
+ aida.cloud1D("frac").fill(fraction);
+ aida.cloud1D("c_CerenEdep_ratio").fill(ratio);
+ c_Ceren_vs_Edep.fill(sumECeren_cor, sumEEdep_cor);
+ aida.cloud2D("c_efrac_ratio").fill(ratio, fraction);
+
+ int bin = (int) Math.rint(ratio / binwidth);
+ if (bin < 0) {
+ bin = 0;
+ }
+ if (bin > maxbin - 1) {
+ bin = maxbin;
+ }
+ slice[bin].fill(fraction);
+ }
+
+ // Run final fits
+ @Override
+ protected void endOfData() {
+ System.out.println("DualCorrection: endOfData");
+ convertandfit(slice, conv_slice);
+ }
+
+ protected void convertandfit(ICloud1D[] slices, IHistogram1D[] conv_slices) {
+ System.out.println("DualCorrection: convert and fit:");
+
+ String fitter = "Chi2";
+ String dpsname = "dps_ratio_" + fitter;
+ dps_ratio = dpsFactory.create(dpsname, "ratio", 2);
+
+ 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();
+ } else {
+ double meanc = slices[i].mean();
+ double rmsc = slices[i].rms();
+ double nsigmas = 3.;
+ int 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 > 300) {
+ System.out.println("Fitter: " + fitter);
+ double meanguess = 0;
+ for (int j = 0; j < conv_slices[i].axis().bins(); j++) {
+ if (conv_slices[i].binHeight(j) == conv_slices[i].maxBinHeight()) {
+ meanguess = conv_slices[i].axis().binCenter(j);
+ }
+ }
+ gauss.setParameter("amplitude", conv_slices[i].maxBinHeight());
+ gauss.setParameter("mean", meanguess);
+ gauss.setParameter("sigma", conv_slices[i].rms() / 2.0);
+ jminuit = fitFactory.createFitter(fitter, "jminuit");
+ try {
+ jminuitResult = jminuit.fit(conv_slices[i], gauss);
+ } catch (RuntimeException ex) {
+ System.out.println("Runtime Exception caught in fit!");
+ System.out.println(ex.toString());
+ } catch (Exception ex) {
+ System.out.println("Exception caught in fit!");
+ System.out.println(ex.toString());
+ }
+ System.out.println("jminuit " + fitter + ": " + jminuitResult.quality());
+ double[] result = jminuitResult.fittedParameters();
+ String functionname = "ratio fitted gauss slices: " + i + " " + fitter;
+ System.out.println(functionname);
+ functionFactory.cloneFunction(functionname, jminuitResult.fittedFunction());
+ System.out.println(result[0] + "," + result[1] + "," + result[2]);
+ IDataPoint dp_cerene = dps_ratio.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]));
+ }
+ }
+ //Now fit the data point sets:
+ System.out.println("Fitter: " + fitter);
+ jminuitResult = jminuit.fit(dps_ratio, poly);
+ String fname = "correction " + fitter + " fitted poly";
+ System.out.println(fname);
+ functionFactory.cloneFunction(fname, jminuitResult.fittedFunction());
+ System.out.println(fitter + ": " + jminuitResult.quality());
+ corfus.addFunction("CS" + Ein_prev.toString() + "gev", jminuitResult.fittedFunction());
+
+ // Now get Edep-only correction
+ double meanc = Edep_frac.mean();
+ double rmsc = Edep_frac.rms();
+ double nsigmas = 3.;
+ int nbins = 100;
+ double minx = meanc - nsigmas * rmsc;
+ double maxx = meanc + nsigmas * rmsc;
+ Edep_frac.setConversionParameters(nbins, minx, maxx);
+ Edep_frac.convertToHistogram();
+ IHistogram1D Edep_conv = Edep_frac.histogram();
+ gauss.setParameter("amplitude", Edep_conv.maxBinHeight());
+ gauss.setParameter("mean", Edep_conv.mean());
+ gauss.setParameter("sigma", Edep_conv.rms());
+ jminuit = fitFactory.createFitter(fitter, "jminuit");
+ try {
+ jminuitResult = jminuit.fit(Edep_conv, gauss);
+ } catch (RuntimeException ex) {
+ System.out.println("Runtime Exception caught in fit!");
+ System.out.println(ex.toString());
+ } catch (Exception ex) {
+ System.out.println("Exception caught in fit!");
+ System.out.println(ex.toString());
+ }
+ fname = "S" + Ein_prev.toString() + "gev";
+ functionFactory.cloneFunction("fname", jminuitResult.fittedFunction());
+ corfus.addFunction(fname, jminuitResult.fittedFunction());
+
+ try {
+ aida.saveAs(AIDAFile);
+ System.out.println("AIDAFile Saved.");
+ } catch (IOException ex) {
+ Logger.getLogger(org.lcsim.mcd.analysis.DRCorrection.DualCorrection.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+
+ public void setCorrFu(CorrectionFunctionSet corrs) {
+ System.out.println("DualCorrection:setCorrFu");
+ this.corfus = corrs;
+ }
+
+ public void setCutConfiguration(cutConfiguration cutco) {
+ System.out.println("DualCorrection:setCutConfiguration");
+ this.cutCon = cutco;
+ }
+
+ public void setMyAIDAFilename(String AIDAFile) {
+ System.out.println("DualCorrection:setMyAIDAFilename");
+ this.AIDAFile = AIDAFile;
+ System.out.println(AIDAFile);
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
diff -N Resolution.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ Resolution.java 18 Apr 2014 20:52:20 -0000 1.1
@@ -0,0 +1,359 @@
+package org.lcsim.mcd.analysis.DRCorrection.alexDRCorr;
+
+/**
+ *
+ * @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 hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+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 Resolution extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+ String AIDAFile;
+ String file_name;
+ IFunctionFactory functionFactory;
+ IFitFactory fitFactory;
+ IDataPointSetFactory dpsFactory;
+ IFunction gauss;
+ static boolean first;
+ static boolean firstEvent;
+ static 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_scor;
+ ICloud1D Eceren_cor;
+ double nsigmas=5.0;
+ int nbins=100;
+ IDataPointSet[] dps_Correctede = null;
+ IDataPointSet[] dps_DRResolution = null;
+ IDataPointSet[] dps_SRResolution = null;
+ IDataPointSet dps_ECorrFrac = null;
+ IDataPointSet dps_CCorrFrac = null;
+ IDataPointSet[] dps_DCorrFrac = null;
+ int point_Correctede[];
+ double[] result;
+ double errors[];
+// String[] Fitters = {"Chi2", "leastsquares"};
+ String[] Fitters = {"Chi2"};
+ // available are "Chi2", "leastsquares", "bml", "cleverchi2"
+ // but "bml", "cleverchi2" seem to produce rubish
+
+ IFitter jminuit;
+ IFitResult jminuitResult;
+ IFunction E_cor;
+ IFunction C_cor;
+ IFunction D_cor;
+ double xval[] = {10.};
+
+
+ cutConfiguration cutCon;
+ CorrectionFunctionSet corfus;
+
+ @Override
+ protected void startOfData() {
+ System.out.println("Start of Data:");
+ dps_Correctede = new IDataPointSet[Fitters.length];
+ dps_DRResolution = new IDataPointSet[Fitters.length];
+ dps_SRResolution = new IDataPointSet[Fitters.length];
+ dps_DCorrFrac = new IDataPointSet[Fitters.length];
+ point_Correctede = new int[Fitters.length];
+ fitFactory = aida.analysisFactory().createFitFactory();
+ functionFactory = aida.analysisFactory().createFunctionFactory(aida.tree());
+ dpsFactory = aida.analysisFactory().createDataPointSetFactory(aida.tree());
+ dps_ECorrFrac = dpsFactory.create("dps_ECorrFrac_", 2);
+ dps_CCorrFrac = dpsFactory.create("dps_CCorrFrac_", 2);
+ for (int i = 0; i < Fitters.length; i++) {
+ String dpsname = "dps_Correctede_" + Fitters[i];
+ point_Correctede[i] = 0;
+ dps_Correctede[i] = dpsFactory.create(dpsname, "electron response mean", 2);
+ dps_DRResolution[i] = dpsFactory.create("dps_DRResolution_"+Fitters[i], 2);
+ dps_SRResolution[i] = dpsFactory.create("dps_SRResolution_"+Fitters[i], 2);
+ dps_DCorrFrac[i] = dpsFactory.create("dps_DCorrFrac_"+Fitters[i], 2);
+ }
+
+ Ein_prev = 0;
+ firstEvent = true;
+ first = true;
+ gauss = functionFactory.createFunctionByName("gauss", "G");
+ }
+
+ @Override
+ protected void process(EventHeader event) {
+ E_in = 0.0;
+ E_kin = 0.0;
+ double Chrg = 0.0;
+ double P_in = 0.0;
+ Hep3Vector origin = new BasicHep3Vector(0,0,0);
+ String DirName = null;
+ List<MCParticle> particles = event.get(MCParticle.class, event.MC_PARTICLES);
+ for (MCParticle particle : particles) {
+ if (particle.getParents().isEmpty()) {
+ E_in = E_in + particle.getEnergy();
+ E_kin = E_kin + particle.getEnergy() - particle.getMass();
+ P_in = particle.getMomentum().magnitude();
+ Particlename = particle.getType().getName();
+ origin = particle.getOrigin();
+ Chrg = particle.getCharge();
+ }
+ break;
+ }
+ if (E_kin < 0.999) {
+ Ein = (int) Math.floor(1000.0*E_kin + 0.5d);
+ } else {
+ Ein = (int) Math.floor(E_kin + 0.5d);
+ }
+ String E_str = Ein.toString();
+ DirName = Particlename.concat(E_str);
+ if (!Ein.equals(Ein_prev)) {
+ if (first) {
+ first = false;
+ } else {
+ System.out.println("E_in: " + E_in);
+ System.out.println("Ein_prev: " + Ein_prev);
+ 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_scor = aida.histogramFactory().createCloud1D("Edep_scor", "single 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);
+ }
+
+ double dif_curv = 0.0;
+ if (Chrg != 0.0) {
+ double rad_curv = P_in * 667.128;
+ double len_curv = rad_curv * 2.0 * Math.asin(0.667 / rad_curv);
+ dif_curv = len_curv - 1.334;
+ }
+ 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.contains("Edep")) {
+ for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+ if (cutCon.T_thresh > 0.0 || cutCon.S_thresh > 0.0) {
+ Hep3Vector diffvec = VecOp.sub(calorimeterHit.getPositionVec(), origin);
+ int ncontribs = calorimeterHit.getMCParticleCount();
+ for (int i=0; i<ncontribs; i++) {
+ double tdep = calorimeterHit.getContributedTime(i);
+ double tofl = tdep - (diffvec.magnitude()+dif_curv)/299.79;
+ double edep = calorimeterHit.getContributedEnergy(i);
+ if ( (tofl < cutCon.T_thresh) && (edep > cutCon.S_thresh/1000.0) ) {
+ sumEEdep += edep;
+ }
+ }
+ }
+ else {
+ sumEEdep = sumEEdep + calorimeterHit.getRawEnergy();
+ }
+ }
+ } // end if Edep
+ if (CollectionName.contains("Opti")) {
+ for (SimCalorimeterHit calorimeterHit : simCalorimeterHits) {
+ if (cutCon.T_thresh > 0.0 || cutCon.C_thresh > 0.0) {
+ Hep3Vector diffvec = VecOp.sub(calorimeterHit.getPositionVec(), origin);
+ int ncontribs = calorimeterHit.getMCParticleCount();
+ for (int i=0; i<ncontribs; i++) {
+ double tdep = calorimeterHit.getContributedTime(i);
+ double tofl = tdep - (diffvec.magnitude()+dif_curv)/299.79;
+ double edep = calorimeterHit.getContributedEnergy(i);
+ if ( (tofl < cutCon.T_thresh) && (edep > cutCon.C_thresh) ) {
+ sumECeren += edep;
+ }
+ }
+ }
+ else {
+ sumECeren = sumECeren + calorimeterHit.getRawEnergy();
+ }
+ }
+ }
+ }
+ Edep.fill(sumEEdep);
+ ECeren.fill(sumECeren);
+
+ xval[0] = sumEEdep;
+ double sumEEdep_cor = corfus.getFunction("SC").value(xval);
+ xval[0] = sumECeren;
+ double sumECeren_cor = corfus.getFunction("CC").value(xval);
+ Edep_cor.fill(sumEEdep_cor);
+ Eceren_cor.fill(sumECeren_cor);
+
+ double ratio = sumECeren_cor / sumEEdep_cor;
+ double fraction = sumEEdep_cor / E_in;
+ xval[0] = ratio;
+ String cfname = "CS" + Ein.toString() + "gev";
+ double cfac = corfus.getFunction(cfname).value(xval);
+ Edep_dcor.fill(sumEEdep_cor / cfac);
+ String sfname = "S" + Ein.toString() + "gev";
+ double sfac = corfus.getFunction(sfname).parameter("mean");
+ Edep_scor.fill(sumEEdep_cor / sfac);
+ }
+
+ @Override
+ protected void endOfData() {
+ System.out.println("End of Data:");
+ convertandfit();
+ try {
+ aida.saveAs(AIDAFile);
+ } catch (IOException ex) {
+ Logger.getLogger(Resolution.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_conv;
+ if (Edep_dcor.isConverted()) {
+ Edep_dcor_conv = Edep_dcor.histogram();
+ } else {
+ double meanc = Edep_dcor.mean();
+ double rmsc = Edep_dcor.rms();
+ double minx = meanc - nsigmas * rmsc;
+ double maxx = meanc + nsigmas * rmsc;
+ Edep_dcor.setConversionParameters(nbins, minx, maxx);
+ Edep_dcor.convertToHistogram();
+ Edep_dcor_conv = Edep_dcor.histogram();
+ }
+ gauss.setParameter("amplitude", Edep_dcor_conv.maxBinHeight());
+ gauss.setParameter("mean", Edep_dcor_conv.mean());
+ gauss.setParameter("sigma", Edep_dcor_conv.rms());
+ System.out.print(Particlename + " " + Ein_prev + " GeV\n");
+ System.out.print(Edep_dcor_conv.maxBinHeight() + "," + Edep_dcor_conv.mean() + "," + Edep_dcor_conv.rms() + "\n");
+
+ for (int i = 0; i < Fitters.length; i++) {
+ System.out.println("Fitter: " + Fitters[i]);
+ gauss.setParameter("amplitude", Edep_dcor_conv.maxBinHeight());
+ gauss.setParameter("mean", Edep_dcor_conv.mean());
+ gauss.setParameter("sigma", Edep_dcor_conv.rms());
+ jminuit = fitFactory.createFitter(Fitters[i], "jminuit");
+ jminuitResult = jminuit.fit(Edep_dcor_conv, gauss);
+ System.out.println("jminuit " + Fitters[i] + ": " + jminuitResult.quality());
+ result = jminuitResult.fittedParameters();
+ errors = jminuitResult.errors();
+ double[] errsplus = jminuitResult.errorsPlus();
+ double[] errsminu = jminuitResult.errorsMinus();
+ String functionname = "Corrected fitted gauss " + Fitters[i];
+ System.out.println(functionname);
+ functionFactory.cloneFunction(functionname, jminuitResult.fittedFunction());
+ 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]++;
+ IDataPoint dp_Resolution = dps_DRResolution[i].addPoint();
+ double sigma = Math.abs(result[2]);
+ dp_Resolution.coordinate(0).setValue(1.0/Math.sqrt(Ein_prev));
+ dp_Resolution.coordinate(1).setValue(Math.abs(sigma)/Ein_prev);
+ dp_Resolution.coordinate(1).setErrorMinus(errsplus[2]/Ein_prev);
+ dp_Resolution.coordinate(1).setErrorPlus(errsminu[2]/Ein_prev);
+ IDataPoint dp_DCorrFrac = dps_DCorrFrac[i].addPoint();
+ double mean = result[1];
+ dp_DCorrFrac.coordinate(0).setValue(Ein_prev);
+ dp_DCorrFrac.coordinate(1).setValue(mean/Ein_prev);
+ dp_DCorrFrac.coordinate(1).setErrorMinus(sigma);
+ dp_DCorrFrac.coordinate(1).setErrorPlus(sigma);
+
+ IHistogram1D Edep_scor_conv;
+ if (Edep_scor.isConverted()) {
+ Edep_scor_conv = Edep_scor.histogram();
+ } else {
+ double meanc = Edep_scor.mean();
+ double rmsc = Edep_scor.rms();
+ double minx = meanc - nsigmas * rmsc;
+ double maxx = meanc + nsigmas * rmsc;
+ Edep_scor.setConversionParameters(nbins, minx, maxx);
+ Edep_scor.convertToHistogram();
+ Edep_scor_conv = Edep_scor.histogram();
+ }
+ gauss.setParameter("amplitude", Edep_scor_conv.maxBinHeight());
+ gauss.setParameter("mean", Edep_scor_conv.mean());
+ gauss.setParameter("sigma", Edep_scor_conv.rms());
+ System.out.print(Particlename + " " + Ein_prev + " GeV\n");
+ System.out.print(Edep_scor_conv.maxBinHeight() + "," + Edep_scor_conv.mean() + "," + Edep_scor_conv.rms() + "\n");
+ jminuit = fitFactory.createFitter("Chi2", "jminuit");
+ jminuitResult = jminuit.fit(Edep_dcor_conv, gauss);
+ System.out.println("jminuit " + Fitters[i] + ": " + jminuitResult.quality());
+ result = jminuitResult.fittedParameters();
+ errors = jminuitResult.errors();
+ IDataPoint dp_SRResolution = dps_SRResolution[i].addPoint();
+ sigma = Math.abs(result[2]);
+ dp_Resolution.coordinate(0).setValue(1.0/Math.sqrt(Ein_prev));
+ dp_Resolution.coordinate(1).setValue(Math.abs(sigma)/Ein_prev);
+ dp_Resolution.coordinate(1).setErrorMinus(errsplus[2]/Ein_prev);
+ dp_Resolution.coordinate(1).setErrorPlus(errsminu[2]/Ein_prev);
+ }
+
+
+ try {
+ aida.saveAs(AIDAFile);
+ } catch (IOException ex) {
+ Logger.getLogger(Resolution.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+
+ public void setCorrFu(CorrectionFunctionSet corrs) {
+ System.out.println("DualCorrection:setCorrFu");
+ this.corfus = corrs;
+ }
+
+ public void setCutConfiguration(cutConfiguration cutco) {
+ System.out.println("DualCorrection:setCutConfiguration");
+ this.cutCon = cutco;
+ }
+
+ public void setMyAIDAFilename(String AIDAFile) {
+ System.out.println("DualCorrection:setMyAIDAFilename");
+ this.AIDAFile = AIDAFile;
+ System.out.println(AIDAFile);
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
diff -N proton_calibration_steering.xml
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ proton_calibration_steering.xml 18 Apr 2014 20:52:20 -0000 1.1
@@ -0,0 +1,33 @@
+<lcsim xmlns:xs="http://www.w3.org/2001/XMLSchema-instance"
+ xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/lcsim/1.0/lcsim.xsd">
+
+ <classpath>
+ <jar>/home/aconway/NetBeansProjects/mcd-analysis/target/mcd-analysis-1.0-SNAPSHOT.jar</jar>
+ </classpath>
+
+ <inputFiles>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-1gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-2gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-5gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-10gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-25gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/proton/proton-90deg-50gev.slcio</file>
+ </inputFiles>
+ <control>
+ <numberOfEvents>-1</numberOfEvents>
+ </control>
+
+ <execute>
+ <driver name="DRCalibrationDriver"/>
+ </execute>
+ <drivers>
+ <driver name="DRCalibrationDriver"
+ type="org.lcsim.mcd.analysis.DRCorrection.alexDRCorr.DRCalibrationDriver">
+ <myP_name>proton</myP_name>
+ <myS_thresh>0.0</myS_thresh>
+ <myC_thresh>0.0</myC_thresh>
+ <myT_thresh>3.0</myT_thresh>
+ <myTag>default</myTag>
+ </driver>
+ </drivers>
+</lcsim>
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
diff -N cutConfiguration.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ cutConfiguration.java 18 Apr 2014 20:52:20 -0000 1.1
@@ -0,0 +1,83 @@
+/*
+ * To change this license header, choose License Headers in Project Properties.
+ * To change this template file, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.mcd.analysis.DRCorrection.alexDRCorr;
+
+/**
+ *
+ * @author aconway
+ */
+public class cutConfiguration {
+ public Double C_thresh;
+ public Double S_thresh;
+ public Double T_thresh;
+ public String P_name;
+ public String tag;
+
+ public cutConfiguration(Double Ct, Double St, Double Tt, String name) {
+ C_thresh = Ct;
+ S_thresh = St;
+ T_thresh = Tt;
+ P_name = name;
+ tag = "default";
+ }
+
+ public cutConfiguration(Double Ct, Double St, Double Tt, String name, String t) {
+ C_thresh = Ct;
+ S_thresh = St;
+ T_thresh = Tt;
+ P_name = name;
+ tag = t;
+ }
+
+ @Override
+ public boolean equals(Object o) {
+ if ((o instanceof cutConfiguration) &&
+ (((cutConfiguration) o).P_name.equals(this.P_name)) &&
+ (((cutConfiguration) o).tag.equals(this.tag)) &&
+ (((cutConfiguration) o).C_thresh.compareTo(this.C_thresh) == 0) &&
+ (((cutConfiguration) o).S_thresh.compareTo(this.S_thresh) == 0) &&
+ (((cutConfiguration) o).T_thresh.compareTo(this.T_thresh) == 0))
+ {
+ return true;
+ } else {
+ return false;
+ }
+ }
+
+ @Override
+ public int hashCode() {
+ int hash = 7;
+ hash = 53 * hash + (this.P_name != null ? this.P_name.hashCode() : 0);
+ hash = 53 * hash + (this.C_thresh != null ? this.C_thresh.hashCode() : 0);
+ hash = 53 * hash + (this.S_thresh != null ? this.S_thresh.hashCode() : 0);
+ hash = 53 * hash + (this.T_thresh != null ? this.T_thresh.hashCode() : 0);
+ return hash;
+ }
+
+ @Override
+ public String toString() {
+ String rv = P_name + "-" +
+ String.valueOf(this.T_thresh) + "ns-" +
+ String.valueOf(this.S_thresh) + "mev-" +
+ String.valueOf(this.C_thresh);
+ if (!this.tag.equals("default")) {
+ rv = rv + "-" + this.tag;
+ }
+ return rv;
+ }
+
+ public String toStringAnyP() {
+ String rv = "any" + "-" +
+ String.valueOf(this.T_thresh) + "ns-" +
+ String.valueOf(this.S_thresh) + "mev-" +
+ String.valueOf(this.C_thresh);
+ if (!this.tag.equals("default")) {
+ rv = rv + "-" + this.tag;
+ }
+ return rv;
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
diff -N electron_calibration_steering.xml
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ electron_calibration_steering.xml 18 Apr 2014 20:52:20 -0000 1.1
@@ -0,0 +1,32 @@
+<lcsim xmlns:xs="http://www.w3.org/2001/XMLSchema-instance"
+ xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/lcsim/1.0/lcsim.xsd">
+
+ <classpath>
+ <jar>/home/aconway/NetBeansProjects/mcd-analysis/target/mcd-analysis-1.0-SNAPSHOT.jar</jar>
+ </classpath>
+ <inputFiles>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/electron/electron-90deg-1gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/electron/electron-90deg-2gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/electron/electron-90deg-5gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/electron/electron-90deg-10gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/electron/electron-90deg-25gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/electron/electron-90deg-50gev.slcio</file>
+ <file>/home/aconway/fermi/eventdata/slcio/singleparticle/electron/electron-90deg-100gev.slcio</file>
+ </inputFiles>
+ <control>
+ <numberOfEvents>-1</numberOfEvents>
+ </control>
+ <execute>
+ <driver name="ElectronCalibrationDriver"/>
+ </execute>
+ <drivers>
+ <driver name="ElectronCalibrationDriver"
+ type="org.lcsim.mcd.analysis.DRCorrection.alexDRCorr.ElectronCalibrationDriver">
+ <myP_name>any</myP_name>
+ <myS_thresh>0.0</myS_thresh>
+ <myC_thresh>0.0</myC_thresh>
+ <myT_thresh>0.0</myT_thresh>
+ <myTag>default</myTag>
+ </driver>
+ </drivers>
+</lcsim>
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection/alexDRCorr
diff -N DRCalibrationDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ DRCalibrationDriver.java 18 Apr 2014 20:52:20 -0000 1.1
@@ -0,0 +1,65 @@
+package org.lcsim.mcd.analysis.DRCorrection.alexDRCorr;
+
+import org.lcsim.util.Driver;
+
+/**
+ *
+ * @author aconway
+ */
+public class DRCalibrationDriver extends Driver {
+
+ private String P_name = "any";
+ private String tag = "default";
+ private Double C_thresh = 0.0;
+ private Double S_thresh = 0.0;
+ private Double T_thresh = 0.0;
+ private String AIDAFile = null;
+
+ DualCorrection dual;
+
+ public DRCalibrationDriver() {
+ dual = new DualCorrection();
+ add(dual);
+ }
+
+ @Override
+ public void startOfData() {
+ System.out.println("DRCalibrationDriver:startOfData");
+ cutConfiguration cutco;
+ if (tag.equals("default")) {
+ cutco = new cutConfiguration(C_thresh, S_thresh, T_thresh, P_name);
+ } else {
+ cutco = new cutConfiguration(C_thresh, S_thresh, T_thresh, P_name, tag);
+ }
+ AIDAFile = "DR-" + cutco.toString() + ".aida";
+ dual.setMyAIDAFilename(AIDAFile);
+ dual.setCutConfiguration(cutco);
+ dual.setCorrFu(new CorrectionFunctionSet(cutco, true));
+ dual.startOfData();
+ }
+
+ public void setMyP_name(String pname) {
+ System.out.println("ElectronCalibrationDriver:setMyP_name = " + pname);
+ this.P_name = pname;
+ }
+
+ public void setMyTag(String tag) {
+ System.out.println("ElectronCalibrationDriver:setMyTag = " + tag);
+ this.tag = tag;
+ }
+
+ public void setMyC_thresh(double cthresh) {
+ System.out.println("ElectronCalibrationDriver:setMyC_thresh = " + cthresh);
+ this.C_thresh = cthresh;
+ }
+
+ public void setMyS_thresh(double sthresh) {
+ System.out.println("ElectronCalibrationDriver:setMyS_thresh = " + sthresh);
+ this.S_thresh = sthresh;
+ }
+
+ public void setMyT_thresh(double tthresh) {
+ System.out.println("ElectronCalibrationDriver:setMyT_thresh = " + tthresh);
+ this.T_thresh = tthresh;
+ }
+}
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection
diff -u -r1.2 -r1.3
--- DualCorrection.java 7 Apr 2014 15:13:48 -0000 1.2
+++ DualCorrection.java 18 Apr 2014 20:52:20 -0000 1.3
@@ -434,15 +434,24 @@
if (entries > 300) {
for (int ii = 0; ii < Fitters.length; ii++) {
System.out.println("Fitter: " + Fitters[ii]);
+ double meanguess = 0;
+ for (int j=0; j<conv_slices[i].axis().bins(); j++) {
+ if (conv_slices[i].binHeight(j) == conv_slices[i].maxBinHeight()) {
+ meanguess = conv_slices[i].axis().binCenter(j);
+ }
+ }
gauss.setParameter("amplitude", conv_slices[i].maxBinHeight());
- gauss.setParameter("mean", conv_slices[i].mean());
- gauss.setParameter("sigma", conv_slices[i].rms());
+ gauss.setParameter("mean", meanguess);
+ gauss.setParameter("sigma", conv_slices[i].rms()/2.0);
jminuit = fitFactory.createFitter(Fitters[ii], "jminuit");
try {
jminuitResult = jminuit.fit(conv_slices[i], gauss);
} catch (RuntimeException ex) {
System.out.println("Runtime Exception caught in fit!");
System.out.println(ex.toString());
+ } catch (Exception ex) {
+ System.out.println("Exception caught in f!");
+ System.out.println(ex.toString());
}
System.out.println("jminuit " + Fitters[ii] + ": " + jminuitResult.quality());
result = jminuitResult.fittedParameters();
mcd-analysis/src/main/java/org/lcsim/mcd/analysis/DRCorrection
diff -u -r1.2 -r1.3
--- ElectronCorrection.java 7 Apr 2014 15:13:48 -0000 1.2
+++ ElectronCorrection.java 18 Apr 2014 20:52:20 -0000 1.3
@@ -391,4 +391,4 @@
this.file_name = file_name;
System.out.println(file_name);
}
-}
+}
\ No newline at end of file
mcd-analysis/src/main/java/org/lcsim/mcd/analysis
diff -u -r1.1 -r1.2
--- ILCstdhepdriver.java 7 Apr 2014 15:13:49 -0000 1.1
+++ ILCstdhepdriver.java 18 Apr 2014 20:52:20 -0000 1.2
@@ -1,17 +1,15 @@
-/*
- * To change this license header, choose License Headers in Project Properties.
- * To change this template file, choose Tools | Templates
- * and open the template in the editor.
- */
-
//package org.lcsim.mcd.analysis;
+import hep.aida.IAxis;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
import java.util.List;
import java.util.ArrayList;
import org.lcsim.event.EventHeader;
import org.lcsim.event.MCParticle;
-import org.lcsim.event.base.BaseMCParticle;
import org.lcsim.util.Driver;
import org.lcsim.util.aida.AIDA;
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogram2D;
/**
*
@@ -21,6 +19,26 @@
public class ILCstdhepdriver extends Driver {
private AIDA aida = AIDA.defaultInstance();
+ private List<String> hadrons = new ArrayList<String>();
+
+ IHistogram1D allhadsen = aida.histogram1D("sum energy wtd by frac energy",
+ 25,0,50);
+ IHistogram2D allhadsfrac = aida.histogram2D("energy vs frac energy",
+ 20,0,50,
+ 25,0.0,0.5);
+
+ public ILCstdhepdriver() {
+ hadrons.add("p");
+ hadrons.add("p_bar");
+ hadrons.add("pi+");
+ hadrons.add("pi-");
+ hadrons.add("n");
+ hadrons.add("n_bar");
+ hadrons.add("K+");
+ hadrons.add("K-");
+ hadrons.add("K0_L");
+ }
+
protected void process(EventHeader event) {
List<MCParticle> particles = event.get(MCParticle.class,"MCParticle");
List<MCParticle> children = new ArrayList<MCParticle>();
@@ -40,14 +58,105 @@
List<MCParticle> fsChildren = getAllFSChildren(higgs);
String pfix = "/"+decayname+"/";
+ double sumEShower = 0.0;
+ double sumEBarrel = 0.0;
+ double sumEHadron = 0.0;
+ for (MCParticle p : fsChildren) {
+ sumEShower += p.getEnergy();
+ double logen = Math.log10(p.getEnergy());
+ double logp = Math.log10(p.getMomentum().magnitude());
+ if (hitsBarrel(p)) {
+ sumEBarrel += p.getEnergy();
+ String name = p.getType().getName();
+ if (hadrons.contains(name)) {
+ sumEHadron += p.getEnergy();
+ }
+ aida.cloud1D(pfix+"hits/"+p.getType().getName()+"/energy").fill(p.getEnergy());
+ aida.cloud1D(pfix+"hits/"+p.getType().getName()+"/momentum").fill(p.getMomentum().magnitude());
+ aida.cloud1D(pfix+"hits/"+p.getType().getName()+"/log energy").fill(logen);
+ aida.cloud1D(pfix+"hits/"+p.getType().getName()+"/log momentum").fill(logp);
+ } else {
+ aida.cloud1D(pfix+"misses/"+p.getType().getName()+"/energy").fill(p.getEnergy());
+ aida.cloud1D(pfix+"misses/"+p.getType().getName()+"/momentum").fill(p.getMomentum().magnitude());
+ aida.cloud1D(pfix+"misses/"+p.getType().getName()+"/log energy").fill(logen);
+ aida.cloud1D(pfix+"misses/"+p.getType().getName()+"/log momentum").fill(logp);
+ }
+ }
+ aida.cloud1D(pfix+"tot shower energy").fill(sumEShower);
+ aida.cloud1D(pfix+"frac barrel energy").fill(sumEBarrel/sumEShower);
+ List<String> pnames = new ArrayList<String>();
for (MCParticle p : fsChildren) {
String particlename = p.getType().getName();
aida.cloud1D(pfix+particlename+"/momentum").fill(p.getMomentum().magnitude());
aida.cloud1D(pfix+particlename+"/energy").fill(p.getEnergy());
aida.cloud1D(pfix+particlename+"/log momentum").fill(Math.log10(p.getMomentum().magnitude()));
aida.cloud1D(pfix+particlename+"/log energy").fill(Math.log10(p.getEnergy()));
- // TODO: look at relative composition of shower particles and energies. Is most of the energy in an event in low energy hadrons?
-
+ if (hitsBarrel(p)) {
+ aida.cloud1D(pfix+particlename+"/log energy wtd by frac energy").fill(Math.log10(p.getEnergy()), p.getEnergy()/sumEBarrel);
+ aida.cloud1D(pfix+particlename+"/sum log energy wtd by frac energy").fill(Math.log10(p.getEnergy()), p.getEnergy()/sumEBarrel);
+ aida.cloud1D(pfix+"all/sum log energy wtd by frac energy").fill(Math.log10(p.getEnergy()), p.getEnergy()/sumEBarrel);
+ aida.cloud1D(pfix+particlename+"/energy wtd by frac energy").fill(p.getEnergy(), p.getEnergy()/sumEBarrel);
+ aida.cloud1D(pfix+particlename+"/sum energy wtd by frac energy").fill(p.getEnergy(), p.getEnergy()/sumEBarrel);
+ aida.cloud1D(pfix+"all/sum energy wtd by frac energy").fill(p.getEnergy(), p.getEnergy()/sumEBarrel);
+ if (!pnames.contains(particlename)) {
+ pnames.add(particlename);
+ }
+ if (hadrons.contains(particlename)) {
+ aida.cloud1D(pfix+"allhads/sum energy wtd by frac energy").fill(p.getEnergy(), p.getEnergy()/sumEHadron);
+ allhadsen.fill(p.getEnergy(),p.getEnergy()/sumEHadron);
+ }
+ }
+ }
+ for (String pname : pnames) {
+ String opname = pfix+pname+"/energy wtd by frac energy";
+ if (!aida.cloud1D(opname).isConverted()) {
+ aida.cloud1D(opname).convertToHistogram();
+ }
+ IHistogram1D ophist = aida.cloud1D(opname).histogram();
+ IAxis opax = ophist.axis();
+ for (int i=0; i<opax.bins(); i++) {
+ double frac = ophist.binHeight(i);
+ double en = ophist.binMean(i);
+ if (frac > 0.005) {
+ aida.cloud2D(pfix+pname+"/log energy vs frac energy").fill(en,frac);
+ }
+ }
+ aida.cloud1D(opname).reset();
+ }
+ String opname = pfix+"allhads/sum energy wtd by frac energy";
+ if (!aida.cloud1D(opname).isConverted()) {
+ aida.cloud1D(opname).convertToHistogram();
+ }
+ IHistogram1D ophist = aida.cloud1D(opname).histogram();
+ IAxis opax = ophist.axis();
+ for (int i=0; i<opax.bins(); i++) {
+ double frac = ophist.binHeight(i);
+ double en = ophist.axis().binCenter(i);
+ if (frac > 0.00) {
+ aida.cloud2D(pfix+"allhads/log energy vs frac energy").fill(en,frac);
+ }
+ }
+ aida.cloud1D(opname).reset();
+ for (int i=0; i<allhadsen.axis().bins(); i++) {
+ double frac = allhadsen.binHeight(i);
+ double en = allhadsen.axis().binCenter(i);
+ if (frac > 0.00005) {
+ aida.cloud2D("cloud energy vs frac energy").fill(en, frac);
+ allhadsfrac.fill(en,frac);
+ }
+ }
+ allhadsen.reset();
+ }
+
+ public boolean hitsBarrel (MCParticle mcp) {
+ if (mcp.getCharge() == 0.0) {
+ return true;
+ } else {
+ Hep3Vector origin = new BasicHep3Vector(mcp.getOriginX(),mcp.getOriginY(),0.0);
+ Hep3Vector pt = new BasicHep3Vector(mcp.getMomentum().x(),mcp.getMomentum().y(),0.0);
+ // diameter of curvature in 5T field in mm
+ double diam = pt.magnitude()*1334.256;
+ return diam > 1350.0 - origin.magnitude();
}
}
CVSspam 0.2.12