diff --git a/include/FileHandler.h b/include/FileHandler.h new file mode 100644 index 0000000..9e14143 --- /dev/null +++ b/include/FileHandler.h @@ -0,0 +1,95 @@ +#pragma once + +#ifndef file_handler_h +#define file_handler_h + +#include "GaussFit.h" +#include + +using std::string; +using std::to_string; + +class FileHandler { +public: + FileHandler(); + FileHandler(string, int n_ = 6); + ~FileHandler(); + +public: + int n = 6, m = 8, pX; + string file; + double adcValue[6][8][2]; + TH1F Left[6][8]; + TH1F Right[6][8]; + +public: + double getADC(TH1F hist); + void solve(); + void save(); + void save(string); + +private: + void init(); +}; + +FileHandler::FileHandler() {} + +FileHandler::FileHandler(string file_, int n_) { + file = file_; + n = n_; +} + +FileHandler::~FileHandler() {} + +void FileHandler::init() { + string L, R; + for (int i = 0; i < n; i++) + for (int j = 0; j < m; 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(), CHANNEL_NUMBER, 0, CHANNEL_NUMBER); + Right[i][j] = TH1F(R.c_str(), R.c_str(), CHANNEL_NUMBER, 0, CHANNEL_NUMBER); + } +} + +double FileHandler::getADC(TH1F hist) { + int n, cnt = 0; + double *parma = new double[3]; + GaussFit *GF = new GaussFit(); + for (int k = 10; k < CHANNEL_NUMBER; k++) { + n = hist.GetBinContent(k); + if (n == 0) continue; + GF->addData(k, n); + // std::cout << k << " " << n << std::endl; + } + parma = GF->fit(); + + return parma[1]; +} + +void FileHandler::solve() { + init(); + readROOTData(file.c_str(), Left, Right, n); + 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() { + string path = rmString(file, ".root") + ".csv"; + save(path); +} + +void FileHandler::save(string path) { + std::ofstream ofs(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] + << std::endl; +} + +#endif diff --git a/include/GaussFit.h b/include/GaussFit.h index 9727fe4..61075ac 100644 --- a/include/GaussFit.h +++ b/include/GaussFit.h @@ -1,16 +1,18 @@ #pragma once -#ifndef utils_h -#define utils_h +#ifndef gauss_fit_h +#define gauss_fit_h +#include "GaussNewton.h" #include "LevenbergMarquardt.h" +#include "clip.h" +#include "utils.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])); + return p[0] * std::exp(-(x - p[1]) * (x - p[1]) / (2 * p[2] * p[2])); } double* GaussianJacobian(double x, double* p) { @@ -28,7 +30,7 @@ public: public: void addData(double x, double y); - double* fit(); + double* fit(int type_ = 0); public: double* parma = new double[3]; @@ -40,31 +42,38 @@ GaussFit::~GaussFit() {} void GaussFit::addData(double x, double y) { data.push_back(Eigen::Vector2d(x, y)); } -double* GaussFit::fit() { +double* GaussFit::fit(int type_) { int n = data.size(); double x, y, sigma; - Eigen::VectorXd mY(n); - Eigen::MatrixXd mX(n, 3); - Eigen::MatrixXd mB(3, 1); + SigmaClip* SC = new SigmaClip(); + data = SC->clip(data); + parma[0] = dataMax(data); + parma[1] = dataAvg(data); + parma[2] = dataStd(data); - 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; + // 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; + // } + // std::cout << parma[0] << " " << parma[1] << ", " << parma[2] << std::endl; + + if (type_ == 0) { + LevenbergMarquardt LM(3, parma, Gaussian, GaussianJacobian); + LM.data = data; + parma = LM.solve(); + } else if (type_ == 1) { + LevenbergMarquardt LM(3, parma, Gaussian, GaussianJacobian, type_ = type_); + LM.data = data; + parma = LM.solve(); + } else { + GaussNewton GN(3, parma, Gaussian, GaussianJacobian); + GN.data = data; + parma = GN.solve(); } - 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(); - + // std::cout << parma[0] << " " << parma[1] << ", " << parma[2] << std::endl; return parma; } diff --git a/include/GaussNewton.h b/include/GaussNewton.h index 222e7a3..b4e1940 100644 --- a/include/GaussNewton.h +++ b/include/GaussNewton.h @@ -1,7 +1,7 @@ #pragma once -#ifndef GNM_h -#define GNM_h +#ifndef gauss_newton_h +#define gauss_newton_h #include #include diff --git a/include/LevenbergMarquardt.h b/include/LevenbergMarquardt.h index 660a1dc..1550a20 100644 --- a/include/LevenbergMarquardt.h +++ b/include/LevenbergMarquardt.h @@ -1,10 +1,9 @@ #pragma once -#ifndef LMM_h -#define LMM_h +#ifndef levenberg_marquardt_h +#define levenberg_marquardt_h #include -#include #include class LevenbergMarquardt { @@ -22,7 +21,7 @@ public: double* parma; double (*Func)(double, double*); double* (*Gunc)(double, double*); - int L, type_, maxIter = 100; + int L, type_, maxIter = 300; std::vector data; private: @@ -77,7 +76,7 @@ double LevenbergMarquardt::calMse(double* parma_) { 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); + mse += std::pow(y - (*Func)(x, parma_), 2); } return mse; @@ -112,7 +111,7 @@ double* LevenbergMarquardt::solve() { 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)); + mu = mu * std::max(1. / 3, 1 - std::pow(2 * cost - 1, 3)); v = 2; } else { mu *= v; diff --git a/include/clip.h b/include/clip.h index 96bf93c..858a4fb 100644 --- a/include/clip.h +++ b/include/clip.h @@ -3,90 +3,62 @@ #ifndef clip_h #define clip_h -#include +#include "utils.h" -#define INF 1e9 +#include class SigmaClip { private: - int maxiters; - double sigma, minValue, maxValue; - double (*cenF)(double *, int) = nullptr; - double (*stdF)(double *, int) = nullptr; + int maxiters = 5; + double sigma = 3, minValue = INF, maxValue = -INF; + double (*cenF)(std::vector data) = nullptr; + double (*stdF)(std::vector data) = nullptr; public: - SigmaClip(double sigma = 3, int maxiters = 5, double (*cenF)(double *, int) = nullptr, - double (*stdF)(double *, int) = nullptr); + SigmaClip(double sigma = 3, int maxiters = 5, + double (*cenF)(std::vector data) = nullptr, + double (*stdF)(std::vector data) = nullptr); ~SigmaClip(); - double *clip(double *data, int L); + std::vector clip(std::vector data); private: - void computeBound(double *data, int L); + void computeBound(std::vector data); }; -double GetAvg(double *data, int L) { - double sum = 0; - for (int i = 0; i < L; i++) sum += data[i]; - return sum / L; -} - -double GetStd(double *data, int L) { - double sum = 0; - double mean = GetAvg(data, L); - for (int i = 0; i < L; i++) sum += (data[i] - mean) * (data[i] - mean); - return sqrt(sum / L); -} - -SigmaClip::SigmaClip(double sigma, int maxiters, double (*cenF)(double *, int), double (*stdF)(double *, int)) { - this->sigma = sigma; - this->maxiters = maxiters; - this->minValue = INF; - this->maxValue = -INF; - - this->cenF = cenF == nullptr ? GetAvg : cenF; - this->stdF = stdF == nullptr ? GetStd : stdF; +SigmaClip::SigmaClip(double sigma_, int maxiters_, + double (*cenF_)(std::vector data), + double (*stdF_)(std::vector data)) { + sigma = sigma_; + maxiters = maxiters_; + cenF = cenF_ == nullptr ? dataAvg : cenF_; + stdF = stdF_ == nullptr ? dataStd : stdF_; } SigmaClip::~SigmaClip() {} -void SigmaClip::computeBound(double *data, int L) { - double std = (*stdF)(data, L); - double mean = (*cenF)(data, L); +void SigmaClip::computeBound(std::vector data) { + double std = (*stdF)(data); + double mean = (*cenF)(data); minValue = mean - (std * sigma); maxValue = mean + (std * sigma); } -double *SigmaClip::clip(double *data, int L) { - int tmp, cnt; - double *value; +std::vector SigmaClip::clip(std::vector data) { + std::vector::iterator itor; + minValue = INF, maxValue = -INF; - for (int i = 1; i <= maxiters; i++) { - computeBound(data, L); - - cnt = 0; - for (int j = 0; j < L; j++) - if (data[j] >= minValue && data[j] <= maxValue) cnt += 1; - - tmp = L, L = cnt, cnt = 0; - value = new double[L]; - for (int j = 0; j < tmp; j++) - if (data[j] >= minValue && data[j] <= maxValue) { - value[cnt] = data[j]; - cnt += 1; - } - - delete data; - data = new double[L]; - for (int j = 0; j < L; j++) data[j] = value[j]; - delete[] value; + for (int k = 1; k <= maxiters; k++) { + computeBound(data); + for (itor = data.begin(); itor != data.end();) { + if ((*itor)(0) < minValue || (*itor)(0) > maxValue) + data.erase(itor); + else + itor++; + } } - value = new double[L + 1]; - value[0] = L; - for (int j = 0; j < L; j++) value[j + 1] = data[j]; - - return value; + return data; } #endif diff --git a/include/csvReader.h b/include/csvReader.h index 1a5f5cb..9ed74d8 100644 --- a/include/csvReader.h +++ b/include/csvReader.h @@ -12,10 +12,10 @@ using std::string; using std::vector; -class csvReader { +class CsvReader { public: - csvReader(string); - ~csvReader(); + CsvReader(string); + ~CsvReader(); public: string file; @@ -24,16 +24,16 @@ public: public: int rows(); void readData(); - string& operator()(int i, int j) const; + string operator()(int i, int j) const; }; -csvReader::csvReader(string file_) { file = file_; } +CsvReader::CsvReader(string file_) { file = file_; } -csvReader::~csvReader() {} +CsvReader::~CsvReader() {} -int csvReader::rows() { return strArray.size(); } +int CsvReader::rows() { return strArray.size(); } -void csvReader::readData() { +void CsvReader::readData() { std::ifstream inFile(file); string lineStr; while (std::getline(inFile, lineStr)) { @@ -45,7 +45,7 @@ void csvReader::readData() { } } -string& csvReader::operator()(int i, int j) const { +string CsvReader::operator()(int i, int j) const { vector lineArray; int n = strArray.size(), m; @@ -57,7 +57,7 @@ string& csvReader::operator()(int i, int j) const { return lineArray.at(j); } -std::ostream& operator<<(std::ostream& cout, csvReader& cR) { +std::ostream& operator<<(std::ostream& cout, CsvReader& cR) { vector lineArray; int n = cR.strArray.size(), m; diff --git a/include/getADC.h b/include/getADC.h deleted file mode 100644 index 620daee..0000000 --- a/include/getADC.h +++ /dev/null @@ -1,20 +0,0 @@ -#include "GaussFit.h" -#include - -#include - -#define CHANNEL_NUMBER 4096 - -double getADC(TH1F *hist) { - int n, cnt = 0; - double *parma = new double[3]; - GaussFit *GF = new GaussFit(); - for (int k = 0; k < CHANNEL_NUMBER; k++) { - n = hist->GetBinContent(k); - if (n == 0) continue; - GF->addData(k, n); - } - parma = GF->fit(); - - return parma[1]; -} diff --git a/include/utils.h b/include/utils.h new file mode 100644 index 0000000..84c9190 --- /dev/null +++ b/include/utils.h @@ -0,0 +1,96 @@ +#pragma once + +#ifndef utils_h +#define utils_h + +#include +#include +#include + +#include +#include + +#define INF 1e9 +#define CHANNEL_NUMBER 4096 + +using std::string; +using std::to_string; + +double dataMax(std::vector data) { + double m = -INF; + for (int i = 0; i < data.size(); i++) { + Eigen::Vector2d &point = data.at(i); + m = std::max(m, point(1)); + } + return m; +} + +double dataAvg(std::vector data) { + int n = 0; + double m = 0; + for (int i = 0; i < data.size(); i++) { + Eigen::Vector2d &point = data.at(i); + n += point(1); + m += point(0) * point(1); + } + return m / n; +} + +double dataStd(std::vector data) { + int n = 0; + double m = 0; + double mu = dataAvg(data); + for (int i = 0; i < data.size(); i++) { + Eigen::Vector2d &point = data.at(i); + n += point(1); + m += std::pow(point(0) - mu, 2) * point(1); + } + 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) { + // 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[6][16]; + + // read adc data + int na, nc; + string adc; + for (int i = 0; i < n; i++) + for (int j = 0; j < m; j++) { + na = i / 2; + nc = j + 2 * m * (i % 2); + adc = "adc" + to_string(na) + "ch" + to_string(nc); + t->SetBranchAddress(adc.c_str(), &dataArray[i][j]); + adc = "adc" + to_string(na) + "ch" + to_string(nc + m); + t->SetBranchAddress(adc.c_str(), &dataArray[i][j + m]); + } + + for (int i = 0; i < ntot; i++) { + t->GetEntry(i); + for (int j = 0; j < n; j++) + for (int k = 0; k < m; k++) { + Left[j][k].Fill(dataArray[j][k]); + Right[j][k].Fill(dataArray[j][k + m]); + } + } +} + +string rmString(string str, string substr) { + int pos; + int len = substr.length(); + while (true) { + pos = str.find(substr); + if (pos < 0) break; + str.erase(pos, len); + } + return str; +} + +#endif diff --git a/main.cpp b/main.cpp index 9dc9730..b49fcf0 100644 --- a/main.cpp +++ b/main.cpp @@ -1,145 +1,32 @@ -#include "csvReader.h" -// #include "getADC.h" -#include -#include -#include +#include "CsvReader.h" +#include "FileHandler.h" #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() { - // TH1F Left[5][8]; - // TH1F Right[5][8]; + int n; + string run; + FileHandler *FH; + CsvReader cR("config2.csv"); - // 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(), CHANNEL_NUMBER, 0, CHANNEL_NUMBER); - // Right[i][j] = TH1F(R.c_str(), R.c_str(), CHANNEL_NUMBER, 0, CHANNEL_NUMBER); - // } - - // readData("F:/NuclearAstroPhy/Q3D-Calibration/2016Q3D/root/raw/201609Q3D1002.root", Left, - // Right); double x = getADC(&Left[0][0]); cout << x << endl; - - csvReader cR("config.csv"); cR.readData(); + n = cR.rows(); + FH = new FileHandler[n - 1]; + + for (int i = 1; i < n; i++) { + run = cR(i, 0); + FH[i - 1].n = 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"); + } + + // FileHandler FH("2016Q3D/root/raw/201609Q3D1270.root"); + // FH.solve(); + // cout << FH.getADC(FH.Left[0][0]) << endl; return 0; }