From 0a8ef2e2527c5e68c033239d9c7cf6df04fa1edd Mon Sep 17 00:00:00 2001 From: YiHui Liu Date: Wed, 6 Jul 2022 23:33:12 +0800 Subject: [PATCH] add: debug mode, GaussFit Func, statistical function; --- CMakeLists.txt | 4 +- include/FileHandler.h | 36 ++++++++------ include/GaussFit.h | 91 +++++++++++++++++++++++++++++------ include/GaussNewton.h | Bin 1897 -> 3910 bytes include/LevenbergMarquardt.h | Bin 3113 -> 6426 bytes include/clip.h | 26 ++++++++-- include/utils.h | 64 ++++++++++++++++++++++-- main.cpp | 15 +++--- 8 files changed, 189 insertions(+), 47 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index fb16bf4..095bb2a 100644 --- a/CMakeLists.txt +++ b/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) diff --git a/include/FileHandler.h b/include/FileHandler.h index 9e14143..fbb36fd 100644 --- a/include/FileHandler.h +++ b/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); - // std::cout << k << " " << n << std::endl; + GF.addData(hist.GetBinCenter(k), n); } - parma = GF->fit(); + parma = GF.fit(); + if (DEBUG) GF.draw(); 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]); - } + readROOTData(file.c_str(), Left, Right, n, m, thMin, thMax); + if (DEBUG) + std::cout << getADC(Left[3][2]) << std::endl; + else + 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] - << std::endl; + 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 61075ac..b2bc160 100644 --- a/include/GaussFit.h +++ b/include/GaussFit.h @@ -7,6 +7,8 @@ #include "LevenbergMarquardt.h" #include "clip.h" #include "utils.h" +#include +#include #include #include @@ -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, sigma; + double x, y; - SigmaClip* SC = new SigmaClip(); - data = SC->clip(data); - parma[0] = dataMax(data); - parma[1] = dataAvg(data); - parma[2] = dataStd(data); + SigmaClip* SC1 = new SigmaClip(1, 5, dataMax2DInd, dataStd2DSQRT); + SigmaClip* SC2 = new SigmaClip(); + data = SC1->clipN(data); + data = SC2->clip(data); - // 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 (getTotal() <= 400) { + for (int i = 0; i < 3; i++) parma[i] = 0; + return parma; + } + + 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 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 diff --git a/include/GaussNewton.h b/include/GaussNewton.h index b4e1940df7090b1f798c6d3b5bb53e4d2ac78cba..6fccc603187698acf6bb3aa8c96efb68441a7493 100644 GIT binary patch literal 3910 zcmd5<%}!H66h3Pc7rwzHB&HQ9qG5x`3Z){N7>ywzq^Y#vrTjE)Ek%NB6IX6Xd?1&` zjW6Lw-$4DonR_~O=l1r;m1!=QnK^ULcmB`)^LtL((vi9x$T9vJIh7OHlcrpkYxbFw zME0@Qz}`O24?y)+y7<+w+mtiu**i6PlhQlN#GX#3ryRd?hek~{g6m(6hUsUjx5J>R+dk^F;>{(6{@|D7$=g?BV zcR6PNbocdUSQqwft!59>%AY_AW0qCkV3$)?=#vnUqnMX@po*=e`i=M3P26kSuJ8H| zC422@1tL@A%2VG1Q1gc8Py%jesq+zMONE@&KaO8G39*N|dBf*u5uUs?4^tG*nq1ME zTh_LRu$!|SV}N?n!-v+E)Jie!U~e~yFSUT{ezcH7POTU4gfrtb8o2D9Eyp!i!p<5; z!Ld|RWG>D*dDk&xbgj=gtFw0-I|aLb8b&9(lGAI!cH0W#q4>V*gS_wci1V6seLS|b;e`$^mEXOU&OQnfgc5WV|Xj1&pkX_~emE>6=^NkSAjqDWO`ZS76FsXrpyF$q!w4CrI*_pR9^XAQLZ4$H7fI$?x+%Rfl6ncCFr)-)ePRL)T zQRti*HMEJ4ubUAby8hJTu>VAy@^I&nhY3H>%pywCn6u!(FkGJ{33Ms&a~T&djLG!G z7p`ppewxP=VG0AM0bYavUWW;b16H2FtaqlN>sY4pmpH2-BT=;CDoJrU_#}Zg?D6f#DsOQ~GSp2u$xV(r zsE#P&bstLhu>k>_JxzH`iFeiYB=u}N<8B(oN;=zqAx{sz1MrY0bs!cQOSw`S)6S`lWChU)%o2 zbi;+5h~{LPoH-}VJvVc)iit6C7v`4Lf;%Ixh9>8iZvtIXj4kv~pri4z z>Jnwg)f42H!*=;TtcAZM{WtX`71WmI%!3}r!8WHc3zJEd;EM%yG^k`*7romXZrnHu zUkFy2p#0g(!&jP0yE-&UysALwm;~9MBm4Y`!b4S@)%2{ou;P*jl&u*OpI+1V{M59X zs+z2;HC>Fr%!bV{iUZ`j-=??N!XBs7I8+L0;Yx+HAwq@Z=twn>#4&8Okr)0^rD0>J OkZM3Q)_CZNk?|L>b2xVZ diff --git a/include/LevenbergMarquardt.h b/include/LevenbergMarquardt.h index 1550a208a20ad512ce969026dcbd7edc2f8ae813..ead604fc72acbdd1ce581b80265009d76e8ff244 100644 GIT binary patch literal 6426 zcmd^E&2Ah;5Uw*42i{;LOK2VAO?E9f1e24&OK?C)7K&t9*2cRI6YnNwceA#!%nNYh zl2_nKICJ3{khs7b0N>X&<^Gu;&&DT&tns*~tE;Q3zxt|X{LkNOawanw%5$0GcP6jo zr5wvx-j#RUwI;co;MoYzPB1^geTnO(Jj3r9)*ZQ(Gx<(lL&FS`3q1R_Sgq@F$d6ZE zO&hpBEt)_(b!D!LbRn1UY7Bcb`2^Pk*!^*2)ekLW zXY1%u{W&ZV(HwEf5W^9Ec5q+dier_@MJfInEVVr>QfVya&^~c!>Dq;@Ox;AFAqMjJ_h~S7=zmU)I zl&qq0*aAu(BvD`i8yN>O7D2+u?9U?f93{`&FH!}s=o{i)uH zDjmC#j^7V|{)+3b@sxL7ZTOgxY9sgG;B7-?-91ELn$vBoNfocJ15H;EqwnDf`D_o_ zpeJHB`TC|3ew<2Lol`j*rBYUxRLVxF)B;^ej$GRWA9(iBHGZs?UayRx_VtT#^7=h{ zw80>;`gRqKyLuMvy4pEMM2^tO48fMCV4yhf!;B+~c8Y4r`ba*??&L)kucinJv1+yY zGQm&^>4lvY_`?|$^F5efNt9?NiYwQ{OC8E|i*VdD>+O4;hw7^!MxJt!} zY;`cloc2PDB}d!1tWL!isl%?;fT@PD@GPMVt{_q?iVhW z^vY~Hz_Y^yUlswa{fv-WPOmp{h23}?4>a@dlK+-DVV99&Vk(2ePDMStfFb&?xjRAa zH|Pm7Mtgz~t4nnH>gA@J?Yb_aZoSrWuJ>MrS=sKbh|IQ1Ka_eL@2EJsbmu&U*L-ud zTuXVOj3wrxN&kN`)8e$@eD|F;bN(DSZ;N<#5zMcCi32=7EUTD#3%m`@fK4KiCzV#< z@MMYCwskUoOKZle)d+S62h%G%m8GHpy5IeqUlyiD;XZ_=KNd&xkX?s`w|>r^GGeq}$JID6KGcYHk~-~m(ogn8aUJ*bY7gHY z`2B(C*2pz=Cq#*zr=FTx%u}&*Y04jOJHG<)yW?%SbN9GgY?*(l_f2s&bo%du-gjS* zm7_Y~fiCnjj$4lHd-#jd_mWn*`*_NH@Ufs=&#l@D;%NJ)cV8W=@o^1pW8{zaSy$0l rfK-=x#MmrrveVXZ-+L~P`Zbcu3?^UZiRP^H@@uf4|K(st6HWgDBaU+^ literal 3113 zcmd5;KX2<)6yN|KfAa3q3Q3`=4!jFU(>{!r8A zODB$Dmbp32yuiJ-T3;@`xupRoZ>U+#$G_He981N=4?q6?+fUjS^RFVsvZ}KG9@%?9 z&nfe7kZvFo>HG|vJ+oz(7-G1Vl#?=I8;Y28QbKS-36n`ml4~z+4Ko7qe_z|y*QJ(L zaF5!Ik!sm2>ji?rJ!7!FTF2M+r8B?PRqr#5>PW|pVf5i82L)PHQ4})itjsVYVI<-y zHk=CQPusZiK1bvt3)=iF*`= z5v)QY?m!-L65Uf?N%;=tysvjKo=v{$mUJi&+{k7FCdxf8Y2+c^pdN;JgCWlFj{@_j zx@W?TfjeaA;^wui5_w0?d$h5ERz9Hi%x?)B&vu)dA*Y0DhCe}bSnVWof_N3;s4$>` z0_igy<8G_6jdUh2s*!c5svkR50N_-N1^U2w`2Xk~h0`1W59U?)R~{ydG-%KD zj3bpfeTKyz>9c}GISNZeK@{%MBiE|K>L|ymy3wc@=minXK<&zGRF!g5x7>2LMcAub zu@ASqtSrVx0sjKMDt+?iEkIY(5pVV~u{oAAbmdPO(`AqsMoB9VtfmsBY~nXjqWfnu zyAqvMdhCz+G=V04b0CpMg`@PGLYcIe?lRGhw!}d?F4>_}((`r@Mn2N;ctD>kMusg$ zLMiW;V$M)=;k`zmyrS18;sgJSkO*(6pewqvqJgueyhiXD1{`u^yHi5v6frHdJH~53 zLXxULiLmV9 wIHLIIc=j;k3egtC5~+W`qeET>kQQ2WI3-!HkZt+M?i^6PLg9f3TyLTM4QNMrl>h($ diff --git a/include/clip.h b/include/clip.h index 858a4fb..0983791 100644 --- a/include/clip.h +++ b/include/clip.h @@ -20,6 +20,7 @@ public: double (*stdF)(std::vector data) = nullptr); ~SigmaClip(); std::vector clip(std::vector data); + std::vector clipN(std::vector data); private: void computeBound(std::vector data); @@ -30,8 +31,8 @@ SigmaClip::SigmaClip(double sigma_, int maxiters_, double (*stdF_)(std::vector data)) { sigma = sigma_; maxiters = maxiters_; - cenF = cenF_ == nullptr ? dataAvg : cenF_; - stdF = stdF_ == nullptr ? dataStd : stdF_; + cenF = cenF_ == nullptr ? dataAvg2D : cenF_; + stdF = stdF_ == nullptr ? dataStd2D : stdF_; } SigmaClip::~SigmaClip() {} @@ -46,8 +47,8 @@ void SigmaClip::computeBound(std::vector data) { std::vector SigmaClip::clip(std::vector data) { std::vector::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 SigmaClip::clip(std::vector data) return data; } +std::vector SigmaClip::clipN(std::vector data) { + std::vector dataN; + std::vector::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 diff --git a/include/utils.h b/include/utils.h index 84c9190..dd50262 100644 --- a/include/utils.h +++ b/include/utils.h @@ -8,15 +8,40 @@ #include #include +#include #include +#define DEBUG 0 #define INF 1e9 #define CHANNEL_NUMBER 4096 +using namespace std; using std::string; using std::to_string; -double dataMax(std::vector data) { +double dataMax(std::vector 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 data) { + double m = 0; + for (int i = 0; i < data.size(); i++) m += data.at(i); + return m / data.size(); +} + +double dataStd(std::vector 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 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 data) { return m; } -double dataAvg(std::vector data) { +double dataAvg2D(std::vector data) { int n = 0; double m = 0; for (int i = 0; i < data.size(); i++) { @@ -36,10 +61,10 @@ double dataAvg(std::vector data) { return m / n; } -double dataStd(std::vector data) { +double dataStd2D(std::vector 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 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 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 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) { diff --git a/main.cpp b/main.cpp index b49fcf0..02aedfe 100644 --- a/main.cpp +++ b/main.cpp @@ -16,17 +16,16 @@ int main() { FH = new FileHandler[n - 1]; for (int i = 1; i < n; i++) { - run = cR(i, 0); - FH[i - 1].n = 5; + if (DEBUG) + 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; }