diff --git a/numass-main/src/main/groovy/inr/numass/scripts/Simulate.groovy b/numass-main/src/main/groovy/inr/numass/scripts/Simulate.groovy index 0a0003eb..33ac1e0c 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/Simulate.groovy +++ b/numass-main/src/main/groovy/inr/numass/scripts/Simulate.groovy @@ -13,88 +13,91 @@ * 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.points.ListPointSet; -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.ModularTritiumSpectrum; -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; - -/** - * - * @author Darksnake - */ - -setDefault(Locale.US); -new MINUITPlugin().startGlobal(); -// global.loadModule(new MINUITModule()); - -FitManager fm = new FitManager(); - -ModularTritiumSpectrum beta = new ModularTritiumSpectrum(8.3e-5, 14390d, 19001d, null); -beta.setCaching(false); - -NBkgSpectrum spectrum = new NBkgSpectrum(beta); -XYModel model = new XYModel("tritium", spectrum, new SpectrumDataAdapter()); - -ParamSet allPars = new ParamSet(); - -allPars.setParValue("N", 3e5); -//значение 6е-6 соответствует полной интенстивности 6е7 распадов в секунду -//Проблема была в переполнении счетчика событий в генераторе. Заменил на long. Возможно стоит поставить туда число с плавающей точкой -allPars.setParError("N", 6); -allPars.setParDomain("N", 0d, Double.POSITIVE_INFINITY); -allPars.setParValue("bkg", 2d); -allPars.setParError("bkg", 0.03); -allPars.setParValue("E0", 18575.0); -allPars.setParError("E0", 2); -allPars.setParValue("mnu2", 0d); -allPars.setParError("mnu2", 1d); -allPars.setParValue("msterile2", 1000 * 1000); -allPars.setParValue("U2", 0); -allPars.setParError("U2", 1e-4); -allPars.setParDomain("U2", -1d, 1d); -allPars.setParValue("X", 0); -allPars.setParError("X", 0.01); -allPars.setParDomain("X", 0d, Double.POSITIVE_INFINITY); -allPars.setParValue("trap", 0); -allPars.setParError("trap", 0.01d); -allPars.setParDomain("trap", 0d, Double.POSITIVE_INFINITY); - -// PrintNamed.printSpectrum(GlobalContext.out(), spectrum, allPars, 0.0, 18700.0, 600); -//String fileName = "d:\\PlayGround\\merge\\scans.out"; -// String configName = "d:\\PlayGround\\SCAN.CFG"; -// ListPointSet config = OldDataReader.readConfig(configName); -SpectrumGenerator generator = new SpectrumGenerator(model, allPars, 12316); - -ListPointSet data = generator.generateData(DataModelUtils.getUniformSpectrumConfiguration(13500d, 18200, 1e6, 60)); - -// data = data.filter("X", Value.of(15510.0), Value.of(18610.0)); -// allPars.setParValue("X", 0.4); -FitState state = FitManager.buildState(data, model, allPars); -new PlotFitResultAction(GlobalContext.instance(), null).runOne(state); - - -FitState res = fm.runTask(state, "QOW", FitTask.TASK_RUN, "N", "bkg", "E0", "U2"); - - - -res.print(out()); - +package inr.numass.scripts; + +import hep.dataforge.context.GlobalContext; +import static hep.dataforge.context.GlobalContext.out; +import hep.dataforge.points.ListPointSet; +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.ModularTritiumSpectrum; +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; + +/** + * + * @author Darksnake + */ + +setDefault(Locale.US); +new MINUITPlugin().startGlobal(); +// global.loadModule(new MINUITModule()); + +FitManager fm = new FitManager(); + +ModularTritiumSpectrum beta = new ModularTritiumSpectrum(8.3e-5, 13990d, 18600d, null); +//beta.setCaching(false); + +NBkgSpectrum spectrum = new NBkgSpectrum(beta); +XYModel model = new XYModel("tritium", spectrum, new SpectrumDataAdapter()); + +ParamSet allPars = new ParamSet(); + +allPars.setParValue("N", 9e5); +//значение 6е-6 соответствует полной интенстивности 6е7 распадов в секунду +//Проблема была в переполнении счетчика событий в генераторе. Заменил на long. Возможно стоит поставить туда число с плавающей точкой +allPars.setParError("N", 6); +allPars.setParDomain("N", 0d, Double.POSITIVE_INFINITY); +allPars.setParValue("bkg", 2d); +allPars.setParError("bkg", 0.03); +allPars.setParValue("E0", 18575.0); +allPars.setParError("E0", 2); +allPars.setParValue("mnu2", 0d); +allPars.setParError("mnu2", 1d); +allPars.setParValue("msterile2", 1000 * 1000); +allPars.setParValue("U2", 0); +allPars.setParError("U2", 1e-4); +allPars.setParDomain("U2", -1d, 1d); +allPars.setParValue("X", 0); +allPars.setParError("X", 0.01); +allPars.setParDomain("X", 0d, Double.POSITIVE_INFINITY); +allPars.setParValue("trap", 1d); +allPars.setParError("trap", 0.01d); +allPars.setParDomain("trap", 0d, Double.POSITIVE_INFINITY); + +// PrintNamed.printSpectrum(GlobalContext.out(), spectrum, allPars, 0.0, 18700.0, 600); +//String fileName = "d:\\PlayGround\\merge\\scans.out"; +// String configName = "d:\\PlayGround\\SCAN.CFG"; +// ListPointSet config = OldDataReader.readConfig(configName); +SpectrumGenerator generator = new SpectrumGenerator(model, allPars, 12316); + +ListPointSet 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/utils/TritiumUtils.java b/numass-main/src/main/java/inr/numass/utils/TritiumUtils.java index c424f5da..f0fbac2e 100644 --- a/numass-main/src/main/java/inr/numass/utils/TritiumUtils.java +++ b/numass-main/src/main/java/inr/numass/utils/TritiumUtils.java @@ -34,7 +34,6 @@ import static java.lang.Math.abs; * * @author Darksnake */ - public class TritiumUtils { // /** @@ -59,6 +58,9 @@ public class TritiumUtils { // return res; // // } + public static ListPointSet correctForDeadTime(ListPointSet data, double dtime) { + return correctForDeadTime(data, adapter(), dtime); + } /** * Коррекция на мертвое время в секундах @@ -67,17 +69,16 @@ public class TritiumUtils { * @param dtime * @return */ - public static ListPointSet correctForDeadTime(ListPointSet data, double dtime) { - SpectrumDataAdapter reader = adapter(); + public static ListPointSet correctForDeadTime(ListPointSet data, SpectrumDataAdapter adapter, double dtime) { +// SpectrumDataAdapter adapter = adapter(); ListPointSet res = new ListPointSet(data.getFormat()); for (DataPoint dp : data) { - double corrFactor = 1 / (1 - dtime * reader.getCount(dp) /reader.getTime(dp)); - res.add(reader.buildSpectrumDataPoint(reader.getX(dp).doubleValue(), (long) (reader.getCount(dp)*corrFactor),reader.getTime(dp))); + double corrFactor = 1 / (1 - dtime * adapter.getCount(dp) / adapter.getTime(dp)); + res.add(adapter.buildSpectrumDataPoint(adapter.getX(dp).doubleValue(), (long) (adapter.getCount(dp) * corrFactor), adapter.getTime(dp))); } return res; } - /** * Поправка масштаба высокого. * @@ -90,15 +91,15 @@ public class TritiumUtils { ListPointSet res = new ListPointSet(data.getFormat()); for (DataPoint dp : data) { double corrFactor = 1 + beta; - res.add(reader.buildSpectrumDataPoint(reader.getX(dp).doubleValue()*corrFactor, reader.getCount(dp), reader.getTime(dp))); + res.add(reader.buildSpectrumDataPoint(reader.getX(dp).doubleValue() * corrFactor, reader.getCount(dp), reader.getTime(dp))); } return res; } - - public static SpectrumDataAdapter adapter(){ + + public static SpectrumDataAdapter adapter() { return new SpectrumDataAdapter("Uset", "CR", "CRerr", "Time"); } - + /** * Integral beta spectrum background with given amplitude (total count rate * from) @@ -131,5 +132,5 @@ public class TritiumUtils { double Fermi = Fn * (1.002037 - 0.001427 * ve); double res = Fermi * pe * Etot; return res * 1E-23; - } + } }