refactor numass resolution
This commit is contained in:
parent
57e3439220
commit
d41ff5fefc
@ -136,9 +136,7 @@ class ListTable @JvmOverloads constructor(override val format: TableFormat, poin
|
||||
* @throws NamingException
|
||||
*/
|
||||
@Throws(NamingException::class)
|
||||
fun row(vararg values: Any): Builder {
|
||||
return row(ValueMap.of(format.namesAsArray(), *values))
|
||||
}
|
||||
fun row(vararg values: Any): Builder = row(ValueMap.of(format.namesAsArray(), *values))
|
||||
|
||||
fun row(values: ValueProvider): Builder {
|
||||
val names = format.namesAsArray()
|
||||
@ -146,17 +144,12 @@ class ListTable @JvmOverloads constructor(override val format: TableFormat, poin
|
||||
return row(ValueMap(map))
|
||||
}
|
||||
|
||||
fun row(vararg values: NamedValue): Builder {
|
||||
return row(ValueMap.of(*values))
|
||||
}
|
||||
fun row(vararg values: NamedValue): Builder = row(ValueMap.of(*values))
|
||||
|
||||
fun row(vararg values: Pair<String, Any>): Builder {
|
||||
return row(ValueMap.of(values.map { NamedValue.of(it.first, it.second) }))
|
||||
}
|
||||
fun row(vararg values: Pair<String, Any>): Builder =
|
||||
row(ValueMap.of(values.map { NamedValue.of(it.first, it.second) }))
|
||||
|
||||
fun row(map: Map<String, Any>): Builder {
|
||||
return row(ValueMap.ofMap(map))
|
||||
}
|
||||
fun row(map: Map<String, Any>): Builder = row(ValueMap.ofMap(map))
|
||||
|
||||
fun rows(points: Iterable<Values>): Builder {
|
||||
for (point in points) {
|
||||
@ -170,16 +163,12 @@ class ListTable @JvmOverloads constructor(override val format: TableFormat, poin
|
||||
return this
|
||||
}
|
||||
|
||||
fun build(): Table {
|
||||
return ListTable(format, points)
|
||||
}
|
||||
fun build(): ListTable = ListTable(format, points)
|
||||
|
||||
/**
|
||||
* Build table without points name check
|
||||
*/
|
||||
fun buildUnsafe(): Table {
|
||||
return ListTable(format, points, true)
|
||||
}
|
||||
fun buildUnsafe(): Table = ListTable(format, points, true)
|
||||
}
|
||||
|
||||
companion object {
|
||||
|
||||
@ -52,8 +52,6 @@ interface Table : NavigableValuesSource, MetaMorph {
|
||||
*/
|
||||
fun getColumn(name: String): Column
|
||||
|
||||
|
||||
@JvmDefault
|
||||
override fun toMeta(): Meta {
|
||||
val res = MetaBuilder("table")
|
||||
res.putNode("format", format.toMeta())
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@ -1,77 +0,0 @@
|
||||
package inr.numass.scripts
|
||||
|
||||
import hep.dataforge.context.Context
|
||||
import hep.dataforge.context.Global
|
||||
import hep.dataforge.grind.GrindShell
|
||||
import hep.dataforge.grind.helpers.PlotHelper
|
||||
import hep.dataforge.stat.fit.ParamSet
|
||||
import inr.numass.NumassPlugin
|
||||
import inr.numass.models.sterile.NumassResolution
|
||||
import inr.numass.models.sterile.SterileNeutrinoSpectrum
|
||||
|
||||
import static hep.dataforge.grind.Grind.buildMeta
|
||||
|
||||
Context ctx = Global.instance()
|
||||
ctx.getPlugins().load(FXPlotManager)
|
||||
ctx.getPlugins().load(NumassPlugin.class)
|
||||
|
||||
GrindShell shell = new GrindShell(ctx)
|
||||
|
||||
shell.eval {
|
||||
PlotHelper plot = plots
|
||||
|
||||
|
||||
ParamSet params = new ParamSet(buildMeta {
|
||||
N(value: 2.7e+06, 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: 1000**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)
|
||||
})
|
||||
|
||||
def meta1 = buildMeta {
|
||||
resolution(width: 8.3e-5, tail: "(0.99797 - 3.05346E-7*D - 5.45738E-10 * D**2 - 6.36105E-14 * D**3)")
|
||||
}
|
||||
|
||||
def meta2 = buildMeta {
|
||||
resolution(width: 8.3e-5, tail: "(0.99797 - 3.05346E-7*D - 5.45738E-10 * D**2 - 6.36105E-14 * D**3)*(1-5e-3*sqrt(E/1000))")
|
||||
}
|
||||
|
||||
def resolution1 = new NumassResolution(
|
||||
ctx,
|
||||
meta1.getMeta("resolution")
|
||||
)
|
||||
|
||||
def resolution2 = new NumassResolution(
|
||||
ctx,
|
||||
meta2.getMeta("resolution")
|
||||
)
|
||||
|
||||
plot.plot(frame: "resolution", from: 13500, to: 19000) { x ->
|
||||
resolution1.value(x, 14000, params)
|
||||
}
|
||||
|
||||
plot.plot(frame: "resolution", from: 13500, to: 19000) { x ->
|
||||
resolution2.value(x, 14000, params)
|
||||
}
|
||||
|
||||
def spectrum1 = new SterileNeutrinoSpectrum(ctx, meta1)
|
||||
def spectrum2 = new SterileNeutrinoSpectrum(ctx, meta2)
|
||||
|
||||
def x = []
|
||||
def y1 = []
|
||||
def y2 = []
|
||||
(13500..19000).step(100).each {
|
||||
x << it
|
||||
y1 << spectrum1.value(it, params)
|
||||
y2 << spectrum2.value(it, params)
|
||||
}
|
||||
|
||||
plot.plot(x, y1, "spectrum1", "spectrum")
|
||||
plot.plot(x, y2, "spectrum2", "spectrum")
|
||||
|
||||
}
|
||||
|
||||
@ -26,7 +26,9 @@ import hep.dataforge.values.Values
|
||||
*
|
||||
* @author Darksnake
|
||||
*/
|
||||
open class NBkgSpectrum(private val source: ParametricFunction) : AbstractParametricFunction(*combineNamesWithEquals(source.namesAsArray(), *list)) {
|
||||
open class NBkgSpectrum(
|
||||
private val source: ParametricFunction,
|
||||
) : AbstractParametricFunction(*combineNamesWithEquals(source.namesAsArray(), *list)) {
|
||||
|
||||
var counter = MultiCounter(this.javaClass.name)
|
||||
|
||||
@ -49,7 +51,7 @@ open class NBkgSpectrum(private val source: ParametricFunction) : AbstractParame
|
||||
|
||||
override fun providesDeriv(name: String): Boolean {
|
||||
return when (name) {
|
||||
"N","bkg" -> true
|
||||
"N", "bkg" -> true
|
||||
else -> this.source.providesDeriv(name)
|
||||
}
|
||||
}
|
||||
|
||||
@ -13,40 +13,16 @@ import hep.dataforge.values.Values
|
||||
import inr.numass.models.ResolutionFunction
|
||||
import inr.numass.utils.ExpressionUtils
|
||||
import org.apache.commons.math3.analysis.BivariateFunction
|
||||
import java.lang.Math.sqrt
|
||||
import java.util.*
|
||||
import kotlin.math.sqrt
|
||||
|
||||
/**
|
||||
* @author [Alexander Nozik](mailto:altavir@gmail.com)
|
||||
*/
|
||||
class NumassResolution(context: Context, meta: Meta) : AbstractParametricBiFunction(list) {
|
||||
|
||||
private val resA: Double = meta.getDouble("A", 8.3e-5)
|
||||
private val resB = meta.getDouble("B", 0.0)
|
||||
private val tailFunction: BivariateFunction = when {
|
||||
meta.hasValue("tail") -> {
|
||||
val tailFunctionStr = meta.getString("tail")
|
||||
if (tailFunctionStr.startsWith("function::")) {
|
||||
FunctionLibrary.buildFrom(context).buildBivariateFunction(tailFunctionStr.substring(10))
|
||||
} else {
|
||||
BivariateFunction { E, U ->
|
||||
val binding = HashMap<String, Any>()
|
||||
binding["E"] = E
|
||||
binding["U"] = U
|
||||
binding["D"] = E - U
|
||||
ExpressionUtils.function(tailFunctionStr, binding)
|
||||
}
|
||||
}
|
||||
}
|
||||
meta.hasValue("tailAlpha") -> {
|
||||
//add polynomial function here
|
||||
val alpha = meta.getDouble("tailAlpha")
|
||||
val beta = meta.getDouble("tailBeta", 0.0)
|
||||
BivariateFunction { E: Double, U: Double -> 1 - (E - U) * (alpha + E / 1000.0 * beta) / 1000.0 }
|
||||
|
||||
}
|
||||
else -> ResolutionFunction.getConstantTail()
|
||||
}
|
||||
class NumassResolution(
|
||||
val resA: Double = 8.3e-5,
|
||||
val resB: Double = 0.0,
|
||||
val tailFunction: BivariateFunction = ResolutionFunction.getConstantTail(),
|
||||
) : AbstractParametricBiFunction(list) {
|
||||
|
||||
override fun derivValue(parName: String, x: Double, y: Double, set: Values): Double {
|
||||
return 0.0
|
||||
@ -82,6 +58,35 @@ class NumassResolution(context: Context, meta: Meta) : AbstractParametricBiFunct
|
||||
companion object {
|
||||
|
||||
private val list = arrayOf<String>() //leaving
|
||||
|
||||
internal fun fromMeta(context: Context, meta: Meta) = NumassResolution(
|
||||
meta.getDouble("A", 8.3e-5),
|
||||
meta.getDouble("B", 0.0),
|
||||
when {
|
||||
meta.hasValue("tail") -> {
|
||||
val tailFunctionStr = meta.getString("tail")
|
||||
if (tailFunctionStr.startsWith("function::")) {
|
||||
FunctionLibrary.buildFrom(context).buildBivariateFunction(tailFunctionStr.substring(10))
|
||||
} else {
|
||||
BivariateFunction { E, U ->
|
||||
val binding = HashMap<String, Any>()
|
||||
binding["E"] = E
|
||||
binding["U"] = U
|
||||
binding["D"] = E - U
|
||||
ExpressionUtils.function(tailFunctionStr, binding)
|
||||
}
|
||||
}
|
||||
}
|
||||
meta.hasValue("tailAlpha") -> {
|
||||
//add polynomial function here
|
||||
val alpha = meta.getDouble("tailAlpha")
|
||||
val beta = meta.getDouble("tailBeta", 0.0)
|
||||
BivariateFunction { E: Double, U: Double -> 1 - (E - U) * (alpha + E / 1000.0 * beta) / 1000.0 }
|
||||
|
||||
}
|
||||
else -> ResolutionFunction.getConstantTail()
|
||||
}
|
||||
)
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
@ -48,7 +48,7 @@ class SterileNeutrinoSpectrum @JvmOverloads constructor(
|
||||
configuration: Meta,
|
||||
val source: ParametricBiFunction = NumassBeta(),
|
||||
val transmission: ParametricBiFunction = NumassTransmission(context, configuration.getMetaOrEmpty("transmission")),
|
||||
val resolution: ParametricBiFunction = NumassResolution(context, configuration.getMeta("resolution", Meta.empty()))
|
||||
val resolution: ParametricBiFunction = NumassResolution.fromMeta(context, configuration.getMeta("resolution", Meta.empty()))
|
||||
) : AbstractParametricFunction(*list) {
|
||||
|
||||
|
||||
|
||||
@ -0,0 +1,155 @@
|
||||
/*
|
||||
* Copyright 2018 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.models
|
||||
|
||||
import hep.dataforge.buildContext
|
||||
import hep.dataforge.configure
|
||||
import hep.dataforge.fx.FXPlugin
|
||||
import hep.dataforge.fx.output.FXOutputManager
|
||||
import hep.dataforge.meta.Meta
|
||||
import hep.dataforge.plots.Plot
|
||||
import hep.dataforge.plots.data.DataPlot
|
||||
import hep.dataforge.plots.jfreechart.JFreeChartPlugin
|
||||
import hep.dataforge.plots.output.plotFrame
|
||||
import hep.dataforge.plots.plotFunction
|
||||
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.step
|
||||
import hep.dataforge.tables.Adapters.X_AXIS
|
||||
import hep.dataforge.tables.Table
|
||||
import hep.dataforge.values.ValueMap
|
||||
import inr.numass.NumassPlugin
|
||||
import inr.numass.data.SpectrumAdapter
|
||||
import inr.numass.data.SpectrumGenerator
|
||||
import inr.numass.models.NBkgSpectrum
|
||||
import inr.numass.models.sterile.NumassResolution
|
||||
import inr.numass.models.sterile.SterileNeutrinoSpectrum
|
||||
import org.apache.commons.math3.analysis.interpolation.SplineInterpolator
|
||||
import java.io.PrintWriter
|
||||
|
||||
private fun getCustomResolution(): NumassResolution {
|
||||
val correctionDataString = """
|
||||
10 0.376892
|
||||
12 0.1615868
|
||||
13 0.1009648
|
||||
14 0.0677539
|
||||
15 0.0487972
|
||||
16 0.037275
|
||||
17 0.02958922
|
||||
18 0.02511696
|
||||
19 0.02247986
|
||||
""".trimIndent()
|
||||
|
||||
val correctionData = correctionDataString.lines().map { line ->
|
||||
val (u, cor) = line.split("\t")
|
||||
u.toDouble() * 1000 to cor.toDouble()
|
||||
}
|
||||
|
||||
val correctionInterpolation = SplineInterpolator().interpolate(
|
||||
correctionData.map { it.first }.toDoubleArray(),
|
||||
correctionData.map { it.second }.toDoubleArray()
|
||||
)
|
||||
return NumassResolution(tailFunction = { e, _ -> 1.0 - correctionInterpolation.value(e) })
|
||||
}
|
||||
|
||||
private fun XYModel.plot(name: String, params: ParamSet): Plot {
|
||||
val x = (14000.0..19000.0).step(100.0).toList()
|
||||
val y = x.map { spectrum.value(it, params) }
|
||||
return DataPlot.plot(name, x.toDoubleArray(), y.toDoubleArray())
|
||||
}
|
||||
|
||||
fun main() {
|
||||
|
||||
val context = buildContext(
|
||||
"NUMASS",
|
||||
NumassPlugin::class.java,
|
||||
FXPlugin::class.java,
|
||||
JFreeChartPlugin::class.java
|
||||
) {
|
||||
output = FXOutputManager()
|
||||
}
|
||||
|
||||
val params = ParamSet().apply {
|
||||
setPar("N", 8e5, 6.0, 0.0, Double.POSITIVE_INFINITY)
|
||||
setPar("bkg", 2.0, 0.03)
|
||||
setPar("E0", 18575.0, 1.0)
|
||||
setPar("mnu2", 0.0, 1.0)
|
||||
setParValue("msterile2", (1000 * 1000).toDouble())
|
||||
setPar("U2", 0.0, 1e-3)
|
||||
setPar("X", 0.0, 0.01)
|
||||
setPar("trap", 1.0, 0.01)
|
||||
}
|
||||
|
||||
val customResolution = getCustomResolution()
|
||||
|
||||
context.plotFrame("fit", stage = "resolution") {
|
||||
plots.configure {
|
||||
"showLine" to true
|
||||
"showSymbol" to false
|
||||
"showErrors" to false
|
||||
"thickness" to 2.0
|
||||
}
|
||||
plots.setType<DataPlot>()
|
||||
|
||||
plotFunction("custom", 14000.0, 18000.0, 1000) { x ->
|
||||
customResolution.value(x, 14100.0, params)
|
||||
}
|
||||
|
||||
val basicResolution = NumassResolution()
|
||||
plotFunction("basic", 14000.0, 18000.0, 1000) { x ->
|
||||
basicResolution.value(x, 14100.0, params)
|
||||
}
|
||||
}
|
||||
|
||||
val dataSpectrum = NBkgSpectrum(SterileNeutrinoSpectrum(context, Meta.empty(), resolution = customResolution))
|
||||
|
||||
val t = 30 * 50 // time in seconds per point
|
||||
|
||||
val adapter = SpectrumAdapter(Meta.empty())
|
||||
val fm = context.getOrLoad(FitManager::class.java)
|
||||
val x = (14000.0..18500.0).step(100.0).toList()
|
||||
val dataModel = XYModel(Meta.empty(), adapter, dataSpectrum)
|
||||
|
||||
val generator = SpectrumGenerator(dataModel, params, 12316)
|
||||
|
||||
val configuration = x.map { ValueMap.ofPairs(X_AXIS to it, "time" to t) }
|
||||
val data: Table = generator.generateData(configuration)
|
||||
|
||||
val modelSpectrum = NBkgSpectrum(SterileNeutrinoSpectrum(context, Meta.empty(), resolution = NumassResolution()))
|
||||
val fitModel = XYModel(Meta.empty(), adapter, modelSpectrum)
|
||||
|
||||
context.plotFrame("fit", stage = "plots") {
|
||||
plots.configure {
|
||||
"showLine" to true
|
||||
"showSymbol" to false
|
||||
"showErrors" to false
|
||||
"thickness" to 4.0
|
||||
}
|
||||
plots.setType<DataPlot>()
|
||||
+dataModel.plot("Data", params)
|
||||
+fitModel.plot("Fit-start", params)
|
||||
}
|
||||
|
||||
val state = FitState(data, fitModel, params)
|
||||
val res = fm.runStage(state, "QOW", FitStage.TASK_RUN, "N", "E0", "bkg", "trap")
|
||||
res.printState(PrintWriter(System.out))
|
||||
val resU2 = fm.runStage(res.optState().get(), "QOW", FitStage.TASK_RUN, "N", "E0", "bkg", "trap", "U2")
|
||||
resU2.printState(PrintWriter(System.out))
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user