feat: implement NumassBetaWithBump model

numass-workspace/build.gradle.kts

@@ -39,11 +39,6 @@ graalvmNative {
 
     // ← добавляем блок с настройками для нескольких бинарников
     binaries {
-        named("main") {
-            imageName = "mass-app"
-            mainClass = "ru.inr.mass.scripts.ApplicationKt"
-        }
-
         // создаём дополнительные бинарники
         register("grid-fit") {
             imageName = "grid-fit"
@@ -78,20 +73,12 @@ tasks.register<CreateStartScripts>("fitCustomStartScripts") {
     outputDir = layout.buildDirectory.dir("scripts-fit").get().asFile
 }
 
-tasks.register<CreateStartScripts>("appStartScripts") {
-    applicationName = "mass-app"
-    mainClass.set("ru.inr.mass.scripts.ApplicationKt")
-    classpath = files(tasks.named<Jar>("jar"), configurations.runtimeClasspath)
-    outputDir = layout.buildDirectory.dir("scripts-app").get().asFile
-}
-
 // Добавляем все скрипты в дистрибутив
 distributions.named("main") {
     contents {
         into("bin") {
             from(tasks.named("gridFitStartScripts"))
             from(tasks.named("fitCustomStartScripts"))
-            from(tasks.named("appStartScripts"))
 //            fileMode = "0755"
             duplicatesStrategy = DuplicatesStrategy.EXCLUDE
         }

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

@@ -0,0 +1,96 @@
+@file:Suppress("PARAMETER_NAME_CHANGED_ON_OVERRIDE")
+
+package ru.inr.mass.scripts
+
+import ru.inr.mass.models.DifferentiableKernel
+import ru.inr.mass.models.Kernel
+import ru.inr.mass.models.NumassBeta
+import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.expressions.symbol
+
+/**
+ * Кастомный бета-спектр Трития:
+ *
+ * NumassBetaWithBump = NumassBeta(args) + A · NumassBeta(args с заменой mnu2 → B)
+ *
+ * Все частные производные реализованы корректно для фита.
+ */
+ class NumassBetaWithBump : DifferentiableKernel {
+    override val x: Symbol = NumassBeta.x
+    override val y: Symbol = NumassBeta.y
+
+    override fun invoke(fs: Double, eIn: Double, arguments: Map<Symbol, Double>): Double {
+        val args1 = arguments
+
+        // Второй вызов — переопределяем только mnu2 на B
+        val bVal = arguments.getValue(B)
+        val args2: Map<Symbol, Double> = arguments.toMutableMap().apply {
+            put(NumassBeta.mnu2, bVal)
+        }
+
+        val term1 = NumassBeta(fs, eIn, args1)
+        val term2 = NumassBeta(fs, eIn, args2)
+
+        return term1 + arguments.getValue(A) * term2
+    }
+
+    override fun derivativeOrNull(symbols: List<Symbol>): Kernel? = when (symbols.size) {
+        0 -> this
+        1 -> {
+            val sym = symbols.first()
+            Kernel { fs, eIn, arguments ->
+                val aVal = arguments.getValue(A)
+                val bVal = arguments.getValue(B)
+
+                val args1 = arguments
+                val args2 = arguments.toMutableMap().apply { this[NumassBeta.mnu2] = bVal }
+
+                when (sym) {
+                    // ∂/∂A = второй спектр
+                    A -> NumassBeta(fs, eIn, args2)
+
+                    // ∂/∂B = A · (∂NumassBeta/∂mnu2) при mnu2 = B
+                    B -> {
+                        val derivMnu2 = NumassBeta.derivativeOrNull(listOf(NumassBeta.mnu2))!!
+                        aVal * derivMnu2(fs, eIn, args2)
+                    }
+
+                    // Параметры первого компонента (обычный NumassBeta)
+                    NumassBeta.e0 -> {
+                        val d1 = NumassBeta.derivativeOrNull(listOf(NumassBeta.e0))!!.invoke(fs, eIn, args1)
+                        val d2 = NumassBeta.derivativeOrNull(listOf(NumassBeta.e0))!!.invoke(fs, eIn, args2)
+                        d1 + aVal * d2
+                    }
+
+                    NumassBeta.mnu2 -> {
+                        // Второй компонент не зависит от исходного mnu2
+                        NumassBeta.derivativeOrNull(listOf(NumassBeta.mnu2))!!.invoke(fs, eIn, args1)
+                    }
+
+                    NumassBeta.msterile2 -> {
+                        val d1 = NumassBeta.derivativeOrNull(listOf(NumassBeta.msterile2))!!.invoke(fs, eIn, args1)
+                        val d2 = NumassBeta.derivativeOrNull(listOf(NumassBeta.msterile2))!!.invoke(fs, eIn, args2)
+                        d1 + aVal * d2
+                    }
+
+                    NumassBeta.u2 -> {
+                        val d1 = NumassBeta.derivativeOrNull(listOf(NumassBeta.u2))!!.invoke(fs, eIn, args1)
+                        val d2 = NumassBeta.derivativeOrNull(listOf(NumassBeta.u2))!!.invoke(fs, eIn, args2)
+                        d1 + aVal * d2
+                    }
+
+                    else -> 0.0
+                }
+            }
+        }
+        else -> null
+    }
+
+    companion object {
+        /** Амплитуда дополнительного компонента */
+        val A: Symbol by symbol
+
+        /** Фиксированное значение mnu² для дополнительного компонента */
+        val B: Symbol by symbol
+    }
+}

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

@@ -0,0 +1,80 @@
+package ru.inr.mass.scripts
+
+import ru.inr.mass.models.DifferentiableKernel
+import ru.inr.mass.models.DifferentiableSpectrum
+import ru.inr.mass.models.NBkgSpectrum
+import ru.inr.mass.models.NumassBeta
+import ru.inr.mass.models.NumassResolution
+import ru.inr.mass.models.NumassTransmission
+import ru.inr.mass.models.Spectrum
+import ru.inr.mass.models.SterileNeutrinoSpectrum
+import ru.inr.mass.models.SterileNeutrinoSpectrumFast
+import ru.inr.mass.workspace.fitWith
+import ru.inr.mass.workspace.freeParameters
+import ru.inr.mass.workspace.iterations
+import space.kscience.kmath.UnstableKMathAPI
+import space.kscience.kmath.data.XYErrorColumnarData
+import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.expressions.symbol
+import space.kscience.kmath.misc.Loggable
+import space.kscience.kmath.optimization.OptimizationLog
+import space.kscience.kmath.optimization.QowOptimizer
+import space.kscience.kmath.optimization.XYFit
+
+class SterileNeutrinoSpectrumRWall(
+    private val wallModel: WallModel,
+    private val trapModel: ITrap,
+    private val adiabaticityModel: AdiabaticityModel,
+    source: DifferentiableKernel = NumassBeta,
+    parallel: Boolean = true,
+) : DifferentiableSpectrum {
+
+    val spectrum = when (parallel) {
+        true -> SterileNeutrinoSpectrumFast(
+            source = source, transmission = NumassTransmission(
+                trapFunc = { ei, ef, _ -> trapModel.value(ei, ei - ef) },
+                adjustX = false,
+            ), resolution = NumassResolution(1.7e-4, tailFunction = { e, u -> adiabaticityModel(e, u) }), fss = null
+        )
+
+        false -> SterileNeutrinoSpectrum(
+            source = source, transmission = NumassTransmission(
+                trapFunc = { ei, ef, _ -> trapModel.value(ei, ei - ef) },
+                adjustX = false,
+            ), resolution = NumassResolution(1.7e-4, tailFunction = { e, u -> adiabaticityModel(e, u) }), fss = null
+        )
+    }
+
+    override fun invoke(u: Double, arguments: Map<Symbol, Double>): Double {
+        val rw = arguments[rearWall]!!
+        return spectrum.invoke(u, arguments) * (1.0 + rw * wallModel(u))
+    }
+
+    override fun derivativeOrNull(symbols: List<Symbol>): Spectrum? {
+        if (symbols.isEmpty()) return this
+        return when (symbols.singleOrNull() ?: TODO("First derivatives only")) {
+            rearWall -> Spectrum { u, args -> wallModel(u) * spectrum.invoke(u, args) }
+            else -> spectrum.derivativeOrNull(symbols)
+        }
+    }
+
+    companion object {
+        val rearWall: Symbol by symbol
+
+        @OptIn(UnstableKMathAPI::class)
+        suspend fun fit(
+            spectrum: NBkgSpectrum,
+            data: XYErrorColumnarData<Double, Double, Double>,
+            fitParams: List<Symbol>,
+            initial: Map<Symbol, Double>,
+            logger: Loggable
+        ): XYFit {
+            return data.fitWith(
+                optimizer = QowOptimizer, modelExpression = spectrum, startingPoint = initial, attributesBuilder = {
+                    freeParameters(*fitParams.toTypedArray())
+                    iterations(20)
+                    OptimizationLog(logger)
+                })
+        }
+    }
+}

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

@@ -1,101 +0,0 @@
-package ru.inr.mass.scripts
-
-import com.charleskorn.kaml.Yaml
-import com.charleskorn.kaml.decodeFromStream
-import kotlinx.coroutines.runBlocking
-import kotlinx.serialization.Serializable
-import ru.inr.mass.models.NBkgSpectrum
-import ru.inr.mass.models.NumassBeta
-import ru.inr.mass.models.NumassTransmission
-import space.kscience.kmath.expressions.Symbol
-import java.io.File
-import kotlin.io.path.Path
-import kotlin.io.path.exists
-
-// TODO: sync params with fit-custom
-@Serializable
-data class FitParams(
-    val norm: Double = 254159.0,
-    val bkg: Double = 0.2,
-    val mnu2: Double = 0.0,
-    val e0: Double = 18569.0,
-    val msterile2: Double = 36e6,
-    val u2: Double = 0.0,
-    val thickness: Double = 0.3,
-    val trap: Double = 1.0,
-    val rearWall: Double = 0.01
-)
-
-@Serializable
-data class FitCustomParams(
-    val spectrum: String = "TABLE_SLICE.tsv",
-    val full: String = "TABLE_FULL.tsv",
-    val postfix: String = "sample",
-    val fitParams: FitParams = FitParams(),
-    val fixBkg: Boolean = false,
-    val fixU2: Boolean = false,
-    val fixTrap: Boolean = false,
-    val trapFunc: TrapModel = TrapModel.FINE_5DAY,
-    val fixRWall: Boolean = false,
-    val fixE0: Boolean = false,
-    val wallFunc: WallModel = WallModel.TSV_2025_12_01,
-    val noAdiabacity: Boolean = false,
-    val adiabacityFunc: AdiabaticityModel = AdiabaticityModel.`2024_11`
-)
-
-fun main (args: Array<String>) {
-    val usage = """
-        Usage: 
-           application CONFIG_PATH # start fitting with configuration from CONFIG_PATH
-           application sample # prints configuration template
-        """.trimIndent()
-
-    if(args.size != 1) {
-        println(usage)
-        return
-    }
-    if(args[0] == "sample") {
-        val yaml = Yaml.default.encodeToString(FitCustomParams.serializer(), FitCustomParams())
-        println(yaml)
-        return
-    }
-    if (!Path(args[0]).exists()) {
-        println(usage)
-        return
-    }
-
-    val config = Yaml.default.decodeFromStream(FitCustomParams.serializer(), File(args[0]).inputStream())
-
-    println("using wall=${config.wallFunc}")
-    println("using trap=${config.trapFunc}")
-    println("using adiabacity=${config.adiabacityFunc}")
-
-    val fitParams: Map<Symbol, Double> = mapOf(
-        NBkgSpectrum.norm to config.fitParams.norm,
-        NBkgSpectrum.bkg to config.fitParams.bkg,
-        NumassBeta.mnu2 to config.fitParams.mnu2,
-        NumassBeta.e0 to config.fitParams.e0,
-        NumassBeta.msterile2 to config.fitParams.msterile2,
-        NumassBeta.u2 to config.fitParams.u2,
-        NumassTransmission.thickness to config.fitParams.thickness,
-        NumassTransmission.trap to config.fitParams.trap,
-        rearWall to config.fitParams.rearWall
-    )
-
-    runBlocking {
-        processCustom(
-            spectrumFile =  config.spectrum,
-            full =config.full,
-            postfix = config.postfix,
-            fitParams = fitParams,
-            fixBkg = config.fixBkg,
-            fixU2 = config.fixU2,
-            fixTrap =  config.fixTrap,
-            trapFunc = config.trapFunc,
-            fixRwall = config.fixRWall,
-            fixE0 = config.fixE0,
-            wallFunc = config.wallFunc,
-            adiabacityFunc = config.adiabacityFunc
-        )
-    }
-}

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

@@ -0,0 +1,107 @@
+package ru.inr.mass.scripts
+
+import ru.inr.mass.scripts.NumassBetaWithBump.Companion.A
+import ru.inr.mass.scripts.NumassBetaWithBump.Companion.B
+import space.kscience.plotly.models.AxisType
+import kotlinx.html.br
+import kotlinx.html.p
+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.scripts.SterileNeutrinoSpectrumRWall.Companion.rearWall
+import ru.inr.mass.workspace.buffer
+import space.kscience.kmath.UnstableKMathAPI
+import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.real.step
+import space.kscience.plotly.*
+import space.kscience.plotly.models.ScatterMode
+import kotlin.io.path.Path
+
+@UnstableKMathAPI
+fun main() {
+    val range = 2000.0..18600.0 step 50.0
+
+    // Параметры взяты из случайного (00.5- 3_17e-05) фита grid-fit
+    val trapModel = TrapModel.LEGACY_2024_07_12
+    val wallModel = WallModel.GEANT_2025_10_09
+    val adiabaticityModel = AdiabaticityModel.vasya
+    val args: Map<Symbol, Double> = mapOf(
+        norm to 403098.17380392,
+        bkg to 0.6671743523307374,
+        mnu2 to 0.0,
+        e0 to 18562.073816431806,
+        msterile2 to 250000.0,
+        u2 to 1.261914688960386E-4,
+        thickness to 0.1,
+        trap to 1.0532058581987842,
+        rearWall to 0.049531505450610504
+    )
+
+    val argsBump = args.toMutableMap()
+    argsBump[A] = 7e-2
+    argsBump[B] = 10.4
+
+    val spectrum = constructSpectrum(trapModel, wallModel, adiabaticityModel)
+    val spectrumWithBump = constructSpectrum(trapModel, wallModel, adiabaticityModel, NumassBetaWithBump())
+
+    Plotly.page {
+        plot {
+            scatter {
+                name = "original"
+                mode = ScatterMode.lines
+                x.buffer = range
+                y.numbers = x.doubles.map { spectrum(it, args) }
+            }
+            scatter {
+                name = "with bump"
+                mode = ScatterMode.lines
+                x.buffer = range
+                y.numbers = x.doubles.map { spectrumWithBump(it, argsBump) }
+            }
+            layout {
+                title = "Rear Wall Contribution"
+                yaxis.type = AxisType.log
+            }
+        }
+        p {
+            +"command args = $args"
+        }
+        br()
+        p {
+            +"trapping model = $trapModel"
+        }
+        br()
+        p {
+            +"wall model = $wallModel"
+        }
+        br()
+    }.makeFile(Path("./beta-comparison.html"))
+
+    println("HV\toriginal\twith bump")
+    range.iterator().forEach {
+        println(
+            "${it}\t${spectrum(it, args)}\t${spectrumWithBump(it, argsBump)}"
+        )
+    }
+}
+
+private fun constructSpectrum(
+    trapModel: TrapModel,
+    wallModel: WallModel,
+    adiabaticityModel: AdiabaticityModel,
+    source: DifferentiableKernel = NumassBeta
+): NBkgSpectrum {
+    return SterileNeutrinoSpectrumRWall(
+        trapModel = trapModel,
+        wallModel = wallModel,
+        adiabaticityModel = adiabaticityModel,
+        source = source,
+        parallel = true
+    ).withNBkg()
+}

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

@@ -12,110 +12,68 @@ import com.github.ajalt.clikt.parameters.options.option
 import com.github.ajalt.clikt.parameters.types.double
 import com.github.ajalt.clikt.parameters.types.enum
 import kotlinx.coroutines.runBlocking
-import kotlinx.html.*
+import kotlinx.html.h3
+import kotlinx.html.h4
+import kotlinx.html.p
+import kotlinx.html.pre
 import ru.inr.mass.models.*
+import ru.inr.mass.scripts.SterileNeutrinoSpectrumRWall.Companion.rearWall
 import ru.inr.mass.workspace.buffer
-import ru.inr.mass.workspace.fitWith
-import ru.inr.mass.workspace.freeParameters
-import ru.inr.mass.workspace.iterations
 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.data.XYErrorColumnarData
 import space.kscience.kmath.data.indices
 import space.kscience.kmath.expressions.Symbol
-import space.kscience.kmath.expressions.symbol
 import space.kscience.kmath.misc.Loggable
+import space.kscience.kmath.operations.toTypedArray
 import space.kscience.kmath.optimization.*
-import space.kscience.kmath.structures.toDoubleArray
 import space.kscience.plotly.*
 import space.kscience.plotly.models.ScatterMode
-import java.io.File
-import java.nio.file.Path
 import kotlin.io.path.Path
 import kotlin.io.path.nameWithoutExtension
-import kotlin.math.pow
 
-val rearWall: Symbol by symbol
-
-
-class SterileNeutrinoSpectrumRWallMP(
-    source: DifferentiableKernel = NumassBeta,
-    transmission: DifferentiableKernel = NumassTransmission(),
-    resolution: DifferentiableKernel = NumassResolution(),
-    val wallFunc: WallModel,
-    fss: FSS? = null
-) : SterileNeutrinoSpectrumFast(source, transmission, resolution, fss) {
-
-    override fun invoke(u: Double, arguments: Map<Symbol, Double>): Double {
-        val rearWall = arguments[rearWall]!!
-        return super.invoke(u, arguments) * (1.0 + rearWall * wallFunc(u))
-    }
-
-    override fun derivativeOrNull(symbols: List<Symbol>): Spectrum? {
-        if (symbols.isEmpty()) return this
-        return when (symbols.singleOrNull() ?: TODO("First derivatives only")) {
-            rearWall -> Spectrum { u, arguments ->
-                wallFunc(u) * super.invoke(u, arguments)
-            }
-
-            else -> super.derivativeOrNull(symbols)
-        }
-    }
-}
-
-private class CustomArgs : CliktCommand() {
+private class FitCustom : CliktCommand(name = "fit-custom") {
 
     init {
-        context {
-            helpFormatter = { MordantHelpFormatter(it, showDefaultValues = true) }
-        }
+        context { helpFormatter = { MordantHelpFormatter(it, showDefaultValues = true) } }
     }
 
-    val spectrum: String by argument().help("path to spectrum TSV file")
-    val full: String? by option().help("path to full spectrum (for test)")
-    val postfix: String? by option().help("save file postfix")
+    val spectrum by argument().help("path to spectrum TSV file")
+    val full by option().help("path to full spectrum (for residuals/test)")
+    val postfix by option().help("save file postfix")
 
-    val norm: Double by option().double().default(8.311751921319484E8).help("NBkgSpectrum.norm")
+    val norm by option().double().default(8.311751921319484E8).help("NBkgSpectrum.norm")
+    val bkg by option().double().default(1612.626961212948).help("NBkgSpectrum.bkg")
+    val fixBkg by option().flag().help("Don't fit bkg")
 
-    val bkg: Double by option().double().default(1612.626961212948).help("NBkgSpectrum.bkg")
-    val fixBkg: Boolean by option().flag().help("Don't fit bkg variable")
+    val mnu2 by option().double().default(0.0).help("NumassBeta.mnu2")
 
-    val mnu2: Double by option().double().default(0.0).help("NumassBeta.mnu2")
+    val e0 by option().double().default(18572.0).help("NumassBeta.e0")
+    val fixE0 by option().flag().help("Don't fit e0")
 
-    val e0: Double by option().double().default(18572.0).help("NumassBeta.e0")
-    val fixE0: Boolean by option().flag().help("Don't fit e0 variable")
+    val msterile2 by option().double().default(0.0).help("NumassBeta.msterile2")
 
-    val msterile2: Double by option().double().default(0.0.pow(2)).help("NumassBeta.msterile2")
+    val u2 by option().double().default(1e-2).help("NumassBeta.u2")
+    val fixU2 by option().flag().help("Don't fit u2")
 
-    val u2: Double by option().double().default(1e-2).help("NumassBeta.u2")
-    val fixU2: Boolean by option().flag().help("Don't fit u2 variable")
+    val thickness by option().double().default(0.3).help("NumassTransmission.thickness")
 
-    val thickness: Double by option().double().default(0.3).help("NumassTransmission.thickness")
+    val trap by option().double().default(1.0).help("NumassTransmission.trap")
+    val trapFunc by option().enum<TrapModel>().default(TrapModel.FINE_5DAY).help("trap model")
+    val fixTrap by option().flag().help("Don't fit trap")
 
-    val trap: Double by option().double().default(1.0).help("NumassTransmission.trap")
-    val trapFunc: TrapModel by option().enum<TrapModel>().default(TrapModel.FINE_5DAY)
-        .help(
-            "trap model"
-        )
-    val fixTrap: Boolean by option().flag().help("Don't fit trap variable")
+    val rwall by option().double().default(0.17).help("rear wall factor")
+    val wallFunc by option().enum<WallModel>().default(WallModel.TSV_2025_12_01).help("wall model")
+    val fixRwall by option().flag().help("Don't fit rwall")
 
-
-    val rwall: Double by option().double().default(0.17).help("rear wall")
-    val wallFunc: WallModel by option().enum<WallModel>().default(WallModel.TSV_2025_12_01)
-        .help(
-            "wall model"
-        )
-    val fixRwall: Boolean by option().flag().help("Don't fit rwall variable")
-
-    // TODO: convert to enum
-    val adiabacityFunc: AdiabaticityModel by option().enum<AdiabaticityModel>().default(AdiabaticityModel.`2024_11`)
-        .help(
-            "adiabacity model"
-        )
+    val adiabacityFunc by option().enum<AdiabaticityModel>().default(AdiabaticityModel.`2024_11`)
+        .help("adiabacity tail model")
 
     @UnstableKMathAPI
-    override fun run() {
+    override fun run() = runBlocking {
+        val data = parse(spectrum)
+        val testData = full?.let { parse(it) } ?: data
 
-        val fitParams: Map<Symbol, Double> = mapOf(
+        val initialParams = mapOf(
             NBkgSpectrum.norm to norm,
             NBkgSpectrum.bkg to bkg,
             NumassBeta.mnu2 to mnu2,
@@ -127,276 +85,110 @@ private class CustomArgs : CliktCommand() {
             rearWall to rwall
         )
 
-        runBlocking {
-            processCustom(
-                spectrum,
-                full,
-                postfix,
-                fitParams,
-                fixBkg,
-                fixU2,
-                fixTrap,
-                fixE0,
-                trapFunc,
-                fixRwall,
-                wallFunc,
-                cliArgs = this@CustomArgs,
-                adiabacityFunc = adiabacityFunc
-            )
+        val spectrumModel = SterileNeutrinoSpectrumRWall(
+            wallModel = wallFunc, trapModel = trapFunc, adiabaticityModel = adiabacityFunc, parallel = true
+        ).withNBkg()
+
+        val fitVars = mutableListOf(NBkgSpectrum.norm)
+        if (!fixE0) fitVars += NumassBeta.e0
+        if (!fixBkg) fitVars += NBkgSpectrum.bkg
+        if (!fixU2) fitVars += NumassBeta.u2
+        if (!fixTrap) fitVars += NumassTransmission.trap
+        if (!fixRwall) fitVars += rearWall
+
+        val logMessages = mutableListOf<String>()
+        val logger = Loggable { tag, msg ->
+            val line = "[$tag] ${msg()}"
+            logMessages += line
+            println(line)
         }
+
+        val fit = SterileNeutrinoSpectrumRWall.fit(
+            spectrumModel,
+            data,
+            fitVars,
+            initialParams,
+            logger
+        )
+
+        val dataPath = Path(spectrum)
+        val outName = dataPath.nameWithoutExtension + (postfix?.let { "-$it" } ?: "") + ".html"
+        val reportPath = dataPath.parent.resolve(outName)
+
+        buildReportPage(
+            spectrumFile = spectrum,
+            postfix = postfix,
+            data = data,
+            testData = testData,
+            model = spectrumModel,
+            initial = initialParams,
+            fit = fit,
+            logMessages = logMessages,
+            cliArgs = this@FitCustom
+        ).makeFile(reportPath, show = false)
+
+        println("Report saved → $reportPath")
+        openInBrowser("file:///$reportPath")
     }
 }
 
-suspend fun fitNumassSpectrumRWall(
-    spectrum: NBkgSpectrum,
-    data: XYErrorColumnarData<Double, Double, Double>,
-    fitParams: List<Symbol>,
-    initial: Map<Symbol, Double>,
-    logger: Loggable
-): XYFit {
-    return data.fitWith(
-        optimizer = QowOptimizer, modelExpression = spectrum, startingPoint = initial, attributesBuilder = {
-            freeParameters(*fitParams.toTypedArray())
-            iterations(20)
-            OptimizationLog(logger)
-        })
-}
-
 @UnstableKMathAPI
-suspend fun processCustom(
+fun buildReportPage(
     spectrumFile: String,
-    full: String?,
     postfix: String?,
-    fitParams: Map<Symbol, Double>,
-    fixBkg: Boolean,
-    fixU2: Boolean,
-    fixTrap: Boolean,
-    fixE0: Boolean,
-    trapFunc: ITrap,
-    fixRwall: Boolean,
-    wallFunc: WallModel,
-    cliArgs: Any? = null,
-    adiabacityFunc: AdiabaticityModel
-) {
-    println(cliArgs)
-
-    val postfix = if (postfix != null) "-$postfix" else ""
-
-    val data = parse(spectrumFile)
-    val testData = if (full != null) parse(full) else data
-
-    val spectrum: NBkgSpectrum = SterileNeutrinoSpectrumRWallMP(
-//        fss = FSS.default,
-        wallFunc = wallFunc,
-        transmission = NumassTransmission(
-//            NumassTransmission.trapFunction,
-            trapFunc = { ei, ef, _ ->
-                val delta = ei - ef
-                trapFunc.value(ei, delta)
-            },
-            adjustX = false,
-        ), resolution = NumassResolution(1.7e-4, tailFunction = { e, u ->
-            adiabacityFunc(e, u)
-        })
-    ).withNBkg()
-
-    val logMessages = emptyList<String>().toMutableList()
-    val dumpLogger = Loggable { tag, block ->
-        logMessages += "[$tag] ${block()}"
-        println("[$tag] ${block()}")
-    }
-
-    val fitVars = mutableListOf(
-        NBkgSpectrum.norm,
-//        NBkgSpectrum.bkg,
-//        rearWall,
-//        NumassBeta.u2,
-//        NumassBeta.mnu2,
-//        NumassBeta.msterile2,
-//        NumassTransmission.trap,
-    )
-
-    if (!fixE0) {
-        fitVars += NumassBeta.e0
-    }
-
-    if (!fixBkg) {
-        fitVars += NBkgSpectrum.bkg
-    }
-
-    if (!fixU2) {
-        fitVars += NumassBeta.u2
-    }
-
-    if (!fixTrap) {
-        fitVars += NumassTransmission.trap
-    }
-
-    if (!fixRwall) {
-        fitVars += rearWall
-    }
-
-    val fit = fitNumassSpectrumRWall(spectrum, data, fitVars, fitParams, dumpLogger)
-
-    val dataPath = Path(spectrumFile)
-    val reportPath = Path(
-        dataPath.parent.toString(), dataPath.nameWithoutExtension + postfix + ".html"
-    )
-
-    val report =
-        reportPage(spectrumFile, postfix, data, spectrum, fitParams, fit, testData, dataPath, cliArgs, logMessages)
-
-    report.makeFile(path = reportPath, show = false)
-    openInBrowser("file:///$reportPath")
-}
-
-@UnstableKMathAPI
-fun reportPage(
-    spectrumFile: String,
-    postfix: String,
     data: XYErrorColumnarData<Double, Double, Double>,
-    spectrum: NBkgSpectrum,
-    fitParams: Map<Symbol, Double>,
-    fit: XYFit,
     testData: XYErrorColumnarData<Double, Double, Double>,
-    dataPath: Path,
-    cliArgs: Any?,
-    logMessages: MutableList<String>
-): PlotlyPage {
-    val report = Plotly.page {
-        h3 {
-            +"Fit for $spectrumFile${postfix.drop(1).takeIf { it.isNotEmpty() }?.let { " ($it)" } ?: ""}"
-        }
-        plot {
-            scatter {
-                name = "Data"
-                mode = ScatterMode.markers
-                x.buffer = data.x
-                y.buffer = data.y
-            }
-            scatter {
-                name = "Initial"
-                mode = ScatterMode.lines
-                x.buffer = data.x
-                y.numbers = x.doubles.map { spectrum(it, fitParams + fit.result) }
-            }
-            scatter {
-                name = "Fit"
-                mode = ScatterMode.lines
-                x.buffer = data.x
-                y.numbers = data.x.toDoubleArray().map { spectrum(it, fitParams + fit.result) }
-            }
-        }
-        plot {
-            layout {
-                title = "Residuals"
-            }
-            scatter {
-                name = "Residuals"
-                mode = ScatterMode.markers
-                x.buffer = testData.x
-                y.numbers = testData.indices.map {
-                    val value = spectrum(testData.x[it], fitParams + fit.result)
-                    val dif = testData.y[it] - value
-                    dif / testData.yErr[it]
-                }
+    model: NBkgSpectrum,
+    initial: Map<Symbol, Double>,
+    fit: XYFit,
+    logMessages: List<String>,
+    cliArgs: Any?
+): PlotlyPage = Plotly.page {
+    h3 { +"Fit for $spectrumFile${postfix?.let { " ($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 {
+        scatter { name = "Data"; mode = ScatterMode.markers; x.buffer = data.x; y.buffer = data.y }
+        scatter {
+            name = "Initial"; mode = ScatterMode.lines; x.buffer = data.x; y.numbers =
+            data.x.toTypedArray().map { model(it, initial) }
         }
-        plot {
-            layout {
-                title = "Residuals distribution"
-            }
-            histogram {
-                x.numbers = data.indices.map {
-                    val value = spectrum(data.x[it], fitParams + fit.result)
-                    val dif = data.y[it] - value
-                    dif / data.yErr[it]
-                }
-                name = "Res histo"
-            }
+        scatter {
+            name = "Fit"; mode = ScatterMode.lines; x.buffer = data.x; y.numbers =
+            data.x.toTypedArray().map { model(it, initial + fit.result) }
         }
-        p {
-            +"command args = $cliArgs"
-        }
-        br()
-        p {
-            +"initial params = $fitParams"
-        }
-        br()
-        p {
-            +"fit result = ${fit.result}"
-        }
-        br()
-        p {
-            +"chi^2 = ${fit.chiSquaredOrNull}"
-        }
-        br()
-        p {
-            +"dof = ${fit.dof}"
-        }
-        h4 {
-            +"Fit log:"
-        }
-        pre {
-            logMessages.forEach {
-                +it
-                +"\n"
-            }
-        }
-        h4 {
-            +"Residuals:"
-        }
-        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
-                        +(dif / testData.yErr[it]).toString()
-                    }
-                }
-            }
-        }
-        h4 {
-            +"Data:"
-        }
-        table {
-            tr {
-                th { +"U_sp" }
-                th { +"counts" }
-                th { +"residials" }
-            }
-            data.indices.forEach {
-                tr {
-                    td { +data.x[it].toString() }
-                    td { +data.y[it].toString() }
-                    td { +data.yErr[it].toString() }
-                }
+    }
+
+    plot {
+        layout.title = "Residuals (σ)"
+        scatter {
+            name = "Residuals"
+            mode = ScatterMode.markers
+            x.buffer = testData.x
+            y.numbers = testData.indices.map {
+                val pred = model(testData.x[it], initial + fit.result)
+                (testData.y[it] - pred) / testData.yErr[it]
             }
         }
     }
-    return report
+
+    plot {
+        layout.title = "Residuals distribution"
+        histogram {
+            name = "Res histo"
+            x.numbers = data.indices.map {
+                val pred = model(data.x[it], initial + fit.result)
+                (data.y[it] - pred) / data.yErr[it]
+            }
+        }
+    }
+
+    p { +"CLI: $cliArgs" }
+    p { +"Initial: $initial" }
+    p { +"Result: ${fit.result}" }
+    p { +"χ² = ${fit.chiSquaredOrNull ?: "—"},  dof = ${fit.dof}" }
+
+    h4 { +"Log" }
+    pre { logMessages.forEach { +"$it\n" } }
 }
 
-fun main(args: Array<String>) = CustomArgs().main(args)
+fun main(args: Array<String>) = FitCustom().main(args)

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

@@ -12,6 +12,7 @@ import com.github.ajalt.clikt.parameters.types.enum
 import kotlinx.coroutines.Dispatchers
 import kotlinx.coroutines.runBlocking
 import ru.inr.mass.models.*
+import ru.inr.mass.scripts.SterileNeutrinoSpectrumRWall.Companion.rearWall
 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.operations.toTypedArray
@@ -55,12 +56,6 @@ private class CalcSpectrum : CliktCommand() {
             "wall model"
         )
 
-    val fixRwall: Boolean by option().flag().help("Don't fit rwall variable")
-
-    val noAdiabacity: Boolean by option().flag().help("Use const 1 instead of adiabacity function")
-
-    val fixE0: Boolean by option().flag().help("Don't fit e0 variable")
-
     // TODO: convert to enum
     val adiabacityFunc: AdiabaticityModel by option().enum<AdiabaticityModel>().default(AdiabaticityModel.`2024_11`)
         .help(
@@ -94,22 +89,15 @@ suspend fun calcSpectrum(fitParams: Map<Symbol, Double>, trapFunc: ITrap, wallFu
 
     runBlocking(Dispatchers.Default) {
 
-        val spectrum: NBkgSpectrum = SterileNeutrinoSpectrumRWallMP(
-            wallFunc = wallFunc,
-            transmission = NumassTransmission(
-                trapFunc = { ei, ef, _ ->
-                    val delta = ei - ef
-                    trapFunc.value(ei, delta)
-                },
-                adjustX = false,
-            ),
-            resolution = NumassResolution(1.7e-4, tailFunction = { e, u -> AdiabaticityModel.`2024_11`(e, u) })
+        val spectrum: NBkgSpectrum = SterileNeutrinoSpectrumRWall(
+            wallModel = wallFunc,
+            trapFunc,
+            AdiabaticityModel.`2024_11`,
+            parallel = true,
         ).withNBkg()
 
         val out = range.map {
 //            async(Dispatchers.Default) {
-
-
             val res = spectrum.invoke(it, fitParams)
 //                return@async res
 //            }

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

@@ -12,10 +12,9 @@ import com.github.ajalt.clikt.parameters.options.option
 import com.github.ajalt.clikt.parameters.options.multiple
 import com.github.ajalt.clikt.parameters.types.double
 import com.github.ajalt.clikt.parameters.types.enum
-import kotlinx.coroutines.coroutineScope
-import kotlinx.coroutines.launch
 import kotlinx.coroutines.runBlocking
 import ru.inr.mass.models.*
+import ru.inr.mass.scripts.SterileNeutrinoSpectrumRWall.Companion.rearWall
 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.misc.Loggable
@@ -27,33 +26,6 @@ import kotlin.io.path.createDirectories
 import kotlin.io.path.nameWithoutExtension
 import kotlin.math.log10
 import kotlin.math.pow
-import kotlinx.coroutines.Dispatchers
-import kotlinx.coroutines.withContext
-
-class SterileNeutrinoSpectrumRWall(
-    source: DifferentiableKernel = NumassBeta,
-    transmission: DifferentiableKernel = NumassTransmission(),
-    resolution: DifferentiableKernel = NumassResolution(),
-    val wallFunc: WallModel,
-    fss: FSS? = null
-) : SterileNeutrinoSpectrum(source, transmission, resolution, fss) {
-
-    override fun invoke(u: Double, arguments: Map<Symbol, Double>): Double {
-        val rearWall = arguments[rearWall]!!
-        return super.invoke(u, arguments) * (1.0 + rearWall * wallFunc(u))
-    }
-
-    override fun derivativeOrNull(symbols: List<Symbol>): Spectrum? {
-        if (symbols.isEmpty()) return this
-        return when (symbols.singleOrNull() ?: TODO("First derivatives only")) {
-            rearWall -> Spectrum { u, arguments ->
-                wallFunc(u) * super.invoke(u, arguments)
-            }
-
-            else -> super.derivativeOrNull(symbols)
-        }
-    }
-}
 
 private class GridArgs : CliktCommand() {
 
@@ -122,16 +94,8 @@ private class GridArgs : CliktCommand() {
             rearWall to rwall
         )
 
-        val spectrumModel: NBkgSpectrum = SterileNeutrinoSpectrumRWallMP(
-            wallFunc = wallFunc, transmission = NumassTransmission(
-                trapFunc = { ei, ef, _ ->
-                    val delta = ei - ef
-                    trapFunc.value(ei, delta)
-                },
-                adjustX = false,
-            ), resolution = NumassResolution(1.7e-4, tailFunction = { e, u ->
-                adiabacityFunc(e, u)
-            })
+        val spectrumModel: NBkgSpectrum = SterileNeutrinoSpectrumRWall(
+            wallModel = wallFunc, trapModel = trapFunc, adiabaticityModel = adiabacityFunc, parallel = true
         ).withNBkg()
 
         val fitVars = mutableListOf(NBkgSpectrum.norm)
@@ -231,7 +195,7 @@ private class GridArgs : CliktCommand() {
                         println("$current/$total [$taskId] $msg")
                     }
 
-                    val fit = fitNumassSpectrumRWall(
+                    val fit = SterileNeutrinoSpectrumRWall.fit(
                         spectrumModel, data, fitVars, fitParams, dumpLogger
                     )
 
@@ -243,17 +207,16 @@ private class GridArgs : CliktCommand() {
                         val reportFileName = "${mStr}-${u2Str}.html"
                         val reportPath = reportsDir.resolve(reportFileName)
 
-                        reportPage(
+                        buildReportPage(
                             spectrum,
                             postfix,
                             data,
+                            testData = dataForView,
                             spectrumModel,
                             fitParams,
                             fit,
-                            dataForView,
-                            dataPath,
+                            logMessages,
                             this@GridArgs,
-                            logMessages
                         ).makeFile(path = reportPath, show = false)
                     }
 

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

@@ -11,6 +11,7 @@ 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.scripts.SterileNeutrinoSpectrumRWall.Companion.rearWall
 import ru.inr.mass.workspace.buffer
 import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.real.step
@@ -19,8 +20,9 @@ import space.kscience.plotly.models.ScatterMode
 
 fun main() {
 
-    val trapInterpolator = TrapModel.FINE_5DAY
-    val wallInterpolator = WallModel.TSV_2025_12_01
+    val trapModel = TrapModel.FINE_5DAY
+    val wallModel = WallModel.TSV_2025_12_01
+    val adiabacityModel = AdiabaticityModel.vasya
 
     val range = 2000.0..18600.0 step 50.0
 
@@ -45,19 +47,8 @@ fun main() {
     val rwall3 = args.toMutableMap();
     rwall3[rearWall] = 3.0
 
-    val spectrumC: NBkgSpectrum = SterileNeutrinoSpectrumRWallMP(
-        fss = null,
-        transmission = NumassTransmission(
-            trapFunc = { ei, ef, _ ->
-                val delta = ei - ef
-                trapInterpolator.value(ei, delta)
-            },
-            adjustX = false,
-        ),
-        wallFunc = wallInterpolator,
-        resolution = NumassResolution(1.7e-4, tailFunction = {
-                e,u -> AdiabaticityModel.vasya(e,u)
-        })
+    val spectrumC: NBkgSpectrum = SterileNeutrinoSpectrumRWall(
+        wallModel, trapModel, adiabacityModel
     ).withNBkg()
 
     Plotly.page {
@@ -86,11 +77,11 @@ fun main() {
         }
         br()
         p {
-            +"trapping model = $trapInterpolator"
+            +"trapping model = $trapModel"
         }
         br()
         p {
-            +"wall model = $wallInterpolator"
+            +"wall model = $wallModel"
         }
         br()
     }.makeFile()
@@ -103,6 +94,12 @@ fun main() {
         val r2 = spectrumC(it, rwall2)
         val r3 = spectrumC(it, rwall3)
 
-        println("${it}\t${spectrumC(it, args)}\t${spectrumC(it, rwall1)}\t${spectrumC(it, rwall2)}\t${spectrumC(it, rwall3)}")
+        println(
+            "${it}\t${spectrumC(it, args)}\t${spectrumC(it, rwall1)}\t${spectrumC(it, rwall2)}\t${
+                spectrumC(
+                    it, rwall3
+                )
+            }"
+        )
     }
 }

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

@@ -9,6 +9,7 @@ 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.scripts.SterileNeutrinoSpectrumRWall.Companion.rearWall
 import ru.inr.mass.workspace.buffer
 import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.real.step
@@ -16,8 +17,6 @@ import space.kscience.plotly.*
 import space.kscience.plotly.models.ScatterMode
 
 fun main() {
-
-    val trapInterpolator = TrapModel.FINE_5DAY
     val range = 2000.0..18600.0 step 50.0
 
     val args: Map<Symbol, Double> = mapOf(
@@ -41,19 +40,8 @@ fun main() {
     val argsTrap3 = args.toMutableMap();
     argsTrap3[trap] = 3.0
 
-    val spectrumC: NBkgSpectrum = SterileNeutrinoSpectrumRWallMP(
-        fss = null,
-        transmission = NumassTransmission(
-            trapFunc = { ei, ef, _ ->
-                val delta = ei - ef
-                trapInterpolator.value(ei, delta)
-            },
-            adjustX = false,
-        ),
-        wallFunc = WallModel.TSV_2025_12_01,
-        resolution = NumassResolution(1.7e-4, tailFunction = {
-                e,u -> AdiabaticityModel.vasya(e,u)
-        })
+    val spectrumC: NBkgSpectrum = SterileNeutrinoSpectrumRWall(
+        WallModel.TSV_2025_12_01, TrapModel.FINE_5DAY, AdiabaticityModel.vasya
     ).withNBkg()
 
     Plotly.page {
@@ -81,6 +69,12 @@ fun main() {
 
     println("HV\ttrap=${args[trap]}\ttrap=${argsTrap1[trap]}\ttrap=${argsTrap2[trap]}\ttrap=${argsTrap3[trap]}")
     range.iterator().forEach {
-        println("${it}\t${spectrumC(it, args)}\t${spectrumC(it, argsTrap1)}\t${spectrumC(it, argsTrap2)}\t${spectrumC(it, argsTrap3)}")
+        println(
+            "${it}\t${spectrumC(it, args)}\t${spectrumC(it, argsTrap1)}\t${
+                spectrumC(
+                    it, argsTrap2
+                )
+            }\t${spectrumC(it, argsTrap3)}"
+        )
     }
 }

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

@@ -9,6 +9,7 @@ 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.scripts.SterileNeutrinoSpectrumRWall.Companion.rearWall
 import ru.inr.mass.workspace.buffer
 import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.real.step
@@ -35,19 +36,8 @@ fun main() {
     val argsRWall = args.toMutableMap();
     argsRWall[rearWall] = 0.04
 
-    val spectrumC: NBkgSpectrum = SterileNeutrinoSpectrumRWallMP(
-        fss = null,
-        transmission = NumassTransmission(
-            trapFunc = { ei, ef, _ ->
-                val delta = ei - ef
-                trapInterpolator.value(ei, delta)
-            },
-            adjustX = false,
-        ),
-        wallFunc = WallModel.GEANT_2025_03_10,
-        resolution = NumassResolution(1.7e-4, tailFunction = {
-                e,u -> AdiabaticityModel.`2024_11`(e,u)
-        })
+    val spectrumC: NBkgSpectrum = SterileNeutrinoSpectrumRWall(
+        WallModel.GEANT_2025_03_10, trapInterpolator, AdiabaticityModel.`2024_11`
     ).withNBkg()
 
     Plotly.page {