optimization

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<Double> {
+            //println("-------------------ORDER $order ------------------------")
             return when {
                 order <= 0 -> error("Non-positive loss cache order")
                 order == 1 -> singleScatterFunction

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

@@ -27,8 +27,8 @@ fun parse(filename: String): XYErrorColumnarData<Double, Double, Double> {
         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<Double, Double, Double> {
 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<Symbol, Double> = 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()
 }