Skip to content

Commit dddb506

Browse files
lynnt20lynnt20
committed
first commit of working MC software trigger producer
1 parent d6bffea commit dddb506

2 files changed

Lines changed: 184 additions & 0 deletions

File tree

‎sbndcode/Trigger/PMT/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -39,6 +39,8 @@ cet_build_plugin(pmtArtdaqFragmentProducer art::module SOURCE pmtArtdaqFragmentP
3939
cet_build_plugin(pmtSoftwareTriggerProducer art::module SOURCE pmtSoftwareTriggerProducer_module.cc LIBRARIES ${MODULE_LIBRARIES})
4040
cet_build_plugin(pmtTriggerProducer art::module SOURCE pmtTriggerProducer_module.cc LIBRARIES ${MODULE_LIBRARIES})
4141
cet_build_plugin(PMTMetricProducer art::module SOURCE PMTMetricProducer_module.cc LIBRARIES ${MODULE_LIBRARIES})
42+
cet_build_plugin(PMTMCMetricProducer art::module SOURCE PMTMCMetricProducer_module.cc LIBRARIES ${MODULE_LIBRARIES})
43+
4244

4345
install_headers()
4446
install_fhicl()

‎sbndcode/Trigger/PMT/PMTMCMetricProducer_module.cc‎

Lines changed: 182 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,182 @@
1+
////////////////////////////////////////////////////////////////////////
2+
// Class: PMTMCMetricProducer
3+
// Plugin Type: producer (Unknown Unknown)
4+
// File: PMTMCMetricProducer_module.cc
5+
//
6+
// Generated at Mon Nov 10 11:37:29 2025 by Lynn Tung using cetskelgen
7+
// from cetlib version 3.18.02.
8+
////////////////////////////////////////////////////////////////////////
9+
10+
#include "art/Framework/Core/EDProducer.h"
11+
#include "art/Framework/Core/ModuleMacros.h"
12+
#include "art/Framework/Principal/Event.h"
13+
#include "art/Framework/Principal/Handle.h"
14+
#include "art/Framework/Principal/Run.h"
15+
#include "art/Framework/Principal/SubRun.h"
16+
#include "canvas/Utilities/InputTag.h"
17+
#include "fhiclcpp/ParameterSet.h"
18+
#include "messagefacility/MessageLogger/MessageLogger.h"
19+
20+
#include "lardataobj/RawData/OpDetWaveform.h"
21+
#include "larcore/Geometry/Geometry.h"
22+
#include "sbndcode/OpDetSim/sbndPDMapAlg.hh"
23+
#include "sbndaq-artdaq-core/Obj/SBND/pmtSoftwareTrigger.hh"
24+
25+
#include <algorithm>
26+
#include <memory>
27+
28+
#include <TTree.h>
29+
#include <string.h>
30+
#include "TH1D.h"
31+
#include "TFile.h"
32+
33+
namespace sbnd {
34+
namespace trigger {
35+
class PMTMCMetricProducer;
36+
}
37+
}
38+
39+
class sbnd::trigger::PMTMCMetricProducer : public art::EDProducer {
40+
public:
41+
explicit PMTMCMetricProducer(fhicl::ParameterSet const& p);
42+
// The compiler-generated destructor is fine for non-base
43+
// classes without bare pointers or other resource use.
44+
45+
// Plugins should not be copied or assigned.
46+
PMTMCMetricProducer(PMTMCMetricProducer const&) = delete;
47+
PMTMCMetricProducer(PMTMCMetricProducer&&) = delete;
48+
PMTMCMetricProducer& operator=(PMTMCMetricProducer const&) = delete;
49+
PMTMCMetricProducer& operator=(PMTMCMetricProducer&&) = delete;
50+
51+
// Required functions.
52+
void produce(art::Event& e) override;
53+
54+
private:
55+
56+
art::ServiceHandle<art::TFileService> tfs;
57+
std::vector <std::string> fPDTypes;
58+
std::string OpDetWaveformsLabel;
59+
opdet::sbndPDMapAlg fPDSMap;
60+
float fStartTime;
61+
62+
float us_to_ticks = 500.;
63+
float ticks_to_us = 1./us_to_ticks;
64+
65+
std::vector <int> fChannelsToIgnore;
66+
std::vector<uint32_t> sumWvfms(const std::vector<uint32_t>& v1, const std::vector<int16_t>& v2);
67+
float estimateBaseline(std::vector<uint32_t> wvfm);
68+
69+
TTree* tree;
70+
int _run, _subrun, _event;
71+
float _flash_peakpe;
72+
float _flash_peaktime;
73+
int _npmts;
74+
};
75+
76+
77+
sbnd::trigger::PMTMCMetricProducer::PMTMCMetricProducer(fhicl::ParameterSet const& p)
78+
: EDProducer{p} // ,
79+
// More initializers here.
80+
{
81+
fPDTypes = p.get<std::vector<std::string>>("PDTypes",{"pmt_coated", "pmt_uncoated"});
82+
OpDetWaveformsLabel = p.get< std::string >("OpDetWaveformsLabel","opdaq");
83+
fStartTime = p.get<float>("StartTime",-1);
84+
fChannelsToIgnore = p.get<std::vector<int>>("ChannelsToIgnore",{});
85+
86+
art::ServiceHandle<art::TFileService> fs;
87+
tree = fs->make<TTree>("tree","metrictree");
88+
tree->Branch("run",&_run,"run/I");
89+
tree->Branch("subrun",&_subrun,"subrun/I");
90+
tree->Branch("event",&_event,"event/I");
91+
tree->Branch("flash_peakpe",&_flash_peakpe,"flash_peakpe/F");
92+
tree->Branch("flash_peaktime",&_flash_peaktime,"flash_peaktime/F");
93+
tree->Branch("npmts",&_npmts, "npmts/I");
94+
95+
// Call appropriate produces<>() functions here.
96+
produces< std::vector<sbnd::trigger::pmtSoftwareTrigger>>();
97+
// Call appropriate consumes<>() for any products to be retrieved by this module.
98+
}
99+
100+
void sbnd::trigger::PMTMCMetricProducer::produce(art::Event& e)
101+
{
102+
103+
std::unique_ptr<std::vector<sbnd::trigger::pmtSoftwareTrigger>> trig_metrics_v = std::make_unique<std::vector<sbnd::trigger::pmtSoftwareTrigger>>();
104+
sbnd::trigger::pmtSoftwareTrigger trig_metrics;
105+
106+
art::Handle< std::vector< raw::OpDetWaveform > > wvfHandle;
107+
e.getByLabel(OpDetWaveformsLabel, wvfHandle);
108+
if(!wvfHandle.isValid() || wvfHandle->size() == 0){
109+
std::cout << "RawWaveform with label " << OpDetWaveformsLabel << " not found..." << std::endl;
110+
e.put(std::move(trig_metrics_v));
111+
return;
112+
}
113+
114+
// hardcoded to match match data version
115+
size_t flash_len=5000;
116+
std::vector<uint32_t> flash(flash_len, 0);
117+
118+
float readout_start = fStartTime;
119+
float readout_end = readout_start + flash_len*ticks_to_us;
120+
121+
int counter=0;
122+
for(auto const& wvf : (*wvfHandle)) {
123+
int wvf_ch = wvf.ChannelNumber();
124+
auto wvf_ts = wvf.TimeStamp();
125+
auto wvf_end = wvf_ts+wvf.Waveform().size()*2e-3;
126+
127+
if (wvf_end < readout_start || wvf_ts > readout_end ) continue;
128+
if (std::find(fChannelsToIgnore.begin(), fChannelsToIgnore.end(), wvf_ch) != fChannelsToIgnore.end() ) continue;
129+
if (std::find(fPDTypes.begin(), fPDTypes.end(), fPDSMap.pdType(wvf_ch) ) == fPDTypes.end() ) continue;
130+
131+
// find where readout_start is
132+
if (wvf_ts <= readout_start && wvf_end >= readout_end){
133+
int offset = int((readout_start - wvf_ts)*us_to_ticks);
134+
std::vector<short> wvfm_snippet(wvf.Waveform().begin() + offset, wvf.Waveform().begin() + offset + flash_len);
135+
flash = sumWvfms(flash, wvfm_snippet);
136+
counter++;
137+
}
138+
}
139+
140+
auto flash_baseline = estimateBaseline(flash);
141+
auto flash_peak_it = std::min_element(flash.begin(),flash.end());
142+
_flash_peakpe = (flash_baseline-(*flash_peak_it))/12.5;
143+
_flash_peaktime = ((flash_peak_it-flash.begin()))*ticks_to_us + readout_start;
144+
_npmts = counter;
145+
146+
trig_metrics.foundBeamTrigger = true;
147+
trig_metrics.nAboveThreshold = _npmts;
148+
trig_metrics.promptPE = 0;
149+
trig_metrics.prelimPE = 0;
150+
trig_metrics.peakPE = _flash_peakpe;
151+
trig_metrics.peaktime = _flash_peaktime;
152+
trig_metrics_v->push_back(trig_metrics);
153+
154+
e.put(std::move(trig_metrics_v));
155+
156+
tree->Fill();
157+
}
158+
159+
160+
std::vector<uint32_t> sbnd::trigger::PMTMCMetricProducer::sumWvfms(const std::vector<uint32_t>& v1,
161+
const std::vector<int16_t>& v2)
162+
{
163+
size_t result_len = (v1.size() > v2.size()) ? v1.size() : v2.size();
164+
std::vector<uint32_t> result(result_len,0);
165+
for (size_t i = 0; i < result_len; i++){
166+
auto value1 = (i < v1.size()) ? v1[i] : 1e5;
167+
auto value2 = (i < v2.size()) ? v2[i] : 1e5;
168+
result.at(i) = value1 + uint32_t(value2);
169+
}
170+
return result;
171+
}
172+
173+
174+
float sbnd::trigger::PMTMCMetricProducer::estimateBaseline(std::vector<uint32_t> wvfm){
175+
// use a copy of 'wvfm' because nth_element might do something weird!
176+
const auto median_it = wvfm.begin() + wvfm.size() / 2;
177+
std::nth_element(wvfm.begin(), median_it , wvfm.end());
178+
auto median = *median_it;
179+
return median;
180+
}
181+
182+
DEFINE_ART_MODULE(sbnd::trigger::PMTMCMetricProducer)

0 commit comments

Comments
 (0)