Skip to content

Commit

Permalink
Adding full module spectrum (new attributes and methods) and use it a…
Browse files Browse the repository at this point in the history 
…s reference for DrawGainMap
  • Loading branch information
@AlvaroEzq
AlvaroEzq committed Mar 16, 2024
1 parent 377d727 commit 85d2563
Show file tree
Hide file tree
Showing 2 changed files with 207 additions and 30 deletions.
21 changes: 20 additions & 1 deletion source/framework/analysis/inc/TRestDataSetGainMap.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,6 +88,8 @@ class TRestDataSetGainMap : public TRestMetadata {
Module* GetModule(const int planeID, const int moduleID);
double GetSlopeParameter(const int planeID, const int moduleID, const double x, const double y);
double GetInterceptParameter(const int planeID, const int moduleID, const double x, const double y);
double GetSlopeParameterFullSpc(const int planeID, const int moduleID);
double GetInterceptParameterFullSpc(const int planeID, const int moduleID);

void SetCalibrationFileName(const std::string& fileName) { fCalibFileName = fileName; }
void SetOutputFileName(const std::string& fileName) { fOutputFileName = fileName; }
Expand Down Expand Up @@ -168,6 +170,12 @@ class TRestDataSetGainMap : public TRestMetadata {
/// Array containing the slope of the linear fit for each segment.
std::vector<std::vector<double>> fSlope = {}; //<

/// Slope of the calibration linear fit of whole module
double fFullSlope = 0; //<

/// Intercept of the calibration linear fit of whole module
double fFullIntercept = 0; //<

/// Array containing the intercept of the linear fit for each segment.
std::vector<std::vector<double>> fIntercept = {}; //<

Expand All @@ -181,9 +189,15 @@ class TRestDataSetGainMap : public TRestMetadata {
/// Array containing the observable spectrum for each segment.
std::vector<std::vector<TH1F*>> fSegSpectra = {}; //<

/// Spectrum of the observable for the whole module.
TH1F* fFullSpectrum = nullptr; //<

/// Array containing the calibration linear fit for each segment.
std::vector<std::vector<TGraph*>> fSegLinearFit = {}; //<

/// Calibration linear fit for the whole module.
TGraph* fFullLinearFit = nullptr; //<

public:
void AddPeak(const double& energyPeak, const TVector2& rangePeak = TVector2(0, 0)) {
fEnergyPeaks.push_back(energyPeak);
Expand All @@ -195,6 +209,8 @@ class TRestDataSetGainMap : public TRestMetadata {
std::pair<int, int> GetIndexMatrix(const double x, const double y) const;
double GetSlope(const double x, const double y) const;
double GetIntercept(const double x, const double y) const;
double GetSlopeFullSpc() const { return fFullSlope; };
double GetInterceptFullSpc() const { return fFullIntercept; };

Int_t GetPlaneId() const { return fPlaneId; }
Int_t GetModuleId() const { return fModuleId; }
Expand All @@ -215,7 +231,7 @@ class TRestDataSetGainMap : public TRestMetadata {
TCanvas* c = nullptr);
void DrawSpectrum(const int index_x, const int index_y, bool drawFits = true, int color = -1,
TCanvas* c = nullptr);
void DrawFullSpectrum();
void DrawFullSpectrum(const bool drawFits = true, const int color = -1, TCanvas* c = nullptr);

void DrawLinearFit(TCanvas* c = nullptr);
void DrawLinearFit(const TVector2& position, TCanvas* c = nullptr);
Expand All @@ -225,7 +241,10 @@ class TRestDataSetGainMap : public TRestMetadata {

void Refit(const TVector2& position, const double energy, const TVector2& range);
void Refit(const size_t x, const size_t y, const size_t peakNumber, const TVector2& range);
void RefitFullSpc(const double energy, const TVector2& range);
void RefitFullSpc(const size_t peakNumber, const TVector2& range);
void UpdateCalibrationFits(const size_t x, const size_t y);
void UpdateCalibrationFitsFullSpc();

void SetPlaneId(const Int_t& planeId) { fPlaneId = planeId; }
void SetModuleId(const Int_t& moduleId) { fModuleId = moduleId; }
Expand Down
216 changes: 187 additions & 29 deletions source/framework/analysis/src/TRestDataSetGainMap.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -300,6 +300,16 @@ double TRestDataSetGainMap::GetSlopeParameter(const int planeID, const int modul
return moduleCal->GetSlope(x, y);
}

/////////////////////////////////////////////
/// \brief Function to get the slope parameter of the whole module with
/// planeID and moduleID.
///
double TRestDataSetGainMap::GetSlopeParameterFullSpc(const int planeID, const int moduleID) {
Module* moduleCal = GetModule(planeID, moduleID);
if (moduleCal == nullptr) return 0; // return numeric_limits<double>::quiet_NaN()
return moduleCal->GetSlopeFullSpc();
}

/////////////////////////////////////////////
/// \brief Function to get the intercept parameter of the module with
/// planeID and moduleID at physical position (x,y)
Expand All @@ -312,6 +322,16 @@ double TRestDataSetGainMap::GetInterceptParameter(const int planeID, const int m
return moduleCal->GetIntercept(x, y);
}

/////////////////////////////////////////////
/// \brief Function to get the intercept parameter of the whole module with
/// planeID and moduleID.
///
double TRestDataSetGainMap::GetInterceptParameterFullSpc(const int planeID, const int moduleID) {
Module* moduleCal = GetModule(planeID, moduleID);
if (moduleCal == nullptr) return 0; // return numeric_limits<double>::quiet_NaN()
return moduleCal->GetInterceptFullSpc();
}

/////////////////////////////////////////////
/// \brief Function to get a list (set) of the plane IDs
///
Expand Down Expand Up @@ -671,12 +691,22 @@ void TRestDataSetGainMap::Module::GenerateGainMap() {
<< p->RESTendl;
}

// --- Definition of histogram whole module ---
std::string hModuleName = "hSpc_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId);
TH1F* hModule = new TH1F(hModuleName.c_str(), "", fNBins, fCalibRange.X(), fCalibRange.Y());

// build the spectrum for the whole module
std::string cut = fDefinitionCut;
if (cut.empty()) cut = "1";
auto histoMod = dataSet.GetDataFrame().Filter(cut).Histo1D(
{"tempMod", "", fNBins, fCalibRange.X(), fCalibRange.Y()}, GetObservable());
hModule = (TH1F*)histoMod->Clone(hModuleName.c_str());

//--- Definition of histogram matrix ---
std::vector<std::vector<TH1F*>> h(fNumberOfSegmentsX, std::vector<TH1F*>(fNumberOfSegmentsY, nullptr));
for (size_t i = 0; i < h.size(); i++) {
for (size_t j = 0; j < h.at(0).size(); j++) {
std::string name = "hSpc_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId) + "_" +
std::to_string(i) + "_" + std::to_string(j);
std::string name = hModuleName + std::to_string(i) + "_" + std::to_string(j);
h[i][j] = new TH1F(name.c_str(), "", fNBins, fCalibRange.X(),
fCalibRange.Y()); // h[column][row] equivalent to h[x][y]
}
Expand Down Expand Up @@ -723,12 +753,23 @@ void TRestDataSetGainMap::Module::GenerateGainMap() {
//--- Fit every peak energy for every segment ---
fSegLinearFit = std::vector(h.size(), std::vector<TGraph*>(h.at(0).size(), nullptr));
for (size_t i = 0; i < h.size(); i++) {
for (size_t j = 0; j < h.at(0).size(); j++) {
RESTExtreme << "Segment[" << i << "][" << j << "]" << p->RESTendl;
for (int j = -1; j < (int)h.at(0).size(); j++) // -1 for the whole module
{
TH1F* hSeg = nullptr;
if (i > 0 && j == -1)
continue;
else if (i == 0 && j == -1) {
hSeg = hModule;
RESTExtreme << "Whole module" << p->RESTendl;
} else {
hSeg = h[i][j];
RESTExtreme << "Segment[" << i << "][" << j << "]" << p->RESTendl;
}

// Search for peaks --> peakPos
std::unique_ptr<TSpectrum> s(new TSpectrum(2 * fEnergyPeaks.size() + 1));
std::vector<double> peakPos;
s->Search(h[i][j], 2, "goff", 0.1);
s->Search(hSeg, 2, "goff", 0.1);
for (int k = 0; k < s->GetNPeaks(); k++) peakPos.push_back(s->GetPositionX()[k]);
std::sort(peakPos.begin(), peakPos.end(), std::greater<double>());
const double ratio = peakPos.size() == 0
Expand Down Expand Up @@ -790,10 +831,10 @@ void TRestDataSetGainMap::Module::GenerateGainMap() {
<< DoubleToString(start, "%.3g") << ", " << DoubleToString(end, "%.3g") << ")"
<< p->RESTendl;

if (h[i][j]->GetFunction(name.c_str())) // remove previous fit
h[i][j]->GetListOfFunctions()->Remove(h[i][j]->GetFunction(name.c_str()));
if (hSeg->GetFunction(name.c_str())) // remove previous fit
hSeg->GetListOfFunctions()->Remove(hSeg->GetFunction(name.c_str()));

h[i][j]->Fit(g, "R+Q0"); // use 0 to not draw the fit but save it
hSeg->Fit(g, "R+Q0"); // use 0 to not draw the fit but save it
mu = g->GetParameter(1);
RESTExtreme << "\t\tgaus mean " << DoubleToString(mu, "%g") << p->RESTendl;
} while (fAutoRangePeaks && peakPos.size() > 0 &&
Expand All @@ -814,16 +855,22 @@ void TRestDataSetGainMap::Module::GenerateGainMap() {
std::unique_ptr<TF1> linearFit;
linearFit = std::unique_ptr<TF1>(new TF1("linearFit", "pol1"));
gr->Fit("linearFit", "SQ"); // Q for quiet mode

fSegLinearFit.at(i).at(j) = (TGraph*)gr->Clone();

const double slope = linearFit->GetParameter(1);
const double intercept = linearFit->GetParameter(0);
calParSlope.at(i).at(j) = slope;
calParIntercept.at(i).at(j) = intercept;

if (j == -1) {
fFullLinearFit = (TGraph*)gr->Clone();
fFullSlope = slope;
fFullIntercept = intercept;
} else {
fSegLinearFit.at(i).at(j) = (TGraph*)gr->Clone();
calParSlope.at(i).at(j) = slope;
calParIntercept.at(i).at(j) = intercept;
}
}
}
fSegSpectra = h;
fFullSpectrum = hModule;
fSlope = calParSlope;
fIntercept = calParIntercept;
}
Expand Down Expand Up @@ -887,6 +934,55 @@ void TRestDataSetGainMap::Module::Refit(const size_t x, const size_t y, const si
UpdateCalibrationFits(x, y);
}

/////////////////////////////////////////////
/// \brief Function to fit again manually a peak for the whole module spectrum. The
/// calibration curve is updated after the fit.
///
/// \param energyPeak The energy of the peak to be fitted (in physical units).
/// \param range The range for the fitting of the peak (in the observables corresponding units).
///
void TRestDataSetGainMap::Module::RefitFullSpc(const double energyPeak, const TVector2& range) {
int peakNumber = -1;
for (size_t i = 0; i < fEnergyPeaks.size(); i++)
if (fEnergyPeaks.at(i) == energyPeak) {
peakNumber = i;
break;
}
if (peakNumber == -1) {
RESTError << "Energy " << energyPeak << " not found in the list of energy peaks" << p->RESTendl;
return;
}
RefitFullSpc((size_t)peakNumber, range);
}

/////////////////////////////////////////////
/// \brief Function to fit again manually a peak for the whole module spectrum. The
/// calibration curve is updated after the fit.
///
/// \param peakNumber The index of the peak to be fitted.
/// \param range The range for the fitting of the peak (in the observables corresponding units).
///
void TRestDataSetGainMap::Module::RefitFullSpc(const size_t peakNumber, const TVector2& range) {
if (!fFullSpectrum) {
RESTError << "No gain map found. Use GenerateGainMap() first." << p->RESTendl;
return;
}
if (peakNumber >= fEnergyPeaks.size()) {
RESTError << "Peak with index " << peakNumber << " not found" << p->RESTendl;
return;
}

// Refit the desired peak
std::string name = "g" + std::to_string(peakNumber);
TF1* g = new TF1(name.c_str(), "gaus", range.X(), range.Y());
while (fFullSpectrum->GetFunction(name.c_str())) // clear previous fits for this peakNumber
fFullSpectrum->GetListOfFunctions()->Remove(fFullSpectrum->GetFunction(name.c_str()));
fFullSpectrum->Fit(g, "R+Q0"); // use 0 to not draw the fit but save it

// Change the point of the graph
UpdateCalibrationFitsFullSpc();
}

/////////////////////////////////////////////
/// \brief Function to update the calibration curve for a given segment of the module. The calibration
/// curve is cleared and then the means of the gaussian fits for each energy peak are added. If there are
Expand Down Expand Up @@ -941,6 +1037,51 @@ void TRestDataSetGainMap::Module::UpdateCalibrationFits(const size_t x, const si
fIntercept.at(x).at(y) = lf->GetParameter(0);
}

/////////////////////////////////////////////
/// \brief Function to update the calibration curve for the whole module. The calibration
/// curve is cleared and then the means of the gaussian fits for each energy peak are added. If there are
/// less than 2 fits, zero points are added. Then, the calibration curve is refitted (linearFit).
///
void TRestDataSetGainMap::Module::UpdateCalibrationFitsFullSpc() {
if (!fFullSpectrum) {
RESTError << "No gain map found. Use GenerateGainMap() first." << p->RESTendl;
return;
}

TGraph* gr = fFullLinearFit;

// Clear the points of the graph
for (size_t i = 0; i < fEnergyPeaks.size(); i++) gr->RemovePoint(i);
// Add the new points to the graph
int c = 0;
for (size_t i = 0; i < fEnergyPeaks.size(); i++) {
std::string fitName = (std::string) "g" + std::to_string(i);
TF1* g = fFullSpectrum->GetFunction(fitName.c_str());
if (!g) {
RESTWarning << "No fit ( " << fitName << " ) found for energy peak " << fEnergyPeaks[i]
<< " in whole module" << p->RESTendl;
continue;
}
gr->SetPoint(c++, g->GetParameter(1), fEnergyPeaks[i]);
}

// Add zero points if needed (if there are less than 2 points)
while (gr->GetN() < 2) {
gr->SetPoint(c++, 0, 0);
RESTWarning << "Not enough points for linear fit at whole module. Adding zero point." << p->RESTendl;
}

// Refit the calibration curve
TF1* lf = nullptr;
if (gr->GetFunction("linearFit"))
lf = gr->GetFunction("linearFit");
else
lf = new TF1("linearFit", "pol1");
gr->Fit(lf, "SQ"); // Q for quiet mode
fFullSlope = lf->GetParameter(1);
fFullIntercept = lf->GetParameter(0);
}

/////////////////////////////////////////////
/// \brief Function to read the parameters from the RML element (TiXmlElement)
/// and set those class members.
Expand Down Expand Up @@ -1109,25 +1250,35 @@ void TRestDataSetGainMap::Module::DrawSpectrum(const bool drawFits, const int co
}
}
}
void TRestDataSetGainMap::Module::DrawFullSpectrum() {
if (fSegSpectra.size() == 0) {
RESTError << "Spectra matrix is empty." << p->RESTendl;
void TRestDataSetGainMap::Module::DrawFullSpectrum(const bool drawFits, const int color, TCanvas* c) {
if (!fFullSpectrum) {
RESTError << "Spectrum is empty." << p->RESTendl;
return;
}
TH1F* sumHist =
new TH1F("fullSpc", "", fSegSpectra[0][0]->GetNbinsX(), fSegSpectra[0][0]->GetXaxis()->GetXmin(),
fSegSpectra[0][0]->GetXaxis()->GetXmax());

sumHist->SetTitle(("Full spectrum;" + GetObservable() + ";counts").c_str());
for (size_t i = 0; i < fSegSpectra.size(); i++) {
for (size_t j = 0; j < fSegSpectra.at(0).size(); j++) {
sumHist->Add(fSegSpectra[i][j]);
}
if (!c) {
std::string t = "fullSpc_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId);
c = new TCanvas(t.c_str(), t.c_str());
}
std::string t = "fullSpc_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId);
TCanvas* c = new TCanvas(t.c_str(), t.c_str());
c->cd();
sumHist->Draw();

fFullSpectrum->SetTitle(("Full spectrum;" + GetObservable() + ";counts").c_str());

if (color > 0) fFullSpectrum->SetLineColor(color);
size_t colorT = fFullSpectrum->GetLineColor();
fFullSpectrum->Draw("same");

if (drawFits)
for (size_t c = 0; c < fEnergyPeaks.size(); c++) {
auto fit = fFullSpectrum->GetFunction(("g" + std::to_string(c)).c_str());
if (!fit) RESTWarning << "Fit for energy peak" << fEnergyPeaks[c] << " not found." << p->RESTendl;
if (!fit) continue;
fit->SetLineColor(c + 2 != colorT ? c + 2 : ++colorT); /* does not work with kRed, kBlue, etc.
as they are not defined with the same number as the first 10 basic colors. See
https://root.cern.ch/doc/master/classTColor.html#C01 and
https://root.cern.ch/doc/master/classTColor.html#C02 */
fit->Draw("same");
}
}

void TRestDataSetGainMap::Module::DrawLinearFit(const TVector2& position, TCanvas* c) {
Expand Down Expand Up @@ -1200,6 +1351,10 @@ void TRestDataSetGainMap::Module::DrawGainMap(const int peakNumber) {
RESTError << "Linear fit matrix is empty." << p->RESTendl;
return;
}
if (!fFullLinearFit) {
RESTError << "Full linear fit is empty." << p->RESTendl;
return;
}

double peakEnergy = fEnergyPeaks[peakNumber];
std::string title = "Gain map for energy " + DoubleToString(peakEnergy, "%g") + ";" +
Expand All @@ -1211,8 +1366,7 @@ void TRestDataSetGainMap::Module::DrawGainMap(const int peakNumber) {
TH2F* hGainMap = new TH2F(("h" + t).c_str(), title.c_str(), fNumberOfSegmentsX, fReadoutRange.X(),
fReadoutRange.Y(), fNumberOfSegmentsY, fReadoutRange.X(), fReadoutRange.Y());

const double peakPosRef =
fSegLinearFit[(fNumberOfSegmentsX - 1) / 2][(fNumberOfSegmentsY - 1) / 2]->GetPointX(peakNumber);
const double peakPosRef = fFullLinearFit->GetPointX(peakNumber);

auto itX = fSplitX.begin();
for (size_t i = 0; i < fSegLinearFit.size(); i++) {
Expand Down Expand Up @@ -1313,5 +1467,9 @@ void TRestDataSetGainMap::Module::Print() const {
}
RESTMetadata << p->RESTendl;
}

RESTMetadata << " Full slope: " << DoubleToString(fFullSlope, "%.3e") << p->RESTendl;
RESTMetadata << " Full intercept: " << DoubleToString(fFullIntercept, "%+.3e") << p->RESTendl;

RESTMetadata << "-----------------------------------------------" << p->RESTendl;
}

0 comments on commit 85d2563

Please sign in to comment.