23 #include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h" 24 #include "larevt/CalibrationDBI/Interface/ChannelStatusService.h" 31 #include "art_root_io/TFileService.h" 37 #include "TEfficiency.h" 66 bool isSignalInROI(
int starttick,
int endtick,
int maxtick,
int roistart,
int roiend);
108 std::vector<art::Ptr<recob::Wire>> wires;
117 for (
auto const& channel : *simChannelHandle) {
123 if (ch1 % 1000 == 0)
mf::LogInfo(
"EvaluateROIEFF") << ch1;
124 int view = wireReadoutGeom.View(ch1);
125 auto const& timeSlices = channel.TDCIDEMap();
129 int tdctick_previous = -999;
133 std::vector<int> signal_starttick;
134 std::vector<int> signal_endtick;
135 std::vector<double> signal_energydeposits;
136 std::vector<double> signal_energy_max;
137 std::vector<double> signal_max_tdctick;
139 for (
auto const& timeSlice : timeSlices) {
140 auto const tpctime = timeSlice.first;
141 auto const& energyDeposits = timeSlice.second;
142 int tdctick =
static_cast<int>(clockData.TPCTDC2Tick(
double(tpctime)));
143 if (tdctick < 0 || tdctick >
int(detProp.ReadOutWindowSize()) - 1)
continue;
146 double e_deposit = 0;
147 for (
auto const& energyDeposit : energyDeposits) {
148 e_deposit += energyDeposit.energy;
151 if (tdctick_previous == -999) {
152 signal_starttick.push_back(tdctick);
157 else if (tdctick - tdctick_previous != 1) {
158 signal_starttick.push_back(tdctick);
159 signal_endtick.push_back(tdctick_previous);
160 signal_energydeposits.push_back(totalE);
161 signal_energy_max.push_back(maxE);
162 signal_max_tdctick.push_back(maxEtick);
167 else if (tdctick - tdctick_previous == 1) {
169 if (maxE < e_deposit) {
175 tdctick_previous = tdctick;
179 signal_endtick.push_back(tdctick_previous);
180 signal_energydeposits.push_back(totalE);
181 signal_energy_max.push_back(maxE);
182 signal_max_tdctick.push_back(maxEtick);
184 if (signal_starttick.size() == 0 ||
185 (signal_endtick.size() == 1 && signal_endtick.back() == -999))
188 for (
auto& wire : wires) {
189 if (wire->Channel() != ch1)
continue;
193 std::vector<float> vecADC =
202 for (
size_t s = 0; s < signal_starttick.size(); s++) {
204 bool IsSignalOutside =
true;
205 for (
const auto& range : signalROI.
get_ranges()) {
212 signal_max_tdctick[s],
215 IsSignalOutside =
false;
220 if (IsSignalOutside) {
221 h_energy[view]->Fill(signal_energy_max[s]);
227 for (
const auto& range : signalROI.
get_ranges()) {
233 signal_max_tdctick[s],
237 double maxadc_sig = 0;
238 int maxadc_tick = -99;
239 for (
int k = roiFirstBinTick; k < roiLastBinTick; k++) {
240 if (vecADC[k] > maxadc_sig) {
241 maxadc_sig = vecADC[k];
246 double maxE_roi = -999.;
247 double maxE_roi_tick = -999.;
249 if (maxE_roi < signal_energy_max[s]) {
250 maxE_roi = signal_energy_max[s];
251 maxE_roi_tick = signal_max_tdctick[s];
255 for (
size_t s2 = s + 1; s2 < signal_starttick.size(); s2++) {
258 signal_max_tdctick[s2],
261 if (maxE_roi < signal_energy_max[s2]) {
262 maxE_roi = signal_energy_max[s2];
263 maxE_roi_tick = signal_max_tdctick[s2];
265 if (s2 == signal_starttick.size() - 1) { s = s2; }
286 double roi_sig[3] = {0., 0., 0.};
287 double roi_total[3] = {0., 0., 0.};
289 for (
auto& wire : wires) {
291 if (
chStatus.IsBad(wirechannel))
continue;
293 int view = wire->View();
299 std::vector<float> vecADC =
302 roi_total[view] += signalROI.
get_ranges().size();
305 for (
const auto& range : signalROI.
get_ranges()) {
306 bool IsSigROI =
false;
311 double maxadc_sig = -99.;
312 for (
int k = roiFirstBinTick; k < roiLastBinTick; k++) {
318 for (
auto const& channel : (*simChannelHandle)) {
319 if (wirechannel != channel.Channel())
continue;
321 int tdctick_previous = -999;
325 std::vector<int> signal_starttick;
326 std::vector<int> signal_endtick;
327 std::vector<double> signal_energydeposits;
328 std::vector<double> signal_energy_max;
329 std::vector<double> signal_max_tdctick;
331 auto const& timeSlices = channel.TDCIDEMap();
333 for (
auto const& timeSlice : timeSlices) {
334 auto const tpctime = timeSlice.first;
335 auto const& energyDeposits = timeSlice.second;
336 int tdctick =
static_cast<int>(clockData.TPCTDC2Tick(
double(tpctime)));
337 if (tdctick < 0 || tdctick >
int(detProp.ReadOutWindowSize()) - 1)
continue;
339 double e_deposit = 0;
340 for (
auto const& energyDeposit : energyDeposits) {
341 e_deposit += energyDeposit.energy;
344 if (tdctick_previous == -999) {
345 signal_starttick.push_back(tdctick);
350 else if (tdctick - tdctick_previous != 1) {
351 signal_starttick.push_back(tdctick);
352 signal_endtick.push_back(tdctick_previous);
353 signal_energydeposits.push_back(totalE);
354 signal_energy_max.push_back(maxE);
355 signal_max_tdctick.push_back(maxEtick);
360 else if (tdctick - tdctick_previous == 1) {
362 if (maxE < e_deposit) {
368 tdctick_previous = tdctick;
372 signal_endtick.push_back(tdctick_previous);
373 signal_energydeposits.push_back(totalE);
374 signal_energy_max.push_back(maxE);
375 signal_max_tdctick.push_back(maxEtick);
377 if (signal_starttick.size() == 0 ||
378 (signal_endtick.size() == 1 && signal_endtick.back() == -999))
382 for (
size_t s = 0; s < signal_starttick.size(); s++) {
385 signal_max_tdctick[s],
394 h_roi[view]->Fill(0);
398 h_roi[view]->Fill(1);
405 for (
int i = 0; i < 3; i++) {
407 double purity = roi_sig[i] / roi_total[i];
408 if (purity == 1) purity = 1. - 1.e-6;
414 if (roi_total[0] + roi_total[1] + roi_total[2]) {
416 (roi_sig[0] + roi_sig[1] + roi_sig[2]) / (roi_total[0] + roi_total[1] + roi_total[2]);
417 if (purity == 1) purity = 1. - 1.e-6;
427 for (
int i = 0; i < 3; ++i) {
429 tfs->make<TH1D>(Form(
"h_energy_%d", i), Form(
"Plane %d; Energy (MeV);", i), 100, 0, 1);
431 tfs->make<TH1D>(Form(
"h_energy_roi_%d", i), Form(
"Plane %d; Energy (MeV);", i), 100, 0, 1);
433 h_purity[i] = tfs->make<TH1D>(Form(
"h_purity_%d", i), Form(
"Plane %d; Purity;", i), 20, 0, 1);
434 h_roi[i] = tfs->make<TH1D>(Form(
"h_roi_%d", i),
435 Form(
"Plane %d; Purity;", i),
441 tfs->make<TH1D>(Form(
"h1_roi_max_%d", i), Form(
"Plane %d; Max adc;", i), 50, 0, 50);
443 tfs->make<TH1D>(Form(
"h1_roi_max_sim_%d", i), Form(
"Plane %d; Max adc;", i), 50, 0, 50);
446 Form(
"Plane %d; tick diff (maxE - max pulse);", i),
452 h_purity_all = tfs->make<TH1D>(
"h_purity_all",
"All Planes; Purity;", 20, 0, 1);
459 for (
int i = 0; i < 3; ++i) {
462 eff_energy[i]->Write(Form(
"eff_energy_%d", i));
477 return roistart <= maxtick && maxtick < roiend;
art::InputTag fSimulationProducerLabel
Utilities related to art service access.
art::InputTag fWireProducerLabel
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
TEfficiency * eff_energy[3]
constexpr auto abs(T v)
Returns the absolute value of the argument.
TH1D * h1_tickdiff_max[3]
const range_list_t & get_ranges() const
Returns the internal list of non-void ranges.
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
EvaluateROIEff(fhicl::ParameterSet const &p)
EDAnalyzer(fhicl::ParameterSet const &pset)
int TDCtick_t
Type representing a TDC tick.
#define DEFINE_ART_MODULE(klass)
void analyze(art::Event const &e) override
EvaluateROIEff & operator=(EvaluateROIEff const &)=delete
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
bool isSignalInROI(int starttick, int endtick, int maxtick, int roistart, int roiend)
object containing MC truth information necessary for making RawDigits and doing back tracking ...
Declaration of basic channel signal object.
unsigned int ChannelID_t
Type representing the ID of a readout channel.
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)