From: Maps Date: Thu, 12 Dec 2024 10:36:56 +0000 (+0100) Subject: pulsing stuff X-Git-Url: https://jspc29.x-matter.uni-frankfurt.de/git/?a=commitdiff_plain;h=d9cb151ce8022f4c412cf0d79cf637ae3f30e4c9;p=mimosis_chain.git pulsing stuff --- 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 <= 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()