REST-for-Physics  v2.3
Rare Event Searches ToolKit for Physics
Loading...
Searching...
No Matches
TRestDataSetGainMap.cxx
1/*************************************************************************
2 * This file is part of the REST software framework. *
3 * *
4 * Copyright (C) 2016 GIFNA/TREX (University of Zaragoza) *
5 * For more information see https://gifna.unizar.es/trex *
6 * *
7 * REST is free software: you can redistribute it and/or modify *
8 * it under the terms of the GNU General Public License as published by *
9 * the Free Software Foundation, either version 3 of the License, or *
10 * (at your option) any later version. *
11 * *
12 * REST is distributed in the hope that it will be useful, *
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15 * GNU General Public License for more details. *
16 * *
17 * You should have a copy of the GNU General Public License along with *
18 * REST in $REST_PATH/LICENSE. *
19 * If not, see https://www.gnu.org/licenses/. *
20 * For the list of contributors see $REST_PATH/CREDITS. *
21 *************************************************************************/
22
139
140#include "TRestDataSetGainMap.h"
141
142ClassImp(TRestDataSetGainMap);
147
162TRestDataSetGainMap::TRestDataSetGainMap(const char* configFilename, std::string name)
163 : TRestMetadata(configFilename) {
166}
167
172
179 this->Initialize();
181
182 // Load Module from RML
183 TiXmlElement* moduleDefinition = GetElement("module");
184 while (moduleDefinition != nullptr) {
185 fModulesCal.push_back(Module(*this));
186 fModulesCal.back().LoadConfigFromTiXmlElement(moduleDefinition);
187
188 moduleDefinition = GetNextElement(moduleDefinition);
189 }
190
191 fCut = (TRestCut*)InstantiateChildMetadata("TRestCut");
192
194}
195
200 for (auto& mod : fModulesCal) {
201 RESTInfo << "Generating gain map of plane " << mod.GetPlaneId() << " module " << mod.GetModuleId()
202 << RESTendl;
203 mod.GenerateGainMap();
205 mod.DrawSpectrum();
206 mod.DrawFullSpectrum();
207 mod.DrawGainMap();
208 }
209 }
210}
211
226void TRestDataSetGainMap::CalibrateDataSet(const std::string& dataSetFileName, std::string outputFileName,
227 std::vector<std::string> excludeColumns) {
228 if (fModulesCal.empty()) {
229 RESTError << "TRestDataSetGainMap::CalibrateDataSet: No modules defined." << RESTendl;
230 return;
231 }
232
233 TRestDataSet dataSet;
234 dataSet.EnableMultiThreading(true);
235
236 if (TRestTools::isDataSet(dataSetFileName)) {
237 dataSet.Import(dataSetFileName);
238 } else {
239 RESTWarning << dataSetFileName << " is not a dataset. Generating a temporal one..." << RESTendl;
240 // generate the dataset with the needed observables
241 dataSet.SetFilePattern(dataSetFileName);
242 dataSet.SetObservablesList({"*"}); // get all observables
243 dataSet.GenerateDataSet();
244 }
245
246 auto dataFrame = dataSet.GetDataFrame();
247
248 // Define a new column with the identifier (pmID) of the module for each row (event)
249 std::string pmIDname = (std::string)GetName() + "_pmID";
250 std::string modCut = fModulesCal[0].GetModuleDefinitionCut();
251 if (modCut.empty()) modCut = "1"; // if no cut is defined, use "1" (all events)
252 int pmID = fModulesCal[0].GetPlaneId() * 10 + fModulesCal[0].GetModuleId();
253
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");
257 else
258 dataFrame = dataFrame.Redefine(pmIDname, modCut + " ? " + std::to_string(pmID) + " : -1");
259
260 for (size_t n = 1; n < fModulesCal.size(); n++) {
261 modCut = fModulesCal[n].GetModuleDefinitionCut();
262 if (modCut.empty()) modCut = "1"; // if no cut is defined, use "1" (all events)
263 pmID = fModulesCal[n].GetPlaneId() * 10 + fModulesCal[n].GetModuleId();
264 dataFrame = dataFrame.Redefine(pmIDname, (modCut + " ? " + std::to_string(pmID) + " : " + pmIDname));
265 }
266
267 // Define a new column with the calibrated observable
268 auto calibrate = [this](double val, double x, double y, int pmID) {
269 for (auto& m : fModulesCal) {
270 if (pmID == m.GetPlaneId() * 10 + m.GetModuleId())
271 return m.GetSlope(x, y) * val + m.GetIntercept(x, y);
272 }
273 // RESTError << "TRestDataSetGainMap::CalibrateDataSet: Module not found for pmID " << pmID <<
274 // RESTendl;
275 return std::numeric_limits<double>::quiet_NaN();
276 };
277 std::string calibObsName = (std::string)GetName() + "_";
278 calibObsName += GetObservable().erase(0, GetObservable().find("_") + 1); // remove the "rawAna_" part
279 dataFrame = dataFrame.Define(calibObsName, calibrate,
281
282 // Define a new column with the calibrated observable for the whole module calibration
283 auto calibrateFullSpc = [this](double val, int pmID) {
284 for (auto& m : fModulesCal) {
285 if (pmID == m.GetPlaneId() * 10 + m.GetModuleId())
286 return m.GetSlopeFullSpc() * val + m.GetInterceptFullSpc();
287 }
288 // RESTError << "TRestDataSetGainMap::CalibrateDataSet: Module not found for pmID " << pmID <<
289 // RESTendl;
290 return std::numeric_limits<double>::quiet_NaN();
291 };
292 std::string calibObsNameFullSpc = (std::string)GetName() + "_";
293 calibObsNameFullSpc +=
294 GetObservable().erase(0, GetObservable().find("_") + 1); // remove the "rawAna_" part
295 calibObsNameFullSpc += "_NoSegmentation";
296 dataFrame = dataFrame.Define(calibObsNameFullSpc, calibrateFullSpc, {fObservable, pmIDname});
297
298 dataSet.SetDataFrame(dataFrame);
299
300 // Format the output file name and export the dataSet
301 if (outputFileName.empty()) outputFileName = dataSetFileName;
302 if (outputFileName == dataSetFileName) { // TRestDataSet cannot be overwritten
303 std::string gmName = GetName();
304 outputFileName = outputFileName.substr(0, outputFileName.find_last_of(".")); // remove extension
305 outputFileName += "_" + gmName + "." + TRestTools::GetFileNameExtension(dataSetFileName);
306 }
307
308 // Export dataset. Exclude columns if requested.
309 auto columns = dataSet.GetDataFrame().GetColumnNames();
310 std::set<std::string> excludeCol = TRestTools::GetMatchingStrings(columns, excludeColumns);
311 // Never exclude the calibObsName, calibObsNameFullSpc and pmIDname
312 excludeCol.erase(calibObsName);
313 excludeCol.erase(calibObsNameFullSpc);
314 excludeCol.erase(pmIDname);
315
316 RESTDebug << "Excluding columns: ";
317 for (auto& c : excludeCol) RESTDebug << c << ", ";
318 RESTDebug << RESTendl;
319
320 dataSet.Export(outputFileName, std::vector<std::string>(excludeCol.begin(), excludeCol.end()));
321
322 // Add this TRestDataSetGainMap metadata to the output file
323 TFile* f = TFile::Open(outputFileName.c_str(), "UPDATE");
324 this->Write();
325 f->Close();
326 delete f;
327}
328
334 if (index < fModulesCal.size()) return &fModulesCal[index];
335
336 RESTError << "No ModuleCalibration with index " << index;
337 if (fModulesCal.empty())
338 RESTError << ". There are no modules defined." << RESTendl;
339 else
340 RESTError << ". Max index is " << fModulesCal.size() - 1 << RESTendl;
341 return nullptr;
342}
343
348TRestDataSetGainMap::Module* TRestDataSetGainMap::GetModule(const int planeID, const int moduleID) {
349 for (auto& i : fModulesCal) {
350 if (i.GetPlaneId() == planeID && i.GetModuleId() == moduleID) return &i;
351 }
352 RESTError << "No ModuleCalibration with planeID " << planeID << " and moduleID " << moduleID << RESTendl;
353 return nullptr;
354}
355
361double TRestDataSetGainMap::GetSlopeParameter(const int planeID, const int moduleID, const double x,
362 const double y) {
363 Module* moduleCal = GetModule(planeID, moduleID);
364 if (moduleCal == nullptr) return 0; // return numeric_limits<double>::quiet_NaN()
365 return moduleCal->GetSlope(x, y);
366}
367
372double TRestDataSetGainMap::GetSlopeParameterFullSpc(const int planeID, const int moduleID) {
373 Module* moduleCal = GetModule(planeID, moduleID);
374 if (moduleCal == nullptr) return 0; // return numeric_limits<double>::quiet_NaN()
375 return moduleCal->GetSlopeFullSpc();
376}
377
383double TRestDataSetGainMap::GetInterceptParameter(const int planeID, const int moduleID, const double x,
384 const double y) {
385 Module* moduleCal = GetModule(planeID, moduleID);
386 if (moduleCal == nullptr) return 0; // return numeric_limits<double>::quiet_NaN()
387 return moduleCal->GetIntercept(x, y);
388}
389
394double TRestDataSetGainMap::GetInterceptParameterFullSpc(const int planeID, const int moduleID) {
395 Module* moduleCal = GetModule(planeID, moduleID);
396 if (moduleCal == nullptr) return 0; // return numeric_limits<double>::quiet_NaN()
397 return moduleCal->GetInterceptFullSpc();
398}
399
404std::set<int> TRestDataSetGainMap::GetPlaneIDs() const {
405 std::set<int> planeIDs;
406 for (const auto& mc : fModulesCal) planeIDs.insert(mc.GetPlaneId());
407 return planeIDs;
408}
409
414std::set<int> TRestDataSetGainMap::GetModuleIDs(const int planeId) const {
415 std::set<int> moduleIDs;
416 for (const auto& mc : fModulesCal)
417 if (mc.GetPlaneId() == planeId) moduleIDs.insert(mc.GetModuleId());
418 return moduleIDs;
419}
420
424std::map<int, std::set<int>> TRestDataSetGainMap::GetModuleIDs() const {
425 std::map<int, std::set<int>> moduleIds;
426 for (const int planeId : GetPlaneIDs())
427 moduleIds.insert(std::pair<int, std::set<int>>(planeId, GetModuleIDs(planeId)));
428 return moduleIds;
429}
430
431TRestDataSetGainMap& TRestDataSetGainMap::operator=(TRestDataSetGainMap& src) {
432 SetName(src.GetName());
433 fOutputFileName = src.GetOutputFileName();
434 fObservable = src.GetObservable();
435 fSpatialObservableX = src.GetSpatialObservableX();
436 fSpatialObservableY = src.GetSpatialObservableY();
437 fCut = src.GetCut();
438 fModulesCal.clear();
439 for (auto pID : src.GetPlaneIDs())
440 for (auto mID : src.GetModuleIDs(pID)) fModulesCal.push_back(*src.GetModule(pID, mID));
441 return *this;
442}
443
449 for (auto& i : fModulesCal) {
450 if (i.GetPlaneId() == moduleCal.GetPlaneId() && i.GetModuleId() == moduleCal.GetModuleId()) {
451 i = moduleCal;
452 return;
453 }
454 }
455 fModulesCal.push_back(moduleCal);
456}
457
462void TRestDataSetGainMap::Import(const std::string& fileName) {
463 if (fileName.empty()) {
464 RESTError << "No input calibration file defined" << RESTendl;
465 return;
466 }
467
468 if (TRestTools::isRootFile(fileName)) {
469 RESTInfo << "Opening " << fileName << RESTendl;
470 TFile* f = TFile::Open(fileName.c_str(), "READ");
471 if (f == nullptr) {
472 RESTError << "Cannot open calibration file " << fileName << RESTendl;
473 return;
474 }
475
476 TRestDataSetGainMap* cal = nullptr;
477 if (f != nullptr) {
478 TIter nextkey(f->GetListOfKeys());
479 TKey* key;
480 while ((key = (TKey*)nextkey())) {
481 std::string kName = key->GetClassName();
482 if (REST_Reflection::GetClassQuick(kName.c_str()) != nullptr &&
483 REST_Reflection::GetClassQuick(kName.c_str())->InheritsFrom("TRestDataSetGainMap")) {
484 cal = f->Get<TRestDataSetGainMap>(key->GetName());
485 *this = *cal;
486 }
487 }
488 }
489 } else
490 RESTError << "File extension not supported for " << fileName << RESTendl;
492}
493
500void TRestDataSetGainMap::Export(const std::string& fileName) {
501 if (!fileName.empty()) fOutputFileName = fileName;
502 if (fOutputFileName.empty()) {
503 RESTError << "No output file defined" << RESTendl;
504 return;
505 }
506
508 TFile* f = TFile::Open(fOutputFileName.c_str(), "UPDATE");
509 this->Write(GetName());
510 f->Close();
511 delete f;
512 RESTInfo << "Calibration saved to " << fOutputFileName << RESTendl;
513 } else
514 RESTError << "File extension for " << fOutputFileName << "is not supported." << RESTendl;
515}
516
522 RESTMetadata << " Calibration dataset: " << fCalibFileName << RESTendl;
523 if (fCut) {
524 RESTMetadata << " Cuts applied: ";
525 /* print only cutStrings and paramCut because
526 TRestDataSet::MakeCut() uses fCut->GetCutStrings() and fCut->GetParamCut() */
527 for (const auto& cut : fCut->GetCutStrings()) RESTMetadata << cut << ", " << RESTendl;
528 for (const auto& cut : fCut->GetParamCut())
529 RESTMetadata << cut.first << " " << cut.second << ", " << RESTendl;
530 }
531 RESTMetadata << " Output file: " << fOutputFileName << RESTendl;
532 RESTMetadata << RESTendl;
533 RESTMetadata << " Number of planes: " << GetNumberOfPlanes() << RESTendl;
534 RESTMetadata << " Number of modules: " << GetNumberOfModules() << RESTendl;
535 RESTMetadata << " Calibration observable: " << fObservable << RESTendl;
536 RESTMetadata << " Spatial observable X: " << fSpatialObservableX << RESTendl;
537 RESTMetadata << " Spatial observable Y: " << fSpatialObservableY << RESTendl;
538 RESTMetadata << "-----------------------------------------------" << RESTendl;
539 for (auto& i : fModulesCal) i.Print();
540 RESTMetadata << "***********************************************" << RESTendl;
541 RESTMetadata << RESTendl;
542}
543
551std::pair<int, int> TRestDataSetGainMap::Module::GetIndexMatrix(const double x, const double y) const {
552 int index_x = -1, index_y = -1;
553
554 if (fSplitX.upper_bound(x) != fSplitX.end()) {
555 index_x = std::distance(fSplitX.begin(), fSplitX.upper_bound(x)) - 1;
556 if (index_x < 0) {
557 RESTWarning << "index_x < 0 for x = " << x << " and fSplitX[0]=" << *fSplitX.begin()
558 << p->RESTendl;
559 index_x = 0;
560 }
561 } else {
562 RESTWarning << "x is out of split for x = " << x << p->RESTendl;
563 index_x = fSplitX.size() - 2;
564 }
565
566 if (fSplitY.upper_bound(y) != fSplitY.end()) {
567 index_y = std::distance(fSplitY.begin(), fSplitY.upper_bound(y)) - 1;
568 if (index_y < 0) {
569 RESTWarning << "index_y < 0 for y = " << y << " and fSplitY[0]=" << *fSplitY.begin()
570 << p->RESTendl;
571 index_y = 0;
572 }
573 } else {
574 RESTWarning << "y is out of split for y = " << y << p->RESTendl;
575 index_y = fSplitY.size() - 2;
576 }
577
578 return std::make_pair(index_x, index_y);
579}
587double TRestDataSetGainMap::Module::GetSlope(const double x, const double y) const {
588 auto [index_x, index_y] = GetIndexMatrix(x, y);
589 if (fSlope.empty()) {
590 RESTError << "Calibration slope matrix is empty. Returning 0" << p->RESTendl;
591 return 0;
592 }
593
594 if (index_x > (int)fSlope.size() || index_y > (int)fSlope.at(0).size()) {
595 RESTError << "Index out of range. Returning 0" << p->RESTendl;
596 return 0;
597 }
598
599 return fSlope[index_x][index_y];
600}
601
609double TRestDataSetGainMap::Module::GetIntercept(const double x, const double y) const {
610 auto [index_x, index_y] = GetIndexMatrix(x, y);
611 if (fIntercept.empty()) {
612 RESTError << "Calibration constant matrix is empty. Returning 0" << p->RESTendl;
613 return 0;
614 }
615
616 if (index_x > (int)fIntercept.size() || index_y > (int)fIntercept.at(0).size()) {
617 RESTError << "Index out of range. Returning 0" << p->RESTendl;
618 return 0;
619 }
620
621 return fIntercept[index_x][index_y];
622}
623
628 SetSplitX();
629 SetSplitY();
630}
638 if (fNumberOfSegmentsX < 1) {
639 RESTError << "SetSplitX: fNumberOfSegmentsX must be >=1." << p->RESTendl;
640 return;
641 }
642 std::set<double> split;
643 for (int i = 0; i <= fNumberOfSegmentsX; i++) { // <= so that the last segment is included
644 double x =
645 fReadoutRange.X() + ((fReadoutRange.Y() - fReadoutRange.X()) / (float)fNumberOfSegmentsX) * i;
646 split.insert(x);
647 }
648 SetSplitX(split);
649}
650
651void TRestDataSetGainMap::Module::SetSplitX(const std::set<double>& splitX) {
652 if (splitX.size() < 2) {
653 RESTError << "SetSplitX: split size must be >=2 (start and end of range must be included)."
654 << p->RESTendl;
655 return;
656 }
657 if (!fSlope.empty())
658 RESTWarning << "SetSplitX: changing split but current gain map and calibration paremeters correspond "
659 "to previous splitting. Use GenerateGainMap() to update them."
660 << p->RESTendl;
661 fSplitX = splitX;
662 fNumberOfSegmentsX = fSplitX.size() - 1;
663}
671 if (fNumberOfSegmentsY < 1) {
672 RESTError << "SetSplitY: fNumberOfSegmentsY must be >=1." << p->RESTendl;
673 return;
674 }
675 std::set<double> split;
676 for (int i = 0; i <= fNumberOfSegmentsY; i++) { // <= so that the last segment is included
677 double y =
678 fReadoutRange.X() + ((fReadoutRange.Y() - fReadoutRange.X()) / (float)fNumberOfSegmentsY) * i;
679 split.insert(y);
680 }
681 SetSplitY(split);
682}
683
684void TRestDataSetGainMap::Module::SetSplitY(const std::set<double>& splitY) {
685 if (splitY.size() < 2) {
686 RESTError << "SetSplitY: split size must be >=2 (start and end of range must be included)."
687 << p->RESTendl;
688 return;
689 }
690 if (!fSlope.empty())
691 RESTWarning << "SetSplitY: changing split but current gain map and calibration paremeters correspond "
692 "to previous splitting. Use GenerateGainMap() to update them."
693 << p->RESTendl;
694 fSplitY = splitY;
695 fNumberOfSegmentsY = fSplitY.size() - 1;
696}
697
715 //--- Initial checks and settings ---
716 std::string dsFileName = fDataSetFileName;
717 if (dsFileName.empty()) dsFileName = p->GetCalibrationFileName();
718 if (dsFileName.empty()) {
719 RESTError << "No calibration file defined" << p->RESTendl;
720 return;
721 }
722
723 TRestDataSet dataSet;
724 dataSet.EnableMultiThreading(true);
725 if (TRestTools::isDataSet(dsFileName)) {
726 dataSet.Import(dsFileName);
727 fDataSetFileName = dsFileName;
728 } else {
729 RESTWarning << dsFileName << " is not a dataset. Generating a temporal one..." << p->RESTendl;
730 // get all the observables needed for the gain map
731 std::vector<std::string> obsList;
732
733 obsList.push_back(p->GetObservable());
734 obsList.push_back(p->GetSpatialObservableX());
735 obsList.push_back(p->GetSpatialObservableY());
736
737 // look for observables (characterized by having a _ in the name) in the definition cut
738 auto modDefCutObs = TRestTools::GetObservablesInString(fDefinitionCut, true);
739 obsList.insert(obsList.end(), modDefCutObs.begin(), modDefCutObs.end());
740
741 // look for observables in the cut
742 for (const auto& cut : p->GetCut()->GetCutStrings()) {
743 auto cutObs = TRestTools::GetObservablesInString(cut, true);
744 obsList.insert(obsList.end(), cutObs.begin(), cutObs.end());
745 }
746 for (const auto& [variable, condition] : p->GetCut()->GetParamCut()) {
747 auto cutObs = TRestTools::GetObservablesInString(variable, true);
748 obsList.insert(obsList.end(), cutObs.begin(), cutObs.end());
749 // not sure if any obs can be in the condition. Just in case...
750 cutObs = TRestTools::GetObservablesInString(condition, true);
751 obsList.insert(obsList.end(), cutObs.begin(), cutObs.end());
752 }
753
754 // remove duplicates
755 std::sort(obsList.begin(), obsList.end());
756 obsList.erase(std::unique(obsList.begin(), obsList.end()), obsList.end());
757
758 // generate the dataset with the needed observables
759 dataSet.SetFilePattern(dsFileName);
760 dataSet.SetObservablesList(obsList);
761 dataSet.GenerateDataSet();
762 fDataSetFileName = dsFileName;
763 }
764
765 dataSet.SetDataFrame(dataSet.MakeCut(p->GetCut()));
766
767 if (fSplitX.empty()) SetSplitX();
768 if (fSplitY.empty()) SetSplitY();
769
770 //--- Get the calibration range if not provided (default is 0,0) ---
771 if (fCalibRange.X() >= fCalibRange.Y()) {
772 // Get spectrum for this file
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()));
777 // double xMin = hpunt->GetXaxis()->GetXmin();
778 double xMax = hpunt->GetXaxis()->GetXmax();
779
780 // Reduce the range to avoid the possible empty (nCounts<1%) end part of the spectrum
781 double fraction = 1, nAtEndSpc = 0, step = 0.66;
782 while (nAtEndSpc * 1. / hpunt->Integral() < 0.01 && fraction > 0.001) {
783 fraction *= step;
784 nAtEndSpc = hpunt->Integral(hpunt->FindFixBin(hpunt->GetXaxis()->GetXmax() * fraction),
785 hpunt->FindFixBin(hpunt->GetXaxis()->GetXmax()));
786 }
787 xMax = hpunt->GetXaxis()->GetXmax() * fraction /
788 step; // previous step is the last one that nAtEndSpc<1%
789 // Set the calibration range if needed
790 // fCalibRange.SetX(xMin);
791 fCalibRange.SetY(xMax);
792 hpunt.reset(); // delete hpunt;
793 RESTDebug << "Calibration range (auto)set to (" << fCalibRange.X() << "," << fCalibRange.Y() << ")"
794 << p->RESTendl;
795 }
796
797 // --- Definition of histogram whole module ---
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());
801
802 // build the spectrum for the whole module
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());
809
810 //--- Definition of histogram matrix ---
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(),
816 fCalibRange.Y()); // h[column][row] equivalent to h[x][y]
817 }
818 }
819
820 // build the spectrum for each segment
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++) {
825 // Get the segment limits from the splits
826 auto xLower = *itX;
827 auto xUpper = *std::next(itX);
828 auto yLower = *itY;
829 auto yUpper = *std::next(itY);
830
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;
841 auto histo = dataSet.GetDataFrame()
842 .Filter(segment_cut)
843 .Histo1D({"temp", "", h[i][j]->GetNbinsX(), h[i][j]->GetXaxis()->GetXmin(),
844 h[i][j]->GetXaxis()->GetXmax()},
845 GetObservable());
846 std::unique_ptr<TH1F> hpunt = std::unique_ptr<TH1F>(static_cast<TH1F*>(histo->Clone()));
847 h[i][j]->Add(hpunt.get());
848 hpunt.reset(); // delete hpunt;
849 itY++;
850 }
851 itX++;
852 }
853
854 //--- Fit every peak energy for every segment ---
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;
866 }
867 fSlope = calParSlope;
868 fIntercept = calParIntercept;
869 fSegSpectra = h;
870
871 //--- Fit every peak energy for the whole module ---
872 delete fFullLinearFit;
873 fFullLinearFit = new TGraph();
874 auto [intercept, slope] = FitPeaks(fFullSpectrum, fFullLinearFit);
875 fFullSlope = slope;
876 fFullIntercept = intercept;
877}
878
879std::pair<double, double> TRestDataSetGainMap::Module::FitPeaks(TH1F* hSeg, TGraph* gr) {
880 if (!hSeg) {
881 RESTError << "No histogram for fitting" << p->RESTendl;
882 return std::make_pair(0, 0);
883 }
884 if (hSeg->Integral() == 0) {
885 RESTError << "Empty spectrum " << hSeg->GetName() << p->RESTendl;
886 return std::make_pair(0, 0);
887 }
888 std::shared_ptr<TGraph> graph = std::shared_ptr<TGraph>(new TGraph());
889 RESTExtreme << "Fitting peaks for " << hSeg->GetName() << p->RESTendl;
890
891 // Search for peaks --> peakPos
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>());
897 const double ratio =
898 peakPos.size() == 0 ? 1 : peakPos.front() / fEnergyPeaks.front(); // to estimate peak position
899
900 // Initialize graph for linear fit
901 graph->SetName("grFit");
902 graph->SetTitle((";" + GetObservable() + ";energy").c_str());
903
904 // Fit every energy peak
905 int c = 0;
906 double mu = 0;
907 for (const auto& energy : fEnergyPeaks) {
908 RESTExtreme << "\t fitting energy " << DoubleToString(energy, "%g") << p->RESTendl;
909 // estimation of the peak position is between start and end
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()) { // if range is defined use it
914 start = fRangePeaks.at(c).X();
915 end = fRangePeaks.at(c).Y();
916 }
917
918 do {
919 if (fAutoRangePeaks) {
920 if (peakPos.size() > 0) {
921 // Find the peak position that is between start and end
922 pos = peakPos.at(0);
923 while (!(start < pos && pos < end)) {
924 // if none of the peak position is
925 // between start and end, use the greater one.
926 if (pos == peakPos.back()) {
927 pos = peakPos.at(0);
928 break;
929 }
930 pos = *std::next(std::find(peakPos.begin(), peakPos.end(),
931 pos)); // get next peak position
932 }
933 peakPos.erase(std::find(peakPos.begin(), peakPos.end(),
934 pos)); // remove this peak position from the list
935 // final estimation of the peak range (idem fitting window) with this peak
936 // position pos
937 start = pos * 0.8;
938 end = pos * 1.2;
939 const double relDist = peakPos.size() > 0 ? (pos - peakPos.front()) / pos : 999;
940 if (relDist < 0.2) { // if the next peak is too close reduce the window width
941 start = pos * (1 - relDist / 2);
942 end = pos * (1 + relDist / 2);
943 }
944 }
945 }
946
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("
950 << DoubleToString(start, "%.3g") << ", " << DoubleToString(end, "%.3g") << ")"
951 << p->RESTendl;
952
953 if (hSeg->GetFunction(name.c_str())) // remove previous fit
954 hSeg->GetListOfFunctions()->Remove(hSeg->GetFunction(name.c_str()));
955
956 hSeg->Fit(g, "R+Q0"); // use 0 to not draw the fit but save it
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)); // avoid small peaks on main peak tail
961 graph->SetPoint(c++, mu, energy);
962 }
963 s.reset(); // delete s;
964
965 if (fZeroPoint) graph->SetPoint(c++, 0, 0);
966 while (graph->GetN() < 2) { // minimun 2 points needed for linear fit
967 graph->SetPoint(c++, 0, 0);
968 SetZeroPoint(true);
969 RESTDebug << "Not enough points for linear fit. Adding and setting zero point to true" << p->RESTendl;
970 }
971
972 // Linear fit
973 std::unique_ptr<TF1> linearFit;
974 linearFit = std::unique_ptr<TF1>(new TF1("linearFit", "pol1"));
975 graph->Fit("linearFit", "SQ"); // Q for quiet mode
976
977 if (gr) *gr = *(TGraph*)graph->Clone(); // if nullptr is passed, do not copy the graph
978 return std::make_pair(linearFit->GetParameter(0), linearFit->GetParameter(1));
979}
980
988void TRestDataSetGainMap::Module::Refit(const TVector2& position, const double energyPeak,
989 const TVector2& range) {
990 auto [index_x, index_y] = GetIndexMatrix(position.X(), position.Y());
991 int peakNumber = -1;
992 for (size_t i = 0; i < fEnergyPeaks.size(); i++)
993 if (fEnergyPeaks.at(i) == energyPeak) {
994 peakNumber = i;
995 break;
996 }
997 if (peakNumber == -1) {
998 RESTError << "Energy " << energyPeak << " not found in the list of energy peaks" << p->RESTendl;
999 return;
1000 }
1001 Refit((size_t)index_x, (size_t)index_y, (size_t)peakNumber, range);
1002}
1003
1013void TRestDataSetGainMap::Module::Refit(const size_t x, const size_t y, const size_t peakNumber,
1014 const TVector2& range) {
1015 if (fSegSpectra.empty()) {
1016 RESTError << "No gain map found. Use GenerateGainMap() first." << p->RESTendl;
1017 return;
1018 }
1019 if (x >= fSegSpectra.size() || y >= fSegSpectra.at(0).size()) {
1020 RESTError << "Segment with index (" << x << ", " << y << ") not found" << p->RESTendl;
1021 return;
1022 }
1023 if (peakNumber >= fEnergyPeaks.size()) {
1024 RESTError << "Peak with index " << peakNumber << " not found" << p->RESTendl;
1025 return;
1026 }
1027
1028 // Refit the desired peak
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())) // clear previous fits for this peakNumber
1033 h->GetListOfFunctions()->Remove(h->GetFunction(name.c_str()));
1034 h->Fit(g, "R+Q0"); // use 0 to not draw the fit but save it
1035
1036 // Change the point of the graph
1037 UpdateCalibrationFits(x, y);
1038}
1039
1047void TRestDataSetGainMap::Module::RefitFullSpc(const double energyPeak, const TVector2& range) {
1048 int peakNumber = -1;
1049 for (size_t i = 0; i < fEnergyPeaks.size(); i++)
1050 if (fEnergyPeaks.at(i) == energyPeak) {
1051 peakNumber = i;
1052 break;
1053 }
1054 if (peakNumber == -1) {
1055 RESTError << "Energy " << energyPeak << " not found in the list of energy peaks" << p->RESTendl;
1056 return;
1057 }
1058 RefitFullSpc((size_t)peakNumber, range);
1059}
1060
1068void TRestDataSetGainMap::Module::RefitFullSpc(const size_t peakNumber, const TVector2& range) {
1069 if (!fFullSpectrum) {
1070 RESTError << "No gain map found. Use GenerateGainMap() first." << p->RESTendl;
1071 return;
1072 }
1073 if (peakNumber >= fEnergyPeaks.size()) {
1074 RESTError << "Peak with index " << peakNumber << " not found" << p->RESTendl;
1075 return;
1076 }
1077
1078 // Refit the desired peak
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())) // clear previous fits for this peakNumber
1082 fFullSpectrum->GetListOfFunctions()->Remove(fFullSpectrum->GetFunction(name.c_str()));
1083 fFullSpectrum->Fit(g, "R+Q0"); // use 0 to not draw the fit but save it
1084
1085 // Change the point of the graph
1086 UpdateCalibrationFitsFullSpc();
1087}
1088
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;
1100 return;
1101 }
1102 if (x >= fSegSpectra.size() || y >= fSegSpectra.at(0).size()) {
1103 RESTError << "Segment with index (" << x << ", " << y << ") not found" << p->RESTendl;
1104 return;
1105 }
1106
1107 TH1F* h = fSegSpectra.at(x).at(y);
1108 TGraph* gr = fSegLinearFit.at(x).at(y);
1109
1110 auto [intercept, slope] = UpdateCalibrationFits(h, gr);
1111 fSlope[x][y] = slope;
1112 fIntercept[x][y] = intercept;
1113}
1114
1115std::pair<double, double> TRestDataSetGainMap::Module::UpdateCalibrationFits(TH1F* h, TGraph* gr) {
1116 if (!h) {
1117 RESTError << "No histogram for updating fits" << p->RESTendl;
1118 return std::make_pair(0, 0);
1119 }
1120 if (!gr) {
1121 RESTError << "No graph for updating fits" << p->RESTendl;
1122 return std::make_pair(0, 0);
1123 }
1124 if (h->Integral() == 0) {
1125 RESTError << "Empty spectrum " << h->GetName() << p->RESTendl;
1126 return std::make_pair(0, 0);
1127 }
1128
1129 // Clear the points of the graph
1130 for (size_t i = 0; i < fEnergyPeaks.size(); i++) gr->RemovePoint(i);
1131 // Add the new points to the graph
1132 int c = 0;
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());
1136 if (!g) {
1137 RESTWarning << "No fit ( " << fitName << " ) found for energy peak " << fEnergyPeaks[i]
1138 << " in histogram " << h->GetName() << p->RESTendl;
1139 continue;
1140 }
1141 gr->SetPoint(c++, g->GetParameter(1), fEnergyPeaks[i]);
1142 }
1143
1144 // Add zero points if needed (if there are less than 2 points)
1145 while (gr->GetN() < 2) {
1146 gr->SetPoint(c++, 0, 0);
1147 }
1148
1149 // Refit the calibration curve
1150 TF1* lf = nullptr;
1151 if (gr->GetFunction("linearFit"))
1152 lf = gr->GetFunction("linearFit");
1153 else
1154 lf = new TF1("linearFit", "pol1");
1155 gr->Fit(lf, "SQ"); // Q for quiet mode
1156
1157 return std::make_pair(lf->GetParameter(0), lf->GetParameter(1));
1158}
1159
1166 auto [intercept, slope] = UpdateCalibrationFits(fFullSpectrum, fFullLinearFit);
1167 fFullSlope = slope;
1168 fFullIntercept = intercept;
1169}
1170
1186 if (module == nullptr) {
1187 RESTError << "TRestDataSetGainMap::Module::LoadConfigFromTiXmlElement: module is nullptr"
1188 << p->RESTendl;
1189 return;
1190 }
1191
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));
1196
1197 el = !module->Attribute("moduleDefinitionCut") ? "Not defined" : module->Attribute("moduleDefinitionCut");
1198 if (!(el.empty() || el == "Not defined")) this->SetModuleDefinitionCut(el);
1199
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));
1206
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));
1211
1212 el = !module->Attribute("dataSetFileName") ? "Not defined" : module->Attribute("dataSetFileName");
1213 if (!(el.empty() || el == "Not defined")) this->SetDataSetFileName(el);
1214
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");
1219
1220 // Get peaks energy and range
1221 TiXmlElement* peakDefinition = (TiXmlElement*)module->FirstChildElement("peak");
1222 while (peakDefinition != nullptr) {
1223 double energy = 0;
1224 TVector2 range = TVector2(0, 0);
1225
1226 std::string ell =
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;
1230 exit(1);
1231 }
1232 energy = StringToDouble(ell);
1233
1234 ell = !peakDefinition->Attribute("range") ? "Not defined" : peakDefinition->Attribute("range");
1235 if (!(ell.empty() || ell == "Not defined")) range = StringTo2DVector(ell);
1236
1237 this->AddPeak(energy, range);
1238 peakDefinition = (TiXmlElement*)peakDefinition->NextSiblingElement();
1239 }
1240}
1241
1242void TRestDataSetGainMap::Module::DrawSpectrum(const TVector2& position, bool drawFits, int color,
1243 TCanvas* c) {
1244 std::pair<size_t, size_t> index = GetIndexMatrix(position.X(), position.Y());
1245 DrawSpectrum(index.first, index.second, drawFits, color, c);
1246}
1247
1248void TRestDataSetGainMap::Module::DrawSpectrum(const int index_x, const int index_y, bool drawFits, int color,
1249 TCanvas* c) {
1250 if (fSegSpectra.size() == 0) {
1251 RESTError << "Spectra matrix is empty." << p->RESTendl;
1252 return;
1253 }
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;
1257 return;
1258 }
1259 if (!fSegSpectra[index_x][index_y]) {
1260 RESTError << "No Spectrum for segment (" << index_x << ", " << index_y << ")." << p->RESTendl;
1261 return;
1262 }
1263
1264 if (!c) {
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());
1268 }
1269
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);
1274 std::string tH = "Spectrum x=[" + DoubleToString(xLower, "%g") + ", " + DoubleToString(xUpper, "%g") +
1275 ") y=[" + DoubleToString(yLower, "%g") + ", " + DoubleToString(yUpper, "%g") + ");" +
1276 GetObservable() + ";counts";
1277 fSegSpectra[index_x][index_y]->SetTitle(tH.c_str());
1278
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");
1282
1283 if (drawFits)
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());
1286 if (!fit)
1287 RESTWarning << "Fit for energy peak " << fEnergyPeaks[c] << " not found." << p->RESTendl;
1288 if (!fit) continue;
1289 fit->SetLineColor(c + 2 != colorT ? c + 2 : ++colorT); /* does not work with kRed, kBlue, etc.
1290 as they are not defined with the same number as the first 10 basic colors. See
1291 https://root.cern.ch/doc/master/classTColor.html#C01 and
1292 https://root.cern.ch/doc/master/classTColor.html#C02 */
1293 fit->Draw("same");
1294 }
1295}
1296
1316void TRestDataSetGainMap::Module::DrawSpectrum(const bool drawFits, const int color, TCanvas* c) {
1317 if (fSegSpectra.size() == 0) {
1318 RESTError << "Spectra matrix is empty." << p->RESTendl;
1319 return;
1320 }
1321 if (!c) {
1322 std::string t = "spectrum_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId);
1323 c = new TCanvas(t.c_str(), t.c_str());
1324 }
1325
1326 size_t nPads = 0;
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;
1332 return;
1333 } else if (nPads == 0)
1334 c->Divide(fSegSpectra.size(), fSegSpectra.at(0).size());
1335
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;
1339 c->cd(pad);
1340 DrawSpectrum(i, j, drawFits, color, c);
1341 }
1342 }
1343}
1344void TRestDataSetGainMap::Module::DrawFullSpectrum(const bool drawFits, const int color, TCanvas* c) {
1345 if (!fFullSpectrum) {
1346 RESTError << "Spectrum is empty." << p->RESTendl;
1347 return;
1348 }
1349
1350 if (!c) {
1351 std::string t = "fullSpc_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId);
1352 c = new TCanvas(t.c_str(), t.c_str());
1353 }
1354 c->cd();
1355
1356 fFullSpectrum->SetTitle(("Full spectrum;" + GetObservable() + ";counts").c_str());
1357
1358 if (color > 0) fFullSpectrum->SetLineColor(color);
1359 size_t colorT = fFullSpectrum->GetLineColor();
1360 fFullSpectrum->Draw("same");
1361
1362 if (drawFits)
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;
1366 if (!fit) continue;
1367 fit->SetLineColor(c + 2 != colorT ? c + 2 : ++colorT); /* does not work with kRed, kBlue, etc.
1368 as they are not defined with the same number as the first 10 basic colors. See
1369 https://root.cern.ch/doc/master/classTColor.html#C01 and
1370 https://root.cern.ch/doc/master/classTColor.html#C02 */
1371 fit->Draw("same");
1372 }
1373}
1374
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);
1378}
1379
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;
1383 return;
1384 }
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;
1388 return;
1389 }
1390 if (!fSegLinearFit[index_x][index_y]) {
1391 RESTError << "No linear fit for segment (" << index_x << ", " << index_y << ")." << p->RESTendl;
1392 return;
1393 }
1394
1395 if (!c) {
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());
1399 }
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);
1404 std::string tH = "Linear Fit x=[" + DoubleToString(xLower, "%g") + ", " + DoubleToString(xUpper, "%g") +
1405 ") y=[" + DoubleToString(yLower, "%g") + ", " + DoubleToString(yUpper, "%g") + ");" +
1406 GetObservable() + ";energy";
1407 fSegLinearFit[index_x][index_y]->SetTitle(tH.c_str());
1408 fSegLinearFit[index_x][index_y]->Draw("AL*");
1409}
1410
1411void TRestDataSetGainMap::Module::DrawLinearFit(TCanvas* c) {
1412 if (fSegLinearFit.size() == 0) {
1413 RESTError << "Spectra matrix is empty." << p->RESTendl;
1414 return;
1415 }
1416 if (!c) {
1417 std::string t = "linearFits_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId);
1418 c = new TCanvas(t.c_str(), t.c_str());
1419 }
1420
1421 size_t nPads = 0;
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;
1427 return;
1428 } else if (nPads == 0)
1429 c->Divide(fSegLinearFit.size(), fSegLinearFit.at(0).size());
1430
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;
1434 c->cd(pad);
1435 DrawLinearFit(i, j, c);
1436 }
1437 }
1438}
1439
1448void TRestDataSetGainMap::Module::DrawGainMap(const int peakNumber, bool fullModuleAsRef) {
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;
1452 return;
1453 }
1454 if (fSegLinearFit.size() == 0) {
1455 RESTError << "Linear fit matrix is empty." << p->RESTendl;
1456 return;
1457 }
1458 if (!fFullLinearFit) {
1459 RESTError << "Full linear fit is empty." << p->RESTendl;
1460 return;
1461 }
1462
1463 double peakEnergy = fEnergyPeaks[peakNumber];
1464 std::string title = "Gain map for energy " + DoubleToString(peakEnergy, "%g") + ";" +
1465 GetSpatialObservableX() + ";" + GetSpatialObservableY(); // + " keV";
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());
1469 gainMap->cd();
1470 TH2F* hGainMap = new TH2F(("h" + t).c_str(), title.c_str(), fNumberOfSegmentsX, fReadoutRange.X(),
1471 fReadoutRange.Y(), fNumberOfSegmentsY, fReadoutRange.X(), fReadoutRange.Y());
1472
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);
1478 }
1479
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++) {
1484 auto xLower = *itX;
1485 auto xUpper = *std::next(itX);
1486 auto yLower = *itY;
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);
1493 itY++;
1494 }
1495 itX++;
1496 }
1497 hGainMap->SetStats(0);
1498 hGainMap->Draw("colz");
1499 hGainMap->SetBarOffset(0.2);
1500 hGainMap->Draw("TEXT SAME");
1501}
1502
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;
1513
1514 RESTMetadata << " Calibration dataset: " << fDataSetFileName << p->RESTendl;
1515 RESTMetadata << p->RESTendl;
1516
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: ";
1524 anyRange = true;
1525 break;
1526 }
1527 if (anyRange)
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() << " )"
1533 << p->RESTendl;
1534 RESTMetadata << " Number of bins: " << fNBins << p->RESTendl;
1535 RESTMetadata << p->RESTendl;
1536
1537 RESTMetadata << " Number of segments X: " << fNumberOfSegmentsX << p->RESTendl;
1538 RESTMetadata << " Number of segments Y: " << fNumberOfSegmentsY << p->RESTendl;
1539 // RESTMetadata << " Draw: " << (fDrawVar ? "true" : "false") << p->RESTendl;
1540 RESTMetadata << " Readout range (" << fReadoutRange.X() << ", " << fReadoutRange.Y() << " )"
1541 << p->RESTendl;
1542 RESTMetadata << "SplitX: ";
1543 for (auto& i : fSplitX) {
1544 RESTMetadata << " " << i;
1545 }
1546 RESTMetadata << p->RESTendl;
1547 RESTMetadata << "SplitY: ";
1548 for (auto& i : fSplitY) {
1549 RESTMetadata << " " << i;
1550 }
1551 RESTMetadata << p->RESTendl;
1552 RESTMetadata << p->RESTendl;
1553
1554 RESTMetadata << " Slope: " << p->RESTendl;
1555 size_t maxSize = 0;
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") << " ";
1562 else
1563 RESTMetadata << " ";
1564 }
1565 RESTMetadata << p->RESTendl;
1566 }
1567 RESTMetadata << " Intercept: " << p->RESTendl;
1568 maxSize = 0;
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") << " ";
1575 else
1576 RESTMetadata << " ";
1577 }
1578 RESTMetadata << p->RESTendl;
1579 }
1580 RESTMetadata << p->RESTendl;
1581 RESTMetadata << " Full slope: " << DoubleToString(fFullSlope, "%.3e") << p->RESTendl;
1582 RESTMetadata << " Full intercept: " << DoubleToString(fFullIntercept, "%+.3e") << p->RESTendl;
1583
1584 RESTMetadata << "-----------------------------------------------" << p->RESTendl;
1585}
A class to help on cuts definitions. To be used with TRestAnalysisTree.
Definition TRestCut.h:31
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...
A base class for any REST metadata class.
virtual void PrintMetadata()
Implemented it in the derived metadata class to print out specific metadata information.
endl_t RESTendl
Termination flag object for TRestStringOutput.
TiXmlElement * GetElement(std::string eleDeclare, TiXmlElement *e=nullptr)
Get an xml element from a given parent element, according to its declaration.
Int_t LoadConfigFromFile(const std::string &configFilename, const std::string &sectionName="")
Give the file name, find out the corresponding section. Then call the main starter.
TRestMetadata * InstantiateChildMetadata(int index, std::string pattern="")
This method will retrieve a new TRestMetadata instance of a child element of the present TRestMetadat...
virtual void InitFromConfigFile()
To make settings from rml file. This method must be implemented in the derived class.
TRestStringOutput::REST_Verbose_Level GetVerboseLevel()
returns the verboselevel in type of REST_Verbose_Level enumerator
void SetSectionName(std::string sName)
set the section name, clear the section content
std::string fConfigFileName
Full name of the rml file.
virtual Int_t Write(const char *name=nullptr, Int_t option=0, Int_t bufsize=0)
overwriting the write() method with fStore considered
TiXmlElement * GetNextElement(TiXmlElement *e)
Get the next sibling xml element of this element, with same eleDeclare.
@ REST_Info
+show most of the information for each steps
@ REST_Debug
+show the defined debug messages
static std::string GetFileNameExtension(const std::string &fullname)
Gets the file extension as the substring found after the latest ".".
static std::set< std::string > GetMatchingStrings(const std::vector< std::string > &stack, const std::vector< std::string > &wantedStrings)
Returns a set of strings that match the wanted strings from the stack. The wanted strings can contain...
static bool isDataSet(const std::string &filename)
It checks if the file contains a dataset object.
static std::vector< std::string > GetObservablesInString(const std::string &observablesStr, bool removeDuplicates=true)
Returns a vector with the observables names found in the input string. The observables names must con...
static bool isRootFile(const std::string &filename)
Returns true if the filename has *.root* extension.
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.