add: debug mode, GaussFit Func, statistical function;

CMakeLists.txt

@@ -4,7 +4,7 @@ project(Q3D)
 set(Eigen3_DIR D:/Microsoft/vcpkg/installed/x64-windows/share/eigen3)
 SET(CMAKE_TOOLCHAIN_FILE D:/Microsoft/vcpkg/scripts/buildsystems/vcpkg.cmake)
 
-find_package(ROOT REQUIRED)
+find_package(ROOT REQUIRED Spectrum)
 find_package(Eigen3 CONFIG REQUIRED)
 
 include(${ROOT_USE_FILE})
@@ -14,6 +14,6 @@ file(GLOB sources ${PROJECT_SOURCE_DIR}/src/*.cpp)
 file(GLOB headers ${PROJECT_SOURCE_DIR}/include/*.h)
 
 add_executable(Q3D main.cpp ${sources} ${headers})
-target_link_libraries(Q3D PRIVATE ${ROOT_LIBRARIES} Eigen3::Eigen)
+target_link_libraries(Q3D ${ROOT_LIBRARIES} Eigen3::Eigen)
 
 install(TARGETS Q3D DESTINATION F:/NuclearAstroPhy/Q3D-Calibration)

include/FileHandler.h

@@ -12,11 +12,12 @@ using std::to_string;
 class FileHandler {
 public:
     FileHandler();
-    FileHandler(string, int n_ = 6);
+    FileHandler(string, int n_ = 6, int thMin_ = 800, int thMax_ = 4000);
     ~FileHandler();
 
 public:
-    int n = 6, m = 8, pX;
+    int n = 6, m = 8;
+    int thMin, thMax, pX;
     string file;
     double adcValue[6][8][2];
     TH1F Left[6][8];
@@ -34,9 +35,11 @@ private:
 
 FileHandler::FileHandler() {}
 
-FileHandler::FileHandler(string file_, int n_) {
+FileHandler::FileHandler(string file_, int n_, int thMin_, int thMax_) {
     file = file_;
     n = n_;
+    thMin = thMin_;
+    thMax = thMax_;
 }
 
 FileHandler::~FileHandler() {}
@@ -55,26 +58,29 @@ void FileHandler::init() {
 double FileHandler::getADC(TH1F hist) {
     int n, cnt = 0;
     double *parma = new double[3];
-    GaussFit *GF = new GaussFit();
+    GaussFit GF = GaussFit();
     for (int k = 10; k < CHANNEL_NUMBER; k++) {
         n = hist.GetBinContent(k);
         if (n == 0) continue;
-        GF->addData(k, n);
+        GF.addData(hist.GetBinCenter(k), n);
-        // std::cout << k << " " << n << std::endl;
     }
-    parma = GF->fit();
+    parma = GF.fit();
+    if (DEBUG) GF.draw();
 
     return parma[1];
 }
 
 void FileHandler::solve() {
     init();
-    readROOTData(file.c_str(), Left, Right, n);
+    readROOTData(file.c_str(), Left, Right, n, m, thMin, thMax);
-    for (int i = 0; i < n; i++)
+    if (DEBUG)
-        for (int j = 0; j < m; j++) {
+        std::cout << getADC(Left[3][2]) << std::endl;
-            adcValue[i][j][0] = getADC(Left[i][j]);
+    else
-            adcValue[i][j][1] = getADC(Right[i][j]);
+        for (int i = 0; i < n; i++)
-        }
+            for (int j = 0; j < m; j++) {
+                adcValue[i][j][0] = getADC(Left[i][j]);
+                adcValue[i][j][1] = getADC(Right[i][j]);
+            }
 }
 
 void FileHandler::save() {
@@ -88,8 +94,8 @@ void FileHandler::save(string path) {
     ofs << "pX,X4,Bind,Left,Right" << std::endl;
     for (int i = 0; i < n; i++)
         for (int j = 0; j < m; j++)
-            ofs << pX << "," << i << "," << j << "," << adcValue[i][j][0] << "," << adcValue[i][j][1]
+            ofs << pX << "," << i << "," << j << "," << adcValue[i][j][0] << ","
-                << std::endl;
+                << adcValue[i][j][1] << std::endl;
 }
 
 #endif

include/GaussFit.h

@@ -7,6 +7,8 @@
 #include "LevenbergMarquardt.h"
 #include "clip.h"
 #include "utils.h"
+#include <TCanvas.h>
+#include <TH1F.h>
 
 #include <Eigen/Dense>
 #include <iostream>
@@ -31,6 +33,9 @@ public:
 public:
     void addData(double x, double y);
     double* fit(int type_ = 0);
+    double RSquare();
+    int getTotal();
+    void draw(std::string title = "./Figure.png");
 
 public:
     double* parma = new double[3];
@@ -38,26 +43,34 @@ public:
 };
 
 GaussFit::GaussFit() {}
+
 GaussFit::~GaussFit() {}
 
 void GaussFit::addData(double x, double y) { data.push_back(Eigen::Vector2d(x, y)); }
 
 double* GaussFit::fit(int type_) {
-    int n = data.size();
+    double x, y;
-    double x, y, sigma;
 
-    SigmaClip* SC = new SigmaClip();
+    SigmaClip* SC1 = new SigmaClip(1, 5, dataMax2DInd, dataStd2DSQRT);
-    data = SC->clip(data);
+    SigmaClip* SC2 = new SigmaClip();
-    parma[0] = dataMax(data);
+    data = SC1->clipN(data);
-    parma[1] = dataAvg(data);
+    data = SC2->clip(data);
-    parma[2] = dataStd(data);
 
-    // for (int i = 0; i < data.size(); i++) {
+    if (getTotal() <= 400) {
-    //     Eigen::Vector2d &point = data.at(i);
+        for (int i = 0; i < 3; i++) parma[i] = 0;
-    //     x = point(0), y = point(1);
+        return parma;
-    //     std::cout << x << " " << y << std::endl;
+    }
-    // }
+
-    // std::cout << parma[0] << " " << parma[1] << ", " << parma[2] << std::endl;
+    parma[0] = dataMax2D(data);
+    parma[1] = dataAvg2D(data);
+    parma[2] = dataStd2D(data);
+
+    if (DEBUG > 1)
+        for (int i = 0; i < data.size(); i++) {
+            Eigen::Vector2d& point = data.at(i);
+            x = point(0), y = point(1);
+            std::cout << x << " " << y << std::endl;
+        }
 
     if (type_ == 0) {
         LevenbergMarquardt LM(3, parma, Gaussian, GaussianJacobian);
@@ -73,8 +86,58 @@ double* GaussFit::fit(int type_) {
         parma = GN.solve();
     }
 
-    // std::cout << parma[0] << " " << parma[1] << ", " << parma[2] << std::endl;
+    if (DEBUG) std::cout << parma[0] << " " << parma[1] << ", " << parma[2] << std::endl;
+    if (DEBUG) std::cout << RSquare() << std::endl;
+
+    if (RSquare() < 0.5)
+        for (int i = 0; i < 3; i++) parma[i] = 0;
+
     return parma;
 }
 
+double GaussFit::RSquare() {
+    double x, y, mu;
+    double RSS = 0, TSS = 0;
+
+    std::vector<double> yData;
+    for (int i = 0; i < data.size(); i++) yData.push_back(data.at(i)(1));
+    mu = dataAvg(yData);
+
+    for (int i = 0; i < data.size(); i++) {
+        Eigen::Vector2d& point = data.at(i);
+        x = point(0), y = point(1);
+        RSS += std::pow(y - Gaussian(x, parma), 2);
+        TSS += std::pow(y - mu, 2);
+    }
+
+    return 1 - RSS / TSS;
+}
+
+int GaussFit::getTotal() {
+    int sum = 0;
+    for (int i = 0; i < data.size(); i++) {
+        Eigen::Vector2d& point = data.at(i);
+        sum += point(1);
+    }
+    return sum;
+}
+
+void GaussFit::draw(std::string title) {
+    TCanvas* c1 = new TCanvas("Gauss Fit", "Gauss Fit", 0, 0, 1600, 1200);
+    TH1F* h1 = new TH1F("", "Raw Data", CHANNEL_NUMBER, 0, CHANNEL_NUMBER);
+    TH1F* h2 = new TH1F("", "Fit Data", CHANNEL_NUMBER, 0, CHANNEL_NUMBER);
+
+    for (int i = 0; i < data.size(); i++) {
+        Eigen::Vector2d& point = data.at(i);
+        h1->SetBinContent(point(0), point(1));
+    }
+    for (int i = 0; i < CHANNEL_NUMBER; i++) h2->SetBinContent(i, Gaussian(i, parma));
+
+    c1->cd(1);
+    h1->Draw("L");
+    h2->SetLineColor(kRed);
+    h2->Draw("SAME");
+    c1->SaveAs(title.c_str());
+}
+
 #endif

include/clip.h

@@ -20,6 +20,7 @@ public:
               double (*stdF)(std::vector<Eigen::Vector2d> data) = nullptr);
     ~SigmaClip();
     std::vector<Eigen::Vector2d> clip(std::vector<Eigen::Vector2d> data);
+    std::vector<Eigen::Vector2d> clipN(std::vector<Eigen::Vector2d> data);
 
 private:
     void computeBound(std::vector<Eigen::Vector2d> data);
@@ -30,8 +31,8 @@ SigmaClip::SigmaClip(double sigma_, int maxiters_,
                      double (*stdF_)(std::vector<Eigen::Vector2d> data)) {
     sigma = sigma_;
     maxiters = maxiters_;
-    cenF = cenF_ == nullptr ? dataAvg : cenF_;
+    cenF = cenF_ == nullptr ? dataAvg2D : cenF_;
-    stdF = stdF_ == nullptr ? dataStd : stdF_;
+    stdF = stdF_ == nullptr ? dataStd2D : stdF_;
 }
 
 SigmaClip::~SigmaClip() {}
@@ -46,8 +47,8 @@ void SigmaClip::computeBound(std::vector<Eigen::Vector2d> data) {
 
 std::vector<Eigen::Vector2d> SigmaClip::clip(std::vector<Eigen::Vector2d> data) {
     std::vector<Eigen::Vector2d>::iterator itor;
+    
     minValue = INF, maxValue = -INF;
-
     for (int k = 1; k <= maxiters; k++) {
         computeBound(data);
         for (itor = data.begin(); itor != data.end();) {
@@ -61,4 +62,23 @@ std::vector<Eigen::Vector2d> SigmaClip::clip(std::vector<Eigen::Vector2d> data)
     return data;
 }
 
+std::vector<Eigen::Vector2d> SigmaClip::clipN(std::vector<Eigen::Vector2d> data) {
+    std::vector<Eigen::Vector2d> dataN;
+    std::vector<Eigen::Vector2d>::iterator itor;
+
+    minValue = INF, maxValue = -INF;
+    for (int k = 1; k <= maxiters; k++) {
+        computeBound(data);
+        for (itor = data.begin(); itor != data.end();) {
+            if ((*itor)(0) < minValue || (*itor)(0) > maxValue) {
+                dataN.push_back(Eigen::Vector2d((*itor)(0), (*itor)(1)));
+                data.erase(itor);
+            } else
+                itor++;
+        }
+    }
+
+    return dataN;
+}
+
 #endif

include/utils.h

@@ -8,15 +8,40 @@
 #include <TTree.h>
 
 #include <Eigen/Dense>
+#include <ios>
 #include <iostream>
 
+#define DEBUG 0
 #define INF 1e9
 #define CHANNEL_NUMBER 4096
 
+using namespace std;
 using std::string;
 using std::to_string;
 
-double dataMax(std::vector<Eigen::Vector2d> data) {
+double dataMax(std::vector<double> data) {
+    double x, m = -INF;
+    for (int i = 0; i < data.size(); i++) {
+        x = data.at(i);
+        m = std::max(m, x);
+    }
+    return m;
+}
+
+double dataAvg(std::vector<double> data) {
+    double m = 0;
+    for (int i = 0; i < data.size(); i++) m += data.at(i);
+    return m / data.size();
+}
+
+double dataStd(std::vector<double> data) {
+    double m = 0;
+    double mu = dataAvg(data);
+    for (int i = 0; i < data.size(); i++) m += std::pow(data.at(i) - mu, 2);
+    return m / (data.size() - 1);
+}
+
+double dataMax2D(std::vector<Eigen::Vector2d> data) {
     double m = -INF;
     for (int i = 0; i < data.size(); i++) {
         Eigen::Vector2d &point = data.at(i);
@@ -25,7 +50,7 @@ double dataMax(std::vector<Eigen::Vector2d> data) {
     return m;
 }
 
-double dataAvg(std::vector<Eigen::Vector2d> data) {
+double dataAvg2D(std::vector<Eigen::Vector2d> data) {
     int n = 0;
     double m = 0;
     for (int i = 0; i < data.size(); i++) {
@@ -36,10 +61,10 @@ double dataAvg(std::vector<Eigen::Vector2d> data) {
     return m / n;
 }
 
-double dataStd(std::vector<Eigen::Vector2d> data) {
+double dataStd2D(std::vector<Eigen::Vector2d> data) {
     int n = 0;
     double m = 0;
-    double mu = dataAvg(data);
+    double mu = dataAvg2D(data);
     for (int i = 0; i < data.size(); i++) {
         Eigen::Vector2d &point = data.at(i);
         n += point(1);
@@ -48,7 +73,22 @@ double dataStd(std::vector<Eigen::Vector2d> data) {
     return std::sqrt(m / (n - 1));
 }
 
-void readROOTData(const char *fin, TH1F Left[6][8], TH1F Right[6][8], int n = 6, int m = 8) {
+double dataMax2DInd(std::vector<Eigen::Vector2d> data) {
+    double x = 0, m = -INF;
+    for (int i = 0; i < data.size(); i++) {
+        Eigen::Vector2d &point = data.at(i);
+        if (point(1) > m) {
+            x = point(0);
+            m = point(1);
+        }
+    }
+    return x;
+}
+
+double dataStd2DSQRT(std::vector<Eigen::Vector2d> data) { return std::sqrt(dataMax2D(data)); }
+
+void readROOTData(const char *fin, TH1F Left[6][8], TH1F Right[6][8], int n = 6, int m = 8,
+                  int thMin = 800, int thMax = 4000) {
     // read file
     TFile *fRun = new TFile(fin);
     TTree *t = (TTree *)fRun->Get("Tree1");
@@ -61,6 +101,7 @@ void readROOTData(const char *fin, TH1F Left[6][8], TH1F Right[6][8], int n = 6,
 
     // read adc data
     int na, nc;
+    double x1, x2;
     string adc;
     for (int i = 0; i < n; i++)
         for (int j = 0; j < m; j++) {
@@ -76,10 +117,23 @@ void readROOTData(const char *fin, TH1F Left[6][8], TH1F Right[6][8], int n = 6,
         t->GetEntry(i);
         for (int j = 0; j < n; j++)
             for (int k = 0; k < m; k++) {
+                x1 = dataArray[j][k];
+                x2 = dataArray[j][k + m];
+                if ((x1 + x2) < thMin || (x1 + x2) > thMax) continue;
                 Left[j][k].Fill(dataArray[j][k]);
                 Right[j][k].Fill(dataArray[j][k + m]);
             }
     }
+
+    // Rebin
+    int nMax;
+    for (int i = 0; i < n; i++)
+        for (int j = 0; j < m; j++) {
+            nMax = Left[i][j].GetMaximum();
+            if (nMax < 50) Left[i][j] = *(TH1F *)(Left[i][j].Rebin(8));
+            nMax = Right[i][j].GetMaximum();
+            if (nMax < 50) Right[i][j] = *(TH1F *)(Right[i][j].Rebin(8));
+        }
 }
 
 string rmString(string str, string substr) {

main.cpp

@@ -16,17 +16,16 @@ int main() {
     FH = new FileHandler[n - 1];
 
     for (int i = 1; i < n; i++) {
-        run = cR(i, 0);
+        if (DEBUG)
-        FH[i - 1].n = 5;
+            run = "1250";
+        else
+            run = cR(i, 0);
+        FH[i - 1] = FileHandler("2016Q3D/root/raw/201609Q3D" + run + ".root", 5);
         FH[i - 1].pX = stoi(cR(i, 5));
-        FH[i - 1].file = "2016Q3D/root/raw/201609Q3D" + run + ".root";
         FH[i - 1].solve();
-        FH[i - 1].save("result/201609Q3D" + run + ".csv");
+        if (!DEBUG) FH[i - 1].save("result/adc/201609Q3D" + run + ".csv");
+        if (DEBUG) break;
     }
 
-    // FileHandler FH("2016Q3D/root/raw/201609Q3D1270.root");
-    // FH.solve();
-    // cout << FH.getADC(FH.Left[0][0]) << endl;
-
     return 0;
 }