[no commit message]

This commit is contained in:
darksnake 2016-03-22 16:57:22 +03:00
parent 793c7cefcc
commit 6d938f3380
9 changed files with 217 additions and 25 deletions

View File

@ -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);

View File

@ -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}")

View File

@ -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");

View File

@ -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();

View File

@ -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<DataPoint> config) {
ListPointSet res = new ListPointSet(adapter.getFormat());
for (Iterator<DataPoint> 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();

View File

@ -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 {

View File

@ -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<DataPoint> data, String maskForX, String maskForY, String maskForYerr, String maskForTime) {
// ListPointSet res = new ListPointSet(XYDataPoint.names);
// for (DataPoint point : data) {

View File

@ -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

View File

@ -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());