feat(tests): add trap and rear wall contribution graphs

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

@@ -229,15 +229,6 @@ suspend fun processCustom(
                 mode = ScatterMode.lines
                 x.buffer = data.x
                 y.numbers = x.doubles.map { spectrum(it, fitParams + fit.result) }
-                File(dataPath.parent.toString(), dataPath.nameWithoutExtension + postfix + ".fit")
-                    .printWriter().use {
-                        out -> out.println("U_sp\tcounts\tresidials")
-                        data.indices.map{
-                            val value = spectrum(data.x[it], fitParams + fit.result)
-                            val dif = data.y[it] - value
-                            out.println("${data.x[it]}\t${data.y[it]}\t${dif/data.yErr[it]}")
-                    }
-                }
             }
             scatter {
                 name = "Fit"
@@ -259,6 +250,17 @@ suspend fun processCustom(
                     val dif = testData.y[it] - value
                     dif/testData.yErr[it]
                 }
+
+                File(dataPath.parent.toString(), dataPath.nameWithoutExtension + postfix + ".fit.tsv")
+                .printWriter().use {
+                    out -> out.println("U_sp\tcounts\tresidials")
+                    testData.indices.forEach{
+                        val value = spectrum(testData.x[it], fitParams + fit.result)
+                        val dif = testData.y[it] - value
+                        val residial = dif/testData.yErr[it]
+                        out.println("${testData.x[it]}\t${value}\t$residial")
+                    }
+                }
             }
         }
         plot {
@@ -296,11 +298,16 @@ suspend fun processCustom(
         table {
             tr {
                 th { +"U_sp" }
+                th { +"Counts" }
                 th { +"residials" }
             }
             testData.indices.forEach {
                 tr {
                     td { +testData.x[it].toString() }
+                    td {
+                        val value = spectrum(testData.x[it], fitParams + fit.result)
+                        +value.toString()
+                    }
                     td {
                         val value = spectrum(testData.x[it], fitParams + fit.result)
                         val dif = testData.y[it] - value

numass-workspace/src/main/kotlin/ru/inr/mass/scripts/trap-compare-spectrums.kt

@@ -0,0 +1,85 @@
+package ru.inr.mass.scripts
+
+import ru.inr.mass.models.*
+import ru.inr.mass.models.NBkgSpectrum.Companion.bkg
+import ru.inr.mass.models.NBkgSpectrum.Companion.norm
+import ru.inr.mass.models.NumassBeta.e0
+import ru.inr.mass.models.NumassBeta.mnu2
+import ru.inr.mass.models.NumassBeta.msterile2
+import ru.inr.mass.models.NumassBeta.u2
+import ru.inr.mass.models.NumassTransmission.Companion.thickness
+import ru.inr.mass.models.NumassTransmission.Companion.trap
+import ru.inr.mass.workspace.buffer
+import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.real.step
+import space.kscience.plotly.*
+import space.kscience.plotly.models.ScatterMode
+
+fun main() {
+
+    val trapInterpolator = TrapInterpolator()
+    val range = 10000.0..18600.0 step 50.0
+
+    val args: Map<Symbol, Double> = mapOf(
+        norm to 231296.31435970013,
+        bkg to 1.156996613764463,
+        mnu2 to 0.0,
+        e0 to 18553.267180007213,
+        msterile2 to 16e6,
+        u2 to -0.006326658267906115,
+        thickness to 0.1,
+        trap to 0.0,
+        rearWall to 0.04192630322035539
+    )
+
+    val argsTrap1 = args.toMutableMap();
+    argsTrap1[trap] = 1.0
+
+    val argsTrap2 = args.toMutableMap();
+    argsTrap2[trap] = 2.0
+
+    val argsTrap3 = args.toMutableMap();
+    argsTrap3[trap] = 3.0
+
+    val spectrumC: NBkgSpectrum = CustomSterileNeutrinoSpectrum(
+        fss = null,
+        transmission = NumassTransmission(
+            trapFunc = { ei, ef, _ ->
+                val delta = ei - ef
+                trapInterpolator.value(ei, delta)
+            },
+            adjustX = false,
+        ),
+        resolution = NumassResolution(1.7e-4, tailFunction = {
+                e,u -> adiabacity(e,u)
+        })
+    ).withNBkg()
+
+    Plotly.page {
+        plot {
+            scatter {
+                name = "CustomSterileNeutrinoSpectrum (trap=${args[trap]})"
+                mode = ScatterMode.lines
+                x.buffer = range
+                y.numbers = x.doubles.map { spectrumC(it, args) }
+            }
+            listOf(argsTrap1, argsTrap2, argsTrap3).forEach { argsTrap ->
+                scatter {
+                    name = "CustomSterileNeutrinoSpectrum (trap=${argsTrap[trap]})"
+                    mode = ScatterMode.lines
+                    x.buffer = range
+                    y.numbers = x.doubles.map { spectrumC(it, argsTrap) }
+                }
+            }
+            layout {
+                title = "Compare Trap"
+//                yaxis.type = AxisType.log
+            }
+        }
+    }.makeFile()
+
+    println("HV\ttrap=${args[trap]}\ttrap=${argsTrap3[trap]}")
+    range.iterator().forEach {
+        println("${it}\t${spectrumC(it, args)}\t${spectrumC(it, argsTrap3)}")
+    }
+}

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

@@ -8,6 +8,7 @@ import space.kscience.plotly.models.TraceType
 import space.kscience.dataforge.meta.Value
 import space.kscience.kmath.structures.asBuffer
 import space.kscience.plotly.models.ScatterMode
+import kotlin.math.exp
 
 fun main() {
     val maxDelta = (WALL_R - WALL_L)
@@ -35,6 +36,18 @@ fun main() {
         }
 
         plot {
+            scatter {
+                name = "Nozik Trap"
+                mode = ScatterMode.lines
+                x.buffer =  (0.0..90.0 step 1.0).toDoubleArray().plus((100.0..maxDelta step 10.0).toDoubleArray()).asBuffer()
+                y.numbers = x.doubles.map {
+                    if (it > 10) {
+                        1.86e-04 * exp(-it / 25.0) + 5.5e-05
+                    } else {
+                        0.0
+                    }
+                }
+            }
             for (HV in (WALL_L..WALL_R step 1000.0)) {
                 scatter {
                     name = "${HV / 1000.0} keV"

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

@@ -60,17 +60,24 @@ class TrapInterpolator {
     private val fineInterpolator = BilinearInterpolator(x, yFine, gridFine)
 
     fun value(ei: Double, delta: Double): Double {
+
+        var eiL = ei;
+
         if (ei < x.first() || ei > x.last())
-            error("trap: ei ($ei) not in [${x.first()}..${x.last()}]")
+            if (ei > x.last() && (ei - x.last()) < 200.0) {
+                eiL = x.last()
+            } else {
+                error("trap: ei ($ei) not in [${x.first()}..${x.last()}]")
+            }
         if (delta < yFine.first())
             error("trap: delta < ${yFine.first()}")
         if (delta > yCoarse.last())
             error("trap: delta > ${yCoarse.last()}")
 
         return if (delta < yFine.last()) {
-            fineInterpolator.value(ei, delta)
+            fineInterpolator.value(eiL, delta)
         } else {
-            coarseInterpolator.value(ei, delta)
+            coarseInterpolator.value(eiL, delta)
         }
     }
 }

numass-workspace/src/main/kotlin/ru/inr/mass/scripts/wall-compare-spectrums.kt

@@ -13,41 +13,27 @@ import ru.inr.mass.workspace.buffer
 import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.real.step
 import space.kscience.plotly.*
-import space.kscience.plotly.models.AxisType
 import space.kscience.plotly.models.ScatterMode
-import kotlin.math.pow
 
 fun main() {
 
-    val trapInterpolator = RearTrapInterpolator()
-    val range = 12000.0..17000.0 step 50.0
-
-    val spectrum: NBkgSpectrum = SterileNeutrinoSpectrum(
-        fss = null,
-        transmission = NumassTransmission(
-            trapFunc = { ei, ef, _ ->
-                val delta = ei - ef
-                trapInterpolator.value(ei, delta)
-            },
-            adjustX = false,
-        ),
-        resolution = NumassResolution(1.7e-4, tailFunction = {
-             e,u -> adiabacity(e,u)
-        })
-    ).withNBkg()
+    val trapInterpolator = TrapInterpolator()
+    val range = 10000.0..18600.0 step 50.0
 
     val args: Map<Symbol, Double> = mapOf(
-        norm to 1.0,
-        bkg to 0.0,
+        norm to 231296.31435970013,
+        bkg to 1.156996613764463,
         mnu2 to 0.0,
-        e0 to 18575.0,
-        msterile2 to 0000.0.pow(2),
-        u2 to 0e-2,
-        thickness to 0.0,
-        trap to 1.0,
-        rearWall to 0.02
+        e0 to 18553.267180007213,
+        msterile2 to 16e6,
+        u2 to -0.006326658267906115,
+        thickness to 0.1,
+        trap to 3.0,
+        rearWall to 0.0
     )
 
+    val argsRWall = args.toMutableMap();
+    argsRWall[rearWall] = 0.04
 
     val spectrumC: NBkgSpectrum = CustomSterileNeutrinoSpectrum(
         fss = null,
@@ -64,28 +50,28 @@ fun main() {
     ).withNBkg()
 
     Plotly.page {
-
         plot {
             scatter {
-                name = "SterileNeutrinoSpectrum"
+                name = "CustomSterileNeutrinoSpectrum (rearWall=${args[rearWall]})"
                 mode = ScatterMode.lines
                 x.buffer = range
-                y.numbers = x.doubles.map {
-                    val value = spectrum(it, args)
-//                    println("${it}\t${value}")
-                    value
-                }
+                y.numbers = x.doubles.map { spectrumC(it, args) }
             }
             scatter {
-                name = "CustomSterileNeutrinoSpectrum"
+                name = "CustomSterileNeutrinoSpectrum (rearWall=${argsRWall[rearWall]})"
                 mode = ScatterMode.lines
                 x.buffer = range
-                y.numbers = x.doubles.map { wallStep(it) * 0.02 }
+                y.numbers = x.doubles.map { spectrumC(it, argsRWall) }
             }
             layout {
-                title = "Sterile neutrino spectrum"
-                yaxis.type = AxisType.log
+                title = "Compare Trap"
+//                yaxis.type = AxisType.log
             }
         }
     }.makeFile()
+
+    println("HV\trearWall=${args[rearWall]}\trearWall=${argsRWall[rearWall]}")
+    range.iterator().forEach {
+        println("${it}\t${spectrumC(it, args)}\t${spectrumC(it, argsRWall)}")
+    }
 }