add: GaussFit, read_data;change: getADC

.gitignore

@@ -8,5 +8,7 @@ build/
 
 # config
 .vscode
+.vs
+*.json
 
 *.code-workspace

include/GaussFit.h

@@ -0,0 +1,71 @@
+#pragma once
+
+#ifndef utils_h
+#define utils_h
+
+#include "LevenbergMarquardt.h"
+
+#include <Eigen/Dense>
+#include <cmath>
+#include <iostream>
+
+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<Eigen::Vector2d> 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

include/GaussNewton.h

@@ -3,9 +3,8 @@
 #ifndef GNM_h
 #define GNM_h
 
-#include "utils.h"
-
 #include <Eigen/Dense>
+#include <iostream>
 
 class GaussNewton {
 public:
@@ -70,8 +69,7 @@ void GaussNewton::calmH_vG() {
 }
 
 double* GaussNewton::solve() {
-    Eigen::VectorXd vX;
+    Eigen::VectorXd vX(L);
-    vX.resize(L);
 
     for (int k = 0; k < maxIter; k++) {
         calmJ_vF();

include/LevenbergMarquardt.h

@@ -3,13 +3,10 @@
 #ifndef LMM_h
 #define LMM_h
 
-#include "utils.h"
-
 #include <Eigen/Dense>
 #include <cmath>
 #include <iostream>
 
-
 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<Eigen::Vector2d> 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::VectorXd vX(L);
-    Eigen::MatrixXd mT, mD = Eigen::MatrixXd::Zero();
+    Eigen::MatrixXd mT, mD(L, L);
-    vX.resize(L);
+    for (int i = 0; i < L; i++) {
-    mD.resize(L, L);
+        for (int j = 0; j < L; j++) mD(i, j) = 0;
-    for (int i = 0; i < L; i++) mD(i, i) = 1;
+        mD(i, i) = 1;
+    }
 
     for (int k = 0; k < maxIter; k++) {
         calmJ_vF();
@@ -111,8 +108,16 @@ 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 = (calMse(parma) - calMse(parmaNew)) / (vX.transpose() * (mu * vX + vG));
-        cost /= 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;

include/getADC.h

@@ -1,192 +1,18 @@
-#include "LevenbergMarquardt.h"
+#include "GaussFit.h"
-#include "utils.h"
-#include <TCanvas.h>
-#include <TF1.h>
-#include <TFile.h>
 #include <TH1F.h>
-#include <TTree.h>
 
-#include <Eigen/Dense>
 #include <iostream>
 
-using namespace std;
+double getADC(TH1F *hist) {
-using Eigen::MatrixXd;
+    int n, cnt = 0;
-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 *parma = new double[3];
-    double sigma, b0, b1, b2;
+    GaussFit *GF = new GaussFit();
-
+    for (int k = 0; k < 300; k++) {
-    cnt = 0;
+        n = hist->GetBinContent(k);
-    for (int i = 2; i < 300; i++) {
-        n = Left->GetBinContent(i);
         if (n == 0) continue;
-        mY(cnt, 0) = log(n);
+        GF->addData(k, n);
-        mX(cnt, 0) = 1, mX(cnt, 1) = i, mX(cnt, 2) = i * i;
-        cnt += 1;
     }
+    parma = GF->fit();
 
-    mB = (mX.transpose() * mX).inverse() * mX.transpose() * mY;
+    return parma[1];
-    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");
 }

include/utils.h

@@ -1,20 +0,0 @@
-#pragma once
-
-#ifndef utils_h
-#define utils_h
-
-#include <cmath>
-
-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

main.cpp

@@ -1,6 +1,142 @@
 #include "getADC.h"
+#include <TFile.h>
+#include <TH1F.h>
+#include <TTree.h>
+
+#include <iostream>
+
+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;
 }