Model for Tristan experiment

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)]
-    ph.plot(data: (2000..19500).step(50).collectEntries { [it, model.value(it, params)] }, name: "spectrum")
+    }, 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());

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 {
+    private fun getShift(pars: ValueProvider): Double = pars.getDouble("shift", 0.0)
-        return pars.getDouble("shift", 0.0)
-    }
 
-    private fun getW(pars: ValueProvider): Double {
+    private fun getW(pars: ValueProvider): Double = pars.getDouble("w")
-        return pars.getDouble("w")
-    }
 
-    override fun providesDeriv(name: String): Boolean {
+    override fun providesDeriv(name: String): Boolean = true
-        return 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")
-    }
 }

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<Double, Double> {
+        val shift = getShift(params)
+        return Pair(shift - cutoff * getTailW(params), shift + cutoff * getW(params))
+    }
+
+    override fun getDerivSupport(parName: String, params: Values): Pair<Double, Double> {
+        val shift = getShift(params)
+        return Pair(shift - cutoff * getTailW(params), shift + cutoff * getW(params))
+    }
+}