From ac3617f60653dc96a67080827a8afcde201b8696 Mon Sep 17 00:00:00 2001 From: Bryan Gonzalez Date: Mon, 8 Jun 2026 19:22:00 -0500 Subject: [PATCH] Added files to PassN/STM --- STM/inspectSTMFile.fcl | 75 ++++++++++++++++++++++++++ STM/plotDummyDigis.C | 77 ++++++++++++++++++++++++++ STM/plotDummyDigis.C~ | 77 ++++++++++++++++++++++++++ STM/plotPH.c | 88 ++++++++++++++++++++++++++++++ STM/plotZeroSuppression.C | 17 ++++++ STM/testSTMChannel.C | 111 ++++++++++++++++++++++++++++++++++++++ 6 files changed, 445 insertions(+) create mode 100644 STM/inspectSTMFile.fcl create mode 100644 STM/plotDummyDigis.C create mode 100644 STM/plotDummyDigis.C~ create mode 100644 STM/plotPH.c create mode 100644 STM/plotZeroSuppression.C create mode 100644 STM/testSTMChannel.C diff --git a/STM/inspectSTMFile.fcl b/STM/inspectSTMFile.fcl new file mode 100644 index 0000000..986d06b --- /dev/null +++ b/STM/inspectSTMFile.fcl @@ -0,0 +1,75 @@ +# Run StrawAndCaloDigisFromFragments module to convert artdaq::Fragment collections +# into TRK and CAL digi collections. +# Usage: mu2e -c DAQ/test/generateDigiFromFragment.fcl -s -n '-1' +# +# +# Contact person G. Pezzullo +process_name : STMDigiFromFragment + +source : { + module_type : RootInput + fileNames : @nil + maxEvents : -1 +} + + +physics : { + + producers : { + makeSTMDigis : { + module_type : STMDigisFromFragments + stmTag : "daq:ContainerSTM" + verbosityLevel : 1 + saveSTMFragSummary : true + saveRawWithHeaderWaveform_HPGe: false + saveRawWaveform_HPGe: true + saveZSWaveform_HPGe: true + saveRawWithHeaderWaveform_LaBr: false + saveRawWaveform_LaBr: true + saveZSWaveform_LaBr: true + } + + makeSTMBinary : { + module_type : STMBinaryDigisFromFragments + stmTag : "daq:ContainerSTM" + verbosityLevel : 1 + rawFile : "raw.bin" + zsFile : "zs.bin" + phFile : "ph.bin" + rawHeaderFile : "rawWithHeader.bin" + eventFile : "event.bin" + } + } + + + analyzers : { + + stmPrint : { + module_type : STMPrintFragments + stmTag : "daq:ContainerSTM" + } + + } + + t1 : [ makeSTMDigis ] + e1 : [ stmOutput ] + #e1 : [ ] + #e1 : [ frag, stmOutput ] + + trigger_paths : [t1] + end_paths : [e1] + +} + +# print summaries +services.scheduler.wantSummary: true +services.TimeTracker.printSummary: true + +outputs : { + stmOutput : { + module_type : RootOutput + fileName : "dig.stm.art" + outputCommands : [ "keep *_*_*_*", "drop artdaq::*_*_*_*" ] +# outputCommands : [ "drop *_*_*_*", "keep mu2e::STMWaveformDigis_*_*_*" ] + } +} diff --git a/STM/plotDummyDigis.C b/STM/plotDummyDigis.C new file mode 100644 index 0000000..2a30acf --- /dev/null +++ b/STM/plotDummyDigis.C @@ -0,0 +1,77 @@ +// +// +void plotDummyDigis(std::string filename = "stmWaveformDigis.root") { + + TFile* file = new TFile(filename.c_str(), "READ"); + + TCanvas* c1 = new TCanvas(); + c1->Divide(2, 2); + + c1->cd(1); + TH1F* hRawWaveform = (TH1F*) file->Get("plotRawWaveformDigisHPGe/evt0_waveform0"); + TString title(hRawWaveform->GetTitle()); + title += " RAW"; + hRawWaveform->SetTitle(title); + hRawWaveform->SetLineWidth(3); + hRawWaveform->SetLineColor(kGreen); + hRawWaveform->Draw("HIST"); + + c1->cd(2); + TH1F* hZSWaveform = (TH1F*) file->Get("plotZSWaveformDigisHPGe/evt0_waveform0"); + title = hZSWaveform->GetTitle(); + title += " ZS"; + hZSWaveform->SetTitle(title); + hZSWaveform->SetLineWidth(3); + hZSWaveform->SetLineColor(kBlue); + hZSWaveform->Draw("HIST"); + + c1->cd(3); + TH1F* hPHWaveform = (TH1F*) file->Get("plotPHWaveformDigisHPGe/evt0_waveform0"); + title = hPHWaveform->GetTitle(); + title += " PH "; + hPHWaveform->SetTitle(title); + hPHWaveform->SetLineWidth(3); + hPHWaveform->SetLineColor(kRed); + hPHWaveform->Draw("HIST"); + +//New Canvas with all the histograms in the same axis +//No need to redefine the histograms from before + TCanvas* c2 = new TCanvas(); + c2->SetTitle("STMWaveformDigis Super Imposed"); + + hRawWaveform->SetLineColor(kGreen); + hRawWaveform->SetLineWidth(3); + + hZSWaveform->SetLineColor(kBlue); + hZSWaveform->SetLineWidth(3); + + hPHWaveform->SetLineColor(kRed); + hPHWaveform->SetLineWidth(3); + + hRawWaveform->GetYaxis()->SetRangeUser(-1800,1800); + hRawWaveform->Draw(); + hZSWaveform->Draw("SAME"); + hPHWaveform->Draw("SAME"); + + auto legend = new TLegend(); + legend->AddEntry(hRawWaveform, "Raw","1"); + legend->AddEntry(hZSWaveform,"ZS","1"); + legend->AddEntry(hPHWaveform,"PH","1"); + //legend->Draw(); + +//New Canvas where the histogram plots the number of bins from the other two + TCanvas* c3 = new TCanvas(); + c3->SetTitle("Number of bins from hPHWaveform and hZSWaveform"); + + TH1F* hDigiSize = new TH1F("hDigiSize","Number of bins in histograms",1000,0,1000); + hDigiSize->Fill(hPHWaveform->GetNbinsX()); + hDigiSize->Fill(hZSWaveform->GetNbinsX()); + hDigiSize->Draw(); + + + //Below is used to create a copy of the previous canvas + // TCanvas* c4 = new TCanvas(); + // c4->DrawClonePad(); + +} + diff --git a/STM/plotDummyDigis.C~ b/STM/plotDummyDigis.C~ new file mode 100644 index 0000000..2a30acf --- /dev/null +++ b/STM/plotDummyDigis.C~ @@ -0,0 +1,77 @@ +// +// +void plotDummyDigis(std::string filename = "stmWaveformDigis.root") { + + TFile* file = new TFile(filename.c_str(), "READ"); + + TCanvas* c1 = new TCanvas(); + c1->Divide(2, 2); + + c1->cd(1); + TH1F* hRawWaveform = (TH1F*) file->Get("plotRawWaveformDigisHPGe/evt0_waveform0"); + TString title(hRawWaveform->GetTitle()); + title += " RAW"; + hRawWaveform->SetTitle(title); + hRawWaveform->SetLineWidth(3); + hRawWaveform->SetLineColor(kGreen); + hRawWaveform->Draw("HIST"); + + c1->cd(2); + TH1F* hZSWaveform = (TH1F*) file->Get("plotZSWaveformDigisHPGe/evt0_waveform0"); + title = hZSWaveform->GetTitle(); + title += " ZS"; + hZSWaveform->SetTitle(title); + hZSWaveform->SetLineWidth(3); + hZSWaveform->SetLineColor(kBlue); + hZSWaveform->Draw("HIST"); + + c1->cd(3); + TH1F* hPHWaveform = (TH1F*) file->Get("plotPHWaveformDigisHPGe/evt0_waveform0"); + title = hPHWaveform->GetTitle(); + title += " PH "; + hPHWaveform->SetTitle(title); + hPHWaveform->SetLineWidth(3); + hPHWaveform->SetLineColor(kRed); + hPHWaveform->Draw("HIST"); + +//New Canvas with all the histograms in the same axis +//No need to redefine the histograms from before + TCanvas* c2 = new TCanvas(); + c2->SetTitle("STMWaveformDigis Super Imposed"); + + hRawWaveform->SetLineColor(kGreen); + hRawWaveform->SetLineWidth(3); + + hZSWaveform->SetLineColor(kBlue); + hZSWaveform->SetLineWidth(3); + + hPHWaveform->SetLineColor(kRed); + hPHWaveform->SetLineWidth(3); + + hRawWaveform->GetYaxis()->SetRangeUser(-1800,1800); + hRawWaveform->Draw(); + hZSWaveform->Draw("SAME"); + hPHWaveform->Draw("SAME"); + + auto legend = new TLegend(); + legend->AddEntry(hRawWaveform, "Raw","1"); + legend->AddEntry(hZSWaveform,"ZS","1"); + legend->AddEntry(hPHWaveform,"PH","1"); + //legend->Draw(); + +//New Canvas where the histogram plots the number of bins from the other two + TCanvas* c3 = new TCanvas(); + c3->SetTitle("Number of bins from hPHWaveform and hZSWaveform"); + + TH1F* hDigiSize = new TH1F("hDigiSize","Number of bins in histograms",1000,0,1000); + hDigiSize->Fill(hPHWaveform->GetNbinsX()); + hDigiSize->Fill(hZSWaveform->GetNbinsX()); + hDigiSize->Draw(); + + + //Below is used to create a copy of the previous canvas + // TCanvas* c4 = new TCanvas(); + // c4->DrawClonePad(); + +} + diff --git a/STM/plotPH.c b/STM/plotPH.c new file mode 100644 index 0000000..849d44c --- /dev/null +++ b/STM/plotPH.c @@ -0,0 +1,88 @@ +// +// This example ROOT macro plots the steps of the Moving Window Deconvolution algorithm +// This ROOT macro runs on the output of the STMMovingWindowDeconvolution module +// with verbosity level set to 5 +// +void plotPH(std::string filename = "ph.root") { + + TFile* file = new TFile(filename.c_str(), "READ"); + + TLegend* leg = new TLegend(0.6, 0.6, 0.8, 0.8); + leg->SetTextSize(0.04); + leg->SetLineColor(kWhite); + + std::string evt_wvf = "evt1_wvf161"; + + TH1F* hRawWaveform = (TH1F*) file->Get(("phLaBr/h_waveform_"+evt_wvf).c_str()); + hRawWaveform->SetStats(false); + hRawWaveform->SetLineColor(kBlack); + hRawWaveform->Draw("HIST"); + leg->AddEntry(hRawWaveform, "raw waveform", "l"); + + TH1F* hDeconvolved = (TH1F*) file->Get(("phLaBr/h_deconvolved_"+evt_wvf).c_str()); + hDeconvolved->SetLineColor(kOrange); + hDeconvolved->Draw("HIST SAME"); + leg->AddEntry(hDeconvolved, "deconvolution (tau parameter)", "l"); + + TH1F* hDifferentiated = (TH1F*) file->Get(("phLaBr/h_differentiated_"+evt_wvf).c_str()); + hDifferentiated->SetLineColor(kSpring); + hDifferentiated->Draw("HIST SAME"); + leg->AddEntry(hDifferentiated, "differentiation (M parameter)", "l"); + + + TH1F* hAveraged = (TH1F*) file->Get(("phLaBr/h_averaged_"+evt_wvf).c_str()); + hAveraged->SetLineColor(kRed); + hAveraged->SetLineWidth(2); + hAveraged->Draw("HIST SAME"); + leg->AddEntry(hAveraged, "average (L parameter)", "l"); + + TH1F* hBaselineMean = (TH1F*) file->Get(("phLaBr/h_baseline_mean_"+evt_wvf).c_str()); + hBaselineMean->SetLineColor(kBlue); + hBaselineMean->SetLineWidth(2); + hBaselineMean->Draw("HIST SAME"); + leg->AddEntry(hBaselineMean, "baseline mean", "l"); + + TH1F* hPeaks = (TH1F*) file->Get(("phLaBr/h_peaks_"+evt_wvf).c_str()); + hPeaks->SetLineColor(kBlue); + hPeaks->SetLineWidth(2); + hPeaks->Draw("HIST SAME"); + leg->AddEntry(hPeaks, "found peaks", "l"); + + TH1F* hPeakThresh = (TH1F*) file->Get(("phLaBr/h_peak_threshold_"+evt_wvf).c_str()); + hPeakThresh->SetLineColor(kRed); + hPeakThresh->SetLineWidth(2); + hPeakThresh->SetLineStyle(kDashed); + hPeakThresh->Draw("HIST SAME"); + leg->AddEntry(hPeakThresh, "peak threshold", "l"); + + double min_y = hRawWaveform->GetMinimum(); + if (hDeconvolved->GetMinimum() < min_y) { + min_y = hDeconvolved->GetMinimum(); + } + if (hDifferentiated->GetMinimum() < min_y) { + min_y = hDifferentiated->GetMinimum(); + } + if (hAveraged->GetMinimum() < min_y) { + min_y = hAveraged->GetMinimum(); + } + if (hPeaks->GetMinimum() < min_y) { + min_y = hPeaks->GetMinimum(); + } + + double max_y = hRawWaveform->GetMaximum(); + if (hDeconvolved->GetMaximum() > max_y) { + max_y = hDeconvolved->GetMaximum(); + } + if (hDifferentiated->GetMaximum() > max_y) { + max_y = hDifferentiated->GetMaximum(); + } + if (hAveraged->GetMaximum() > max_y) { + max_y = hAveraged->GetMaximum(); + } + if (hPeaks->GetMaximum() > max_y) { + max_y = hPeaks->GetMaximum(); + } + hRawWaveform->GetYaxis()->SetRangeUser(min_y*1.1,max_y*1.1); + + leg->Draw(); +} diff --git a/STM/plotZeroSuppression.C b/STM/plotZeroSuppression.C new file mode 100644 index 0000000..b22aef6 --- /dev/null +++ b/STM/plotZeroSuppression.C @@ -0,0 +1,17 @@ +// +// This example ROOT macro plots an unsuppressed STMWaveformDigi (blue) and a handful of zero-suppressed digis on top (red). +// This ROOT macro runs on the output of Offline/STMReco/fcl/plotSTMWaveformDigis.fcl +// +void plotZeroSuppression(std::string filename = "stmWaveformDigis.root", int n_zp_digis_to_draw = 5) { + + TFile* file = new TFile(filename.c_str(), "READ"); + + TH1F* hRawWaveform = (TH1F*) file->Get("plotHPGeWaveformDigis/evt1_waveform0"); + hRawWaveform->Draw("HIST"); + + for (int i_zp_digi = 0; i_zp_digi < n_zp_digis_to_draw; ++i_zp_digi) { + TH1F* hZPWaveform = (TH1F*) file->Get(Form("plotHPGeWaveformDigisZP/evt1_waveform%d", i_zp_digi)); + hZPWaveform->SetLineColor(kRed); + hZPWaveform->Draw("HIST SAME"); + } +} diff --git a/STM/testSTMChannel.C b/STM/testSTMChannel.C new file mode 100644 index 0000000..53616d6 --- /dev/null +++ b/STM/testSTMChannel.C @@ -0,0 +1,111 @@ +#include "DataProducts/inc/STMChannel.hh" +#include "canvas/Utilities/InputTag.h" + +void testSTMChannel() { + mu2e::STMChannel ch_default; + std::cout << "Default constructor should be \"unknown\": " << ch_default << " : "; + if (ch_default.name().find("unknown") != std::string::npos) { + std::cout << "OK" << std::endl; + } + else { + std::cout << "FAILED" << std::endl; + return -1; + } + std::cout << "and isValid() should be false: isValid() = " << std::boolalpha << ch_default.isValid() << " : "; + if (ch_default.isValid() == false) { + std::cout << "OK" << std::endl; + } + else { + std::cout << "FAILED" << std::endl; + return -1; + } + + std::cout << std::endl; + + // Test findByName() + const int n_channels = 4; + std::string names[n_channels] = {"HPGe", "LaBr", "unknown", "a channel that does not exist"}; + mu2e::STMChannel expectedChannels[n_channels] = {mu2e::STMChannel::HPGe, mu2e::STMChannel::LaBr, mu2e::STMChannel::unknown, mu2e::STMChannel::unknown}; + bool expectedValids[n_channels] = {true, true, false, false}; + for (int i_channel = 0; i_channel < n_channels; ++i_channel) { + std::string name = names[i_channel]; + mu2e::STMChannel expectedChannel = expectedChannels[i_channel]; + bool expectedValid = expectedValids[i_channel]; + + std::cout << "Looking for \"" << name << "\"... " << std::endl; + mu2e::STMChannel ch = mu2e::STMChannel::findByName(name); + std::cout << "Expecting to find " << expectedChannel << " and found: " << ch << " : "; + if (ch == expectedChannel) { + std::cout << "OK" << std::endl; + } + else { + std::cout << "FAILED" << std::endl; + return -1; + } + + std::cout << "isValid() should be " << std::boolalpha << expectedValid << ": isValid() = " << ch.isValid() << " : "; + if (ch.isValid() == expectedValid) { + std::cout << "OK" << std::endl; + } + else { + std::cout << "FAILED" << std::endl; + return -1; + } + + std::cout << std::endl; + } + + std::string expected_all_printed = "List of STM channels: \n 0 HPGe\n 1 LaBr\n 2 unknown\n"; + std::stringstream all_printed; + mu2e::STMChannel::printAll(all_printed); + std::cout << "Output of printAll(): " << std::endl; + std::cout << all_printed.str(); + std::cout << "Is it what we expected? : "; + if (all_printed.str().compare(expected_all_printed) == 0) { + std::cout << "YES" << std::endl; + } + else { + std::cout << "NO" << std::endl; + return -1; + } + std::cout << std::endl; + + art::InputTag hpgeInstance("", "HPGe"); + std::cout << "Testing STMChannel::getChannel() with art::InputTag " << hpgeInstance << ": "; + if (mu2e::STMChannel::getChannel(hpgeInstance) == mu2e::STMChannel::HPGe) { + std::cout << "PASSED" << std::endl; + } + else { + std::cout << "FAILED" << std::endl; + return -1; + } + art::InputTag labrInstance("", "LaBr"); + std::cout << "Testing STMChannel::getChannel() with art::InputTag " << labrInstance << ": "; + if (mu2e::STMChannel::getChannel(labrInstance) == mu2e::STMChannel::LaBr) { + std::cout << "PASSED" << std::endl; + } + else { + std::cout << "FAILED" << std::endl; + return -1; + } + art::InputTag hpgeLabel("labelledHPGe", ""); + std::cout << "Testing STMChannel::getChannel() with art::InputTag " << hpgeLabel << ": "; + if (mu2e::STMChannel::getChannel(hpgeLabel) == mu2e::STMChannel::HPGe) { + std::cout << "PASSED" << std::endl; + } + else { + std::cout << "FAILED" << std::endl; + return -1; + } + art::InputTag labrLabel("labelledLaBr", ""); + std::cout << "Testing STMChannel::getChannel() with art::InputTag " << labrLabel << ": "; + if (mu2e::STMChannel::getChannel(labrLabel) == mu2e::STMChannel::LaBr) { + std::cout << "PASSED" << std::endl; + } + else { + std::cout << "FAILED" << std::endl; + return -1; + } + + std::cout << "All tests passed" << std::endl; +}