Author: [log in to unmask]
Date: Thu Sep 1 04:27:30 2016
New Revision: 4478
Log: (empty)
Added:
java/trunk/users/src/main/java/org/hps/users/spaul/AidaHistogramUtil.java
java/trunk/users/src/main/java/org/hps/users/spaul/BestTrack.java
java/trunk/users/src/main/java/org/hps/users/spaul/PeakPositionVsTime.java
java/trunk/users/src/main/java/org/hps/users/spaul/PrintAidaTree.java
java/trunk/users/src/main/java/org/hps/users/spaul/PrintStuff.java
java/trunk/users/src/main/java/org/hps/users/spaul/RunningPedestalPrintDriver.java
java/trunk/users/src/main/java/org/hps/users/spaul/ScalerAttenuationCalculator.java
java/trunk/users/src/main/java/org/hps/users/spaul/ScalerZeropointCalculator.java
java/trunk/users/src/main/java/org/hps/users/spaul/TestAIDA2ROOT.java
java/trunk/users/src/main/java/org/hps/users/spaul/bh/
java/trunk/users/src/main/java/org/hps/users/spaul/bh/BumpHunt.java
java/trunk/users/src/main/java/org/hps/users/spaul/bh/Fitters.java
java/trunk/users/src/main/java/org/hps/users/spaul/bh/GenerateFakeStuff.java
java/trunk/users/src/main/java/org/hps/users/spaul/bh/GetEpsilon2.java
java/trunk/users/src/main/java/org/hps/users/spaul/bh/GetPullDistribution.java
java/trunk/users/src/main/java/org/hps/users/spaul/bh/GetYieldPerSignal.java
java/trunk/users/src/main/java/org/hps/users/spaul/bh/GoodnessOfFitAccumulator.java
java/trunk/users/src/main/java/org/hps/users/spaul/bh/Histogram.java
java/trunk/users/src/main/java/org/hps/users/spaul/bh/IFitResult.java
java/trunk/users/src/main/java/org/hps/users/spaul/bh/LogLikelihood.java
java/trunk/users/src/main/java/org/hps/users/spaul/bh/Polynomial.java
java/trunk/users/src/main/java/org/hps/users/spaul/bh/PolynomialPlusGaussian.java
java/trunk/users/src/main/java/org/hps/users/spaul/bh/WindowResult.java
java/trunk/users/src/main/java/org/hps/users/spaul/bh/WindowSet.java
java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/
java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/Cuts.java
java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/GenerateFakeResolutions.java
java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/InvariantMassCalculator.java
java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/MakeMassHistogram.java
java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/PlotMassResolutionsMC.java
java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/RecoilCalculator.java
java/trunk/users/src/main/java/org/hps/users/spaul/dimuon/
java/trunk/users/src/main/java/org/hps/users/spaul/dimuon/DimuonFilterDriver.java
java/trunk/users/src/main/java/org/hps/users/spaul/pulser/
java/trunk/users/src/main/java/org/hps/users/spaul/pulser/UncorrClusterDriver.java
Added: java/trunk/users/src/main/java/org/hps/users/spaul/AidaHistogramUtil.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/AidaHistogramUtil.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/AidaHistogramUtil.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,31 @@
+package org.hps.users.spaul;
+
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogramFactory;
+import hep.aida.ref.histogram.Histogram1D;
+
+public class AidaHistogramUtil {
+ public static IHistogram1D shift(IHistogram1D h, final int nbins){
+ return new Histogram1D(){
+ public double binHeight(int i){
+ return super.binHeight(i-nbins);
+ }
+ public double binError(int i){
+ return super.binError(i-nbins);
+ }
+ };
+ }
+ public static IHistogram1D shift(IHistogramFactory hf, IHistogram1D h, final int nbins){
+ IHistogram1D copy = hf.createCopy("copy", h);
+ copy.reset();
+ if(nbins > 0)
+ for(int i = 0; i< copy.axis().bins()-nbins; i++){
+ copy.fill(h.binMean(i), h.binHeight(i+nbins));
+ }
+ else
+ for(int i = -nbins; i< copy.axis().bins(); i++){
+ copy.fill(h.binMean(i), h.binHeight(i+nbins));
+ }
+ return copy;
+ }
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/BestTrack.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/BestTrack.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/BestTrack.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,25 @@
+package org.hps.users.spaul;
+
+import java.util.List;
+
+import org.lcsim.event.Track;
+import org.lcsim.event.TrackerHit;
+
+public class BestTrack {
+ static boolean shareHits(Track t1, Track t2){
+ for(TrackerHit hit : t1.getTrackerHits()){
+ for(TrackerHit hit2 : t2.getTrackerHits()){
+ if(equals(hit, hit2))
+ return true;
+ }
+ }
+ return false;
+ }
+
+ static double threshold = 3;
+ static boolean equals(TrackerHit hit1, TrackerHit hit2){
+ return Math.hypot(hit1.getPosition()[0]-hit2.getPosition()[0],
+ hit1.getPosition()[1]-hit2.getPosition()[1])< threshold;
+ }
+
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/PeakPositionVsTime.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/PeakPositionVsTime.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/PeakPositionVsTime.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,103 @@
+package org.hps.users.spaul;
+
+import java.util.List;
+
+import hep.aida.IHistogram1D;
+
+import org.hps.conditions.beam.BeamEnergy.BeamEnergyCollection;
+import org.hps.recon.ecal.EcalRawConverter;
+import org.hps.record.triggerbank.AbstractIntData;
+import org.hps.record.triggerbank.TIData;
+import org.lcsim.conditions.ConditionsManager;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.GenericObject;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+public class PeakPositionVsTime extends Driver{
+
+
+ AIDA aida = AIDA.defaultInstance();
+
+ private int skipEvents = 10000;
+ private int samplesPerPrint = 1000;
+ IHistogram1D radPeak, feePeak;
+ public void setSamplesPerPrint(int samplesPerPrint){
+ this.samplesPerPrint = samplesPerPrint;
+ }
+ public void setSkipEvents(int skipEvents){
+ this.skipEvents = skipEvents;
+ }
+ double beamEnergy = 0;
+ public void startOfData(){
+
+ BeamEnergyCollection beamEnergyCollection =
+ this.getConditionsManager().getCachedConditions(BeamEnergyCollection.class, "beam_energies").getCachedData();
+ if(beamEnergyCollection != null && beamEnergyCollection.size() != 0)
+ beamEnergy = beamEnergyCollection.get(0).getBeamEnergy();
+ radPeak = aida.histogram1D("radpeak", 200, 0, 1.5*beamEnergy);
+ feePeak = aida.histogram1D("feepeak", 200, 0, 1.5*beamEnergy);
+
+ }
+
+
+ @Override
+ public void process(EventHeader event){
+ int eventNum = event.getEventNumber();
+ if(eventNum < skipEvents)
+ return;
+ if(eventNum % samplesPerPrint != 0){
+ double feePeakValue = getMax(feePeak);
+ double radPeakValue = getMax(radPeak);
+ long time = event.getTimeStamp();
+ System.out.println(time + "\t" + eventNum + "\t" + feePeakValue + "\t" + radPeakValue);
+ radPeak.reset();
+ feePeak.reset();
+ }
+
+
+ boolean singles1 = false, pairs1 = false;
+ for (GenericObject gob : event.get(GenericObject.class,"TriggerBank"))
+ {
+ if (!(AbstractIntData.getTag(gob) == TIData.BANK_TAG)) continue;
+ TIData tid = new TIData(gob);
+ if (tid.isSingle1Trigger()) singles1 = true;
+ if (tid.isPair1Trigger()) pairs1 = true;
+ }
+
+
+ if(! singles1 && ! pairs1)
+ return;
+ List<Cluster> clusters = event.get(Cluster.class, "EcalClustersCorr");
+ if(singles1){
+ for(Cluster cluster : clusters){
+ feePeak.fill(cluster.getEnergy());
+ }
+ }
+ if(pairs1){
+ for(Cluster c1 : clusters){
+ for(Cluster c2 : clusters){
+ double e1 = c1.getEnergy();
+ double e2 = c2.getEnergy();
+ if(e1 < .8*beamEnergy && e2 < .8*beamEnergy //avoid FEEs
+ && c1.getPosition()[1]*c2.getPosition()[1]<0 //both on opposite halves of detector
+ ){
+ radPeak.fill(e1+e2);
+ }
+ }
+ }
+ }
+
+ }
+ double getMax(IHistogram1D hist){
+ int peakBin = 0;
+ for(int i = 0; i< hist.axis().bins(); i++){
+ if(hist.binHeight(i) > hist.binHeight(peakBin))
+ peakBin = i;
+ }
+
+
+ return (hist.axis().binLowerEdge(peakBin)+hist.axis().binUpperEdge(peakBin))/2.;
+ }
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/PrintAidaTree.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/PrintAidaTree.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/PrintAidaTree.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,23 @@
+package org.hps.users.spaul;
+
+import java.io.IOException;
+import java.util.Arrays;
+
+import hep.aida.IAnalysisFactory;
+import hep.aida.IHistogramFactory;
+import hep.aida.ITree;
+
+public class PrintAidaTree {
+ public static void main(String arg[]) throws IllegalArgumentException, IOException{
+ String input = arg[0];
+ IAnalysisFactory af = IAnalysisFactory.create();
+ System.out.println("opening file..");
+ ITree tree = af.createTreeFactory().create(input, "xml", true);
+ System.out.println("opened file");
+ String[] names = tree.listObjectNames(".", true);
+ String[] types = tree.listObjectTypes(".", true);
+ for(int i = 0; i<names.length; i++){
+ System.out.println(names[i] + "\t" + types[i]);
+ }
+ }
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/PrintStuff.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/PrintStuff.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/PrintStuff.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,45 @@
+package org.hps.users.spaul;
+
+import java.util.List;
+import java.util.Map;
+
+import org.hps.datacat.DatacatUtilities;
+import org.hps.datacat.EvioDatasetIndex;
+
+import org.srs.datacat.model.DatasetModel;
+import org.srs.datacat.model.DatasetResultSetModel;
+import org.srs.datacat.model.dataset.DatasetWithViewModel;
+
+public class PrintStuff {
+ public static void main(String arg[]){
+ /*for(int run = 7000; run <= 8100; run++){
+ EvioDatasetIndex edi;
+ try{
+ edi = new EvioDatasetIndex(new DatacatUtilities(), run);
+ }catch(Exception e){
+ continue;
+ }
+ for(int fileNumber = 0; fileNumber < edi.getDatasets().getCount(); fileNumber++){
+ DatasetModel dataset = edi.findByFileNumber(fileNumber);
+ DatasetWithViewModel datasetView = (DatasetWithViewModel) dataset;
+ if(datasetView == null)
+ continue;
+ Map<String, Object> metadata = datasetView.getMetadataMap();
+
+ long firstTimestamp = (Long) metadata.get("FIRST_HEAD_TIMESTAMP");
+ long lastTimestamp = (Long) metadata.get("LAST_HEAD_TIMESTAMP");
+
+ long tiTimeOffset = Long.parseLong((String) metadata.get("TI_TIME_MIN_OFFSET"));
+ //long fileNumber = (Long) metadata.get("FILE");
+ long firstPhysicsEvent = (Long) metadata.get("FIRST_PHYSICS_EVENT");
+ long lastPhysicsEvent = (Long) metadata.get("LAST_PHYSICS_EVENT");
+ long nEvents = lastPhysicsEvent- firstPhysicsEvent+1;
+
+ int firstTI = 0;
+ int lastTI = 0;
+ System.out.printf("%d,%d,%d,0,%d,%d,%d,%d,%d\n", run, fileNumber, nEvents, firstTimestamp, lastTimestamp, firstTI, lastTI, tiTimeOffset);
+
+ }
+ }*/
+ }
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/RunningPedestalPrintDriver.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/RunningPedestalPrintDriver.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/RunningPedestalPrintDriver.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,38 @@
+package org.hps.users.spaul;
+
+import org.hps.recon.ecal.EcalRawConverter;
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+public class RunningPedestalPrintDriver extends Driver{
+ private static int NCHANNELS = 422;
+ AIDA aida = AIDA.defaultInstance();
+ private int skipEvents = 10000;
+ private int prescale = 1000;
+ public void setPrescale(int prescale){
+ this.prescale = prescale;
+ }
+ public void setSkipEvents(int skipEvents){
+ this.skipEvents = skipEvents;
+ }
+
+ @Override
+ public void process(EventHeader event){
+ int eventNum = event.getEventNumber();
+ if(eventNum < skipEvents)
+ return;
+ if(eventNum % prescale != 0)
+ return;
+ long time = event.getTimeStamp();
+ EcalRawConverter erc = new EcalRawConverter();
+ double[] pedestals = new double[NCHANNELS];
+ for(int cellID = 1; cellID<= NCHANNELS; cellID++){
+ pedestals[cellID-1] = erc.getSingleSamplePedestal(event, cellID);
+ }
+ String line = time + "\t" + eventNum;
+ for(int i = 0; i< NCHANNELS; i++)
+ line += "\t" + prescale;
+ System.out.println(line);
+ }
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/ScalerAttenuationCalculator.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/ScalerAttenuationCalculator.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/ScalerAttenuationCalculator.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,255 @@
+package org.hps.users.spaul;
+
+import java.io.BufferedReader;
+import java.io.FileNotFoundException;
+import java.io.FileReader;
+import java.io.IOException;
+import java.text.SimpleDateFormat;
+import java.util.Date;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.TimeZone;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+
+import org.apache.commons.cli.CommandLine;
+import org.apache.commons.cli.CommandLineParser;
+import org.apache.commons.cli.HelpFormatter;
+import org.apache.commons.cli.Option;
+import org.apache.commons.cli.Options;
+import org.apache.commons.cli.ParseException;
+import org.apache.commons.cli.PosixParser;
+import org.apache.commons.csv.CSVFormat;
+import org.apache.commons.csv.CSVParser;
+import org.apache.commons.csv.CSVRecord;
+import org.hps.conditions.database.DatabaseConditionsManager;
+import org.hps.conditions.svt.SvtAlignmentConstant;
+import org.hps.conditions.svt.SvtBiasConstant;
+import org.hps.conditions.svt.SvtBiasConstant.SvtBiasConstantCollection;
+import org.hps.conditions.svt.SvtMotorPosition;
+import org.hps.conditions.svt.SvtMotorPosition.SvtMotorPositionCollection;
+import org.hps.conditions.svt.SvtTimingConstants;
+import org.hps.rundb.RunManager;
+import org.hps.users.meeg.SvtChargeIntegrator;
+
+/**
+ * Calculate the zeropoint value for each run from a myadump of scalerS2b
+ * @author Sebouh Paul
+ *
+ */
+public class ScalerAttenuationCalculator {
+ //after the value of scalerS2b drops below this value, (at a beam trip)
+ // take 10 samples, and use the 10th sample is the zero point value
+ static double zp_threshold = 500;
+
+ /**
+ *
+ * @param args the command line arguments (requires a CSV run/file log file
+ * and a MYA dump file.)
+ */
+ public static void main(String[] args) {
+
+ Options options = new Options();
+ options.addOption(new Option("z", true, "use zeropoint value from a file"));
+
+ final CommandLineParser parser = new PosixParser();
+ CommandLine cl = null;
+
+ try {
+ cl = parser.parse(options, args);
+ } catch (ParseException e) {
+ throw new RuntimeException("Cannot parse.", e);
+ }
+
+
+
+ Map<Integer, Long> runErrorMap = new HashMap<Integer, Long>();
+
+
+ if (cl.getArgs().length != 2) {
+ printUsage(options);
+ return;
+ }
+
+ List<CSVRecord> records = null;
+ try {
+ FileReader reader = new FileReader(cl.getArgs()[0]);
+ final CSVParser csvParser = new CSVParser(reader, CSVFormat.DEFAULT);
+
+ records = csvParser.getRecords();
+
+ csvParser.close();
+ } catch (FileNotFoundException ex) {
+ Logger.getLogger(ScalerAttenuationCalculator.class.getName()).log(Level.SEVERE, null, ex);
+ } catch (IOException ex) {
+ Logger.getLogger(ScalerAttenuationCalculator.class.getName()).log(Level.SEVERE, null, ex);
+ }
+
+ SimpleDateFormat dateFormat = new SimpleDateFormat("yyyy-MM-dd HH:mm:ss.SSS");
+ dateFormat.setTimeZone(TimeZone.getTimeZone("America/New_York"));
+
+ Map<Integer, Double> zeropointMap = new HashMap<Integer, Double>();
+ if(cl.hasOption("z")){
+ try {
+ BufferedReader br = new BufferedReader(new FileReader(cl.getOptionValue("z")));
+ String line;
+ System.err.println("zero point file header: " + br.readLine()); //discard the first line
+ double zeropoint = 0;
+ while ((line = br.readLine()) != null) {
+ String arr[] = line.split("[ \t]+");
+ int run = Integer.parseInt(arr[0]);
+
+ //if the zero point is unknown (ie, a very short run),
+ // use the zeropoint of the previous run
+ if(Double.parseDouble(arr[1]) != 0)
+ zeropoint = Double.parseDouble(arr[1]);
+
+ zeropointMap.put(run, zeropoint);
+ // System.out.format("%d %d\n", run, errorTime);
+ }
+ } catch (FileNotFoundException ex) {
+ Logger.getLogger(SvtChargeIntegrator.class.getName()).log(Level.SEVERE, null, ex);
+ } catch (IOException ex) {
+ Logger.getLogger(SvtChargeIntegrator.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+
+ try {
+ BufferedReader br = new BufferedReader(new FileReader(cl.getArgs()[1]));
+ String line;
+ System.err.println("myaData header: " + br.readLine()); //discard the first line
+ /*if (perRun) {
+ if (cl.hasOption("e")) {
+ System.out.println("run_num\tnominal_position\tnEvents\ttotalQ\ttotalQ_withbias\ttotalQ_atnom\ttotalQ_noerror\tgatedQ\tgatedQ_withbias\tgatedQ_atnom\tgatedQ_noerror\tgoodQ\tgoodQ_withbias\tgoodQ_atnom\tgoodQ_noerror");
+ } else {
+ System.out.println("run_num\tnominal_position\tnEvents\ttotalQ\ttotalQ_withbias\ttotalQ_atnom\tgatedQ\tgatedQ_withbias\tgatedQ_atnom\tgoodQ\tgoodQ_withbias\tgoodQ_atnom");
+ }
+ } else {
+ if (cl.hasOption("e")) {
+ System.out.println("run_num\tfile_num\tnominal_position\tnEvents\ttotalQ\ttotalQ_withbias\ttotalQ_atnom\ttotalQ_noerror\tgatedQ\tgatedQ_withbias\tgatedQ_atnom\tgatedQ_noerror\tgoodQ\tgoodQ_withbias\tgoodQ_atnom\tgoodQ_noerror");
+ } else {
+ System.out.println("run_num\tfile_num\tnominal_position\tnEvents\ttotalQ\ttotalQ_withbias\ttotalQ_atnom\tgatedQ\tgatedQ_withbias\tgatedQ_atnom\tgoodQ\tgoodQ_withbias\tgoodQ_atnom");
+ }
+ }*/
+ System.out.println("run_num\tzeropoint\tattenuation");
+
+
+ Date date = null;
+ Date lastDate;
+
+ for (CSVRecord record : records) {
+ int runNum = Integer.parseInt(record.get(0));
+
+
+
+
+ Date startDate, endDate;
+ long firstTime, lastTime;//Unix time from head bank
+ long firstTI, lastTI;//TI timestamp from TI bank
+
+
+ firstTime = Long.parseLong(record.get(7));
+ lastTime = Long.parseLong(record.get(8));
+ firstTI = Long.parseLong(record.get(10));
+ lastTI = Long.parseLong(record.get(11));
+
+
+
+
+ if (firstTime == 0 || lastTime == 0) {
+ continue;
+ }
+ startDate = new Date(firstTime * 1000);
+ endDate = new Date(lastTime * 1000);
+
+
+
+
+
+ br.mark(1000);
+
+
+
+
+ double attenuation_count = 0, attenuation_sum = 0;
+
+ //only change the zeropoint if there
+ double zeropoint = zeropointMap.get(runNum);
+
+ double crossing_time = 0;
+
+ double bpm_threshold = 40;
+
+ while ((line = br.readLine()) != null) {
+ String arr[] = line.split(" +");
+
+ if (arr.length != 5) {
+ throw new RuntimeException("this line is not correct.");
+ }
+ lastDate = date;
+ date = dateFormat.parse(arr[0] + " " + arr[1]);
+ if (date.before(startDate)) { //not in the file's time range yet; keep reading the file
+ continue;
+ }
+
+ double scaler, livetime, bpm;
+ if (arr[2].equals("<undefined>")) {
+ scaler = 0;
+ } else {
+ scaler = Double.parseDouble(arr[2]);
+ }
+
+ if (arr[3].equals("<undefined>")) {
+ bpm = 0;
+ } else {
+ bpm = Double.parseDouble(arr[3]);
+ }
+
+
+
+ if (lastDate != null) {
+ long dtStart = Math.max(startDate.getTime(), lastDate.getTime());
+ long dtEnd = Math.min(date.getTime(), endDate.getTime());
+ double dt = (dtEnd - dtStart) / 1000.0;
+
+ // average only the events when the bpm has been above a threshold for
+ // at least 15 seconds
+ if(bpm > bpm_threshold && crossing_time == 0)
+ crossing_time =date.getTime();
+ else if(bpm < bpm_threshold)
+ crossing_time = 0;
+ else if(crossing_time != 0 && date.getTime() > crossing_time+15000 && bpm > bpm_threshold){
+ attenuation_sum+= (scaler-zeropoint)/bpm;
+ attenuation_count ++;
+ }
+
+
+ }
+ if (date.after(endDate)) {//this is the last interval overlapping the file's time range; backtrack so this line will be read again for the next file
+ date = lastDate;
+ try {
+ br.reset();
+ } catch (IOException e) {
+ }
+ break;
+ }
+ br.mark(1000);
+ //current = nextCurrent;
+ }
+ System.out.printf("%d\t%f\t%f\n", runNum, zeropoint, attenuation_sum/attenuation_count);
+
+ }
+ } catch (Exception ex) {
+ Logger.getLogger(ScalerAttenuationCalculator.class.getName()).log(Level.SEVERE, null, ex);
+ }
+
+ }
+
+ private static void printUsage(Options options) {
+ HelpFormatter formatter = new HelpFormatter();
+ System.err.println("SvtChargeIntegrator <CSV log file> <MYA dump file>");
+ formatter.printHelp("Need to adhere to these options", options);
+ }
+}
+
Added: java/trunk/users/src/main/java/org/hps/users/spaul/ScalerZeropointCalculator.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/ScalerZeropointCalculator.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/ScalerZeropointCalculator.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,209 @@
+package org.hps.users.spaul;
+
+import java.io.BufferedReader;
+import java.io.FileNotFoundException;
+import java.io.FileReader;
+import java.io.IOException;
+import java.text.SimpleDateFormat;
+import java.util.Date;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.TimeZone;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+
+import org.apache.commons.cli.CommandLine;
+import org.apache.commons.cli.CommandLineParser;
+import org.apache.commons.cli.HelpFormatter;
+import org.apache.commons.cli.Option;
+import org.apache.commons.cli.Options;
+import org.apache.commons.cli.ParseException;
+import org.apache.commons.cli.PosixParser;
+import org.apache.commons.csv.CSVFormat;
+import org.apache.commons.csv.CSVParser;
+import org.apache.commons.csv.CSVRecord;
+import org.hps.conditions.database.DatabaseConditionsManager;
+import org.hps.conditions.svt.SvtAlignmentConstant;
+import org.hps.conditions.svt.SvtBiasConstant;
+import org.hps.conditions.svt.SvtBiasConstant.SvtBiasConstantCollection;
+import org.hps.conditions.svt.SvtMotorPosition;
+import org.hps.conditions.svt.SvtMotorPosition.SvtMotorPositionCollection;
+import org.hps.conditions.svt.SvtTimingConstants;
+import org.hps.rundb.RunManager;
+
+/**
+ * Calculate the zeropoint value for each run from a myadump of scalerS2b
+ * @author Sebouh Paul
+ *
+ */
+public class ScalerZeropointCalculator {
+ //after the value of scalerS2b drops below this value, (at a beam trip)
+ // take 10 samples, and use the 10th sample is the zero point value
+ static double threshold = 500;
+
+ /**
+ *
+ * @param args the command line arguments (requires a CSV run/file log file
+ * and a MYA dump file.)
+ */
+ public static void main(String[] args) {
+
+ Options options = new Options();
+
+ final CommandLineParser parser = new PosixParser();
+ CommandLine cl = null;
+ try {
+ cl = parser.parse(options, args);
+ } catch (ParseException e) {
+ throw new RuntimeException("Cannot parse.", e);
+ }
+
+
+
+ Map<Integer, Long> runErrorMap = new HashMap<Integer, Long>();
+
+
+ if (cl.getArgs().length != 2) {
+ printUsage(options);
+ return;
+ }
+
+ List<CSVRecord> records = null;
+ try {
+ FileReader reader = new FileReader(cl.getArgs()[0]);
+ final CSVParser csvParser = new CSVParser(reader, CSVFormat.DEFAULT);
+
+ records = csvParser.getRecords();
+
+ csvParser.close();
+ } catch (FileNotFoundException ex) {
+ Logger.getLogger(ScalerZeropointCalculator.class.getName()).log(Level.SEVERE, null, ex);
+ } catch (IOException ex) {
+ Logger.getLogger(ScalerZeropointCalculator.class.getName()).log(Level.SEVERE, null, ex);
+ }
+
+ SimpleDateFormat dateFormat = new SimpleDateFormat("yyyy-MM-dd HH:mm:ss.SSS");
+ dateFormat.setTimeZone(TimeZone.getTimeZone("America/New_York"));
+
+ try {
+ BufferedReader br = new BufferedReader(new FileReader(cl.getArgs()[1]));
+ String line;
+ System.err.println("myaData header: " + br.readLine()); //discard the first line
+ /*if (perRun) {
+ if (cl.hasOption("e")) {
+ System.out.println("run_num\tnominal_position\tnEvents\ttotalQ\ttotalQ_withbias\ttotalQ_atnom\ttotalQ_noerror\tgatedQ\tgatedQ_withbias\tgatedQ_atnom\tgatedQ_noerror\tgoodQ\tgoodQ_withbias\tgoodQ_atnom\tgoodQ_noerror");
+ } else {
+ System.out.println("run_num\tnominal_position\tnEvents\ttotalQ\ttotalQ_withbias\ttotalQ_atnom\tgatedQ\tgatedQ_withbias\tgatedQ_atnom\tgoodQ\tgoodQ_withbias\tgoodQ_atnom");
+ }
+ } else {
+ if (cl.hasOption("e")) {
+ System.out.println("run_num\tfile_num\tnominal_position\tnEvents\ttotalQ\ttotalQ_withbias\ttotalQ_atnom\ttotalQ_noerror\tgatedQ\tgatedQ_withbias\tgatedQ_atnom\tgatedQ_noerror\tgoodQ\tgoodQ_withbias\tgoodQ_atnom\tgoodQ_noerror");
+ } else {
+ System.out.println("run_num\tfile_num\tnominal_position\tnEvents\ttotalQ\ttotalQ_withbias\ttotalQ_atnom\tgatedQ\tgatedQ_withbias\tgatedQ_atnom\tgoodQ\tgoodQ_withbias\tgoodQ_atnom");
+ }
+ }*/
+ System.out.println("run_num\tzeropoint");
+
+
+ Date date = null;
+ Date lastDate;
+
+ for (CSVRecord record : records) {
+ int runNum = Integer.parseInt(record.get(0));
+
+
+
+
+ Date startDate, endDate;
+ long firstTime, lastTime;//Unix time from head bank
+ long firstTI, lastTI;//TI timestamp from TI bank
+
+
+ firstTime = Long.parseLong(record.get(7));
+ lastTime = Long.parseLong(record.get(8));
+ firstTI = Long.parseLong(record.get(10));
+ lastTI = Long.parseLong(record.get(11));
+
+
+
+
+ if (firstTime == 0 || lastTime == 0) {
+ continue;
+ }
+ startDate = new Date(firstTime * 1000);
+ endDate = new Date(lastTime * 1000);
+
+
+
+
+
+ br.mark(1000);
+
+
+ //find the 10th sample after the first trip begins (10 is an arbitrary number)
+
+ double zeroPoint = 0;
+ int count = 0;
+ while ((line = br.readLine()) != null) {
+ String arr[] = line.split(" +");
+
+ if (arr.length != 4) {
+ throw new RuntimeException("this line is not correct.");
+ }
+ lastDate = date;
+ date = dateFormat.parse(arr[0] + " " + arr[1]);
+ if (date.before(startDate)) { //not in the file's time range yet; keep reading the file
+ continue;
+ }
+
+ double current, livetime;
+ if (arr[2].equals("<undefined>")) {
+ current = 0;
+ } else {
+ current = Double.parseDouble(arr[2]);
+ }
+
+
+
+
+ if (lastDate != null) {
+ long dtStart = Math.max(startDate.getTime(), lastDate.getTime());
+ long dtEnd = Math.min(date.getTime(), endDate.getTime());
+ double dt = (dtEnd - dtStart) / 1000.0;
+
+ if(current < threshold)
+ count++;
+ if(count == 10){
+ zeroPoint = current;
+ break;
+ }
+
+ }
+ if (date.after(endDate)) {//this is the last interval overlapping the file's time range; backtrack so this line will be read again for the next file
+ date = lastDate;
+ try {
+ br.reset();
+ } catch (IOException e) {
+ }
+ break;
+ }
+ br.mark(1000);
+ //current = nextCurrent;
+ }
+ System.out.printf("%d\t%f\n", runNum, zeroPoint);
+
+ }
+ } catch (Exception ex) {
+ Logger.getLogger(ScalerZeropointCalculator.class.getName()).log(Level.SEVERE, null, ex);
+ }
+
+ }
+
+ private static void printUsage(Options options) {
+ HelpFormatter formatter = new HelpFormatter();
+ System.err.println("SvtChargeIntegrator <CSV log file> <MYA dump file>");
+ formatter.printHelp("Need to adhere to these options", options);
+ }
+}
+
Added: java/trunk/users/src/main/java/org/hps/users/spaul/TestAIDA2ROOT.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/TestAIDA2ROOT.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/TestAIDA2ROOT.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,15 @@
+package org.hps.users.spaul;
+
+import java.io.IOException;
+
+import hep.aida.IHistogram1D;
+
+import org.lcsim.util.aida.AIDA;
+
+public class TestAIDA2ROOT {
+ public static void main(String arg[]) throws IOException{
+ AIDA aida = AIDA.defaultInstance();
+ IHistogram1D h = aida.histogram1D("/folder/title;xaxis;yaxis", 100, 0, 10);
+ aida.saveAs("/Users/sebouhpaul/test.root");
+ }
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/bh/BumpHunt.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/bh/BumpHunt.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/bh/BumpHunt.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,104 @@
+package org.hps.users.spaul.bh;
+
+import hep.aida.IAnalysisFactory;
+
+import java.io.FileNotFoundException;
+import java.io.PrintWriter;
+import java.util.ArrayList;
+import java.util.Arrays;
+
+import org.apache.commons.math3.fitting.GaussianCurveFitter;
+import org.apache.commons.math3.fitting.PolynomialCurveFitter;
+import org.apache.commons.math3.fitting.SimpleCurveFitter;
+import org.apache.commons.math3.distribution.ChiSquaredDistribution;
+
+public class BumpHunt {
+ static int degree;
+ static int iter;
+ public static void main(String arg[]) throws FileNotFoundException{
+
+ Histogram h = new Histogram(arg[0]);
+ WindowSet ws = new WindowSet(arg[1]);
+ degree = Integer.parseInt(arg[2]);
+ String outfile = null;
+ if(arg.length > 3)
+ outfile = arg[3];
+ results = new ArrayList<WindowResult>();
+ for(int i = 0; i< ws.N; i++){
+ try{
+ results.add(process(h, ws, i));
+
+ } catch(Exception e){
+
+ }
+ }
+
+ if(outfile != null){
+ PrintWriter pw = new PrintWriter(outfile);
+ for(WindowResult wr : results){
+ pw.printf("%f %f %f %f %f", wr.min, wr.max, wr.mean, wr.yield, wr.pull);
+ }
+ }
+ }
+ static ArrayList<WindowResult> results;
+
+ public static WindowResult process(Histogram h, WindowSet ws, int i) {
+ double min = ws.minMass[i];
+ double max = ws.maxMass[i];
+ double res = ws.resolution[i];
+
+ Polynomial p = Fitters.polyfit(h, min, max, degree);
+
+ //PolynomialPlusGaussian pplusg = Fitters.pplusgfit(h, p, min, max, res);
+ PolynomialPlusGaussian pplusg = Fitters.boxFit(h, p, min, max, res);
+
+ pplusg = Fitters.pplusgfit(h, p, min, max, res, pplusg.N/.80);
+
+ System.out.println(min + " " + max);
+ System.out.println(Arrays.toString(p.p));
+
+ //test statistic for polynomial-only fit.
+ double ll1 = LogLikelihood.getLogLikelihood(h, min, max, p);
+ System.out.println("poly LL: " + ll1);
+ System.out.println(Arrays.toString(pplusg.p) + " " + pplusg.mean + " " + pplusg.sigma + " " + pplusg.N);
+ double ll2 = LogLikelihood.getLogLikelihood(h, min, max, pplusg);
+ System.out.println("gauss poly LL: " + ll2);
+
+ double p_value = new ChiSquaredDistribution(1).cumulativeProbability(ll1-ll2);
+
+ System.out.println();
+
+
+ double yield = pplusg.N;
+ double yieldError = getYieldError(p, h, res, pplusg.mean);
+
+ double pull = yield/yieldError;
+
+ WindowResult wr = new WindowResult();
+ wr.pplusg = pplusg;
+ wr.poly = p;
+ wr.max = max;
+ wr.min = min;
+ wr.mean = pplusg.mean;
+ wr.yield = yield;
+ wr.pull = pull;
+ wr.p_value = p_value;
+
+ return wr;
+ }
+
+
+ final static double twoSqrtPi = 2*Math.sqrt(Math.PI);
+ private static double getYieldError(Polynomial p, Histogram h, double res,
+ double mean) {
+ double bw = h.maxMass[0]-h.minMass[0];
+
+ return Math.sqrt(p.get(mean)*res/bw)*.345;
+ }
+
+
+
+
+
+
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/bh/Fitters.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/bh/Fitters.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/bh/Fitters.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,139 @@
+package org.hps.users.spaul.bh;
+
+import java.util.ArrayList;
+
+import org.apache.commons.math3.analysis.ParametricUnivariateFunction;
+import org.apache.commons.math3.analysis.UnivariateFunction;
+import org.apache.commons.math3.fitting.PolynomialCurveFitter;
+import org.apache.commons.math3.fitting.SimpleCurveFitter;
+import org.apache.commons.math3.fitting.WeightedObservedPoint;
+
+public class Fitters {
+ static Polynomial polyfit(Histogram h, double min, double max, int degree){
+ PolynomialCurveFitter pcf = PolynomialCurveFitter.create(degree);
+
+ ArrayList<WeightedObservedPoint> wops = new ArrayList();
+ for(int j = 0; j< h.Nbins; j++){
+ if(h.minMass[j]>= min && h.maxMass[j]<= max){
+ double w = 1/h.h[j];
+ WeightedObservedPoint wop = new WeightedObservedPoint(w, (h.minMass[j]+h.maxMass[j])/2., h.h[j]);
+ wops.add(wop);
+ }
+ }
+ double fit[] = pcf.fit(wops);
+
+ Polynomial p = new Polynomial(fit.length);
+ p.p = fit;
+ return p;
+ }
+ static PolynomialPlusGaussian pplusgfit(Histogram h, Polynomial p, double min, double max, double res, double yield_guess){
+
+ double start[] = new double[p.k+3];
+ start[start.length-3] = (min+max)/2;
+ start[start.length-2] = res;
+
+ start[start.length-1] = yield_guess;
+
+ SimpleCurveFitter scf = SimpleCurveFitter.create(new ParametricGpG(), start);
+ ArrayList<WeightedObservedPoint> wops = new ArrayList();
+ for(int j = 0; j< h.Nbins; j++){
+ if(h.minMass[j]>= min && h.maxMass[j]<= max){
+ double w = 1/h.h[j];
+ WeightedObservedPoint wop = new WeightedObservedPoint(w, (h.minMass[j]+h.maxMass[j])/2., h.h[j]);
+ wops.add(wop);
+ }
+ }
+ double[] result = scf.fit(wops);
+ PolynomialPlusGaussian ppg = new PolynomialPlusGaussian(p.k);
+ for(int i = 0; i< p.k; i++)
+ ppg.p[i] = result[i];
+ ppg.mean = result[p.k];
+ ppg.sigma = result[p.k+1];
+ ppg.N = result[p.k+2];
+
+ return ppg;
+ }
+
+ static PolynomialPlusGaussian boxFit(Histogram h, Polynomial p, double min, double max, double res){
+
+
+ double boxmin = (max+min)/2. - 1.25*res;
+ double boxmax = (max+min)/2. + 1.25*res;
+ double N = 0;
+ for(int j = 0; j< h.Nbins; j++){
+
+ if(h.minMass[j]>= boxmin && h.maxMass[j]<= boxmax){
+ double binMean = (h.maxMass[j]+ h.minMass[j])/2.;
+ N+= h.h[j]-p.get(binMean);
+ }
+ }
+
+ //N/=.6;
+
+
+ PolynomialPlusGaussian ppg = new PolynomialPlusGaussian(p.k);
+ for(int i = 0; i< p.k; i++)
+ ppg.p[i] = p.p[i];
+ ppg.mean = (max+min)/2.;
+ ppg.sigma = res;
+ ppg.N = N;
+
+ return ppg;
+ }
+
+
+
+
+
+ static class ParametricGpG implements
+ ParametricUnivariateFunction{
+
+
+
+ @Override
+ public double[] gradient(double arg0, double... arg1) {
+
+ double[] ret = new double [arg1.length];
+ double powM = 1;
+ for(int i = 0; i<arg1.length-3; i++){
+ ret[i] = powM;
+ powM *= arg0;
+
+ }
+ double mean = arg1[arg1.length-3];
+ double sigma = arg1[arg1.length-2];
+ double N = arg1[arg1.length-1];
+
+ ret[arg1.length-3] = 0;
+ //-N*(mean-arg0)/Math.pow(sigma,2)*Math.exp(-Math.pow((arg0-mean)/sigma,2)/2)
+ // /(sqrt2pi*sigma);
+
+ ret[arg1.length-2] = 0;
+ //N*Math.exp(-Math.pow((arg0-mean)/sigma,2)/2)/(sqrt2pi)
+ //*(-1/(sigma*sigma) + (arg0-mean)*(arg0-mean)/Math.pow(sigma,4));
+
+ ret[arg1.length-1] = Math.exp(-Math.pow((arg0-mean)/sigma,2)/2)
+ /(sqrt2pi*sigma);
+
+ return ret;
+ }
+
+ @Override
+ public double value(double arg0, double... arg1) {
+
+ double y = 0;
+ double powM = 1;
+ for(int i = 0; i<arg1.length-3; i++){
+ y += powM*arg1[i];
+ powM *= arg0;
+ }
+ double mean = arg1[arg1.length-3];
+ double sigma = arg1[arg1.length-2];
+ double N = arg1[arg1.length-1];
+ return y + N*Math.exp(-Math.pow((arg0-mean)/sigma,2)/2)
+ /(sqrt2pi*sigma);
+ }
+ };
+
+ static final double sqrt2pi = Math.sqrt(2*Math.PI);
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/bh/GenerateFakeStuff.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/bh/GenerateFakeStuff.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/bh/GenerateFakeStuff.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,50 @@
+package org.hps.users.spaul.bh;
+
+import java.io.FileNotFoundException;
+import java.io.PrintWriter;
+import java.util.Random;
+
+public class GenerateFakeStuff {
+ static Random random = new Random();
+ static PolynomialPlusGaussian pplusg;
+ public static void main(String arg[]) throws FileNotFoundException{
+
+ PrintWriter pw1 = new PrintWriter(arg[0]);
+ double signalSize = Double.parseDouble(arg[1]);
+ pw1.println(1000);
+
+
+ double min = 0;
+ double max = 100;
+ pplusg = new PolynomialPlusGaussian(3);
+ //pplusg.p = new double[]{10000, -300, 20};
+
+ pplusg.p = randomPoly(min, max);
+ pplusg.mean = 50;
+ pplusg.sigma = 1;
+ pplusg.N = signalSize*Math.sqrt(pplusg.get(pplusg.mean));
+
+ System.out.println(pplusg);
+ for(int i = 0; i< 1000; i++){
+ double x = i/10.;
+ double y = pplusg.get(x);
+
+ pw1.println(i/10. + " " + (i+1)/10. + " " + Math.round(y + Math.sqrt(y)*random.nextGaussian()));
+ }
+ pw1.close();
+
+ PrintWriter pw2 = new PrintWriter("windows.txt");
+ pw2.println(90);
+ for(int i = 0; i< 90; i++){
+
+ pw2.println(i + " " + (i+10) + " " + 1);
+
+ }
+ pw2.close();
+ }
+ private static double[] randomPoly(double min, double max) {
+
+ double[] p = new double[]{10000*(.7+.6*random.nextDouble()), -30*random.nextGaussian(), .2*random.nextGaussian()};
+ return p;
+ }
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/bh/GetEpsilon2.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/bh/GetEpsilon2.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/bh/GetEpsilon2.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,13 @@
+package org.hps.users.spaul.bh;
+
+import java.io.FileNotFoundException;
+
+public class GetEpsilon2 {
+ public static void main(String arg[]) throws FileNotFoundException{
+ String radHistFile = arg[0];
+ String bumpHuntResults = arg[1];
+ Histogram radHist = new Histogram(radHistFile);
+
+
+ }
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/bh/GetPullDistribution.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/bh/GetPullDistribution.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/bh/GetPullDistribution.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,40 @@
+package org.hps.users.spaul.bh;
+
+import java.io.FileNotFoundException;
+
+import hep.aida.*;
+
+public class GetPullDistribution {
+ public static void main(String arg[]) throws FileNotFoundException{
+ IAnalysisFactory af = IAnalysisFactory.create();
+ ITree tree = af.createTreeFactory().create();
+ IHistogramFactory hf = af.createHistogramFactory(tree);
+
+ IHistogram1D h1 = hf.createHistogram1D("pull", 100, -10, 10);
+ IPlotter p = af.createPlotterFactory().create();
+ p.region(0).plot(h1);
+ p.show();
+
+ IHistogram1D h2 = hf.createHistogram1D("yield", 100, -1000, 1000);
+ p = af.createPlotterFactory().create();
+ p.region(0).plot(h2);
+ p.show();
+
+ IHistogram1D h3 = hf.createHistogram1D("p_value", 100, 0, 1);
+ p = af.createPlotterFactory().create();
+ p.region(0).plot(h3);
+ p.show();
+
+ for(int i = 0; i<100000; i++){
+ GenerateFakeStuff.main(new String[]{"hist.txt", "0"});
+ BumpHunt.main(new String[]{"hist.txt", "windows.txt", "4"});
+ for(WindowResult wr : BumpHunt.results){
+ //if(wr.p_value > .999){
+ h1.fill(wr.pull);
+ h2.fill(wr.yield);
+ //}
+ h3.fill(wr.p_value);
+ }
+ }
+ }
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/bh/GetYieldPerSignal.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/bh/GetYieldPerSignal.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/bh/GetYieldPerSignal.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,82 @@
+package org.hps.users.spaul.bh;
+
+import java.io.FileNotFoundException;
+import java.util.Random;
+
+import org.hps.users.spaul.StyleUtil;
+
+import hep.aida.*;
+
+public class GetYieldPerSignal {
+ public static void main(String arg[]) throws FileNotFoundException{
+ IAnalysisFactory af = IAnalysisFactory.create();
+ ITree tree = af.createTreeFactory().create();
+ IHistogramFactory hf = af.createHistogramFactory(tree);
+
+ IHistogram1D h1 = hf.createHistogram1D("pull", 100, -10, 10);
+ IPlotter p = af.createPlotterFactory().create();
+ p.region(0).plot(h1);
+ p.show();
+
+ IHistogram1D h2 = hf.createHistogram1D("yield", 100, -5000, 5000);
+ p = af.createPlotterFactory().create();
+ p.region(0).plot(h2);
+ p.show();
+
+ IHistogram1D h3 = hf.createHistogram1D("p_value", 100, 0, 1);
+ p = af.createPlotterFactory().create();
+ p.region(0).plot(h3);
+ p.show();
+
+ IHistogram2D h4 = hf.createHistogram2D("yield vs signal", 100, -1000, 4000, 100, -1000, 4000);
+
+ p = af.createPlotterFactory().create();
+ p.region(0).plot(h4);
+ StyleUtil.stylize(p.region(0), "yield", "signal");
+ p.show();
+
+
+ WindowSet ws = new WindowSet("windows.txt");
+ BumpHunt.degree = 2;
+
+ Random r = new Random();
+
+ //variables for a linear fit
+ double sx = 0, sy = 0,sy2 = 0, sx2 = 0, sxy = 0, n = 0;
+ for(int i = 0; i<10000; i++){
+ GenerateFakeStuff.main(new String[]{"hist.txt", "5"});
+
+ Histogram h = new Histogram("hist.txt");
+
+
+ WindowResult wr = BumpHunt.process(h, ws, 45);
+ h1.fill(wr.pull);
+ h2.fill(wr.yield);
+ //}
+ h3.fill(wr.p_value);
+ //}
+ double s = r.nextDouble()*10;
+ GenerateFakeStuff.main(new String[]{"hist.txt", Double.toString(s)});
+ h = new Histogram("hist.txt");
+
+
+ wr = BumpHunt.process(h, ws, 45);
+ h4.fill(GenerateFakeStuff.pplusg.N, wr.yield);
+ double x = GenerateFakeStuff.pplusg.N;
+ double y = wr.yield;
+ sy += y;
+ sy2 += y*y;
+ sx +=x;
+ sx2+= x*x;
+ sxy += x*y;
+ n++;
+ }
+ double m = (sxy*n-sy*sx)/(sx2*n-sx*sx);
+ double b = sy/n-m*sx/n;
+ System.out.println("m = " + m);
+ System.out.println("b = " + b);
+
+
+
+ }
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/bh/GoodnessOfFitAccumulator.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/bh/GoodnessOfFitAccumulator.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/bh/GoodnessOfFitAccumulator.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,60 @@
+package org.hps.users.spaul.bh;
+
+public class GoodnessOfFitAccumulator{
+ private double logLikelinessSum = 0;
+ private double chi2Sum = 0;
+ /**
+ * Accumulate a data point
+ * @param h the measured value
+ * @param f the fit value
+ */
+ public void accumulate(double h, double f) {
+ logLikelinessSum+=2*(f-h+h*Math.log(h/f));
+ chi2Sum+= (f-h)*(f-h)/h;
+ total_dof ++;
+ }
+ /**
+ *@return -2*log(L)
+ */
+ public double getLogLikeliness() {
+ return logLikelinessSum;
+ }
+ /**
+ *@return Pearson's chi^2
+ */
+ public double getChi2() {
+ return chi2Sum;
+ }
+ /**
+ *@return the number of degrees of freedom
+ */
+ public double getDOF(){
+ return total_dof;
+ }
+ /**
+ * @return the number of parameters in the fit
+ */
+ public double getFitParameterCount(){
+ return fitParamCount;
+ }
+
+ /**
+ * resets the sums
+ */
+ public void reset(){
+ logLikelinessSum =0;
+ chi2Sum = 0;
+ total_dof = -fitParamCount;
+ }
+ private int fitParamCount;
+ private int total_dof;
+ /**
+ * Constructs an instance of GoodnessOfFit.
+ * @param fitParamCount the number of parameters in the fit
+ */
+ public GoodnessOfFitAccumulator(int fitParamCount){
+ this.fitParamCount = fitParamCount;
+ this.total_dof = -fitParamCount;
+ }
+
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/bh/Histogram.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/bh/Histogram.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/bh/Histogram.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,26 @@
+package org.hps.users.spaul.bh;
+
+import java.io.File;
+import java.io.FileNotFoundException;
+import java.util.Scanner;
+
+public class Histogram {
+ int Nbins;
+ Histogram(String s) throws FileNotFoundException{
+ this(new File(s));
+ }
+ double[] h, minMass, maxMass;
+ public Histogram(File file) throws FileNotFoundException {
+ Scanner scanner = new Scanner(file);
+ Nbins = scanner.nextInt();
+ minMass = new double[Nbins];
+ maxMass = new double[Nbins];
+ h = new double[Nbins];
+ for(int i = 0; i< Nbins; i++){
+ minMass[i] = scanner.nextDouble();
+ maxMass[i] = scanner.nextDouble();
+ h[i] = scanner.nextDouble();
+ }
+ }
+
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/bh/IFitResult.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/bh/IFitResult.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/bh/IFitResult.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,6 @@
+package org.hps.users.spaul.bh;
+
+public interface IFitResult {
+ double get(double m);
+ int paramCount();
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/bh/LogLikelihood.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/bh/LogLikelihood.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/bh/LogLikelihood.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,15 @@
+package org.hps.users.spaul.bh;
+
+import org.apache.commons.math3.fitting.WeightedObservedPoint;
+
+public class LogLikelihood {
+ static double getLogLikelihood(Histogram h, double min, double max, IFitResult result){
+ GoodnessOfFitAccumulator gofa = new GoodnessOfFitAccumulator(result.paramCount());
+ for(int j = 0; j< h.Nbins; j++){
+ if(h.minMass[j]>= min && h.maxMass[j]<= max){
+ gofa.accumulate(h.h[j], result.get((h.minMass[j]+h.maxMass[j])/2.));
+ }
+ }
+ return gofa.getLogLikeliness();
+ }
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/bh/Polynomial.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/bh/Polynomial.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/bh/Polynomial.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,20 @@
+package org.hps.users.spaul.bh;
+
+public class Polynomial implements IFitResult {
+ double p[];
+ int k;
+ Polynomial(int k){
+ this.p = new double[k];
+ this.k = k;
+ }
+ public double get(double m){
+ double y = 0;
+ double powM = 1;
+ for(int i = 0; i<k; i++){
+ y += powM*p[i];
+ powM *= m;
+ }
+ return y;
+ }
+ public int paramCount(){return k;}
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/bh/PolynomialPlusGaussian.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/bh/PolynomialPlusGaussian.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/bh/PolynomialPlusGaussian.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,38 @@
+package org.hps.users.spaul.bh;
+
+import java.util.Arrays;
+
+public class PolynomialPlusGaussian implements IFitResult{
+
+ /**
+ * The square root of two pi
+ */
+ static final double sqrt2pi = Math.sqrt(2*Math.PI);
+
+
+ //parameters of the fit
+ double p[];
+ int k;
+ double mean, sigma, N;
+
+ PolynomialPlusGaussian(int k){
+ this.p = new double[k];
+ this.k = k;
+ }
+ public double get(double m){
+ double y = 0;
+ double powM = 1;
+ for(int i = 0; i<k; i++){
+ y += powM*p[i];
+ powM *= m;
+ }
+ return y + N*Math.exp(-Math.pow((m-mean)/sigma,2)/2)
+ /(sqrt2pi*sigma);
+ }
+ public int paramCount(){
+ return k + 3;
+ }
+ public String toString(){
+ return "p + g " + Arrays.toString(p) + " " + mean + " " + sigma + " " + N;
+ }
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/bh/WindowResult.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/bh/WindowResult.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/bh/WindowResult.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,14 @@
+package org.hps.users.spaul.bh;
+
+public class WindowResult {
+ double min, max;
+
+ IFitResult poly;
+ IFitResult pplusg;
+
+ double mean;
+ double yield;
+ double pull;
+
+ double p_value;
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/bh/WindowSet.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/bh/WindowSet.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/bh/WindowSet.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,40 @@
+package org.hps.users.spaul.bh;
+
+import java.io.File;
+import java.io.FileNotFoundException;
+import java.util.Scanner;
+
+public class WindowSet {
+ int N;
+
+ double minMass[];
+ double maxMass[];
+ double resolution[];
+
+ WindowSet(String s) throws FileNotFoundException{
+ this(new File(s));
+ }
+ public WindowSet(File file) throws FileNotFoundException {
+ Scanner scanner = new Scanner(file);
+ N = scanner.nextInt();
+ minMass = new double[N];
+ maxMass = new double[N];
+ resolution = new double[N];
+ for(int i = 0; i< N; i++){
+ minMass[i] = scanner.nextDouble();
+ maxMass[i] = scanner.nextDouble();
+ resolution[i] = scanner.nextDouble();
+ }
+ printToLatexTable();
+
+ }
+ void printToLatexTable(){
+ System.out.println("window # & center & $\\sigma_M$ & bottom & top & width \\\\");
+ for(int i = 0; i<N; i++){
+ System.out.printf("%d & %.0f & %.0f & %.0f & %.0f & %.0f \\\\\n",
+ i+1, (minMass[i]+maxMass[i])/2., resolution[i],
+ minMass[i], maxMass[i], maxMass[i]-minMass[i]);
+ }
+ }
+
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/Cuts.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/Cuts.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/Cuts.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,64 @@
+package org.hps.users.spaul.bh.test;
+
+import org.lcsim.event.Vertex;
+
+public class Cuts {
+
+ static double xtarget = 0; //NABO
+
+ static double xmax =2, ymax = 1,
+ vchi2max = 10,
+ v0pzmin = 1.2, v0pzmax = 2.75,
+ minP = .1, maxP = 2.0,
+ tchi2max = 50;
+
+ public static boolean passesCuts(Vertex v) {
+ double x = v.getPosition().x();
+ double y = v.getPosition().y();
+ double z = v.getPosition().z();
+ if(Math.abs(x-xtarget)> 2 || Math.abs(y) > 1)
+ return false;
+ System.out.println("passed position test");
+ if(v.getChi2()> vchi2max)
+ return false;
+ System.out.println("passed vertex chi2 test");
+
+ double px = v.getAssociatedParticle().getMomentum().x();
+ double py = v.getAssociatedParticle().getMomentum().y();
+ double pz = v.getAssociatedParticle().getMomentum().z();
+
+ if(pz < v0pzmin || pz > v0pzmax)
+ return false;
+
+ System.out.println("passed pztot test");
+
+ double px1 = v.getAssociatedParticle().getParticles().get(0).getMomentum().x();
+ double py1 = v.getAssociatedParticle().getParticles().get(0).getMomentum().y();
+ double pz1 = v.getAssociatedParticle().getParticles().get(0).getMomentum().z();
+
+ if(pz1 < minP || pz1> maxP)
+ return false;
+ System.out.println("passed pz1 test");
+
+ double px2 = v.getAssociatedParticle().getParticles().get(1).getMomentum().x();
+ double py2 = v.getAssociatedParticle().getParticles().get(1).getMomentum().y();
+ double pz2 = v.getAssociatedParticle().getParticles().get(1).getMomentum().z();
+
+ if(pz2 < minP || pz2> maxP)
+ return false;
+ System.out.println("passed pz2 test");
+
+ if(v.getAssociatedParticle().getParticles().get(0).getTracks().size() == 0
+ || v.getAssociatedParticle().getParticles().get(1).getTracks().size() == 0)
+ return false;
+
+ double chi2_1 = v.getAssociatedParticle().getParticles().get(0).getTracks().get(0).getChi2();
+
+ double chi2_2 = v.getAssociatedParticle().getParticles().get(1).getTracks().get(0).getChi2();
+ if(chi2_1 > tchi2max || chi2_2 > tchi2max)
+ return false;
+
+ return true;
+ }
+
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/GenerateFakeResolutions.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/GenerateFakeResolutions.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/GenerateFakeResolutions.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,27 @@
+package org.hps.users.spaul.bh.test;
+
+import java.io.IOException;
+import java.util.Random;
+
+import hep.aida.*;
+
+public class GenerateFakeResolutions {
+ public static void main(String arg[]) throws IllegalArgumentException, IOException{
+ IAnalysisFactory af = IAnalysisFactory.create();
+ Random r = new Random();
+ int masses[] = {25, 50, 75, 100, 150, 200, 250};
+ for(int m : masses){
+ String filename = m + "_fake.aida";
+ ITree tree = af.createTreeFactory().create(filename,"xml",false,true);
+ IHistogramFactory hf = af.createHistogramFactory(tree);
+ IHistogram1D h = hf.createHistogram1D("recon mass", 300, 0, 300);
+ double dm$m = .05;
+ for(int i = 0; i< 10000; i++){
+ h.fill(m + m*dm$m*r.nextGaussian());
+ }
+ tree.commit();
+ tree.close();
+ System.out.println("wrote file " + filename);
+ }
+ }
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/InvariantMassCalculator.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/InvariantMassCalculator.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/InvariantMassCalculator.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,87 @@
+package org.hps.users.spaul.bh.test;
+
+import hep.physics.vec.Hep3Vector;
+
+import java.util.List;
+
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.Vertex;
+
+public class InvariantMassCalculator {
+ static double Ebeam = 2.2;
+
+ static String vertexCollectionName = "TargetConstrainedV0Vertices";
+
+ static public double getMass(Vertex v){
+
+
+ return getMassRecon(v.getAssociatedParticle().getParticles().get(0),
+ v.getAssociatedParticle().getParticles().get(1));
+
+
+ }
+ //special case: we have the recoiled particle which is particle #3
+ static double getMassRecon(ReconstructedParticle part1, ReconstructedParticle part2, ReconstructedParticle part3) {
+
+ Hep3Vector p1 = part1.getMomentum();
+ Hep3Vector p2 = part2.getMomentum();
+
+ double E1 = getCombinedEnergy(part1);
+ double E2 = getCombinedEnergy(part2);
+ double E3 = getCombinedEnergy(part3);
+ if(E1 > E2){
+ E1 = Ebeam - E2 - E3;
+ }
+ else if(E1 < E2){
+ E2 = Ebeam - E1 - E3;
+ }
+
+ double denom = Math.sqrt(p1.magnitudeSquared()*p2.magnitudeSquared());
+ double cos = (p1.x()*p2.x()+p1.y()*p2.y()+p1.z()*p2.z())
+ /denom;
+
+ double sin2 =
+ (Math.pow(p1.y()*p2.z()-p1.z()*p2.y(),2) +
+ Math.pow(p1.z()*p2.x()-p1.x()*p2.z(),2) +
+ Math.pow(p1.x()*p2.y()-p1.y()*p2.x(),2)
+ )/denom*denom;
+
+ //this is equivalent to 2*sin(theta/2)*sqrt(E1*E2);
+ return Math.sqrt(2*E1*E2*sin2/(1+cos));
+
+ //return Math.sqrt(2*E1*E2*(1-cos));
+
+ }
+
+
+ static double getMassRecon(ReconstructedParticle part1, ReconstructedParticle part2) {
+
+ Hep3Vector p1 = part1.getMomentum();
+ Hep3Vector p2 = part2.getMomentum();
+
+ double E1 = getCombinedEnergy(part1);
+ double E2 = getCombinedEnergy(part2);
+
+ double denom = Math.sqrt(p1.magnitudeSquared()*p2.magnitudeSquared());
+ double cos = (p1.x()*p2.x()+p1.y()*p2.y()+p1.z()*p2.z())
+ /denom;
+ double sin2 =
+ (Math.pow(p1.y()*p2.z()-p1.z()*p2.y(),2) +
+ Math.pow(p1.z()*p2.x()-p1.x()*p2.z(),2) +
+ Math.pow(p1.x()*p2.y()-p1.y()*p2.x(),2)
+ )/denom*denom;
+
+ //this is equivalent to 2*sin(theta/2)*sqrt(E1*E2);
+ //return Math.sqrt(2*E1*E2*sin2/(1+cos));
+
+ return Math.sqrt(2*E1*E2*(1-cos));
+
+ }
+
+ static double getCombinedEnergy(ReconstructedParticle part){
+ // @TODO figure out how to combine the Ecal measurement with the
+ // SVT measurement
+ return part.getMomentum().magnitude();
+ }
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/MakeMassHistogram.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/MakeMassHistogram.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/MakeMassHistogram.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,110 @@
+package org.hps.users.spaul.bh.test;
+
+import hep.aida.IAnalysisFactory;
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogram2D;
+import hep.aida.IHistogramFactory;
+import hep.aida.ITree;
+
+import java.io.File;
+import java.io.IOException;
+import java.util.List;
+
+import org.hps.conditions.ConditionsDriver;
+import org.hps.users.spaul.feecc.CustomBinning;
+import org.hps.users.spaul.feecc.RemoveDuplicateParticles;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.Vertex;
+import org.lcsim.lcio.LCIOReader;
+
+public class MakeMassHistogram {
+
+ static String vertexCollectionName = "TargetConstrainedV0Vertices";
+ static boolean display = false;
+ static CustomBinning cb;
+ public static void main(String arg[]) throws IllegalArgumentException, IOException{
+
+
+ String input = arg[0];
+ String output = arg[1];
+
+ IAnalysisFactory af = IAnalysisFactory.create();
+ ITree tree = af.createTreeFactory().create(output,"xml",false,true);
+ IHistogramFactory hf = af.createHistogramFactory(tree);
+ setupHistograms(tree, hf);
+
+ ConditionsDriver hack = new ConditionsDriver();
+ //hack.setXmlConfigResource("/u/group/hps/hps_soft/detector-data/detectors/HPS-EngRun2015-Nominal-v3");
+ hack.setDetectorName("HPS-PhysicsRun2016-2pt2-v0");
+ hack.setFreeze(true);
+ hack.setRunNumber(Integer.parseInt(arg[2]));
+ hack.initialize();
+ File file = new File(input);
+ File inputs[] = new File[]{file};
+ if(file.isDirectory()){
+ inputs = file.listFiles();
+ }
+ for(int i = 0; i< inputs.length; i++){
+ LCIOReader reader = new LCIOReader(inputs[i]);
+ //reader.open(input);
+ //reader.
+ EventHeader event = reader.read();
+ int nEvents = 0;
+ try{
+ outer : while(event != null){
+ processEvent(event);
+
+ //System.out.println(Q2);
+
+ event = reader.read();
+ }
+ } catch (Exception e){
+ e.printStackTrace();
+ }
+ }
+ tree.commit();
+ tree.close();
+
+ }
+
+
+
+
+
+ static IHistogram1D h1;
+
+
+ static void setupHistograms(ITree tree, IHistogramFactory hf){
+ //h1 = hf.createHistogram2D("px\\/pz vs py\\/pz", 160, -.16, .24, 160, -.2, .2);
+
+ // 1 MeV. bins.
+ double maxMass = .4;
+ double massBin = .001;
+ h1 = hf.createHistogram1D("recon mass", (int) (maxMass/massBin), 0, maxMass);
+
+
+ }
+
+ private static void processEvent(EventHeader event) {
+ if(event.getEventNumber() %10000 == 0)
+ System.out.println("event number " + event.getEventNumber());
+
+
+
+
+ List<ReconstructedParticle> particles = event.get(ReconstructedParticle.class, "FinalStateParticles");
+ particles = RemoveDuplicateParticles.removeDuplicateParticles(particles);
+ List<Vertex> vertices = event.get(Vertex.class, vertexCollectionName);
+
+ //System.out.println(mA);
+ for(Vertex v : vertices){
+
+ if(Cuts.passesCuts(v))
+ h1.fill(InvariantMassCalculator.getMass(v));
+ }
+ //h9_t.fill(event.get(MCParticle.class, "MCParticles").get(0).getEnergy());
+
+ }
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/PlotMassResolutionsMC.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/PlotMassResolutionsMC.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/PlotMassResolutionsMC.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,53 @@
+package org.hps.users.spaul.bh.test;
+
+import java.io.IOException;
+
+import org.hps.users.spaul.StyleUtil;
+
+import hep.aida.*;
+
+public class PlotMassResolutionsMC {
+ public static void main(String[] arg) throws IllegalArgumentException, IOException{
+ IAnalysisFactory af = IAnalysisFactory.create();
+ ITreeFactory tf = af.createTreeFactory();
+ IPlotterFactory pf = af.createPlotterFactory();
+
+ IPlotter p1 = pf.create();
+ p1.createRegions(1,2);
+ IPlotterRegion r1 = p1.region(0);
+ IPlotterRegion r2 = p1.region(1);
+
+ IDataPointSet dps = af.createDataPointSetFactory(af.createTreeFactory().create()).create("mass resolutions", 2);
+
+ for(int i = 0; i<arg.length/2; i++){
+ double mass = Double.parseDouble(arg[2*i]);
+ String file = arg[2*i+1];
+ ITree tree = tf.create(file);
+ IHistogram1D massRecon = (IHistogram1D) tree.find("recon mass");
+
+ massRecon.setTitle(mass + " MeV");
+ massRecon.scale(1/massRecon.maxBinHeight());
+ double rms = massRecon.rms();
+ double n = massRecon.entries();
+ // http://web.eecs.umich.edu/~fessler/papers/files/tr/stderr.pdf
+ double drms = rms *Math.sqrt(Math.sqrt(2/(n-1)));
+ r1.plot(massRecon);
+
+ IDataPoint dp = dps.addPoint();
+ dp.coordinate(0).setValue(mass);
+ dp.coordinate(1).setValue(rms/mass);
+ dp.coordinate(1).setErrorPlus(drms/mass);
+ dp.coordinate(1).setErrorMinus(drms/mass);
+ }
+ r2.plot(dps);
+
+ StyleUtil.stylize(r1, "Mass Resolution", "A' mass (GeV)", "arbitrary units");
+ StyleUtil.stylize(r2, "Mass Resolution", "A' mass (GeV)", "dm/m");
+ //r2.style().dataStyle().setParameter("showLineBetweenPoints", arg1)
+ r2.setYLimits(0, .2);
+ r1.setXLimits(0, .300);
+ r2.setXLimits(0, .300);
+
+ p1.show();
+ }
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/RecoilCalculator.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/RecoilCalculator.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/bh/test/RecoilCalculator.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,78 @@
+package org.hps.users.spaul.bh.test;
+
+import java.util.Random;
+
+import hep.aida.IAnalysisFactory;
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogramFactory;
+import hep.aida.IPlotter;
+import hep.aida.ITree;
+
+public class RecoilCalculator {
+
+ //some ridiculous equations that came out of mathematica, and then translated into java.
+ //might have some mistakes in there
+ static public double getProductAB(double a, double b, double sum, double wa, double wb, double wsum){
+ return ((a*wsum*wa - b*wsum*wb + sum* wsum* wb + a*wa*wb)*(-a*wsum*wa + sum*wsum*wa +
+ b*wsum*wb + b*wa*wb))/Math.pow(wsum*(wa + wb) + wa*wb,2);
+ }
+ static public double getProductError(double a, double b, double sum, double wa, double wb, double wsum){
+ return Math.sqrt((2*b*sum*wsum*wb*(2*wsum*wsum*(wa - wb) + wsum*wa*(wa - wb) + wa*wa*wb) +
+ b*b*(wsum+ wa)*wb*wb*(4*wsum*wsum + wa*wb + wsum*(wa + wb)) +
+ a*a*wa*wa*(wsum + wb)*(4*wsum*wsum + wa*wb + wsum*(wa + wb)) +
+ sum*sum*wsum*wsum*(wsum*(wa - wb)*(wa-wb)+ wa*wb*(wa + wb)) -
+ 2*a*wsum*wa*(sum*(2*wsum*wsum*(wa - wb) + wsum*(wa - wb)*wb - wa*wb*wb) +
+ b*wb*(4*wsum*wsum + 3*wa*wb + 3*wsum*(wa + wb))))
+ /Math.pow(wa*wb + wsum*(wa + wb),3));
+
+ }
+
+ public static void main(String arg[]){
+ IAnalysisFactory af = IAnalysisFactory.create();
+ ITree tree = af.createTreeFactory().create();
+ IHistogramFactory hf = af.createHistogramFactory(tree);
+
+ IHistogram1D h1 = hf.createHistogram1D("stuff", 100, .4, 1);
+ IPlotter p = af.createPlotterFactory().create();
+ p.region(0).plot(h1);
+ p.show();
+
+ IHistogram1D h2 = hf.createHistogram1D("stuff2", 100, 0, .05);
+ IPlotter p2 = af.createPlotterFactory().create();
+ p2.region(0).plot(h2);
+ p2.show();
+
+ double smean = 2;
+ double amean = 1.5;
+ double bmean = .5;
+
+ double ds = .003;
+ double da = .005;
+ double db = .07;
+
+ System.out.println(getProductAB(amean, bmean, smean, 1/(da*da), 1/(db*db), 1/(ds*ds))
+ -getProductAB(bmean, amean, smean, 1/(db*db), 1/(da*da), 1/(ds*ds)))
+ ;
+ System.out.println(getProductError(amean, bmean, smean, 1/(da*da), 1/(db*db), 1/(ds*ds))
+ -getProductError(bmean, amean, smean, 1/(db*db), 1/(da*da), 1/(ds*ds)))
+ ;
+ Random r =new Random();
+ for(int i = 0; i< 10000; i++){
+
+ double s = r.nextGaussian()*ds + smean;
+ double a = r.nextGaussian()*da + amean;
+ double b = r.nextGaussian()*db + bmean;
+
+ double product = getProductAB(a, b, s, 1/(da*da), 1/(db*db), 1/(ds*ds));
+ double dproduct = getProductError(a,b,s, 1/(da*da), 1/(db*db), 1/(ds*ds));
+
+ h1.fill(product);
+ h2.fill(dproduct);
+ }
+
+
+
+ }
+
+
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/dimuon/DimuonFilterDriver.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/dimuon/DimuonFilterDriver.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/dimuon/DimuonFilterDriver.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,59 @@
+package org.hps.users.spaul.dimuon;
+
+import java.util.List;
+
+import org.hps.record.triggerbank.AbstractIntData;
+import org.hps.record.triggerbank.TIData;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.GenericObject;
+import org.lcsim.event.Track;
+import org.lcsim.util.Driver;
+
+public class DimuonFilterDriver extends Driver{
+ public void process(EventHeader event) {
+
+
+
+ // only keep singles triggers:
+ if (!event.hasCollection(GenericObject.class,"TriggerBank"))
+ throw new Driver.NextEventException();
+ boolean isSingles=false;
+ for (GenericObject gob : event.get(GenericObject.class,"TriggerBank"))
+ {
+ if (!(AbstractIntData.getTag(gob) == TIData.BANK_TAG)) continue;
+ TIData tid = new TIData(gob);
+ if (tid.isSingle0Trigger() || tid.isSingle1Trigger())
+ {
+ isSingles=true;
+ break;
+ }
+ }
+ if (!isSingles) throw new Driver.NextEventException();
+
+ if (!event.hasCollection(Cluster.class, "EcalClusters"))
+ throw new Driver.NextEventException();
+
+
+
+ if (!event.hasCollection(Track.class, "GBLTracks"))
+ throw new Driver.NextEventException();
+
+ List<Track> tracks = event.get(Track.class,"GBLTracks");
+
+ if(tracks.size()<2)
+ throw new Driver.NextEventException();
+
+ boolean foundPlus = false;
+ for(Track t : tracks){
+ if(t.getCharge()==1)
+ foundPlus = true;
+ }
+
+
+
+
+ }
+
+
+}
Added: java/trunk/users/src/main/java/org/hps/users/spaul/pulser/UncorrClusterDriver.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/spaul/pulser/UncorrClusterDriver.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/spaul/pulser/UncorrClusterDriver.java Thu Sep 1 04:27:30 2016
@@ -0,0 +1,92 @@
+package org.hps.users.spaul.pulser;
+
+import hep.aida.IHistogram2D;
+
+import java.util.List;
+
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.EventHeader;
+import org.lcsim.geometry.Detector;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+public class UncorrClusterDriver extends Driver{
+ AIDA aida = AIDA.defaultInstance();
+ public static boolean fid_ECal(double x, double y)
+ {
+ y = Math.abs(y);
+
+ boolean in_fid = false;
+ double x_edge_low = -262.74;
+ double x_edge_high = 347.7;
+ double y_edge_low = 33.54;
+ double y_edge_high = 75.18;
+
+ double x_gap_low = -106.66;
+ double x_gap_high = 42.17;
+ double y_gap_high = 47.18;
+
+ y = Math.abs(y);
+
+ if( x > x_edge_low && x < x_edge_high && y > y_edge_low && y < y_edge_high )
+ {
+ if( !(x > x_gap_low && x < x_gap_high && y > y_edge_low && y < y_gap_high) )
+ {
+ in_fid = true;
+ }
+ }
+
+ return in_fid;
+ }
+ private void processUncorrectedCluster(Cluster c){
+ double time = 0;
+ double seedEnergy = -100;
+ if(Math.abs(c.getPosition()[1])< 47.18)
+ return;
+ for(CalorimeterHit hit : c.getCalorimeterHits()){
+ //if(Math.abs(hit.getPositionVec().y())< 47.18)
+ //continue;
+ if(hit.getCorrectedEnergy() > seedEnergy){
+ seedEnergy = hit.getCorrectedEnergy();
+ time = hit.getTime();
+
+ }
+ }
+ double clusterEnergy = c.getEnergy();
+ //c.get
+ int size = c.getSize();
+
+
+ energyHist.fill(clusterEnergy, time);
+
+ seedHist.fill(seedEnergy, time);
+ sizeHist.fill(size, time);
+
+ if(clusterEnergy>.6 && seedEnergy>.4 && size>=3){
+ energyHistCut.fill(clusterEnergy, time);
+
+ seedHistCut.fill(seedEnergy, time);
+ sizeHistCut.fill(size, time);
+ }
+ }
+
+ IHistogram2D energyHist, seedHist, sizeHist,
+ energyHistCut, seedHistCut, sizeHistCut;
+ @Override
+ public void process(EventHeader event){
+ List<Cluster> clusters = event.get(Cluster.class, "EcalClusters");
+ for(Cluster c : clusters)
+ processUncorrectedCluster(c);
+ }
+ @Override
+ public void detectorChanged(Detector d){
+ energyHist = aida.histogram2D("uncorr/energy vs time", 100, 0, 2, 100, 0, 400);
+ seedHist = aida.histogram2D("uncorr/seed vs time", 100, 0, 2, 100, 0, 400);
+ sizeHist = aida.histogram2D("uncorr/size vs time", 10, 0, 10, 100, 0, 400);
+ energyHistCut = aida.histogram2D("uncorr/energy vs time with cuts", 100, 0, 2, 100, 0, 400);
+ seedHistCut = aida.histogram2D("uncorr/seed vs time with cuts", 100, 0, 2, 100, 0, 400);
+ sizeHistCut = aida.histogram2D("uncorr/size vs time with cuts", 10, 0, 10, 100, 0, 400);
+
+ }
+}
|
