140#include "TRestDataSetGainMap.h"
183 TiXmlElement* moduleDefinition =
GetElement(
"module");
184 while (moduleDefinition !=
nullptr) {
186 fModulesCal.back().LoadConfigFromTiXmlElement(moduleDefinition);
201 RESTInfo <<
"Generating gain map of plane " << mod.GetPlaneId() <<
" module " << mod.GetModuleId()
203 mod.GenerateGainMap();
206 mod.DrawFullSpectrum();
227 std::vector<std::string> excludeColumns) {
229 RESTError <<
"TRestDataSetGainMap::CalibrateDataSet: No modules defined." <<
RESTendl;
234 dataSet.EnableMultiThreading(
true);
237 dataSet.
Import(dataSetFileName);
239 RESTWarning << dataSetFileName <<
" is not a dataset. Generating a temporal one..." <<
RESTendl;
241 dataSet.SetFilePattern(dataSetFileName);
242 dataSet.SetObservablesList({
"*"});
249 std::string pmIDname = (std::string)GetName() +
"_pmID";
250 std::string modCut =
fModulesCal[0].GetModuleDefinitionCut();
251 if (modCut.empty()) modCut =
"1";
254 auto columnList = dataFrame.GetColumnNames();
255 if (std::find(columnList.begin(), columnList.end(), pmIDname) == columnList.end())
256 dataFrame = dataFrame.Define(pmIDname, modCut +
" ? " + std::to_string(pmID) +
" : -1");
258 dataFrame = dataFrame.Redefine(pmIDname, modCut +
" ? " + std::to_string(pmID) +
" : -1");
262 if (modCut.empty()) modCut =
"1";
264 dataFrame = dataFrame.Redefine(pmIDname, (modCut +
" ? " + std::to_string(pmID) +
" : " + pmIDname));
268 auto calibrate = [
this](
double val,
double x,
double y,
int pmID) {
270 if (pmID == m.GetPlaneId() * 10 + m.GetModuleId())
271 return m.GetSlope(x, y) * val + m.GetIntercept(x, y);
275 return std::numeric_limits<double>::quiet_NaN();
277 std::string calibObsName = (std::string)GetName() +
"_";
278 calibObsName += GetObservable().erase(0, GetObservable().find(
"_") + 1);
279 dataFrame = dataFrame.Define(calibObsName, calibrate,
283 auto calibrateFullSpc = [
this](
double val,
int pmID) {
285 if (pmID == m.GetPlaneId() * 10 + m.GetModuleId())
286 return m.GetSlopeFullSpc() * val + m.GetInterceptFullSpc();
290 return std::numeric_limits<double>::quiet_NaN();
292 std::string calibObsNameFullSpc = (std::string)GetName() +
"_";
293 calibObsNameFullSpc +=
294 GetObservable().erase(0, GetObservable().find(
"_") + 1);
295 calibObsNameFullSpc +=
"_NoSegmentation";
296 dataFrame = dataFrame.Define(calibObsNameFullSpc, calibrateFullSpc, {
fObservable, pmIDname});
298 dataSet.SetDataFrame(dataFrame);
301 if (outputFileName.empty()) outputFileName = dataSetFileName;
302 if (outputFileName == dataSetFileName) {
303 std::string gmName = GetName();
304 outputFileName = outputFileName.substr(0, outputFileName.find_last_of(
"."));
312 excludeCol.erase(calibObsName);
313 excludeCol.erase(calibObsNameFullSpc);
314 excludeCol.erase(pmIDname);
316 RESTDebug <<
"Excluding columns: ";
317 for (
auto& c : excludeCol) RESTDebug << c <<
", ";
320 dataSet.
Export(outputFileName, std::vector<std::string>(excludeCol.begin(), excludeCol.end()));
323 TFile* f = TFile::Open(outputFileName.c_str(),
"UPDATE");
336 RESTError <<
"No ModuleCalibration with index " << index;
338 RESTError <<
". There are no modules defined." <<
RESTendl;
350 if (i.GetPlaneId() == planeID && i.GetModuleId() == moduleID)
return &i;
352 RESTError <<
"No ModuleCalibration with planeID " << planeID <<
" and moduleID " << moduleID <<
RESTendl;
364 if (moduleCal ==
nullptr)
return 0;
374 if (moduleCal ==
nullptr)
return 0;
375 return moduleCal->GetSlopeFullSpc();
386 if (moduleCal ==
nullptr)
return 0;
396 if (moduleCal ==
nullptr)
return 0;
397 return moduleCal->GetInterceptFullSpc();
405 std::set<int> planeIDs;
406 for (
const auto& mc :
fModulesCal) planeIDs.insert(mc.GetPlaneId());
415 std::set<int> moduleIDs;
417 if (mc.GetPlaneId() == planeId) moduleIDs.insert(mc.GetModuleId());
425 std::map<int, std::set<int>> moduleIds;
427 moduleIds.insert(std::pair<
int, std::set<int>>(planeId,
GetModuleIDs(planeId)));
432 SetName(src.GetName());
450 if (i.GetPlaneId() == moduleCal.GetPlaneId() && i.GetModuleId() == moduleCal.GetModuleId()) {
463 if (fileName.empty()) {
464 RESTError <<
"No input calibration file defined" <<
RESTendl;
469 RESTInfo <<
"Opening " << fileName <<
RESTendl;
470 TFile* f = TFile::Open(fileName.c_str(),
"READ");
472 RESTError <<
"Cannot open calibration file " << fileName <<
RESTendl;
478 TIter nextkey(f->GetListOfKeys());
480 while ((key = (TKey*)nextkey())) {
481 std::string kName = key->GetClassName();
490 RESTError <<
"File extension not supported for " << fileName <<
RESTendl;
503 RESTError <<
"No output file defined" <<
RESTendl;
509 this->
Write(GetName());
524 RESTMetadata <<
" Cuts applied: ";
527 for (
const auto& cut :
fCut->GetCutStrings()) RESTMetadata << cut <<
", " <<
RESTendl;
528 for (
const auto& cut :
fCut->GetParamCut())
529 RESTMetadata << cut.first <<
" " << cut.second <<
", " <<
RESTendl;
533 RESTMetadata <<
" Number of planes: " << GetNumberOfPlanes() <<
RESTendl;
534 RESTMetadata <<
" Number of modules: " << GetNumberOfModules() <<
RESTendl;
538 RESTMetadata <<
"-----------------------------------------------" <<
RESTendl;
540 RESTMetadata <<
"***********************************************" <<
RESTendl;
552 int index_x = -1, index_y = -1;
555 index_x = std::distance(
fSplitX.begin(),
fSplitX.upper_bound(x)) - 1;
557 RESTWarning <<
"index_x < 0 for x = " << x <<
" and fSplitX[0]=" << *
fSplitX.begin()
562 RESTWarning <<
"x is out of split for x = " << x <<
p->
RESTendl;
567 index_y = std::distance(
fSplitY.begin(),
fSplitY.upper_bound(y)) - 1;
569 RESTWarning <<
"index_y < 0 for y = " << y <<
" and fSplitY[0]=" << *
fSplitY.begin()
574 RESTWarning <<
"y is out of split for y = " << y <<
p->
RESTendl;
578 return std::make_pair(index_x, index_y);
588 auto [index_x, index_y] = GetIndexMatrix(x, y);
589 if (fSlope.empty()) {
590 RESTError <<
"Calibration slope matrix is empty. Returning 0" << p->RESTendl;
594 if (index_x > (
int)fSlope.size() || index_y > (
int)fSlope.at(0).size()) {
595 RESTError <<
"Index out of range. Returning 0" << p->RESTendl;
599 return fSlope[index_x][index_y];
610 auto [index_x, index_y] = GetIndexMatrix(x, y);
611 if (fIntercept.empty()) {
612 RESTError <<
"Calibration constant matrix is empty. Returning 0" << p->RESTendl;
616 if (index_x > (
int)fIntercept.size() || index_y > (
int)fIntercept.at(0).size()) {
617 RESTError <<
"Index out of range. Returning 0" << p->RESTendl;
621 return fIntercept[index_x][index_y];
638 if (fNumberOfSegmentsX < 1) {
639 RESTError <<
"SetSplitX: fNumberOfSegmentsX must be >=1." << p->RESTendl;
642 std::set<double> split;
643 for (
int i = 0; i <= fNumberOfSegmentsX; i++) {
645 fReadoutRange.X() + ((fReadoutRange.Y() - fReadoutRange.X()) / (
float)fNumberOfSegmentsX) * i;
652 if (splitX.size() < 2) {
653 RESTError <<
"SetSplitX: split size must be >=2 (start and end of range must be included)."
658 RESTWarning <<
"SetSplitX: changing split but current gain map and calibration paremeters correspond "
659 "to previous splitting. Use GenerateGainMap() to update them."
662 fNumberOfSegmentsX = fSplitX.size() - 1;
671 if (fNumberOfSegmentsY < 1) {
672 RESTError <<
"SetSplitY: fNumberOfSegmentsY must be >=1." << p->RESTendl;
675 std::set<double> split;
676 for (
int i = 0; i <= fNumberOfSegmentsY; i++) {
678 fReadoutRange.X() + ((fReadoutRange.Y() - fReadoutRange.X()) / (
float)fNumberOfSegmentsY) * i;
685 if (splitY.size() < 2) {
686 RESTError <<
"SetSplitY: split size must be >=2 (start and end of range must be included)."
691 RESTWarning <<
"SetSplitY: changing split but current gain map and calibration paremeters correspond "
692 "to previous splitting. Use GenerateGainMap() to update them."
695 fNumberOfSegmentsY = fSplitY.size() - 1;
716 std::string dsFileName = fDataSetFileName;
717 if (dsFileName.empty()) dsFileName = p->GetCalibrationFileName();
718 if (dsFileName.empty()) {
719 RESTError <<
"No calibration file defined" << p->RESTendl;
724 dataSet.EnableMultiThreading(
true);
726 dataSet.
Import(dsFileName);
727 fDataSetFileName = dsFileName;
729 RESTWarning << dsFileName <<
" is not a dataset. Generating a temporal one..." << p->RESTendl;
731 std::vector<std::string> obsList;
733 obsList.push_back(p->GetObservable());
734 obsList.push_back(p->GetSpatialObservableX());
735 obsList.push_back(p->GetSpatialObservableY());
739 obsList.insert(obsList.end(), modDefCutObs.begin(), modDefCutObs.end());
742 for (
const auto& cut : p->GetCut()->GetCutStrings()) {
744 obsList.insert(obsList.end(), cutObs.begin(), cutObs.end());
746 for (
const auto& [variable, condition] : p->GetCut()->GetParamCut()) {
748 obsList.insert(obsList.end(), cutObs.begin(), cutObs.end());
751 obsList.insert(obsList.end(), cutObs.begin(), cutObs.end());
755 std::sort(obsList.begin(), obsList.end());
756 obsList.erase(std::unique(obsList.begin(), obsList.end()), obsList.end());
759 dataSet.SetFilePattern(dsFileName);
760 dataSet.SetObservablesList(obsList);
762 fDataSetFileName = dsFileName;
765 dataSet.SetDataFrame(dataSet.
MakeCut(p->GetCut()));
767 if (fSplitX.empty()) SetSplitX();
768 if (fSplitY.empty()) SetSplitY();
771 if (fCalibRange.X() >= fCalibRange.Y()) {
773 std::string cut = fDefinitionCut;
774 if (cut.empty()) cut =
"1";
775 auto histo = dataSet.
GetDataFrame().Filter(cut).Histo1D({
"temp",
"", fNBins, 0, 0}, GetObservable());
776 std::unique_ptr<TH1F> hpunt = std::unique_ptr<TH1F>(
static_cast<TH1F*
>(histo->Clone()));
778 double xMax = hpunt->GetXaxis()->GetXmax();
781 double fraction = 1, nAtEndSpc = 0, step = 0.66;
782 while (nAtEndSpc * 1. / hpunt->Integral() < 0.01 && fraction > 0.001) {
784 nAtEndSpc = hpunt->Integral(hpunt->FindFixBin(hpunt->GetXaxis()->GetXmax() * fraction),
785 hpunt->FindFixBin(hpunt->GetXaxis()->GetXmax()));
787 xMax = hpunt->GetXaxis()->GetXmax() * fraction /
791 fCalibRange.SetY(xMax);
793 RESTDebug <<
"Calibration range (auto)set to (" << fCalibRange.X() <<
"," << fCalibRange.Y() <<
")"
798 std::string hModuleName =
"hSpc_" + std::to_string(fPlaneId) +
"_" + std::to_string(fModuleId);
799 delete fFullSpectrum;
800 fFullSpectrum =
new TH1F(hModuleName.c_str(),
"", fNBins, fCalibRange.X(), fCalibRange.Y());
803 std::string cut = fDefinitionCut;
804 if (cut.empty()) cut =
"1";
805 auto histoMod = dataSet.
GetDataFrame().Filter(cut).Histo1D(
806 {
"tempMod",
"", fNBins, fCalibRange.X(), fCalibRange.Y()}, GetObservable());
807 std::unique_ptr<TH1F> hpuntMod = std::unique_ptr<TH1F>(
static_cast<TH1F*
>(histoMod->Clone()));
808 fFullSpectrum->Add(hpuntMod.get());
811 std::vector<std::vector<TH1F*>> h(fNumberOfSegmentsX, std::vector<TH1F*>(fNumberOfSegmentsY,
nullptr));
812 for (
size_t i = 0; i < h.size(); i++) {
813 for (
size_t j = 0; j < h.at(0).size(); j++) {
814 std::string name = hModuleName +
"_" + std::to_string(i) +
"_" + std::to_string(j);
815 h[i][j] =
new TH1F(name.c_str(),
"", fNBins, fCalibRange.X(),
821 auto itX = fSplitX.begin();
822 for (
size_t i = 0; i < h.size(); i++) {
823 auto itY = fSplitY.begin();
824 for (
size_t j = 0; j < h.at(0).size(); j++) {
827 auto xUpper = *std::next(itX);
829 auto yUpper = *std::next(itY);
831 std::string segment_cut =
"";
832 if (!GetSpatialObservableX().empty())
833 segment_cut += GetSpatialObservableX() +
">=" + std::to_string(xLower) +
"&&" +
834 GetSpatialObservableX() +
"<" + std::to_string(xUpper);
835 if (!GetSpatialObservableY().empty())
836 segment_cut +=
"&&" + GetSpatialObservableY() +
">=" + std::to_string(yLower) +
"&&" +
837 GetSpatialObservableY() +
"<" + std::to_string(yUpper);
838 if (!fDefinitionCut.empty()) segment_cut +=
"&&" + fDefinitionCut;
839 if (segment_cut.empty()) segment_cut =
"1";
840 RESTExtreme <<
"Segment[" << i <<
"][" << j <<
"] cut: " << segment_cut << p->RESTendl;
843 .Histo1D({
"temp",
"", h[i][j]->GetNbinsX(), h[i][j]->GetXaxis()->GetXmin(),
844 h[i][j]->GetXaxis()->GetXmax()},
846 std::unique_ptr<TH1F> hpunt = std::unique_ptr<TH1F>(
static_cast<TH1F*
>(histo->Clone()));
847 h[i][j]->Add(hpunt.get());
855 std::vector<std::vector<double>> calParSlope(fNumberOfSegmentsX,
856 std::vector<double>(fNumberOfSegmentsY, 0));
857 std::vector<std::vector<double>> calParIntercept(fNumberOfSegmentsX,
858 std::vector<double>(fNumberOfSegmentsY, 0));
859 fSegLinearFit = std::vector(h.size(), std::vector<TGraph*>(h.at(0).size(),
nullptr));
860 for (
size_t i = 0; i < h.size(); i++)
861 for (
int j = 0; j < (int)h.at(0).size(); j++) {
862 fSegLinearFit[i][j] =
new TGraph();
863 auto [intercept, slope] = FitPeaks(h[i][j], fSegLinearFit[i][j]);
864 calParSlope[i][j] = slope;
865 calParIntercept[i][j] = intercept;
867 fSlope = calParSlope;
868 fIntercept = calParIntercept;
872 delete fFullLinearFit;
873 fFullLinearFit =
new TGraph();
874 auto [intercept, slope] = FitPeaks(fFullSpectrum, fFullLinearFit);
876 fFullIntercept = intercept;
879std::pair<double, double> TRestDataSetGainMap::Module::FitPeaks(TH1F* hSeg, TGraph* gr) {
881 RESTError <<
"No histogram for fitting" << p->RESTendl;
882 return std::make_pair(0, 0);
884 if (hSeg->Integral() == 0) {
885 RESTError <<
"Empty spectrum " << hSeg->GetName() << p->RESTendl;
886 return std::make_pair(0, 0);
888 std::shared_ptr<TGraph> graph = std::shared_ptr<TGraph>(
new TGraph());
889 RESTExtreme <<
"Fitting peaks for " << hSeg->GetName() << p->RESTendl;
892 std::unique_ptr<TSpectrum> s(
new TSpectrum(2 * fEnergyPeaks.size() + 1));
893 std::vector<double> peakPos;
894 s->Search(hSeg, 2,
"goff", 0.1);
895 for (
int k = 0; k < s->GetNPeaks(); k++) peakPos.push_back(s->GetPositionX()[k]);
896 std::sort(peakPos.begin(), peakPos.end(), std::greater<double>());
898 peakPos.size() == 0 ? 1 : peakPos.front() / fEnergyPeaks.front();
901 graph->SetName(
"grFit");
902 graph->SetTitle((
";" + GetObservable() +
";energy").c_str());
907 for (
const auto& energy : fEnergyPeaks) {
908 RESTExtreme <<
"\t fitting energy " <<
DoubleToString(energy,
"%g") << p->RESTendl;
910 double pos = energy * ratio;
911 double start = pos * 0.8;
912 double end = pos * 1.2;
913 if (fRangePeaks.at(c).X() < fRangePeaks.at(c).Y()) {
914 start = fRangePeaks.at(c).X();
915 end = fRangePeaks.at(c).Y();
919 if (fAutoRangePeaks) {
920 if (peakPos.size() > 0) {
923 while (!(start < pos && pos < end)) {
926 if (pos == peakPos.back()) {
930 pos = *std::next(std::find(peakPos.begin(), peakPos.end(),
933 peakPos.erase(std::find(peakPos.begin(), peakPos.end(),
939 const double relDist = peakPos.size() > 0 ? (pos - peakPos.front()) / pos : 999;
941 start = pos * (1 - relDist / 2);
942 end = pos * (1 + relDist / 2);
947 std::string name =
"g" + std::to_string(c);
948 TF1* g =
new TF1(name.c_str(),
"gaus", start, end);
949 RESTExtreme <<
"\t\tat " <<
DoubleToString(pos,
"%.3g") <<
". Range("
953 if (hSeg->GetFunction(name.c_str()))
954 hSeg->GetListOfFunctions()->Remove(hSeg->GetFunction(name.c_str()));
956 hSeg->Fit(g,
"R+Q0");
957 mu = g->GetParameter(1);
958 RESTExtreme <<
"\t\tgaus mean " <<
DoubleToString(mu,
"%g") << p->RESTendl;
959 }
while (fAutoRangePeaks && peakPos.size() > 0 &&
960 !(start < mu && mu < end));
961 graph->SetPoint(c++, mu, energy);
965 if (fZeroPoint) graph->SetPoint(c++, 0, 0);
966 while (graph->GetN() < 2) {
967 graph->SetPoint(c++, 0, 0);
969 RESTDebug <<
"Not enough points for linear fit. Adding and setting zero point to true" << p->RESTendl;
973 std::unique_ptr<TF1> linearFit;
974 linearFit = std::unique_ptr<TF1>(
new TF1(
"linearFit",
"pol1"));
975 graph->Fit(
"linearFit",
"SQ");
977 if (gr) *gr = *(TGraph*)graph->Clone();
978 return std::make_pair(linearFit->GetParameter(0), linearFit->GetParameter(1));
989 const TVector2& range) {
990 auto [index_x, index_y] = GetIndexMatrix(position.X(), position.Y());
992 for (
size_t i = 0; i < fEnergyPeaks.size(); i++)
993 if (fEnergyPeaks.at(i) == energyPeak) {
997 if (peakNumber == -1) {
998 RESTError <<
"Energy " << energyPeak <<
" not found in the list of energy peaks" << p->RESTendl;
1001 Refit((
size_t)index_x, (
size_t)index_y, (
size_t)peakNumber, range);
1014 const TVector2& range) {
1015 if (fSegSpectra.empty()) {
1016 RESTError <<
"No gain map found. Use GenerateGainMap() first." << p->RESTendl;
1019 if (x >= fSegSpectra.size() || y >= fSegSpectra.at(0).size()) {
1020 RESTError <<
"Segment with index (" << x <<
", " << y <<
") not found" << p->RESTendl;
1023 if (peakNumber >= fEnergyPeaks.size()) {
1024 RESTError <<
"Peak with index " << peakNumber <<
" not found" << p->RESTendl;
1029 std::string name =
"g" + std::to_string(peakNumber);
1030 TF1* g =
new TF1(name.c_str(),
"gaus", range.X(), range.Y());
1031 TH1F* h = fSegSpectra.at(x).at(y);
1032 while (h->GetFunction(name.c_str()))
1033 h->GetListOfFunctions()->Remove(h->GetFunction(name.c_str()));
1037 UpdateCalibrationFits(x, y);
1048 int peakNumber = -1;
1049 for (
size_t i = 0; i < fEnergyPeaks.size(); i++)
1050 if (fEnergyPeaks.at(i) == energyPeak) {
1054 if (peakNumber == -1) {
1055 RESTError <<
"Energy " << energyPeak <<
" not found in the list of energy peaks" << p->RESTendl;
1058 RefitFullSpc((
size_t)peakNumber, range);
1069 if (!fFullSpectrum) {
1070 RESTError <<
"No gain map found. Use GenerateGainMap() first." << p->RESTendl;
1073 if (peakNumber >= fEnergyPeaks.size()) {
1074 RESTError <<
"Peak with index " << peakNumber <<
" not found" << p->RESTendl;
1079 std::string name =
"g" + std::to_string(peakNumber);
1080 TF1* g =
new TF1(name.c_str(),
"gaus", range.X(), range.Y());
1081 while (fFullSpectrum->GetFunction(name.c_str()))
1082 fFullSpectrum->GetListOfFunctions()->Remove(fFullSpectrum->GetFunction(name.c_str()));
1083 fFullSpectrum->Fit(g,
"R+Q0");
1086 UpdateCalibrationFitsFullSpc();
1097void TRestDataSetGainMap::Module::UpdateCalibrationFits(
const size_t x,
const size_t y) {
1098 if (fSegSpectra.empty()) {
1099 RESTError <<
"No gain map found. Use GenerateGainMap() first." << p->RESTendl;
1102 if (x >= fSegSpectra.size() || y >= fSegSpectra.at(0).size()) {
1103 RESTError <<
"Segment with index (" << x <<
", " << y <<
") not found" << p->RESTendl;
1107 TH1F* h = fSegSpectra.at(x).at(y);
1108 TGraph* gr = fSegLinearFit.at(x).at(y);
1110 auto [intercept, slope] = UpdateCalibrationFits(h, gr);
1111 fSlope[x][y] = slope;
1112 fIntercept[x][y] = intercept;
1115std::pair<double, double> TRestDataSetGainMap::Module::UpdateCalibrationFits(TH1F* h, TGraph* gr) {
1117 RESTError <<
"No histogram for updating fits" << p->RESTendl;
1118 return std::make_pair(0, 0);
1121 RESTError <<
"No graph for updating fits" << p->RESTendl;
1122 return std::make_pair(0, 0);
1124 if (h->Integral() == 0) {
1125 RESTError <<
"Empty spectrum " << h->GetName() << p->RESTendl;
1126 return std::make_pair(0, 0);
1130 for (
size_t i = 0; i < fEnergyPeaks.size(); i++) gr->RemovePoint(i);
1133 for (
size_t i = 0; i < fEnergyPeaks.size(); i++) {
1134 std::string fitName = (std::string)
"g" + std::to_string(i);
1135 TF1* g = h->GetFunction(fitName.c_str());
1137 RESTWarning <<
"No fit ( " << fitName <<
" ) found for energy peak " << fEnergyPeaks[i]
1138 <<
" in histogram " << h->GetName() << p->RESTendl;
1141 gr->SetPoint(c++, g->GetParameter(1), fEnergyPeaks[i]);
1145 while (gr->GetN() < 2) {
1146 gr->SetPoint(c++, 0, 0);
1151 if (gr->GetFunction(
"linearFit"))
1152 lf = gr->GetFunction(
"linearFit");
1154 lf =
new TF1(
"linearFit",
"pol1");
1157 return std::make_pair(lf->GetParameter(0), lf->GetParameter(1));
1166 auto [intercept, slope] = UpdateCalibrationFits(fFullSpectrum, fFullLinearFit);
1168 fFullIntercept = intercept;
1186 if (module ==
nullptr) {
1187 RESTError <<
"TRestDataSetGainMap::Module::LoadConfigFromTiXmlElement: module is nullptr"
1192 std::string el = !
module->Attribute("planeId") ? "Not defined" : module->Attribute("planeId");
1193 if (!(el.empty() || el ==
"Not defined")) this->SetPlaneId(
StringToInteger(el));
1194 el = !
module->Attribute("moduleId") ? "Not defined" : module->Attribute("moduleId");
1195 if (!(el.empty() || el ==
"Not defined")) this->SetModuleId(
StringToInteger(el));
1197 el = !
module->Attribute("moduleDefinitionCut") ? "Not defined" : module->Attribute("moduleDefinitionCut");
1198 if (!(el.empty() || el ==
"Not defined")) this->SetModuleDefinitionCut(el);
1200 el = !
module->Attribute("numberOfSegmentsX") ? "Not defined" : module->Attribute("numberOfSegmentsX");
1201 if (!(el.empty() || el ==
"Not defined")) this->SetNumberOfSegmentsX(
StringToInteger(el));
1202 el = !
module->Attribute("numberOfSegmentsY") ? "Not defined" : module->Attribute("numberOfSegmentsY");
1203 if (!(el.empty() || el ==
"Not defined")) this->SetNumberOfSegmentsY(
StringToInteger(el));
1204 el = !
module->Attribute("readoutRange") ? "Not defined" : module->Attribute("readoutRange");
1205 if (!(el.empty() || el ==
"Not defined")) this->SetReadoutRange(
StringTo2DVector(el));
1207 el = !
module->Attribute("calibRange") ? "Not defined" : module->Attribute("calibRange");
1208 if (!(el.empty() || el ==
"Not defined")) this->SetCalibrationRange(
StringTo2DVector(el));
1209 el = !
module->Attribute("nBins") ? "Not defined" : module->Attribute("nBins");
1210 if (!(el.empty() || el ==
"Not defined")) this->SetNBins(
StringToInteger(el));
1212 el = !
module->Attribute("dataSetFileName") ? "Not defined" : module->Attribute("dataSetFileName");
1213 if (!(el.empty() || el ==
"Not defined")) this->SetDataSetFileName(el);
1215 el = !
module->Attribute("zeroPoint") ? "Not defined" : module->Attribute("zeroPoint");
1216 if (!(el.empty() || el ==
"Not defined")) this->SetZeroPoint(
ToLower(el) ==
"true");
1217 el = !
module->Attribute("autoRangePeaks") ? "Not defined" : module->Attribute("autoRangePeaks");
1218 if (!(el.empty() || el ==
"Not defined")) this->SetAutoRangePeaks(
ToLower(el) ==
"true");
1221 TiXmlElement* peakDefinition = (TiXmlElement*)module->FirstChildElement(
"peak");
1222 while (peakDefinition !=
nullptr) {
1224 TVector2 range = TVector2(0, 0);
1227 !peakDefinition->Attribute(
"energy") ?
"Not defined" : peakDefinition->Attribute(
"energy");
1228 if (ell.empty() || ell ==
"Not defined") {
1229 RESTError <<
"< peak variable key does not contain energy!" << p->RESTendl;
1234 ell = !peakDefinition->Attribute(
"range") ?
"Not defined" : peakDefinition->Attribute(
"range");
1237 this->AddPeak(energy, range);
1238 peakDefinition = (TiXmlElement*)peakDefinition->NextSiblingElement();
1244 std::pair<size_t, size_t> index = GetIndexMatrix(position.X(), position.Y());
1245 DrawSpectrum(index.first, index.second, drawFits, color, c);
1250 if (fSegSpectra.size() == 0) {
1251 RESTError <<
"Spectra matrix is empty." << p->RESTendl;
1254 if (index_x < 0 || index_y < 0 || index_x >= (
int)fSegSpectra.size() ||
1255 index_y >= (
int)fSegSpectra.at(index_x).size()) {
1256 RESTError <<
"Index out of range." << p->RESTendl;
1259 if (!fSegSpectra[index_x][index_y]) {
1260 RESTError <<
"No Spectrum for segment (" << index_x <<
", " << index_y <<
")." << p->RESTendl;
1265 std::string t =
"spectrum_" + std::to_string(fPlaneId) +
"_" + std::to_string(fModuleId) +
"_" +
1266 std::to_string(index_x) +
"_" + std::to_string(index_y);
1267 c =
new TCanvas(t.c_str(), t.c_str());
1270 auto xLower = *std::next(fSplitX.begin(), index_x);
1271 auto xUpper = *std::next(fSplitX.begin(), index_x + 1);
1272 auto yLower = *std::next(fSplitY.begin(), index_y);
1273 auto yUpper = *std::next(fSplitY.begin(), index_y + 1);
1276 GetObservable() +
";counts";
1277 fSegSpectra[index_x][index_y]->SetTitle(tH.c_str());
1279 if (color > 0) fSegSpectra[index_x][index_y]->SetLineColor(color);
1280 size_t colorT = fSegSpectra[index_x][index_y]->GetLineColor();
1281 fSegSpectra[index_x][index_y]->Draw(
"same");
1284 for (
size_t c = 0; c < fEnergyPeaks.size(); c++) {
1285 auto fit = fSegSpectra[index_x][index_y]->GetFunction((
"g" + std::to_string(c)).c_str());
1287 RESTWarning <<
"Fit for energy peak " << fEnergyPeaks[c] <<
" not found." << p->RESTendl;
1289 fit->SetLineColor(c + 2 != colorT ? c + 2 : ++colorT);
1317 if (fSegSpectra.size() == 0) {
1318 RESTError <<
"Spectra matrix is empty." << p->RESTendl;
1322 std::string t =
"spectrum_" + std::to_string(fPlaneId) +
"_" + std::to_string(fModuleId);
1323 c =
new TCanvas(t.c_str(), t.c_str());
1327 for (
const auto&
object : *c->GetListOfPrimitives())
1328 if (object->InheritsFrom(TVirtualPad::Class())) ++nPads;
1329 if (nPads != 0 && nPads != fSegSpectra.size() * fSegSpectra.at(0).size()) {
1330 RESTError <<
"Canvas " << c->GetName() <<
" has " << nPads <<
" pads, but "
1331 << fSegSpectra.size() * fSegSpectra.at(0).size() <<
" are needed." << p->RESTendl;
1333 }
else if (nPads == 0)
1334 c->Divide(fSegSpectra.size(), fSegSpectra.at(0).size());
1336 for (
size_t i = 0; i < fSegSpectra.size(); i++) {
1337 for (
size_t j = 0; j < fSegSpectra[i].size(); j++) {
1338 int pad = fSegSpectra.size() * (fSegSpectra[i].size() - 1) + 1 + i - fSegSpectra.size() * j;
1340 DrawSpectrum(i, j, drawFits, color, c);
1344void TRestDataSetGainMap::Module::DrawFullSpectrum(
const bool drawFits,
const int color, TCanvas* c) {
1345 if (!fFullSpectrum) {
1346 RESTError <<
"Spectrum is empty." << p->RESTendl;
1351 std::string t =
"fullSpc_" + std::to_string(fPlaneId) +
"_" + std::to_string(fModuleId);
1352 c =
new TCanvas(t.c_str(), t.c_str());
1356 fFullSpectrum->SetTitle((
"Full spectrum;" + GetObservable() +
";counts").c_str());
1358 if (color > 0) fFullSpectrum->SetLineColor(color);
1359 size_t colorT = fFullSpectrum->GetLineColor();
1360 fFullSpectrum->Draw(
"same");
1363 for (
size_t c = 0; c < fEnergyPeaks.size(); c++) {
1364 auto fit = fFullSpectrum->GetFunction((
"g" + std::to_string(c)).c_str());
1365 if (!fit) RESTWarning <<
"Fit for energy peak" << fEnergyPeaks[c] <<
" not found." << p->RESTendl;
1367 fit->SetLineColor(c + 2 != colorT ? c + 2 : ++colorT);
1375void TRestDataSetGainMap::Module::DrawLinearFit(
const TVector2& position, TCanvas* c) {
1376 std::pair<size_t, size_t> index = GetIndexMatrix(position.X(), position.Y());
1377 DrawLinearFit(index.first, index.second, c);
1380void TRestDataSetGainMap::Module::DrawLinearFit(
const int index_x,
const int index_y, TCanvas* c) {
1381 if (fSegLinearFit.size() == 0) {
1382 RESTError <<
"Spectra matrix is empty." << p->RESTendl;
1385 if (index_x < 0 || index_y < 0 || index_x >= (
int)fSegLinearFit.size() ||
1386 index_y >= (
int)fSegLinearFit.at(index_x).size()) {
1387 RESTError <<
"Index out of range." << p->RESTendl;
1390 if (!fSegLinearFit[index_x][index_y]) {
1391 RESTError <<
"No linear fit for segment (" << index_x <<
", " << index_y <<
")." << p->RESTendl;
1396 std::string t =
"linearFit_" + std::to_string(fPlaneId) +
"_" + std::to_string(fModuleId) +
"_" +
1397 std::to_string(index_x) +
"_" + std::to_string(index_y);
1398 c =
new TCanvas(t.c_str(), t.c_str());
1400 auto xLower = *std::next(fSplitX.begin(), index_x);
1401 auto xUpper = *std::next(fSplitX.begin(), index_x + 1);
1402 auto yLower = *std::next(fSplitY.begin(), index_y);
1403 auto yUpper = *std::next(fSplitY.begin(), index_y + 1);
1406 GetObservable() +
";energy";
1407 fSegLinearFit[index_x][index_y]->SetTitle(tH.c_str());
1408 fSegLinearFit[index_x][index_y]->Draw(
"AL*");
1411void TRestDataSetGainMap::Module::DrawLinearFit(TCanvas* c) {
1412 if (fSegLinearFit.size() == 0) {
1413 RESTError <<
"Spectra matrix is empty." << p->RESTendl;
1417 std::string t =
"linearFits_" + std::to_string(fPlaneId) +
"_" + std::to_string(fModuleId);
1418 c =
new TCanvas(t.c_str(), t.c_str());
1422 for (
const auto&
object : *c->GetListOfPrimitives())
1423 if (object->InheritsFrom(TVirtualPad::Class())) ++nPads;
1424 if (nPads != 0 && nPads != fSegLinearFit.size() * fSegLinearFit.at(0).size()) {
1425 RESTError <<
"Canvas " << c->GetName() <<
" has " << nPads <<
" pads, but "
1426 << fSegLinearFit.size() * fSegLinearFit.at(0).size() <<
" are needed." << p->RESTendl;
1428 }
else if (nPads == 0)
1429 c->Divide(fSegLinearFit.size(), fSegLinearFit.at(0).size());
1431 for (
size_t i = 0; i < fSegLinearFit.size(); i++) {
1432 for (
size_t j = 0; j < fSegLinearFit[i].size(); j++) {
1433 int pad = fSegLinearFit.size() * (fSegLinearFit[i].size() - 1) + 1 + i - fSegLinearFit.size() * j;
1435 DrawLinearFit(i, j, c);
1449 if (peakNumber < 0 || peakNumber >= (
int)fEnergyPeaks.size()) {
1450 RESTError <<
"Peak number out of range (peakNumber should be between 0 and "
1451 << fEnergyPeaks.size() - 1 <<
" )" << p->RESTendl;
1454 if (fSegLinearFit.size() == 0) {
1455 RESTError <<
"Linear fit matrix is empty." << p->RESTendl;
1458 if (!fFullLinearFit) {
1459 RESTError <<
"Full linear fit is empty." << p->RESTendl;
1463 double peakEnergy = fEnergyPeaks[peakNumber];
1464 std::string title =
"Gain map for energy " +
DoubleToString(peakEnergy,
"%g") +
";" +
1465 GetSpatialObservableX() +
";" + GetSpatialObservableY();
1466 std::string t =
"gainMap" + std::to_string(peakNumber) +
"_" + std::to_string(fPlaneId) +
"_" +
1467 std::to_string(fModuleId);
1468 TCanvas* gainMap =
new TCanvas(t.c_str(), t.c_str());
1470 TH2F* hGainMap =
new TH2F((
"h" + t).c_str(), title.c_str(), fNumberOfSegmentsX, fReadoutRange.X(),
1471 fReadoutRange.Y(), fNumberOfSegmentsY, fReadoutRange.X(), fReadoutRange.Y());
1473 double peakPosRef = fFullLinearFit->GetPointX(peakNumber);
1474 if (!fullModuleAsRef) {
1475 int index_x = fNumberOfSegmentsX > 0 ? (fNumberOfSegmentsX - 1) / 2 : 0;
1476 int index_y = fNumberOfSegmentsY > 0 ? (fNumberOfSegmentsY - 1) / 2 : 0;
1477 peakPosRef = fSegLinearFit[index_x][index_y]->GetPointX(peakNumber);
1480 auto itX = fSplitX.begin();
1481 for (
size_t i = 0; i < fSegLinearFit.size(); i++) {
1482 auto itY = fSplitY.begin();
1483 for (
size_t j = 0; j < fSegLinearFit.at(0).size(); j++) {
1485 auto xUpper = *std::next(itX);
1487 auto yUpper = *std::next(itY);
1488 float xMean = (xUpper + xLower) / 2.;
1489 float yMean = (yUpper + yLower) / 2.;
1490 auto [index_x, index_y] = GetIndexMatrix(xMean, yMean);
1491 if (!fSegLinearFit[index_x][index_y])
continue;
1492 hGainMap->Fill(xMean, yMean, fSegLinearFit[index_x][index_y]->GetPointX(peakNumber) / peakPosRef);
1497 hGainMap->SetStats(0);
1498 hGainMap->Draw(
"colz");
1499 hGainMap->SetBarOffset(0.2);
1500 hGainMap->Draw(
"TEXT SAME");
1508 RESTMetadata <<
"-----------------------------------------------" << p->RESTendl;
1509 RESTMetadata <<
" Plane ID: " << fPlaneId << p->RESTendl;
1510 RESTMetadata <<
" Module ID: " << fModuleId << p->RESTendl;
1511 RESTMetadata <<
" Definition cut: " << fDefinitionCut << p->RESTendl;
1512 RESTMetadata << p->RESTendl;
1514 RESTMetadata <<
" Calibration dataset: " << fDataSetFileName << p->RESTendl;
1515 RESTMetadata << p->RESTendl;
1517 RESTMetadata <<
" Energy peaks: ";
1518 for (
const auto& peak : fEnergyPeaks) RESTMetadata << peak <<
", ";
1519 RESTMetadata << p->RESTendl;
1520 bool anyRange =
false;
1521 for (
const auto& r : fRangePeaks)
1522 if (r.X() < r.Y()) {
1523 RESTMetadata <<
" Range peaks: ";
1528 for (
const auto& r : fRangePeaks) RESTMetadata <<
"(" << r.X() <<
", " << r.Y() <<
") ";
1529 if (anyRange) RESTMetadata << p->RESTendl;
1530 RESTMetadata <<
" Auto range peaks: " << (fAutoRangePeaks ?
"true" :
"false") << p->RESTendl;
1531 RESTMetadata <<
" Zero point: " << (fZeroPoint ?
"true" :
"false") << p->RESTendl;
1532 RESTMetadata <<
" Calibration range: (" << fCalibRange.X() <<
", " << fCalibRange.Y() <<
" )"
1534 RESTMetadata <<
" Number of bins: " << fNBins << p->RESTendl;
1535 RESTMetadata << p->RESTendl;
1537 RESTMetadata <<
" Number of segments X: " << fNumberOfSegmentsX << p->RESTendl;
1538 RESTMetadata <<
" Number of segments Y: " << fNumberOfSegmentsY << p->RESTendl;
1540 RESTMetadata <<
" Readout range (" << fReadoutRange.X() <<
", " << fReadoutRange.Y() <<
" )"
1542 RESTMetadata <<
"SplitX: ";
1543 for (
auto& i : fSplitX) {
1544 RESTMetadata <<
" " << i;
1546 RESTMetadata << p->RESTendl;
1547 RESTMetadata <<
"SplitY: ";
1548 for (
auto& i : fSplitY) {
1549 RESTMetadata <<
" " << i;
1551 RESTMetadata << p->RESTendl;
1552 RESTMetadata << p->RESTendl;
1554 RESTMetadata <<
" Slope: " << p->RESTendl;
1556 for (
auto& x : fSlope)
1557 if (maxSize < x.size()) maxSize = x.size();
1558 for (
size_t j = 0; j < maxSize; j++) {
1559 for (
size_t k = 0; k < fSlope.size(); k++) {
1560 if (j < fSlope[k].size())
1561 RESTMetadata <<
DoubleToString(fSlope[k][fSlope[k].size() - 1 - j],
"%.3e") <<
" ";
1563 RESTMetadata <<
" ";
1565 RESTMetadata << p->RESTendl;
1567 RESTMetadata <<
" Intercept: " << p->RESTendl;
1569 for (
auto& x : fIntercept)
1570 if (maxSize < x.size()) maxSize = x.size();
1571 for (
size_t j = 0; j < maxSize; j++) {
1572 for (
size_t k = 0; k < fIntercept.size(); k++) {
1573 if (j < fIntercept[k].size())
1574 RESTMetadata <<
DoubleToString(fIntercept[k][fIntercept[k].size() - 1 - j],
"%+.3e") <<
" ";
1576 RESTMetadata <<
" ";
1578 RESTMetadata << p->RESTendl;
1580 RESTMetadata << p->RESTendl;
1581 RESTMetadata <<
" Full slope: " <<
DoubleToString(fFullSlope,
"%.3e") << p->RESTendl;
1582 RESTMetadata <<
" Full intercept: " <<
DoubleToString(fFullIntercept,
"%+.3e") << p->RESTendl;
1584 RESTMetadata <<
"-----------------------------------------------" << p->RESTendl;
A class to help on cuts definitions. To be used with TRestAnalysisTree.
void SetSplitY()
Function to set the class members for segmentation of the detector plane along the Y axis.
void SetSplits()
Function to set the class members for segmentation of the detector plane along the X and Y axis.
void SetSplitX()
Function to set the class members for segmentation of the detector plane along the X axis.
void DrawSpectrum(const bool drawFits=true, const int color=-1, TCanvas *c=nullptr)
Function to draw the spectrum for each segment of the module on the same canvas. The canvas is divide...
void LoadConfigFromTiXmlElement(const TiXmlElement *module)
Function to read the parameters from the RML element (TiXmlElement) and set those class members.
void GenerateGainMap()
Function that calculates the calibration parameters for each segment defined at fSplitX and fSplitY a...
double GetSlope(const double x, const double y) const
Function to get the calibration parameter slope for a given x and y position on the detector plane.
std::set< double > fSplitX
Split points in the x direction.
void Print() const
Prints on screen the information about the members of Module.
void Refit(const TVector2 &position, const double energy, const TVector2 &range)
Function to fit again manually a peak for a given segment of the module.
std::set< double > fSplitY
Split points in the y direction.
void RefitFullSpc(const double energy, const TVector2 &range)
Function to fit again manually a peak for the whole module spectrum. The calibration curve is updated...
const TRestDataSetGainMap * p
Pointer to the parent class.
void UpdateCalibrationFitsFullSpc()
Function to update the calibration curve for the whole module. The calibration curve is cleared and t...
void DrawGainMap(const int peakNumber=0, const bool fullModuleAsRef=true)
Function to draw the relative gain map for a given energy peak of the module.
double GetIntercept(const double x, const double y) const
Function to get the calibration parameter intercept for a given x and y position on the detector plan...
std::pair< int, int > GetIndexMatrix(const double x, const double y) const
Function to get the index of the matrix of calibration parameters for a given x and y position on the...
Metadata class to calculate,store and apply the gain corrected calibration of a group of detectors.
double GetInterceptParameter(const int planeID, const int moduleID, const double x, const double y)
Function to get the intercept parameter of the module with planeID and moduleID at physical position ...
double GetSlopeParameter(const int planeID, const int moduleID, const double x, const double y)
Function to get the slope parameter of the module with planeID and moduleID at physical position (x,...
std::string fObservable
Observable that will be used to calculate the gain map.
double GetInterceptParameterFullSpc(const int planeID, const int moduleID)
Function to get the intercept parameter of the whole module with planeID and moduleID.
void CalibrateDataSet(const std::string &dataSetFileName, std::string outputFileName="", std::vector< std::string > excludeColumns={})
Function to calibrate a dataset with this gain map.
~TRestDataSetGainMap()
Default destructor.
void GenerateGainMap()
Function to calculate the calibration parameters of all modules.
TRestDataSetGainMap()
Default constructor.
void Initialize() override
Making default settings.
std::vector< Module > fModulesCal
List of modules.
void SetModule(const Module &moduleCal)
Function to set a module calibration. If the module calibration already exists (same planeId and modu...
void InitFromConfigFile() override
Initialization of TRestDataSetGainMap members through a RML file.
double GetSlopeParameterFullSpc(const int planeID, const int moduleID)
Function to get the slope parameter of the whole module with planeID and moduleID.
void Export(const std::string &fileName="")
Function to export the calibration to the file fileName.
void PrintMetadata() override
Prints on screen the information about the metadata members.
std::map< int, std::set< int > > GetModuleIDs() const
Function to get the map of the module IDs for each plane ID.
std::string fSpatialObservableY
Observable that will be used to segmentize the gain map in the y direction.
std::string fCalibFileName
Name of the file that contains the calibration data.
TRestCut * fCut
Cut to be applied to the calibration data.
std::string fOutputFileName
Name of the file where the gain map was (or will be) exported.
std::set< int > GetPlaneIDs() const
Function to get a list (set) of the plane IDs.
std::string fSpatialObservableX
Observable that will be used to segmentize the gain map in the x direction.
Module * GetModule(const size_t index=0)
Function to retrieve the module calibration by index. Default is 0.
void Import(const std::string &fileName)
Function to import the calibration parameters from the root file fileName.
It allows to group a number of runs that satisfy given metadata conditions.
void Import(const std::string &fileName)
This function imports metadata from a root file it import metadata info from the previous dataSet whi...
ROOT::RDF::RNode GetDataFrame() const
Gives access to the RDataFrame.
ROOT::RDF::RNode MakeCut(const TRestCut *cut)
This function applies a TRestCut to the dataframe and returns a dataframe with the applied cuts....
void GenerateDataSet()
This function generates the data frame with the filelist and column names (or observables) that have ...
void Export(const std::string &filename, std::vector< std::string > excludeColumns={})
It will generate an output file with the dataset compilation. Only the selected branches and the file...
@ REST_Info
+show most of the information for each steps
@ REST_Debug
+show the defined debug messages
TClass * GetClassQuick()
Get the type of a "class" object, returning the wrapped type identifier "TClass".
Double_t StringToDouble(std::string in)
Gets a double from a string.
Int_t StringToInteger(std::string in)
Gets an integer from a string.
std::string DoubleToString(Double_t d, std::string format="%8.6e")
Gets a string from a double.
TVector2 StringTo2DVector(std::string in)
Gets a 2D-vector from a string.
std::string ToLower(std::string in)
Convert string to its lower case. Alternative of TString::ToLower.