NPM/numass-framework
4
0
Fork 1
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Fixing fit

Browse Source
This commit is contained in:
Alexander Nozik 2018-07-05 11:14:21 +03:00
parent 1d524f46c8
commit 0e7600359d
5 changed files with 39 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
numass-main/src/main/groovy/inr/numass/scripts/Simulate.groovy
View File

@ -16,9 +16,10 @@
package inr.numass.scripts package inr.numass.scripts
import hep.dataforge.context.Global import hep.dataforge.context.Global
import hep.dataforge.io.ColumnedDataWriter
import hep.dataforge.meta.Meta import hep.dataforge.meta.Meta
import hep.dataforge.stat.fit.FitManager 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.fit.ParamSet
import hep.dataforge.stat.models.XYModel import hep.dataforge.stat.models.XYModel
import inr.numass.NumassPlugin import inr.numass.NumassPlugin
@ -37,8 +38,6 @@ import static java.util.Locale.setDefault
setDefault(Locale.US); setDefault(Locale.US);
new NumassPlugin().startGlobal() new NumassPlugin().startGlobal()
FitManager fm = Global.INSTANCE.get(FitManager)
SterileNeutrinoSpectrum sp = new SterileNeutrinoSpectrum(Global.INSTANCE, Meta.empty()); SterileNeutrinoSpectrum sp = new SterileNeutrinoSpectrum(Global.INSTANCE, Meta.empty());
//beta.setCaching(false); //beta.setCaching(false);
@ -83,14 +82,20 @@ def data = generator.generateData(DataModelUtils.getUniformSpectrumConfiguration
// allPars.setParValue("X", 0.4); // allPars.setParValue("X", 0.4);
ColumnedDataWriter.writeTable(System.out, data, "--- DATA ---"); //ColumnedDataWriter.writeTable(System.out, data, "--- DATA ---");
//FitState state = new FitState(data, model, allPars); //FitState state = new FitState(data, model, allPars);
////new PlotFitResultAction().eval(state); ////new PlotFitResultAction().eval(state);
// //
// //
//def res = fm.runStage(state, "QOW", FitStage.TASK_RUN, "N", "bkg", "E0", "U2"); //def res = fm.runStage(state, "QOW", FitStage.TASK_RUN, "N", "bkg", "E0", "U2");
//
// FitResult res = new FitHelper(Global.INSTANCE).fit(data)
// .model(model)
//res.print(out()); .params(allPars)
.stage("MINUIT", FitStage.TASK_RUN, "N")
.run();
//
//
//
res.printState(new PrintWriter(System.out));

2
numass-main/src/main/kotlin/inr/numass/models/sterile/NumassBeta.kt
View File

@ -218,7 +218,7 @@ class NumassBeta : AbstractParametricBiFunction(list) {
companion object { companion object {
private val K = 1E-23 private const val K = 1E-23
private val list = arrayOf("E0", "mnu2", "msterile2", "U2") private val list = arrayOf("E0", "mnu2", "msterile2", "U2")
} }

26
numass-main/src/main/kotlin/inr/numass/models/sterile/NumassResolution.kt
View File

@ -31,9 +31,9 @@ class NumassResolution(context: Context, meta: Meta) : AbstractParametricBiFunct
} else { } else {
BivariateFunction { E, U -> BivariateFunction { E, U ->
val binding = HashMap<String, Any>() val binding = HashMap<String, Any>()
binding.put("E", E) binding["E"] = E
binding.put("U", U) binding["U"] = U
binding.put("D", E - U) binding["D"] = E - U
ExpressionUtils.function(tailFunctionStr, binding) ExpressionUtils.function(tailFunctionStr, binding)
} }
} }
@ -54,12 +54,10 @@ class NumassResolution(context: Context, meta: Meta) : AbstractParametricBiFunct
private fun getValueFast(E: Double, U: Double): Double { private fun getValueFast(E: Double, U: Double): Double {
val delta = resA * E val delta = resA * E
return if (E - U < 0) { return when {
0.0 E - U < 0 -> 0.0
} else if (E - U > delta) { E - U > delta -> tailFunction.value(E, U)
tailFunction.value(E, U) else -> (E - U) / delta
} else {
(E - U) / delta
} }
} }
@ -74,12 +72,10 @@ class NumassResolution(context: Context, meta: Meta) : AbstractParametricBiFunct
} }
assert(resB > 0) assert(resB > 0)
val delta = resA * E val delta = resA * E
return if (E - U < 0) { return when {
0.0 E - U < 0 -> 0.0
} else if (E - U > delta) { E - U > delta -> tailFunction.value(E, U)
tailFunction.value(E, U) else -> (1 - sqrt(1 - (E - U) / E * resB)) / (1 - sqrt(1 - resA * resB))
} else {
(1 - sqrt(1 - (E - U) / E * resB)) / (1 - sqrt(1 - resA * resB))
} }
} }

8
numass-main/src/main/kotlin/inr/numass/models/sterile/NumassTransmission.kt
View File

@ -34,8 +34,8 @@ class NumassTransmission(context: Context, meta: Meta) : AbstractParametricBiFun
} else { } else {
BivariateFunction { Ei: Double, Ef: Double -> BivariateFunction { Ei: Double, Ef: Double ->
val binding = HashMap<String, Any>() val binding = HashMap<String, Any>()
binding.put("Ei", Ei) binding["Ei"] = Ei
binding.put("Ef", Ef) binding["Ef"] = Ef
ExpressionUtils.function(trapFuncStr, binding) ExpressionUtils.function(trapFuncStr, binding)
} }
} }
@ -45,7 +45,7 @@ class NumassTransmission(context: Context, meta: Meta) : AbstractParametricBiFun
} }
} }
fun getX(eIn: Double, set: Values): Double { private fun getX(eIn: Double, set: Values): Double {
return if (adjustX) { return if (adjustX) {
//From our article //From our article
set.getDouble("X") * Math.log(eIn / ION_POTENTIAL) * eIn * ION_POTENTIAL / 1.9580741410115568e6 set.getDouble("X") * Math.log(eIn / ION_POTENTIAL) * eIn * ION_POTENTIAL / 1.9580741410115568e6
@ -95,7 +95,7 @@ class NumassTransmission(context: Context, meta: Meta) : AbstractParametricBiFun
companion object { companion object {
private val list = arrayOf("trap", "X") private val list = arrayOf("trap", "X")
private val ION_POTENTIAL = 15.4//eV private const val ION_POTENTIAL = 15.4//eV
} }
} }

16
numass-main/src/main/kotlin/inr/numass/models/sterile/SterileNeutrinoSpectrum.kt
View File

@ -96,12 +96,6 @@ class SterileNeutrinoSpectrum @JvmOverloads constructor(
return 0.0 return 0.0
} }
val fsSource: (Double) -> Double = fss?.let { fss ->
{ eIn: Double -> (0 until fss.size()).sumByDouble { fss.getP(it) * sourceFunction.value(fss.getE(it), eIn, set) } }
} ?: { eIn: Double -> sourceFunction.value(0.0, eIn, set) }
val integrator: UnivariateIntegrator<*> = if (fast) { val integrator: UnivariateIntegrator<*> = if (fast) {
when { when {
eMax - u < 300 -> NumassIntegrator.getFastInterator() eMax - u < 300 -> NumassIntegrator.getFastInterator()
@ -113,7 +107,15 @@ class SterileNeutrinoSpectrum @JvmOverloads constructor(
NumassIntegrator.getHighDensityIntegrator() NumassIntegrator.getHighDensityIntegrator()
} }
return integrator.integrate(u, eMax) { eIn -> fsSource(eIn) * transResFunction.value(eIn, u, set) } 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")) { private inner class TransRes : AbstractParametricBiFunction(arrayOf("X", "trap")) {