Update dependencies and build
This commit is contained in:
parent
c20ac1ffb8
commit
8e45062957
GPG Key ID:
B10A55DCC7B9AEBB
@ -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)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -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")
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -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"))
|
||||
|
||||
@ -19,7 +19,6 @@ application{
|
||||
|
||||
compileKotlin {
|
||||
kotlinOptions {
|
||||
jvmTarget = "1.8"
|
||||
javaParameters = true
|
||||
}
|
||||
}
|
||||
|
||||
@ -1,5 +1,4 @@
|
||||
|
||||
description = 'dataforge-minuit'
|
||||
dependencies {
|
||||
api project(':dataforge-stat')
|
||||
api project(':dataforge-stat')
|
||||
}
|
||||
2
gradle/wrapper/gradle-wrapper.properties
vendored
2
gradle/wrapper/gradle-wrapper.properties
vendored
@ -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
|
||||
|
||||
@ -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")
|
||||
|
||||
@ -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"
|
||||
}
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@ -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
|
||||
|
||||
@ -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"]}"
|
||||
}
|
||||
|
||||
@ -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";
|
||||
@ -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;
|
||||
@ -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))
|
||||
@ -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);
|
||||
@ -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)
|
||||
)
|
||||
}
|
||||
@ -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,"")
|
||||
@ -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()
|
||||
@ -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);
|
||||
|
||||
@ -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));
|
||||
|
||||
@ -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());
|
||||
|
||||
@ -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")
|
||||
//}
|
||||
@ -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());
|
||||
//
|
||||
@ -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}")
|
||||
@ -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())
|
||||
@ -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));
|
||||
@ -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));
|
||||
@ -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);
|
||||
|
||||
@ -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")
|
||||
}
|
||||
@ -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")
|
||||
}
|
||||
@ -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"
|
||||
)
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -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
|
||||
//
|
||||
|
||||
@ -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)
|
||||
|
||||
@ -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
|
||||
}
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
Loading…
Reference in New Issue
Block a user