NPM/numass-framework
4
0
Fork 1
Code Issues Pull Requests Packages Projects Releases Wiki Activity

[no commit message]

Browse Source
This commit is contained in:
darksnake 2016-04-25 15:32:23 +03:00
parent 73908683b8
commit dc03dfebf8
2 changed files with 100 additions and 96 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

173
numass-main/src/main/groovy/inr/numass/scripts/Simulate.groovy
View File

@ -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());

23
numass-main/src/main/java/inr/numass/utils/TritiumUtils.java
View File

@ -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;
}
}
}