WIP porting model from legacy
This commit is contained in:
parent
c6beaf3f83
commit
f4d26197d9
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-6.8.1-bin.zip
|
||||
distributionUrl=https\://services.gradle.org/distributions/gradle-7.0-bin.zip
|
||||
zipStoreBase=GRADLE_USER_HOME
|
||||
zipStorePath=wrapper/dists
|
||||
|
||||
@ -0,0 +1,48 @@
|
||||
/*
|
||||
* 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 hep.dataforge.stat.parametric
|
||||
|
||||
import hep.dataforge.exceptions.NameNotFoundException
|
||||
|
||||
abstract class AbstractParametric : AbstractNamedSet {
|
||||
constructor(names: NameList?) : super(names) {}
|
||||
constructor(list: Array<String?>?) : super(list) {}
|
||||
constructor(set: NameSetContainer?) : super(set) {}
|
||||
|
||||
/**
|
||||
* Provide default value for parameter with name `name`. Throws
|
||||
* NameNotFound if no default found for given parameter.
|
||||
*
|
||||
* @param name
|
||||
* @return
|
||||
*/
|
||||
protected fun getDefaultParameter(name: String?): Double {
|
||||
throw NameNotFoundException(name)
|
||||
}
|
||||
|
||||
/**
|
||||
* Extract value from input vector using default value if required parameter
|
||||
* not found
|
||||
*
|
||||
* @param name
|
||||
* @param set
|
||||
* @return
|
||||
*/
|
||||
protected fun getParameter(name: String?, set: Values): Double {
|
||||
//FIXME add default value
|
||||
return set.getDouble(name)
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,22 @@
|
||||
/*
|
||||
* 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 hep.dataforge.stat.parametric
|
||||
|
||||
import hep.dataforge.exceptions.NotDefinedException
|
||||
|
||||
abstract class AbstractParametricBiFunction : AbstractParametric, ParametricBiFunction {
|
||||
constructor(names: NameList?) : super(names) {}
|
||||
constructor(list: Array<String?>?) : super(list) {}
|
||||
constructor(set: NameSetContainer?) : super(set) {}
|
||||
|
||||
fun derivValue(parName: String?, x: Double, y: Double, set: Values?): Double {
|
||||
return if (!getNames().contains(parName)) {
|
||||
0
|
||||
} else {
|
||||
throw NotDefinedException()
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,34 @@
|
||||
/*
|
||||
* 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 hep.dataforge.stat.parametric
|
||||
|
||||
import hep.dataforge.names.NameList
|
||||
|
||||
typealias NameList = List<String>
|
||||
|
||||
/**
|
||||
*
|
||||
*
|
||||
* Abstract AbstractNamedSpectrum class.
|
||||
*
|
||||
* @author Alexander Nozik
|
||||
* @version $Id: $Id
|
||||
*/
|
||||
abstract class AbstractParametricFunction : AbstractParametric, ParametricFunction {
|
||||
constructor(names: NameList?) : super(names) {}
|
||||
constructor(vararg list: String?) : super(list) {}
|
||||
constructor(set: NameSetContainer?) : super(set) {}
|
||||
}
|
||||
@ -0,0 +1,14 @@
|
||||
/*
|
||||
* 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 hep.dataforge.stat.parametric
|
||||
|
||||
import hep.dataforge.names.NameList
|
||||
|
||||
abstract class AbstractParametricValue : AbstractParametric, ParametricValue {
|
||||
constructor(names: NameList?) : super(names) {}
|
||||
constructor(list: Array<String?>?) : super(list) {}
|
||||
constructor(set: NameSetContainer?) : super(set) {}
|
||||
}
|
||||
@ -0,0 +1,108 @@
|
||||
/*
|
||||
* 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 hep.dataforge.stat.parametric
|
||||
|
||||
import hep.dataforge.stat.fit.Param
|
||||
|
||||
/**
|
||||
*
|
||||
* DerivativeCalculator class.
|
||||
*
|
||||
* @author Alexander Nozik
|
||||
* @version $Id: $Id
|
||||
*/
|
||||
object DerivativeCalculator {
|
||||
private const val numPoints = 3
|
||||
|
||||
/**
|
||||
* Calculates finite differences derivative via 3 points differentiator.
|
||||
*
|
||||
* @param function a [hep.dataforge.stat.parametric.ParametricValue] object.
|
||||
* @param point a [hep.dataforge.stat.fit.ParamSet] object.
|
||||
* @param parName a [String] object.
|
||||
* @return a double.
|
||||
*/
|
||||
fun calculateDerivative(function: ParametricValue?, point: ParamSet, parName: String?): Double {
|
||||
val projection: UnivariateFunction = ParametricUtils.getNamedProjection(function, parName, point)
|
||||
val par: Param = point.getByName(parName)
|
||||
val diff =
|
||||
FiniteDifferencesDifferentiator(numPoints, par.getErr() / 2.0, par.getLowerBound(), par.getUpperBound())
|
||||
val derivative: UnivariateDifferentiableFunction = diff.differentiate(projection)
|
||||
val x = DerivativeStructure(1, 1, 0, point.getDouble(parName))
|
||||
val y: DerivativeStructure = derivative.value(x)
|
||||
return y.getPartialDerivative(1)
|
||||
}
|
||||
|
||||
/**
|
||||
* Calculates finite differences derivative via 3 points differentiator.
|
||||
*
|
||||
* @param function a [org.apache.commons.math3.analysis.UnivariateFunction] object.
|
||||
* @param point a double.
|
||||
* @param step a double.
|
||||
* @return a double.
|
||||
*/
|
||||
fun calculateDerivative(function: UnivariateFunction?, point: Double, step: Double): Double {
|
||||
val diff = FiniteDifferencesDifferentiator(numPoints, step)
|
||||
val derivative: UnivariateDifferentiableFunction = diff.differentiate(function)
|
||||
val x = DerivativeStructure(1, 1, 0, point)
|
||||
val y: DerivativeStructure = derivative.value(x)
|
||||
return y.getPartialDerivative(1)
|
||||
}
|
||||
|
||||
/**
|
||||
*
|
||||
* providesValidDerivative.
|
||||
*
|
||||
* @param function a [hep.dataforge.stat.parametric.ParametricValue] object.
|
||||
* @param point a [hep.dataforge.stat.fit.ParamSet] object.
|
||||
* @param tolerance a double.
|
||||
* @param parName a [String] object.
|
||||
* @return a boolean.
|
||||
*/
|
||||
fun providesValidDerivative(
|
||||
function: ParametricValue,
|
||||
point: ParamSet,
|
||||
tolerance: Double,
|
||||
parName: String?
|
||||
): Boolean {
|
||||
if (!function.providesDeriv(parName)) {
|
||||
return false
|
||||
}
|
||||
val calculatedDeriv = calculateDerivative(function, point, parName)
|
||||
val providedDeriv: Double = function.derivValue(parName, point)
|
||||
return safeRelativeDifference(calculatedDeriv, providedDeriv) <= tolerance
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns safe from (no devision by zero) relative difference between two
|
||||
* input values
|
||||
*
|
||||
* @param val1
|
||||
* @param val2
|
||||
* @return
|
||||
*/
|
||||
private fun safeRelativeDifference(val1: Double, val2: Double): Double {
|
||||
if (Precision.equals(val1, val2, Precision.EPSILON)) {
|
||||
return 0
|
||||
}
|
||||
val average: Double = Math.abs(val1 + val2) / 2
|
||||
return if (average > Precision.EPSILON) {
|
||||
Math.abs(val1 - val2) / average
|
||||
} else {
|
||||
Double.POSITIVE_INFINITY
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,75 @@
|
||||
/*
|
||||
* 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 ru.inr.mass.models
|
||||
|
||||
import hep.dataforge.names.NamesUtils.combineNamesWithEquals
|
||||
import hep.dataforge.stat.parametric.AbstractParametricFunction
|
||||
import hep.dataforge.stat.parametric.ParametricFunction
|
||||
import hep.dataforge.utils.MultiCounter
|
||||
import hep.dataforge.values.ValueProvider
|
||||
import hep.dataforge.values.Values
|
||||
import space.kscience.kmath.expressions.Symbol
|
||||
|
||||
typealias Values = Map<Symbol, Double>
|
||||
|
||||
/**
|
||||
*
|
||||
* @author Darksnake
|
||||
*/
|
||||
open class NBkgSpectrum(private val source: ParametricFunction) : AbstractParametricFunction(*combineNamesWithEquals(source.namesAsArray(), *list)) {
|
||||
|
||||
override fun derivValue(parName: String, x: Double, set: Values): Double {
|
||||
return when (parName) {
|
||||
"N" -> source.value(x, set)
|
||||
"bkg" -> 1.0
|
||||
else -> getN(set) * source.derivValue(parName, x, set)
|
||||
}
|
||||
}
|
||||
|
||||
private fun getBkg(set: ValueProvider): Double {
|
||||
return set.getDouble("bkg")
|
||||
}
|
||||
|
||||
private fun getN(set: ValueProvider): Double {
|
||||
return set.getDouble("N")
|
||||
}
|
||||
|
||||
override fun providesDeriv(name: String): Boolean {
|
||||
return when (name) {
|
||||
"N","bkg" -> true
|
||||
else -> this.source.providesDeriv(name)
|
||||
}
|
||||
}
|
||||
|
||||
override fun value(x: Double, set: Values): Double {
|
||||
this.counter.increase("value")
|
||||
return getN(set) * source.value(x, set) + getBkg(set)
|
||||
}
|
||||
|
||||
override fun getDefaultParameter(name: String): Double {
|
||||
return when (name) {
|
||||
"bkg" -> 0.0
|
||||
"N" -> 1.0
|
||||
else -> super.getDefaultParameter(name)
|
||||
}
|
||||
}
|
||||
|
||||
companion object {
|
||||
|
||||
private val list = arrayOf("N", "bkg")
|
||||
}
|
||||
|
||||
}
|
||||
@ -0,0 +1,225 @@
|
||||
/*
|
||||
* 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.models.sterile
|
||||
|
||||
import hep.dataforge.exceptions.NotDefinedException
|
||||
import hep.dataforge.stat.parametric.AbstractParametricBiFunction
|
||||
import hep.dataforge.stat.parametric.AbstractParametricFunction
|
||||
import hep.dataforge.stat.parametric.ParametricFunction
|
||||
import hep.dataforge.values.Values
|
||||
|
||||
import java.lang.Math.*
|
||||
|
||||
/**
|
||||
* A bi-function for beta-spectrum calculation taking energy and final state as
|
||||
* input.
|
||||
*
|
||||
*
|
||||
* dissertation p.33
|
||||
*
|
||||
*
|
||||
* @author [Alexander Nozik](mailto:altavir@gmail.com)
|
||||
*/
|
||||
class NumassBeta : AbstractParametricBiFunction(list) {
|
||||
|
||||
/**
|
||||
* Beta spectrum derivative
|
||||
*
|
||||
* @param n parameter number
|
||||
* @param E0
|
||||
* @param mnu2
|
||||
* @param E
|
||||
* @return
|
||||
* @throws NotDefinedException
|
||||
*/
|
||||
@Throws(NotDefinedException::class)
|
||||
private fun derivRoot(n: Int, E0: Double, mnu2: Double, E: Double): Double {
|
||||
val D = E0 - E//E0-E
|
||||
if (D == 0.0) {
|
||||
return 0.0
|
||||
}
|
||||
|
||||
return if (mnu2 >= 0) {
|
||||
if (E >= E0 - sqrt(mnu2)) {
|
||||
0.0
|
||||
} else {
|
||||
val bare = sqrt(D * D - mnu2)
|
||||
when (n) {
|
||||
0 -> factor(E) * (2.0 * D * D - mnu2) / bare
|
||||
1 -> -factor(E) * 0.5 * D / bare
|
||||
else -> 0.0
|
||||
}
|
||||
}
|
||||
} else {
|
||||
val mu = sqrt(-0.66 * mnu2)
|
||||
if (E >= E0 + mu) {
|
||||
0.0
|
||||
} else {
|
||||
val root = sqrt(Math.max(D * D - mnu2, 0.0))
|
||||
val exp = exp(-1 - D / mu)
|
||||
when (n) {
|
||||
0 -> factor(E) * (D * (D + mu * exp) / root + root * (1 - exp))
|
||||
1 -> factor(E) * (-(D + mu * exp) / root * 0.5 - root * exp * (1 + D / mu) / 3.0 / mu)
|
||||
else -> 0.0
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Derivative of spectrum with sterile neutrinos
|
||||
*
|
||||
* @param name
|
||||
* @param E
|
||||
* @param E0
|
||||
* @param pars
|
||||
* @return
|
||||
* @throws NotDefinedException
|
||||
*/
|
||||
@Throws(NotDefinedException::class)
|
||||
private fun derivRootsterile(name: String, E: Double, E0: Double, pars: Values): Double {
|
||||
val mnu2 = getParameter("mnu2", pars)
|
||||
val mst2 = getParameter("msterile2", pars)
|
||||
val u2 = getParameter("U2", pars)
|
||||
|
||||
return when (name) {
|
||||
"E0" -> {
|
||||
if (u2 == 0.0) {
|
||||
derivRoot(0, E0, mnu2, E)
|
||||
} else {
|
||||
u2 * derivRoot(0, E0, mst2, E) + (1 - u2) * derivRoot(0, E0, mnu2, E)
|
||||
}
|
||||
}
|
||||
"mnu2" -> (1 - u2) * derivRoot(1, E0, mnu2, E)
|
||||
"msterile2" -> {
|
||||
if (u2 == 0.0) {
|
||||
0.0
|
||||
} else {
|
||||
u2 * derivRoot(1, E0, mst2, E)
|
||||
}
|
||||
}
|
||||
"U2" -> root(E0, mst2, E) - root(E0, mnu2, E)
|
||||
else -> 0.0
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
* The part independent of neutrino mass. Includes global normalization
|
||||
* constant, momentum and Fermi correction
|
||||
*
|
||||
* @param E
|
||||
* @return
|
||||
*/
|
||||
private fun factor(E: Double): Double {
|
||||
val me = 0.511006E6
|
||||
val eTot = E + me
|
||||
val pe = sqrt(E * (E + 2.0 * me))
|
||||
val ve = pe / eTot
|
||||
val yfactor = 2.0 * 2.0 * 1.0 / 137.039 * Math.PI
|
||||
val y = yfactor / ve
|
||||
val fn = y / abs(1.0 - exp(-y))
|
||||
val fermi = fn * (1.002037 - 0.001427 * ve)
|
||||
val res = fermi * pe * eTot
|
||||
return K * res
|
||||
}
|
||||
|
||||
override fun providesDeriv(name: String): Boolean {
|
||||
return true
|
||||
}
|
||||
|
||||
/**
|
||||
* Bare beta spectrum with Mainz negative mass correction
|
||||
*
|
||||
* @param E0
|
||||
* @param mnu2
|
||||
* @param E
|
||||
* @return
|
||||
*/
|
||||
private fun root(E0: Double, mnu2: Double, E: Double): Double {
|
||||
//bare beta-spectrum
|
||||
val delta = E0 - E
|
||||
val bare = factor(E) * delta * sqrt(Math.max(delta * delta - mnu2, 0.0))
|
||||
return when {
|
||||
mnu2 >= 0 -> Math.max(bare, 0.0)
|
||||
delta == 0.0 -> 0.0
|
||||
delta + 0.812 * sqrt(-mnu2) <= 0 -> 0.0 //sqrt(0.66)
|
||||
else -> {
|
||||
val aux = sqrt(-mnu2 * 0.66) / delta
|
||||
Math.max(bare * (1 + aux * exp(-1 - 1 / aux)), 0.0)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* beta-spectrum with sterile neutrinos
|
||||
*
|
||||
* @param E
|
||||
* @param E0
|
||||
* @param pars
|
||||
* @return
|
||||
*/
|
||||
private fun rootsterile(E: Double, E0: Double, pars: Values): Double {
|
||||
val mnu2 = getParameter("mnu2", pars)
|
||||
val mst2 = getParameter("msterile2", pars)
|
||||
val u2 = getParameter("U2", pars)
|
||||
|
||||
return if (u2 == 0.0) {
|
||||
root(E0, mnu2, E)
|
||||
} else {
|
||||
u2 * root(E0, mst2, E) + (1 - u2) * root(E0, mnu2, E)
|
||||
}
|
||||
// P(rootsterile)+ (1-P)root
|
||||
}
|
||||
|
||||
override fun getDefaultParameter(name: String): Double {
|
||||
return when (name) {
|
||||
"mnu2", "U2", "msterile2" -> 0.0
|
||||
else -> super.getDefaultParameter(name)
|
||||
}
|
||||
}
|
||||
|
||||
override fun derivValue(parName: String, fs: Double, eIn: Double, pars: Values): Double {
|
||||
val e0 = getParameter("E0", pars)
|
||||
return derivRootsterile(parName, eIn, e0 - fs, pars)
|
||||
}
|
||||
|
||||
override fun value(fs: Double, eIn: Double, pars: Values): Double {
|
||||
val e0 = getParameter("E0", pars)
|
||||
return rootsterile(eIn, e0 - fs, pars)
|
||||
}
|
||||
|
||||
/**
|
||||
* Get univariate spectrum with given final state
|
||||
*/
|
||||
fun getSpectrum(fs: Double = 0.0): ParametricFunction {
|
||||
return BetaSpectrum(fs);
|
||||
}
|
||||
|
||||
inner class BetaSpectrum(val fs: Double) : AbstractParametricFunction(*list) {
|
||||
|
||||
override fun providesDeriv(name: String): Boolean {
|
||||
return this@NumassBeta.providesDeriv(name)
|
||||
}
|
||||
|
||||
override fun derivValue(parName: String, x: Double, set: Values): Double {
|
||||
return this@NumassBeta.derivValue(parName, fs, x, set)
|
||||
}
|
||||
|
||||
override fun value(x: Double, set: Values): Double {
|
||||
return this@NumassBeta.value(fs, x, set)
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
companion object {
|
||||
|
||||
private const val K = 1E-23
|
||||
private val list = arrayOf("E0", "mnu2", "msterile2", "U2")
|
||||
}
|
||||
|
||||
}
|
||||
@ -0,0 +1,87 @@
|
||||
/*
|
||||
* 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.models.sterile
|
||||
|
||||
import hep.dataforge.context.Context
|
||||
import hep.dataforge.maths.functions.FunctionLibrary
|
||||
import hep.dataforge.meta.Meta
|
||||
import hep.dataforge.stat.parametric.AbstractParametricBiFunction
|
||||
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.*
|
||||
|
||||
/**
|
||||
* @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()
|
||||
}
|
||||
|
||||
override fun derivValue(parName: String, x: Double, y: Double, set: Values): Double {
|
||||
return 0.0
|
||||
}
|
||||
|
||||
private fun getValueFast(E: Double, U: Double): Double {
|
||||
val delta = resA * E
|
||||
return when {
|
||||
E - U < 0 -> 0.0
|
||||
E - U > delta -> tailFunction.value(E, U)
|
||||
else -> (E - U) / delta
|
||||
}
|
||||
}
|
||||
|
||||
override fun providesDeriv(name: String): Boolean {
|
||||
return true
|
||||
}
|
||||
|
||||
override fun value(E: Double, U: Double, set: Values): Double {
|
||||
assert(resA > 0)
|
||||
if (resB <= 0) {
|
||||
return this.getValueFast(E, U)
|
||||
}
|
||||
assert(resB > 0)
|
||||
val delta = resA * E
|
||||
return when {
|
||||
E - U < 0 -> 0.0
|
||||
E - U > delta -> tailFunction.value(E, U)
|
||||
else -> (1 - sqrt(1 - (E - U) / E * resB)) / (1 - sqrt(1 - resA * resB))
|
||||
}
|
||||
}
|
||||
|
||||
companion object {
|
||||
|
||||
private val list = arrayOf<String>() //leaving
|
||||
}
|
||||
|
||||
}
|
||||
@ -0,0 +1,86 @@
|
||||
/*
|
||||
* 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.models.sterile
|
||||
|
||||
import hep.dataforge.context.Context
|
||||
import hep.dataforge.maths.functions.FunctionLibrary
|
||||
import hep.dataforge.meta.Meta
|
||||
import hep.dataforge.stat.parametric.AbstractParametricBiFunction
|
||||
import hep.dataforge.values.Values
|
||||
import inr.numass.models.misc.LossCalculator
|
||||
import inr.numass.utils.ExpressionUtils
|
||||
import org.apache.commons.math3.analysis.BivariateFunction
|
||||
import org.slf4j.LoggerFactory
|
||||
import java.util.*
|
||||
|
||||
|
||||
|
||||
/**
|
||||
* @author [Alexander Nozik](mailto:altavir@gmail.com)
|
||||
*/
|
||||
class NumassTransmission(context: Context, meta: Meta) : AbstractParametricBiFunction(list) {
|
||||
private val trapFunc: BivariateFunction
|
||||
//private val lossCache = HashMap<Double, UnivariateFunction>()
|
||||
|
||||
private val adjustX: Boolean = meta.getBoolean("adjustX", false)
|
||||
|
||||
init {
|
||||
if (meta.hasValue("trapping")) {
|
||||
val trapFuncStr = meta.getString("trapping")
|
||||
trapFunc = if (trapFuncStr.startsWith("function::")) {
|
||||
FunctionLibrary.buildFrom(context).buildBivariateFunction(trapFuncStr.substring(10))
|
||||
} else {
|
||||
BivariateFunction { Ei: Double, Ef: Double ->
|
||||
val binding = HashMap<String, Any>()
|
||||
binding["Ei"] = Ei
|
||||
binding["Ef"] = Ef
|
||||
ExpressionUtils.function(trapFuncStr, binding)
|
||||
}
|
||||
}
|
||||
} else {
|
||||
LoggerFactory.getLogger(javaClass).warn("Trapping function not defined. Using default")
|
||||
trapFunc = FunctionLibrary.buildFrom(context).buildBivariateFunction("numass.trap.nominal")
|
||||
}
|
||||
}
|
||||
|
||||
override fun derivValue(parName: String, eIn: Double, eOut: Double, set: Values): Double {
|
||||
return when (parName) {
|
||||
"trap" -> trapFunc.value(eIn, eOut)
|
||||
"X" -> LossCalculator.getTotalLossDeriv(set, eIn, eOut)
|
||||
else -> super.derivValue(parName, eIn, eOut, set)
|
||||
}
|
||||
}
|
||||
|
||||
override fun providesDeriv(name: String): Boolean {
|
||||
return true
|
||||
}
|
||||
|
||||
override fun value(eIn: Double, eOut: Double, set: Values): Double {
|
||||
// loss part
|
||||
val loss = LossCalculator.getTotalLossValue(set, eIn, eOut)
|
||||
// double loss;
|
||||
//
|
||||
// if(eIn-eOut >= 300){
|
||||
// loss = 0;
|
||||
// } else {
|
||||
// UnivariateFunction lossFunction = this.lossCache.computeIfAbsent(X, theX ->
|
||||
// FunctionCaching.cacheUnivariateFunction(0, 300, 400, x -> calculator.getTotalLossValue(theX, eIn, eIn - x))
|
||||
// );
|
||||
//
|
||||
// loss = lossFunction.value(eIn - eOut);
|
||||
// }
|
||||
|
||||
//trapping part
|
||||
val trap = getParameter("trap", set) * trapFunc.value(eIn, eOut)
|
||||
return loss + trap
|
||||
}
|
||||
|
||||
companion object {
|
||||
|
||||
private val list = arrayOf("trap", "X")
|
||||
}
|
||||
|
||||
}
|
||||
@ -0,0 +1,40 @@
|
||||
/*
|
||||
* 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 hep.dataforge.stat.parametric
|
||||
|
||||
import hep.dataforge.exceptions.NotDefinedException
|
||||
|
||||
/**
|
||||
*
|
||||
* @author Alexander Nozik
|
||||
*/
|
||||
interface ParametricBiFunction : NameSetContainer {
|
||||
fun derivValue(parName: String?, x: Double, y: Double, set: Values?): Double
|
||||
fun value(x: Double, y: Double, set: Values?): Double
|
||||
fun providesDeriv(name: String?): Boolean
|
||||
fun derivative(parName: String?): ParametricBiFunction? {
|
||||
return if (providesDeriv(parName)) {
|
||||
object : ParametricBiFunction {
|
||||
override fun derivValue(parName: String?, x: Double, y: Double, set: Values?): Double {
|
||||
throw NotDefinedException()
|
||||
}
|
||||
|
||||
override fun value(x: Double, y: Double, set: Values?): Double {
|
||||
return this@ParametricBiFunction.derivValue(parName, x, y, set)
|
||||
}
|
||||
|
||||
override fun providesDeriv(name: String?): Boolean {
|
||||
return !names.contains(name)
|
||||
}
|
||||
|
||||
val names: NameList
|
||||
get() = this@ParametricBiFunction.getNames()
|
||||
}
|
||||
} else {
|
||||
throw NotDefinedException()
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,39 @@
|
||||
/*
|
||||
* 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.models.sterile
|
||||
|
||||
import hep.dataforge.names.NameList
|
||||
import hep.dataforge.stat.parametric.ParametricBiFunction
|
||||
import hep.dataforge.values.Values
|
||||
|
||||
class ParametricBiFunctionCache(val function: ParametricBiFunction): ParametricBiFunction {
|
||||
override fun derivValue(parName: String?, x: Double, y: Double, set: Values?): Double {
|
||||
TODO("not implemented") //To change body of created functions use File | Settings | File Templates.
|
||||
}
|
||||
|
||||
override fun getNames(): NameList {
|
||||
TODO("not implemented") //To change body of created functions use File | Settings | File Templates.
|
||||
}
|
||||
|
||||
override fun value(x: Double, y: Double, set: Values?): Double {
|
||||
TODO("not implemented") //To change body of created functions use File | Settings | File Templates.
|
||||
}
|
||||
|
||||
override fun providesDeriv(name: String?): Boolean {
|
||||
TODO("not implemented") //To change body of created functions use File | Settings | File Templates.
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,58 @@
|
||||
/*
|
||||
* 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 hep.dataforge.stat.parametric
|
||||
|
||||
import hep.dataforge.exceptions.NotDefinedException
|
||||
|
||||
/**
|
||||
*
|
||||
*
|
||||
* NamedSpectrum interface.
|
||||
*
|
||||
* @author Alexander Nozik
|
||||
* @version $Id: $Id
|
||||
*/
|
||||
interface ParametricFunction : NameSetContainer {
|
||||
fun derivValue(parName: String?, x: Double, set: Values?): Double
|
||||
fun value(x: Double, set: Values?): Double
|
||||
fun providesDeriv(name: String?): Boolean
|
||||
fun derivative(parName: String?): ParametricFunction? {
|
||||
return if (providesDeriv(parName)) {
|
||||
object : ParametricFunction {
|
||||
override fun derivValue(parName: String?, x: Double, set: Values?): Double {
|
||||
return if (names.contains(parName)) {
|
||||
throw NotDefinedException()
|
||||
} else {
|
||||
0
|
||||
}
|
||||
}
|
||||
|
||||
override fun value(x: Double, set: Values?): Double {
|
||||
return this@ParametricFunction.derivValue(parName, x, set)
|
||||
}
|
||||
|
||||
override fun providesDeriv(name: String?): Boolean {
|
||||
return !names.contains(name)
|
||||
}
|
||||
|
||||
val names: NameList
|
||||
get() = this@ParametricFunction.getNames()
|
||||
}
|
||||
} else {
|
||||
throw NotDefinedException()
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,130 @@
|
||||
/*
|
||||
* 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 hep.dataforge.stat.parametric
|
||||
|
||||
import hep.dataforge.exceptions.NotDefinedException
|
||||
|
||||
/**
|
||||
* Универсальная обертка, которая объединяет именованную и обычную функцию.
|
||||
*
|
||||
* @author Alexander Nozik
|
||||
* @version $Id: $Id
|
||||
*/
|
||||
class ParametricMultiFunctionWrapper : hep.dataforge.stat.parametric.ParametricValue, MultiFunction {
|
||||
private val multiFunc: MultiFunction?
|
||||
private val nFunc: hep.dataforge.stat.parametric.ParametricValue?
|
||||
private val names: NameList
|
||||
|
||||
constructor(names: NameList, multiFunc: MultiFunction?) {
|
||||
this.names = names
|
||||
nFunc = null
|
||||
this.multiFunc = multiFunc
|
||||
}
|
||||
|
||||
constructor(nFunc: hep.dataforge.stat.parametric.ParametricValue) {
|
||||
names = nFunc.getNames()
|
||||
this.nFunc = nFunc
|
||||
multiFunc = null
|
||||
}
|
||||
|
||||
/**
|
||||
* {@inheritDoc}
|
||||
*/
|
||||
override fun derivValue(parName: String?, pars: Values): Double {
|
||||
return if (nFunc != null) {
|
||||
nFunc.derivValue(parName, pars)
|
||||
} else {
|
||||
require(pars.getNames().contains(names.asArray())) { "Wrong parameter set." }
|
||||
require(names.contains(parName)) { "Wrong derivative parameter name." }
|
||||
multiFunc.derivValue(getNumberByName(parName),
|
||||
MathUtils.getDoubleArray(pars, getNames().asArray()))
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* {@inheritDoc}
|
||||
*/
|
||||
@Throws(NotDefinedException::class)
|
||||
fun derivValue(n: Int, vector: DoubleArray?): Double {
|
||||
return if (multiFunc != null) {
|
||||
multiFunc.derivValue(n, vector)
|
||||
} else {
|
||||
val set = NamedVector(names.asArray(), vector)
|
||||
nFunc!!.derivValue(names.asArray().get(n), set)
|
||||
}
|
||||
}
|
||||
|
||||
val dimension: Int
|
||||
get() = getNames().size()
|
||||
|
||||
/**
|
||||
* {@inheritDoc}
|
||||
*/
|
||||
fun getNames(): NameList {
|
||||
return names
|
||||
}
|
||||
|
||||
private fun getNumberByName(name: String?): Int {
|
||||
return getNames().asList().indexOf(name)
|
||||
}
|
||||
|
||||
/**
|
||||
* {@inheritDoc}
|
||||
*
|
||||
*
|
||||
* Подразумевается, что аналитически заданы все(!) производные
|
||||
*/
|
||||
fun providesDeriv(n: Int): Boolean {
|
||||
return if (nFunc != null && nFunc.providesDeriv(getNames().asArray().get(n))) {
|
||||
true
|
||||
} else multiFunc != null && multiFunc.providesDeriv(n)
|
||||
}
|
||||
|
||||
/**
|
||||
* {@inheritDoc}
|
||||
*/
|
||||
override fun providesDeriv(name: String?): Boolean {
|
||||
return if (nFunc != null) {
|
||||
nFunc.providesDeriv(name)
|
||||
} else {
|
||||
multiFunc.providesDeriv(getNumberByName(name))
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* {@inheritDoc}
|
||||
*/
|
||||
override fun value(pars: Values): Double {
|
||||
return if (nFunc != null) {
|
||||
nFunc.value(pars)
|
||||
} else {
|
||||
require(pars.getNames().contains(names.asArray())) { "Wrong parameter set." }
|
||||
this.value(MathUtils.getDoubleArray(pars, getNames().asArray()))
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* {@inheritDoc}
|
||||
*/
|
||||
override fun value(vector: DoubleArray?): Double {
|
||||
return if (multiFunc != null) {
|
||||
multiFunc.value(vector)
|
||||
} else {
|
||||
val set = NamedVector(names.asArray(), vector)
|
||||
nFunc!!.value(set)
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,151 @@
|
||||
/*
|
||||
* 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 hep.dataforge.stat.parametric
|
||||
|
||||
import hep.dataforge.exceptions.NameNotFoundException
|
||||
|
||||
/**
|
||||
*
|
||||
* ParametricUtils class.
|
||||
*
|
||||
* @author Alexander Nozik
|
||||
* @version $Id: $Id
|
||||
*/
|
||||
object ParametricUtils {
|
||||
/**
|
||||
* Создает одномерное сечение многомерной именованной функции
|
||||
*
|
||||
* @param nFunc - исходная именованная функция
|
||||
* @param parName - имя параметра, по которому делается сечение
|
||||
* @param pars - Точка, вкоторой вычеслено сечение
|
||||
* @return a [UniFunction] object.
|
||||
*/
|
||||
fun getNamedProjection(nFunc: ParametricValue, parName: String?, pars: Values?): UniFunction {
|
||||
return object : UniFunction() {
|
||||
var curPars: NamedVector = NamedVector(pars)
|
||||
@Throws(NotDefinedException::class)
|
||||
fun derivValue(x: Double): Double {
|
||||
curPars.setValue(parName, x)
|
||||
return nFunc.derivValue(parName, curPars)
|
||||
}
|
||||
|
||||
fun providesDeriv(): Boolean {
|
||||
return nFunc.providesDeriv(parName)
|
||||
}
|
||||
|
||||
fun value(x: Double): Double {
|
||||
curPars.setValue(parName, x)
|
||||
return nFunc.value(curPars)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
*
|
||||
* getNamedProjectionDerivative.
|
||||
*
|
||||
* @param nFunc a [hep.dataforge.stat.parametric.ParametricValue] object.
|
||||
* @param parName a [String] object.
|
||||
* @param derivativeName a [String] object.
|
||||
* @param pars
|
||||
* @return a [org.apache.commons.math3.analysis.UnivariateFunction] object.
|
||||
*/
|
||||
fun getNamedProjectionDerivative(
|
||||
nFunc: ParametricValue,
|
||||
parName: String?, derivativeName: String?, pars: Values?
|
||||
): UnivariateFunction {
|
||||
return object : UnivariateFunction() {
|
||||
var curPars: NamedVector = NamedVector(pars)
|
||||
fun value(x: Double): Double {
|
||||
curPars.setValue(parName, x)
|
||||
return nFunc.derivValue(derivativeName, curPars)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Функция, которая запоминает исходную точку, и при нехватке параметров
|
||||
* берет значения оттуда.
|
||||
*
|
||||
* @param func a [hep.dataforge.stat.parametric.ParametricValue] object.
|
||||
* @param initPars
|
||||
* @param freePars - Описывает, каким параметрам можно будет изменяться.
|
||||
* Если null, то разрешено изменение всех параметров.
|
||||
* @return a [hep.dataforge.stat.parametric.ParametricValue] object.
|
||||
*/
|
||||
fun getNamedSubFunction(func: ParametricValue, initPars: Values, vararg freePars: String?): ParametricValue {
|
||||
require(initPars.getNames().contains(func.namesAsArray())) { "InitPars does not cover all of func parameters." }
|
||||
val names: NameList
|
||||
if (freePars.size > 0) {
|
||||
names = NameList(freePars)
|
||||
} else {
|
||||
names = initPars.getNames()
|
||||
}
|
||||
return object : AbstractParametricValue(names) {
|
||||
private val allPars: NamedVector = NamedVector(initPars)
|
||||
private val f: ParametricValue = func
|
||||
override fun derivValue(derivParName: String?, pars: Values): Double {
|
||||
if (!pars.getNames().contains(this.namesAsArray())) {
|
||||
throw NameNotFoundException()
|
||||
}
|
||||
for (name in this.getNames()) {
|
||||
allPars.setValue(name, pars.getDouble(name))
|
||||
}
|
||||
return f.derivValue(derivParName, allPars)
|
||||
}
|
||||
|
||||
override fun providesDeriv(name: String?): Boolean {
|
||||
return f.providesDeriv(name) && this.getNames().contains(name)
|
||||
}
|
||||
|
||||
override fun value(pars: Values): Double {
|
||||
if (!pars.getNames().contains(this.namesAsArray())) {
|
||||
throw NameNotFoundException()
|
||||
}
|
||||
for (name in this.getNames()) {
|
||||
allPars.setValue(name, pars.getDouble(name))
|
||||
}
|
||||
return f.value(allPars)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
fun getSpectrumDerivativeFunction(name: String?, s: ParametricFunction, pars: Values?): UnivariateFunction {
|
||||
return UnivariateFunction { x: Double -> s.derivValue(name, x, pars) }
|
||||
}
|
||||
|
||||
fun getSpectrumFunction(s: ParametricFunction, pars: Values?): UnivariateFunction {
|
||||
return UnivariateFunction { x: Double -> s.value(x, pars) }
|
||||
}
|
||||
|
||||
fun getSpectrumPointFunction(s: ParametricFunction, x: Double): ParametricValue {
|
||||
return object : AbstractParametricValue(s) {
|
||||
@Throws(NotDefinedException::class, NamingException::class)
|
||||
override fun derivValue(derivParName: String?, pars: Values?): Double {
|
||||
return s.derivValue(derivParName, x, pars)
|
||||
}
|
||||
|
||||
override fun providesDeriv(name: String?): Boolean {
|
||||
return s.providesDeriv(name)
|
||||
}
|
||||
|
||||
@Throws(NamingException::class)
|
||||
override fun value(pars: Values?): Double {
|
||||
return s.value(x, pars)
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,52 @@
|
||||
/*
|
||||
* 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 hep.dataforge.stat.parametric
|
||||
|
||||
import hep.dataforge.exceptions.NotDefinedException
|
||||
|
||||
/**
|
||||
* A function mapping parameter set to real value
|
||||
*
|
||||
* @author Alexander Nozik
|
||||
*/
|
||||
interface ParametricValue : NameSetContainer {
|
||||
/**
|
||||
* Value
|
||||
* @param pars
|
||||
* @return
|
||||
*/
|
||||
fun value(pars: Values?): Double
|
||||
|
||||
/**
|
||||
* Partial derivative value for given parameter
|
||||
* @param derivParName
|
||||
* @param pars
|
||||
* @return
|
||||
*/
|
||||
fun derivValue(derivParName: String?, pars: Values?): Double {
|
||||
throw NotDefinedException()
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns true if this object provides explicit analytical value derivative for given parameter
|
||||
*
|
||||
* @param name a [String] object.
|
||||
* @return a boolean.
|
||||
*/
|
||||
fun providesDeriv(name: String?): Boolean {
|
||||
return false
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,169 @@
|
||||
/*
|
||||
* 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.models.sterile
|
||||
|
||||
import hep.dataforge.context.Context
|
||||
import hep.dataforge.description.NodeDef
|
||||
import hep.dataforge.description.NodeDefs
|
||||
import hep.dataforge.description.ValueDef
|
||||
import hep.dataforge.description.ValueDefs
|
||||
import hep.dataforge.exceptions.NotDefinedException
|
||||
import hep.dataforge.maths.integration.UnivariateIntegrator
|
||||
import hep.dataforge.meta.Meta
|
||||
import hep.dataforge.stat.parametric.AbstractParametricBiFunction
|
||||
import hep.dataforge.stat.parametric.AbstractParametricFunction
|
||||
import hep.dataforge.stat.parametric.ParametricBiFunction
|
||||
import hep.dataforge.values.ValueType.BOOLEAN
|
||||
import hep.dataforge.values.Values
|
||||
import inr.numass.getFSS
|
||||
import inr.numass.models.FSS
|
||||
import inr.numass.models.misc.LossCalculator
|
||||
import inr.numass.utils.NumassIntegrator
|
||||
|
||||
/**
|
||||
* Compact all-in-one model for sterile neutrino spectrum
|
||||
*
|
||||
* @author Alexander Nozik
|
||||
*/
|
||||
@NodeDefs(
|
||||
NodeDef(key = "resolution"),
|
||||
NodeDef(key = "transmission")
|
||||
)
|
||||
@ValueDefs(
|
||||
ValueDef(key = "fssFile", info = "The name for external FSS file. By default internal FSS file is used"),
|
||||
ValueDef(key = "useFSS", type = arrayOf(BOOLEAN))
|
||||
)
|
||||
|
||||
/**
|
||||
* @param source variables:Eo offset,Ein; parameters: "mnu2", "msterile2", "U2"
|
||||
* @param transmission variables:Ein,Eout; parameters: "A"
|
||||
* @param resolution variables:Eout,U; parameters: "X", "trap"
|
||||
*/
|
||||
|
||||
class SterileNeutrinoSpectrum @JvmOverloads constructor(
|
||||
context: Context,
|
||||
configuration: Meta,
|
||||
val source: ParametricBiFunction = NumassBeta(),
|
||||
val transmission: ParametricBiFunction = NumassTransmission(context, configuration.getMetaOrEmpty("transmission")),
|
||||
val resolution: ParametricBiFunction = NumassResolution(context, configuration.getMeta("resolution", Meta.empty()))
|
||||
) : AbstractParametricFunction(*list) {
|
||||
|
||||
|
||||
/**
|
||||
* auxiliary function for trans-res convolution
|
||||
*/
|
||||
private val transRes: ParametricBiFunction = TransRes()
|
||||
private val fss: FSS? = getFSS(context, configuration)
|
||||
// private boolean useMC;
|
||||
private val fast: Boolean = configuration.getBoolean("fast", true)
|
||||
|
||||
override fun derivValue(parName: String, u: Double, set: Values): Double {
|
||||
return when (parName) {
|
||||
"U2", "msterile2", "mnu2", "E0" -> integrate(u, source.derivative(parName), transRes, set)
|
||||
"X", "trap" -> integrate(u, source, transRes.derivative(parName), set)
|
||||
else -> throw NotDefinedException()
|
||||
}
|
||||
}
|
||||
|
||||
override fun value(u: Double, set: Values): Double {
|
||||
return integrate(u, source, transRes, set)
|
||||
}
|
||||
|
||||
override fun providesDeriv(name: String): Boolean {
|
||||
return source.providesDeriv(name) && transmission.providesDeriv(name) && resolution.providesDeriv(name)
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Direct Gauss-Legendre integration
|
||||
*
|
||||
* @param u
|
||||
* @param sourceFunction
|
||||
* @param transResFunction
|
||||
* @param set
|
||||
* @return
|
||||
*/
|
||||
private fun integrate(
|
||||
u: Double,
|
||||
sourceFunction: ParametricBiFunction,
|
||||
transResFunction: ParametricBiFunction,
|
||||
set: Values): Double {
|
||||
|
||||
val eMax = set.getDouble("E0") + 5.0
|
||||
|
||||
if (u >= eMax) {
|
||||
return 0.0
|
||||
}
|
||||
|
||||
val integrator: UnivariateIntegrator<*> = if (fast) {
|
||||
when {
|
||||
eMax - u < 300 -> NumassIntegrator.getFastInterator()
|
||||
eMax - u > 2000 -> NumassIntegrator.getHighDensityIntegrator()
|
||||
else -> NumassIntegrator.getDefaultIntegrator()
|
||||
}
|
||||
|
||||
} else {
|
||||
NumassIntegrator.getHighDensityIntegrator()
|
||||
}
|
||||
|
||||
return integrator.integrate(u, eMax) { eIn -> sumByFSS(eIn, sourceFunction, set) * transResFunction.value(eIn, u, set) }
|
||||
}
|
||||
|
||||
private fun sumByFSS(eIn: Double, sourceFunction: ParametricBiFunction, set: Values): Double {
|
||||
return if (fss == null) {
|
||||
sourceFunction.value(0.0, eIn, set)
|
||||
} else {
|
||||
(0 until fss.size()).sumByDouble { fss.getP(it) * sourceFunction.value(fss.getE(it), eIn, set) }
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
private inner class TransRes : AbstractParametricBiFunction(arrayOf("X", "trap")) {
|
||||
|
||||
override fun providesDeriv(name: String): Boolean {
|
||||
return true
|
||||
}
|
||||
|
||||
override fun derivValue(parName: String, eIn: Double, u: Double, set: Values): Double {
|
||||
return when (parName) {
|
||||
"X" -> throw NotDefinedException()//TODO implement p0 derivative
|
||||
"trap" -> lossRes(transmission.derivative(parName), eIn, u, set)
|
||||
else -> super.derivValue(parName, eIn, u, set)
|
||||
}
|
||||
}
|
||||
|
||||
override fun value(eIn: Double, u: Double, set: Values): Double {
|
||||
|
||||
val p0 = LossCalculator.p0(set, eIn)
|
||||
return p0 * resolution.value(eIn, u, set) + lossRes(transmission, eIn, u, set)
|
||||
}
|
||||
|
||||
private fun lossRes(transFunc: ParametricBiFunction, eIn: Double, u: Double, set: Values): Double {
|
||||
val integrand = { eOut: Double -> transFunc.value(eIn, eOut, set) * resolution.value(eOut, u, set) }
|
||||
|
||||
val border = u + 30
|
||||
val firstPart = NumassIntegrator.getFastInterator().integrate(u, Math.min(eIn, border), integrand)
|
||||
val secondPart: Double = if (eIn > border) {
|
||||
if (fast) {
|
||||
NumassIntegrator.getDefaultIntegrator().integrate(border, eIn, integrand)
|
||||
} else {
|
||||
NumassIntegrator.getHighDensityIntegrator().integrate(border, eIn, integrand)
|
||||
}
|
||||
} else {
|
||||
0.0
|
||||
}
|
||||
return firstPart + secondPart
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
companion object {
|
||||
|
||||
private val list = arrayOf("X", "trap", "E0", "mnu2", "msterile2", "U2")
|
||||
}
|
||||
|
||||
}
|
||||
@ -6,8 +6,8 @@ pluginManagement {
|
||||
gradlePluginPortal()
|
||||
}
|
||||
|
||||
val toolsVersion = "0.9.3"
|
||||
val kotlinVersion = "1.4.32"
|
||||
val toolsVersion = "0.9.5-dev"
|
||||
val kotlinVersion = "1.5.0-M2"
|
||||
|
||||
plugins {
|
||||
id("ru.mipt.npm.gradle.project") version toolsVersion
|
||||
|
||||
Loading…
Reference in New Issue
Block a user