feat: implement NumassBetaWithBump model
This commit is contained in:
@@ -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
|
||||
}
|
||||
}
|
||||
@@ -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()
|
||||
}
|
||||
Reference in New Issue
Block a user