All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SimTestPulseAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: SimTestPulseAna
3 // Module Type: analyzer
4 // File: SimTestPulseAna_module.cc
5 //
6 // Generated at Tue May 23 13:08:14 2017 by Kazuhiro Terao using artmod
7 // from cetpkgsupport v1_11_00.
8 ////////////////////////////////////////////////////////////////////////
9 
10 #include "art/Framework/Core/EDAnalyzer.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 "art_root_io/TFileService.h"
17 #include "fhiclcpp/ParameterSet.h"
18 #include "messagefacility/MessageLogger/MessageLogger.h"
19 #include "larcore/Geometry/Geometry.h"
20 #include "larcorealg/Geometry/GeometryCore.h"
21 #include "lardataobj/RawData/RawDigit.h"
22 #include "lardataobj/RecoBase/Wire.h"
23 #include "lardataobj/RecoBase/Hit.h"
24 
25 #include "ParamHolder.h"
26 #include <TFile.h>
27 #include <TTree.h>
28 
29 class SimTestPulseAna;
30 
31 class SimTestPulseAna : public art::EDAnalyzer
32 {
33 public:
34  explicit SimTestPulseAna(fhicl::ParameterSet const & p);
35  // The destructor generated by the compiler is fine for classes
36  // without bare pointers or other resource use.
37 
38  // Plugins should not be copied or assigned.
39  SimTestPulseAna(SimTestPulseAna const &) = delete;
40  SimTestPulseAna(SimTestPulseAna &&) = delete;
41  SimTestPulseAna & operator = (SimTestPulseAna const &) = delete;
42  SimTestPulseAna & operator = (SimTestPulseAna &&) = delete;
43 
44  // Required functions.
45  void analyze(art::Event const & e) override;
46 
47  void beginJob() override;
48  void endJob() override;
49  void compute_params(const std::vector<float>& wf,
50  float& mean, float& rms, float& min, float& max,
51  int& argmin, int& argmax);
52 
53  TTree* CreateTree(std::string name);
54  TTree* CreateHitTree(std::string name);
55 
56 private:
57  TTree* _hit_tree;
58  TTree* _wire_tree;
59  TTree* _raw_digit_tree;
60  TTree* _filtered_digit_tree;
61  art::InputTag _hit_producer;
62  art::InputTag _wire_producer;
63  art::InputTag _raw_digit_producer;
64  art::InputTag _filtered_digit_producer;
65  size_t _num_samples;
66  size_t _reco_tick_offset;
67  bool _verbose;
68  // TTree variables
69  int _run, _subrun, _event;
70  int _ch, _wire;
71  int _argmax, _argmin;
72  int _tick_offset;
73  int _plane;
74  int _signal_id;
75  int _start_tick;
76  std::vector<float> _wf;
77  float _mean;
78  float _std_dev;
79  float _max;
80  float _min;
81  // Following are specific to hits
82  float _peakTime;
83  float _peakAmplitude;
84  float _rms;
85  float _summedADC;
86  float _integral;
87  float _chisquare;
88  float _baseline;
89 };
90 
91 void SimTestPulseAna::compute_params(const std::vector<float>& wf,
92  float& mean, float& std_dev, float& min, float& max,
93  int& argmin, int& argmax)
94 {
95  mean = std_dev = 0;
96  min = std::numeric_limits<float>::max();
97  max = std::numeric_limits<float>::min();
98  for(size_t idx=0; idx<wf.size(); ++idx)
99  {
100  auto const& v = wf[idx];
101  if(v < min) { min = v; argmin=idx; }
102  if(v > max) { max = v; argmax=idx; }
103  mean += v;
104  }
105  mean /= (float)(wf.size());
106  for(auto const& v : wf)
107  std_dev += pow(v - mean,2);
108  std_dev = sqrt(std_dev / (float)(wf.size()));
109 }
110 
111 
112 SimTestPulseAna::SimTestPulseAna(fhicl::ParameterSet const & p)
113  : EDAnalyzer(p)
114  , _wire_tree(nullptr)
115  , _raw_digit_tree(nullptr)
116  , _filtered_digit_tree(nullptr)
117 {
118  _verbose = p.get<bool> ("Verbose",false);
119  _hit_producer = p.get<art::InputTag>("HitProducer","");
120  _wire_producer = p.get<art::InputTag>("WireProducer","");
121  _raw_digit_producer = p.get<art::InputTag>("RawDigitProducer","");
122  _filtered_digit_producer = p.get<art::InputTag>("FilteredDigitProducer","");
123  _reco_tick_offset = p.get<size_t> ("RecoTickOffset",2400);
124  _num_samples = p.get<size_t> ("NumSample",10);
125 
126  if(_num_samples<1)
127  {
128  std::cerr << "\033[93m[ERROR]\033[00m waveform size will be less than 1 tick ..." << std::endl;
129  throw std::exception();
130  }
131 }
132 
133 TTree* SimTestPulseAna::CreateTree(std::string name)
134 {
135  art::ServiceHandle<art::TFileService> tfs;
136 
137  TTree* tree = tfs->make<TTree>(name.c_str(),"");
138 
139  tree->Branch("run",&_run,"run/I");
140  tree->Branch("subrun",&_subrun,"subrun/I");
141  tree->Branch("event",&_event,"event/I");
142  tree->Branch("ch",&_ch,"ch/I");
143  tree->Branch("wire",&_wire,"wire/I");
144  tree->Branch("plane",&_plane,"plane/I");
145  tree->Branch("signal_id",&_signal_id,"signal_id/I");
146  tree->Branch("start_tick",&_start_tick,"start_tick/I");
147  tree->Branch("tick_offset",&_tick_offset,"tick_offset/I");
148  tree->Branch("wf","std::vector<float>",&_wf);
149  tree->Branch("mean",&_mean,"mean/F");
150  tree->Branch("std",&_std_dev,"std/F");
151  tree->Branch("min",&_min,"min/F");
152  tree->Branch("max",&_max,"max/F");
153  tree->Branch("argmax",&_argmax,"argmax/I");
154  tree->Branch("argmin",&_argmin,"argmin/I");
155  return tree;
156 }
157 
158 TTree* SimTestPulseAna::CreateHitTree(std::string name)
159 {
160  art::ServiceHandle<art::TFileService> tfs;
161 
162  TTree* tree = tfs->make<TTree>(name.c_str(),"");
163 
164  tree->Branch("run",&_run,"run/I");
165  tree->Branch("subrun",&_subrun,"subrun/I");
166  tree->Branch("event",&_event,"event/I");
167  tree->Branch("ch",&_ch,"ch/I");
168  tree->Branch("wire",&_wire,"wire/I");
169  tree->Branch("plane",&_plane,"plane/I");
170  tree->Branch("signal_id",&_signal_id,"signal_id/I");
171  tree->Branch("peakTime",&_peakTime,"peakTime/F");
172  tree->Branch("peakAmplitude",&_peakAmplitude,"peakAmplitude/F");
173  tree->Branch("rms",&_rms,"rms/F");
174  tree->Branch("summedADC",&_summedADC,"summedADC/F");
175  tree->Branch("integral",&_integral,"integral/F");
176  tree->Branch("chisquare",&_chisquare,"chisquare/F");
177  tree->Branch("baseline",&_baseline,"baseline/F");
178 
179  return tree;
180 }
181 
182 void SimTestPulseAna::beginJob()
183 {
184  _hit_tree = ( _hit_producer.empty() ? nullptr : this->CreateHitTree(_hit_producer.label()));
185  _wire_tree = ( _wire_producer.empty() ? nullptr : this->CreateTree(_wire_producer.label()));
186  _raw_digit_tree = ( _raw_digit_producer.empty() ? nullptr : this->CreateTree(_raw_digit_producer.label()) );
187  _filtered_digit_tree = ( _filtered_digit_producer.empty() ? nullptr : this->CreateTree(_filtered_digit_producer.label()+"_f") );
188 }
189 
190 void SimTestPulseAna::endJob()
191 {
192 }
193 
194 void SimTestPulseAna::analyze(art::Event const & e)
195 {
196  art::ServiceHandle<geo::Geometry> geo;
197  _run = e.id().run();
198  _subrun = e.id().subRun();
199  _event = e.id().event();
200 
201  auto const& hit_array = alternative::ParamHolder::get().TruthHitArray();
202 
203  for(auto const& hit : hit_array)
204  {
205  if(hit.tick < _reco_tick_offset) {
206  if(_verbose)
207  std::cout << "[IMBECILE!] Skipping truth charge deposition ID " << hit.signal_id << " @ tick " << hit.tick << std::endl;
208  continue;
209  }
210 
211  _signal_id = hit.signal_id;
212 
213  //
214  // RawDigit by daq
215  //
216  if(!_raw_digit_producer.empty()) {
217 
218  _tick_offset = 0;
219  art::Handle<std::vector<raw::RawDigit> > digit_h;
220  e.getByLabel(_raw_digit_producer,digit_h);
221  if(!digit_h.isValid()) std::cerr << "[IMBECILE!] Failed to fetch RawDigit with label " << _raw_digit_producer << std::endl;
222 
223  for(auto const& ch : hit.channel_list) {
224  _ch = ch;
225  auto const wid = geo->ChannelToWire(ch).front();
226  _wire = (int)(wid.Wire);
227  _plane = wid.Plane;
228  bool found=false;
229  for(auto const& digit : *digit_h) {
230  if((int)(digit.Channel()) != ch) continue;
231 
232  auto const& adcs = digit.ADCs();
233  size_t start = (hit.tick > _num_samples ? hit.tick - _num_samples : 0);
234  size_t end = (adcs.size() > (hit.tick + _num_samples + 1) ? (hit.tick + _num_samples + 1) : adcs.size() -1);
235  _start_tick = start;
236  _wf.resize(end-start+1);
237  for(size_t index=start; index<=end; ++index) _wf[index-start] = adcs[index];
238  this->compute_params(_wf, _mean, _std_dev, _min, _max, _argmin, _argmax);
239 
240  if(_verbose)
241  std::cout << "[IMBECILE!] Storing " << _raw_digit_producer << " raw::RawDigit ... ch=" << _ch
242  << " plane=" << _plane
243  << " wire=" << _wire
244  << " ... " << adcs.size() << " ADC samples start @ tick=" << _start_tick
245  << " ... mean=" << _mean
246  << " std=" << _std_dev
247  << " ... min=" << _min << " @ tick=" << _start_tick + _argmin
248  << " ... max=" << _max << " @ tick=" << _start_tick + _argmax
249  << std::endl;
250 
251  _raw_digit_tree->Fill();
252  found=true;
253  break;
254  }
255  if(!found) std::cout << "[IMBECILE!] Could not find target channel " << _ch
256  << " or wire " << _wire
257  << " for producer " << _raw_digit_producer << std::endl;
258  }
259  }
260 
261  //
262  // RawDigit by noise filter
263  //
264  if(!_filtered_digit_producer.empty()) {
265 
266  _tick_offset = _reco_tick_offset;
267  int signal_tick = hit.tick - _tick_offset;
268 
269  art::Handle<std::vector<raw::RawDigit> > digit_h;
270  e.getByLabel(_filtered_digit_producer,digit_h);
271  if(!digit_h.isValid()) std::cerr << "[IMBECILE!] Failed to fetch RawDigit with label " << _filtered_digit_producer << std::endl;
272 
273  for(auto const& ch : hit.channel_list) {
274  _ch = ch;
275  auto const wid = geo->ChannelToWire(ch).front();
276  _wire = wid.Wire;
277  _plane = wid.Plane;
278  bool found=false;
279  for(auto const& digit : *digit_h) {
280  if((int)(digit.Channel()) != ch) continue;
281 
282  auto const& adcs = digit.ADCs();
283  size_t start = ( signal_tick > (int)(_num_samples) ? signal_tick - _num_samples : 0);
284  size_t end = (adcs.size() > (signal_tick + _num_samples + 1) ? (signal_tick + _num_samples + 1) : adcs.size() -1);
285  _wf.resize(end-start+1);
286  _start_tick = start;
287  for(size_t index=start; index<=end; ++index) _wf[index-start] = adcs[index];
288  this->compute_params(_wf, _mean, _std_dev, _min, _max, _argmin, _argmax);
289 
290  if(_verbose)
291  std::cout << "[IMBECILE!] Storing " << _filtered_digit_producer << " raw::RawDigit ... ch=" << _ch
292  << " plane=" << _plane
293  << " wire=" << _wire
294  << " ... " << adcs.size() << " ADC samples start @ tick=" << _start_tick
295  << " ... mean=" << _mean
296  << " std=" << _std_dev
297  << " ... min=" << _min << " @ tick=" << _start_tick + _argmin
298  << " ... max=" << _max << " @ tick=" << _start_tick + _argmax
299  << std::endl;
300 
301  _filtered_digit_tree->Fill();
302  found=true;
303  break;
304  }
305  if(!found) std::cout << "[IMBECILE!] Could not find target channel " << _ch
306  << " or wire " << _wire
307  << " for producer " << _filtered_digit_producer << std::endl;
308  }
309  }
310 
311  //
312  // Wire
313  //
314  if(!_wire_producer.empty())
315  {
316  _tick_offset = _reco_tick_offset;
317  int signal_tick = hit.tick - _tick_offset;
318 
319  art::Handle<std::vector<recob::Wire> > wire_h;
320  e.getByLabel(_wire_producer,wire_h);
321  if(!wire_h.isValid()) std::cerr << "[IMBECILE!] Failed to fetch Wire with label " << _wire_producer << std::endl;
322 
323  for(auto const& ch : hit.channel_list)
324  {
325  _ch = ch;
326  auto const wid = geo->ChannelToWire(ch).front();
327  _wire = wid.Wire;
328  _plane = wid.Plane;
329  bool found=false;
330  for(auto const& wire : *wire_h)
331  {
332  if((int)(wire.Channel()) != ch) continue;
333 
334  auto const& signalROI = wire.SignalROI();
335 
336  for(const auto& range : signalROI.get_ranges())
337  {
338  // check if this range is relevant
339  if(signal_tick < (int)(range.begin_index()) || signal_tick > (int)(range.begin_index() + range.data().size()))
340  continue;
341 
342  _wf.clear();
343  _wf.reserve(range.data().size());
344  for(auto const& v : range.data()) _wf.push_back(v);
345  _start_tick = range.begin_index();
346  this->compute_params(_wf, _mean, _std_dev, _min, _max, _argmin, _argmax);
347 
348  if(_verbose)
349  std::cout << "[IMBECILE!] Storing " << _wire_producer << " recob::Wire ... ch=" << _ch
350  << " plane=" << _plane
351  << " wire=" << _wire
352  << " ... " << range.data().size() << " ADC samples start @ tick=" << _start_tick
353  << " ... mean=" << _mean
354  << " std=" << _std_dev
355  << " ... min=" << _min << " @ tick=" << _start_tick + _argmin
356  << " ... max=" << _max << " @ tick=" << _start_tick + _argmax
357  << std::endl;
358 
359  _wire_tree->Fill();
360  found=true;
361  break;
362  }
363  break;
364  }
365  if(!found) std::cout << "[IMBECILE!] Could not find target channel " << _ch
366  << " or wire " << _wire
367  << " for producer " << _wire_producer << std::endl;
368  }
369  }
370 
371  // Recover hits
372  if(!_hit_producer.empty())
373  {
374  _tick_offset = _reco_tick_offset;
375  int signal_tick = hit.tick - _tick_offset;
376 
377  art::Handle<std::vector<recob::Hit> > hit_h;
378  e.getByLabel(_hit_producer,hit_h);
379  if(!hit_h.isValid()) std::cerr << "[IMBECILE!] Failed to fetch Wire with label " << _hit_producer << std::endl;
380 
381  for(auto const& ch : hit.channel_list)
382  {
383  _ch = ch;
384  auto const wid = geo->ChannelToWire(ch).front();
385  _wire = wid.Wire;
386  _plane = wid.Plane;
387  bool found=false;
388 
389  float diffPeakTime(1000.);
390 
391  for(auto const& recoHit : *hit_h)
392  {
393  // Look for hit to match channel
394  if((int)(recoHit.Channel()) != ch) continue;
395 
396  // Look for hit range to match tick range
397  if(signal_tick < (int)(recoHit.PeakTime() - 3.*recoHit.RMS()) || signal_tick > (int)(recoHit.PeakTime() + 3.*recoHit.RMS())) continue;
398 
399  if (std::abs(recoHit.PeakTime() - signal_tick) < diffPeakTime)
400  {
401  found = true;
402 
403  _peakTime = recoHit.PeakTime();
404  _peakAmplitude = recoHit.PeakAmplitude();
405  _rms = recoHit.RMS();
406  _summedADC = recoHit.SummedADC();
407  _integral = recoHit.Integral();
408  _chisquare = recoHit.GoodnessOfFit();
409  _baseline = _wf[recoHit.StartTick()-_start_tick];
410 
411  diffPeakTime = std::abs(_peakTime - signal_tick);
412  }
413  }
414 
415  if (found) _hit_tree->Fill();
416 
417  if(!found) std::cout << "[IMBECILE!] Could not find target channel " << _ch
418  << " or wire " << _wire
419  << " for producer " << _hit_producer << std::endl;
420  }
421  }
422  }
423 
424  return;
425 }
426 
427 DEFINE_ART_MODULE(SimTestPulseAna)
SimTestPulseAna::_raw_digit_producer
art::InputTag _raw_digit_producer
Definition: SimTestPulseAna_module.cc:63
SimTestPulseAna::_hit_producer
art::InputTag _hit_producer
Definition: SimTestPulseAna_module.cc:61
alternative::ParamHolder::get
static ParamHolder & get()
Definition: ParamHolder.h:23
cerr
BEGIN_PROLOG could also be cerr
Definition: messageservice.fcl:10
Hit.h
Declaration of signal hit object.
p
pdgs p
Definition: selectors.fcl:22
RawDigit.h
Definition of basic raw digits.
SimTestPulseAna::CreateTree
TTree * CreateTree(std::string name)
Definition: SimTestPulseAna_module.cc:133
hit
process_name hit
Definition: cheaterreco.fcl:51
SimTestPulseAna::SimTestPulseAna
SimTestPulseAna(fhicl::ParameterSet const &p)
Definition: SimTestPulseAna_module.cc:112
GeometryCore.h
Access the description of detector geometry.
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
SimTestPulseAna::compute_params
void compute_params(const std::vector< float > &wf, float &mean, float &rms, float &min, float &max, int &argmin, int &argmax)
Definition: SimTestPulseAna_module.cc:91
SimTestPulseAna::operator=
SimTestPulseAna & operator=(SimTestPulseAna const &)=delete
icarus::ns::util::end
auto end(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:585
SimTestPulseAna::CreateHitTree
TTree * CreateHitTree(std::string name)
Definition: SimTestPulseAna_module.cc:158
alternative::ParamHolder::TruthHitArray
const std::vector< alternative::TruthHit > & TruthHitArray() const
Definition: ParamHolder.cxx:13
SimTestPulseAna::analyze
void analyze(art::Event const &e) override
Definition: SimTestPulseAna_module.cc:194
pmtana::mean
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:13
SimTestPulseAna::_wf
std::vector< float > _wf
Definition: SimTestPulseAna_module.cc:76
e
do i e
Definition: updateVersionList.sh:68
Wire.h
Declaration of basic channel signal object.
name
then echo fcl name
Definition: sbndpoms_genfclwithrunnumber_maker.sh:220
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
SimTestPulseAna::_filtered_digit_producer
art::InputTag _filtered_digit_producer
Definition: SimTestPulseAna_module.cc:64
SimTestPulseAna::_wire_producer
art::InputTag _wire_producer
Definition: SimTestPulseAna_module.cc:62
Geometry.h
art framework interface to geometry description
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
SimTestPulseAna::beginJob
void beginJob() override
Definition: SimTestPulseAna_module.cc:182
Generated by   doxygen 1.8.5