Skip to content
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

Jump to bottom

Fixing problems with TRestDetectorSignal fitting methods #119

Merged
merged 16 commits into from
Jun 4, 2024
Merged

Fixing problems with TRestDetectorSignal fitting methods #119

Changes from 12 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
85 changes: 46 additions & 39 deletions src/TRestDetectorSignal.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -291,16 +291,34 @@ 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 where 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 where 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 h1("h1", "h1", GetNumberOfPoints(), GetTime(0), GetTime(GetNumberOfPoints() - 1));
Copy link
Contributor Author

Choose a reason for hiding this comment

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

We are maintaining the fit with the signal stored as histogram but as Juanan commented, it is not necessary as we can get the signal as TGraph with TRestDetectorSignal::GetGraph(Int_t color). What do you think?

Copy link
Member

Choose a reason for hiding this comment

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

I think TGraph is a better choice but on practice I don't think it makes a difference since the histogram has one bin per signal time bin, so the information content of both representations is the same.

You can change them to use TGraph if you want but I would do it in a different PR (unless it's trivial to do so).


// copying the signal peak to a histogram
for (int i = 0; i < GetNumberOfPoints(); i++) {
h1->Fill(GetTime(i), GetData(i));
h1.SetBinContent(i + 1, GetData(i));
}
/*
TCanvas* c = new TCanvas("c", "Signal fit", 200, 10, 1280, 720);
Expand All @@ -310,12 +328,12 @@ TRestDetectorSignal::GetMaxGauss() // returns a 2vector with the time of the pe
*/

TFitResultPtr fitResult =
h1->Fit(gaus, "QNRS"); // Q = quiet, no info in screen; N = no plot; R = fit in the function range; S
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

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

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

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,8 +342,8 @@ 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);
Expand All @@ -338,9 +356,6 @@ TRestDetectorSignal::GetMaxGauss() // returns a 2vector with the time of the pe

TVector2 fitParam(time, energy);

delete h1;
delete gaus;

return fitParam;
}

Expand All @@ -356,21 +371,20 @@ TRestDetectorSignal::GetMaxLandau() // returns a 2vector with the time of the p
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.
TF1 landau("landau", "landau", lowerLimit, upperLimit);
TH1F h1("h1", "h1", 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));
h1.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;
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
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,8 +393,8 @@ 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);
Expand All @@ -393,9 +407,6 @@ TRestDetectorSignal::GetMaxLandau() // returns a 2vector with the time of the p

TVector2 fitParam(time, energy);

delete h1;
delete landau;

return fitParam;
}

Expand Down Expand Up @@ -424,23 +435,22 @@ TRestDetectorSignal::GetMaxAget() // returns a 2vector with the time of the pea
Double_t lowerLimit = maxRawTime - 0.2; // us
Double_t upperLimit = maxRawTime + 0.7; // us
Copy link
Contributor Author

Choose a reason for hiding this comment

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

I haven't changed this limits (also in the landau fit). Maybe we should adjust them too

Copy link
Member

Choose a reason for hiding this comment

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

Yes, I think we should. The new dynamic time window should work better in most cases than what we have now.


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);
TF1 aget("aget", agetResponseFunction, lowerLimit, upperLimit, 3); //
TH1F h1("h1", "h1", 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));
h1.SetBinContent(i + 1, GetData(i));
}

TFitResultPtr fitResult =
h1->Fit(aget, "QNRS"); // Q = quiet, no info in screen; N = no plot; R = fit in
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

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,16 +459,13 @@ 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;
}
Double_t TRestDetectorSignal::GetMaxPeakTime(Int_t from, Int_t to) { return GetTime(GetMaxIndex(from, to)); }
Expand Down
Loading