Skip to content
GitLab
Commit c418962d authored by Jesper Halkjær Jensen's avatar Jesper Halkjær Jensen
Browse files

Hacky [ci skip]

parent 967247d5
Hide whitespace changes
Inline Side-by-side
include/ausa/calibration/geometry/IsotropicGeometryFitter.h
View file @ c418962d
......@@ -53,7 +53,8 @@ namespace AUSA {
*
* @pre The algorithm assumes that bin (i, j) contain the hits in pixel (i, j). No index shifting is applied!
*/
IsotropicGeometryFitter(TH2 &h, std::shared_ptr<DoubleSidedDetector> d, std::string minimizer = "Minuit", std::string algorithm = "", FitType fittype = FitType::NEYMAN);
IsotropicGeometryFitter(std::vector<TH2*>& h, std::vector<std::shared_ptr<DoubleSidedDetector>>& d,
std::string minimizer = "Minuit", std::string algorithm = "", FitType fittype = FitType::NEYMAN);
/**
* Wrapper object containing the result of the fitting.
......@@ -64,6 +65,8 @@ namespace AUSA {
TMatrixD localCorrelation, localCovariance, globalCovariance;
double chi2;
int status;
TH2* h;
std::shared_ptr<DoubleSidedDetector> d;
};
/**
......@@ -76,25 +79,26 @@ namespace AUSA {
*
* @return The optimal center position, scaling constant and the errors on those.
*/
Result fit();
std::vector<Result> fit();
void setdNLocked(bool lock_z=true);
TH2D makeHistogram(Result* result = nullptr) const;
TH2D makeHistogram(double A , const TVector3& v) const;
TH2D makeHistogram(double A , const TVector3& v, TH2 &h, std::shared_ptr<DoubleSidedDetector> d) const;
TH2D histogramDifference(Result *result = nullptr) const;
TH2D histogramDifference(double A , const TVector3& v) const;
TH2D histogramDifference(double A , const TVector3& v, TH2 &h, std::shared_ptr<DoubleSidedDetector> d) const;
private:
TH2& h;
std::shared_ptr<DoubleSidedDetector> d;
std::unique_ptr<Result> last;
std::vector<std::vector<bool>> exclude;
std::vector<TH2*>& hists;
std::vector<std::shared_ptr<DoubleSidedDetector>>& detectors;
double estimateA();
double estimateA(TH2 &h, std::shared_ptr<DoubleSidedDetector> d);
};
}
......
source/calibration/geometry/IsotropicGeometryFitter.cpp
View file @ c418962d
......@@ -22,8 +22,8 @@ using Chi2Calc = std::function<double(double, double)>;
namespace {
class Chi2Function : public ROOT::Math::IMultiGenFunction {
public:
Chi2Function(TH2 &h, std::shared_ptr<DSD> detector, const std::vector<std::vector<bool>>& exclude, Chi2Calc calc)
: h(h), d(detector), exclude(exclude), calc(calc)
Chi2Function(TH2 &h, std::shared_ptr<DSD> detector, /*const std::vector<std::vector<bool>>& exclude,*/ Chi2Calc calc)
: h(h), d(detector), /*exclude(exclude), */calc(calc)
{
up = detector->getOrientation();
cross = detector->getNormal().Cross(up);
......@@ -33,7 +33,7 @@ namespace {
}
virtual ROOT::Math::IBaseFunctionMultiDim *Clone() const {
return new Chi2Function(h, d, exclude, calc);
return new Chi2Function(h, d, /*exclude, */calc);
}
virtual unsigned int NDim() const {
......@@ -49,7 +49,7 @@ namespace {
double chi2 = 0;
for (UInt_t f = 1; f <= d->frontStripCount(); f++) {
for (UInt_t b = 1; b <= d->backStripCount(); b++) {
if (exclude[f-1][b-1]) continue;
// if (exclude[f-1][b-1]) continue;
auto ref = h.GetBinContent(h.GetBin(f, b));
auto val = A*d->getPixelSolidAngle(f, b);
chi2 += calc(ref, val);
......@@ -62,36 +62,73 @@ namespace {
Chi2Calc calc;
std::shared_ptr<DSD> d;
TVector3 up, cross, norm, center;
const std::vector<std::vector<bool>>& exclude;
// const std::vector<std::vector<bool>>& exclude;
};
}
class MultipleChi2Function : public ROOT::Math::IMultiGenFunction {
public:
MultipleChi2Function(std::vector<Chi2Function> funcs): funcs(funcs)
{
IsotropicGeometryFitter::IsotropicGeometryFitter(TH2 &h, std::shared_ptr<DoubleSidedDetector> d, std::string minimizer, std::string algorithm, FitType fittype)
: DSDGeometryFitter(minimizer, algorithm, std::vector<std::string>{"A", "dU", "dC", "dN"}),
h(h), d(d), last(nullptr)
{
setFitType(fittype);
setCutOff(10);
if (h.GetNbinsX() != d->frontStripCount() || h.GetNbinsY() != d->backStripCount()) {
throw std::invalid_argument("Histogram must have " + std::to_string(d->frontStripCount())
+ " x bins and " + std::to_string(d->backStripCount()) + " y bins.");
}
}
virtual ROOT::Math::IBaseFunctionMultiDim *Clone() const {
return new MultipleChi2Function(funcs);
}
virtual unsigned int NDim() const {
return funcs.size()*3 + 1;
}
virtual double DoEval(const double *par) const {
double parameters[4];
parameters[0] = par[0]; // A
double chi2 = 0;
for(size_t i=0; i<funcs.size(); i++) {
parameters[1] = par[i * 3 + 1];
parameters[2] = par[i * 3 + 2];
parameters[3] = par[i * 3 + 3];
chi2 += funcs[i].DoEval(parameters);
}
return chi2;
}
std::vector<Chi2Function> funcs;
};
}
std::vector<std::string> getNames(int size) {
std::vector<std::string> result;
result.emplace_back("A");
for(int i=0; i<size; i++) {
result.emplace_back("dU" + std::to_string(i));
result.emplace_back("dC" + std::to_string(i));
result.emplace_back("dN" + std::to_string(i));
}
return result;
}
IsotropicGeometryFitter::Result IsotropicGeometryFitter::fit() {
IsotropicGeometryFitter::IsotropicGeometryFitter(std::vector<TH2 *>& h,
std::vector<std::shared_ptr<DoubleSidedDetector>>& d,
std::string minimizer, std::string algorithm, FitType fittype) :
DSDGeometryFitter(minimizer, algorithm, getNames(h.size())), last(nullptr), hists(h), detectors(d)
{
setFitType(fittype);
setCutOff(10);
exclude = generateExclusionMap(h, *d);
}
// Calculate scaling factor for max bin.
auto A = estimateA();
auto original = d->getCenter();
d->setCenter(original);
std::vector<IsotropicGeometryFitter::Result> IsotropicGeometryFitter::fit()
{
std::vector<Chi2Function> funcs;
Chi2Calc calc;
if (fitType == Stat::Chi2Type::LOGLIKELIHOOD)
......@@ -101,27 +138,39 @@ IsotropicGeometryFitter::Result IsotropicGeometryFitter::fit() {
else
calc = Stat::chi2Pearson;
Chi2Function func{h,d, exclude, calc};
// Calculate scaling factor for max bin.
auto A = estimateA(*(hists[0]), detectors[0]);
auto r = doFit(func, {A, 0, 0, 0}, {0.1*A, 0.1, 0.1, 0.1});
for(int i=0; i<hists.size(); i++) {
auto& h = hists[i];
auto& d = detectors[i];
// Covariance matrix without the A estimate
TMatrixD V(3,3);
for(int i=0; i<3; i++) {
// We omit the A estimate
V[i][0] = r.cov[i+1][1];
V[i][1] = r.cov[i+1][2];
V[i][2] = r.cov[i+1][3];
auto original = d->getCenter();
d->setCenter(original);
funcs.emplace_back(Chi2Function{*h, d, /*exclude, */calc});
}
TVector3 newCenter{r.values+1};
auto func = MultipleChi2Function(funcs);
auto r = doFit(func, {A, 0, 0, 0}, {0.1*A, 0.1, 0.1, 0.1});
std::vector<IsotropicGeometryFitter::Result> results;
TMatrixD cov(3,3);
std::tie(newCenter, cov) = Stat::transformCoordinates(func.center, func.up, func.cross, func.norm, newCenter, V);
TMatrixD V(3,3);
last = std::make_unique<Result>(Result{r.values[0], r.err[0], newCenter,
{sqrt(cov[0][0]), sqrt(cov[1][1]), sqrt(cov[2][2])}, r.cor, r.cov, cov, r.chi2, r.status});
for(int i=0; i<hists.size(); i++) {
TVector3 newCenter{r.values+1 + 3 * i};
std::tie(newCenter, cov) = Stat::transformCoordinates(func.funcs[i].center, func.funcs[i].up, func.funcs[i].cross, func.funcs[i].norm, newCenter, V);
results.emplace_back(Result{r.values[0], r.err[0], newCenter,
{0,0,0}, r.cor, r.cov, cov, r.chi2, r.status, hists[i], detectors[i]});
}
return *last;
return results;
}
TH2D IsotropicGeometryFitter::makeHistogram(IsotropicGeometryFitter::Result *result) const {
......@@ -132,10 +181,10 @@ TH2D IsotropicGeometryFitter::makeHistogram(IsotropicGeometryFitter::Result *res
result = last.get();
}
return makeHistogram(result -> A, result->center);
return makeHistogram(result -> A, result->center, *result->h, result->d);
}
TH2D IsotropicGeometryFitter::makeHistogram(double A, const TVector3 &v) const {
TH2D IsotropicGeometryFitter::makeHistogram(double A, const TVector3 &v, TH2 &h, std::shared_ptr<DoubleSidedDetector> d) const {
auto original = d->getCenter();
d->setCenter(v);
......@@ -169,11 +218,11 @@ TH2D IsotropicGeometryFitter::histogramDifference(Result *result) const {
result = last.get();
}
return histogramDifference(result -> A, result -> center);
return histogramDifference(result -> A, result -> center, *result->h, result->d);
}
TH2D IsotropicGeometryFitter::histogramDifference(double A, const TVector3 &v) const {
auto hist = makeHistogram(A, v);
TH2D IsotropicGeometryFitter::histogramDifference(double A, const TVector3 &v, TH2 &h, std::shared_ptr<DoubleSidedDetector> d) const {
auto hist = makeHistogram(A, v, h, d);
hist.Add(&h, -1);
UInt_t nF = d->frontStripCount();
......@@ -181,9 +230,9 @@ TH2D IsotropicGeometryFitter::histogramDifference(double A, const TVector3 &v) c
for (UInt_t f = 1; f <= nF; f++) {
for (UInt_t b = 1; b <= nB; b++) {
if (exclude[f-1][b-1]) {
hist.SetBinContent(f, b, 0);
}
// if (exclude[f-1][b-1]) {
// hist.SetBinContent(f, b, 0);
// }
}
}
auto name = "fit_diff_" + d->getName();
......@@ -192,13 +241,13 @@ TH2D IsotropicGeometryFitter::histogramDifference(double A, const TVector3 &v) c
return hist;
}
double IsotropicGeometryFitter::estimateA() {
double IsotropicGeometryFitter::estimateA(TH2 &h, std::shared_ptr<DoubleSidedDetector> d) {
double A = 0;
int count = 0;
for (UInt_t f = 1; f <= d->frontStripCount(); f++) {
for (UInt_t b = 1; b <= d->backStripCount(); b++) {
if (exclude[f-1][b-1]) continue;
// if (exclude[f-1][b-1]) continue;
auto N = h.GetBinContent(f, b);
auto solidAngle = d->getPixelSolidAngle(f, b);
......@@ -213,6 +262,14 @@ double IsotropicGeometryFitter::estimateA() {
}
void IsotropicGeometryFitter::setdNLocked(bool lock_z) {
if(lock_z) fixParameter("dN", 0);
else freeParameter("dN");
if(lock_z) {
for(int i=0; i<hists.size(); i++) {
fixParameter("dN" + std::to_string(i), 0);
}
}
else {
for(int i=0; i<hists.size(); i++) {
freeParameter("dN" + std::to_string(i));
}
}
}
\ No newline at end of file
Supports Markdown
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