Commit in mcd-analysis/src/main/java/org/lcsim/mcd/analysis on MAIN | |||
DRCorrection/alexDRCorr/proton_resolution_steering.xml | +33 | added 1.1 | |
/CorrectionFunctionSet.java | +171 | added 1.1 | |
/correctionFunctions.java | +65 | added 1.1 | |
/electronCorrection.java | +326 | added 1.1 | |
/ElectronCalibrationDriver.java | +64 | added 1.1 | |
/DRResolutionDriver.java | +65 | added 1.1 | |
/DualCorrection.java | +353 | added 1.1 | |
/Resolution.java | +359 | added 1.1 | |
/proton_calibration_steering.xml | +33 | added 1.1 | |
/cutConfiguration.java | +83 | added 1.1 | |
/electron_calibration_steering.xml | +32 | added 1.1 | |
/DRCalibrationDriver.java | +65 | added 1.1 | |
DRCorrection/DualCorrection.java | +11 | -2 | 1.2 -> 1.3 |
/ElectronCorrection.java | +1 | -1 | 1.2 -> 1.3 |
ILCstdhepdriver.java | +118 | -9 | 1.1 -> 1.2 |
+1779 | -12 |
Added new package, alexDRCorr, to implement certain changes to the analysis. C/S correction functions are now polynomials and stored as files and the 'DetectorConfiguration' is replaced by a simpler, more appropriate interface. I also use an S-only correction to compare against C/S.
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>
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); + } + } +}
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); + } +}
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); + } +}
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; + } +}
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; + } +}
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); + } +}
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); + } +}
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>
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; + } +}
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>
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; + } +}
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();
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
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();
} }
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1