diff --git a/macros/REST_RefitGainMap.C b/macros/REST_RefitGainMap.C new file mode 100644 index 000000000..b1532168b --- /dev/null +++ b/macros/REST_RefitGainMap.C @@ -0,0 +1,497 @@ +////////////////////////////////////////////////////////// +/// Macro to refit the gain map with a simple 'GUI' to allow +/// the user to add or delete peaks and update the fits with +/// the FIT PANEL. This macro is meant to be used with the +/// FIT PANEL. The GUI consists of a dialog canvas and two +/// graphical canvases, one with all the segments and another +/// one with the segment selected alone. The dialog canvas +/// contains a button for each module. This will draw the +/// segments of the module in the first canvas (cAll). Then +/// the user can select the segment to be shown in the second +/// canvas by clicking on the first one. When the spectrum of +/// the segment is shown in the second canvas (cAlone), the +/// user can add or delete peaks with the buttons 'Add peak' +/// and 'Delete peak'. The user can also update the fits +/// using the FIT PANEL by selecting the fit type 'Prev. Fit' +/// and fit function named 'g'+peakNumber (in energy descending +/// order starting at 0). It is important to check the option +/// 'Add to list' in the FIT PANEL (the draw option 'SAME' is +/// optional but it can be helpful), if not all the fits will +/// be deleted. After fitting with the fit panel, click on the +/// button 'UpdateFits' to update the fits in the canvas and the +/// calibration curve of that segment. +/// Finally, export the gain map with the button 'Export' on the +/// dialog canvas. The user can also export the gain map to a +/// root file with the command 'gm.Export();' in the terminal. +/// +/// \author: Álvaro Ezquerro aezquerro@unizar.es +/// +/// Known bugs: +/// - Sometimes when clicking the drawAlone button, the segment +/// spectrum displayed in cAlone is not the correct one. I have +/// not been able to reproduce it systematically. +/// - Crash when deleting all the peaks of a segment. +/// - Crash (seg fault) sometimes when deleting one of the peaks. +/// Although you can still continue using the macro (but it crashes +/// totally when exiting the root terminal with .q) +/// +/// Possible improvements: +/// - Encapsule the canvases (and fit panel maybe?) in a proper +/// ROOT GUI. +/// +////////////////////////////////////////////////////////// + +// Forward declaration of TFitEditor2 +class TFitEditor2; +typedef std::multimap::iterator fPrevFitIter; + +// Class definition for GainMapRefitter +class GainMapRefitter { + public: + // Member variables (formerly global variables) + TRestDataSetGainMap gm; + TRestDataSetGainMap::Module* m; + std::map meansAux; + std::string fitNameAux; + const int screenWidth; + const int screenHeight; + const int cAllHeight; + const int cAllWidth; + const int cAloneHeight; + const int cAloneWidth; + TCanvas* cAll; + TCanvas* cAlone; + TFitEditor* fitEditor; + TFitEditor2* fitEditor2; + + public: + GainMapRefitter() + : m(nullptr), + screenWidth(gClient->GetDisplayWidth()), + screenHeight(gClient->GetDisplayHeight()), + cAllHeight(2.0 / 3 * screenHeight), + cAllWidth(1.25 * cAllHeight), + cAloneHeight(0.45 * screenHeight), + cAloneWidth(0.8 * cAllHeight), + cAll(new TCanvas("cAll", "cAll", cAllWidth, cAllHeight)), + cAlone(new TCanvas("cAlone", "cAlone", cAloneWidth, cAloneHeight)), + fitEditor(nullptr), + fitEditor2(nullptr) {} + + ~GainMapRefitter() { + delete cAll; + delete cAlone; + // Note: fitEditor and fitEditor2 are managed by ROOT and not deleted here + } + + // Method declarations + void drawAll(); + void drawAlone(const int x, const int y); + void drawWithinAll(int x, int y); + void clearCanvas(TCanvas* c, size_t n_subPad = 0); + void highlightDrawnAlonePad(const int x, const int y); + void DeletePeak(const int x, const int y, const int peakNumber); + void AddPeak(const int x, const int y, const int peakNumber); + void UpdateFits(const int x, const int y); + void changeModule(int plane, int module); + TDialogCanvas* createDialog(); + void ImportGainMap(const std::string& gainMapFile) { gm.Import(gainMapFile); } + + // Getter for gm (if needed externally) + TRestDataSetGainMap& GetGainMap() { return gm; } +}; + +// Class that inherits from TFitEditor to be able to access the protected methods of TFitEditor +class TFitEditor2 : public TFitEditor { + public: + TFitEditor2(TVirtualPad* pad, TObject* obj) : TFitEditor(pad, obj){}; + virtual ~TFitEditor2() {} + + // making public the protected method GetFitObjectListOfFunctions() + TList* GetFitObjectListOfFunctions() { return TFitEditor::GetFitObjectListOfFunctions(); } + TF1* GetFitFunction() { return TFitEditor::GetFitFunction(); } + void CleanFitFunctionList(std::string containingString = "Prev") { + /*// print fPrevFit + std::cout << "Before fPrevFit:" << std::endl; + for (auto it = fPrevFit.begin(); it != fPrevFit.end(); ++it) { + std::cout << it->first->GetName() << " => " << it->second->GetName() << '\n'; + } //*/ + std::pair look = fPrevFit.equal_range(fFitObject); + for (fPrevFitIter it = look.first; it != look.second; ++it) { + std::string name = it->second->GetName(); + if (name.find(containingString) != std::string::npos) { + fPrevFit.erase(it); + break; + } + } + TFitEditor::FillFunctionList(); // update the list of functions in the fit panel + } + + void ToggleNecessaryOptions() { + if (fAdd2FuncList) { + if (fAdd2FuncList->GetState() == kButtonUp) fAdd2FuncList->SetState(kButtonDown); + } + + if (fDrawSame) { + if (fDrawSame->GetState() == kButtonUp) fDrawSame->SetState(kButtonDown); + } + } + + TGComboBox* GetFunctionListComboBox() { return fFuncList; } + + std::string GetSelectedFunctionName() { + if (!fFuncList) return ""; + return fFuncList->GetSelectedEntry()->GetTitle(); + } + + void SelectFunctionByName(const char* name) { + if (!fFuncList) return; + if (fFuncList->FindEntry(name)) + if (fFuncList->GetSelected() != fFuncList->FindEntry(name)->EntryId()) + fFuncList->Select(fFuncList->FindEntry(name)->EntryId()); + } +}; + +// Method definitions +void GainMapRefitter::clearCanvas(TCanvas* c, size_t n_subPad) { + c->cd(n_subPad); + if (n_subPad == 0) { + c->GetListOfPrimitives()->Clear(); + c->Modified(); + c->Range(0, 0, 1, 1); + } else { + TPad* subpad = (TPad*)c->GetPad(n_subPad); + subpad->SetFillColor(0); + subpad->GetListOfPrimitives()->Clear(); + subpad->Modified(); + subpad->Range(0, 0, 1, 1); + } +} + +void GainMapRefitter::drawWithinAll(int x, int y) { + size_t nSubPad = x + 1 + m->GetNumberOfSegmentsX() * (m->GetNumberOfSegmentsY() - y - 1); + if (cAll->GetCanvasImp()) { + clearCanvas(cAll, nSubPad); + } else { + cAll = new TCanvas("cAll", "cAll", cAllWidth, cAllHeight); + cAll->Divide(m->GetNumberOfSegmentsX(), m->GetNumberOfSegmentsY()); + } + cAll->cd(nSubPad); + m->DrawSpectrum((size_t)x, (size_t)y, true, -1, cAll); + + std::string action = (std::string) "gRefitter->drawAlone(" + std::to_string(x) + (std::string) "," + + std::to_string(m->fSegSpectra[x].size() - 1 - y) + (std::string) ");"; + TButton* but = new TButton("Draw alone", action.c_str(), .5, .8, .8, .88); + but->Draw(); + cAll->Update(); +} + +void GainMapRefitter::drawAlone(const int x, const int y) { + for (TObject* obj : *m->fSegSpectra.at(x).at(y)->GetListOfFunctions()) { + if (obj && obj->InheritsFrom(TF1::Class())) { + TF1* g = (TF1*)obj; + meansAux[(std::string)g->GetName()] = g->GetParameter(1); + } + } + + if (cAlone->GetCanvasImp()) { + clearCanvas(cAlone); + } else { + cAlone = new TCanvas("cAlone", "cAlone", cAloneWidth, cAloneHeight); + } + + cAlone->cd(); + std::string canvasTitle = "hSpc_" + std::to_string(m->GetPlaneId()) + "_" + + std::to_string(m->GetModuleId()) + "_" + std::to_string(x) + "_" + + std::to_string(y); + cAlone->SetTitle(canvasTitle.c_str()); + m->DrawSpectrum(x, y, true, -1, cAlone); + + std::string action = "gRefitter->UpdateFits(" + std::to_string(x) + "," + std::to_string(y) + ")"; + TButton* butAlone = new TButton("UpdateFits", action.c_str(), .7, .75, .9, .825); + butAlone->Draw(); + + for (size_t n = 0; n < m->fEnergyPeaks.size(); n++) { + bool peakHasFit = false; + std::string objName = "g" + std::to_string(n); + TF1* f = m->fSegSpectra[x][y]->GetFunction(objName.c_str()); + if (f) { + std::string action = "gRefitter->DeletePeak(" + std::to_string(x) + "," + std::to_string(y) + + "," + std::to_string(n) + ")"; + std::string name = "Delete energy " + DoubleToString(m->fEnergyPeaks.at(n), "%g"); + TButton* but = new TButton(name.c_str(), action.c_str(), .7, .65 - n * .05, .95, .7 - n * .05); + but->Draw(); + peakHasFit = true; + } + if (!peakHasFit) { + std::string action = "gRefitter->AddPeak(" + std::to_string(x) + "," + std::to_string(y) + "," + + std::to_string(n) + ")"; + std::string name = "Add energy " + std::to_string(m->fEnergyPeaks.at(n)); + TButton* but = new TButton(name.c_str(), action.c_str(), .7, .65 - n * .05, .95, .7 - n * .05); + but->Draw(); + } + } + cAlone->Update(); + + highlightDrawnAlonePad(x, y); + if (!fitEditor) { + fitEditor = TFitEditor::GetInstance(cAlone, m->fSegSpectra.at(x).at(y)); + fitEditor2 = (TFitEditor2*)fitEditor; + } + fitEditor2->ToggleNecessaryOptions(); + fitEditor2->CleanFitFunctionList("Prev"); + fitEditor->Show(cAlone, m->fSegSpectra.at(x).at(y)); + fitEditor2->SelectFunctionByName(fitNameAux.c_str()); +} + +void GainMapRefitter::highlightDrawnAlonePad(const int x, const int y) { + uint COLOR = 38; + size_t nSubPad = x + 1 + m->GetNumberOfSegmentsX() * (m->GetNumberOfSegmentsY() - y - 1); + size_t nPads = 0; + for (const auto& object : *cAll->GetListOfPrimitives()) + if (object->InheritsFrom(TVirtualPad::Class())) ++nPads; + + for (size_t i = 0; i < nPads; i++) { + if (i == nSubPad) continue; + TVirtualPad* pad = (TVirtualPad*)cAll->cd(i + 1); + if (pad->GetFillColor() != 0) { + pad->SetFillColor(0); + pad->Modified(); + pad->Update(); + } + } + + if (nSubPad > nPads) { + std::cout << "Error: the number of pads is " << nPads << " and the selected pad is " << nSubPad + << std::endl; + return; + } + TVirtualPad* pad = (TVirtualPad*)cAll->cd(nSubPad); + if (pad->GetFillColor() == COLOR) return; + pad->SetFillColor(COLOR); + pad->Modified(); + pad->Update(); +} + +void GainMapRefitter::DeletePeak(const int x, const int y, const int peakNumber) { + TH1F* h = m->fSegSpectra.at(x).at(y); + TGraph* gr = m->fSegLinearFit.at(x).at(y); + std::string objName = "g" + std::to_string(peakNumber); + TF1* f = h->GetFunction(objName.c_str()); + if (f) { + h->GetListOfFunctions()->Remove(f); + gr->RemovePoint(peakNumber); + fitEditor2->CleanFitFunctionList(objName.c_str()); + } + + drawWithinAll(x, y); + drawAlone(x, y); +} + +void GainMapRefitter::AddPeak(const int x, const int y, const int peakNumber) { + TH1F* h = m->fSegSpectra.at(x).at(y); + std::string objName = "g" + std::to_string(peakNumber); + TF1* g = new TF1(objName.c_str(), "gaus", meansAux[objName] * 0.8, meansAux[objName] * 1.2); + h->Fit(g, "R+Q0"); + + drawWithinAll(x, y); + drawAlone(x, y); +} + +void GainMapRefitter::UpdateFits(const int x, const int y) { + TH1F* h = m->fSegSpectra.at(x).at(y); + TList* list = h->GetListOfFunctions(); + for (size_t n = 0; n < m->fEnergyPeaks.size(); n++) { + std::string objName = "g" + std::to_string(n); + TF1* firstFit = h->GetFunction(objName.c_str()); + TF1* lastFit = nullptr; + for (int i = list->GetSize() - 1; i >= 0; i--) { + TObject* obj = list->At(i); + if (obj && obj->InheritsFrom(TF1::Class()) && obj->GetName() == objName) { + if (!lastFit) + lastFit = (TF1*)list->Remove(obj); + else + list->Remove(obj); + } + } + + if (firstFit && lastFit && firstFit != lastFit) lastFit->Copy(*firstFit); + if (firstFit) h->GetListOfFunctions()->Add(firstFit); + + if (h->GetFunction(objName.c_str())) m->UpdateCalibrationFits(x, y); + } + + drawWithinAll(x, y); + fitNameAux = fitEditor2->GetSelectedFunctionName(); + drawAlone(x, y); + + std::cout << std::endl; + std::cout << "Segment " << x << ", " << y << std::endl; + for (TObject* obj : *h->GetListOfFunctions()) { + if (obj && obj->InheritsFrom(TF1::Class())) { + TF1* f = (TF1*)obj; + std::cout << "Function " << f->GetName() << std::endl; + std::cout << "\tmean : " << f->GetParameter(1) << std::endl; + std::cout << "\tsigma : " << f->GetParameter(2) << std::endl; + } + } + + TGraph* gr = m->fSegLinearFit.at(x).at(y); + if (!gr) return; + std::cout << "Graph " << gr->GetName() << std::endl; + for (int i = 0; i < gr->GetN(); i++) { + double x, y; + gr->GetPoint(i, x, y); + std::cout << "\tPoint " << i << ": " << x << ", " << y << std::endl; + } + for (TObject* obj : *gr->GetListOfFunctions()) { + if (obj && obj->InheritsFrom(TF1::Class())) { + TF1* f = (TF1*)obj; + std::cout << "Function " << f->GetName() << std::endl; + std::cout << "\tslope : " << f->GetParameter(1) << std::endl; + std::cout << "\tintercept : " << f->GetParameter(0) << std::endl; + } + } + std::cout << std::endl; +} + +void GainMapRefitter::drawAll() { + cAll->cd(); + if (cAll->GetCanvasImp()) { + clearCanvas(cAll); + } else { + cAll = new TCanvas("cAll", "cAll", 900, 700); + } + + m->DrawSpectrum(true, -1, cAll); + for (size_t i = 0; i < m->fSegSpectra.size(); i++) { + for (size_t j = 0; j < m->fSegSpectra[i].size(); j++) { + cAll->cd(i + 1 + m->fSegSpectra[i].size() * j); + std::string action = (std::string) "gRefitter->drawAlone(" + std::to_string(i) + + (std::string) "," + std::to_string(m->fSegSpectra[i].size() - 1 - j) + + (std::string) ");"; + TButton* but = new TButton("Draw alone", action.c_str(), .5, .8, .8, .88); + but->Draw(); + } + } + cAll->Update(); +} + +void GainMapRefitter::changeModule(int plane, int module) { + if (m) { + for (size_t i = 0; i < m->fSegSpectra.size(); i++) + for (size_t j = 0; j < m->fSegSpectra[i].size(); j++) m->UpdateCalibrationFits(i, j); + } + m = gm.GetModule(plane, module); + if (!m) return; + + if (cAlone->GetCanvasImp()) { + clearCanvas(cAlone); + cAlone->SetTitle("cAlone"); + cAlone->Update(); + } + + drawAll(); + cAll->SetTitle( + ((std::string) "Plane " + std::to_string(plane) + (std::string) ", Module " + std::to_string(module)) + .c_str()); + cAll->Update(); +} + +TDialogCanvas* GainMapRefitter::createDialog() { + double width = 300, height = 100; + int nPlanes = gm.GetNumberOfPlanes(); + int nModules = gm.GetModuleIDs(*gm.GetPlaneIDs().begin()).size(); + width = width * nPlanes; + height = height * nModules; + TDialogCanvas* dialog = new TDialogCanvas("Module selection", "", width, height); + dialog->GetCanvasImp()->SetWindowPosition(25, 50); + + int i = 0; + for (auto pm : gm.GetPlaneIDs()) { + int j = 0; + for (auto mm : gm.GetModuleIDs(pm)) { + std::string action = + "gRefitter->changeModule(" + std::to_string(pm) + "," + std::to_string(mm) + ");"; + std::string name = "Plane " + std::to_string(pm) + ", Module " + std::to_string(mm); + TButton* but = + new TButton(name.c_str(), action.c_str(), .1 + i * .8 / nPlanes, .4 + j * .4 / nModules, + .1 + (i + 1) * .8 / nPlanes, .2 + (j + 1) * .4 / nModules); + but->Draw(); + j++; + } + i++; + } + + TButton* butExport = new TButton("Export", "gRefitter->gm.Export();", .4, .1, .6, .3); + butExport->Draw(); + return dialog; +} + +// Global pointer to GainMapRefitter instance. It needs to be global because ROOT's button actions +// are executed in the global scope by CINT/Cling. +GainMapRefitter* gRefitter = nullptr; +// Main function to initialize the class +void REST_RefitGainMap(std::string gainMapFile) { + gRefitter = new GainMapRefitter(); + gRefitter->ImportGainMap(gainMapFile); + + std::cout << std::endl; + std::cout << " ************************************************" << std::endl; + std::cout << " *************** REFITTING MACRO ****************" << std::endl; + std::cout << " * This macro is meant to be used with the FIT *" << std::endl; + std::cout << " * PANEL. First steps: *" << std::endl; + std::cout << " * *" << std::endl; + std::cout << " * 1. Select the module in the dialog canvas. *" << std::endl; + std::cout << " * *" << std::endl; + std::cout << " * 2. Click on the button 'Draw alone' of the *" << std::endl; + std::cout << " * segment you want to refit. The segment *" << std::endl; + std::cout << " * spectrum will be shown in the small canvas *" << std::endl; + std::cout << " * and the Fit Panel will open. It should auto- *" << std::endl; + std::cout << " * matically select the histogram of the corres-*" << std::endl; + std::cout << " * ponding segment. Check that the fit type is *" << std::endl; + std::cout << " * 'Prev. Fit'. The fit function are named *" << std::endl; + std::cout << " * 'g'+peak number (in energy descending order *" << std::endl; + std::cout << " * starting at 0). For example: *" << std::endl; + std::cout << " * - g0 is 22.5keV *" << std::endl; + std::cout << " * - g1 is 8.0keV *" << std::endl; + std::cout << " * *" << std::endl; + std::cout << " * It's important to check option 'Add to list' *" << std::endl; + std::cout << " * in the FIT PANEL (the draw option 'SAME' is *" << std::endl; + std::cout << " * optional but it is helpful). *" << std::endl; + std::cout << " * *" << std::endl; + std::cout << " * 3. After selecting the desired peak function *" << std::endl; + std::cout << " * ('g0', 'g1'...), select the range with the *" << std::endl; + std::cout << " * bottom 'X' named scroller and click on the *" << std::endl; + std::cout << " * button 'Fit'. *" << std::endl; + std::cout << " * *" << std::endl; + std::cout << " * 4. After fitting with the fit panel, click *" << std::endl; + std::cout << " * on the button 'UpdateFits' to update the *" << std::endl; + std::cout << " * fits in the canvas and the calibration curve *" << std::endl; + std::cout << " * of that segment. *" << std::endl; + std::cout << " * *" << std::endl; + std::cout << " * 5. Repeat steps 2-4 for all the segments of *" << std::endl; + std::cout << " * the module and step 1 to change the module. *" << std::endl; + std::cout << " * *" << std::endl; + std::cout << " * 6. Export the gain map with the button *" << std::endl; + std::cout << " * 'Export' on the dialog canvas. *" << std::endl; + std::cout << " ************************************************" << std::endl; + std::cout << std::endl; + + auto dialog = gRefitter->createDialog(); + + Int_t x = 0, y = 0; + UInt_t w = 0, h = 0; + dialog->GetCanvasImp()->GetWindowGeometry(x, y, w, h); + + gRefitter->cAll->SetWindowPosition( + x, + y + 35 + + 100 * + gRefitter->GetGainMap().GetModuleIDs(*gRefitter->GetGainMap().GetPlaneIDs().begin()).size()); + Int_t x2 = 0, y2 = 0; + UInt_t w2 = 0, h2 = 0; + gRefitter->cAll->GetCanvasImp()->GetWindowGeometry(x2, y2, w2, h2); + + gRefitter->cAlone->SetWindowPosition(x2 + w2, y2); +} diff --git a/source/framework/analysis/inc/TRestDataSetGainMap.h b/source/framework/analysis/inc/TRestDataSetGainMap.h index 617c5f0aa..11b2d0d07 100644 --- a/source/framework/analysis/inc/TRestDataSetGainMap.h +++ b/source/framework/analysis/inc/TRestDataSetGainMap.h @@ -150,8 +150,8 @@ class TRestDataSetGainMap : public TRestMetadata { /// Number of bins for the spectrum histograms. Int_t fNBins = 100; //< - /// Cut that defines which events are from this module. - std::string fDefinitionCut = ""; //< + /// Cut that defines which events are from this module. Default is "1" (all events). + std::string fDefinitionCut = "1"; //< /// Number of segments in the x direction. Int_t fNumberOfSegmentsX = 1; //< diff --git a/source/framework/analysis/src/TRestDataSetGainMap.cxx b/source/framework/analysis/src/TRestDataSetGainMap.cxx index eab96f424..499199cdc 100644 --- a/source/framework/analysis/src/TRestDataSetGainMap.cxx +++ b/source/framework/analysis/src/TRestDataSetGainMap.cxx @@ -25,36 +25,34 @@ /// parameters for a given detector with multiple (or just one) modules. /// The modules are defined using the Module class (defined internally). /// It performs a gain correction based on a spatial segmentation of the -/// detector module. This is useful forbig modules such as the ones used -/// in TREX-DM experiment. The input data files for this methods are -/// TRestDataSet for both calculating the calibration parameters and -/// performing the calibration of the desired events. +/// detector module. This is useful for big modules such as the ones used +/// in TREX-DM experiment. The input data for this methods can be a +/// TRestDataSet, a TRestRun file or a file pattern for several TRestRun +/// files for calculating the calibration parameters. The input data file +/// to be used for performing the calibration of the desired events must +/// be a TRestDataSet. /// /// To correct the gain inhomogeneities, the module is divided in /// fNumberOfSegmentsX*fNumberOfSegmentsY segments. The energy peaks provided /// are fitted in each segment. The events are associated to each segment based on /// the observables fSpatialObservableX and fSpatialObservablesY. This results in -/// the following plot that can be obtain with the function -/// TRestDataSetGainMap::Module::DrawSpectrum() +/// the following plot: /// -/// \htmlonly \endhtmlonly -/// ![Peak fitting of each segment. Plot obtain with -/// TRestDataSetGainMap::Module::DrawSpectrum()](drawSpectrum.png) +/// \image html drawSpectrum.png "Plotted by TRestDataSetGainMap::Module::DrawSpectrum()" width=700px /// -/// Also, the peak position provides a gain map that can be plotted with the function -/// TRestDataSetGainMap::Module::DrawGainMap(peakNumber) +/// Also, the peak position provides a gain map: /// -/// \htmlonly \endhtmlonly -/// ![Gain map. Plot obtain with TRestDataSetGainMap::Module::DrawGainMap()](drawGainMap.png) +/// \image html drawGainMap.png "Plotted by TRestDataSetGainMap::Module::DrawGainMap()" width=500px /// /// The result is a better energy resolution with the gain corrected /// calibration (red) than the plain calibration (blue). /// -/// \htmlonly \endhtmlonly -/// ![Gain correction comparison.](gainCorrectionComparison.png) - +/// \image html gainCorrectionComparison.png "Gain correction comparison." width=500px +/// +/// /// ### Parameters -/// * **calibFileName**: name of the file to use for the calibration. It should be a TRestDataSet +/// * **calibFileName**: name of the file to use for the calibration. It can be a TRestDataSet file, +/// a TRestRun file or a file pattern for several TRestRun files. /// * **outputFileName**: name of the file to save this calibration metadata /// * **observable**: name of the observable to be calibrated. It must be a branch of the calibration /// TRestDataSet @@ -100,6 +98,11 @@ /// h->Draw(); /// \endcode /// +/// Example to refit the peaks with the *REST_RefitGainMap* macro (using restRootMacros): +/// \code +/// REST_RefitGainMap("myCalibration.root"); +/// \endcode +/// /// Example to refit manually the peaks of the gain map if any of them is not well fitted /// (using restRoot): /// \code @@ -126,11 +129,10 @@ /// /// History of developments: /// -/// 2023-September: First implementation of TRestDataSetGainMap -/// Álvaro Ezquerro +/// 2023-September: First implementation of TRestDataSetGainMap, Álvaro Ezquerro /// /// \class TRestDataSetGainMap -/// \author: Álvaro Ezquerro aezquerro@unizar.es +/// \author Álvaro Ezquerro aezquerro@unizar.es /// ///
/// @@ -234,6 +236,7 @@ void TRestDataSetGainMap::CalibrateDataSet(const std::string& dataSetFileName, s // Define a new column with the identifier (pmID) of the module for each row (event) std::string pmIDname = (std::string)GetName() + "_pmID"; std::string modCut = fModulesCal[0].GetModuleDefinitionCut(); + if (modCut.empty()) modCut = "1"; // if no cut is defined, use "1" (all events) int pmID = fModulesCal[0].GetPlaneId() * 10 + fModulesCal[0].GetModuleId(); auto columnList = dataFrame.GetColumnNames(); @@ -244,6 +247,7 @@ void TRestDataSetGainMap::CalibrateDataSet(const std::string& dataSetFileName, s for (size_t n = 1; n < fModulesCal.size(); n++) { modCut = fModulesCal[n].GetModuleDefinitionCut(); + if (modCut.empty()) modCut = "1"; // if no cut is defined, use "1" (all events) pmID = fModulesCal[n].GetPlaneId() * 10 + fModulesCal[n].GetModuleId(); dataFrame = dataFrame.Redefine(pmIDname, (modCut + " ? " + std::to_string(pmID) + " : " + pmIDname)); } @@ -286,8 +290,7 @@ void TRestDataSetGainMap::CalibrateDataSet(const std::string& dataSetFileName, s if (outputFileName == dataSetFileName) { // TRestDataSet cannot be overwritten std::string gmName = GetName(); outputFileName = outputFileName.substr(0, outputFileName.find_last_of(".")); // remove extension - outputFileName += "_" + gmName; - outputFileName += TRestTools::GetFileNameExtension(dataSetFileName); + outputFileName += "_" + gmName + "." + TRestTools::GetFileNameExtension(dataSetFileName); } // Export dataset. Exclude columns if requested. @@ -713,15 +716,47 @@ void TRestDataSetGainMap::Module::GenerateGainMap() { return; } - if (!TRestTools::fileExists(dsFileName)) { - RESTError << "Calibration file " << dsFileName << " does not exist." << p->RESTendl; - return; - } - if (!TRestTools::isDataSet(dsFileName)) RESTWarning << dsFileName << " is not a dataset." << p->RESTendl; TRestDataSet dataSet; dataSet.EnableMultiThreading(true); - dataSet.Import(dsFileName); - fDataSetFileName = dsFileName; + if (TRestTools::isDataSet(dsFileName)) { + dataSet.Import(dsFileName); + fDataSetFileName = dsFileName; + } else { + RESTWarning << dsFileName << " is not a dataset. Generating a temporal one..." << p->RESTendl; + // get all the observables needed for the gain map + std::vector obsList; + + obsList.push_back(p->GetObservable()); + obsList.push_back(p->GetSpatialObservableX()); + obsList.push_back(p->GetSpatialObservableY()); + + // look for observables (characterized by having a _ in the name) in the definition cut + auto modDefCutObs = TRestTools::GetObservablesInString(fDefinitionCut, true); + obsList.insert(obsList.end(), modDefCutObs.begin(), modDefCutObs.end()); + + // look for observables in the cut + for (const auto& cut : p->GetCut()->GetCutStrings()) { + auto cutObs = TRestTools::GetObservablesInString(cut, true); + obsList.insert(obsList.end(), cutObs.begin(), cutObs.end()); + } + for (const auto& [variable, condition] : p->GetCut()->GetParamCut()) { + auto cutObs = TRestTools::GetObservablesInString(variable, true); + obsList.insert(obsList.end(), cutObs.begin(), cutObs.end()); + // not sure if any obs can be in the condition. Just in case... + cutObs = TRestTools::GetObservablesInString(condition, true); + obsList.insert(obsList.end(), cutObs.begin(), cutObs.end()); + } + + // remove duplicates + std::sort(obsList.begin(), obsList.end()); + obsList.erase(std::unique(obsList.begin(), obsList.end()), obsList.end()); + + // generate the dataset with the needed observables + dataSet.SetFilePattern(dsFileName); + dataSet.SetObservablesList(obsList); + dataSet.GenerateDataSet(); + fDataSetFileName = dsFileName; + } dataSet.SetDataFrame(dataSet.MakeCut(p->GetCut())); diff --git a/source/framework/tools/inc/TRestTools.h b/source/framework/tools/inc/TRestTools.h index 68f2fa9e3..1c2540896 100644 --- a/source/framework/tools/inc/TRestTools.h +++ b/source/framework/tools/inc/TRestTools.h @@ -83,6 +83,8 @@ class TRestTools { static std::string GetFileNameExtension(const std::string& fullname); static std::string GetFileNameRoot(const std::string& fullname); + static std::vector GetObservablesInString(const std::string& observablesStr, + bool removeDuplicates = true); static int GetBinaryFileColumns(std::string fname); diff --git a/source/framework/tools/src/TRestTools.cxx b/source/framework/tools/src/TRestTools.cxx index 5c868261c..84f37a7db 100644 --- a/source/framework/tools/src/TRestTools.cxx +++ b/source/framework/tools/src/TRestTools.cxx @@ -754,6 +754,8 @@ bool TRestTools::isRunFile(const std::string& filename) { bool TRestTools::isDataSet(const std::string& filename) { if (!isRootFile(filename)) return false; + if (!TRestTools::fileExists(filename)) return false; + TFile* f = TFile::Open((TString)filename); TIter nextkey(f->GetListOfKeys()); @@ -835,6 +837,30 @@ string TRestTools::GetFileNameRoot(const string& fullname) { return filesystem::path(fullname).stem().string(); } +////////////////////////////////////////////////// +/// \brief Returns a vector with the observables names found in the input string. +/// The observables names must contain the character "_" to be identified as such. +/// e.g. Input: "x1_x2 + x3 - x4*y_z". Output: {"x1_x2", "y_z"} +/// Input: "hitsAna_xMean*hitsAna_xMean+hitsAna_yMean*hitsAna_yMean<100" and true, +// Output: {"hitsAna_xMean", "hitsAna_yMean"}. +/// Input: "hitsAna_xMean*hitsAna_xMean+hitsAna_yMean*hitsAna_yMean<100" and false, +/// Output: {"hitsAna_xMean", "hitsAna_xMean", "hitsAna_yMean", "hitsAna_yMean"} +/// +std::vector TRestTools::GetObservablesInString(const std::string& observablesStr, + bool removeDuplicates) { + std::vector obsList; + std::string text = observablesStr; + while (text.find("_") != std::string::npos) { + size_t pos_ = text.find("_"); + size_t beginning = text.find_last_of(" -+*/)(^%<>", pos_) + 1; + size_t end = text.find_first_of(" -+*/)(^%<>", pos_); + std::string obs = text.substr(beginning, end - beginning); + text = Replace(text, obs, "1", 0, removeDuplicates ? 0 : 1); + obsList.push_back(obs); + } + return obsList; +} + /////////////////////////////////////////////// /// \brief Returns the input string but without multiple slashes ("/") ///