Update spectrum computation

This commit is contained in:
Alexander Nozik 2021-05-18 18:19:48 +03:00
parent 014628f2fd
commit 0a5812f0fa
6 changed files with 563 additions and 36 deletions

View File

@ -19,12 +19,14 @@ import space.kscience.kmath.data.ColumnarData
import space.kscience.kmath.misc.Symbol
import space.kscience.kmath.misc.UnstableKMathAPI
import space.kscience.kmath.misc.symbol
import space.kscience.kmath.real.div
import space.kscience.kmath.real.sum
import space.kscience.kmath.structures.Buffer
import space.kscience.kmath.structures.DoubleBuffer
@OptIn(UnstableKMathAPI::class)
public class FSS(public val es: DoubleBuffer, public val ps: DoubleBuffer) : ColumnarData<Double> {
public class FSS(public val es: Buffer<Double>, public val ps: Buffer<Double>) : ColumnarData<Double> {
override val size: Int get() = ps.size
@ -47,7 +49,7 @@ public class FSS(public val es: DoubleBuffer, public val ps: DoubleBuffer) : Col
}.toList()
val es = DoubleBuffer(data.size) { data[it].first }
val ps = DoubleBuffer(data.size) { data[it].second }
FSS(es, ps)
FSS(es, ps / ps.sum())
}
}
}

View File

@ -42,5 +42,9 @@ public class NBkgSpectrum(public val source: Spectrum) : DifferentiableSpectrum
public val norm: Symbol by symbol
public val bkg: Symbol by symbol
}
}
/**
* Apply transformation adding norming-factor and the background
*/
public fun Spectrum.withNBkg(): NBkgSpectrum = NBkgSpectrum(this)

View File

@ -67,10 +67,10 @@ public class NumassTransmission(
else -> null
}
override fun invoke(x: Double, y: Double, arguments: Map<Symbol, Double>): Double {
override fun invoke(ei: Double, ef: Double, arguments: Map<Symbol, Double>): Double {
// loss part
val thickness = arguments[thickness] ?: 0.0
val loss = getTotalLossValue(thickness, x, y)
val loss = getTotalLossValue(thickness, ei, ef)
// double loss;
//
// if(eIn-eOut >= 300){
@ -84,7 +84,7 @@ public class NumassTransmission(
// }
//trapping part
val trap = arguments[trap] ?: 1.0 * trapFunc(x, y, arguments)
val trap = (arguments[trap] ?: 1.0) * trapFunc(ei, ef, arguments)
return loss + trap
}
@ -155,13 +155,11 @@ public class NumassTransmission(
* @param order
* @return
*/
private fun getLoss(order: Int): UnivariateFunction<Double> {
return getCachedSpectrum(order)
}
private fun getLoss(order: Int): UnivariateFunction<Double> = getCachedSpectrum(order)
private fun getLossProbDerivs(x: Double): List<Double> {
val res = ArrayList<Double>()
val probs = getLossProbabilities(x)
val probs = lossProbabilities(x)
var delta = exp(-x)
res.add((delta - probs[0]) / x)
@ -209,7 +207,7 @@ public class NumassTransmission(
return res
}
private fun getLossProbabilities(x: Double): List<Double> =
public fun lossProbabilities(x: Double): List<Double> =
lossProbCache.getOrPut(x) { calculateLossProbabilities(x) }
private fun getLossProbability(order: Int, X: Double): Double {
@ -220,7 +218,7 @@ public class NumassTransmission(
1.0
}
}
val probs = getLossProbabilities(X)
val probs = lossProbabilities(X)
return if (order >= probs.size) {
0.0
} else {
@ -228,11 +226,11 @@ public class NumassTransmission(
}
}
private fun getLossValue(order: Int, Ei: Double, Ef: Double): Double {
private fun getLossValue(order: Int, ei: Double, ef: Double): Double {
return when {
Ei - Ef < 5.0 -> 0.0
Ei - Ef >= getMargin(order) -> 0.0
else -> getLoss(order).invoke(Ei - Ef)
ei - ef < 5.0 -> 0.0
ei - ef >= getMargin(order) -> 0.0
else -> getLoss(order).invoke(ei - ef)
}
}
@ -303,17 +301,17 @@ public class NumassTransmission(
/**
* Значение полной функции потерь с учетом всех неисчезающих порядков
*
* @param x
* @param thickness
* @param Ei
* @param Ef
* @return
*/
private fun getTotalLossValue(x: Double, Ei: Double, Ef: Double): Double {
return if (x == 0.0) {
private fun getTotalLossValue(thickness: Double, Ei: Double, Ef: Double): Double {
return if (thickness == 0.0) {
0.0
} else {
val probs = getLossProbabilities(x)
(1 until probs.size).sumByDouble { i ->
val probs = lossProbabilities(thickness)
(1 until probs.size).sumOf { i ->
probs[i] * getLossValue(i, Ei, Ef)
}
}
@ -468,9 +466,8 @@ public class NumassTransmission(
//
// }
internal val defaultTrapping: Kernel = Kernel { e, u, arguments ->
val d = e - u
0.99797 - 3.05346E-7 * d - 5.45738E-10 * d.pow(2.0) - 6.36105E-14 * d.pow(3.0)
internal val defaultTrapping: Kernel = Kernel { ei, ef, _ ->
1.2e-4 - 4.5e-9 * ei
}
}

View File

@ -33,7 +33,6 @@ public class SterileNeutrinoSpectrum(
public val fss: FSS? = FSS.default,
) : DifferentiableSpectrum {
/**
* auxiliary function for trans-res convolution
*/

View File

@ -4,30 +4,94 @@ import inr.numass.models.sterile.NumassBeta.e0
import inr.numass.models.sterile.NumassBeta.mnu2
import inr.numass.models.sterile.NumassBeta.msterile2
import inr.numass.models.sterile.NumassBeta.u2
import inr.numass.models.sterile.NumassResolution
import inr.numass.models.sterile.NumassTransmission
import inr.numass.models.sterile.NumassTransmission.Companion.thickness
import inr.numass.models.sterile.NumassTransmission.Companion.trap
import inr.numass.models.sterile.SterileNeutrinoSpectrum
import kotlinx.html.code
import ru.inr.mass.models.NBkgSpectrum.Companion.bkg
import ru.inr.mass.models.NBkgSpectrum.Companion.norm
import ru.inr.mass.models.withNBkg
import ru.inr.mass.workspace.buffer
import space.kscience.kmath.misc.Symbol
import space.kscience.kmath.real.step
import space.kscience.plotly.Plotly
import space.kscience.plotly.makeFile
import space.kscience.plotly.*
import space.kscience.plotly.models.ScatterMode
import space.kscience.plotly.scatter
import space.kscience.plotly.models.appendXY
import kotlin.math.pow
fun main() {
val spectrum = SterileNeutrinoSpectrum()
val spectrum = SterileNeutrinoSpectrum().withNBkg()
val args: Map<Symbol, Double> = mapOf(
norm to 8e5,
bkg to 2.0,
mnu2 to 0.0,
e0 to 18575.0,
msterile2 to 1e6,
msterile2 to 1000.0.pow(2),
u2 to 1e-2,
thickness to 0.2
thickness to 0.1,
trap to 1.0
)
Plotly.plot {
Plotly.page {
plot {
scatter {
name = "Computed spectrum"
mode = ScatterMode.lines
x.buffer = 14000.0..18600.0 step 10.0
y.numbers = x.numbers.map { spectrum(it.toDouble(), args) }
y.numbers = x.doubles.map { spectrum(it, args) }
}
scatter {
name = "Old spectrum"
mode = ScatterMode.markers
javaClass.getResource("/old-spectrum.dat").readText().lines().map {
val (u, w) = it.split("\t")
appendXY(u.toDouble(), w.toDouble())
}
}
layout {
title = "Sterile neutrino spectrum"
}
}
plot {
val resolution = NumassResolution()
scatter {
name = "resolution"
x.buffer = 14000.0..14015.0 step 0.1
y.numbers = x.doubles.map { resolution(it, 14005.0, args) }
}
layout {
title = "Resolution, U = 14005.0"
}
}
plot {
val transmission = NumassTransmission()
scatter {
name = "transmission"
x.buffer = 14000.0..14100.0 step 0.2
y.numbers = x.doubles.map { transmission(it, 14005.0, args) }
}
layout {
title = "Resolution, U = 14005.0"
}
}
code {
+"""
norm to 8e5,
bkg to 2.0,
mnu2 to 0.0,
e0 to 18575.0,
msterile2 to 1000.0.pow(2),
u2 to 1e-2,
thickness to 0.1,
trap to 1.0
""".trimIndent()
}
}.makeFile()
println(spectrum(14000.0, args))
}

View File

@ -0,0 +1,461 @@
14000.0 21144.06826503362
14010.0 21006.572782462954
14020.0 20869.73215534316
14030.0 20733.405474277584
14040.0 20597.678076756514
14050.0 20462.59068330559
14060.0 20328.060309100107
14070.0 20194.08435615353
14080.0 20060.5857702678
14090.0 19927.897515335742
14100.0 19795.744180130758
14110.0 19664.120723613014
14120.0 19533.02778705386
14130.0 19402.632381360218
14140.0 19272.699152362227
14150.0 19143.347627762996
14160.0 19014.531805682913
14170.0 18886.408790202804
14180.0 18758.932989789675
14190.0 18631.846229720144
14200.0 18505.26601159422
14210.0 18379.392629602586
14220.0 18254.10223375718
14230.0 18129.363510394855
14240.0 18005.14584358382
14250.0 17881.62029183643
14260.0 17758.66138973419
14270.0 17635.98360171343
14280.0 17513.983186686626
14290.0 17392.578925822396
14300.0 17271.760548864386
14310.0 17151.482688239
14320.0 17031.779417997725
14330.0 16912.505887381285
14340.0 16793.96906902998
14350.0 16675.96921303285
14360.0 16558.41743144571
14370.0 16441.432931201125
14380.0 16324.953623636095
14390.0 16209.104211058508
14400.0 16093.73021604293
14410.0 15978.940653610127
14420.0 15864.736533819127
14430.0 15751.067774220024
14440.0 15637.866760895515
14450.0 15525.184953929003
14460.0 15413.09295555013
14470.0 15301.508741199292
14480.0 15190.496880751798
14490.0 15080.039018348876
14500.0 14970.124261418177
14510.0 14860.7425238067
14520.0 14751.695201569893
14530.0 14643.356106912928
14540.0 14535.488170477993
14550.0 14428.215148489058
14560.0 14321.411043220709
14570.0 14215.171022243501
14580.0 14109.36523470313
14590.0 14004.14388378637
14600.0 13899.542222301901
14610.0 13795.291516338752
14620.0 13691.556823690074
14630.0 13588.464699901404
14640.0 13485.832954556325
14650.0 13383.644657513441
14660.0 13282.035247171907
14670.0 13180.87709231583
14680.0 13080.311002851942
14690.0 12980.244679012389
14700.0 12880.747932860657
14710.0 12781.647374930337
14720.0 12682.989706972592
14730.0 12584.94741175745
14740.0 12487.397543672685
14750.0 12390.359015628186
14760.0 12293.780735136428
14770.0 12197.70496416098
14780.0 12102.166403721063
14790.0 12007.015454031693
14800.0 11912.479014703355
14810.0 11818.412654395914
14820.0 11724.755853920242
14830.0 11631.690434050815
14840.0 11539.09461085014
14850.0 11446.934809590526
14860.0 11355.29826945811
14870.0 11264.103206419837
14880.0 11173.391414671089
14890.0 11083.191157977213
14900.0 10993.542243060661
14910.0 10904.3437625975
14920.0 10815.57991970103
14930.0 10727.343250432023
14940.0 10639.51381282399
14950.0 10552.124239914427
14960.0 10465.302372109274
14970.0 10379.001774243998
14980.0 10293.08219246836
14990.0 10207.61593622077
15000.0 10122.661048674467
15010.0 10038.272286139583
15020.0 9954.240722704115
15030.0 9870.630464233276
15040.0 9787.596581873242
15050.0 9704.92878000812
15060.0 9622.769864774036
15070.0 9541.08694072584
15080.0 9459.82229559507
15090.0 9379.053465373388
15100.0 9298.772441033982
15110.0 9218.879315513399
15120.0 9139.468822719997
15130.0 9060.509843593685
15140.0 8982.055742086564
15150.0 8903.99444334693
15160.0 8826.317501558422
15170.0 8749.187182236928
15180.0 8672.453689772115
15190.0 8596.246591749144
15200.0 8520.427754118471
15210.0 8445.087555717617
15220.0 8370.11748545368
15230.0 8295.647637731907
15240.0 8221.547649216913
15250.0 8147.94083274912
15260.0 8074.7804536478025
15270.0 8002.059066003284
15280.0 7929.752465260882
15290.0 7857.8918938278075
15300.0 7786.432964345517
15310.0 7715.42864418354
15320.0 7644.86753765072
15330.0 7574.717041508431
15340.0 7504.95708969505
15350.0 7435.691961473507
15360.0 7366.800829538794
15370.0 7298.3870564224835
15380.0 7230.332137707168
15390.0 7162.762258781331
15400.0 7095.543205027884
15410.0 7028.816020258866
15420.0 6962.4441705414565
15430.0 6896.48494424743
15440.0 6830.93648543435
15450.0 6765.884954538622
15460.0 6701.179718966927
15470.0 6636.9212378276525
15480.0 6573.027429324024
15490.0 6509.553847185331
15500.0 6446.471412141977
15510.0 6383.843741433191
15520.0 6321.607038216404
15530.0 6259.747821228214
15540.0 6198.295937294384
15550.0 6137.242641018659
15560.0 6076.626172822268
15570.0 6016.366457084789
15580.0 5956.510619644465
15590.0 5897.062530707578
15600.0 5837.968553009778
15610.0 5779.315758519532
15620.0 5721.031250356002
15630.0 5663.187996894608
15640.0 5605.664107554607
15650.0 5548.540712406761
15660.0 5491.813140382814
15670.0 5435.489532444286
15680.0 5379.53988312815
15690.0 5323.957501315221
15700.0 5268.784529428275
15710.0 5213.949789397878
15720.0 5159.511366420885
15730.0 5105.412684951784
15740.0 5051.750839397989
15750.0 4998.457610290456
15760.0 4945.515385186621
15770.0 4892.961760871097
15780.0 4840.756414526476
15790.0 4788.944773061411
15800.0 4737.498617468057
15810.0 4686.3909452158
15820.0 4635.694502745357
15830.0 4585.326432683856
15840.0 4535.337006496239
15850.0 4485.73808232071
15860.0 4436.422695192254
15870.0 4387.535162584088
15880.0 4338.971648277426
15890.0 4290.778543741371
15900.0 4242.964756654218
15910.0 4195.485498285614
15920.0 4148.340792808015
15930.0 4101.574722284884
15940.0 4055.1298072593986
15950.0 4009.0525281291575
15960.0 3963.3736738704315
15970.0 3917.9677718938033
15980.0 3872.9398448260504
15990.0 3828.2164200868797
16000.0 3783.895173928012
16010.0 3739.886718632015
16020.0 3696.2261734572267
16030.0 3652.8814126004845
16040.0 3609.910341273943
16050.0 3567.2495172963345
16060.0 3524.9272329072337
16070.0 3482.955227529935
16080.0 3441.2881310607927
16090.0 3399.989789171778
16100.0 3358.990279803416
16110.0 3318.3354584720028
16120.0 3277.987694893713
16130.0 3237.99206058647
16140.0 3198.3194118107995
16150.0 3158.9648599415045
16160.0 3119.9291430710155
16170.0 3081.1934964628285
16180.0 3042.823324903815
16190.0 3004.7344280439656
16200.0 2966.9851740278386
16210.0 2929.5430695917976
16220.0 2892.3987181568477
16230.0 2855.6100275832555
16240.0 2819.085845132773
16250.0 2782.877962791771
16260.0 2747.0291572078468
16270.0 2711.4448475331446
16280.0 2676.1773173268366
16290.0 2641.2040797156797
16300.0 2606.549844158465
16310.0 2572.189980185393
16320.0 2538.1307408394996
16330.0 2504.376545819634
16340.0 2470.9082731773897
16350.0 2437.7523686766176
16360.0 2404.8884406197394
16370.0 2372.3267433090205
16380.0 2340.041679868535
16390.0 2308.0566681265714
16400.0 2276.364900483522
16410.0 2244.9608235478695
16420.0 2213.851181038246
16430.0 2183.036130062953
16440.0 2152.4898830631255
16450.0 2122.222946832527
16460.0 2092.2739581641
16470.0 2062.5670871504676
16480.0 2033.1582544737187
16490.0 2004.0295712007196
16500.0 1975.1870752561688
16510.0 1946.6174732521065
16520.0 1918.3181762061895
16530.0 1890.2873908718332
16540.0 1862.5411704472265
16550.0 1835.0590968312558
16560.0 1807.8538348100612
16570.0 1780.911853576155
16580.0 1754.2308050172048
16590.0 1727.8325389404906
16600.0 1701.698442698838
16610.0 1675.8157795027485
16620.0 1650.1758848467434
16630.0 1624.863638082829
16640.0 1599.7502526032201
16650.0 1574.9037795340814
16660.0 1550.3385739398202
16670.0 1526.0249879111723
16680.0 1501.9609190174958
16690.0 1478.1371551736843
16700.0 1454.562266355945
16710.0 1431.2629710458543
16720.0 1408.2020306676277
16730.0 1385.3967626619628
16740.0 1362.8228152351494
16750.0 1340.4990745445132
16760.0 1318.4219161315934
16770.0 1296.6010672641116
16780.0 1275.0088762832552
16790.0 1253.6453650584724
16800.0 1232.5442493232742
16810.0 1211.658289353798
16820.0 1191.0176647821295
16830.0 1170.6303036088184
16840.0 1150.4556908576478
16850.0 1130.5155785149711
16860.0 1110.8069491517203
16870.0 1091.3309731180668
16880.0 1072.0845257414933
16890.0 1053.0627415365889
16900.0 1034.2608134619222
16910.0 1015.6883432426487
16920.0 997.3451227461114
16930.0 979.2156990505905
16940.0 961.3166816443793
16950.0 943.6242449695036
16960.0 926.1527819466919
16970.0 908.9005704757275
16980.0 891.875990713501
16990.0 875.0380298220415
17000.0 858.4315340935602
17010.0 842.0345556237725
17020.0 825.839835611
17030.0 809.8541972689023
17040.0 794.0760340237487
17050.0 778.4974453963638
17060.0 763.1325187902158
17070.0 747.9695606512495
17080.0 733.0108899888257
17090.0 718.2378941774687
17100.0 703.6802707006024
17110.0 689.3200740592449
17120.0 675.1484180538769
17130.0 661.1772283643826
17140.0 647.4029267672108
17150.0 633.8133144520649
17160.0 620.4231905966275
17170.0 607.2137164548266
17180.0 594.2020648590451
17190.0 581.3783456474185
17200.0 568.7274210368327
17210.0 556.2733857970691
17220.0 544.0045977017528
17230.0 531.9117824415157
17240.0 519.9971020875178
17250.0 508.2678056184506
17260.0 496.71467779586277
17270.0 485.3427232164403
17280.0 474.1374682331266
17290.0 463.1120302522859
17300.0 452.26178665747074
17310.0 441.575171561372
17320.0 431.06841514926606
17330.0 420.7248016629219
17340.0 410.5455962743998
17350.0 400.5381056021479
17360.0 390.692606578034
17370.0 381.01256962550866
17380.0 371.49059203665655
17390.0 362.1385795412321
17400.0 352.9374563180437
17410.0 343.89363070121476
17420.0 335.0149122363285
17430.0 326.2810786848798
17440.0 317.7041270409196
17450.0 309.27957718766265
17460.0 301.0107264179578
17470.0 292.88782336642765
17480.0 284.9186052303832
17490.0 277.0936223153505
17500.0 269.41820031707056
17510.0 261.8883409551077
17520.0 254.50167234406157
17530.0 247.25811424081263
17540.0 240.15673078197244
17550.0 233.1964162491122
17560.0 226.37615287644778
17570.0 219.6924947596777
17580.0 213.14795710242234
17590.0 206.7333310548883
17600.0 200.45259075229302
17610.0 194.30102616977774
17620.0 188.27566991731828
17630.0 182.3796313106596
17640.0 176.60863758080657
17650.0 170.9619608607341
17660.0 165.43719009680015
17670.0 160.03533971145148
17680.0 154.75268002752233
17690.0 149.59013583552016
17700.0 144.5407096228448
17710.0 139.60946073903887
17720.0 134.79274745016315
17730.0 130.09050970754674
17740.0 125.49983265538104
17750.0 121.02105847914751
17760.0 116.65005220556819
17770.0 112.38900206712604
17780.0 108.23370450262375
17790.0 104.18525080537835
17800.0 100.24013083467196
17810.0 96.39772145483639
17820.0 92.65717123562793
17830.0 89.0162460869758
17840.0 85.47462006423214
17850.0 82.03090410010446
17860.0 78.68169683311878
17870.0 75.42879564467147
17880.0 72.26661207443473
17890.0 69.19740085849467
17900.0 66.2188166204438
17910.0 63.329476362960214
17920.0 60.52817559123287
17930.0 57.81382840461139
17940.0 55.18459677046772
17950.0 52.639207385766255
17960.0 50.17628063513036
17970.0 47.79485685506756
17980.0 45.49349261460641
17990.0 43.27070393645204
18000.0 41.124820497329566
18010.0 39.05416297660028
18020.0 37.05763103192207
18030.0 35.13403855083936
18040.0 33.28241970446307
18050.0 31.501328016008742
18060.0 29.789232968602683
18070.0 28.14480287113967
18080.0 26.566892964175935
18090.0 25.053956614987822
18100.0 23.604395867135885
18110.0 22.21670656050556
18120.0 20.889529413788292
18130.0 19.62201632664601
18140.0 18.412882964633525
18150.0 17.260422663377845
18160.0 16.16336469143388
18170.0 15.120474319163025
18180.0 14.130233108077414
18190.0 13.190878829930073
18200.0 12.301470618491864
18210.0 11.460817915538579
18220.0 10.667432266646426
18230.0 9.91985795035206
18240.0 9.21683299142454
18250.0 8.556612650396238
18260.0 7.938160755458398
18270.0 7.3601775328248
18280.0 6.8213004386287
18290.0 6.320297242431397
18300.0 5.855529641786286
18310.0 5.425186892341504
18320.0 5.028036943022659
18330.0 4.6631960393091365
18340.0 4.32930799884052
18350.0 4.0249367725373055
18360.0 3.7486423441053933
18370.0 3.499083287193633
18380.0 3.2748421288122014
18390.0 3.0745715651847014
18400.0 2.8964113339074795
18410.0 2.7395339740716547
18420.0 2.60254765486773
18430.0 2.4840384601692094
18440.0 2.38261999444518
18450.0 2.2966842827978713
18460.0 2.225103820486433
18470.0 2.166504319747381
18480.0 2.119417872453323
18490.0 2.082513315540327
18500.0 2.0544548933160414
18510.0 2.033862847336914
18520.0 2.019465429438452
18530.0 2.010044135140901
18540.0 2.004413085375449
18550.0 2.0014591003363216
18560.0 2.0002524725947946
18570.0 2.0000024904161795
18580.0 2.0
18590.0 2.0
18600.0 2.0