Skip to content

Fixing problems with TRestDetectorSignal fitting methods #119

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 16 commits into from
Jun 4, 2024
Merged

Fixing problems with TRestDetectorSignal fitting methods #119

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
f77da54
first commit
mariajmz May 30, 2024
1075cf3
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] May 30, 2024
f910715
gaussfit
mariajmz May 30, 2024
cb7136c
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] May 31, 2024
1798f30
fixing errors
mariajmz May 31, 2024
693b74a
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] May 31, 2024
eeb0276
delete std::cout
mariajmz Jun 3, 2024
bbec7c3
landau fit
mariajmz Jun 3, 2024
e25384c
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Jun 3, 2024
5829873
aget fit
mariajmz Jun 3, 2024
cff964f
Merge branch 'mariajmz_fits' of github.com:rest-for-physics/detectorl…
mariajmz Jun 3, 2024
818533a
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Jun 3, 2024
12be284
dynamic fit limits for landau and aget
mariajmz Jun 4, 2024
06826f3
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Jun 4, 2024
ef37ea2
code review
lobis Jun 4, 2024
ea80a78
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Jun 4, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions inc/TRestDetectorSignal.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -138,7 +138,7 @@ class TRestDetectorSignal {
void ExponentialConvolution(Double_t fromTime, Double_t decayTime, Double_t offset = 0);
void SignalAddition(TRestDetectorSignal* inSgnl);

Bool_t isSorted();
Bool_t isSorted() const;
void Sort();

void GetDifferentialSignal(TRestDetectorSignal* diffSgnl, Int_t smearPoints = 5);
Expand All @@ -157,7 +157,7 @@ class TRestDetectorSignal {
fSignalCharge.clear();
}

void WriteSignalToTextFile(const TString& filename);
void WriteSignalToTextFile(const TString& filename) const;
void Print() const;

TGraph* GetGraph(Int_t color = 1);
Expand Down
167 changes: 102 additions & 65 deletions src/TRestDetectorSignal.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -291,31 +291,48 @@ TRestDetectorSignal::GetMaxGauss() // returns a 2vector with the time of the pe
Double_t maxRawTime =
GetTime(maxRaw); // The time of the bin where the maximum of the raw signal is found
Double_t energy = 0, time = 0;
Double_t lowerLimit = maxRawTime - 0.2; // us
Double_t upperLimit = maxRawTime + 0.4; // us

TF1* gaus = new TF1("gaus", "gaus", lowerLimit, upperLimit);
TH1F* h1 = new TH1F("h1", "h1", 1000, 0,
10); // Histogram to store the signal. For now the number of bins is fixed.
// Define fit limits
Double_t threshold = GetData(maxRaw) * 0.9; // 90% of the maximum value

Double_t lowerLimit = maxRawTime, upperLimit = maxRawTime;

// Find the lower limit: time when signal drops to 90% of the max before the peak
for (int i = maxRaw; i >= 0; --i) {
if (GetData(i) <= threshold) {
lowerLimit = GetTime(i);
break;
}
}

// Find the upper limit: time when signal drops to 90% of the max after the peak
for (int i = maxRaw; i < GetNumberOfPoints(); ++i) {
if (GetData(i) <= threshold) {
upperLimit = GetTime(i);
break;
}
}

TF1 gaus("gaus", "gaus", lowerLimit, upperLimit);
TH1F h("h", "h", GetNumberOfPoints(), GetTime(0), GetTime(GetNumberOfPoints() - 1));

// copying the signal peak to a histogram
for (int i = 0; i < GetNumberOfPoints(); i++) {
h1->Fill(GetTime(i), GetData(i));
h.SetBinContent(i + 1, GetData(i));
}
/*
TCanvas* c = new TCanvas("c", "Signal fit", 200, 10, 1280, 720);
h1->GetXaxis()->SetTitle("Time (us)");
h1->GetYaxis()->SetTitle("Amplitude");
h1->Draw();
h->GetXaxis()->SetTitle("Time (us)");
h->GetYaxis()->SetTitle("Amplitude");
h->Draw();
*/

TFitResultPtr fitResult =
h1->Fit(gaus, "QNRS"); // Q = quiet, no info in screen; N = no plot; R = fit in the function range; S
// = save and return the fit result
TFitResultPtr fitResult = h.Fit(&gaus, "QNRS"); // Q = quiet, no info in screen; N = no plot; R = fit in
// the function range; S = save and return the fit result

if (fitResult->IsValid()) {
energy = gaus->GetParameter(0);
time = gaus->GetParameter(1);
energy = gaus.GetParameter(0);
time = gaus.GetParameter(1);
} else {
// the fit failed, return -1 to indicate failure
energy = -1;
Copy link
Member

@lobis lobis May 31, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We can take this opportunity to rethink if we want to return energy -1 or do something else.

Copy link
Member

@lobis lobis May 31, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I made #120 to address this.

Expand All @@ -324,24 +341,19 @@ TRestDetectorSignal::GetMaxGauss() // returns a 2vector with the time of the pe
<< "WARNING: bad fit to signal with ID " << GetID() << " with maximum at time = " << maxRawTime
<< " ns "
<< "\n"
<< "Failed fit parameters = " << gaus->GetParameter(0) << " || " << gaus->GetParameter(1)
<< " || " << gaus->GetParameter(2) << "\n"
<< "Failed fit parameters = " << gaus.GetParameter(0) << " || " << gaus.GetParameter(1) << " || "
<< gaus.GetParameter(2) << "\n"
<< "Assigned fit parameters : energy = " << energy << ", time = " << time << endl;
/*
TCanvas* c2 = new TCanvas("c2", "Signal fit", 200, 10, 1280, 720);
h1->Draw();
h->Draw();
c2->Update();
getchar();
delete c2;
*/
}

TVector2 fitParam(time, energy);

delete h1;
delete gaus;

return fitParam;
return {time, energy};
}

// z position by landau fit
Expand All @@ -353,24 +365,42 @@ TRestDetectorSignal::GetMaxLandau() // returns a 2vector with the time of the p
Double_t maxRawTime =
GetTime(maxRaw); // The time of the bin where the maximum of the raw signal is found
Double_t energy = 0, time = 0;
Double_t lowerLimit = maxRawTime - 0.2; // us
Double_t upperLimit = maxRawTime + 0.4; // us

TF1* landau = new TF1("landau", "landau", lowerLimit, upperLimit);
TH1F* h1 = new TH1F("h1", "h1", 1000, 0,
10); // Histogram to store the signal. For now the number of bins is fixed.
// Define fit limits
Double_t threshold = GetData(maxRaw) * 0.9; // 90% of the maximum value

Double_t lowerLimit = maxRawTime, upperLimit = maxRawTime;

// Find the lower limit: time when signal drops to 90% of the max before the peak
for (int i = maxRaw; i >= 0; --i) {
if (GetData(i) <= threshold) {
lowerLimit = GetTime(i);
break;
}
}

// Find the upper limit: time when signal drops to 90% of the max after the peak
for (int i = maxRaw; i < GetNumberOfPoints(); ++i) {
if (GetData(i) <= threshold) {
upperLimit = GetTime(i);
break;
}
}

TF1 landau("landau", "landau", lowerLimit, upperLimit);
TH1F h("h", "h", GetNumberOfPoints(), GetTime(0), GetTime(GetNumberOfPoints() - 1));

// copying the signal peak to a histogram
for (int i = 0; i < GetNumberOfPoints(); i++) {
h1->Fill(GetTime(i), GetData(i));
h.SetBinContent(i + 1, GetData(i));
}

TFitResultPtr fitResult =
h1->Fit(landau, "QNRS"); // Q = quiet, no info in screen; N = no plot; R = fit in the function range;
// S = save and return the fit result
h.Fit(&landau, "QNRS"); // Q = quiet, no info in screen; N = no plot; R = fit in the function range;
// S = save and return the fit result
if (fitResult->IsValid()) {
energy = landau->GetParameter(0);
time = landau->GetParameter(1);
energy = landau.GetParameter(0);
time = landau.GetParameter(1);
} else {
// the fit failed, return -1 to indicate failure
energy = -1;
Expand All @@ -379,24 +409,19 @@ TRestDetectorSignal::GetMaxLandau() // returns a 2vector with the time of the p
<< "WARNING: bad fit to signal with ID " << GetID() << " with maximum at time = " << maxRawTime
<< " ns "
<< "\n"
<< "Failed fit parameters = " << landau->GetParameter(0) << " || " << landau->GetParameter(1)
<< " || " << landau->GetParameter(2) << "\n"
<< "Failed fit parameters = " << landau.GetParameter(0) << " || " << landau.GetParameter(1)
<< " || " << landau.GetParameter(2) << "\n"
<< "Assigned fit parameters : energy = " << energy << ", time = " << time << endl;
/*
TCanvas* c2 = new TCanvas("c2", "Signal fit", 200, 10, 1280, 720);
h1->Draw();
h->Draw();
c2->Update();
getchar();
delete c2;
*/
}

TVector2 fitParam(time, energy);

delete h1;
delete landau;

return fitParam;
return {time, energy};
}

// z position by aget fit
Expand All @@ -420,27 +445,43 @@ TRestDetectorSignal::GetMaxAget() // returns a 2vector with the time of the pea
Double_t maxRawTime =
GetTime(maxRaw); // The time of the bin where the maximum of the raw signal is found
Double_t energy = 0, time = 0;
// The intervals below are small because otherwise the function doesn't fit anymore.
Double_t lowerLimit = maxRawTime - 0.2; // us
Double_t upperLimit = maxRawTime + 0.7; // us

TF1* aget = new TF1("aget", agetResponseFunction, lowerLimit, upperLimit, 3); //
TH1F* h1 = new TH1F("h1", "h1", 1000, 0,
10); // Histogram to store the signal. For now the number of bins is fixed.
aget->SetParameters(500, maxRawTime, 1.2);
// Define fit limits
Double_t threshold = GetData(maxRaw) * 0.9; // 90% of the maximum value

Double_t lowerLimit = maxRawTime, upperLimit = maxRawTime;

// Find the lower limit: time when signal drops to 90% of the max before the peak
for (int i = maxRaw; i >= 0; --i) {
if (GetData(i) <= threshold) {
lowerLimit = GetTime(i);
break;
}
}

// Find the upper limit: time when signal drops to 90% of the max after the peak
for (int i = maxRaw; i < GetNumberOfPoints(); ++i) {
if (GetData(i) <= threshold) {
upperLimit = GetTime(i);
break;
}
}

TF1 aget("aget", agetResponseFunction, lowerLimit, upperLimit, 3); //
TH1F h("h", "h", GetNumberOfPoints(), GetTime(0), GetTime(GetNumberOfPoints() - 1));
aget.SetParameters(500, maxRawTime, 1.2);

// copying the signal peak to a histogram
for (int i = 0; i < GetNumberOfPoints(); i++) {
h1->Fill(GetTime(i), GetData(i));
h.SetBinContent(i + 1, GetData(i));
}

TFitResultPtr fitResult =
h1->Fit(aget, "QNRS"); // Q = quiet, no info in screen; N = no plot; R = fit in
// the function range; S = save and return the fit result
TFitResultPtr fitResult = h.Fit(&aget, "QNRS"); // Q = quiet, no info in screen; N = no plot; R = fit in
// the function range; S = save and return the fit result

if (fitResult->IsValid()) {
energy = aget->GetParameter(0);
time = aget->GetParameter(1);
energy = aget.GetParameter(0);
time = aget.GetParameter(1);
} else {
// the fit failed, return -1 to indicate failure
energy = -1;
Expand All @@ -449,18 +490,14 @@ TRestDetectorSignal::GetMaxAget() // returns a 2vector with the time of the pea
<< "WARNING: bad fit to signal with ID " << GetID() << " with maximum at time = " << maxRawTime
<< " ns "
<< "\n"
<< "Failed fit parameters = " << aget->GetParameter(0) << " || " << aget->GetParameter(1)
<< " || " << aget->GetParameter(2) << "\n"
<< "Failed fit parameters = " << aget.GetParameter(0) << " || " << aget.GetParameter(1) << " || "
<< aget.GetParameter(2) << "\n"
<< "Assigned fit parameters : energy = " << energy << ", time = " << time << endl;
}

TVector2 fitParam(time, energy);

delete h1;
delete aget;

return fitParam;
return {time, energy};
}

Double_t TRestDetectorSignal::GetMaxPeakTime(Int_t from, Int_t to) { return GetTime(GetMaxIndex(from, to)); }

Double_t TRestDetectorSignal::GetMinPeakValue() { return GetData(GetMinIndex()); }
Expand Down Expand Up @@ -518,7 +555,7 @@ Int_t TRestDetectorSignal::GetTimeIndex(Double_t t) {
return -1;
}

Bool_t TRestDetectorSignal::isSorted() {
Bool_t TRestDetectorSignal::isSorted() const {
for (int i = 0; i < GetNumberOfPoints() - 1; i++) {
if (GetTime(i + 1) < GetTime(i)) {
return false;
Expand Down Expand Up @@ -718,7 +755,7 @@ void TRestDetectorSignal::GetSignalGaussianConvolution(TRestDetectorSignal* conv
cout << "Final charge of the pulse " << totChargeFinal << endl;
}

void TRestDetectorSignal::WriteSignalToTextFile(const TString& filename) {
void TRestDetectorSignal::WriteSignalToTextFile(const TString& filename) const {
FILE* fff = fopen(filename.Data(), "w");
for (int i = 0; i < GetNumberOfPoints(); i++) {
fprintf(fff, "%e\t%e\n", GetTime(i), GetData(i));
Expand Down