Restored kodex. A lot of minor fixes.

This commit is contained in:
Alexander Nozik 2018-01-10 17:01:33 +03:00
parent d5ea5190fa
commit 01c910682e
14 changed files with 178 additions and 133 deletions

View File

@ -28,7 +28,6 @@ import inr.numass.models.BetaSpectrum
import inr.numass.models.ModularSpectrum import inr.numass.models.ModularSpectrum
import inr.numass.models.NBkgSpectrum import inr.numass.models.NBkgSpectrum
import static hep.dataforge.maths.RandomUtils.setSeed
import static inr.numass.utils.DataModelUtils.getUniformSpectrumConfiguration import static inr.numass.utils.DataModelUtils.getUniformSpectrumConfiguration
PrintWriter out = Global.out(); PrintWriter out = Global.out();

View File

@ -36,26 +36,32 @@ new GrindShell(ctx).eval {
def response = new ModGauss(5.0) def response = new ModGauss(5.0)
ParametricFunction spectrum = NumassModelsKt.convolute(beta, response) ParametricFunction spectrum = NumassModelsKt.convolute(beta, response)
def model = new XYModel(Meta.empty(), new SpectrumAdapter(Meta.empty()), new NBkgSpectrum(spectrum)); def model = new XYModel(Meta.empty(), new SpectrumAdapter(Meta.empty()), new NBkgSpectrum(spectrum))
ParamSet params = morph(ParamSet, [:], "params") { ParamSet params = morph(ParamSet, [:], "params") {
N(value: 1e+14, err: 30, lower: 0) N(value: 1e+14, err: 30, lower: 0)
bkg(value: 5.0, err: 0.1) bkg(value: 5.0, err: 0.1)
E0(value: 18575.0, err: 0.1) E0(value: 18575.0, err: 0.1)
mnu2(value: 0, err: 0.01) mnu2(value: 0, err: 0.01)
msterile2(value: 1000**2, err: 1) msterile2(value: 7000**2, err: 1)
U2(value: 0.0, err: 1e-3) U2(value: 0.0, err: 1e-3)
//X(value: 0.0, err: 0.01, lower: 0) //X(value: 0.0, err: 0.01, lower: 0)
//trap(value: 1.0, err: 0.05) //trap(value: 1.0, err: 0.05)
w(value: 150, err: 5) w(value: 150, err: 5)
//shift(value: 1, err: 1e-2) //shift(value: 1, err: 1e-2)
tailAmp(value: 0.01, err: 1e-2) //tailAmp(value: 0.01, err: 1e-2)
tailW(value: 300, err: 1) tailW(value: 300, err: 1)
} }
ph.plotFunction(-2000d, 500d, 400, "actual", "response") { double x -> // double norm = NumassIntegrator.defaultIntegrator.integrate(1000d, 18500d) {
response.value(x, params) // model.value(it, params)
} // }
// println("The total number of events is $norm")
//
// ph.plotFunction(-2000d, 500d, 400, "actual", "response") { double x ->
// response.value(x, params)
// }
SpectrumGenerator generator = new SpectrumGenerator(model, params, 12316); SpectrumGenerator generator = new SpectrumGenerator(model, params, 12316);
@ -65,7 +71,7 @@ new GrindShell(ctx).eval {
.configure(showLine: true, showSymbol: false, showErrors: false, thickness: 2, connectionType: "spline", color: "red") .configure(showLine: true, showSymbol: false, showErrors: false, thickness: 2, connectionType: "spline", color: "red")
Table data = generator.generateData(DataModelUtils.getUniformSpectrumConfiguration(10000, 19500, 1, 950)); Table data = generator.generateData(DataModelUtils.getUniformSpectrumConfiguration(7000, 19500, 1, 1000));
//params.setParValue("w", 151) //params.setParValue("w", 151)
//params.setParValue("tailAmp", 0.011) //params.setParValue("tailAmp", 0.011)

View File

@ -48,7 +48,7 @@ public class Numass {
.ln() .ln()
.text("\t") .text("\t")
.content( .content(
MarkupUtils.markupDescriptor(DescriptorUtils.buildDescriptor("method::hep.dataforge.data.DataManager.read")) MarkupUtils.INSTANCE.markupDescriptor(DescriptorUtils.buildDescriptor("method::hep.dataforge.data.DataManager.read"))
) )
.ln() .ln()
.text("***Allowed actions***", "red") .text("***Allowed actions***", "red")
@ -60,7 +60,7 @@ public class Numass {
am.getAllActions() am.getAllActions()
.map(name -> am.optAction(name).get()) .map(name -> am.optAction(name).get())
.map(ActionDescriptor::build).forEach(descriptor -> .map(ActionDescriptor::build).forEach(descriptor ->
builder.text("\t").content(MarkupUtils.markupDescriptor(descriptor)) builder.text("\t").content(MarkupUtils.INSTANCE.markupDescriptor(descriptor))
); );
builder.text("***End of actions list***", "red"); builder.text("***End of actions list***", "red");

View File

@ -16,6 +16,7 @@
package inr.numass.data; package inr.numass.data;
import hep.dataforge.meta.Meta; import hep.dataforge.meta.Meta;
import hep.dataforge.stat.RandomKt;
import hep.dataforge.stat.fit.ParamSet; import hep.dataforge.stat.fit.ParamSet;
import hep.dataforge.stat.models.Generator; import hep.dataforge.stat.models.Generator;
import hep.dataforge.stat.models.XYModel; import hep.dataforge.stat.models.XYModel;
@ -27,7 +28,6 @@ import org.apache.commons.math3.random.JDKRandomGenerator;
import org.apache.commons.math3.random.RandomDataGenerator; import org.apache.commons.math3.random.RandomDataGenerator;
import org.apache.commons.math3.random.RandomGenerator; import org.apache.commons.math3.random.RandomGenerator;
import static hep.dataforge.maths.RandomUtils.getDefaultRandomGenerator;
import static java.lang.Double.isNaN; import static java.lang.Double.isNaN;
import static java.lang.Math.sqrt; import static java.lang.Math.sqrt;
@ -62,7 +62,7 @@ public class SpectrumGenerator implements Generator {
} }
public SpectrumGenerator(XYModel source, ParamSet params) { public SpectrumGenerator(XYModel source, ParamSet params) {
this(source, params, INSTANCE.getDefaultRandomGenerator()); this(source, params, RandomKt.getDefaultGenerator());
} }
@Override @Override

View File

@ -180,12 +180,12 @@ public class LossCalculator {
final LossCalculator loss = LossCalculator.instance; final LossCalculator loss = LossCalculator.instance;
final List<Double> probs = loss.getGunLossProbabilities(set.getDouble("X")); final List<Double> probs = loss.getGunLossProbabilities(set.getDouble("X"));
UnivariateFunction single = (double e) -> probs.get(1) * scatterFunction.value(e); UnivariateFunction single = (double e) -> probs.get(1) * scatterFunction.value(e);
frame.add(XYFunctionPlot.plotFunction("Single scattering", single::value, 0, 100, 1000)); frame.add(XYFunctionPlot.plotFunction("Single scattering", 0, 100, 1000, single::value));
for (int i = 2; i < probs.size(); i++) { for (int i = 2; i < probs.size(); i++) {
final int j = i; final int j = i;
UnivariateFunction scatter = (double e) -> probs.get(j) * loss.getLossValue(j, e, 0d); UnivariateFunction scatter = (double e) -> probs.get(j) * loss.getLossValue(j, e, 0d);
frame.add(XYFunctionPlot.plotFunction(j + " scattering", scatter::value, 0, 100, 1000)); frame.add(XYFunctionPlot.plotFunction(j + " scattering", 0, 100, 1000, scatter::value));
} }
UnivariateFunction total = (eps) -> { UnivariateFunction total = (eps) -> {
@ -199,11 +199,11 @@ public class LossCalculator {
return sum; return sum;
}; };
frame.add(XYFunctionPlot.plotFunction("Total loss", total::value, 0, 100, 1000)); frame.add(XYFunctionPlot.plotFunction("Total loss", 0, 100, 1000, total::value));
} else { } else {
frame.add(XYFunctionPlot.plotFunction("Differential crosssection", scatterFunction::value, 0, 100, 2000)); frame.add(XYFunctionPlot.plotFunction("Differential crosssection", 0, 100, 2000, scatterFunction::value));
} }
} }

View File

@ -3,13 +3,13 @@ package inr.numass.data
import hep.dataforge.maths.chain.Chain import hep.dataforge.maths.chain.Chain
import hep.dataforge.maths.chain.MarkovChain import hep.dataforge.maths.chain.MarkovChain
import hep.dataforge.maths.chain.StatefulChain import hep.dataforge.maths.chain.StatefulChain
import hep.dataforge.stat.defaultGenerator
import inr.numass.data.api.NumassBlock import inr.numass.data.api.NumassBlock
import inr.numass.data.api.NumassEvent import inr.numass.data.api.NumassEvent
import inr.numass.data.api.SimpleBlock import inr.numass.data.api.SimpleBlock
import kotlinx.coroutines.experimental.channels.takeWhile import kotlinx.coroutines.experimental.channels.takeWhile
import kotlinx.coroutines.experimental.channels.toList import kotlinx.coroutines.experimental.channels.toList
import kotlinx.coroutines.experimental.runBlocking import kotlinx.coroutines.experimental.runBlocking
import org.apache.commons.math3.random.JDKRandomGenerator
import org.apache.commons.math3.random.RandomGenerator import org.apache.commons.math3.random.RandomGenerator
import java.time.Duration import java.time.Duration
import java.time.Instant import java.time.Instant
@ -24,12 +24,10 @@ private fun RandomGenerator.nextDeltaTime(cr: Double): Long {
fun generateBlock(start: Instant, length: Long, chain: Chain<NumassEvent>): NumassBlock { fun generateBlock(start: Instant, length: Long, chain: Chain<NumassEvent>): NumassBlock {
val events = runBlocking { chain.asChannel().takeWhile { it.timeOffset < length }.toList()} val events = runBlocking { chain.channel.takeWhile { it.timeOffset < length }.toList()}
return SimpleBlock(start, Duration.ofNanos(length), events) return SimpleBlock(start, Duration.ofNanos(length), events)
} }
internal val defaultGenerator = JDKRandomGenerator()
internal val defaultAmplitudeGenerator: RandomGenerator.(NumassEvent?, Long) -> Short = { _, _ -> ((nextDouble() + 2.0) * 100).toShort() } internal val defaultAmplitudeGenerator: RandomGenerator.(NumassEvent?, Long) -> Short = { _, _ -> ((nextDouble() + 2.0) * 100).toShort() }
fun buildSimpleEventChain( fun buildSimpleEventChain(

View File

@ -0,0 +1,55 @@
package inr.numass.models.mc
import hep.dataforge.fx.plots.PlotManager
import hep.dataforge.kodex.buildMeta
import hep.dataforge.kodex.global
import hep.dataforge.maths.chain.Chain
import hep.dataforge.plots.data.XYFunctionPlot
import hep.dataforge.stat.PolynomialDistribution
import hep.dataforge.stat.fit.ParamSet
import inr.numass.NumassPlugin
import inr.numass.models.sterile.SterileNeutrinoSpectrum
fun sampleBeta(params: ParamSet): Chain<Double> {
TODO()
}
fun main(args: Array<String>) {
NumassPlugin().startGlobal()
val pm = PlotManager().apply { startGlobal() }
val meta = buildMeta("model") {
"fast" to true
node("resolution") {
"width" to 1.22e-4
"tailAlpha" to 1e-2
}
}
val allPars = ParamSet()
.setPar("N", 7e+05, 1.8e+03, 0.0, java.lang.Double.POSITIVE_INFINITY)
.setPar("bkg", 1.0, 0.050)
.setPar("E0", 18575.0, 1.4)
.setPar("mnu2", 0.0, 1.0)
.setPar("msterile2", 1000.0 * 1000.0, 0.0)
.setPar("U2", 0.0, 1e-4, -1.0, 1.0)
.setPar("X", 0.0, 0.01, 0.0, java.lang.Double.POSITIVE_INFINITY)
.setPar("trap", 1.0, 0.01, 0.0, java.lang.Double.POSITIVE_INFINITY)
val sp = SterileNeutrinoSpectrum(global, meta)
val spectrumPlot = XYFunctionPlot.plotFunction("spectrum", 14000.0, 18600.0, 500) {
sp.value(it, allPars)
}
val distribution = PolynomialDistribution(0.0, 5000.0, 3.0);
val distributionPlot = XYFunctionPlot.plotFunction("distribution", 14000.0, 18500.0, 500) {
50*distribution.density(18575.0 - it)
}
pm.getPlotFrame("beta").apply {
add(spectrumPlot)
add(distributionPlot)
}
}

View File

@ -1,6 +1,6 @@
package inr.numass.models.misc package inr.numass.models.misc
import hep.dataforge.maths.HyperSquareDomain import hep.dataforge.maths.domains.HyperSquareDomain
import hep.dataforge.values.Values import hep.dataforge.values.Values
interface FunctionSupport { interface FunctionSupport {

View File

@ -1,104 +0,0 @@
/*
* To change this license header, choose License Headers in Project Properties.
* To change this template file, choose Tools | Templates
* and open the template in the editor.
*/
package inr.numass.models.sterile
import hep.dataforge.context.Context
import hep.dataforge.maths.MathPlugin
import hep.dataforge.meta.Meta
import hep.dataforge.meta.MetaBuilder
import hep.dataforge.stat.fit.ParamSet
import hep.dataforge.stat.parametric.ParametricFunction
import inr.numass.Numass
import inr.numass.models.BetaSpectrum
import inr.numass.models.ModularSpectrum
import inr.numass.models.NBkgSpectrum
import inr.numass.models.ResolutionFunction
import org.apache.commons.math3.analysis.BivariateFunction
/**
*
* @author Alexander Nozik
*/
object TestModels {
/**
* @param args the command line arguments
*/
@JvmStatic
fun main(args: Array<String>) {
val context = Numass.buildContext()
/*
<model modelName="sterile" fssFile="FS.txt">
<resolution width = "1.22e-4" tailAlpha="1e-2"/>
<transmission trapping="numass.trap2016"/>
</model>
*/
val meta = MetaBuilder("model")
.putNode(MetaBuilder("resolution")
.putValue("width", 1.22e-4)
.putValue("tailAlpha", 1e-2)
)
.putNode(MetaBuilder("transmission")
.putValue("trapping", "numass.trap2016")
)
val oldFunc = oldModel(context, meta)
val newFunc = newModel(context, meta)
val allPars = ParamSet()
.setPar("N", 7e+05, 1.8e+03, 0.0, java.lang.Double.POSITIVE_INFINITY)
.setPar("bkg", 1.0, 0.050)
.setPar("E0", 18575.0, 1.4)
.setPar("mnu2", 0.0, 1.0)
.setPar("msterile2", 1000.0 * 1000.0, 0.0)
.setPar("U2", 0.0, 1e-4, -1.0, 1.0)
.setPar("X", 0.04, 0.01, 0.0, java.lang.Double.POSITIVE_INFINITY)
.setPar("trap", 1.0, 0.01, 0.0, java.lang.Double.POSITIVE_INFINITY)
var u = 14000.0
while (u < 18600) {
// double oldVal = oldFunc.value(u, allPars);
// double newVal = newFunc.value(u, allPars);
val oldVal = oldFunc.derivValue("trap", u, allPars)
val newVal = newFunc.derivValue("trap", u, allPars)
System.out.printf("%f\t%g\t%g\t%g%n", u, oldVal, newVal, 1.0 - oldVal / newVal)
u += 100.0
}
}
private fun oldModel(context: Context, meta: Meta): ParametricFunction {
val A = meta.getDouble("resolution", meta.getDouble("resolution.width", 8.3e-5)!!)!!//8.3e-5
val from = meta.getDouble("from", 13900.0)!!
val to = meta.getDouble("to", 18700.0)!!
context.chronicle.report("Setting up tritium model with real transmission function")
val resolutionTail: BivariateFunction
if (meta.hasValue("resolution.tailAlpha")) {
resolutionTail = ResolutionFunction.getAngledTail(meta.getDouble("resolution.tailAlpha")!!, meta.getDouble("resolution.tailBeta", 0.0)!!)
} else {
resolutionTail = ResolutionFunction.getRealTail()
}
//RangedNamedSetSpectrum beta = new BetaSpectrum(context.io().getFile("FS.txt"));
val beta = BetaSpectrum()
val sp = ModularSpectrum(beta, ResolutionFunction(A, resolutionTail), from, to)
if (meta.getBoolean("caching", false)!!) {
context.chronicle.report("Caching turned on")
sp.setCaching(true)
}
//Adding trapping energy dependence
if (meta.hasValue("transmission.trapping")) {
val trap = MathPlugin.buildFrom(context).buildBivariateFunction(meta.getString("transmission.trapping"))
sp.setTrappingFunction(trap)
}
return NBkgSpectrum(sp)
}
private fun newModel(context: Context, meta: Meta): ParametricFunction {
val sp = SterileNeutrinoSpectrum(context, meta)
return NBkgSpectrum(sp)
}
}

View File

@ -1,4 +1,4 @@
package inr.numass.scripts package inr.numass.scripts.utils
import hep.dataforge.io.XMLMetaWriter import hep.dataforge.io.XMLMetaWriter
import hep.dataforge.kodex.buildMeta import hep.dataforge.kodex.buildMeta
@ -74,7 +74,7 @@ fun main(args: Array<String>) {
val directory = if (args.isNotEmpty()) { val directory = if (args.isNotEmpty()) {
args.first() args.first()
} else { } else {
"." ""
} }
val path = Paths.get(directory) val path = Paths.get(directory)

View File

@ -13,13 +13,11 @@
* See the License for the specific language governing permissions and * See the License for the specific language governing permissions and
* limitations under the License. * limitations under the License.
*/ */
package inr.numass.run; package inr.numass.models;
import hep.dataforge.exceptions.NamingException; import hep.dataforge.exceptions.NamingException;
import hep.dataforge.stat.fit.MINUITPlugin; import hep.dataforge.stat.fit.MINUITPlugin;
import hep.dataforge.stat.fit.ParamSet; import hep.dataforge.stat.fit.ParamSet;
import inr.numass.models.BetaSpectrum;
import inr.numass.models.ModularSpectrum;
import java.io.File; import java.io.File;
import java.io.FileNotFoundException; import java.io.FileNotFoundException;

View File

@ -0,0 +1,93 @@
/*
* To change this license header, choose License Headers in Project Properties.
* To change this template file, choose Tools | Templates
* and open the template in the editor.
*/
package inr.numass.models
import hep.dataforge.context.Context
import hep.dataforge.kodex.step
import hep.dataforge.maths.MathPlugin
import hep.dataforge.meta.Meta
import hep.dataforge.meta.MetaBuilder
import hep.dataforge.stat.fit.ParamSet
import hep.dataforge.stat.parametric.ParametricFunction
import inr.numass.Numass
import inr.numass.models.sterile.SterileNeutrinoSpectrum
import org.apache.commons.math3.analysis.BivariateFunction
/**
* @param args the command line arguments
*/
fun main(args: Array<String>) {
val context = Numass.buildContext()
/*
<model modelName="sterile" fssFile="FS.txt">
<resolution width = "1.22e-4" tailAlpha="1e-2"/>
<transmission trapping="numass.trap2016"/>
</model>
*/
val meta = MetaBuilder("model")
.putNode(MetaBuilder("resolution")
.putValue("width", 1.22e-4)
.putValue("tailAlpha", 1e-2)
)
.putNode(MetaBuilder("transmission")
// .putValue("trapping", "function::numass.trap.nominal")
)
val oldFunc = oldModel(context, meta)
val newFunc = newModel(context, meta)
val allPars = ParamSet()
.setPar("N", 7e+05, 1.8e+03, 0.0, java.lang.Double.POSITIVE_INFINITY)
.setPar("bkg", 1.0, 0.050)
.setPar("E0", 18575.0, 1.4)
.setPar("mnu2", 0.0, 1.0)
.setPar("msterile2", 1000.0 * 1000.0, 0.0)
.setPar("U2", 0.0, 1e-4, -1.0, 1.0)
.setPar("X", 0.0, 0.01, 0.0, java.lang.Double.POSITIVE_INFINITY)
.setPar("trap", 1.0, 0.01, 0.0, java.lang.Double.POSITIVE_INFINITY)
for (u in 14000.0..18600.0 step 100.0) {
val oldVal = oldFunc.value(u, allPars);
val newVal = newFunc.value(u, allPars);
// val oldVal = oldFunc.derivValue("trap", u, allPars)
// val newVal = newFunc.derivValue("trap", u, allPars)
System.out.printf("%f\t%g\t%g\t%g%n", u, oldVal, newVal, 1.0 - oldVal / newVal)
}
}
private fun oldModel(context: Context, meta: Meta): ParametricFunction {
val A = meta.getDouble("resolution", meta.getDouble("resolution.width", 8.3e-5))//8.3e-5
val from = meta.getDouble("from", 13900.0)
val to = meta.getDouble("to", 18700.0)
context.chronicle.report("Setting up tritium model with real transmission function")
val resolutionTail: BivariateFunction = if (meta.hasValue("resolution.tailAlpha")) {
ResolutionFunction.getAngledTail(meta.getDouble("resolution.tailAlpha"), meta.getDouble("resolution.tailBeta", 0.0))
} else {
ResolutionFunction.getRealTail()
}
//RangedNamedSetSpectrum beta = new BetaSpectrum(context.io().getFile("FS.txt"));
val sp = ModularSpectrum(BetaSpectrum(), ResolutionFunction(A, resolutionTail), from, to)
if (meta.getBoolean("caching", false)) {
context.chronicle.report("Caching turned on")
sp.setCaching(true)
}
//Adding trapping energy dependence
if (meta.hasValue("transmission.trapping")) {
val trap = MathPlugin.buildFrom(context).buildBivariateFunction(meta.getString("transmission.trapping"))
sp.setTrappingFunction(trap)
}
return NBkgSpectrum(sp)
}
private fun newModel(context: Context, meta: Meta): ParametricFunction {
val sp = SterileNeutrinoSpectrum(context, meta)
return NBkgSpectrum(sp)
}

View File

@ -41,8 +41,8 @@ public class TestNeLossParametrisation {
System.out.println(norm); System.out.println(norm);
frame.add(XYFunctionPlot.plotFunction("old", oldFunction::value, 0, 30, 300)); frame.add(XYFunctionPlot.plotFunction("old", 0, 30, 300, oldFunction::value));
frame.add(XYFunctionPlot.plotFunction("new", newFunction::value, 0, 30, 300)); frame.add(XYFunctionPlot.plotFunction("new", 0, 30, 300, newFunction::value));
} }
public static UnivariateFunction getSingleScatterFunction( public static UnivariateFunction getSingleScatterFunction(

View File

@ -33,7 +33,7 @@ public class TransmissionInterpolatorTest {
TransmissionInterpolator interpolator = TransmissionInterpolator.fromFile(Global.instance(), TransmissionInterpolator interpolator = TransmissionInterpolator.fromFile(Global.instance(),
"d:\\sterile-new\\loss2014-11\\.dataforge\\merge\\empty_sum.onComplete", "Uset", "CR", 15, 0.8, 19002d); "d:\\sterile-new\\loss2014-11\\.dataforge\\merge\\empty_sum.onComplete", "Uset", "CR", 15, 0.8, 19002d);
frame.add(DataPlot.plot("data", interpolator.getX(), interpolator.getY())); frame.add(DataPlot.plot("data", interpolator.getX(), interpolator.getY()));
frame.add(XYFunctionPlot.plotFunction("interpolated", interpolator::value, interpolator.getXmin(), interpolator.getXmax(), 2000)); frame.add(XYFunctionPlot.plotFunction("interpolated", interpolator.getXmin(), interpolator.getXmax(), 2000, interpolator::value));
// PrintFunction.printFuntionSimple(new PrintWriter(System.onComplete), interpolator, interpolator.getXmin(), interpolator.getXmax(), 500); // PrintFunction.printFuntionSimple(new PrintWriter(System.onComplete), interpolator, interpolator.getXmin(), interpolator.getXmax(), 500);
} }