diff --git a/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassTransmission.kt b/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassTransmission.kt index 075ce99..d6e72e5 100644 --- a/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassTransmission.kt +++ b/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassTransmission.kt @@ -20,7 +20,7 @@ import kotlin.math.* * @author [Alexander Nozik](mailto:altavir@gmail.com) */ public class NumassTransmission( - public val trapFunc: Kernel = defaultTrapping, + public val trapFunc: Kernel = trapFunction, // private val adjustX: Boolean = false, ) : DifferentiableKernel { // private val trapFunc: Kernel = if (meta.hasValue("trapping")) { @@ -80,8 +80,8 @@ public class NumassTransmission( //trapping part // - //val trap = (arguments[trap] ?: 1.0) * trapFunc(ei, ef, arguments) - val trap = 0.0 + val trap = (arguments[trap] ?: 1.0) * trapFunc(ei, ef, arguments) + //val trap = 0.0 return loss + trap } @@ -136,6 +136,7 @@ public class NumassTransmission( } private fun getCachedSpectrum(order: Int): Function1D { + //println("-------------------ORDER $order ------------------------") return when { order <= 0 -> error("Non-positive loss cache order") order == 1 -> singleScatterFunction diff --git a/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/fit.kt b/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/fit.kt index fa2b0b2..2532f9d 100644 --- a/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/fit.kt +++ b/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/fit.kt @@ -27,8 +27,8 @@ fun parse(filename: String): XYErrorColumnarData { val array = it.split("\t").map { it.toDouble() } x.add(array[0]) y.add(array[1]) - //errors.add(array[2]) - errors.add(1.0) + errors.add(array[2]) + //errors.add(array[1] / 100) } return XYErrorColumnarData.of( x.asBuffer(), @@ -41,16 +41,16 @@ fun parse(filename: String): XYErrorColumnarData { suspend fun main() { val spectrum: NBkgSpectrum = SterileNeutrinoSpectrum( // fss = FSS.default, - transmission = NumassTransmission(NumassTransmission.trapFunction), - resolution = NumassResolution(1.7e-4) + transmission = NumassTransmission(NumassTransmission.trapFunction), +// resolution = NumassResolution(1.7e-4) ).withNBkg() val args: Map = mapOf( - NBkgSpectrum.norm to 319034.0, - NBkgSpectrum.bkg to 0.029, + NBkgSpectrum.norm to 8.84e8, + NBkgSpectrum.bkg to -1.702, NumassBeta.mnu2 to 0.0, - NumassBeta.e0 to 18575.0, - NumassBeta.msterile2 to 1000.0.pow(2), + NumassBeta.e0 to 18572.0, + NumassBeta.msterile2 to 0.0.pow(2), NumassBeta.u2 to 1e-2, NumassTransmission.thickness to 0.3, NumassTransmission.trap to 1.0 @@ -64,40 +64,41 @@ suspend fun main() { // // val strategy = (12000.0..19000.0 step 100.0).asSequence().associateWith { timePerPoint } - val data = parse("/home/sabina/Numass/Fit/tr5wobkg") + val data = parse("/home/sabina/Numass/Fit/Data/03.10.2023") - //val fit: XYFit = data.fitWith( - // optimizer = QowOptimizer, - // modelExpression = spectrum, - // startingPoint = args, - // OptimizationParameters(NBkgSpectrum.norm, NBkgSpectrum.bkg), - // OptimizationIterations(20) - //) + val fit: XYFit = data.fitWith( + optimizer = QowOptimizer, + modelExpression = spectrum, + startingPoint = args, + OptimizationParameters(NumassBeta.e0, NBkgSpectrum.norm, NBkgSpectrum.bkg), + OptimizationIterations(20) + ) //println("Chi squared/dof: ${fit.chiSquaredOrNull}/${fit.dof}") Plotly.page { plot { -/* scatter { + + scatter { name = "Data" mode = ScatterMode.markers x.buffer = data.x y.buffer = data.y - }*/ + } + scatter { name = "Initial" mode = ScatterMode.lines x.buffer = 12000.0..18600.0 step 50.0 - y.numbers = x.doubles.map { spectrum(it, args) } - File("/home/sabina/Numass/Fit/output").printWriter().use { out -> + y.numbers = x.doubles.map { spectrum(it, args + fit.resultPoint) } + File("/home/sabina/Numass/Fit/trap1").printWriter().use { out -> y.numbers.forEach { out.println(it) } } } - } -/* scatter { + scatter { name = "Fit" mode = ScatterMode.lines x.buffer = 12000.0..18600.0 step 10.0 @@ -125,6 +126,18 @@ suspend fun main() { } name = "Res histo" } - }*/ + } + plot { + scatter { + name = "Residuals_lines" + mode = ScatterMode.lines + x.buffer = data.x + y.numbers = data.indices.map { + val value = spectrum(data.x[it], args + fit.resultPoint) + val dif = data.y[it] - value + dif / data.yErr[it] + } + } + } }.makeFile() } \ No newline at end of file