Commit 204be0d5 authored by Iain's avatar Iain
Browse files

adding untracked files

parent 4f1563a1
Pipeline #49097 passed with stages
in 4 minutes and 26 seconds
#include "Spline3Refit.h"
#include "../WaveformStatistic.h"
#include "TROOT.h"
#include "TSpline.h"
#include "fstream"
#include "sstream"
void addPoint(TF1 *fit, double *x, double *y, int n, int &i, int len) {
if (y[n] == 1.0) {
fit->FixParameter(i, x[n]);
fit->FixParameter(i + 1, y[n]);
} else {
fit->SetParameter(i, x[n]);
fit->SetParameter(i + 1, y[n]);
double cx = 0.1;
double cy = 0.1;
double dx1, dx2 = 0.0;
if (n == len - 1)
dx1 = (x[n] - x[n - 1]) / 2;
else if (n == 0)
dx2 = (x[n + 1] - x[n]) / 2;
fit->SetParLimits(i, x[n] - dx1 * cx, x[n] + dx2 * cx);
fit->SetParLimits(i + 1, y[n] * (1.0 - cy), y[n] * (1.0 + cy));
}
i += 2;
}
void addPointFixed(TF1 *fit, double *x, double *y, int n, int &i, int len) {
fit->FixParameter(i, x[n]);
fit->FixParameter(i + 1, y[n]);
i += 2;
}
void Spline3Refit::loadModel(std::string path_name) {
std::ifstream myfile;
myfile.open(path_name);
x = std::vector<double>();
y = std::vector<double>();
double xx, yy;
int row = 0;
while (!myfile.eof()) {
std::string str;
std::getline(myfile, str);
std::stringstream ss(str);
ss >> xx;
ss >> yy;
x.push_back(xx);
y.push_back(yy);
row++;
}
num_points = row;
myfile.close();
}
/*
* Sets the pulse strategies fitting parameters
* In this implementation the pulse rise, fall and baseline values are left
* fitting parameters
*/
void Spline3Refit::setFitParemeters(TF1 *fit, Pulse &p, Json::Value cfg) {
double riseTimeSigma = cfg["RiseTimeSigma"].asDouble();
double fallTimeTau = cfg["FallTimeTau"].asInt();
double time2Frac = cfg["Time2Frac"].asDouble();
double fallTime2Tau = cfg["FallTime2Tau"].asInt();
if (!num_points) {
loadModel(cfg["SplineModelPath"].asString());
}
fit->ReleaseParameter(0);
fit->SetParameter(0, p.baseline);
fit->FixParameter(1, 0); // first derivitive start point
fit->SetParameter(2, 0); // first derivitive end point
fit->SetParameter(3, 0);
// fit->SetParameter(4, 0);
fit->FixParameter(4, num_points); // num knots
int i = 6;
for (int n = 0; n < num_points; n++)
if (cfg["FloatingKnots"].asBool())
addPoint(fit, x.data(), y.data(), n, i, num_points);
else
addPointFixed(fit, x.data(), y.data(), n, i, num_points);
}
std::string Spline3Refit::getFitFunctionName(int chID) {
return "FSpline3Refit_" + std::to_string(chID);
};
static double SplineFit(double *x_s, double *pars) {
double x = x_s[0];
double baseline = pars[0];
double b1 = pars[1];
double e1 = pars[2];
double b2 = pars[3];
double e2 = 0.0; // pars[4];// TODO: sort this out
const int num_knots = pars[4];
double knots_x[num_knots];
double knots_y[num_knots];
double x_offset = pars[2 * num_knots + 7];
if (x <= pars[6] + x_offset)
return baseline;
for (int i = 0; i < num_knots; i++) {
knots_x[i] = pars[2 * i + 6];
knots_y[i] = pars[2 * i + 7]; // making height zero
}
double a = knots_x[0];
double b = knots_x[num_knots - 1];
if (x >= b + x_offset)
return baseline;
TSpline5 spline5("TestSpline", knots_x, knots_y, num_knots, "b1e1b2e2", b1,
e1, b2, e2);
return baseline + pars[2 * num_knots + 6] * spline5.Eval(x - x_offset);
}
static double SplineMutli(double *x_s, double *pars) {
double x = x_s[0];
double baseline = pars[0];
double b1 = pars[1];
double e1 = pars[2];
double b2 = pars[3];
double e2 = 0.0; // pars[4];// TODO: sort this out
const int num_knots = pars[4];
double knots_x[num_knots];
double knots_y[num_knots];
for (int i = 0; i < num_knots; i++) {
knots_x[i] = pars[2 * i + 6];
knots_y[i] = pars[2 * i + 7]; // making height zero
}
double a = knots_x[0];
double b = knots_x[num_knots - 1];
// TSpline5
// spline("TestSpline",knots_x,knots_y,num_knots,"b1e1b2e2",b1,e1,b2,e2);
TSpline3 spline("TestSpline", knots_x, knots_y, num_knots, "e1", e1);
double val = baseline;
for (int iPulse = 0; iPulse < pars[5]; iPulse++) {
double x_offset = pars[2 * (num_knots + iPulse) + 7];
if (x <= a + x_offset)
continue;
if (x >= b + x_offset)
continue;
val += pars[2 * (num_knots + iPulse) + 6] * spline.Eval(x - x_offset);
}
return val;
}
TF1 *Spline3Refit::createFitFunction(Json::Value cfg, int chID) {
TString fit_name(getFitFunctionName(chID));
int mNFitPulseMax = cfg["NumPulsesFit"].asInt();
TF1 *fitFunction = (TF1 *)gROOT->FindObject(fit_name);
if (fitFunction && gROOT->FindObjectAny(fit_name)) {
// fitFunction->Delete();
} else {
fitFunction =
new TF1(fit_name, SplineMutli, 0., 0.5e-6, 32 + 2 * mNFitPulseMax);
}
return fitFunction;
}
void Spline3Refit::packPulseFit(Pulse *pulse, TF1 *fitFunction,
int fitPulseIndex, double fitLowLimit,
double fitUppLimit, Json::Value cfg) {
if (pulse) {
double riseTimeSigma = cfg["RiseTimeSigma"].asDouble();
int riseTimeSigmaCoefficient = cfg["RiseTimeSigmaCoefficient"].asInt();
double fallTimeTau = cfg["FallTimeTau"].asDouble();
double fallTime2Tau = cfg["FallTime2Tau"].asDouble();
int polarity = cfg["Polarity"].asInt();
// fit->SetParameter(i, -p.amp);
// fit->SetParameter(i+1, p.time);
fitFunction->ReleaseParameter(6 + 2 * (fitPulseIndex + num_points));
fitFunction->FixParameter(6 + 2 * (fitPulseIndex + num_points),
pulse->amp * polarity);
if (polarity > 0)
fitFunction->SetParLimits(6 + 2 * (fitPulseIndex + num_points), 0, 1e10);
else // polarity cant be zero otherwise it breaks other things
fitFunction->SetParLimits(6 + 2 * (fitPulseIndex + num_points), -1e10, 0);
fitFunction->ReleaseParameter(7 + 2 * (fitPulseIndex + num_points));
fitFunction->FixParameter(
7 + 2 * (fitPulseIndex + num_points),
pulse->time); //- riseTimeSigma * riseTimeSigmaCoefficient);
fitFunction->SetParLimits(7 + 2 * (fitPulseIndex + num_points), fitLowLimit,
fitUppLimit);
} else {
fitFunction->FixParameter(6 + 2 * (fitPulseIndex + num_points), 0);
fitFunction->FixParameter(7 + 2 * (fitPulseIndex + num_points), 0);
}
}
void Spline3Refit::packPulseRefit(Pulse *pulse, TF1 *fitFunction,
int fitPulseIndex, double fitLowLimit,
double fitUppLimit, Json::Value cfg) {
if (pulse) {
fitFunction->ReleaseParameter(6 + 2 * (fitPulseIndex + num_points));
fitFunction->FixParameter(6 + 2 * (fitPulseIndex + num_points),
pulse->fit.amp);
fitFunction->ReleaseParameter(7 + 2 * (fitPulseIndex + num_points));
fitFunction->SetParameter(7 + 2 * (fitPulseIndex + num_points),
pulse->fit.time);
} else {
fitFunction->FixParameter(6 + 2 * (fitPulseIndex + num_points), 0);
fitFunction->FixParameter(7 + 2 * (fitPulseIndex + num_points), 0);
}
}
void Spline3Refit::unpackPulseFit(Pulse *pulse, TF1 *fitFunction,
int fitPulseIndex) {
pulse->fit.amp =
fitFunction->GetParameter(6 + 2 * (fitPulseIndex + num_points));
pulse->fit.time =
fitFunction->GetParameter(7 + 2 * (fitPulseIndex + num_points));
pulse->fit.baseline = fitFunction->GetParameter(0);
for (int i = 0; i < 2 * num_points; i++) {
pulse->fit.knots[i] = fitFunction->GetParameter(i + 6);
}
pulse->fit.chi2 = fitFunction->GetChisquare();
pulse->fit.NDF = fitFunction->GetNDF();
pulse->fit.charge = fitFunction->Integral(
pulse->fit.lowLimit,
pulse->fit.highLimit); // TODO: integral is the same for each
// pulse in the pulse group
}
\ No newline at end of file
#ifndef __SPLINE3_REFIT_H_
#define __SPLINE3_REFIT_H_
#include "AdditiveChi2Refit.h"
/**
* An PulseFittingStrategy1 implementation used to analyze and fit Pulses in a
* Waveform.
*
* Overides the PulseFittingStrategy1's 'setFitParemeters' to allow all the
* parameters to be fit
*/
class Spline3Refit : public IAdditiveChi2Refit {
public:
void setFitParemeters(TF1 *fit, Pulse &p, Json::Value cfg) override;
std::string getFitFunctionName(int chID) override;
void unpackPulseFit(Pulse *pulse, TF1 *fitFunction,
int fitPulseIndex) override;
void packPulseFit(Pulse *pulse, TF1 *fitFunction, int fitPulseIndex,
double fitLowLimit, double fitUppLimit,
Json::Value cfg) override;
void packPulseRefit(Pulse *pulse, TF1 *fitFunction, int fitPulseIndex,
double fitLowLimit, double fitUppLimit,
Json::Value cfg) override;
TF1 *createFitFunction(Json::Value cfg, int chID) override;
private:
void loadModel(std::string path_name);
int num_points;
std::vector<double> x;
std::vector<double> y;
};
#endif // __SPLINE3_REFIT_H_
#include "SplineWriter.h"
#include <string>
/**
* SplineWriter constructor.
* Sets up the TFile to write to and creates the branches for each Channel.
*/
SplineWriter::SplineWriter(Json::Value cfg, DataSource dataSource)
: written(false) {
mProcDir = gDirectory;
num_knots = 8;
std::string inputFilePath = dataSource.filePath;
// TODO: this part isn't OS agnostic
std::string run_name = inputFilePath.substr(
inputFilePath.find_last_of("/") + 1, inputFilePath.find(".mid.gz"));
std::string out_file_name =
inputFilePath.substr(0, inputFilePath.find_last_of("/"));
char tmp_name[80];
sprintf(tmp_name, "%s/%s.Spline.TTree.root", out_file_name.c_str(),
run_name.c_str());
mOutputFile = new TFile(std::string(tmp_name).c_str(), "RECREATE");
mTree = new TTree("KnotTree", "Data Results");
mTree->Branch("x", &x, "x/D");
mTree->Branch("y", &y, "y/D");
mTree->Branch("id", &id, "id/I");
mProcDir->cd();
}
/**
* Writes the TTree to disk.
*/
void SplineWriter::write() {
if (!written) {
written = true;
mOutputFile->cd();
mTree->Write();
mOutputFile->Close();
mProcDir->cd();
}
}
/**
* Updates the TTreeBranch for each Channel at each Event.
*/
bool SplineWriter::push(int indexEvent, Channel &ch, int chID, Waveform &wf,
Json::Value cfg) {
bool containsPulses = true;
for (auto pulse = ch.getPulses()->begin(); pulse != ch.getPulses()->end();
pulse++) {
for (int i = 0; i < num_knots; i++) {
x = pulse->fit.knots[i * 2];
y = pulse->fit.knots[i * 2 + 1];
id = i;
mTree->Fill();
}
}
return containsPulses;
}
/**
* Checks the runtime parameters and returns a boolean indicating whether or not
* they contain every value that the algorithm depends on. It's a bit ugly, but
* pre-checking those parameters lets the program fail fast.
*/
bool SplineWriter::checkParams(Json::Value cfg) { return true; }
/**
* SplineWriter deconstructor.
*/
SplineWriter::~SplineWriter() {
delete mProcDir;
delete mOutputFile;
delete mTree;
}
\ No newline at end of file
#ifndef SPLINE_WRITER_H
#define SPLINE_WRITER_H
#include "../reader/DataFile.h"
#include "PulseWritingStrategy.h"
#include "TFile.h"
#include "TROOT.h"
#include "TTree.h"
#include "TTreeBranch.h"
#include <string>
/**
* IPulseWritingStrategy implementation of a TTree.
* A single tree is created with a branch for each Channel, a branch for eventID
* and triggerTime. Each branch for a Channel is a TTreeBranch.
*
* TTree Structure:
* \startuml
* !include ../docs/ttree.puml
* \enduml
*/
class SplineWriter : public IPulseWritingStrategy {
public:
SplineWriter(Json::Value cfg, DataSource dataSource);
void write();
bool push(int indexEvent, Channel &ch, int chID, Waveform &wf,
Json::Value cfg);
bool checkParams(Json::Value cfg);
inline TTree *getTree() { return mTree; };
~SplineWriter();
private:
TDirectory *mProcDir;
TFile *mOutputFile;
TTree *mTree;
double x;
double y;
int id;
bool written;
int num_knots;
};
#endif // SPLINE_WRITER_H
{
"EventLimit": 100,
"DataSource": {
"FilePath": "./ntp/LOLX/run01099_000.mid",
"Device": "V1740",
"Channels": 30,
"Variables": {}
},
"ChannelConfigs": [
{
"PulseFindingStrategy": "PFS",
"PulseAnalysisStrategy": "AdditiveChi2RefitSpline3",
"BaselineStrategy": "MostProbable",
"EventProcessor": "LOLX",
"WritingStrategy": "SplineWriter",
"Variables": {
"__Comment0__": "Pulse Shape",
"Noise": 1.1,
"RiseTimeSigma": 1.854,
"FallTimeTau": 9.11,
"FallTime2Tau": 0.0,
"Polarity": -1,
"Time2Frac": 0,
"RiseTimeSigmaCoefficient": 5.0,
"FallTimeTauCoefficient": 3,
"FallTime2TauCoefficient": 1,
"SplineModelPath": "./experiments/SplineModels/LOLX",
"FloatingKnots": 0,
"__Comment1__": "Baseline Extraction",
"WaveformAmplitudeBinning": 1,
"MaxAmplitude": 16384,
"MinAmplitude": 1,
"RMSCutNoiseCoefficient": 15,
"__Comment2__": "Pulse Finding",
"NumBin": 0,
"PulseBaselineCalcBeginOffset": 15,
"PulseBaselineCalcEndOffset": 5,
"PulseBaselineDifferenceCutCoefficient": 4,
"PulseChargeCutCoefficient": 6,
"PulseChargeCalcOffset": 2,
"PulseStartCalcOffset": 2,
"PulseStartNoiseThresholdCoefficient": 1,
"PulseEndNoiseThresholdCoefficient": 2,
"PulseOppositePolarityAmpCutCoefficient": 0.5,
"MinPulseWidth": 7,
"__Comment3__": "Pulse Fitting",
"FitOption": "QR0",
"TemplateCheck": 0,
"doFit": 1,
"refit": 0,
"NumPulsesFit": 15,
"MinTimeDifference": 0.5e-9,
"MinChiSquareForRefit": 5.5,
"RiseTimeSigmaUpperLimitCoefficient": 9,
"PulseFitLowerBoundOffsetCoefficient": 7,
"PulseFitLowerBoundOffset": 1
}
}
]
}
-6.36796 0.0261292
-2.20208 0.536544
0.827654 1
12.5679 0.536544
18.2486 0.252753
33.3973 0.0261292
75.0561 -0.014704
165.569 0.00775425
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment