diff --git a/numass-main/src/main/groovy/inr/numass/scripts/OldTest.groovy b/numass-main/src/main/groovy/inr/numass/scripts/OldTest.groovy index d4021b0c..4979d4d3 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/OldTest.groovy +++ b/numass-main/src/main/groovy/inr/numass/scripts/OldTest.groovy @@ -50,7 +50,7 @@ FitManager fm = new FitManager(); File fssfile = new File("c:\\Users\\Darksnake\\Dropbox\\PlayGround\\FS.txt"); BivariateFunction resolution = new ResolutionFunction(2.28e-4); -resolution.setTailFunction(ResolutionFunction.getRealTail()) +//resolution.setTailFunction(ResolutionFunction.getRealTail()) ModularTritiumSpectrum sp = new ModularTritiumSpectrum(resolution, 18395d, 18580d, fssfile); sp.setCaching(false); @@ -87,11 +87,11 @@ allPars.setParDomain("trap", 0d, Double.POSITIVE_INFINITY); ListPointSet data = readData("c:\\Users\\Darksnake\\Dropbox\\PlayGround\\RUN23.DAT", 18400d); -FitState state = fm.buildState(data, model, allPars); +FitState state = new FitState(data, model, allPars); FitState res = fm.runDefaultTask(state, "E0", "N", "bkg"); -res = fm.runDefaultTask(res, "E0", "N", "bkg", "U2"); +res = fm.runDefaultTask(res, "E0", "N", "bkg", "mnu2"); res.print(out); diff --git a/numass-main/src/main/groovy/inr/numass/scripts/SystTransmission.groovy b/numass-main/src/main/groovy/inr/numass/scripts/SystTransmission.groovy new file mode 100644 index 00000000..7d7fe106 --- /dev/null +++ b/numass-main/src/main/groovy/inr/numass/scripts/SystTransmission.groovy @@ -0,0 +1,112 @@ +/* + * Copyright 2015 Alexander Nozik. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package inr.numass.scripts; + +import hep.dataforge.context.GlobalContext; +import static hep.dataforge.context.GlobalContext.out; +import hep.dataforge.points.ListPointSet; +import hep.dataforge.datafitter.FitManager; +import hep.dataforge.datafitter.FitState; +import hep.dataforge.datafitter.FitTask; +import hep.dataforge.datafitter.MINUITPlugin + +import hep.dataforge.datafitter.ParamSet; +import hep.dataforge.datafitter.models.XYModel; +import hep.dataforge.exceptions.NamingException; +import hep.dataforge.exceptions.PackFormatException; +import inr.numass.data.SpectrumDataAdapter; +import inr.numass.data.SpectrumGenerator; +import inr.numass.models.ModularTritiumSpectrum; +import inr.numass.models.NBkgSpectrum; +import inr.numass.models.ResolutionFunction +import inr.numass.utils.DataModelUtils; +import hep.dataforge.plotfit.PlotFitResultAction; +import hep.dataforge.plots.PlotFrame +import hep.dataforge.plots.data.PlottableFunction +import hep.dataforge.plots.jfreechart.JFreeChartFrame +import java.io.FileNotFoundException; +import java.util.Locale; +import org.apache.commons.math3.analysis.BivariateFunction + +import static java.util.Locale.setDefault; + +/** + * + * @author Darksnake + */ + +setDefault(Locale.US); +new MINUITPlugin().startGlobal(); + +FitManager fm = new FitManager(); + +ResolutionFunction resolution = new ResolutionFunction(8.3e-5); +resolution.setTailFunction(ResolutionFunction.getRealTail()); +ModularTritiumSpectrum beta = new ModularTritiumSpectrum(resolution, 18395d, 18580d, null); +beta.setCaching(false); + +NBkgSpectrum spectrum = new NBkgSpectrum(beta); +XYModel model = new XYModel("tritium", spectrum, new SpectrumDataAdapter()); + +ParamSet allPars = new ParamSet(); + + +allPars.setPar("N", 6e7, 1e5, 0, Double.POSITIVE_INFINITY); + +allPars.setPar("bkg", 2, 0.1 ); + +allPars.setPar("E0", 18575.0, 0.1 ); + +allPars.setPar("mnu2", 0, 2); + +def mster = 3000;// Mass of sterile neutrino in eV + +allPars.setPar("msterile2", mster**2, 1); + +allPars.setPar("U2", 0, 1e-4); + +allPars.setPar("X", 0, 0.05, 0d, Double.POSITIVE_INFINITY); + +allPars.setPar("trap", 1, 0.01, 0d, Double.POSITIVE_INFINITY); + +int seed = 12316 +SpectrumGenerator generator = new SpectrumGenerator(model, allPars, seed); + +def config = DataModelUtils.getUniformSpectrumConfiguration(18400d, 18580, 1e7, 60) +//def config = DataModelUtils.getSpectrumConfigurationFromResource("/data/run23.cfg") + +ListPointSet data = generator.generateExactData(config); + +FitState state = new FitState(data, model, allPars); + +println("Simulating data with real tail. Seed = ${seed}") + +println("Fitting data with real parameters") + +FitState res = fm.runTask(state, "QOW", FitTask.TASK_RUN, "N", "bkg","E0", "mnu2"); +res.print(out()); + +def mnu2 = res.getParameters().getValue("mnu2"); + +println("Setting constant tail and fitting") +resolution.setTailFunction(ResolutionFunction.getConstantTail()); + +res = fm.runTask(state, "QOW", FitTask.TASK_RUN, "N", "bkg","E0", "mnu2"); +res.print(out()); + +def diff = res.getParameters().getValue("mnu2") - mnu2; + +println("\n\nSquared mass difference: ${diff}") \ No newline at end of file diff --git a/numass-main/src/main/groovy/inr/numass/scripts/Systematics.groovy b/numass-main/src/main/groovy/inr/numass/scripts/Systematics.groovy index 3e9d794b..a9c33afe 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/Systematics.groovy +++ b/numass-main/src/main/groovy/inr/numass/scripts/Systematics.groovy @@ -88,7 +88,7 @@ ListPointSet data = generator.generateData(DataModelUtils.getUniformSpectrumConf // data = data.filter("X", Value.of(15510.0), Value.of(18610.0)); allPars.setParValue("U2", 0); -FitState state = FitManager.buildState(data, model, allPars); +FitState state = new FitState(data, model, allPars); //new PlotFitResultAction(GlobalContext.instance(), null).runOne(state); //double delta = 4e-6; @@ -99,8 +99,8 @@ FitState state = FitManager.buildState(data, model, allPars); resolution.setTailFunction(ResolutionFunction.getRealTail()) -PlotFrame frame = JFreeChartFrame.drawFrame("Transmission function", null); -frame.add(new PlottableFunction("transmission",null, {U -> resolution.value(18500,U)},13500,18505,500)); +//PlotFrame frame = JFreeChartFrame.drawFrame("Transmission function", null); +//frame.add(new PlottableFunction("transmission",null, {U -> resolution.value(18500,U)},13500,18505,500)); FitState res = fm.runTask(state, "QOW", FitTask.TASK_RUN, "N", "bkg", "E0", "U2", "trap"); diff --git a/numass-main/src/main/java/inr/numass/data/SpectrumDataAdapter.java b/numass-main/src/main/java/inr/numass/data/SpectrumDataAdapter.java index 901e78e2..d848f4cd 100644 --- a/numass-main/src/main/java/inr/numass/data/SpectrumDataAdapter.java +++ b/numass-main/src/main/java/inr/numass/data/SpectrumDataAdapter.java @@ -15,15 +15,15 @@ */ package inr.numass.data; -import hep.dataforge.points.DataPoint; -import hep.dataforge.points.MapPoint; -import hep.dataforge.points.XYAdapter; import hep.dataforge.exceptions.DataFormatException; import hep.dataforge.exceptions.NameNotFoundException; import hep.dataforge.meta.Meta; import hep.dataforge.meta.MetaBuilder; -import hep.dataforge.values.Value; +import hep.dataforge.points.DataPoint; +import hep.dataforge.points.MapPoint; import hep.dataforge.points.PointAdapter; +import hep.dataforge.points.XYAdapter; +import hep.dataforge.values.Value; /** * @@ -82,7 +82,7 @@ public class SpectrumDataAdapter extends XYAdapter { @Override public Value getYerr(DataPoint point) throws NameNotFoundException { - if (providesYError(point)) { + if (super.providesYError(point)) { return Value.of(super.getYerr(point).doubleValue() / getTime(point)); } else { double y = super.getY(point).doubleValue(); diff --git a/numass-main/src/main/java/inr/numass/data/SpectrumGenerator.java b/numass-main/src/main/java/inr/numass/data/SpectrumGenerator.java index 06b08cab..061a4649 100644 --- a/numass-main/src/main/java/inr/numass/data/SpectrumGenerator.java +++ b/numass-main/src/main/java/inr/numass/data/SpectrumGenerator.java @@ -15,21 +15,18 @@ */ package inr.numass.data; -import hep.dataforge.points.DataPoint; -import hep.dataforge.points.ListPointSet; import hep.dataforge.datafitter.ParamSet; import hep.dataforge.datafitter.models.Generator; import hep.dataforge.datafitter.models.XYModel; import static hep.dataforge.maths.RandomUtils.getDefaultRandomGenerator; +import hep.dataforge.points.DataPoint; +import hep.dataforge.points.ListPointSet; import static java.lang.Double.isNaN; import static java.lang.Math.sqrt; import java.util.Iterator; import org.apache.commons.math3.random.JDKRandomGenerator; import org.apache.commons.math3.random.RandomDataGenerator; import org.apache.commons.math3.random.RandomGenerator; -import static java.lang.Double.isNaN; -import static java.lang.Double.isNaN; -import static java.lang.Double.isNaN; /** * Генератор наборов данных для спектров. На входе требуется набор данных, @@ -74,6 +71,26 @@ public class SpectrumGenerator implements Generator { return res; } + /** + * Generate spectrum points with error derived from configuration but with + * zero spread (exactly the same as model provides) + * + * @param config + * @return + */ + public ListPointSet generateExactData(Iterable config) { + ListPointSet res = new ListPointSet(adapter.getFormat()); + for (Iterator it = config.iterator(); it.hasNext();) { + res.add(this.generateExactDataPoint(it.next())); + } + return res; + } + + public DataPoint generateExactDataPoint(DataPoint configPoint) { + double mu = this.getMu(configPoint); + return adapter.buildSpectrumDataPoint(this.getX(configPoint), (long) mu, this.getTime(configPoint)); + } + @Override public DataPoint generateDataPoint(DataPoint configPoint) { double mu = this.getMu(configPoint); @@ -110,7 +127,7 @@ public class SpectrumGenerator implements Generator { double time = this.getTime(configPoint); - return adapter.buildSpectrumDataPoint(this.getX(configPoint), (long)y, time); + return adapter.buildSpectrumDataPoint(this.getX(configPoint), (long) y, time); } @Override @@ -128,9 +145,8 @@ public class SpectrumGenerator implements Generator { // } // return sqrt(this.getMu(point)); // } - private double getTime(DataPoint point) { - + return adapter.getTime(point); // if (point.containsName("time")) { // return point.getValue("time").doubleValue(); diff --git a/numass-main/src/main/java/inr/numass/models/ResolutionFunction.java b/numass-main/src/main/java/inr/numass/models/ResolutionFunction.java index 20c26733..9a796b5b 100644 --- a/numass-main/src/main/java/inr/numass/models/ResolutionFunction.java +++ b/numass-main/src/main/java/inr/numass/models/ResolutionFunction.java @@ -103,6 +103,10 @@ public class ResolutionFunction implements BivariateFunction { return (double x, double y) -> f.value(x - y); } + + public static BivariateFunction getConstantTail(){ + return new ConstantTailFunction(); + } private static class ConstantTailFunction implements BivariateFunction { diff --git a/numass-main/src/main/java/inr/numass/utils/DataModelUtils.java b/numass-main/src/main/java/inr/numass/utils/DataModelUtils.java index 2c312d05..17bf364f 100644 --- a/numass-main/src/main/java/inr/numass/utils/DataModelUtils.java +++ b/numass-main/src/main/java/inr/numass/utils/DataModelUtils.java @@ -18,7 +18,7 @@ package inr.numass.utils; import hep.dataforge.points.DataPoint; import hep.dataforge.points.ListPointSet; import hep.dataforge.points.MapPoint; - +import java.util.Scanner; /** * @@ -34,13 +34,25 @@ public class DataModelUtils { for (int i = 0; i < numpoints; i++) { // формула работает даже в том случае когда порядок точек обратный double x = from + (to - from) / (numpoints - 1) * i; - DataPoint point = new MapPoint(list, x,time); + DataPoint point = new MapPoint(list, x, time); res.add(point); } return res; } - + + public static ListPointSet getSpectrumConfigurationFromResource(String resource) { + final String[] list = {"x", "time"}; + ListPointSet res = new ListPointSet(list); + Scanner scan = new Scanner(DataModelUtils.class.getResourceAsStream(resource)); + while (scan.hasNextLine()) { + double x = scan.nextDouble(); + int time = scan.nextInt(); + res.add(new MapPoint(list, x, time)); + } + return res; + } + // public static ListPointSet maskDataSet(Iterable data, String maskForX, String maskForY, String maskForYerr, String maskForTime) { // ListPointSet res = new ListPointSet(XYDataPoint.names); // for (DataPoint point : data) { diff --git a/numass-main/src/main/resources/data/run23.cfg b/numass-main/src/main/resources/data/run23.cfg new file mode 100644 index 00000000..d79ff862 --- /dev/null +++ b/numass-main/src/main/resources/data/run23.cfg @@ -0,0 +1,48 @@ +18400.0 7300 +18425.0 7300 +18450.0 7300 +18475.0 7300 +18500.0 14600 +18505.0 14600 +18510.0 14600 +18515.0 14600 +18520.0 14600 +18525.0 14600 +18530.0 14600 +18535.0 14600 +18540.0 14600 +18542.0 14600 +18544.0 14600 +18546.0 14600 +18548.0 14600 +18550.0 14600 +18552.0 14600 +18553.0 14600 +18554.0 14600 +18555.0 20500 +18556.0 20500 +18557.0 20500 +18558.0 20500 +18559.0 20500 +18560.0 29200 +18561.0 29200 +18562.0 29200 +18563.0 29200 +18564.0 29200 +18565.0 29200 +18566.0 29200 +18567.0 29200 +18568.0 29200 +18569.0 29200 +18570.0 29200 +18571.0 29200 +18572.0 23300 +18573.0 23300 +18574.0 23300 +18575.0 14600 +18577.0 14600 +18580.0 14600 +18585.0 14600 +18590.0 14600 +18670.0 14600 +18770.0 14600 \ No newline at end of file diff --git a/numass-prop/src/main/java/inr/numass/prop/PropTest.java b/numass-prop/src/main/java/inr/numass/prop/PropTest.java index 93af1e50..4d81abc8 100644 --- a/numass-prop/src/main/java/inr/numass/prop/PropTest.java +++ b/numass-prop/src/main/java/inr/numass/prop/PropTest.java @@ -17,7 +17,6 @@ package inr.numass.prop; import hep.dataforge.context.GlobalContext; import static hep.dataforge.context.GlobalContext.out; -import hep.dataforge.points.DataPoint; import hep.dataforge.datafitter.FitManager; import hep.dataforge.datafitter.FitState; import hep.dataforge.datafitter.ParamSet; @@ -26,10 +25,11 @@ import hep.dataforge.datafitter.models.HistogramModel; import hep.dataforge.functions.ParametricFunction; import hep.dataforge.maths.MatrixOperations; import hep.dataforge.maths.RandomUtils; +import hep.dataforge.points.DataPoint; +import hep.dataforge.points.PointSet; import inr.numass.models.BetaSpectrum; import inr.numass.models.NBkgSpectrum; import java.io.FileNotFoundException; -import hep.dataforge.points.PointSet; /** * Hello world! @@ -84,7 +84,7 @@ public class PropTest { // allPars.setParValue("base", 1e-3); // allPars.setParValue("w", 5.470); allPars.setParValue("dw", 2e-2); - FitState state = FitManager.buildState(data, model, allPars); + FitState state = new FitState(data, model, allPars); FitState res = fm.runDefaultTask(state, "U2", "E0", "N"); res.print(out());