diff --git a/.gitignore b/.gitignore index 2c45ad5..eaac86d 100644 --- a/.gitignore +++ b/.gitignore @@ -8,5 +8,7 @@ build/ # config .vscode +.vs +*.json *.code-workspace diff --git a/Q3D.exe b/Q3D.exe index 84a4999..29f843b 100644 Binary files a/Q3D.exe and b/Q3D.exe differ diff --git a/include/GaussFit.h b/include/GaussFit.h new file mode 100644 index 0000000..9727fe4 --- /dev/null +++ b/include/GaussFit.h @@ -0,0 +1,71 @@ +#pragma once + +#ifndef utils_h +#define utils_h + +#include "LevenbergMarquardt.h" + +#include +#include +#include + +double Gaussian(double x, double* p) { + return p[0] * exp(-(x - p[1]) * (x - p[1]) / (2 * p[2] * p[2])); +} + +double* GaussianJacobian(double x, double* p) { + double* resJ = new double[3]; + resJ[0] = -Gaussian(x, p) / p[0]; + resJ[1] = -(x - p[1]) * Gaussian(x, p) / (p[2] * p[2]); + resJ[2] = -(x - p[1]) * (x - p[1]) * Gaussian(x, p) / (p[2] * p[2] * p[2]); + return resJ; +} + +class GaussFit { +public: + GaussFit(); + ~GaussFit(); + +public: + void addData(double x, double y); + double* fit(); + +public: + double* parma = new double[3]; + std::vector data; +}; + +GaussFit::GaussFit() {} +GaussFit::~GaussFit() {} + +void GaussFit::addData(double x, double y) { data.push_back(Eigen::Vector2d(x, y)); } + +double* GaussFit::fit() { + int n = data.size(); + double x, y, sigma; + + Eigen::VectorXd mY(n); + Eigen::MatrixXd mX(n, 3); + Eigen::MatrixXd mB(3, 1); + + for (int i = 0; i < n; i++) { + Eigen::Vector2d& point = data.at(i); + x = point(0), y = point(1); + mY(i, 0) = log(y); + mX(i, 0) = 1, mX(i, 1) = x, mX(i, 2) = x * x; + } + + mB = (mX.transpose() * mX).inverse() * mX.transpose() * mY; + sigma = -1 / mB(2, 0); + parma[1] = mB(1, 0) * sigma / 2; + parma[0] = exp(mB(0, 0) + parma[1] * parma[1] / sigma); + parma[2] = sqrt(sigma); + + LevenbergMarquardt LM(3, parma, Gaussian, GaussianJacobian); + LM.data = data; + parma = LM.solve(); + + return parma; +} + +#endif diff --git a/include/GaussNewton.h b/include/GaussNewton.h index ce17b48..222e7a3 100644 --- a/include/GaussNewton.h +++ b/include/GaussNewton.h @@ -3,9 +3,8 @@ #ifndef GNM_h #define GNM_h -#include "utils.h" - #include +#include class GaussNewton { public: @@ -70,8 +69,7 @@ void GaussNewton::calmH_vG() { } double* GaussNewton::solve() { - Eigen::VectorXd vX; - vX.resize(L); + Eigen::VectorXd vX(L); for (int k = 0; k < maxIter; k++) { calmJ_vF(); @@ -86,4 +84,4 @@ double* GaussNewton::solve() { return parma; } -#endif \ No newline at end of file +#endif diff --git a/include/LevenbergMarquardt.h b/include/LevenbergMarquardt.h index 25bbbb8..660a1dc 100644 --- a/include/LevenbergMarquardt.h +++ b/include/LevenbergMarquardt.h @@ -3,13 +3,10 @@ #ifndef LMM_h #define LMM_h -#include "utils.h" - #include #include #include - class LevenbergMarquardt { public: LevenbergMarquardt(int L_, double* parma_, double (*Func_)(double, double*), @@ -21,11 +18,11 @@ public: double* solve(); public: - double mu = 1., eps = 1e-5; + double mu = 1., eps = 1e-6; double* parma; double (*Func)(double, double*); double* (*Gunc)(double, double*); - int L, type_, maxIter = 20; + int L, type_, maxIter = 100; std::vector data; private: @@ -75,27 +72,27 @@ void LevenbergMarquardt::calmH_vG() { } double LevenbergMarquardt::calMse(double* parma_) { - int n = data.size(); double x, y, mse = 0; - for (int i = 0; i < n; i++) { + for (int i = 0; i < data.size(); i++) { Eigen::Vector2d& point = data.at(i); x = point(0), y = point(1); mse += pow(y - (*Func)(x, parma_), 2); } - return mse / n; + return mse; } double* LevenbergMarquardt::solve() { double v = 2, cost; - double* parmaNew; + double* parmaNew = new double[L]; - Eigen::VectorXd vX; - Eigen::MatrixXd mT, mD = Eigen::MatrixXd::Zero(); - vX.resize(L); - mD.resize(L, L); - for (int i = 0; i < L; i++) mD(i, i) = 1; + Eigen::VectorXd vX(L); + Eigen::MatrixXd mT, mD(L, L); + for (int i = 0; i < L; i++) { + for (int j = 0; j < L; j++) mD(i, j) = 0; + mD(i, i) = 1; + } for (int k = 0; k < maxIter; k++) { calmJ_vF(); @@ -111,11 +108,19 @@ double* LevenbergMarquardt::solve() { if (vX.norm() <= eps) return parma; for (int i = 0; i < L; i++) parmaNew[i] = parma[i] + vX[i]; - cost = calMse(parma) - calMse(parmaNew); - cost /= vX.transpose() * (mu * vX + vG); + cost = (calMse(parma) - calMse(parmaNew)) / (vX.transpose() * (mu * vX + vG)); + + if (cost > 0) { + for (int i = 0; i < L; i++) parma[i] = parmaNew[i]; + mu = mu * std::max(1. / 3, 1 - pow(2 * cost - 1, 3)); + v = 2; + } else { + mu *= v; + v *= 2; + } } return parma; } -#endif \ No newline at end of file +#endif diff --git a/include/getADC.h b/include/getADC.h index 80d778a..2633e83 100644 --- a/include/getADC.h +++ b/include/getADC.h @@ -1,192 +1,18 @@ -#include "LevenbergMarquardt.h" -#include "utils.h" -#include -#include -#include +#include "GaussFit.h" #include -#include -#include #include -using namespace std; -using Eigen::MatrixXd; -using Eigen::VectorXd; - -void getADC(const char *fin) { - // read file - TFile *fRun = new TFile(fin); - TTree *t = (TTree *)fRun->Get("Tree1"); - - // data numbers - Long64_t ntot = t->GetEntriesFast(); - - // five X-4 block - UInt_t bArray[16]; - UInt_t cArray[16]; - UInt_t dArray[16]; - UInt_t eArray[16]; - UInt_t fArray[16]; - - // calibration coefficient - Float_t CalC[5][8]; - Float_t Calk1[5][8]; - Float_t Calk2[5][8]; - - // read adc data - t->SetBranchAddress("adc0ch0", &bArray[0]); - t->SetBranchAddress("adc0ch1", &bArray[1]); - t->SetBranchAddress("adc0ch2", &bArray[2]); - t->SetBranchAddress("adc0ch3", &bArray[3]); - t->SetBranchAddress("adc0ch4", &bArray[4]); - t->SetBranchAddress("adc0ch5", &bArray[5]); - t->SetBranchAddress("adc0ch6", &bArray[6]); - t->SetBranchAddress("adc0ch7", &bArray[7]); - t->SetBranchAddress("adc0ch8", &bArray[8]); - t->SetBranchAddress("adc0ch9", &bArray[9]); - t->SetBranchAddress("adc0ch10", &bArray[10]); - t->SetBranchAddress("adc0ch11", &bArray[11]); - t->SetBranchAddress("adc0ch12", &bArray[12]); - t->SetBranchAddress("adc0ch13", &bArray[13]); - t->SetBranchAddress("adc0ch14", &bArray[14]); - t->SetBranchAddress("adc0ch15", &bArray[15]); - - t->SetBranchAddress("adc0ch16", &cArray[0]); - t->SetBranchAddress("adc0ch17", &cArray[1]); - t->SetBranchAddress("adc0ch18", &cArray[2]); - t->SetBranchAddress("adc0ch19", &cArray[3]); - t->SetBranchAddress("adc0ch20", &cArray[4]); - t->SetBranchAddress("adc0ch21", &cArray[5]); - t->SetBranchAddress("adc0ch22", &cArray[6]); - t->SetBranchAddress("adc0ch23", &cArray[7]); - t->SetBranchAddress("adc0ch24", &cArray[8]); - t->SetBranchAddress("adc0ch25", &cArray[9]); - t->SetBranchAddress("adc0ch26", &cArray[10]); - t->SetBranchAddress("adc0ch27", &cArray[11]); - t->SetBranchAddress("adc0ch28", &cArray[12]); - t->SetBranchAddress("adc0ch29", &cArray[13]); - t->SetBranchAddress("adc0ch30", &cArray[14]); - t->SetBranchAddress("adc0ch31", &cArray[15]); - - t->SetBranchAddress("adc1ch0", &dArray[0]); - t->SetBranchAddress("adc1ch1", &dArray[1]); - t->SetBranchAddress("adc1ch2", &dArray[2]); - t->SetBranchAddress("adc1ch3", &dArray[3]); - t->SetBranchAddress("adc1ch4", &dArray[4]); - t->SetBranchAddress("adc1ch5", &dArray[5]); - t->SetBranchAddress("adc1ch6", &dArray[6]); - t->SetBranchAddress("adc1ch7", &dArray[7]); - t->SetBranchAddress("adc1ch8", &dArray[8]); - t->SetBranchAddress("adc1ch9", &dArray[9]); - t->SetBranchAddress("adc1ch10", &dArray[10]); - t->SetBranchAddress("adc1ch11", &dArray[11]); - t->SetBranchAddress("adc1ch12", &dArray[12]); - t->SetBranchAddress("adc1ch13", &dArray[13]); - t->SetBranchAddress("adc1ch14", &dArray[14]); - t->SetBranchAddress("adc1ch15", &dArray[15]); - - t->SetBranchAddress("adc1ch16", &eArray[0]); - t->SetBranchAddress("adc1ch17", &eArray[1]); - t->SetBranchAddress("adc1ch18", &eArray[2]); - t->SetBranchAddress("adc1ch19", &eArray[3]); - t->SetBranchAddress("adc1ch20", &eArray[4]); - t->SetBranchAddress("adc1ch21", &eArray[5]); - t->SetBranchAddress("adc1ch22", &eArray[6]); - t->SetBranchAddress("adc1ch23", &eArray[7]); - t->SetBranchAddress("adc1ch24", &eArray[8]); - t->SetBranchAddress("adc1ch25", &eArray[9]); - t->SetBranchAddress("adc1ch26", &eArray[10]); - t->SetBranchAddress("adc1ch27", &eArray[11]); - t->SetBranchAddress("adc1ch28", &eArray[12]); - t->SetBranchAddress("adc1ch29", &eArray[13]); - t->SetBranchAddress("adc1ch30", &eArray[14]); - t->SetBranchAddress("adc1ch31", &eArray[15]); - - t->SetBranchAddress("adc2ch0", &fArray[0]); - t->SetBranchAddress("adc2ch1", &fArray[1]); - t->SetBranchAddress("adc2ch2", &fArray[2]); - t->SetBranchAddress("adc2ch3", &fArray[3]); - t->SetBranchAddress("adc2ch4", &fArray[4]); - t->SetBranchAddress("adc2ch5", &fArray[5]); - t->SetBranchAddress("adc2ch6", &fArray[6]); - t->SetBranchAddress("adc2ch7", &fArray[7]); - t->SetBranchAddress("adc2ch8", &fArray[8]); - t->SetBranchAddress("adc2ch9", &fArray[9]); - t->SetBranchAddress("adc2ch10", &fArray[10]); - t->SetBranchAddress("adc2ch11", &fArray[11]); - t->SetBranchAddress("adc2ch12", &fArray[12]); - t->SetBranchAddress("adc2ch13", &fArray[13]); - t->SetBranchAddress("adc2ch14", &fArray[14]); - t->SetBranchAddress("adc2ch15", &fArray[15]); - - TH1F *Left = new TH1F("Left", "Left", 300, 0, 300); - TH1F *Right = new TH1F("Right", "Right", 300, 0, 300); - - for (int i = 0; i < ntot; i++) { - t->GetEntry(i); - Left->Fill(bArray[0]); - Right->Fill(bArray[8]); - } - - // use matrix and least square method to get coefficient - // SigmaClip *clip = new SigmaClip(5); - int cnt = 0; - double vX[300]; - for (int i = 2; i < 300; i++) - if (Left->GetBinContent(i) > 0) { - vX[cnt] = i; - cnt += 1; - } - - VectorXd mY(cnt); - MatrixXd mX(cnt, 3); - MatrixXd mB(3, 1); - double n; +double getADC(TH1F *hist) { + int n, cnt = 0; double *parma = new double[3]; - double sigma, b0, b1, b2; - - cnt = 0; - for (int i = 2; i < 300; i++) { - n = Left->GetBinContent(i); + GaussFit *GF = new GaussFit(); + for (int k = 0; k < 300; k++) { + n = hist->GetBinContent(k); if (n == 0) continue; - mY(cnt, 0) = log(n); - mX(cnt, 0) = 1, mX(cnt, 1) = i, mX(cnt, 2) = i * i; - cnt += 1; + GF->addData(k, n); } + parma = GF->fit(); - mB = (mX.transpose() * mX).inverse() * mX.transpose() * mY; - b0 = mB(0, 0), b1 = mB(1, 0), b2 = mB(2, 0); - sigma = -1 / b2; - parma[1] = b1 * sigma / 2; - parma[0] = exp(b0 + parma[1] * parma[1] / sigma); - parma[2] = sqrt(sigma); - - cout << parma[0] << ", " << parma[1] << ", " << parma[2] << endl; - - LevenbergMarquardt LM(3, parma, Gaussian, GaussianJacobian); - for (int i = 2; i < 300; i++) { - n = Left->GetBinContent(i); - if (n == 0) continue; - LM.addData(i, n); - } - parma = LM.solve(); - cout << parma[0] << ", " << parma[1] << ", " << parma[2] << endl; - - // TCanvas *c1 = new TCanvas("c1", "Fitting with Gaussian function"); - // c1->Divide(1, 2); - - // c1->cd(1); - // c1->SetGrid(); - // Left->Draw(); - - // c1->cd(2); - // c1->SetGrid(); - // Right->Draw(); - - // TF1 *GausFunc = new TF1("GausFunc", "[0]*exp(-(x-[1])*(x-[1])/2.0/[2]/[2])", 0, 300); - // GausFunc->SetParLimits(0, 1000, 2400); - // GausFunc->SetParLimits(1, 100, 200); - // GausFunc->SetParLimits(2, 0, 15); - // Left->Fit(GausFunc, "R"); - // Right->Fit(GausFunc, "R"); + return parma[1]; } diff --git a/include/utils.h b/include/utils.h deleted file mode 100644 index 06f3942..0000000 --- a/include/utils.h +++ /dev/null @@ -1,20 +0,0 @@ -#pragma once - -#ifndef utils_h -#define utils_h - -#include - -double Gaussian(double x, double* p) { - return p[0] * exp(-(x - p[1]) * (x - p[1]) / (2 * p[2] * p[2])); -} - -double* GaussianJacobian(double x, double *p) { - double* resJ = new double[3]; - resJ[0] = -Gaussian(x, p) / p[0]; - resJ[1] = -(x - p[1]) * Gaussian(x, p) / (p[2] * p[2]); - resJ[2] = -(x - p[1]) * (x - p[1]) * Gaussian(x, p) / (p[2] * p[2] * p[2]); - return resJ; -} - -#endif diff --git a/main.cpp b/main.cpp index 02eeceb..81210e2 100644 --- a/main.cpp +++ b/main.cpp @@ -1,6 +1,142 @@ #include "getADC.h" +#include +#include +#include + +#include + +using namespace std; +using std::string; +using std::to_string; + +void readData(const char *fin, TH1F Left[5][8], TH1F Right[5][8]) { + // read file + TFile *fRun = new TFile(fin); + TTree *t = (TTree *)fRun->Get("Tree1"); + + // data numbers + Long64_t ntot = t->GetEntriesFast(); + + // five X-4 block + UInt_t dataArray[5][16]; + + // calibration coefficient + Float_t CalC[5][8]; + Float_t Calk1[5][8]; + Float_t Calk2[5][8]; + + // read adc data + t->SetBranchAddress("adc0ch0", &dataArray[0][0]); + t->SetBranchAddress("adc0ch1", &dataArray[0][1]); + t->SetBranchAddress("adc0ch2", &dataArray[0][2]); + t->SetBranchAddress("adc0ch3", &dataArray[0][3]); + t->SetBranchAddress("adc0ch4", &dataArray[0][4]); + t->SetBranchAddress("adc0ch5", &dataArray[0][5]); + t->SetBranchAddress("adc0ch6", &dataArray[0][6]); + t->SetBranchAddress("adc0ch7", &dataArray[0][7]); + t->SetBranchAddress("adc0ch8", &dataArray[0][8]); + t->SetBranchAddress("adc0ch9", &dataArray[0][9]); + t->SetBranchAddress("adc0ch10", &dataArray[0][10]); + t->SetBranchAddress("adc0ch11", &dataArray[0][11]); + t->SetBranchAddress("adc0ch12", &dataArray[0][12]); + t->SetBranchAddress("adc0ch13", &dataArray[0][13]); + t->SetBranchAddress("adc0ch14", &dataArray[0][14]); + t->SetBranchAddress("adc0ch15", &dataArray[0][15]); + + t->SetBranchAddress("adc0ch16", &dataArray[1][0]); + t->SetBranchAddress("adc0ch17", &dataArray[1][1]); + t->SetBranchAddress("adc0ch18", &dataArray[1][2]); + t->SetBranchAddress("adc0ch19", &dataArray[1][3]); + t->SetBranchAddress("adc0ch20", &dataArray[1][4]); + t->SetBranchAddress("adc0ch21", &dataArray[1][5]); + t->SetBranchAddress("adc0ch22", &dataArray[1][6]); + t->SetBranchAddress("adc0ch23", &dataArray[1][7]); + t->SetBranchAddress("adc0ch24", &dataArray[1][8]); + t->SetBranchAddress("adc0ch25", &dataArray[1][9]); + t->SetBranchAddress("adc0ch26", &dataArray[1][10]); + t->SetBranchAddress("adc0ch27", &dataArray[1][11]); + t->SetBranchAddress("adc0ch28", &dataArray[1][12]); + t->SetBranchAddress("adc0ch29", &dataArray[1][13]); + t->SetBranchAddress("adc0ch30", &dataArray[1][14]); + t->SetBranchAddress("adc0ch31", &dataArray[1][15]); + + t->SetBranchAddress("adc1ch0", &dataArray[2][0]); + t->SetBranchAddress("adc1ch1", &dataArray[2][1]); + t->SetBranchAddress("adc1ch2", &dataArray[2][2]); + t->SetBranchAddress("adc1ch3", &dataArray[2][3]); + t->SetBranchAddress("adc1ch4", &dataArray[2][4]); + t->SetBranchAddress("adc1ch5", &dataArray[2][5]); + t->SetBranchAddress("adc1ch6", &dataArray[2][6]); + t->SetBranchAddress("adc1ch7", &dataArray[2][7]); + t->SetBranchAddress("adc1ch8", &dataArray[2][8]); + t->SetBranchAddress("adc1ch9", &dataArray[2][9]); + t->SetBranchAddress("adc1ch10", &dataArray[2][10]); + t->SetBranchAddress("adc1ch11", &dataArray[2][11]); + t->SetBranchAddress("adc1ch12", &dataArray[2][12]); + t->SetBranchAddress("adc1ch13", &dataArray[2][13]); + t->SetBranchAddress("adc1ch14", &dataArray[2][14]); + t->SetBranchAddress("adc1ch15", &dataArray[2][15]); + + t->SetBranchAddress("adc1ch16", &dataArray[3][0]); + t->SetBranchAddress("adc1ch17", &dataArray[3][1]); + t->SetBranchAddress("adc1ch18", &dataArray[3][2]); + t->SetBranchAddress("adc1ch19", &dataArray[3][3]); + t->SetBranchAddress("adc1ch20", &dataArray[3][4]); + t->SetBranchAddress("adc1ch21", &dataArray[3][5]); + t->SetBranchAddress("adc1ch22", &dataArray[3][6]); + t->SetBranchAddress("adc1ch23", &dataArray[3][7]); + t->SetBranchAddress("adc1ch24", &dataArray[3][8]); + t->SetBranchAddress("adc1ch25", &dataArray[3][9]); + t->SetBranchAddress("adc1ch26", &dataArray[3][10]); + t->SetBranchAddress("adc1ch27", &dataArray[3][11]); + t->SetBranchAddress("adc1ch28", &dataArray[3][12]); + t->SetBranchAddress("adc1ch29", &dataArray[3][13]); + t->SetBranchAddress("adc1ch30", &dataArray[3][14]); + t->SetBranchAddress("adc1ch31", &dataArray[3][15]); + + t->SetBranchAddress("adc2ch0", &dataArray[4][0]); + t->SetBranchAddress("adc2ch1", &dataArray[4][1]); + t->SetBranchAddress("adc2ch2", &dataArray[4][2]); + t->SetBranchAddress("adc2ch3", &dataArray[4][3]); + t->SetBranchAddress("adc2ch4", &dataArray[4][4]); + t->SetBranchAddress("adc2ch5", &dataArray[4][5]); + t->SetBranchAddress("adc2ch6", &dataArray[4][6]); + t->SetBranchAddress("adc2ch7", &dataArray[4][7]); + t->SetBranchAddress("adc2ch8", &dataArray[4][8]); + t->SetBranchAddress("adc2ch9", &dataArray[4][9]); + t->SetBranchAddress("adc2ch10", &dataArray[4][10]); + t->SetBranchAddress("adc2ch11", &dataArray[4][11]); + t->SetBranchAddress("adc2ch12", &dataArray[4][12]); + t->SetBranchAddress("adc2ch13", &dataArray[4][13]); + t->SetBranchAddress("adc2ch14", &dataArray[4][14]); + t->SetBranchAddress("adc2ch15", &dataArray[4][15]); + + for (int i = 0; i < ntot; i++) { + t->GetEntry(i); + for (int j = 0; j < 5; j++) + for (int k = 0; k < 8; k++) { + Left[j][k].Fill(dataArray[j][k]); + Right[j][k].Fill(dataArray[j][k + 8]); + } + } +} int main() { - getADC("2016Q3D/root/raw/201609Q3D1002.root"); + TH1F Left[5][8]; + TH1F Right[5][8]; + + string L, R; + for (int i = 0; i < 5; i++) + for (int j = 0; j < 8; j++) { + L = "Left" + to_string(i) + to_string(j); + R = "Right" + to_string(i) + to_string(j); + Left[i][j] = TH1F(L.c_str(), L.c_str(), 300, 0, 300); + Right[i][j] = TH1F(R.c_str(), R.c_str(), 300, 0, 300); + } + + readData("F:/NuclearAstroPhy/Q3D-Calibration/2016Q3D/root/raw/201609Q3D1002.root", Left, Right); + double x = getADC(&Left[0][0]); + cout << x << endl; + return 0; }