Merge remote-tracking branch 'spc/master'

This commit is contained in:
Alexander Nozik 2023-11-08 11:46:53 +03:00
commit 9a1e9dc996
Signed by: altavir
GPG Key ID: B10A55DCC7B9AEBB
35 changed files with 548 additions and 2338 deletions

View File

@ -1,6 +1,9 @@
import org.jetbrains.kotlin.gradle.dsl.KotlinJvmProjectExtension
import org.jetbrains.kotlin.gradle.tasks.KotlinCompile
plugins {
kotlin("jvm") version "1.5.31"
id("org.openjfx.javafxplugin") version "0.0.10" apply false
kotlin("jvm") version "1.9.20" apply false
id("org.openjfx.javafxplugin") version "0.1.0" apply false
id("com.github.johnrengelman.shadow") version "7.1.0" apply false
}
@ -15,47 +18,44 @@ allprojects {
maven("https://oss.sonatype.org/content/repositories/snapshots")
}
extensions.findByType<KotlinJvmProjectExtension>()?.jvmToolchain(16)
dependencies {
api(kotlin("reflect"))
api("org.jetbrains:annotations:23.0.0")
testImplementation("junit:junit:4.13.2")
add("api", kotlin("reflect"))
add("api", "org.jetbrains:annotations:24.0.1")
add("testImplementation", "junit:junit:4.13.2")
//Spock dependencies. To be removed later
// https://mvnrepository.com/artifact/org.spockframework/spock-core
testImplementation("org.spockframework:spock-core:2.0-groovy-3.0")
add("testImplementation", "org.spockframework:spock-core:2.3-groovy-4.0")
}
tasks {
compileJava{
withType<Jar> {
duplicatesStrategy = DuplicatesStrategy.EXCLUDE
}
withType<JavaCompile> {
options.encoding = "UTF-8"
}
compileTestJava{
options.encoding = "UTF-8"
}
compileKotlin {
withType<KotlinCompile> {
kotlinOptions {
jvmTarget = "16"
javaParameters = true
freeCompilerArgs = freeCompilerArgs + listOf(
"-Xjvm-default=all",
"-progressive",
"-Xuse-experimental=kotlin.Experimental"
"-Xjvm-default=all"
)
}
}
compileTestKotlin {
kotlinOptions {
jvmTarget = "16"
javaParameters = true
freeCompilerArgs = freeCompilerArgs + listOf(
"-Xjvm-default=all",
"-progressive",
"-Xuse-experimental=kotlin.Experimental"
)
withType<GroovyCompile>{
val compileKotlinTask = getByName<KotlinCompile>("compileKotlin")
dependsOn(compileKotlinTask)
afterEvaluate {
classpath += files(compileKotlinTask.destinationDirectory)
}
}
}

View File

@ -12,9 +12,7 @@ description = "A tornadofx based kotlin library"
dependencies {
api(project(":dataforge-plots"))
api(project(":dataforge-gui:dataforge-html"))
api("no.tornado:tornadofx:1.7.20"){
exclude("org.jetbrains.kotlin")
}
api("no.tornado:tornadofx:1.7.20")
api("org.controlsfx:controlsfx:11.1.0")
api("no.tornado:tornadofx-controlsfx:0.1.1")
api("org.fxmisc.richtext:richtextfx:0.10.7")
@ -24,4 +22,3 @@ dependencies {
//compileOnly project(":dataforge-plots:plots-jfc")
}

View File

@ -10,9 +10,9 @@ javafx {
description = "jFreeChart plugin"
dependencies {
api("org.jfree:jfreesvg:3.4.2")
api("org.jfree:org.jfree.svg:5.0.5")
// https://mvnrepository.com/artifact/org.jfree/org.jfree.chart.fx
api(group= "org.jfree", name= "jfreechart-fx", version= "1.0.1")
api("org.jfree:org.jfree.chart.fx:2.0.1")
api(project(":dataforge-plots"))

View File

@ -19,7 +19,6 @@ application{
compileKotlin {
kotlinOptions {
jvmTarget = "1.8"
javaParameters = true
}
}

View File

@ -1,4 +1,3 @@
description = 'dataforge-minuit'
dependencies {
api project(':dataforge-stat')

View File

@ -1,5 +1,5 @@
distributionBase=GRADLE_USER_HOME
distributionPath=wrapper/dists
distributionUrl=https\://services.gradle.org/distributions/gradle-7.3-bin.zip
distributionUrl=https\://services.gradle.org/distributions/gradle-8.4-bin.zip
zipStoreBase=GRADLE_USER_HOME
zipStorePath=wrapper/dists

View File

@ -2,13 +2,10 @@ apply plugin: 'groovy'
description = 'The GRIND (GRoovy INteractive Dataforge) environment'
compileGroovy.dependsOn(compileKotlin)
compileGroovy.classpath += files(compileKotlin.destinationDir)
dependencies {
api project(":dataforge-core")
// https://mvnrepository.com/artifact/org.codehaus.groovy/groovy-all
api 'org.codehaus.groovy:groovy-all:3.0.9'
api 'org.apache.groovy:groovy-all:4.0.15'
testImplementation project(":dataforge-gui")

View File

@ -1,11 +1,9 @@
import com.google.protobuf.gradle.protobuf
import com.google.protobuf.gradle.protoc
import org.jetbrains.kotlin.gradle.tasks.KotlinCompile
plugins {
idea
kotlin("jvm")
id("com.google.protobuf") version "0.8.16"
id("com.google.protobuf") version "0.9.4"
}
@ -14,13 +12,12 @@ repositories {
}
dependencies {
api("com.google.protobuf:protobuf-java:3.17.1")
api("com.google.protobuf:protobuf-java:3.25.0")
api(project(":numass-core:numass-data-api"))
api(project(":dataforge-storage"))
}
tasks.withType<KotlinCompile> {
kotlinOptions.jvmTarget = "1.8"
dependsOn(":numass-core:numass-data-proto:generateProto")
}
@ -36,7 +33,7 @@ protobuf {
// Configure the protoc executable
protoc {
// Download from repositories
artifact = "com.google.protobuf:protoc:3.17.1"
artifact = "com.google.protobuf:protoc:3.25.0"
}
generatedFilesBaseDir = "$projectDir/gen"
}

View File

@ -1,5 +1,4 @@
plugins {
id 'groovy'
id 'application'
id "org.openjfx.javafxplugin"
}
@ -15,18 +14,16 @@ mainClassName = 'inr.numass.LaunchGrindShell'
description = "Main numass project"
compileGroovy.dependsOn(compileKotlin)
compileGroovy.classpath += files(compileKotlin.destinationDir)
dependencies {
api group: 'commons-cli', name: 'commons-cli', version: '1.+'
api group: 'commons-io', name: 'commons-io', version: '2.+'
api project(':numass-core')
api project(':numass-core:numass-signal-processing')
compileOnly "org.jetbrains.kotlin:kotlin-main-kts:1.3.21"
compileOnly "org.jetbrains.kotlin:kotlin-main-kts:1.9.20"
api project(':dataforge-stat:dataforge-minuit')
api project(':grind:grind-terminal')
api project(":dataforge-gui")
api project(':dataforge-plots:plots-jfc')
//api "hep.dataforge:dataforge-html"
// https://mvnrepository.com/artifact/org.ehcache/ehcache

View File

@ -1,23 +0,0 @@
package inr.numass.scripts
import hep.dataforge.tables.Table
import hep.dataforge.workspace.FileBasedWorkspace
import hep.dataforge.workspace.Workspace
import java.nio.file.Paths
/**
* Created by darksnake on 26-Dec-16.
*/
Workspace numass = FileBasedWorkspace.build(Paths.get("D:/Work/Numass/sterile2016_10/workspace.groovy"))
numass.runTask("prepare", "fill_1_all").forEachData(Table) {
Table table = it.get();
def dp18 = table.find { it["Uset"] == 18000 }
def dp17 = table.find { it["Uset"] == 17000 }
println "${it.name}\t${dp18["CR"]}\t${dp18["CRerr"]}\t${dp17["CR"]}\t${dp17["CRerr"]}"
}

View File

@ -1,109 +0,0 @@
/*
* 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.maths.integration.UnivariateIntegrator
import hep.dataforge.plots.PlotFrame
import hep.dataforge.plots.data.XYFunctionPlot
import hep.dataforge.plots.jfreechart.JFreeChartFrame
import hep.dataforge.stat.fit.ParamSet
import inr.numass.models.misc.LossCalculator
import org.apache.commons.math3.analysis.UnivariateFunction
import org.apache.commons.math3.analysis.solvers.BisectionSolver
ParamSet params = new ParamSet()
.setParValue("exPos", 12.76)
.setParValue("ionPos", 13.95)
.setParValue("exW", 1.2)
.setParValue("ionW", 13.5)
.setParValue("exIonRatio", 4.55)
UnivariateFunction scatterFunction = LossCalculator.getSingleScatterFunction(params);
PlotFrame frame = JFreeChartFrame.drawFrame("Differential scatter function", null);
frame.add(XYFunctionPlot.plotFunction("differential", scatterFunction, 0, 100, 400));
UnivariateIntegrator integrator = NumassContext.defaultIntegrator;
double border = 13.6;
UnivariateFunction ratioFunction = { e -> integrator.integrate(0, e, scatterFunction) / integrator.integrate(e, 100, scatterFunction) }
double ratio = ratioFunction.value(border);
println "The true excitation to ionization ratio with border energy $border is $ratio";
double resolution = 1.5d;
def X = 0.527;
LossCalculator calculator = LossCalculator.INSTANCE;
List<Double> lossProbs = calculator.getGunLossProbabilities(X);
UnivariateFunction newScatterFunction = { double d ->
double res = scatterFunction.value(d);
for (i = 1; i < lossProbs.size(); i++) {
res += lossProbs.get(i) * calculator.getLossValue(i, d, 0);
}
return res;
}
UnivariateFunction resolutionValue = { double e ->
if (e <= 0d) {
return 0d;
} else if (e >= resolution) {
return 1d;
} else {
return e / resolution;
}
};
UnivariateFunction integral = { double u ->
if (u <= 0d) {
return 0d;
} else {
UnivariateFunction integrand = { double e -> resolutionValue.value(u - e) * newScatterFunction.value(e) };
return integrator.integrate(0d, u, integrand)
}
}
frame.add(XYFunctionPlot.plotFunction("integral", integral, 0, 100, 800));
BisectionSolver solver = new BisectionSolver(1e-3);
UnivariateFunction integralShifted = { u ->
def integr = integral.value(u);
return integr / (1 - integr) - ratio;
}
double integralBorder = solver.solve(400, integralShifted, 10d, 20d);
println "The integral border is $integralBorder";
double newBorder = 14.43
double integralValue = integral.value(newBorder);
double err = Math.abs(integralValue / (1 - integralValue) / ratio - 1d)
println "The relative error ic case of using $newBorder instead of real one is $err";

View File

@ -1,46 +0,0 @@
/*
* To change this license header, choose License Headers in Project Properties.
* To change this template file, choose Tools | Templates
* and open the template in the editor.
*/
package inr.numass.scripts
import hep.dataforge.maths.integration.GaussRuleIntegrator
import hep.dataforge.maths.integration.UnivariateIntegrator
import org.apache.commons.math3.analysis.UnivariateFunction
UnivariateIntegrator integrator = new GaussRuleIntegrator(400);
def exPos = 12.695;
def ionPos = 13.29;
def exW = 1.22;
def ionW = 11.99;
def exIonRatio = 3.6;
def cutoff = 25d
UnivariateFunction func = {double eps ->
if (eps <= 0) {
return 0;
}
double z1 = eps - exPos;
double ex = exIonRatio * Math.exp(-2 * z1 * z1 / exW / exW);
double z = 4 * (eps - ionPos) * (eps - ionPos);
double ion = 1 / (1 + z / ionW / ionW);
double res;
if (eps < exPos) {
res = ex;
} else {
res = Math.max(ex, ion);
}
return res;
};
//caclulating lorentz integral analythically
double tailNorm = (Math.atan((ionPos - cutoff) * 2d / ionW) + 0.5 * Math.PI) * ionW / 2d;
final double norm = integrator.integrate(0d, cutoff, func) + tailNorm;
println 1/norm;

View File

@ -1,41 +0,0 @@
/*
* To change this license header, choose License Headers in Project Properties.
* To change this template file, choose Tools | Templates
* and open the template in the editor.
*/
package inr.numass.scripts
import hep.dataforge.maths.integration.GaussRuleIntegrator
import hep.dataforge.maths.integration.UnivariateIntegrator
import org.apache.commons.math3.analysis.UnivariateFunction
UnivariateIntegrator integrator = new GaussRuleIntegrator(400);
def exPos = 12.587;
def ionPos = 11.11;
def exW = 1.20;
def ionW = 11.02;
def exIonRatio = 2.43;
def cutoff = 20d
UnivariateFunction loss = LossCalculator.getSingleScatterFunction(exPos, ionPos, exW, ionW, exIonRatio);
println integrator.integrate(0, 600, loss);
println integrator.integrate(0, cutoff, loss);
println integrator.integrate(cutoff, 600d, loss);
println (integrator.integrate(0, cutoff, loss) + integrator.integrate(cutoff, 3000d, loss));
//double tailValue = (Math.atan((ionPos-cutoff)*2d/ionW) + 0.5*Math.PI)*ionW/2;
//println tailValue
//println integrator.integrate(loss,0,100);
//println integrator.integrate(loss,100,600);
//def lorentz = {eps->
// double z = 4 * (eps - ionPos) * (eps - ionPos);
// 1 / (1 + z / ionW / ionW);
//}
//
//println(integrator.integrate(lorentz, cutoff, 800))

View File

@ -1,104 +0,0 @@
/*
* 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.Global
import hep.dataforge.stat.fit.FitManager
import hep.dataforge.stat.fit.FitState
import hep.dataforge.stat.fit.MINUITPlugin
import hep.dataforge.stat.fit.ParamSet
import hep.dataforge.stat.models.XYModel
import hep.dataforge.tables.ListTable
import inr.numass.data.SpectrumAdapter
import inr.numass.models.BetaSpectrum
import inr.numass.models.ModularSpectrum
import inr.numass.models.NBkgSpectrum
import inr.numass.models.ResolutionFunction
import org.apache.commons.math3.analysis.BivariateFunction
import static inr.numass.utils.OldDataReader.readData
PrintWriter out = Global.out();
Locale.setDefault(Locale.US);
new MINUITPlugin().startGlobal();
FitManager fm = new FitManager();
// setSeed(543982);
File fssfile = new File("c:\\Users\\Darksnake\\Dropbox\\PlayGround\\FS.txt");
BivariateFunction resolution = new ResolutionFunction(2.28e-4);
//resolution.setTailFunction(ResolutionFunction.getRealTail())
ModularSpectrum sp = new ModularSpectrum(new BetaSpectrum(fssfile), resolution, 18395d, 18580d);
sp.setCaching(false);
//RangedNamedSetSpectrum beta = new BetaSpectrum(fssfile);
//ModularSpectrum sp = new ModularSpectrum(beta, 2.28e-4, 18395d, 18580d);
// ModularTritiumSpectrum beta = new ModularTritiumSpectrum(2.28e-4, 18395d, 18580d, "d:\\PlayGround\\FS.txt");
NBkgSpectrum spectrum = new NBkgSpectrum(sp);
XYModel model = new XYModel("tritium", spectrum, new SpectrumAdapter());
ParamSet allPars = new ParamSet();
allPars.setParValue("N", 602533.94);
//значение 6е-6 соответствует полной интенстивности 6е7 распадов в секунду
//Проблема была в переполнении счетчика событий в генераторе. Заменил на long. Возможно стоит поставить туда число с плавающей точкой
allPars.setParError("N", 1000);
allPars.setParDomain("N", 0d, Double.POSITIVE_INFINITY);
allPars.setParValue("bkg", 0.012497889);
allPars.setParError("bkg", 1e-4);
allPars.setParValue("E0", 18575.986);
allPars.setParError("E0", 0.05);
allPars.setParValue("mnu2", 0d);
allPars.setParError("mnu2", 1d);
allPars.setParValue("msterile2", 50 * 50);
allPars.setParValue("U2", 0);
allPars.setParError("U2", 1e-2);
allPars.setParDomain("U2", -1d, 1d);
allPars.setParValue("X", 0.47);
allPars.setParError("X", 0.014);
allPars.setParDomain("X", 0d, Double.POSITIVE_INFINITY);
allPars.setParValue("trap", 1d);
allPars.setParError("trap", 0.2d);
allPars.setParDomain("trap", 0d, Double.POSITIVE_INFINITY);
ListTable data = readData("c:\\Users\\Darksnake\\Dropbox\\PlayGround\\RUN23.DAT", 18400d);
FitState state = new FitState(data, model, allPars);
FitState res = fm.runDefaultStage(state, "E0", "N", "bkg");
res = fm.runDefaultStage(res, "E0", "N", "bkg", "mnu2");
res.print(out);
//spectrum.counter.print(onComplete);
//
//// fm.setPriorProb(new Gaussian("X", 0.47, 0.47*0.03));
//// fm.setPriorProb(new MultivariateGaussianPrior(allPars.getSubSet("X","trap")));
//res = fm.runStage(res, "MINUIT", "run", "E0", "N", "bkg", "mnu2");
////
//res.print(onComplete);
//sp.setCaching(true);
//sp.setSuppressWarnings(true);
//
//BayesianConfidenceLimit bm = new BayesianConfidenceLimit();
//bm.printMarginalLikelihood(onComplete, "U2", res, ["E0", "N", "bkg", "U2", "X"], 10000);
// PrintNamed.printLike2D(Out.onComplete, "like", res, "N", "E0", 30, 60, 2);

View File

@ -1,61 +0,0 @@
/*
* To change this license header, choose License Headers in Project Properties.
* To change this template file, choose Tools | Templates
* and open the template in the editor.
*/
package inr.numass.scripts
import inr.numass.models.misc.LossCalculator
LossCalculator loss = LossCalculator.INSTANCE
def X = 0.36
def lossProbs = loss.getGunLossProbabilities(X);
printf("%8s\t%8s\t%8s\t%8s\t%n",
"eps",
"p1",
"p2",
"p3"
)
/*
'exPos' = 12.587 ± 0.049
'ionPos' = 11.11 ± 0.50
'exW' = 1.20 ± 0.12
'ionW' = 11.02 ± 0.68
'exIonRatio' = 2.43 ± 0.42
*/
def singleScatter = loss.getSingleScatterFunction(
12.860,
16.62,
1.71,
12.09,
4.59
);
for (double d = 0; d < 30; d += 0.3) {
double ei = 18500;
double ef = ei - d;
printf("%8f\t%8f\t%8f\t%8f\t%n",
d,
lossProbs[1] * singleScatter.value(ei - ef),
lossProbs[2] * loss.getLossValue(2, ei, ef),
lossProbs[3] * loss.getLossValue(3, ei, ef)
)
}
for (double d = 30; d < 100; d += 1) {
double ei = 18500;
double ef = ei - d;
printf("%8f\t%8f\t%8f\t%8f\t%n",
d,
lossProbs[1] * singleScatter.value(ei - ef),
lossProbs[2] * loss.getLossValue(2, ei, ef),
lossProbs[3] * loss.getLossValue(3, ei, ef)
)
}

View File

@ -1,15 +0,0 @@
/*
* To change this license header, choose License Headers in Project Properties.
* To change this template file, choose Tools | Templates
* and open the template in the editor.
*/
package inr.numass.scripts
import hep.dataforge.io.ColumnedDataReader
import hep.dataforge.io.ColumnedDataWriter
import hep.dataforge.tables.Table
File file = new File("D:\\Work\\Numass\\sterile2016\\empty.dat" )
Table referenceTable = new ColumnedDataReader(file).toTable();
ColumnedDataWriter.writeTable(System.out, referenceTable,"")

View File

@ -1,65 +0,0 @@
package inr.numass.scripts
import hep.dataforge.context.Global
import hep.dataforge.grind.Grind
import hep.dataforge.meta.Meta
import hep.dataforge.stat.fit.FitManager
import hep.dataforge.stat.fit.ParamSet
import hep.dataforge.stat.models.XYModel
import inr.numass.NumassPlugin
/**
* Created by darksnake on 14-Dec-16.
*/
Locale.setDefault(Locale.US);
new NumassPlugin().startGlobal();
def fm = Global.instance().provide("fitting", FitManager.class).getFitManager();
def mm = fm.modelManager
Meta modelMeta = Grind.buildMeta(modelName: "sterile") {
resolution(width: 8.3e-5, tailAlpha: 3e-3)
transmission(trapping: "function::numass.trap.nominal") // numass.trap.nominal = 1.2e-4 - 4.5e-9 * Ei
}
/*
'N' = 2.76515e+06 ± 2.4e+03 (0.00000,Infinity)
'bkg' = 41.195 ± 0.053
'E0' = 18576.35 ± 0.32
'mnu2' = 0.00 ± 0.010
'msterile2' = 1000000.00 ± 1.0
'U2' = 0.00314 ± 0.0010
'X' = 0.12000 ± 0.010 (0.00000,Infinity)
'trap' = 1.089 ± 0.026
*/
Meta paramMeta = Grind.buildMeta("params") {
N(value: 2.76515e+06, err: 30, lower: 0)
bkg(value: 41.195, err: 0.1)
E0(value: 18576.35, err: 0.1)
mnu2(value: 0, err: 0.01)
msterile2(value: 1000**2, err: 1)
U2(value: 0.00314, err: 1e-3)
X(value: 0.12000, err: 0.01, lower: 0)
trap(value: 1.089, err: 0.05)
}
XYModel model = mm.getModel(modelMeta)
ParamSet allPars = ParamSet.fromMeta(paramMeta);
def a = 16000;
def b = 18600;
def step = 50;
for (double x = a; x < b; x += step) {
println "${x}\t${model.value(x, allPars)}"
}
Global.terminate()

View File

@ -1,110 +0,0 @@
/*
* 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.Global
import hep.dataforge.meta.MetaBuilder
import hep.dataforge.stat.fit.ParamSet
import inr.numass.data.SpectrumInformation
import inr.numass.models.BetaSpectrum
import inr.numass.models.ModularSpectrum
import inr.numass.models.NBkgSpectrum
import inr.numass.models.ResolutionFunction
import org.apache.commons.math3.analysis.UnivariateFunction
import static java.util.Locale.setDefault
setDefault(Locale.US);
Global global = Global.instance();
// global.loadModule(new MINUIT());
// FitManager fm = new FitManager("data 2013");
UnivariateFunction reolutionTail = {x ->
if (x > 1500) {
return 0.98;
} else //Intercept = 1.00051, Slope = -1.3552E-5
{
return 1.00051 - 1.3552E-5 * x;
}
};
ModularSpectrum beta = new ModularSpectrum(new BetaSpectrum(),
new ResolutionFunction(8.3e-5, reolutionTail), 14490d, 19001d);
beta.setCaching(false);
NBkgSpectrum spectrum = new NBkgSpectrum(beta);
// XYModel model = new XYModel("tritium", spectrum);
ParamSet allPars = new ParamSet();
allPars.setParValue("N", 3090.1458);
//значение 6е-6 соответствует полной интенстивности 6е7 распадов в секунду
//Проблема была в переполнении счетчика событий в генераторе. Заменил на long. Возможно стоит поставить туда число с плавающей точкой
allPars.setParError("N", 6);
allPars.setParDomain("N", 0d, Double.POSITIVE_INFINITY);
allPars.setParValue("bkg", 2.2110028);
allPars.setParError("bkg", 0.03);
allPars.setParValue("E0", 18580.742);
allPars.setParError("E0", 2);
allPars.setParValue("mnu2", 0d);
allPars.setParError("mnu2", 1d);
allPars.setParValue("msterile2", 1000 * 1000);
allPars.setParValue("U2", 0);
allPars.setParError("U2", 1e-4);
allPars.setParDomain("U2", -1d, 1d);
allPars.setParValue("X", 1.0);
allPars.setParError("X", 0.01);
allPars.setParDomain("X", 0d, Double.POSITIVE_INFINITY);
allPars.setParValue("trap", 1.0d);
allPars.setParError("trap", 0.01d);
allPars.setParDomain("trap", 0d, Double.POSITIVE_INFINITY);
SpectrumInformation sign = new SpectrumInformation(spectrum);
// double Elow = 14000d;
// double Eup = 18600d;
// int numpoints = (int) ((Eup - Elow) / 50);
// double time = 1e6 / numpoints;
// DataSet config = getUniformSpectrumConfiguration(Elow, Eup, time, numpoints);
// NamedMatrix infoMatrix = sign.getInformationMatrix(allPars, config,"U2","E0","N");
//
// PrintNamed.printNamedMatrix(Out.out, infoMatrix);
// NamedMatrix cov = sign.getExpetedCovariance(allPars, config,"U2","E0","N");
//
// PrintWriter onComplete = Global.onComplete();
//
// printNamedMatrix(out, cov);
//
// cov = sign.getExpetedCovariance(allPars, config,"U2","E0","N","X");
//
// printNamedMatrix(out, cov);
//PlotPlugin pm = new PlotPlugin();
Map<String, UnivariateFunction> functions = new HashMap<>();
functions.put("U2", sign.getSignificanceFunction(allPars, "U2", "U2"));
// functions.put("UX", sign.getSignificanceFunction(allPars, "U2", "X"));
functions.put("X", sign.getSignificanceFunction(allPars, "X", "X"));
functions.put("trap", sign.getSignificanceFunction(allPars, "trap", "trap"));
functions.put("E0", sign.getSignificanceFunction(allPars, "E0", "E0"));
MetaBuilder builder = new MetaBuilder("significance");
builder.putValue("from", 14000d);
builder.putValue("to", 18500d);
pm.plotFunction(builder.build(), functions);
// printFuntionSimple(out(), func, 14000d, 18600d, 200);

View File

@ -1,102 +0,0 @@
/*
* 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.Global
import hep.dataforge.io.ColumnedDataWriter
import hep.dataforge.meta.Meta
import hep.dataforge.stat.fit.FitHelper
import hep.dataforge.stat.fit.FitResult
import hep.dataforge.stat.fit.FitStage
import hep.dataforge.stat.fit.ParamSet
import hep.dataforge.stat.models.XYModel
import inr.numass.NumassPlugin
import inr.numass.data.SpectrumAdapter
import inr.numass.data.SpectrumGenerator
import inr.numass.models.NBkgSpectrum
import inr.numass.models.sterile.SterileNeutrinoSpectrum
import inr.numass.utils.DataModelUtils
import static java.util.Locale.setDefault
/**
*
* @author Darksnake
*/
setDefault(Locale.US);
new NumassPlugin().startGlobal()
SterileNeutrinoSpectrum sp = new SterileNeutrinoSpectrum(Global.INSTANCE, Meta.empty());
//beta.setCaching(false);
NBkgSpectrum spectrum = new NBkgSpectrum(sp);
XYModel model = new XYModel(Meta.empty(), new SpectrumAdapter(Meta.empty()), spectrum);
ParamSet allPars = new ParamSet();
allPars.setParValue("N", 2e6 / 100);
//значение 6е-6 соответствует полной интенстивности 6е7 распадов в секунду
//Проблема была в переполнении счетчика событий в генераторе. Заменил на long. Возможно стоит поставить туда число с плавающей точкой
allPars.setParError("N", 6);
allPars.setParDomain("N", 0d, Double.POSITIVE_INFINITY);
allPars.setParValue("bkg", 2d);
allPars.setParError("bkg", 0.03);
allPars.setParValue("E0", 18575.0);
allPars.setParError("E0", 2);
allPars.setParValue("mnu2", 0d);
allPars.setParError("mnu2", 1d);
allPars.setParValue("msterile2", 8000 * 8000);
allPars.setParValue("U2", 0);
allPars.setParError("U2", 1e-4);
allPars.setParDomain("U2", -1d, 1d);
allPars.setParValue("X", 0);
allPars.setParError("X", 0.01);
allPars.setParDomain("X", 0d, Double.POSITIVE_INFINITY);
allPars.setParValue("trap", 0d);
allPars.setParError("trap", 0.01d);
allPars.setParDomain("trap", 0d, Double.POSITIVE_INFINITY);
// PrintNamed.printSpectrum(Global.onComplete(), spectrum, allPars, 0.0, 18700.0, 600);
//String fileName = "d:\\PlayGround\\merge\\scans.onComplete";
// String configName = "d:\\PlayGround\\SCAN.CFG";
// ListTable config = OldDataReader.readConfig(configName);
SpectrumGenerator generator = new SpectrumGenerator(model, allPars, 12316);
def data = generator.generateData(DataModelUtils.getUniformSpectrumConfiguration(12000, 18500, 604800 / 100 * 100, 130));
//data = TritiumUtils.correctForDeadTime(data, new SpectrumAdapter(), 10e-9);
// data = data.filter("X", Value.of(15510.0), Value.of(18610.0));
// allPars.setParValue("X", 0.4);
ColumnedDataWriter.writeTable(System.out, data, "--- DATA ---");
//FitState state = new FitState(data, model, allPars);
////new PlotFitResultAction().eval(state);
//
//
//def res = fm.runStage(state, "QOW", FitStage.TASK_RUN, "N", "bkg", "E0", "U2");
FitResult res = new FitHelper(Global.INSTANCE).fit(data)
.model(model)
.params(allPars)
.stage("QOW", FitStage.TASK_RUN, "N", "E0")
.run();
//
//
//
res.printState(new PrintWriter(System.out));

View File

@ -1,67 +0,0 @@
/*
* 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.Global
import hep.dataforge.io.FittingIOUtils
import hep.dataforge.stat.fit.FitManager
import hep.dataforge.stat.fit.ParamSet
import hep.dataforge.stat.models.XYModel
import inr.numass.data.SpectrumAdapter
import inr.numass.models.GunSpectrum
import inr.numass.models.NBkgSpectrum
import static java.util.Locale.setDefault
setDefault(Locale.US);
Global global = Global.instance();
// global.loadModule(new MINUITModule());
FitManager fm = new FitManager();
GunSpectrum gsp = new GunSpectrum();
NBkgSpectrum spectrum = new NBkgSpectrum(gsp);
XYModel model = new XYModel("gun", spectrum, new SpectrumAdapter());
ParamSet allPars = new ParamSet()
.setPar("N", 1e3, 1e2)
.setPar("pos", 18500, 0.1)
.setPar("bkg", 50, 1)
.setPar("resA", 5.3e-5, 1e-5)
.setPar("sigma", 0.3, 0.03);
PrintNamed.printSpectrum(new PrintWriter(System.out), spectrum, allPars, 18495, 18505, 100);
allPars.setParValue("sigma", 0.6);
FittingIOUtils.printSpectrum(new PrintWriter(System.out), spectrum, allPars, 18495, 18505, 100);
// //String fileName = "d:\\PlayGround\\merge\\scans.onComplete";
//// String configName = "d:\\PlayGround\\SCAN.CFG";
//// ListTable config = OldDataReader.readConfig(configName);
// SpectrumGenerator generator = new SpectrumGenerator(model, allPars, 12316);
//
// ListTable data = generator.generateData(DataModelUtils.getUniformSpectrumConfiguration(18495, 18505, 20, 20));
//
//// data = data.filter("X", Value.of(15510.0), Value.of(18610.0));
//// allPars.setParValue("X", 0.4);
// FitState state = FitTaskManager.buildState(data, model, allPars);
//
// FitState res = fm.runStage(state, "QOW", FitTask.TASK_RUN, "N", "bkg", "pos", "sigma");
//
// res.print(onComplete());

View File

@ -1,157 +0,0 @@
///*
// * To change this license header, choose License Headers in Project Properties.
// * To change this template file, choose Tools | Templates
// * and open the template in the editor.
// */
//
//package inr.numass.scripts
//
//import hep.dataforge.grind.Grind
//import hep.dataforge.values.Values
//import inr.numass.NumassUtils
//import inr.numass.data.api.NumassPoint
//import inr.numass.data.storage.NumassDataLoader
//
//import inr.numass.data.PileUpSimulator
//import inr.numass.utils.UnderflowCorrection
//import org.apache.commons.math3.random.JDKRandomGenerator
//
//rnd = new JDKRandomGenerator();
//
//////Loading data
//File dataDir = new File("D:\\Work\\Numass\\data\\2016_10\\Fill_1\\set_28")
////File dataDir = new File("D:\\Work\\Numass\\data\\2016_10\\Fill_2_wide\\set_7")
//if (!dataDir.exists()) {
// println "dataDir directory does not exist"
//}
//def data = NumassDataLoader.fromLocalDir(null, dataDir).getNMPoints()
//
////File rootDir = new File("D:\\Work\\Numass\\data\\2016_10\\Fill_1")
//////File rootDir = new File("D:\\Work\\Numass\\data\\2016_10\\Fill_2_wide")
//////File rootDir = new File("D:\\Work\\Numass\\data\\2017_01\\Fill_2_wide")
////
////NumassStorage storage = NumassStorage.buildLocalNumassRoot(rootDir, true);
////
////Collection<NMPoint> data = NumassDataUtils.joinSpectra(
//// StorageUtils.loaderStream(storage)
//// .filter { it.key.matches("set_3.") }
//// .map {
//// println "loading ${it.key}"
//// it.value
//// }
////)
//
////Simulation process
//Map<String, List<NumassPoint>> res = [:]
//
//List<NumassPoint> generated = new ArrayList<>();
//List<NumassPoint> registered = new ArrayList<>();
//List<NumassPoint> firstIteration = new ArrayList<>();
//List<NumassPoint> secondIteration = new ArrayList<>();
//List<NumassPoint> pileup = new ArrayList<>();
//
//lowerChannel = 400;
//upperChannel = 1800;
//
//PileUpSimulator buildSimulator(NumassPoint point, double cr, NumassPoint reference = null, boolean extrapolate = true, double scale = 1d) {
// def cfg = Grind.buildMeta(cr: cr) {
// pulser(mean: 3450, sigma: 86.45, freq: 66.43)
// }
// //NMEventGeneratorWithPulser generator = new NMEventGeneratorWithPulser(rnd, cfg)
//
// def generator = b
//
// if (extrapolate) {
// double[] chanels = new double[RawNMPoint.MAX_CHANEL];
// double[] values = new double[RawNMPoint.MAX_CHANEL];
// Values fitResult = new UnderflowCorrection().fitPoint(point, 400, 600, 1800, 20);
//
// def amp = fitResult.getDouble("amp")
// def sigma = fitResult.getDouble("expConst")
// if (sigma > 0) {
//
// for (int i = 0; i < upperChannel; i++) {
// chanels[i] = i;
// if (i < lowerChannel) {
// values[i] = point.getLength()*amp * Math.exp((i as double) / sigma)
// } else {
// values[i] = Math.max(0, point.getCount(i) - (reference == null ? 0 : reference.getCount(i)) as int);
// }
// }
// generator.loadSpectrum(chanels, values)
// } else {
// generator.loadSpectrum(point, reference, lowerChannel, upperChannel);
// }
// } else {
// generator.loadSpectrum(point, reference, lowerChannel, upperChannel);
// }
//
// return new PileUpSimulator(point.length * scale, rnd, generator).withUset(point.voltage).generate();
//}
//
//double adjustCountRate(PileUpSimulator simulator, NumassPoint point) {
// double generatedInChannel = simulator.generated().getCountInWindow(lowerChannel, upperChannel);
// double registeredInChannel = simulator.registered().getCountInWindow(lowerChannel, upperChannel);
// return (generatedInChannel / registeredInChannel) * (point.getCountInWindow(lowerChannel, upperChannel) / point.getLength());
//}
//
//data.forEach { point ->
// double cr = NumassUtils.countRateWithDeadTime(point, lowerChannel, upperChannel, 6.55e-6);
//
// PileUpSimulator simulator = buildSimulator(point, cr);
//
// //second iteration to exclude pileup overlap
// NumassPoint pileupPoint = simulator.pileup();
// firstIteration.add(simulator.registered());
//
// //updating count rate
// cr = adjustCountRate(simulator, point);
// simulator = buildSimulator(point, cr, pileupPoint);
//
// pileupPoint = simulator.pileup();
// secondIteration.add(simulator.registered());
//
// cr = adjustCountRate(simulator, point);
// simulator = buildSimulator(point, cr, pileupPoint);
//
// generated.add(simulator.generated());
// registered.add(simulator.registered());
// pileup.add(simulator.pileup());
//}
//res.put("original", data);
//res.put("generated", generated);
//res.put("registered", registered);
//// res.put("firstIteration", new SimulatedPoint("firstIteration", firstIteration));
//// res.put("secondIteration", new SimulatedPoint("secondIteration", secondIteration));
//res.put("pileup", pileup);
//
//def keys = res.keySet();
//
////print spectra for selected point
//double u = 16500d;
//
//List<Map> points = res.values().collect { it.find { it.voltage == u }.getMap(20, true) }
//
//println "\n Spectrum example for U = ${u}\n"
//
//print "channel\t"
//println keys.join("\t")
//
//points.first().keySet().each {
// print "${it}\t"
// println points.collect { map -> map[it] }.join("\t")
//}
//
////printing count rate in window
//print "U\tLength\t"
//print keys.collect { it + "[total]" }.join("\t") + "\t"
//print keys.collect { it + "[pulse]" }.join("\t") + "\t"
//println keys.join("\t")
//
//for (int i = 0; i < data.size(); i++) {
// print "${data.get(i).getVoltage()}\t"
// print "${data.get(i).getLength()}\t"
// print keys.collect { res[it].get(i).getTotalCount() }.join("\t") + "\t"
// print keys.collect { res[it].get(i).getCountInWindow(3100, 3800) }.join("\t") + "\t"
// println keys.collect { res[it].get(i).getCountInWindow(400, 3100) }.join("\t")
//}

View File

@ -1,81 +0,0 @@
/*
* 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.Global
import hep.dataforge.grind.GrindMetaBuilder
import hep.dataforge.io.FittingIOUtils
import hep.dataforge.meta.Meta
import hep.dataforge.stat.fit.ParamSet
import hep.dataforge.stat.models.XYModel
import hep.dataforge.stat.parametric.ParametricFunction
import inr.numass.data.SpectrumAdapter
import inr.numass.models.NBkgSpectrum
import inr.numass.models.sterile.SterileNeutrinoSpectrum
import static java.util.Locale.setDefault
/**
*
* @author Darksnake
*/
setDefault(Locale.US);
//ParametricFunction beta = new BetaSpectrum();
//ModularSpectrum beta = new ModularSpectrum(new BetaSpectrum(), 8.3e-5, 13990d, 18600d);
//beta.setCaching(false)
Meta cfg = new GrindMetaBuilder().meta() {
resolution(width: 8.3e-5)
}.build();
ParametricFunction beta = new SterileNeutrinoSpectrum(Global.instance(), cfg);
NBkgSpectrum spectrum = new NBkgSpectrum(beta);
XYModel model = new XYModel(spectrum, new SpectrumAdapter());
ParamSet allPars = new ParamSet();
allPars.setPar("N", 6.6579e+05, 1.8e+03, 0d, Double.POSITIVE_INFINITY);
allPars.setPar("bkg", 0.5387, 0.050);
allPars.setPar("E0", 18574.94, 1.4);
allPars.setPar("mnu2", 0d, 1d);
allPars.setPar("msterile2", 1000d * 1000d, 0);
allPars.setPar("U2", 0.0, 1e-4, -1d, 1d);
allPars.setPar("X", 0.04, 0.01, 0d, Double.POSITIVE_INFINITY);
allPars.setPar("trap", 1.634, 0.01, 0d, Double.POSITIVE_INFINITY);
FittingIOUtils.printSpectrum(Global.out(), spectrum, allPars, 14000, 18600.0, 400);
//SpectrumGenerator generator = new SpectrumGenerator(model, allPars, 12316);
//
//ListTable data = generator.generateData(DataModelUtils.getUniformSpectrumConfiguration(14000d, 18500, 2000, 90));
//
//data = NumassUtils.correctForDeadTime(data, new SpectrumAdapter(), 1e-8);
//// data = data.filter("X", Value.of(15510.0), Value.of(18610.0));
//// allPars.setParValue("X", 0.4);
//FitState state = new FitState(data, model, allPars);
////new PlotFitResultAction().eval(state);
//
//
//FitState res = fm.runStage(state, "QOW", FitTask.TASK_RUN, "N", "bkg", "E0", "U2", "trap");
//
//
//
//res.print(onComplete());
//

View File

@ -1,101 +0,0 @@
/*
* 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.stat.fit.FitManager
import hep.dataforge.stat.fit.FitState
import hep.dataforge.stat.fit.MINUITPlugin
import hep.dataforge.stat.fit.ParamSet
import hep.dataforge.stat.models.XYModel
import hep.dataforge.tables.ListTable
import inr.numass.data.SpectrumAdapter
import inr.numass.data.SpectrumGenerator
import inr.numass.models.BetaSpectrum
import inr.numass.models.ModularSpectrum
import inr.numass.models.NBkgSpectrum
import inr.numass.models.ResolutionFunction
import inr.numass.utils.DataModelUtils
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());
resolution.setTailFunction(ResolutionFunction.getAngledTail(0.00325));
ModularSpectrum beta = new ModularSpectrum(new BetaSpectrum(), resolution, 18395d, 18580d);
beta.setCaching(false);
NBkgSpectrum spectrum = new NBkgSpectrum(beta);
XYModel model = new XYModel("tritium", spectrum, new SpectrumAdapter());
ParamSet allPars = new ParamSet();
allPars.setPar("N", 6e9, 1e5, 0, Double.POSITIVE_INFINITY);
allPars.setPar("bkg", 0.002, 0.005 );
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(18530d, 18580, 1e7, 60)
//def config = DataModelUtils.getSpectrumConfigurationFromResource("/data/run23.cfg")
ListTable data = generator.generateExactData(config);
FitState state = new FitState(data, model, allPars);
println("Simulating data with real tail")
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

@ -1,99 +0,0 @@
/*
* 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.Global
import hep.dataforge.io.FittingIOUtils
import hep.dataforge.meta.Meta
import hep.dataforge.stat.fit.*
import hep.dataforge.stat.models.XYModel
import inr.numass.data.SpectrumAdapter
import inr.numass.data.SpectrumGenerator
import inr.numass.models.BetaSpectrum
import inr.numass.models.ModularSpectrum
import inr.numass.models.NBkgSpectrum
import inr.numass.models.ResolutionFunction
import inr.numass.utils.DataModelUtils
import org.apache.commons.math3.analysis.BivariateFunction
import static java.util.Locale.setDefault
/**
*
* @author Darksnake
*/
setDefault(Locale.US);
new MINUITPlugin().startGlobal();
FitManager fm = Global.INSTANCE.get(FitManager)
BivariateFunction resolution = new ResolutionFunction(8.3e-5);
ModularSpectrum beta = new ModularSpectrum(new BetaSpectrum(), resolution, 13490d, 18575d);
beta.setCaching(false);
NBkgSpectrum spectrum = new NBkgSpectrum(beta);
XYModel model = new XYModel(Meta.empty(), new SpectrumAdapter(Meta.empty()), spectrum);
ParamSet allPars = new ParamSet();
allPars.setPar("N", 6e5, 10, 0, Double.POSITIVE_INFINITY);
allPars.setPar("bkg", 2d, 0.1);
allPars.setPar("E0", 18575.0, 0.05);
allPars.setPar("mnu2", 0, 1);
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", 0, 0.01, 0d, Double.POSITIVE_INFINITY);
SpectrumGenerator generator = new SpectrumGenerator(model, allPars, 12316);
def data = generator.generateData(DataModelUtils.getUniformSpectrumConfiguration(14000d, 18200, 1e6, 60));
// data = data.filter("X", Value.of(15510.0), Value.of(18610.0));
allPars.setParValue("U2", 0);
FitState state = new FitState(data, model, allPars);
//new PlotFitResultAction(Global.instance(), null).runOne(state);
//double delta = 4e-6;
//resolution.setTailFunction{double E, double U ->
// 1-delta*(E-U);
//}
//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));
def res = fm.runStage(state, "QOW", FitStage.TASK_RUN, "N", "bkg", "E0", "U2", "trap");
res.printState(new PrintWriter(System.out))
FittingIOUtils.printResiduals(new PrintWriter(System.out), res.optState().get())

View File

@ -1,68 +0,0 @@
/*
* 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.io.PrintFunction
import hep.dataforge.maths.integration.UnivariateIntegrator
import hep.dataforge.plots.data.XYFunctionPlot
import hep.dataforge.plots.jfreechart.JFreeChartFrame
import hep.dataforge.stat.fit.ParamSet
import inr.numass.models.ExperimentalVariableLossSpectrum
import org.apache.commons.math3.analysis.UnivariateFunction
//double exPos = 12.94
//double exW = 1.31
//double ionPos = 14.13
//double ionW = 12.79
//double exIonRatio = 0.6059
ParamSet params = new ParamSet()
.setParValue("shift",0)
.setParValue("X", 0.4)
.setParValue("exPos", 12.94)
.setParValue("ionPos", 15.6)
.setParValue("exW", 1.31)
.setParValue("ionW", 12.79)
.setParValue("exIonRatio", 0.6059)
ExperimentalVariableLossSpectrum lsp = new ExperimentalVariableLossSpectrum(19005, 8e-5, 19010,0.2);
JFreeChartFrame frame = JFreeChartFrame.drawFrame("Experimental Loss Test", null);
UnivariateIntegrator integrator = NumassContext.defaultIntegrator
UnivariateFunction exFunc = lsp.excitation(params.getValue("exPos"), params.getValue("exW"));
frame.add(XYFunctionPlot.plotFunction("ex", exFunc, 0d, 50d, 500));
println "excitation norm factor " + integrator.integrate(0, 50, exFunc)
UnivariateFunction ionFunc = lsp.ionization(params.getValue("ionPos"), params.getValue("ionW"));
frame.add(XYFunctionPlot.plotFunction("ion", ionFunc, 0d, 50d, 500));
println "ionization norm factor " + integrator.integrate(0, 200, ionFunc)
UnivariateFunction sumFunc = lsp.singleScatterFunction(params);
frame.add(XYFunctionPlot.plotFunction("sum", sumFunc, 0d, 50d, 500));
println "sum norm factor " + integrator.integrate(0, 100, sumFunc)
PrintFunction.printFunctionSimple(new PrintWriter(System.out), sumFunc, 0d, 50d, 100)
JFreeChartFrame integerFrame = JFreeChartFrame.drawFrame("Experimental Loss Test", null);
UnivariateFunction integr = { d-> lsp.value(d,params)}
integerFrame.add(XYFunctionPlot.plotFunction("integr", integr, 18950d, 19005d, 500));

View File

@ -1,46 +0,0 @@
/*
* 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.plots.data.XYFunctionPlot
import hep.dataforge.plots.jfreechart.JFreeChartFrame
import org.apache.commons.math3.analysis.UnivariateFunction
def lorenz = {x, x0, gama -> 1/(3.14*gama*(1+(x-x0)*(x-x0)/gama/gama))}
def excitationSpectrum = {Map<Double,Double> lines, double gama ->
UnivariateFunction function = {x->
double res = 0;
lines.each{k,v -> res += lorenz(x,k,gama)*v};
return res;
}
return function;
}
def lines =
[
12.6:0.5,
12.4:0.3,
12.2:0.2
]
UnivariateFunction excitation = excitationSpectrum(lines,0.08)
JFreeChartFrame frame = JFreeChartFrame.drawFrame("theoretical loss spectrum", null);
frame.add(XYFunctionPlot.plotFunction("excitation", excitation, 0d, 20d, 500));

View File

@ -1,119 +0,0 @@
/*
* 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.Global
import hep.dataforge.data.DataSet
import hep.dataforge.stat.fit.FitManager
import hep.dataforge.stat.fit.FitState
import hep.dataforge.stat.fit.ParamSet
import hep.dataforge.stat.likelihood.BayesianConfidenceLimit
import hep.dataforge.stat.models.XYModel
import hep.dataforge.tables.ListTable
import inr.numass.data.SpectrumGenerator
import inr.numass.models.BetaSpectrum
import inr.numass.models.ModularSpectrum
import inr.numass.models.NBkgSpectrum
import static inr.numass.utils.DataModelUtils.getUniformSpectrumConfiguration
PrintWriter out = Global.out();
FitManager fm = new FitManager();
setSeed(543982);
// TritiumSpectrum beta = new TritiumSpectrum(2e-4, 13995d, 18580d);
File fssfile = new File("c:\\Users\\Darksnake\\Dropbox\\PlayGround\\FS.txt");
ModularSpectrum beta = new ModularSpectrum(new BetaSpectrum(),8.3e-5, 14400d, 19010d);
beta.setCaching(false);
NBkgSpectrum spectrum = new NBkgSpectrum(beta);
XYModel model = new XYModel("tritium", spectrum);
ParamSet allPars = new ParamSet();
allPars.setParValue("N", 6e5);
//значение 6е-6 соответствует полной интенстивности 6е7 распадов в секунду
//Проблема была в переполнении счетчика событий в генераторе. Заменил на long. Возможно стоит поставить туда число с плавающей точкой
allPars.setParError("N", 25);
allPars.setParDomain("N", 0d, Double.POSITIVE_INFINITY);
allPars.setParValue("bkg", 5);
allPars.setParError("bkg", 1e-3);
allPars.setParValue("E0", 18575d);
allPars.setParError("E0", 0.1);
allPars.setParValue("mnu2", 0d);
allPars.setParError("mnu2", 1d);
allPars.setParValue("msterile2", 1000 * 1000);
allPars.setParValue("U2", 0);
allPars.setParError("U2", 1e-4);
allPars.setParDomain("U2", 0d, 1d);
allPars.setParValue("X", 0.0);
allPars.setParDomain("X", 0d, Double.POSITIVE_INFINITY);
allPars.setParValue("trap", 1d);
allPars.setParError("trap", 0.01d);
allPars.setParDomain("trap", 0d, Double.POSITIVE_INFINITY);
// PlotPlugin pm = new PlotPlugin();
// String plotTitle = "Tritium spectrum";
// pm.plotFunction(ParametricUtils.getSpectrumFunction(spectrum, allPars), 14000, 18600, 500,plotTitle, null);
// PrintNamed.printSpectrum(Out.onComplete, beta.trapping, allPars, 14000d, 18600d, 500);
// double e = 18570d;
// trans.alpha = 1e-4;
// trans.plotTransmission(System.onComplete, allPars, e, e-1000d, e+100d, 200);
SpectrumGenerator generator = new SpectrumGenerator(model, allPars);
// ColumnedDataFile file = new ColumnedDataFile("d:\\PlayGround\\RUN36.cfg");
// ListTable config = file.getPoints("time","X");
double Elow = 14000d;
double Eup = 18600d;
int numpoints = (int) ((Eup - Elow) / 50);
double time = 1e6 / numpoints; // 3600 / numpoints;
DataSet config = getUniformSpectrumConfiguration(Elow, Eup, time, numpoints);
// config.addAll(DataModelUtils.getUniformSpectrumConfiguration(Eup, Elow, time, numpoints));// в обратную сторону
ListTable data = generator.generateData(config);
// plotTitle = "Generated tritium spectrum data";
// pm.plotXYScatter(data, "X", "Y",plotTitle, null);
// bareBeta.setFSS("D:\\PlayGround\\FSS.dat");
// data = tritiumUtils.applyDrift(data, 2.8e-6);
FitState state = fm.buildState(data, model, allPars);
// fm.checkDerivs();
// res.print(Out.onComplete);
// fm.checkFitDerivatives();
FitState res = fm.runDefaultStage(state, "U2", "N", "trap");
res.print(out);
// res = fm.runFrom(res);
// res = fm.generateErrorsFrom(res);
beta.setCaching(true);
beta.setSuppressWarnings(true);
BayesianConfidenceLimit bm = new BayesianConfidenceLimit();
// bm.setPriorProb(new OneSidedUniformPrior("trap", 0, true));
// bm.setPriorProb(new Gaussian("trap", 1d, 0.002));
// bm.printMarginalLikelihood(Out.onComplete,"U2", res);
FitState conf = bm.getConfidenceInterval("U2", res, ["U2", "N", "trap"]);
// plotTitle = String.format("Marginal likelihood for parameter \'%s\'", "U2");
// pm.plotFunction(bm.getMarginalLikelihood("U2", res), 0, 2e-3, 40,plotTitle, null);
conf.print(out);
// PrintNamed.printLogProbRandom(Out.onComplete, res, 5000,0.5d, "E0","N");
spectrum.counter.print(out);

View File

@ -1,100 +0,0 @@
package inr.numass.scripts.models
import hep.dataforge.context.Context
import hep.dataforge.context.Global
import hep.dataforge.grind.GrindShell
import hep.dataforge.grind.helpers.PlotHelper
import hep.dataforge.meta.Meta
import hep.dataforge.stat.fit.FitManager
import hep.dataforge.stat.fit.FitStage
import hep.dataforge.stat.fit.FitState
import hep.dataforge.stat.fit.ParamSet
import hep.dataforge.stat.models.XYModel
import hep.dataforge.stat.parametric.ParametricFunction
import hep.dataforge.tables.Table
import inr.numass.NumassPlugin
import inr.numass.data.SpectrumAdapter
import inr.numass.data.SpectrumGenerator
import inr.numass.models.NBkgSpectrum
import inr.numass.models.NumassModelsKt
import inr.numass.models.misc.ModGauss
import inr.numass.models.sterile.NumassBeta
import inr.numass.utils.DataModelUtils
Context ctx = Global.instance()
ctx.getPlugins().load(NumassPlugin)
new GrindShell(ctx).eval {
PlotHelper ph = plots
def beta = new NumassBeta().getSpectrum(0)
def response = new ModGauss(5.0)
ParametricFunction spectrum = NumassModelsKt.convolute(beta, response)
def model = new XYModel(Meta.empty(), new SpectrumAdapter(Meta.empty()), new NBkgSpectrum(spectrum))
ParamSet params = morph(ParamSet, [:], "params") {
N(value: 1e+14, err: 30, lower: 0)
bkg(value: 5.0, err: 0.1)
E0(value: 18575.0, err: 0.1)
mnu2(value: 0, err: 0.01)
msterile2(value: 7000**2, err: 1)
U2(value: 0.0, err: 1e-3)
//X(value: 0.0, err: 0.01, lower: 0)
//trap(value: 1.0, err: 0.05)
w(value: 150, err: 5)
//shift(value: 1, err: 1e-2)
//tailAmp(value: 0.01, err: 1e-2)
tailW(value: 300, err: 1)
}
// double norm = NumassIntegrator.defaultIntegrator.integrate(1000d, 18500d) {
// model.value(it, params)
// }
// println("The total number of events is $norm")
//
// ph.plotFunction(-2000d, 500d, 400, "actual", "response") { double x ->
// response.value(x, params)
// }
SpectrumGenerator generator = new SpectrumGenerator(model, params, 12316);
ph.plot(data: (2000..19500).step(50).collectEntries {
[it, model.value(it, params)]
}, name: "spectrum", frame: "test")
.configure(showLine: true, showSymbol: false, showErrors: false, thickness: 2, connectionType: "spline", color: "red")
Table data = generator.generateData(DataModelUtils.getUniformSpectrumConfiguration(7000, 19500, 1, 1000));
//params.setParValue("w", 151)
//params.setParValue("tailAmp", 0.011)
//params.setParValue("X", 0.01)
//params.setParValue("trap", 0.01)
//params.setParValue("mnu2", 4)
ph.plotFunction(-2000d, 500d, 400, "supposed", "response") { double x ->
response.value(x, params)
}
ph.plot(data: (2000..19500).step(50).collectEntries {
[it, model.value(it, params)]
}, name: "spectrum-mod", frame: "test")
.configure(showLine: true, showSymbol: false, showErrors: false, thickness: 2, connectionType: "spline", color: "green")
ph.plot(data: data, frame: "test", adapter: new SpectrumAdapter(Meta.empty()))
.configure(color: "blue")
FitState state = new FitState(data, model, params);
def fm = ctx.get(FitManager)
def res = fm.runStage(state, "MINUIT", FitStage.TASK_RUN, "N", "bkg", "E0", "U2");
res.printState(ctx.getOutput.out().newPrintWriter());
NumassIOKt.display(res, ctx, "fit")
}

View File

@ -1,66 +0,0 @@
package inr.numass.scripts.underflow
import hep.dataforge.cache.CachePlugin
import hep.dataforge.context.Context
import hep.dataforge.context.Global
import hep.dataforge.data.DataNode
import hep.dataforge.grind.GrindShell
import hep.dataforge.io.ColumnedDataWriter
import hep.dataforge.meta.Meta
import hep.dataforge.tables.ColumnTable
import hep.dataforge.tables.Table
import inr.numass.NumassPlugin
import inr.numass.data.NumassDataUtils
import static hep.dataforge.grind.Grind.buildMeta
Context ctx = Global.instance()
ctx.getPlugins().load(FXPlotManager)
ctx.getPlugins().load(NumassPlugin.class)
ctx.getPlugins().load(CachePlugin.class)
Meta meta = buildMeta {
data(dir: "D:\\Work\\Numass\\data\\2017_05\\Fill_2", mask: "set_.{1,3}")
generate(t0: 3e4, sort: false)
}
def shell = new GrindShell(ctx);
DataNode<Table> spectra = UnderflowUtils.getSpectraMap(shell, meta);
shell.eval {
def columns = [:];
Map<Integer, Table> binned = [:]
(14500..17500).step(500).each {
Table up = binned.computeIfAbsent(it) { key ->
NumassDataUtils.spectrumWithBinning(spectra.optData(key as String).get().get(), 20, 400, 3100);
}
Table lo = binned.computeIfAbsent(it - 500) { key ->
NumassDataUtils.spectrumWithBinning(spectra.optData(key as String).get().get(), 20, 400, 3100);
}
columns << [channel: up.channel]
columns << [(it as String): NumassDataUtils.subtractSpectrum(lo, up).getColumn("cr")]
}
ColumnedDataWriter.writeTable(System.out, ColumnTable.of(columns), "Response function")
// println()
// println()
//
// columns.clear()
//
// binned.each { key, table ->
// columns << [channel: table.channel]
// columns << [(key as String): table.cr]
// }
//
// ColumnedDataWriter.writeTable(System.out,
// ColumnTable.of(columns),
// "Spectra")
}

View File

@ -1,139 +0,0 @@
/*
* To change this license header, choose License Headers in Project Properties.
* To change this template file, choose Tools | Templates
* and open the template in the editor.
*/
package inr.numass.scripts.underflow
import hep.dataforge.cache.CachePlugin
import hep.dataforge.context.Context
import hep.dataforge.context.Global
import hep.dataforge.data.DataNode
import hep.dataforge.grind.GrindShell
import hep.dataforge.grind.helpers.PlotHelper
import hep.dataforge.io.ColumnedDataWriter
import hep.dataforge.meta.Meta
import hep.dataforge.plots.PlotGroup
import hep.dataforge.plots.data.DataPlot
import hep.dataforge.tables.Adapters
import hep.dataforge.tables.Table
import inr.numass.NumassPlugin
import inr.numass.data.analyzers.NumassAnalyzerKt
import inr.numass.subthreshold.Threshold
import javafx.application.Platform
import static hep.dataforge.grind.Grind.buildMeta
import static inr.numass.data.analyzers.NumassAnalyzer.CHANNEL_KEY
import static inr.numass.data.analyzers.NumassAnalyzer.COUNT_RATE_KEY
Context ctx = Global.instance()
ctx.getPlugins().load(NumassPlugin)
ctx.getPlugins().load(CachePlugin)
Meta meta = buildMeta(t0: 3e4) {
data(dir: "D:\\Work\\Numass\\data\\2017_11\\Fill_2", mask: "set_3.")
subtract(reference: 18500)
fit(xLow: 400, xHigh: 600, upper: 3000, binning: 20, method: "pow")
scan(
hi: [500, 600, 700, 800, 900],
lo: [350, 400, 450],
def: 700
)
window(lo: 300, up: 3000)
}
def shell = new GrindShell(ctx);
DataNode<Table> spectra = Threshold.INSTANCE.getSpectraMap(ctx, meta).computeAll();
shell.eval {
//subtracting reference point
Map<Double, Table> spectraMap
if (meta.hasValue("subtract.reference")) {
String referenceVoltage = meta["subtract.reference"]
println "subtracting reference point ${referenceVoltage}"
def referencePoint = spectra.get(referenceVoltage)
spectraMap = spectra
.findAll { it.name != referenceVoltage }
.collectEntries {
[(it.meta["voltage"]): NumassAnalyzerKt.subtractAmplitudeSpectrum(it.get(), referencePoint)]
}
} else {
spectraMap = spectra.collectEntries { return [(it.meta["voltage"]): it.get()] }
}
//Showing selected points
def showPoints = { Map points, int binning = 20, int loChannel = 300, int upChannel = 2000 ->
def plotGroup = new PlotGroup("points");
def adapter = Adapters.buildXYAdapter(CHANNEL_KEY, COUNT_RATE_KEY)
points.each {
plotGroup.set(
DataPlot.plot(
it.key as String,
adapter,
NumassAnalyzerKt.withBinning(it.value as Table, binning)
)
)
}
//configuring and plotting
plotGroup.configure(showLine: true, showSymbol: false, showErrors: false, connectionType: "step")
def frame = (plots as PlotHelper).getManager().getPlotFrame("Spectra")
frame.configureValue("yAxis.type", "log")
frame.add(plotGroup)
}
showPoints(spectraMap.findAll { it.key in [14100d, 14200d, 14300d, 14400d, 14800d, 15000d, 15200d, 15700d] })
meta["scan.hi"].each { xHigh ->
println "Caclculate correctuion for upper linearity bound: ${xHigh}"
Table correctionTable = TableTransform.filter(
Threshold.INSTANCE.calculateSubThreshold(
spectraMap,
meta.getMeta("fit").builder.setValue("xHigh", xHigh)
),
"correction",
0,
2
)
if (xHigh == meta["scan.def"]) {
ColumnedDataWriter.writeTable(System.out, correctionTable, "underflow parameters")
}
Platform.runLater {
(plots as PlotHelper).plot(
data: correctionTable,
adapter: Adapters.buildXYAdapter("U", "correction"),
name: "upper_${xHigh}",
frame: "upper"
)
}
}
meta["scan.lo"].each { xLow ->
println "Caclculate correctuion for lower linearity bound: ${xLow}"
Table correctionTable = TableTransform.filter(
Threshold.INSTANCE.calculateSubThreshold(
spectraMap,
meta.getMeta("fit").builder.setValue("xLow", xLow)
),
"correction",
0,
2
)
Platform.runLater {
(plots as PlotHelper).plot(
data: correctionTable,
adapter: Adapters.buildXYAdapter("U", "correction"),
name: "lower_${xLow}",
frame: "lower"
)
}
}
}

View File

@ -1,17 +1,8 @@
package inr.numass.scripts.analysis
import hep.dataforge.context.Global
import hep.dataforge.context.plugin
import hep.dataforge.fx.output.FXOutputManager
import hep.dataforge.grind.workspace.GroovyWorkspaceParser
import hep.dataforge.plots.plotData
import hep.dataforge.tables.Adapters
import hep.dataforge.tables.Table
import hep.dataforge.workspace.FileBasedWorkspace
import hep.dataforge.workspace.Workspace
import hep.dataforge.workspace.context
import hep.dataforge.workspace.data
import inr.numass.displayChart
import java.io.File
fun main() {
@ -32,7 +23,7 @@ fun main() {
//workspace.context.setValue("cache.enabled", false)
workspace.runTask("fit","Fill_4")
workspace.runTask("fit", "Fill_4")
// val result = workspace.runTask("dif", "adiab_19").first().get() as Table
//

View File

@ -14,7 +14,7 @@ import kotlinx.coroutines.runBlocking
import java.io.File
private suspend fun createSummaryNode(storage: Storage): MetaBuilder {
Global.logger.info("Reading content of shelf {}", storage.fullName)
Global.logger.info("Reading content of shelf ${storage.fullName}")
val builder = MetaBuilder("shelf")
.setValue("name", storage.name)
@ -24,7 +24,7 @@ private suspend fun createSummaryNode(storage: Storage): MetaBuilder {
if(element is Storage && element.name.startsWith("Fill")){
builder.putNode(createSummaryNode(element))
} else if(element is NumassDataLoader){
Global.logger.info("Reading content of set {}", element.fullName)
Global.logger.info("Reading content of set ${element.fullName}")
val setBuilder = MetaBuilder("set")
.setValue("name", element.name)

View File

@ -2,7 +2,7 @@ plugins {
kotlin("jvm")
id("org.openjfx.javafxplugin")
//id("com.github.johnrengelman.shadow")
id("org.beryx.runtime") version "1.12.7"
id("org.beryx.runtime") version "1.13.0"
application
}

View File

@ -1,5 +1,9 @@
rootProject.name = "numass"
plugins {
id("org.gradle.toolchains.foojay-resolver-convention") version "0.7.0"
}
include(":dataforge-core") // core classes and data structures
include(":dataforge-core:dataforge-json") // json meta converter
//include(":dataforge-core:osgi") // osgi framework for easy plugin delivery