model correction

numass-data-proto/src/gen/kotlin/ru/inr/mass/data/proto/Point.kt

@@ -2,10 +2,29 @@
 // Source: ru.inr.mass.data.proto.Point in numass-proto.proto
 package ru.inr.mass.`data`.proto
 
-import com.squareup.wire.*
+import com.squareup.wire.FieldEncoding
+import com.squareup.wire.Message
+import com.squareup.wire.ProtoAdapter
+import com.squareup.wire.ProtoReader
+import com.squareup.wire.ProtoWriter
+import com.squareup.wire.ReverseProtoWriter
+import com.squareup.wire.Syntax
 import com.squareup.wire.Syntax.PROTO_3
-import com.squareup.wire.internal.immutableCopyOf
+import com.squareup.wire.WireField
-import com.squareup.wire.internal.redactElements
+import com.squareup.wire.`internal`.immutableCopyOf
+import com.squareup.wire.`internal`.redactElements
+import kotlin.Any
+import kotlin.AssertionError
+import kotlin.Boolean
+import kotlin.Deprecated
+import kotlin.DeprecationLevel
+import kotlin.Int
+import kotlin.Long
+import kotlin.Nothing
+import kotlin.String
+import kotlin.Unit
+import kotlin.collections.List
+import kotlin.jvm.JvmField
 import okio.ByteString
 
 public class Point(

numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassResolution.kt

@@ -15,9 +15,10 @@ import kotlin.math.sqrt
 
 @Suppress("PARAMETER_NAME_CHANGED_ON_OVERRIDE")
 public class NumassResolution(
-    public val resA: Double = 8.3e-5,
+    public val resA: Double = 1.7e-4,
     public val resB: Double = 0.0,
-    public val tailFunction: (Double, Double) -> Double = { _, _ -> 1.0 },
+    // Recent formula for adiabacity
+    public val tailFunction: (Double, Double) -> Double = { e, u -> -1.75559e-13 * (e - u) * (e - u) * (e - u) - 5.97479e-11 * (e - u) * (e - u) + 3.26473e-7 * (e - u) + 1.0 },
 ) : DifferentiableKernel {
 
 //    private val tailFunction: Kernel = when {

numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassTransmission.kt

@@ -54,8 +54,10 @@ public class NumassTransmission(
                 }
                 sum
             }
+
             else -> null
         }
+
         else -> null
     }
 
@@ -63,6 +65,7 @@ public class NumassTransmission(
         // loss part
         val thickness = arguments[thickness] ?: 0.0
         val loss = getTotalLossValue(thickness, ei, ef)
+        //val loss = 0.0
         //        double loss;
         //
         //        if(eIn-eOut >= 300){
@@ -76,7 +79,9 @@ public class NumassTransmission(
         //        }
 
         //trapping part
-        val trap = (arguments[trap] ?: 1.0) * trapFunc(ei, ef, arguments)
+        //
+        //val trap = (arguments[trap] ?: 1.0) * trapFunc(ei, ef, arguments)
+        val trap = 0.0
         return loss + trap
     }
 
@@ -86,7 +91,7 @@ public class NumassTransmission(
 
         private val cache = HashMap<Int, Function1D<Double>>()
 
-        private const val ION_POTENTIAL = 15.4//eV
+        private const val ION_POTENTIAL = 13.6//eV
 
 
         private fun getX(arguments: Map<Symbol, Double>, eIn: Double, adjustX: Boolean = false): Double {
@@ -331,6 +336,7 @@ public class NumassTransmission(
                     val z = eps - pos1
                     A1 * exp(-2.0 * z * z / w1 / w1)
                 }
+
                 else -> {
                     val z = 4.0 * (eps - pos2) * (eps - pos2)
                     A2 / (1 + z / w2 / w2)
@@ -399,7 +405,7 @@ public class NumassTransmission(
 //            return getSingleScatterFunction(exPos, ionPos, exW, ionW, exIonRatio)
 //        }
 
-        public val trapFunction: (Double, Double) -> Double = { ei: Double, ef: Double ->
+        public val trapFunction: Kernel = Kernel { ei: Double, ef: Double, _ ->
             val eps = ei - ef
             if (eps > 10) {
                 1.86e-04 * exp(-eps / 25.0) + 5.5e-05

numass-workspace/src/main/kotlin/ru/inr/mass/scripts/fit.kt

@@ -3,28 +3,51 @@ package ru.inr.mass.scripts
 import ru.inr.mass.models.*
 import ru.inr.mass.workspace.buffer
 import ru.inr.mass.workspace.fitWith
-import ru.inr.mass.workspace.generate
+import space.kscience.kmath.data.XYColumnarData
+import space.kscience.kmath.data.XYErrorColumnarData
+import space.kscience.kmath.data.indices
 import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.expressions.derivative
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.asSequence
 import space.kscience.kmath.optimization.*
 import space.kscience.kmath.real.step
+import space.kscience.kmath.structures.asBuffer
 import space.kscience.plotly.*
 import space.kscience.plotly.models.ScatterMode
+import java.io.File
 import kotlin.math.pow
 
+fun parse(filename: String): XYErrorColumnarData<Double, Double, Double> {
+    val x: MutableList<Double> = mutableListOf()
+    val y: MutableList<Double> = mutableListOf()
+    val errors: MutableList<Double> = mutableListOf()
+
+    File(filename).forEachLine {
+        val array = it.split("\t").map { it.toDouble() }
+        x.add(array[0])
+        y.add(array[1])
+        //errors.add(array[2])
+        errors.add(1.0)
+    }
+    return XYErrorColumnarData.of(
+        x.asBuffer(),
+        y.asBuffer(),
+        errors.map { it }.asBuffer()
+    )
+}
 
 @OptIn(UnstableKMathAPI::class)
 suspend fun main() {
     val spectrum: NBkgSpectrum = SterileNeutrinoSpectrum(
-        fss = FSS.default,
+//        fss = FSS.default,
-//        resolution = NumassResolution(8.2e-5)
+        transmission = NumassTransmission(NumassTransmission.trapFunction),
+        resolution = NumassResolution(1.7e-4)
     ).withNBkg()
 
     val args: Map<Symbol, Double> = mapOf(
-        NBkgSpectrum.norm to 8e5,
+        NBkgSpectrum.norm to 319034.0,
-        NBkgSpectrum.bkg to 2.0,
+        NBkgSpectrum.bkg to 0.029,
         NumassBeta.mnu2 to 0.0,
         NumassBeta.e0 to 18575.0,
         NumassBeta.msterile2 to 1000.0.pow(2),
@@ -33,46 +56,75 @@ suspend fun main() {
         NumassTransmission.trap to 1.0
     )
 
-    listOf(NBkgSpectrum.norm, NBkgSpectrum.bkg, NumassBeta.e0).forEach {
+//    listOf(NBkgSpectrum.norm, NBkgSpectrum.bkg, NumassBeta.e0).forEach {
-        println("$it: ${spectrum.derivative(it).invoke(14000.0, args)}")
+//        println("$it: ${spectrum.derivative(it).invoke(14000.0, args)}")
-    }
+//    }
 
-    val timePerPoint = 30.0
+//    val timePerPoint = 30.0
+//
+//    val strategy = (12000.0..19000.0 step 100.0).asSequence().associateWith { timePerPoint }
 
-    val strategy = (12000.0..19000.0 step 100.0).asSequence().associateWith { timePerPoint }
+    val data = parse("/home/sabina/Numass/Fit/tr5wobkg")
 
-    val generatedData = spectrum.generate(strategy, args)
+    //val fit: XYFit = data.fitWith(
+    //    optimizer = QowOptimizer,
+    //    modelExpression = spectrum,
+    //    startingPoint = args,
+    //    OptimizationParameters(NBkgSpectrum.norm, NBkgSpectrum.bkg),
+    //    OptimizationIterations(20)
+    //)
 
-    val fit: XYFit = generatedData.fitWith(
+    //println("Chi squared/dof: ${fit.chiSquaredOrNull}/${fit.dof}")
-        optimizer = QowOptimizer,
-        modelExpression = spectrum,
-        startingPoint = args + mapOf(NBkgSpectrum.norm to 8.1e5),
-        OptimizationParameters(NBkgSpectrum.norm, NBkgSpectrum.bkg, NumassBeta.e0),
-        OptimizationIterations(20)
-    )
 
-    println("Chi squared/dof: ${fit.chiSquaredOrNull}/${fit.dof}")
+    Plotly.page {
+        plot {
+/*            scatter {
+                name = "Data"
+                mode = ScatterMode.markers
+                x.buffer = data.x
+                y.buffer = data.y
+            }*/
 
-    Plotly.plot {
+            scatter {
-        scatter {
+                name = "Initial"
-            name = "Generated"
+                mode = ScatterMode.lines
-            mode = ScatterMode.markers
+                x.buffer = 12000.0..18600.0 step 50.0
-            x.buffer = generatedData.x
+                y.numbers = x.doubles.map { spectrum(it, args) }
-            y.buffer = generatedData.y
+                File("/home/sabina/Numass/Fit/output").printWriter().use { out ->
+                    y.numbers.forEach {
+                        out.println(it)
+                    }
+                }
+            }
         }
-
+/*            scatter {
-        scatter {
+                name = "Fit"
-            name = "Initial"
+                mode = ScatterMode.lines
-            mode = ScatterMode.lines
+                x.buffer = 12000.0..18600.0 step 10.0
-            x.buffer = 12000.0..18600.0 step 50.0
+                y.numbers = x.doubles.map { spectrum(it, args + fit.resultPoint) }
-            y.numbers = x.doubles.map { spectrum(it, args) }
+            }
         }
-
+        plot{
-        scatter {
+            scatter {
-            name = "Fit"
+                name = "Residuals"
-            mode = ScatterMode.lines
+                mode = ScatterMode.markers
-            x.buffer = 12000.0..18600.0 step 10.0
+                x.buffer = data.x
-            y.numbers = x.doubles.map { spectrum(it, args + fit.resultPoint) }
+                y.numbers = data.indices.map{
+                    val value = spectrum(data.x[it], args + fit.resultPoint)
+                    val dif = data.y[it] - value
+                    dif/data.yErr[it]
+                }
+            }
         }
+        plot{
+            histogram {
+                    x.numbers = data.indices.map{
+                        val value = spectrum(data.x[it], args + fit.resultPoint)
+                        val dif = data.y[it] - value
+                        dif/data.yErr[it]
+                    }
+                    name = "Res histo"
+            }
+        }*/
     }.makeFile()
 }