From d9bb3ada0a17d7954a93fbd51ead292c316da386 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Wed, 29 Nov 2017 16:07:24 +0300 Subject: [PATCH] Model for Tristan experiment --- .../numass/scripts/models/TristanModel.groovy | 37 +++++--- .../kotlin/inr/numass/models/misc/Gauss.kt | 20 +---- .../kotlin/inr/numass/models/misc/ModGauss.kt | 84 +++++++++++++++++++ 3 files changed, 114 insertions(+), 27 deletions(-) create mode 100644 numass-main/src/main/kotlin/inr/numass/models/misc/ModGauss.kt diff --git a/numass-main/src/main/groovy/inr/numass/scripts/models/TristanModel.groovy b/numass-main/src/main/groovy/inr/numass/scripts/models/TristanModel.groovy index 0ec10081..350c822b 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/models/TristanModel.groovy +++ b/numass-main/src/main/groovy/inr/numass/scripts/models/TristanModel.groovy @@ -19,7 +19,7 @@ import inr.numass.data.SpectrumAdapter import inr.numass.data.SpectrumGenerator import inr.numass.models.NBkgSpectrum import inr.numass.models.NumassModelsKt -import inr.numass.models.misc.Gauss +import inr.numass.models.misc.ModGauss import inr.numass.models.sterile.NumassBeta import inr.numass.utils.DataModelUtils @@ -30,14 +30,16 @@ ctx.getPluginManager().load(PlotManager) ctx.getPluginManager().load(NumassPlugin) new GrindShell(ctx).eval { + PlotHelper ph = plots + def beta = new NumassBeta().getSpectrum(0) - def response = new Gauss(5.0) + def response = new ModGauss(5.0) ParametricFunction spectrum = NumassModelsKt.convolute(beta, response) def model = new XYModel(Meta.empty(), new SpectrumAdapter(Meta.empty()), new NBkgSpectrum(spectrum)); ParamSet params = morph(ParamSet, [:], "params") { - N(value: 1e+12, err: 30, lower: 0) + N(value: 1e+14, err: 30, lower: 0) bkg(value: 5.0, err: 0.1) E0(value: 18575.0, err: 0.1) mnu2(value: 0, err: 0.01) @@ -47,35 +49,48 @@ new GrindShell(ctx).eval { //trap(value: 1.0, err: 0.05) w(value: 150, err: 5) //shift(value: 1, err: 1e-2) - tail(value: 1e-4, err: 1e-5) + tailAmp(value: 0.01, err: 1e-2) + tailW(value: 300, err: 1) + } + + ph.plotFunction(-2000d, 500d, 400, "actual", "response") { double x -> + response.value(x, params) } SpectrumGenerator generator = new SpectrumGenerator(model, params, 12316); - PlotHelper ph = plots - - ph.plot(data: (2000..19500).step(50).collectEntries { [it, model.value(it, params)] }, name: "spectrum") + ph.plot(data: (2000..19500).step(50).collectEntries { + [it, model.value(it, params)] + }, name: "spectrum", frame: "test") .configure(showLine: true, showSymbol: false, showErrors: false, thickness: 2, connectionType: "spline", color: "red") Table data = generator.generateData(DataModelUtils.getUniformSpectrumConfiguration(10000, 19500, 1, 950)); - //params.setParValue("w", 151) + params.setParValue("w", 151) + //params.setParValue("tailAmp", 0.011) //params.setParValue("X", 0.01) //params.setParValue("trap", 0.01) //params.setParValue("mnu2", 4) - ph.plot(data: (2000..19500).step(50).collectEntries { [it, model.value(it, params)] }, name: "spectrum-mod") + + ph.plotFunction(-2000d, 500d, 400, "supposed", "response") { double x -> + response.value(x, params) + } + + ph.plot(data: (2000..19500).step(50).collectEntries { + [it, model.value(it, params)] + }, name: "spectrum-mod", frame: "test") .configure(showLine: true, showSymbol: false, showErrors: false, thickness: 2, connectionType: "spline", color: "green") - ph.plot(data: data, adapter: new SpectrumAdapter(Meta.empty())) + ph.plot(data: data, frame: "test", adapter: new SpectrumAdapter(Meta.empty())) .configure(color: "blue") FitState state = new FitState(data, model, params); def fm = ctx.getFeature(FitManager) - def res = fm.runStage(state, "QOW", FitStage.TASK_RUN, "N", "bkg", "E0", "U2"); + def res = fm.runStage(state, "MINUIT", FitStage.TASK_RUN, "N", "bkg", "E0", "U2"); res.printState(ctx.io.out().newPrintWriter()); diff --git a/numass-main/src/main/kotlin/inr/numass/models/misc/Gauss.kt b/numass-main/src/main/kotlin/inr/numass/models/misc/Gauss.kt index dddc1c6f..785f5f33 100644 --- a/numass-main/src/main/kotlin/inr/numass/models/misc/Gauss.kt +++ b/numass-main/src/main/kotlin/inr/numass/models/misc/Gauss.kt @@ -24,19 +24,13 @@ import java.lang.Math.* /** * @author Darksnake */ -class Gauss(private val cutoff: Double = 4.0) : AbstractParametricFunction(*list), FunctionSupport { +class Gauss(private val cutoff: Double = 4.0) : AbstractParametricFunction("w", "shift"), FunctionSupport { - private fun getShift(pars: ValueProvider): Double { - return pars.getDouble("shift", 0.0) - } + private fun getShift(pars: ValueProvider): Double = pars.getDouble("shift", 0.0) - private fun getW(pars: ValueProvider): Double { - return pars.getDouble("w") - } + private fun getW(pars: ValueProvider): Double = pars.getDouble("w") - override fun providesDeriv(name: String): Boolean { - return true - } + override fun providesDeriv(name: String): Boolean = true override fun value(d: Double, pars: Values): Double { @@ -72,10 +66,4 @@ class Gauss(private val cutoff: Double = 4.0) : AbstractParametricFunction(*list val w = getW(params) return Pair(shift - cutoff * w, shift + cutoff * w) } - - - companion object { - - private val list = arrayOf("w", "shift") - } } diff --git a/numass-main/src/main/kotlin/inr/numass/models/misc/ModGauss.kt b/numass-main/src/main/kotlin/inr/numass/models/misc/ModGauss.kt new file mode 100644 index 00000000..0714ebb7 --- /dev/null +++ b/numass-main/src/main/kotlin/inr/numass/models/misc/ModGauss.kt @@ -0,0 +1,84 @@ +/* + * Copyright 2015 Alexander Nozik. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package inr.numass.models.misc + +import hep.dataforge.stat.parametric.AbstractParametricFunction +import hep.dataforge.values.ValueProvider +import hep.dataforge.values.Values + +import java.lang.Math.* + +/** + * @author Darksnake + */ +class ModGauss(private val cutoff: Double = 4.0) : AbstractParametricFunction("w", "shift", "tailAmp", "tailW"), FunctionSupport { + + private fun getShift(pars: ValueProvider): Double = pars.getDouble("shift", 0.0) + + private fun getTailAmp(pars: ValueProvider): Double = pars.getDouble("tailAmp", 0.0) + + private fun getTailW(pars: ValueProvider): Double = pars.getDouble("tailW", 100.0) + + private fun getW(pars: ValueProvider): Double = pars.getDouble("w") + + override fun providesDeriv(name: String): Boolean = true + + + override fun value(d: Double, pars: Values): Double { + val shift = getShift(pars) + if (d - shift > cutoff * getW(pars)) { + return 0.0 + } + val aux = (d - shift) / getW(pars) + val tail = if (d > getShift(pars)) { + 0.0 + } else { + val tailW = getTailW(pars) + getTailAmp(pars) / tailW * Math.exp((d - shift) / tailW) + } + return exp(-aux * aux / 2) / getW(pars) / sqrt(2 * Math.PI) + tail + } + + override fun derivValue(parName: String, d: Double, pars: Values): Double { + if (abs(d - getShift(pars)) > cutoff * getW(pars)) { + return 0.0 + } + val pos = getShift(pars) + val w = getW(pars) + val tailW = getTailW(pars) + + return when (parName) { + "shift" -> this.value(d, pars) * (d - pos) / w / w + "w" -> this.value(d, pars) * ((d - pos) * (d - pos) / w / w / w - 1 / w) + "tailAmp" -> if (d > pos) { + 0.0 + } else { + Math.exp((d - pos) / tailW) / tailW + } + else -> return 0.0; + } + } + + override fun getSupport(params: Values): Pair { + val shift = getShift(params) + return Pair(shift - cutoff * getTailW(params), shift + cutoff * getW(params)) + } + + override fun getDerivSupport(parName: String, params: Values): Pair { + val shift = getShift(params) + return Pair(shift - cutoff * getTailW(params), shift + cutoff * getW(params)) + } +}