[no commit message]

This commit is contained in:
Alexander Nozik 2016-07-07 20:40:43 +03:00
parent 1a1c17a1ae
commit 3d6d5f7688
13 changed files with 201 additions and 114 deletions

View File

@ -20,7 +20,7 @@ import hep.dataforge.datafitter.FitManager;
import hep.dataforge.datafitter.ParamSet; import hep.dataforge.datafitter.ParamSet;
import hep.dataforge.datafitter.models.XYModel; import hep.dataforge.datafitter.models.XYModel;
import hep.dataforge.exceptions.NamingException; import hep.dataforge.exceptions.NamingException;
import hep.dataforge.io.PrintNamed; import hep.dataforge.io.FittingIOUtils;
import inr.numass.data.SpectrumDataAdapter; import inr.numass.data.SpectrumDataAdapter;
import inr.numass.models.GunSpectrum; import inr.numass.models.GunSpectrum;
import inr.numass.models.NBkgSpectrum; import inr.numass.models.NBkgSpectrum;
@ -52,7 +52,7 @@ PrintNamed.printSpectrum(new PrintWriter(System.out), spectrum, allPars, 18495,
allPars.setParValue("sigma", 0.6); allPars.setParValue("sigma", 0.6);
PrintNamed.printSpectrum(new PrintWriter(System.out), spectrum, allPars, 18495, 18505, 100); FittingIOUtils.printSpectrum(new PrintWriter(System.out), spectrum, allPars, 18495, 18505, 100);
// //String fileName = "d:\\PlayGround\\merge\\scans.out"; // //String fileName = "d:\\PlayGround\\merge\\scans.out";
//// String configName = "d:\\PlayGround\\SCAN.CFG"; //// String configName = "d:\\PlayGround\\SCAN.CFG";

View File

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

View File

@ -106,7 +106,6 @@ public class NumassPlugin extends BasicPlugin {
double to = an.getDouble("to", 19010d); double to = an.getDouble("to", 19010d);
RangedNamedSetSpectrum beta = new BetaSpectrum(context.io().getFile("FS.txt")); RangedNamedSetSpectrum beta = new BetaSpectrum(context.io().getFile("FS.txt"));
ModularSpectrum sp = new ModularSpectrum(beta, A, from, to); ModularSpectrum sp = new ModularSpectrum(beta, A, from, to);
sp.setCaching(false);
NBkgSpectrum spectrum = new NBkgSpectrum(sp); NBkgSpectrum spectrum = new NBkgSpectrum(sp);
return new XYModel(spectrum, getAdapter(an)); return new XYModel(spectrum, getAdapter(an));
@ -125,7 +124,6 @@ public class NumassPlugin extends BasicPlugin {
} }
NBkgSpectrum spectrum = new NBkgSpectrum(sp); NBkgSpectrum spectrum = new NBkgSpectrum(sp);
sp.setCaching(false);
return new XYModel(spectrum, getAdapter(an)); return new XYModel(spectrum, getAdapter(an));
}); });
@ -212,9 +210,9 @@ public class NumassPlugin extends BasicPlugin {
BivariateFunction resolutionTail = ResolutionFunction.getRealTail(); BivariateFunction resolutionTail = ResolutionFunction.getRealTail();
RangedNamedSetSpectrum beta = new BetaSpectrum(context.io().getFile("FS.txt")); RangedNamedSetSpectrum beta = new BetaSpectrum(context.io().getFile("FS.txt"));
ModularSpectrum sp = new ModularSpectrum(beta, new ResolutionFunction(A, resolutionTail), from, to); ModularSpectrum sp = new ModularSpectrum(beta, new ResolutionFunction(A, resolutionTail), from, to);
if (!an.getBoolean("caching", false)) { if (an.getBoolean("caching", false)) {
context.getReport().report("Caching turned off"); context.getReport().report("Caching turned on");
sp.setCaching(false); sp.setCaching(true);
} }
//Adding trapping energy dependence //Adding trapping energy dependence
//Intercept = 4.95745, B1 = -0.36879, B2 = 0.00827 //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")); RangedNamedSetSpectrum beta = new BetaSpectrum(context.io().getFile("FS.txt"));
ModularSpectrum sp = new ModularSpectrum(beta, new ResolutionFunction(A, reolutionTail), from, to); ModularSpectrum sp = new ModularSpectrum(beta, new ResolutionFunction(A, reolutionTail), from, to);
sp.setCaching(false);
NBkgSpectrum spectrum = new NBkgSpectrum(sp); NBkgSpectrum spectrum = new NBkgSpectrum(sp);
return new XYModel(spectrum, getAdapter(an)); return new XYModel(spectrum, getAdapter(an));

View File

@ -15,6 +15,7 @@ import inr.numass.storage.NMPoint;
import inr.numass.storage.NumassData; import inr.numass.storage.NumassData;
import inr.numass.storage.RawNMPoint; import inr.numass.storage.RawNMPoint;
import inr.numass.utils.PileUpSimulator; import inr.numass.utils.PileUpSimulator;
import inr.numass.utils.TritiumUtils;
import java.time.Instant; import java.time.Instant;
import java.util.ArrayList; import java.util.ArrayList;
import java.util.LinkedHashMap; import java.util.LinkedHashMap;
@ -32,7 +33,7 @@ public class PileupSimulationAction extends OneToOneAction<NumassData, Map<Strin
@Override @Override
protected Map<String, NumassData> execute(Context context, Reportable log, String name, Laminate inputMeta, NumassData input) { protected Map<String, NumassData> execute(Context context, Reportable log, String name, Laminate inputMeta, NumassData input) {
int lowerChannel = inputMeta.getInt("lowerChannel", 1); 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<NMPoint> generated = new ArrayList<>(); List<NMPoint> generated = new ArrayList<>();
List<NMPoint> registered = new ArrayList<>(); List<NMPoint> registered = new ArrayList<>();
@ -44,7 +45,7 @@ public class PileupSimulationAction extends OneToOneAction<NumassData, Map<Strin
input.getNMPoints().forEach(point -> { input.getNMPoints().forEach(point -> {
double length = point.getLength() * scale; 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) PileUpSimulator simulator = new PileUpSimulator(cr, length)
.withGenerator(point, null, lowerChannel, upperChannel) .withGenerator(point, null, lowerChannel, upperChannel)

View File

@ -34,6 +34,7 @@ import inr.numass.storage.NMPoint;
import inr.numass.storage.NumassData; import inr.numass.storage.NumassData;
import inr.numass.storage.RawNMPoint; import inr.numass.storage.RawNMPoint;
import inr.numass.utils.ExpressionUtils; import inr.numass.utils.ExpressionUtils;
import inr.numass.utils.TritiumUtils;
import java.io.OutputStream; import java.io.OutputStream;
import java.time.Instant; import java.time.Instant;
import java.util.ArrayList; import java.util.ArrayList;
@ -103,9 +104,9 @@ public class PrepareDataAction extends OneToOneAction<NumassData, Table> {
long wind = point.getCountInWindow(a, b); long wind = point.getCountInWindow(a, b);
// count rate after all corrections // 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 // 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); double correctionFactor = correction(log, point, meta);

View File

@ -40,7 +40,7 @@ public class ModularSpectrum extends AbstractParametricFunction {
BivariateFunction resolution; BivariateFunction resolution;
RangedNamedSetSpectrum sourceSpectrum; RangedNamedSetSpectrum sourceSpectrum;
BivariateFunction trappingFunction; BivariateFunction trappingFunction;
boolean caching = true; boolean caching = false;
double cacheMin; double cacheMin;
double cacheMax; double cacheMax;
@ -67,7 +67,6 @@ public class ModularSpectrum extends AbstractParametricFunction {
public ModularSpectrum(RangedNamedSetSpectrum source, BivariateFunction resolution) { public ModularSpectrum(RangedNamedSetSpectrum source, BivariateFunction resolution) {
this(source, resolution, Double.NaN, Double.NaN); this(source, resolution, Double.NaN, Double.NaN);
setCaching(false);
} }
/** /**

View File

@ -122,7 +122,7 @@ public class NMEventGenerator {
chanel = 1600; 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) { public double nextExpDecay(double mean) {

View File

@ -103,44 +103,48 @@ public class PileUpSimulator {
private boolean random(double prob) { private boolean random(double prob) {
double r = generator.nextUniform(); double r = generator.nextUniform();
return r < prob; return r <= prob;
} }
public synchronized PileUpSimulator generate() { 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; boolean pileupFlag = false;
while (true) { while (true) {
NMEvent next = generator.nextEvent(current); NMEvent next = generator.nextEvent(last);
if (next.getTime() > pointLength) { if (next.getTime() > pointLength) {
break; break;
} }
generated.add(next.clone()); generated.add(next);
//flag that shows that previous event was pileup
//not counting double pileups //not counting double pileups
if (current != null) { if (last != null) {
double delay = (next.getTime() - current.getTime()) / us; //time between events in microseconds double delay = (next.getTime() - lastRegisteredTime) / us; //time between events in microseconds
if (nextEventRegistered(delay)) { if (nextEventRegistered(delay)) {
//just register new event //just register new event
registred.add(next.clone()); registred.add(next);
lastRegisteredTime = next.getTime();
pileupFlag = false; pileupFlag = false;
} else if (pileup(delay) && !pileupFlag) { } else if (pileup(delay) && !pileupFlag) {
//pileup event //pileup event
short newChannel = pileupChannel(delay, current.getChanel(), next.getChanel()); short newChannel = pileupChannel(delay, last.getChanel(), next.getChanel());
NMEvent newEvent = new NMEvent(newChannel, current.getTime()); NMEvent newEvent = new NMEvent(newChannel, last.getTime());
//replace already registered event by event with new channel //replace already registered event by event with new channel
registred.remove(registred.size() - 1); registred.remove(registred.size() - 1);
registred.add(newEvent.clone()); registred.add(newEvent);
pileup.add(newEvent.clone()); pileup.add(newEvent);
//do not change DAQ close time
pileupFlag = true; // up the flag to avoid secondary pileup pileupFlag = true; // up the flag to avoid secondary pileup
} else { } else {
// second event not registered // second event not registered, DAQ closed
pileupFlag = false; pileupFlag = false;
} }
} else { } else {
//register first event //register first event
registred.add(next.clone()); registred.add(next);
lastRegisteredTime = next.getTime();
} }
current = next; last = next;
} }
return this; return this;
} }

View File

@ -19,6 +19,7 @@ import hep.dataforge.tables.DataPoint;
import hep.dataforge.tables.ListTable; import hep.dataforge.tables.ListTable;
import hep.dataforge.tables.Table; import hep.dataforge.tables.Table;
import inr.numass.data.SpectrumDataAdapter; import inr.numass.data.SpectrumDataAdapter;
import inr.numass.storage.NMPoint;
import static java.lang.Math.abs; import static java.lang.Math.abs;
import static java.lang.Math.exp; import static java.lang.Math.exp;
import static java.lang.Math.sqrt; import static java.lang.Math.sqrt;
@ -30,7 +31,6 @@ import org.apache.commons.math3.analysis.UnivariateFunction;
*/ */
public class TritiumUtils { public class TritiumUtils {
public static Table correctForDeadTime(ListTable data, double dtime) { public static Table correctForDeadTime(ListTable data, double dtime) {
return correctForDeadTime(data, adapter(), dtime); return correctForDeadTime(data, adapter(), dtime);
} }
@ -106,4 +106,21 @@ public class TritiumUtils {
double res = Fermi * pe * Etot; double res = Fermi * pe * Etot;
return res * 1E-23; 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());
}
} }

View File

@ -19,7 +19,7 @@ package inr.numass.storage;
* *
* @author Darksnake * @author Darksnake
*/ */
public class NMEvent implements Cloneable{ public class NMEvent{
protected final short chanel; protected final short chanel;
protected final double time; protected final double time;
@ -28,10 +28,10 @@ public class NMEvent implements Cloneable{
this.time = time; this.time = time;
} }
@Override // @Override
public NMEvent clone() { // public NMEvent clone() {
return new NMEvent(chanel, time); // return new NMEvent(chanel, time);
} // }
/** /**
* @return the chanel * @return the chanel

View File

@ -129,25 +129,6 @@ public class NMPoint {
return res; 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<DataPoint> getData() { public List<DataPoint> getData() {
List<DataPoint> data = new ArrayList<>(); List<DataPoint> data = new ArrayList<>();
for (int i = 0; i < RawNMPoint.MAX_CHANEL; i++) { for (int i = 0; i < RawNMPoint.MAX_CHANEL; i++) {

View File

@ -69,7 +69,7 @@ public class RawNMPoint implements Cloneable {
public RawNMPoint clone() { public RawNMPoint clone() {
ArrayList<NMEvent> newevents = new ArrayList<>(); ArrayList<NMEvent> newevents = new ArrayList<>();
for (NMEvent event : this.getEvents()) { for (NMEvent event : this.getEvents()) {
newevents.add(event.clone()); newevents.add(event);
} }
return new RawNMPoint(getUset(), getUread(), newevents, getLength()); return new RawNMPoint(getUset(), getUread(), newevents, getLength());
} }

View File

@ -42,6 +42,7 @@ import hep.dataforge.tables.Table;
import hep.dataforge.tables.XYAdapter; import hep.dataforge.tables.XYAdapter;
import inr.numass.storage.NMPoint; import inr.numass.storage.NMPoint;
import inr.numass.storage.NumassData; import inr.numass.storage.NumassData;
import inr.numass.utils.TritiumUtils;
import java.io.File; import java.io.File;
import java.io.IOException; import java.io.IOException;
import java.net.URL; 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) { private DataPoint getSpectrumPoint(NMPoint point, int lowChannel, int upChannel, double dTime) {
double u = point.getUread(); double u = point.getUread();
return new MapPoint(new String[]{"x", "y", "yErr"}, u, return new MapPoint(new String[]{"x", "y", "yErr"}, u,
point.getCountRate(lowChannel, upChannel, dTime), TritiumUtils.countRateWithDeadTime(point,lowChannel, upChannel, dTime),
point.getCountRateErr(lowChannel, upChannel, dTime)); TritiumUtils.countRateWithDeadTimeErr(point,lowChannel, upChannel, dTime));
} }
/** /**
@ -402,8 +403,8 @@ public class NumassLoaderViewComponent extends AnchorPane implements Initializab
point.getLength(), point.getLength(),
point.getEventsCount(), point.getEventsCount(),
point.getCountInWindow(loChannel, upChannel), point.getCountInWindow(loChannel, upChannel),
point.getCountRate(loChannel, upChannel, dTime), TritiumUtils.countRateWithDeadTime(point,loChannel, upChannel, dTime),
point.getCountRateErr(loChannel, upChannel, dTime), TritiumUtils.countRateWithDeadTimeErr(point,loChannel, upChannel, dTime),
point.getStartTime() point.getStartTime()
} }
)); ));