Skip to content

Commit

Permalink
Merge pull request #119 from rest-for-physics/mariajmz_fits
Browse files Browse the repository at this point in the history 
Fixing problems with TRestDetectorSignal fitting methods
  • Loading branch information
@lobis
lobis authored Jun 4, 2024
2 parents e47e9a5 + ea80a78 commit 9df8c2f
Show file tree
Hide file tree
Showing 2 changed files with 104 additions and 67 deletions.
4 changes: 2 additions & 2 deletions inc/TRestDetectorSignal.h
View file
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
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;
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

0 comments on commit 9df8c2f

Please sign in to comment.