diff --git a/numass-data-proto/src/gen/kotlin/ru/inr/mass/data/proto/Point.kt b/numass-data-proto/src/gen/kotlin/ru/inr/mass/data/proto/Point.kt index 1ccfaa2..d6266c8 100644 --- a/numass-data-proto/src/gen/kotlin/ru/inr/mass/data/proto/Point.kt +++ b/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.internal.redactElements +import com.squareup.wire.WireField +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( diff --git a/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassResolution.kt b/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassResolution.kt index 16ab8f4..3b572c0 100644 --- a/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassResolution.kt +++ b/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 { diff --git a/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassTransmission.kt b/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassTransmission.kt index 6323441..075ce99 100644 --- a/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassTransmission.kt +++ b/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>() - private const val ION_POTENTIAL = 15.4//eV + private const val ION_POTENTIAL = 13.6//eV private fun getX(arguments: Map, 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 diff --git a/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/fit.kt b/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/fit.kt index 45a9c4d..fa2b0b2 100644 --- a/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/fit.kt +++ b/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 { + val x: MutableList = mutableListOf() + val y: MutableList = mutableListOf() + val errors: MutableList = 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, -// resolution = NumassResolution(8.2e-5) +// fss = FSS.default, + transmission = NumassTransmission(NumassTransmission.trapFunction), + resolution = NumassResolution(1.7e-4) ).withNBkg() val args: Map = mapOf( - NBkgSpectrum.norm to 8e5, - NBkgSpectrum.bkg to 2.0, + NBkgSpectrum.norm to 319034.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 { - println("$it: ${spectrum.derivative(it).invoke(14000.0, args)}") - } +// listOf(NBkgSpectrum.norm, NBkgSpectrum.bkg, NumassBeta.e0).forEach { +// 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( - 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}") - 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 { - name = "Generated" - mode = ScatterMode.markers - x.buffer = generatedData.x - y.buffer = generatedData.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.forEach { + out.println(it) + } + } + } } - - scatter { - name = "Initial" - mode = ScatterMode.lines - x.buffer = 12000.0..18600.0 step 50.0 - y.numbers = x.doubles.map { spectrum(it, args) } +/* scatter { + name = "Fit" + mode = ScatterMode.lines + x.buffer = 12000.0..18600.0 step 10.0 + y.numbers = x.doubles.map { spectrum(it, args + fit.resultPoint) } + } } - - scatter { - name = "Fit" - mode = ScatterMode.lines - x.buffer = 12000.0..18600.0 step 10.0 - y.numbers = x.doubles.map { spectrum(it, args + fit.resultPoint) } + plot{ + scatter { + name = "Residuals" + mode = ScatterMode.markers + 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] + } + } } + 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() } \ No newline at end of file