From d9cb151ce8022f4c412cf0d79cf637ae3f30e4c9 Mon Sep 17 00:00:00 2001 From: Maps Date: Thu, 12 Dec 2024 11:36:56 +0100 Subject: [PATCH] pulsing stuff --- scripts/pulse/FitSCurves.cpp | 101 ++++++++++++---- scripts/pulse/FitSCurvesVariableStep.cpp | 141 +++++++++++++++++++++++ scripts/pulse/FitSCurves_4Jun.cpp | 84 ++++++++++++++ scripts/pulse/demo.root | Bin 0 -> 9919 bytes scripts/pulse/fit-raw.py | 121 +++++++++---------- scripts/pulse/plot-raw.py | 6 +- 6 files changed, 366 insertions(+), 87 deletions(-) create mode 100644 scripts/pulse/FitSCurvesVariableStep.cpp create mode 100644 scripts/pulse/FitSCurves_4Jun.cpp create mode 100644 scripts/pulse/demo.root diff --git a/scripts/pulse/FitSCurves.cpp b/scripts/pulse/FitSCurves.cpp index 1d3ce17..8c3a379 100644 --- a/scripts/pulse/FitSCurves.cpp +++ b/scripts/pulse/FitSCurves.cpp @@ -4,9 +4,14 @@ #include "TMath.h" #include "TString.h" #include "TGraph.h" +#include "TH1I.h" #include "TH1D.h" #include #include +#include "TF1.h" +#include "TFile.h" +#include "TText.h" +#include "TCanvas.h" using namespace std; @@ -14,12 +19,10 @@ TH1I* split(string text, char delim) { string line; stringstream ss(text); - TH1I* h = new TH1I("h", "h",255,0,254); + TH1I* h = new TH1I("h", "h",255,0,255); Int_t counter = 0; while(getline(ss, line, delim)) { TString lline (line); -// cout << counter << endl; -// cout << lline.Atoi() << endl; h->SetBinContent(counter, lline.Atoi()); ++counter; } @@ -31,20 +34,41 @@ TH1I* split(string text, char delim) { int FitSCurves(TString mat = "A", TString f = "VCASNA-116.csv", TString dir = "VCASNA-100", Float_t offset = 0., Float_t slopex = 1.0) { vector v_x(255); // Make x axis for hists - iota(begin(v_x), end(v_x), 0); - - TH1I* Fpn = new TH1I("FPNHist", "FPNHist", 255, 0, 254); - TH1I* ThN = new TH1I("ThermNoiseHist", "ThermNoiseHist", 255, 0, 254); - - + + TH1D* Fpn = new TH1D("FPNHist", "FPNHist", 255, 0, 255); + TH1D* ThN = new TH1D("ThermNoiseHist", "ThermNoiseHist", 40, 0, 20); + cout << "Reading data from CSV..." << endl; ifstream ifs(f); string line; Int_t linecounter = 0; - Int_t failedFit = 0; + Int_t failedFit = 0; + Int_t noisy = 0; + Int_t dead = 0; + + Int_t noisyCondition = (255 - 50) * 500; // Pixel noisy if it responds 100 % of time above 50 electrons + Int_t deadCondition = (255 - 200) * 500; // Pixel dead if it responds less than 100 % of times after 205 electrons + while (getline(ifs, line)) { + + TH1I* SCurve = split(line, '\t'); + + if (SCurve->Integral() > noisyCondition) { // new + // new + ++noisy; // new + break; // new + // new + } // new + // new + if (SCurve->Integral() < deadCondition) { // new + // new + ++dead; // new + break; // new + // new + } // new + TF1 *ferrf = new TF1("ferrf", "0.5*[0]*(1+TMath::Erf((x-[1])/(TMath::Sqrt2()*[2])))",0,250); ferrf->SetParameter(0,500); ferrf->SetParameter(1, 150); @@ -55,21 +79,13 @@ int FitSCurves(TString mat = "A", TString f = "VCASNA-116.csv", TString dir = "V ferrf->SetParLimits(1,80,220); ferrf->SetParLimits(2,0.1,25); - - - TH1I* SCurve = split(line, '\t'); - - try{SCurve->Fit("ferrf");} + try{ + SCurve->Fit("ferrf"); + Fpn->Fill(ferrf->GetParameter(1)); + ThN->Fill(ferrf->GetParameter(2)); + } catch (...) {++failedFit;} - //catch(++failedFit;) - - //Float_t height = ferrf->GetParameter(0); - //Float_t mean = ferrf->GetParameter(1); - //Float_t tn = ferrf->GetParameter(2); - //cout << height <SetParameter(1, 150); + gausFpn->SetParameter(2, 10); + gausThn->SetParameter(1, 5); + gausThn->SetParameter(2, 4); + ThN->Draw(); + ThN->Fit("gausThn"); + TString thnText = "Thermal noise: " + to_string(gausThn->GetParameter(1) * slopex); + TText* mylatexThn = new TText(0,0,thnText); + mylatexThn->Draw(); + t1->SaveAs("test.png"); + + TCanvas *t2 = new TCanvas(); + Fpn->Draw(); + Fpn->Fit("gausFpn"); + TString fpnText = "Mean: " + to_string(gausFpn->GetParameter(1)*slopex + offset) + ", FPN: " + to_string(gausFpn->GetParameter(2)*slopex) + ", Dead: " + to_string(dead) + ", Noisy: " + to_string(noisy); + //TString thnText = "Thermal noise: " + to_string(gausThn->GetParameter(1) * slopex); + TText* mylatexFpn = new TText(2,6,fpnText); + mylatexFpn->SetTextSize(0.02); + mylatexFpn->Draw(); + t2->SaveAs("FPN.png"); + + cout << gausThn->GetParameter(0) << endl; + cout << gausThn->GetParameter(1) << endl; + + cout << "Thermal Noise: " << gausThn->GetParameter(1) * slopex << endl; + cout << "Mean Threshold: " << gausFpn->GetParameter(1) * slopex + offset << endl; + cout << "Fixed Pattern Noise: " << gausFpn->GetParameter(2) * slopex << endl; + cout << "Dead: " << dead << endl; + cout << "Noisy: " << noisy << endl; cout << "Data read. Fitting now..." << endl; + + TFile* t = new TFile("demo.root", "recreate"); + t->WriteObject(Fpn, "FPN"); + t->WriteObject(ThN, "THN"); + t->Close(); return 0; diff --git a/scripts/pulse/FitSCurvesVariableStep.cpp b/scripts/pulse/FitSCurvesVariableStep.cpp new file mode 100644 index 0000000..7326a3b --- /dev/null +++ b/scripts/pulse/FitSCurvesVariableStep.cpp @@ -0,0 +1,141 @@ +#include +#include +#include +#include "TMath.h" +#include "TString.h" +#include "TGraph.h" +#include "TH1I.h" +#include "TH1D.h" +#include +#include +#include "TF1.h" +#include "TFile.h" +#include "TText.h" +#include "TCanvas.h" + +using namespace std; + +TH1I* split(string text, char delim, Int_t stepsize = 5) { + string line; + + stringstream ss(text); + TH1I* h = new TH1I("h", "h",255/stepsize,0,255); + Int_t counter = 0; + while(getline(ss, line, delim)) { + TString lline (line); + if ((counter+3) % stepsize == 0){ + cout << lline << endl; + h->SetBinContent(counter/stepsize, lline.Atoi());} + counter++; + } + return h; + delete h; +} + + +int FitSCurves(TString mat = "A", TString f = "VCASNA-116.csv", TString dir = "VCASNA-100", Float_t offset = 0., Float_t slopex = 1.0, Int_t stepsize=5) { + + vector v_x(255/stepsize); // Make x axis for hists + + TH1D* Fpn = new TH1D("FPNHist", "FPNHist", 255, 0, 255); + TH1D* ThN = new TH1D("ThermNoiseHist", "ThermNoiseHist", 40, 0, 20); + + cout << "Reading data from CSV..." << endl; + + ifstream ifs(f); + + string line; + Int_t linecounter = 0; + Int_t failedFit = 0; + Int_t noisy = 0; + Int_t dead = 0; + Float_t max = 4000.; + + Int_t noisyCondition = (255 - 50.) /stepsize * max; // Pixel noisy if it responds 100 % of time above 50 electrons + Int_t deadCondition = (255 - 200.) /stepsize * max; // Pixel dead if it responds less than 100 % of times after 205 electrons + + while (getline(ifs, line)) { + + TH1I* SCurve = split(line, '\t', stepsize); + + if (SCurve->Integral() > noisyCondition) { // new + // new + ++noisy; // new + break; // new + // new + } // new + // new + if (SCurve->Integral() < deadCondition) { // new + // new + ++dead; // new + break; // new + // new + } // new + + TF1 *ferrf = new TF1("ferrf", "0.5*[0]*(1+TMath::Erf((x-[1])/(TMath::Sqrt2()*[2])))",0,250); + ferrf->SetParameter(0,max); + ferrf->SetParameter(1, 150); + ferrf->SetParameter(2, 5); + + ferrf->SetParameters(max,150.,5.); + ferrf->SetParLimits(0,0,max*1.1); + ferrf->SetParLimits(1,80,220); + ferrf->SetParLimits(2,0.1,25); + + try{ + SCurve->Fit("ferrf"); + Fpn->Fill(ferrf->GetParameter(1)); + ThN->Fill(ferrf->GetParameter(2)); + } + catch (...) {++failedFit;} + + + // Fill thermal noise hist + delete SCurve; + delete ferrf; + + ++linecounter; + } + TCanvas *t1 = new TCanvas(); + + TF1 *gausThn = new TF1("gausThn", "gaus", 1, 20); + TF1 *gausFpn = new TF1("gausFpn", "gaus", 30, 255); + gausFpn->SetParameter(1, 150); + gausFpn->SetParameter(2, 10); + gausThn->SetParameter(1, 5); + gausThn->SetParameter(2, 4); + ThN->Draw(); + ThN->Fit("gausThn"); + TString thnText = "Thermal noise: " + to_string(gausThn->GetParameter(1) * slopex); + TText* mylatexThn = new TText(0,0,thnText); + mylatexThn->Draw(); + t1->SaveAs("test.png"); + + TCanvas *t2 = new TCanvas(); + Fpn->Draw(); + Fpn->Fit("gausFpn"); + TString fpnText = "Mean: " + to_string(gausFpn->GetParameter(1)*slopex + offset) + ", FPN: " + to_string(gausFpn->GetParameter(2)*slopex) + ", Dead: " + to_string(dead) + ", Noisy: " + to_string(noisy); + //TString thnText = "Thermal noise: " + to_string(gausThn->GetParameter(1) * slopex); + TText* mylatexFpn = new TText(2,6,fpnText); + mylatexFpn->SetTextSize(0.02); + mylatexFpn->Draw(); + t2->SaveAs("FPN.png"); + + cout << gausThn->GetParameter(0) << endl; + cout << gausThn->GetParameter(1) << endl; + + cout << "Thermal Noise: " << gausThn->GetParameter(1) * slopex << endl; + cout << "Mean Threshold: " << gausFpn->GetParameter(1) * slopex + offset << endl; + cout << "Fixed Pattern Noise: " << gausFpn->GetParameter(2) * slopex << endl; + cout << "Dead: " << dead << endl; + cout << "Noisy: " << noisy << endl; + cout << "Data read. Fitting now..." << endl; + + TFile* t = new TFile("demo.root", "recreate"); + t->WriteObject(Fpn, "FPN"); + t->WriteObject(ThN, "THN"); + t->Close(); + + return 0; + +} diff --git a/scripts/pulse/FitSCurves_4Jun.cpp b/scripts/pulse/FitSCurves_4Jun.cpp new file mode 100644 index 0000000..9d8c815 --- /dev/null +++ b/scripts/pulse/FitSCurves_4Jun.cpp @@ -0,0 +1,84 @@ +#include +#include +#include +#include "TMath.h" +#include "TString.h" +#include "TGraph.h" +#include "TH1D.h" +#include +#include +#include "TF1.h" + +using namespace std; + +TH1I* split(string text, char delim) { + string line; + + stringstream ss(text); + TH1I* h = new TH1I("h", "h",255,0,254); + Int_t counter = 0; + while(getline(ss, line, delim)) { + TString lline (line); +// cout << counter << endl; +// cout << lline.Atoi() << endl; + h->SetBinContent(counter, lline.Atoi()); + ++counter; + } + return h; + delete h; +} + + +int FitSCurves(TString mat = "A", TString f = "VCASNA-116.csv", TString dir = "VCASNA-100", Float_t offset = 0., Float_t slopex = 1.0) { + + vector v_x(255); // Make x axis for hists + + TH1I* Fpn = new TH1I("FPNHist", "FPNHist", 255, 0, 254); + TH1I* ThN = new TH1I("ThermNoiseHist", "ThermNoiseHist", 255, 0, 254); + + cout << "Reading data from CSV..." << endl; + + ifstream ifs(f); + + string line; + Int_t linecounter = 0; + Int_t failedFit = 0; + while (getline(ifs, line)) { + TF1 *ferrf = new TF1("ferrf", "0.5*[0]*(1+TMath::Erf((x-[1])/(TMath::Sqrt2()*[2])))",0,250); + ferrf->SetParameter(0,500); + ferrf->SetParameter(1, 150); + ferrf->SetParameter(2, 5); + + ferrf->SetParameters(500.,150.,5.); + ferrf->SetParLimits(0,0,550); + ferrf->SetParLimits(1,80,220); + ferrf->SetParLimits(2,0.1,25); + + + + TH1I* SCurve = split(line, '\t'); + + try{SCurve->Fit("ferrf");} + catch (...) {++failedFit;} + //catch(++failedFit;) + + //Float_t height = ferrf->GetParameter(0); + //Float_t mean = ferrf->GetParameter(1); + //Float_t tn = ferrf->GetParameter(2); + + //cout << height <rs)Fu2Rsn~d8+eXEyYT0{{TP1OOoK{cBNw_iTR+%V{huO27Be}TkZ3!jPrH2)(1%?1E~tB6_Jm_nPF+S=1Q{a4BVXA1!T z@AN+p03hz~VfMfF&kFz`d;33?0+IJ$MihpB{*3=JuKef!O+|u91WZgp4$j`l%GB6} z(8bcl##BoN$_xU)fCBIUxm5R5ND@c$_aQ}6I7DHLAj(D&Ak791tir`g5hle4i@*dy zP8md!tCpalYu0-JttiT$$#+sYQ^5s_jA&H{h}TndO=uJses3mZNhx33Vt8Rxs=AEv zZ(`vId)n~5Huc}}?%dh<7MOr)K|J>x6cR)u0OE775V~(ZRYYt`YFvde;sZCo3=<`I zR7+s54WDOWn8XIv9zcsdlm^qGfm2xxaiODWCq-`1NLS~uhIHw2&+Q_;s_$Pfz&h9U zgH~W4^;(1L*I$4gqkb#=@%sGFJ&dO!@=?yJz@*D^4SU^zp zO#T}A8FYE%3D0Q?Yd-Tira8#zAi9_qefk~ZG_pgeZw~|q11q}Z44LwqP2;M@k#vYh z+}eQw59n^sx&XQb#bg6B8I~n&&|2gyX!rGV(Y+=kL_@oM4tTz(?b}U#p(knzKI5oG zSl1b3dm=mJQ6QEGfM>2BIdi|@vp)TQQq&B^@Gc^0@E_l7qzA3Bqb@rsCEpT2O?hg28K(6UFo!;8w z5BGX`&qpX#8*iI zF;y4xx@s*4OllUPdD0*F%`p@cpk@vRv870rXLd&n!&RjJt5;BXie775@|&&3RM7ao zrXzPOr=CCNjkMTg_BygPSVO}~Q>JF}^V72F$iuAuY44%x@6Wf(I{E+ZYoDPfYG^xozdcv}94 ztF!aF8Y_n(zsPy4vm0dBlGm`Pb8vfK-*fZ~InvD69mcoK+D}v7f8@RCw&tX7767{& zz6uP#B2`Tfcrci*Lw|W!vl0c+OyDQd%TH_3k$mj3o)~_#{t$UcBr&H_tsCvZCv$R| z@b0ima8sS{HjnEV^EIly0L#w5Z;}b6l?Bnz*pqBvjls6Dlya(oiz^GzASaw*qT(1d z%&d$qEeaOCYSKW#p-9J<4j+LQV}ltpieP0~a9R2WvMhZJq?(jVq&9(Z`l09^@WYi? z0m^{J!iX1EKWm}}CO`PNSUB-brJXVq2dkZqsZcdZz@DWpKKf_pHCBvuRG>3Rk+x_; z^Nc~-^R0k0)vou6qx2+OQ18BU!fs=!s1LJ(o8o}zA7n7^3?Aa|5kDr22)EV)Zput} zVu>-RF<-))oA$#?PHYBlwqjA|J(zs3%Jwy(1Zvvv2oxW(~hNy_n*s2QD7>HTf_tgeEg~$UjPio3MXi=NVU8Ybbqz+{E>#>?4Y>Zui zizvFop>#}_|7{K^hagF~Wy!e)zYTW>qbbEiF?jwF!J5hWoO$UtZru6^o^)(0QYce3!&cx2}XV3+QHb zwbovGNMR(bzfK5pvZpq=S={~HuiSex;sP!a3{dNk!QueV!NAY>H2}|Bv5ajEy`OHX z9!}Oz4w>{q4n&QrlzSLlZ7^7M*ZnaF@vJy4 zk+$2Q(STJ?G%6~wsM!c>ShsHmD@l-Jfyc~Cp0CAIxCnIf5$=JWS8d7dus-<&SI?aQ z)*qfe9^;rg@WucK0ZVz}t`-Eg9vqYPI#CpWdwLMQ&3;E$=ZQ6}=t$nW}YugiE}BzvF#ugn%2l$PFGD`vzrafy<}j zQFfN|($P88DFqxg1iJNR<#(>(>4);)=q4JovljiFK7}|HI(412o?(D%zP*fB=U+;d zob};eunXlmt-ti&)w^+A&(hz5uNoefe(bX?Z0-|Q>}c(|F_|-Dcy-fQtlvowRl10I z{GsMR>Ec!#b@vk~b0$C(Cukst{rWSW<|{`+N`CJ7vxYrG#2z}3@2#2C=gHBtNg5&n z*OEBVl5OacPBluu5UgGz4M%weFq(H(zp`|D&83vA z*`K@lu(r`wx#ly&iLPtxpfu0n>OTC)%B?Re5=d= zk>>t|E9d+1nC)AWVyMm+%_T^kIglvrO|J&2qwR29xwf|nRs^fm3;pRQK2Y<*QK5_H zU*?idz;R=xU322=Ip{K9F!iYPqSNVHJ@5cs@bd|l+p3|+6zu6t9D?9sJhfVmy4Hd* zEE{eokI{OKJ;eBq4`)HW}< z#ckeepEB#mSD6bHKVOI3I&vFn`+$*MO{|@gth2Y=d{STi!}zHUV=(V<8QzXn-{D~n z6^>ebha}aB-WYVZsTO|kb033C__RlsEPj|8P407(VXW(^>C6lq4DJ+b(OsEt#=7R? zySvk(_5IwaVj8(p#8QvhHK(-R&tn?u#H&+EmFaT>|7_0cWh?cEt-P?3qrAAMJ00)U zZL+hvWw3)f-OK*63jY$i2HCE7C8hIt3(usHu2phij}2W3b_2b7TC1AuYB0+U&z(LU zX=)nJQS&o5g%LOOSMDK0A^&CZ)!nfb`;hHV?;mmH2seeRzfa)LaKQ*2_nXdelhrnV zYGQG*#cE8bLJ}=>PLPxGS!7#uMAy?YvEk<`gT7}5!pqx_|B+N=D?JJKkG4hwEA_EaC#Xt zqeqSNEB~+#wZ!9j&BFt6ffbYL4vWg^HY#iTOvsgcqtks|0m1;)a%ooen#lGjrBq5BiFYK}>pl`~5Jp&($AB^Yg;%1h))BB>vP zKA&3TuX^TpmJOr$M>78cqg?1eU=;X|xBW*w0Vd!oGM3IRaLO)DriQkrPLg(J_9!Cu zu0}SVgf^CT)}|(Z7iX8hJdP9t5GwZ{kJGWX!CiIj;XX{<7fz;>ZlGJg)rh5wD#S0= zPqFzW9c4r-P7HJ=Q%)tw+{z&u&CJo9o{CDQTen_~+#;7vS_01;5%^OoFt2aD85 z^*Io$vw~qm@PzI~OqF4euKNO0FaU+a#+_H|`Oo3t)oFR?fk$0AVaPJlKGzHbvtGjO zmHhq#dPZnQwa)^tbgjojqH&98hs)tdxPoJss03#n>+>5t(agbnCvEQ@=C>+jWuXM< z3Wy!CdHvCs2epFB!U|ZSu2@i!N0=mC(}uhVYmtoHecWGF%vw6^l3Z^2+PSqj9n$iy zwf7C02oJ~|IE%vEE>pf9r3KpBQqyKTc&o}Wf(Y75Srx4 z&>T>y!m2vbGZ_+}-I-@tu;>h;ipguJ_QRw}( z3Tx*cszF-Oo}=ODz2~+`pUKbLV66=x2A%66vhE@Sbk(@!E?e;}gJ6Ydsfzyw;I_mT zOO)VMV=7ZuZ^7uEbGQ%nO5vJ}vSgqOC!Isc2L*RcD3D1#5tb~|gJyG-hoNax3Fj4} z*(R;3Hh>S~ovoufC^b-bFE${HE{e{%J)x|F3GN;mm2=EU>` z5wl>jL^$UN`oi5n-#OZ+BnT-Dr_di~re@AT$81+Tv5yy>>-WS^?pn8y`_Gw$mK~Rn z>_pIP4JF`QJgXEUDb;E7ku~X`y%PP?Td;jkm^TO>X(+<3G@;akV@pKi%2QO3)sV0Q zm}nz*Nr7HAZTCqLA%~kj)H8TIy1l-6vJ8#+-UsB#kuifdU^jj`-;15R5ruWj(TXxp zq?H(~5N5?{4l4E)dV2hl2cLoLD&6Eg#)bNnpldZ`8?pm!y}2WCPc91RrdMZne1D}7 zO#o6DUvjcwZP11*pWq-Y7d(RcJ&H0oD1xzm+PWU6P>fPLDaWA&*t0EJ;A&`CCyk7l zl>>f5iLE;k>%c(|Vs)E3br5cbtvccj0n{SBu?&;%&L2rpc9+^P5Wo8q$E6r8-U^J} zYbvE}Aof+yRB^VPl*t$|ml6L}0XC=4f`2QnapXN1SN`f&)G5IovAevyQbuIQ2B|9* zW102UG=#jiQ|(QYskaPvmm=FZl2vBB@F}sB@lhi8I+ruiY5XHYL7!{3K8o&lDRbs6 z_~?dFxW<{2*{N9!V6RYwCfu1?CbQB8LeJ?ylJyvU;4^ zi7E=J6y_2FFJ-J7(%0fI1ndCCJgAgyPNBh^VUbcVqdeDttaZ%O1%$PTngA$$=o(?_C3(EC^*cJ!s&Goaup8XeKfbA(%JSn7l z`q5UqG;Mp1408K0v6KIFX&sN-E3TTwzW2ki|iPtCA;(L32;a65;3@ff0 zV2B67Z%If01!~_^oKHjQ(_TW>P2NeEF(6t%5UHFCh*%m3Ph54vNrG>|zQy2MU!Wgurg20gz9f=G_ z4H;y)aAs;KMWuHh3f9f>@46GIh?IU~=2!|BlxbuP>~)HsA`);gj`#O$ZDd|fn<%Yh z+Cn2#L(UYN*k|Vq5%2u@jMzk9!7E~fMG?RFI2E~~`fLgYAp)z8yAj90<6146F6sMf z)#zN*EZVrIKbTE{-c_7YO7Lu@?=CdoidzNW%G#C^g*2QS@>x`;{n?3ZwREEYY><%p{FAswM0W^(!(6>l)ZqJFhaC`904`LtM2< z>$8O0i`(i$)db$-Rq25^r-u4QM&<@aDB6M8A#~eAgdY2W;KwZ+3Ku;h3C-?k7QzWZmP<$A!p*-YeRq9ENAj@# zKv>uJnSe?y@~kY1Ux}2M)0;4bHbgQ8 zZhb}EVIp;WmADFE0*Z#{ens<+ns`r>nuba!4ENYTqJgdLuqNNKAS~nKyAmutokw8B zb)dF6j8XyOo2`YSWQ3KXAd^_{wTCd?z`0)}Y|)MI#`jEd_ch(H;g8u`*OHpi^pyCO z$jt@LAcbusL5#_Wl5HePgtSIkOk}e~1*~E5Jox-k1un(Hlee&vX~$1K)jXvv{lwCW z140#sRThj&X{pF@mh#4~u?h!e_YD$$EEm$U=J)iKkOE)Hl$Cls7F^N>eqY`VMFoa? z0iz9zf$QtP6kx}rjAm#}Tf3=*g1eF}Aad~R2g%RVJV^U(6EN*jU|}O@Cs3mAJkft^ zG=wKm_T!wt+SZSh&NWzfSI;HFek%t$4`Wfi+t#n_z@-Nmy(2(Za?QcoAAqlkI1VbYZ|{Pl!{~?yvZ?R>!EM`Q7oNOQAA`~*^^UC%p|;q% z=zVVFm*b1(S>l<$+X^H7l>g=V#A*Wd;DsOB8#^HKtgijTQ{*j-U*y*IVutTaR>#4;Y z$=t*`Ayv@w7DHZsWIv6j(a;#{fLTgb%YR~sG{?8VEne# zG2m*?-=>?X2%o>On9faxnw;U4awD>juDMRD4X$+W>h6ddv&O_9;N8bVW^9!H>wS)N z_)&A~JI>)|ANYCj?{wC9C63yQLDMx-%dj+{Vms-t3|Q)X)+C#T6d0{*?m$Y6)LTWZ zW77)<_~xd-KqpO5VT8JV3e>>ByQi%$ds175b@Wf3kwB|i;IZZdD`C8MmdqsfD(U`c zW!kUA(}2(nNAD1GcZW~l}^_-LLqAis8mqby-m`oOw75)ZQHt&nY_5c zTp8soH#M)~GC)|DI%8V9$bV|a8Pif0uyf=SLzj3SZDl?`#3^6C3Eh{(mhUCp9f(j> zTcIc#Stji9J8#KaN$_UDW|F)?VzREJ#^9sba08oVt|nL}k4vG5CN+ErFFykl<80-pmT+u+@hWRVUv ziZXhv>Oz(Rtz0GydD0=CO`#*mMs~7{nIipDT>R8Sj+l)(R&;+YCqDMA3gaA^m=bId|K>1nfd?8-hmCXy45iwN?!X&oiJR&C|j#SneQ z!MUm!QBCn=j1%uzJME6Zs<(7@tz2f;Vzwg2wlN!(ZPmMHl*t3OE4+FNfx8H=(W^>z z&QIZ$V(7^{_C7(gL8)j=Yr7}aAHNRnHDN=zA>3VP&$_10y5dBrv2i==<7BC z_`3WA-1i`;88(QzdM?wiw&@;Ppi=NDcM|@BceZ>+(#O0=l!Iu7PsWhX z3ZtCNZ&iJvr2c8idZ(28Z3E}*Jis-T7Z>{xKa4Z{#Fy0jT%)6sBq3tq(BwfSpeoUx zzQ6D7M7Bp(krPFl-8WQm!ndSfjOvLL=jWEAb05^t>GJ0d`pd~OY&HH?cm~Cxir!`B z1X4vBHnuEqEHzq37>al=1Cn0-Yn+AM0qjOT|1eWgmReTNa!3 zN)(o#zi2wszI-PR3USoKiU%*)yF#;G=9M&bc>Gg^_ecYYoF*jDw1BX`^9U7`9xjTlcTzBnZb;V4pjY=1=iTJGD@GZETWJiuK zMkR-_DojMjSunDZuYAM)Z@Wq*Z|wW)X@kg8prc}m_@zM4O64lSx{;(ysIvB%-G>+> z-~uf~k1i8&xz7hm3$PTa#Z|FNmUd?x7rEq=xdFWb}z0@@& zX*&Y#OqkNFNcj=O<@(h`REqNTt6~pq2*JV+AH3#3#{}L#v(s`BOCaNFzp3?CxbLNW z^7J(=PgLS<>XP|A;p0c58DUEKkN3LU_6*}rjL`?|P<2bOXz!{!LepPv2@k@LHNnuX z5JS!}9#_U9I@t#J6ntX%Sht(iLU5|w_=INvK|5HDoqH)tv0J!gJ9Rry>Wv=9xE!1M2UW;k4vMP|;4G-Ds7 zM;9evS9WW;9tV&~}QX50%ei5`x3? zBWm|qtDH3<;VE+-r^8}Y-XW~TCUIjoA^|!fzJ#I~dZE zh^ewBU>2_L82`I_Wne*ZyP_V<(Uj5V7PzbOFLfla^scESTABi0&65e>j206^Uv4@<a38z)3HnWvKzDQUBGpP)6iYF61RPJX}V zP~I~jeU6?!-f|L6JHxUNe}%1+{2jkHl0sQ}lDD303cm!j_FlxXR*NydOloQ7BX8L7 zY0K9<#qDs?LH?BSYcp>L(ig02Io9jnv3lG;S(eGaVf#N>)?Z=OzhV3TCd>N!5flLU of74X`zg}|x&B6As9GL$p=bsks@9h77FX!*I0C}%FTQcB(0m-$hS^xk5 literal 0 HcmV?d00001 diff --git a/scripts/pulse/fit-raw.py b/scripts/pulse/fit-raw.py index 0182d78..731fc77 100755 --- a/scripts/pulse/fit-raw.py +++ b/scripts/pulse/fit-raw.py @@ -18,8 +18,12 @@ def gauss_func(x, mean, sig, ampl): f = sys.argv[1] dirName = sys.argv[2] -offset = float(sys.argv[3]) -slopex = float(sys.argv[4]) +# offset = float(sys.argv[3]) +# slopex = float(sys.argv[4]) + +offset = 10 +slopex = 2 + os.chdir(dirName) @@ -43,31 +47,28 @@ verbosityPercentage = 0.01 verbosityPercentageIncrement = 0.01 -print("Status: 0 \%") +for i in range(len(df.loc[:])): + hist = df.loc[i][2:257] -for i in range(len(df.loc[:])): - - if sum(df.loc[i][0:255]) <= maxDead: + if sum(hist) <= maxDead: arrDead.append(i) continue - elif sum(df.loc[i][0:255]) >= 200*maxCounts: + elif sum(hist) >= 200*maxCounts: arrNoisy.append(i) continue - - xS = [] - for j in range(255): - if df.loc[i][j] != 0: - xS.append(j) else: - # hist = df.loc[i][0:255] - params = [] try: - params, tmp = curve_fit(error_func, x, df.loc[i][0:255], guess) # guess before hist - # print(params) + params, tmp = curve_fit( + error_func, + x, + hist, + guess + ) + except RuntimeError: continue errorcounter += 1 @@ -79,63 +80,63 @@ for i in range(len(df.loc[:])): if mean < 255 and mean >= 0: fpn[mean] += 1 - # plt.bar(x,hist,width=1) - # plt.plot(x,error_func(x,params[0],params[1],params[2]),color="red") - # plt.savefig("raw-fit-" + str(i) + ".png") - # plt.clf() - # plt.cla() - # plt.close() - - if (i / len(df.loc[:])) > verbosityPercentage: - print("\rStatus:" + str(verbosityPercentage*100) +"\%") - verbosityPercentage += verbosityPercentageIncrement + plt.bar(x,hist,width=1) + plt.plot(x,error_func(x,params[0],params[1],params[2]),color="red") + plt.savefig("raw-fit-" + str(i) + ".png") + plt.clf() + plt.cla() + plt.close() + + # if (i / len(df.loc[:])) > verbosityPercentage: + # print("\rStatus:" + str(verbosityPercentage*100) +"\%") + # verbosityPercentage += verbosityPercentageIncrement + +# fpnMeanAndSigma = np.zeros(255) -fpnMeanAndSigma = np.zeros(255) +# for i in range(0, 255): -for i in irange(0,255): +# fpnMeanAndSigma[i] = i * fpn[i] +# calcmeansum += i * fpn[i] +# calcmeanheight += fpn[i] - fpnMeanAndSigma[i] = i * fpn[i] - calcmeansum += i * fpn[i] - calcmeanheight += fpn[i] +# stddev = 0 +# means = calcmeansum / calcmeanheight -stddev = 0 -means = calcmeansum / calcmeanheight +# for i in range(0,255): +# stddev += (i - means)**2 * (fpn[i] / calcmeanheight) -for i in range(0,255): - stddev += (i - means)**2 * (fpn[i] / calcmeanheight) +# means = offset + slopex * (calcmeansum / calcmeanheight) +# stddev = slopex * np.sqrt(stddev) -means = offset + slopex * (calcmeansum / calcmeanheight) -stddev = slopex * np.sqrt(stddev) +# # print(means) +# # print(stddev) -# print(means) -# print(stddev) +# guess = [fpn.argmax(), 10, max(fpn)] +# params = [] +# try: +# params, tmp = curve_fit(gauss_func, x, fpn, guess) -guess = [fpn.argmax(), 10, max(fpn)] -params = [] -try: - params, tmp = curve_fit(gauss_func, x, fpn, guess) - -except RuntimeError: - errorcounter += 1 - pass -except RuntimeWarning: - pass +# except RuntimeError: +# errorcounter += 1 +# pass +# except RuntimeWarning: +# pass -# print("\n\nFPNFIT: " + str(params[1]) + "\n\n") -# print("FIT:" + params + "\n\n") +# # print("\n\nFPNFIT: " + str(params[1]) + "\n\n") +# # print("FIT:" + params + "\n\n") -# print("\n\nFPNCALC: " + str(np.std(fpn)) + "\n\n") -# print("CALC:" + str(np.mean(fpn)) + " " + str(np.std(fpn)) + "\n\n") -# print() -# print() +# # print("\n\nFPNCALC: " + str(np.std(fpn)) + "\n\n") +# # print("CALC:" + str(np.mean(fpn)) + " " + str(np.std(fpn)) + "\n\n") +# # print() +# # print() -stringPlot = "GMDT: " + str(means) + "\nFPN: " + str(stddev) + "\nERRS: " + str(errorcounter) -plt.bar(x,fpn,width=1) -plt.plot(x,gauss_func(x,means/slopex - offset,stddev/slopex,max(fpn)),color="red") # plt.plot(x,gauss_func(x,params[0],params[1],params[2]),color="red") -plt.text(0,0,stringPlot) -plt.savefig("fpn-" + f + ".png") +# stringPlot = "GMDT: " + str(means) + "\nFPN: " + str(stddev) + "\nERRS: " + str(errorcounter) +# plt.bar(x,fpn,width=1) +# plt.plot(x,gauss_func(x,means/slopex - offset,stddev/slopex,max(fpn)),color="red") # plt.plot(x,gauss_func(x,params[0],params[1],params[2]),color="red") +# plt.text(0,0,stringPlot) +# plt.savefig("fpn-" + f + ".png") diff --git a/scripts/pulse/plot-raw.py b/scripts/pulse/plot-raw.py index 82fc7e5..016ef59 100755 --- a/scripts/pulse/plot-raw.py +++ b/scripts/pulse/plot-raw.py @@ -11,10 +11,10 @@ df = pd.read_csv(f, delimiter='\t',header=None) x = np.linspace(0, 255, num=255) for i in range(len(df.loc[:])): - if sum(df.loc[i][0:255]) != 0: - hist = df.loc[i][0:255] + if sum(df.loc[i][2:257]) != 0: + hist = df.loc[i][2:257] plt.bar(x,hist,width=1) - plt.savefig("fpn-" + str(i) + ".png") + plt.savefig("fpn-" + str(int(df.loc[i][0])) + "-" + str(int(df.loc[i][1])) + ".png") plt.clf() plt.cla() plt.close() -- 2.43.0