diff --git a/numass-main/src/main/groovy/inr/numass/scripts/SimulateGun.groovy b/numass-main/src/main/groovy/inr/numass/scripts/SimulateGun.groovy index 3f43363c..e494b2a3 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/SimulateGun.groovy +++ b/numass-main/src/main/groovy/inr/numass/scripts/SimulateGun.groovy @@ -13,59 +13,59 @@ * See the License for the specific language governing permissions and * limitations under the License. */ -package inr.numass.scripts; - -import hep.dataforge.context.GlobalContext; -import hep.dataforge.datafitter.FitManager; -import hep.dataforge.datafitter.ParamSet; -import hep.dataforge.datafitter.models.XYModel; -import hep.dataforge.exceptions.NamingException; -import hep.dataforge.io.PrintNamed; -import inr.numass.data.SpectrumDataAdapter; -import inr.numass.models.GunSpectrum; -import inr.numass.models.NBkgSpectrum; -import java.io.FileNotFoundException; -import java.io.PrintWriter; -import java.util.Locale; -import static java.util.Locale.setDefault; - - -setDefault(Locale.US); -GlobalContext global = GlobalContext.instance(); -// global.loadModule(new MINUITModule()); - -FitManager fm = new FitManager(); - -GunSpectrum gsp = new GunSpectrum(); -NBkgSpectrum spectrum = new NBkgSpectrum(gsp); - -XYModel model = new XYModel("gun", spectrum, new SpectrumDataAdapter()); - -ParamSet allPars = new ParamSet() -.setPar("N", 1e3, 1e2) -.setPar("pos", 18500, 0.1) -.setPar("bkg", 50, 1) -.setPar("resA", 5.3e-5, 1e-5) -.setPar("sigma", 0.3, 0.03); - -PrintNamed.printSpectrum(new PrintWriter(System.out), spectrum, allPars, 18495, 18505, 100); - -allPars.setParValue("sigma", 0.6); - -PrintNamed.printSpectrum(new PrintWriter(System.out), spectrum, allPars, 18495, 18505, 100); - -// //String fileName = "d:\\PlayGround\\merge\\scans.out"; -//// String configName = "d:\\PlayGround\\SCAN.CFG"; -//// ListTable config = OldDataReader.readConfig(configName); -// SpectrumGenerator generator = new SpectrumGenerator(model, allPars, 12316); -// -// ListTable data = generator.generateData(DataModelUtils.getUniformSpectrumConfiguration(18495, 18505, 20, 20)); -// -//// data = data.filter("X", Value.of(15510.0), Value.of(18610.0)); -//// allPars.setParValue("X", 0.4); -// FitState state = FitTaskManager.buildState(data, model, allPars); -// -// FitState res = fm.runTask(state, "QOW", FitTask.TASK_RUN, "N", "bkg", "pos", "sigma"); -// -// res.print(out()); - +package inr.numass.scripts; + +import hep.dataforge.context.GlobalContext; +import hep.dataforge.datafitter.FitManager; +import hep.dataforge.datafitter.ParamSet; +import hep.dataforge.datafitter.models.XYModel; +import hep.dataforge.exceptions.NamingException; +import hep.dataforge.io.FittingIOUtils; +import inr.numass.data.SpectrumDataAdapter; +import inr.numass.models.GunSpectrum; +import inr.numass.models.NBkgSpectrum; +import java.io.FileNotFoundException; +import java.io.PrintWriter; +import java.util.Locale; +import static java.util.Locale.setDefault; + + +setDefault(Locale.US); +GlobalContext global = GlobalContext.instance(); +// global.loadModule(new MINUITModule()); + +FitManager fm = new FitManager(); + +GunSpectrum gsp = new GunSpectrum(); +NBkgSpectrum spectrum = new NBkgSpectrum(gsp); + +XYModel model = new XYModel("gun", spectrum, new SpectrumDataAdapter()); + +ParamSet allPars = new ParamSet() +.setPar("N", 1e3, 1e2) +.setPar("pos", 18500, 0.1) +.setPar("bkg", 50, 1) +.setPar("resA", 5.3e-5, 1e-5) +.setPar("sigma", 0.3, 0.03); + +PrintNamed.printSpectrum(new PrintWriter(System.out), spectrum, allPars, 18495, 18505, 100); + +allPars.setParValue("sigma", 0.6); + +FittingIOUtils.printSpectrum(new PrintWriter(System.out), spectrum, allPars, 18495, 18505, 100); + +// //String fileName = "d:\\PlayGround\\merge\\scans.out"; +//// String configName = "d:\\PlayGround\\SCAN.CFG"; +//// ListTable config = OldDataReader.readConfig(configName); +// SpectrumGenerator generator = new SpectrumGenerator(model, allPars, 12316); +// +// ListTable data = generator.generateData(DataModelUtils.getUniformSpectrumConfiguration(18495, 18505, 20, 20)); +// +//// data = data.filter("X", Value.of(15510.0), Value.of(18610.0)); +//// allPars.setParValue("X", 0.4); +// FitState state = FitTaskManager.buildState(data, model, allPars); +// +// FitState res = fm.runTask(state, "QOW", FitTask.TASK_RUN, "N", "bkg", "pos", "sigma"); +// +// res.print(out()); + diff --git a/numass-main/src/main/groovy/inr/numass/scripts/SterileDemo.groovy b/numass-main/src/main/groovy/inr/numass/scripts/SterileDemo.groovy new file mode 100644 index 00000000..fe2464f6 --- /dev/null +++ b/numass-main/src/main/groovy/inr/numass/scripts/SterileDemo.groovy @@ -0,0 +1,86 @@ +/* + * 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.scripts; + +import hep.dataforge.context.GlobalContext; +import static hep.dataforge.context.GlobalContext.out; +import hep.dataforge.tables.ListTable; +import hep.dataforge.datafitter.FitManager; +import hep.dataforge.datafitter.FitState; +import hep.dataforge.datafitter.FitTask; +import hep.dataforge.datafitter.MINUITPlugin + +import hep.dataforge.datafitter.ParamSet; +import hep.dataforge.datafitter.models.XYModel; +import hep.dataforge.exceptions.NamingException; +import hep.dataforge.exceptions.PackFormatException; +import inr.numass.data.SpectrumDataAdapter; +import inr.numass.data.SpectrumGenerator; +import inr.numass.models.BetaSpectrum +import inr.numass.models.ModularSpectrum; +import inr.numass.models.NBkgSpectrum; +import inr.numass.utils.DataModelUtils; +import hep.dataforge.plotfit.PlotFitResultAction; +import java.io.FileNotFoundException; +import java.util.Locale; +import static java.util.Locale.setDefault; +import inr.numass.utils.TritiumUtils; +import inr.numass.data.SpectrumDataAdapter; +import hep.dataforge.io.FittingIOUtils + +/** + * + * @author Darksnake + */ + +setDefault(Locale.US); + +ModularSpectrum beta = new ModularSpectrum(new BetaSpectrum(), 8.3e-5, 13990d, 18600d); +beta.setCaching(false); + +NBkgSpectrum spectrum = new NBkgSpectrum(beta); +XYModel model = new XYModel(spectrum, new SpectrumDataAdapter()); + +ParamSet allPars = new ParamSet(); + +allPars.setPar("N", 6.6579e+05, 1.8e+03, 0d, Double.POSITIVE_INFINITY); +allPars.setPar("bkg", 0.5387, 0.050); +allPars.setPar("E0", 18574.94, 1.4); +allPars.setPar("mnu2", 0d, 1d); +allPars.setPar("msterile2", 1000d * 1000d,0); +allPars.setPar("U2", 0.0, 1e-4, -1d, 1d); +allPars.setPar("X", 0.04000, 0.01, 0d, Double.POSITIVE_INFINITY); +allPars.setPar("trap", 1.634, 0.01,0d, Double.POSITIVE_INFINITY); + +FittingIOUtils.printSpectrum(GlobalContext.out(), spectrum, allPars, 14000.0, 18600.0, 600); + +//SpectrumGenerator generator = new SpectrumGenerator(model, allPars, 12316); +// +//ListTable data = generator.generateData(DataModelUtils.getUniformSpectrumConfiguration(14000d, 18500, 2000, 90)); +// +//data = TritiumUtils.correctForDeadTime(data, new SpectrumDataAdapter(), 1e-8); +//// data = data.filter("X", Value.of(15510.0), Value.of(18610.0)); +//// allPars.setParValue("X", 0.4); +//FitState state = new FitState(data, model, allPars); +////new PlotFitResultAction().eval(state); +// +// +//FitState res = fm.runTask(state, "QOW", FitTask.TASK_RUN, "N", "bkg", "E0", "U2", "trap"); +// +// +// +//res.print(out()); +// diff --git a/numass-main/src/main/java/inr/numass/NumassPlugin.java b/numass-main/src/main/java/inr/numass/NumassPlugin.java index 39028153..8b594dc9 100644 --- a/numass-main/src/main/java/inr/numass/NumassPlugin.java +++ b/numass-main/src/main/java/inr/numass/NumassPlugin.java @@ -106,7 +106,6 @@ public class NumassPlugin extends BasicPlugin { double to = an.getDouble("to", 19010d); RangedNamedSetSpectrum beta = new BetaSpectrum(context.io().getFile("FS.txt")); ModularSpectrum sp = new ModularSpectrum(beta, A, from, to); - sp.setCaching(false); NBkgSpectrum spectrum = new NBkgSpectrum(sp); return new XYModel(spectrum, getAdapter(an)); @@ -125,7 +124,6 @@ public class NumassPlugin extends BasicPlugin { } NBkgSpectrum spectrum = new NBkgSpectrum(sp); - sp.setCaching(false); return new XYModel(spectrum, getAdapter(an)); }); @@ -212,9 +210,9 @@ public class NumassPlugin extends BasicPlugin { BivariateFunction resolutionTail = ResolutionFunction.getRealTail(); RangedNamedSetSpectrum beta = new BetaSpectrum(context.io().getFile("FS.txt")); ModularSpectrum sp = new ModularSpectrum(beta, new ResolutionFunction(A, resolutionTail), from, to); - if (!an.getBoolean("caching", false)) { - context.getReport().report("Caching turned off"); - sp.setCaching(false); + if (an.getBoolean("caching", false)) { + context.getReport().report("Caching turned on"); + sp.setCaching(true); } //Adding trapping energy dependence //Intercept = 4.95745, B1 = -0.36879, B2 = 0.00827 @@ -244,7 +242,6 @@ public class NumassPlugin extends BasicPlugin { }; RangedNamedSetSpectrum beta = new BetaSpectrum(context.io().getFile("FS.txt")); ModularSpectrum sp = new ModularSpectrum(beta, new ResolutionFunction(A, reolutionTail), from, to); - sp.setCaching(false); NBkgSpectrum spectrum = new NBkgSpectrum(sp); return new XYModel(spectrum, getAdapter(an)); diff --git a/numass-main/src/main/java/inr/numass/actions/PileupSimulationAction.java b/numass-main/src/main/java/inr/numass/actions/PileupSimulationAction.java index 042e0386..c939a53a 100644 --- a/numass-main/src/main/java/inr/numass/actions/PileupSimulationAction.java +++ b/numass-main/src/main/java/inr/numass/actions/PileupSimulationAction.java @@ -15,6 +15,7 @@ import inr.numass.storage.NMPoint; import inr.numass.storage.NumassData; import inr.numass.storage.RawNMPoint; import inr.numass.utils.PileUpSimulator; +import inr.numass.utils.TritiumUtils; import java.time.Instant; import java.util.ArrayList; import java.util.LinkedHashMap; @@ -32,7 +33,7 @@ public class PileupSimulationAction extends OneToOneAction execute(Context context, Reportable log, String name, Laminate inputMeta, NumassData input) { int lowerChannel = inputMeta.getInt("lowerChannel", 1); - int upperChannel = inputMeta.getInt("upperChannel", RawNMPoint.MAX_CHANEL-1); + int upperChannel = inputMeta.getInt("upperChannel", RawNMPoint.MAX_CHANEL - 1); List generated = new ArrayList<>(); List registered = new ArrayList<>(); @@ -44,7 +45,7 @@ public class PileupSimulationAction extends OneToOneAction { double length = point.getLength() * scale; - double cr = point.getCountRate(lowerChannel, upperChannel, 6.4e-6); + double cr = TritiumUtils.countRateWithDeadTime(point, lowerChannel, upperChannel, 6.4e-6); PileUpSimulator simulator = new PileUpSimulator(cr, length) .withGenerator(point, null, lowerChannel, upperChannel) diff --git a/numass-main/src/main/java/inr/numass/actions/PrepareDataAction.java b/numass-main/src/main/java/inr/numass/actions/PrepareDataAction.java index fecd2027..7667c784 100644 --- a/numass-main/src/main/java/inr/numass/actions/PrepareDataAction.java +++ b/numass-main/src/main/java/inr/numass/actions/PrepareDataAction.java @@ -34,6 +34,7 @@ import inr.numass.storage.NMPoint; import inr.numass.storage.NumassData; import inr.numass.storage.RawNMPoint; import inr.numass.utils.ExpressionUtils; +import inr.numass.utils.TritiumUtils; import java.io.OutputStream; import java.time.Instant; import java.util.ArrayList; @@ -103,9 +104,9 @@ public class PrepareDataAction extends OneToOneAction { long wind = point.getCountInWindow(a, b); // count rate after all corrections - double cr = point.getCountRate(a, b, deadTimeFunction.apply(point)); + double cr = TritiumUtils.countRateWithDeadTime(point, a, b, deadTimeFunction.apply(point)); // count rate error after all corrections - double crErr = point.getCountRateErr(a, b, deadTimeFunction.apply(point)); + double crErr = TritiumUtils.countRateWithDeadTimeErr(point, a, b, deadTimeFunction.apply(point)); double correctionFactor = correction(log, point, meta); diff --git a/numass-main/src/main/java/inr/numass/models/ModularSpectrum.java b/numass-main/src/main/java/inr/numass/models/ModularSpectrum.java index 4bf123a0..59417966 100644 --- a/numass-main/src/main/java/inr/numass/models/ModularSpectrum.java +++ b/numass-main/src/main/java/inr/numass/models/ModularSpectrum.java @@ -40,7 +40,7 @@ public class ModularSpectrum extends AbstractParametricFunction { BivariateFunction resolution; RangedNamedSetSpectrum sourceSpectrum; BivariateFunction trappingFunction; - boolean caching = true; + boolean caching = false; double cacheMin; double cacheMax; @@ -67,7 +67,6 @@ public class ModularSpectrum extends AbstractParametricFunction { public ModularSpectrum(RangedNamedSetSpectrum source, BivariateFunction resolution) { this(source, resolution, Double.NaN, Double.NaN); - setCaching(false); } /** diff --git a/numass-main/src/main/java/inr/numass/utils/NMEventGenerator.java b/numass-main/src/main/java/inr/numass/utils/NMEventGenerator.java index 014919bc..85db4225 100644 --- a/numass-main/src/main/java/inr/numass/utils/NMEventGenerator.java +++ b/numass-main/src/main/java/inr/numass/utils/NMEventGenerator.java @@ -122,7 +122,7 @@ public class NMEventGenerator { chanel = 1600; } - return new NMEvent(chanel, prev == null ? 0 : prev.getTime() + nextExpDecay(1 / cr)); + return new NMEvent(chanel, prev == null ? 0 : prev.getTime() + nextExpDecay(1d / cr)); } public double nextExpDecay(double mean) { diff --git a/numass-main/src/main/java/inr/numass/utils/PileUpSimulator.java b/numass-main/src/main/java/inr/numass/utils/PileUpSimulator.java index 7ed78e6f..504eb95a 100644 --- a/numass-main/src/main/java/inr/numass/utils/PileUpSimulator.java +++ b/numass-main/src/main/java/inr/numass/utils/PileUpSimulator.java @@ -103,44 +103,48 @@ public class PileUpSimulator { private boolean random(double prob) { double r = generator.nextUniform(); - return r < prob; + return r <= prob; } public synchronized PileUpSimulator generate() { - NMEvent current = null; + NMEvent last = null;// last event + double lastRegisteredTime = 0; // Time of DAQ closing + //flag that shows that previous event was pileup boolean pileupFlag = false; while (true) { - NMEvent next = generator.nextEvent(current); + NMEvent next = generator.nextEvent(last); if (next.getTime() > pointLength) { break; } - generated.add(next.clone()); - //flag that shows that previous event was pileup + generated.add(next); //not counting double pileups - if (current != null) { - double delay = (next.getTime() - current.getTime()) / us; //time between events in microseconds + if (last != null) { + double delay = (next.getTime() - lastRegisteredTime) / us; //time between events in microseconds if (nextEventRegistered(delay)) { //just register new event - registred.add(next.clone()); + registred.add(next); + lastRegisteredTime = next.getTime(); pileupFlag = false; } else if (pileup(delay) && !pileupFlag) { //pileup event - short newChannel = pileupChannel(delay, current.getChanel(), next.getChanel()); - NMEvent newEvent = new NMEvent(newChannel, current.getTime()); + short newChannel = pileupChannel(delay, last.getChanel(), next.getChanel()); + NMEvent newEvent = new NMEvent(newChannel, last.getTime()); //replace already registered event by event with new channel registred.remove(registred.size() - 1); - registred.add(newEvent.clone()); - pileup.add(newEvent.clone()); + registred.add(newEvent); + pileup.add(newEvent); + //do not change DAQ close time pileupFlag = true; // up the flag to avoid secondary pileup } else { - // second event not registered + // second event not registered, DAQ closed pileupFlag = false; } } else { //register first event - registred.add(next.clone()); + registred.add(next); + lastRegisteredTime = next.getTime(); } - current = next; + last = next; } return this; } diff --git a/numass-main/src/main/java/inr/numass/utils/TritiumUtils.java b/numass-main/src/main/java/inr/numass/utils/TritiumUtils.java index 82c8ce35..a86e3345 100644 --- a/numass-main/src/main/java/inr/numass/utils/TritiumUtils.java +++ b/numass-main/src/main/java/inr/numass/utils/TritiumUtils.java @@ -19,6 +19,7 @@ import hep.dataforge.tables.DataPoint; import hep.dataforge.tables.ListTable; import hep.dataforge.tables.Table; import inr.numass.data.SpectrumDataAdapter; +import inr.numass.storage.NMPoint; import static java.lang.Math.abs; import static java.lang.Math.exp; import static java.lang.Math.sqrt; @@ -30,7 +31,6 @@ import org.apache.commons.math3.analysis.UnivariateFunction; */ public class TritiumUtils { - public static Table correctForDeadTime(ListTable data, double dtime) { return correctForDeadTime(data, adapter(), dtime); } @@ -106,4 +106,21 @@ public class TritiumUtils { double res = Fermi * pe * Etot; return res * 1E-23; } + + public static double countRateWithDeadTime(NMPoint p, int from, int to, double deadTime) { + double wind = p.getCountInWindow(from, to) / p.getLength(); + double res; + if (deadTime > 0) { + double total = p.getEventsCount(); + double time = p.getLength(); + res = wind / (1 - total * deadTime / time); + } else { + res = wind; + } + return res; + } + + public static double countRateWithDeadTimeErr(NMPoint p, int from, int to, double deadTime) { + return Math.sqrt(countRateWithDeadTime(p,from, to, deadTime) / p.getLength()); + } } diff --git a/numass-storage/src/main/java/inr/numass/storage/NMEvent.java b/numass-storage/src/main/java/inr/numass/storage/NMEvent.java index e56c60ce..e7b89851 100644 --- a/numass-storage/src/main/java/inr/numass/storage/NMEvent.java +++ b/numass-storage/src/main/java/inr/numass/storage/NMEvent.java @@ -19,7 +19,7 @@ package inr.numass.storage; * * @author Darksnake */ -public class NMEvent implements Cloneable{ +public class NMEvent{ protected final short chanel; protected final double time; @@ -28,10 +28,10 @@ public class NMEvent implements Cloneable{ this.time = time; } - @Override - public NMEvent clone() { - return new NMEvent(chanel, time); - } +// @Override +// public NMEvent clone() { +// return new NMEvent(chanel, time); +// } /** * @return the chanel diff --git a/numass-storage/src/main/java/inr/numass/storage/NMPoint.java b/numass-storage/src/main/java/inr/numass/storage/NMPoint.java index 448b3e3a..fad39faa 100644 --- a/numass-storage/src/main/java/inr/numass/storage/NMPoint.java +++ b/numass-storage/src/main/java/inr/numass/storage/NMPoint.java @@ -129,25 +129,6 @@ public class NMPoint { return res; } - - //TODO move dead time out of here! - public double getCountRate(int from, int to, double deadTime) { - double wind = getCountInWindow(from, to) / getLength(); - double res; - if (deadTime > 0) { - double total = getEventsCount(); - double time = getLength(); - res = wind / (1 - total * deadTime / time); - } else { - res = wind; - } - return res; - } - - public double getCountRateErr(int from, int to, double deadTime) { - return Math.sqrt(getCountRate(from, to, deadTime) / getLength()); - } - public List getData() { List data = new ArrayList<>(); for (int i = 0; i < RawNMPoint.MAX_CHANEL; i++) { diff --git a/numass-storage/src/main/java/inr/numass/storage/RawNMPoint.java b/numass-storage/src/main/java/inr/numass/storage/RawNMPoint.java index be45ef7f..da652534 100644 --- a/numass-storage/src/main/java/inr/numass/storage/RawNMPoint.java +++ b/numass-storage/src/main/java/inr/numass/storage/RawNMPoint.java @@ -69,7 +69,7 @@ public class RawNMPoint implements Cloneable { public RawNMPoint clone() { ArrayList newevents = new ArrayList<>(); for (NMEvent event : this.getEvents()) { - newevents.add(event.clone()); + newevents.add(event); } return new RawNMPoint(getUset(), getUread(), newevents, getLength()); } diff --git a/numass-viewer/src/main/java/inr/numass/viewer/NumassLoaderViewComponent.java b/numass-viewer/src/main/java/inr/numass/viewer/NumassLoaderViewComponent.java index 1927659a..90f9c8b7 100644 --- a/numass-viewer/src/main/java/inr/numass/viewer/NumassLoaderViewComponent.java +++ b/numass-viewer/src/main/java/inr/numass/viewer/NumassLoaderViewComponent.java @@ -42,6 +42,7 @@ import hep.dataforge.tables.Table; import hep.dataforge.tables.XYAdapter; import inr.numass.storage.NMPoint; import inr.numass.storage.NumassData; +import inr.numass.utils.TritiumUtils; import java.io.File; import java.io.IOException; import java.net.URL; @@ -325,8 +326,8 @@ public class NumassLoaderViewComponent extends AnchorPane implements Initializab private DataPoint getSpectrumPoint(NMPoint point, int lowChannel, int upChannel, double dTime) { double u = point.getUread(); return new MapPoint(new String[]{"x", "y", "yErr"}, u, - point.getCountRate(lowChannel, upChannel, dTime), - point.getCountRateErr(lowChannel, upChannel, dTime)); + TritiumUtils.countRateWithDeadTime(point,lowChannel, upChannel, dTime), + TritiumUtils.countRateWithDeadTimeErr(point,lowChannel, upChannel, dTime)); } /** @@ -402,8 +403,8 @@ public class NumassLoaderViewComponent extends AnchorPane implements Initializab point.getLength(), point.getEventsCount(), point.getCountInWindow(loChannel, upChannel), - point.getCountRate(loChannel, upChannel, dTime), - point.getCountRateErr(loChannel, upChannel, dTime), + TritiumUtils.countRateWithDeadTime(point,loChannel, upChannel, dTime), + TritiumUtils.countRateWithDeadTimeErr(point,loChannel, upChannel, dTime), point.getStartTime() } ));