All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Public Attributes | Private Member Functions | Private Attributes | List of all members
opdet::opHitFinderSBND Class Reference

#include <opHitFinderSBND.hh>

Inheritance diagram for opdet::opHitFinderSBND:

Public Member Functions

 opHitFinderSBND (fhicl::ParameterSet const &p, const detinfo::DetectorClocks *timeService)
 
std::vector< recob::OpHit > MakeHits (const std::vector< raw::OpDetWaveform > &waveforms)
 
 opHitFinderSBND (fhicl::ParameterSet const &p)
 
 opHitFinderSBND (opHitFinderSBND const &)=delete
 
 opHitFinderSBND (opHitFinderSBND &&)=delete
 
opHitFinderSBND & operator= (opHitFinderSBND const &)=delete
 
opHitFinderSBND & operator= (opHitFinderSBND &&)=delete
 
void produce (art::Event &e) override
 

Public Attributes

opdet::sbndPDMapAlg map
 

Private Member Functions

void subtractBaseline (TH1D *hist, std::string pdtype, double &rms)
 
bool findPeak (TH1D *h, size_t &time, double &Area, double rms, double &amplitude, std::string type)
 
void denoise (TH1D *h)
 
void TV1D_denoise (float *input, float *&output, const int width, const float lambda)
 
void subtractBaseline (std::vector< double > &waveform, std::string pdtype, std::string electronicsType, double &rms)
 
bool findAndSuppressPeak (std::vector< double > &waveform, size_t &timebin, double &Area, double &amplitude, const int &threshold, const std::string &opdetType, const std::string &electronicsType)
 
void denoise (std::vector< double > &waveform, std::vector< double > &outwaveform)
 
bool TV1D_denoise (std::vector< double > &waveform, std::vector< double > &outwaveform, const double lambda)
 
void TV1D_denoise_v2 (std::vector< double > &input, std::vector< double > &output, unsigned int width, const double lambda)
 

Private Attributes

std::string fInputModuleName
 
double fSampling
 
double fBaselineSample
 
double fUseDenoising
 
double fPulsePolarityPMT
 
double fPulsePolarityArapuca
 
double fSaturation
 
double fArea1pePMT
 
double fArea1peSiPM
 
int fThresholdPMT
 
int fThresholdArapuca
 
int fEvNumber
 
int fChNumber
 
TH1D * wvfHist
 
std::stringstream histname
 
double fSampling_Daphne
 
bool fUseDenoising
 
std::string opdetType
 
std::string electronicsType
 
int threshold
 
std::vector< double > fwaveform
 
std::vector< double > outwvform
 

Detailed Description

Definition at line 42 of file opHitFinderSBND.hh.

Constructor & Destructor Documentation

opdet::opHitFinderSBND::opHitFinderSBND ( fhicl::ParameterSet const &  p,
const detinfo::DetectorClocks *  timeService 
)

Definition at line 13 of file opHitFinderSBND.cc.

15  {
16  fInputModuleName = p.get< std::string >("InputModule" );
17  fBaselineSample = p.get< int >("BaselineSample"); //in ticks
18  fSaturation = p.get< double >("Saturation" ); //in number of p.e.
19  fArea1pePMT = p.get< double >("Area1pePMT" ); //in ADC*ns for PMTs
20  fArea1peSiPM = p.get< double >("Area1peSiPM" ); //in ADC*ns for SiPMs
21  fThresholdPMT = p.get< double >("ThresholdPMT" ); //in ADC
22  fThresholdArapuca = p.get< double >("ThresholdArapuca"); //in ADC
23  fUseDenoising = p.get< int >("UseDenoising");
24  fPulsePolarityPMT = p.get< int >("PulsePolarityPMT");
25  fPulsePolarityArapuca = p.get< int >("PulsePolarityArapuca");
26 
27  fSampling = (timeService->OpticalClock().Frequency())*500./64; //in GHz. This number is wrong! Therefore the hard coded value
28  }
p
pdgs p
Definition: selectors.fcl:22
opdet::opHitFinderSBND::opHitFinderSBND ( fhicl::ParameterSet const &  p)
explicit

Definition at line 109 of file opHitFinderSBND_module.cc.

110  : EDProducer{p}
111  // Initialize member data here.
112  {
113  fInputModuleName = p.get< std::string >("InputModule" );
114  fBaselineSample = p.get< int >("BaselineSample"); //in ticks
115  fSaturation = p.get< double >("Saturation" ); //in number of p.e.
116  fArea1pePMT = p.get< double >("Area1pePMT" ); //in ADC*ns for PMTs
117  fArea1peSiPM = p.get< double >("Area1peSiPM" ); //in ADC*ns for SiPMs
118  fThresholdPMT = p.get< double >("ThresholdPMT" ); //in ADC
119  fThresholdArapuca = p.get< double>("ThresholdArapuca"); //in ADC
120  fPulsePolarityPMT = p.get< int >("PulsePolarityPMT");
121  fPulsePolarityArapuca = p.get<int>("PulsePolarityArapuca");
122  fUseDenoising = p.get< bool >("UseDenoising");
123 
124  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();
125  fSampling = clockData.OpticalClock().Frequency(); // MHz
126  fSampling_Daphne = p.get<double>("DaphneFrequency");
127 
128  // Call appropriate produces<>() functions here.
129  produces<std::vector<recob::OpHit>>();
130  }
p
pdgs p
Definition: selectors.fcl:22
opdet::opHitFinderSBND::opHitFinderSBND ( opHitFinderSBND const &  )
delete
opdet::opHitFinderSBND::opHitFinderSBND ( opHitFinderSBND &&  )
delete

Member Function Documentation

void opdet::opHitFinderSBND::denoise ( TH1D *  h)
private

Definition at line 139 of file opHitFinderSBND.cc.

139  {
140 
141  int wavelength = (int)h->GetNbinsX();
142  float lambda = 10.0;
143  float* input;
144  float* output;
145 
146  //std::cout <<"denoise wavl:" << wavelength << std::endl;
147 
148  input = (float*)malloc(sizeof(*input)*wavelength);
149  output = (float*)malloc(sizeof(*output)*wavelength);
150 
151  for(int i=0; i<wavelength; i++){
152  input[i] = h->GetBinContent(i);
153  output[i] = input[i];
154  }
155 
156  //std::cout <<"denoise full array" << std::endl;
157 
158  TV1D_denoise(input, output, (const int)wavelength, (const float)lambda);
159 
160  //std::cout <<"denoise denoise" << std::endl;
161 
162  for(int i=0; i<wavelength; i++){
163  if(output[i])h->SetBinContent(i,output[i]);
164  }
165 
166  }
evgen::ldm::lambda
double lambda(double a, double b, double c)
Definition: HNLMakeDecay_tool.cc:134
h
while getopts h
Definition: getStatistics.sh:13
opdet::opHitFinderSBND::TV1D_denoise
void TV1D_denoise(float *input, float *&output, const int width, const float lambda)
Definition: opHitFinderSBND.cc:168
output
BEGIN_PROLOG sequence::SlidingWindowTriggerPatternsOppositeWindows END_PROLOG simSlidingORM6O6 effSlidingORW output
Definition: triggersim_eastmodule_icarus.fcl:134
void opdet::opHitFinderSBND::denoise ( std::vector< double > &  waveform,
std::vector< double > &  outwaveform 
)
private

Definition at line 303 of file opHitFinderSBND_module.cc.

304  {
305 
306  int wavelength = waveform.size();
307  outwaveform = waveform; // copy
308  double lambda = 10.0;
309  const uint retries = 5; uint try_ = 0;
310  if (wavelength > 0) {
311  while (try_ <= retries) {
312  if (TV1D_denoise(waveform, outwaveform, lambda)) break;
313  try_++;
314  mf::LogInfo("opHitFinder") << try_ << "/" << retries
315  << " Coming out of TV1D_denoise() unsuccessfully, "
316  << "using lambda: " << lambda;
317  lambda += 0.1 * lambda;
318  if (try_ == retries) mf::LogWarning("opHitFinder") << "Couldn't denoise!";
319  }
320  }
321  // if (wavelength > 0) TV1D_denoise_v2(waveform, outwaveform, wavelength, lambda);
322 
323  // TODO: fairly certain this for is completely redundant,
324  // and if not a std::move or swap would be better. ~icaza
325  for(int i = 0; i < wavelength; i++) {
326  if(outwaveform[i]) waveform[i] = outwaveform[i];
327  }
328  } // void opHitFinderSBND::denoise()
evgen::ldm::lambda
double lambda(double a, double b, double c)
Definition: HNLMakeDecay_tool.cc:134
opdet::opHitFinderSBND::TV1D_denoise
void TV1D_denoise(float *input, float *&output, const int width, const float lambda)
Definition: opHitFinderSBND.cc:168
bool opdet::opHitFinderSBND::findAndSuppressPeak ( std::vector< double > &  waveform,
size_t &  timebin,
double &  Area,
double &  amplitude,
const int &  threshold,
const std::string &  opdetType,
const std::string &  electronicsType 
)
private

Definition at line 261 of file opHitFinderSBND_module.cc.

266  {
267 
268  std::vector<double>::iterator max_element_it = std::max_element(waveform.begin(), waveform.end());
269  amplitude = *max_element_it;
270  if(amplitude < threshold) return false; // stop if there's no more peaks
271  timebin = std::distance(waveform.begin(), max_element_it);
272 
273  // it_e contains the iterator to the last element in the peak
274  // where waveform is above threshold
275  auto it_e = std::find_if(max_element_it,
276  waveform.end(),
277  [threshold](const double& x)->bool
278  {return x < threshold;} );
279  // it_s contains the iterator to the first element in the peak
280  // where waveform is above threshold
281  auto it_s = std::find_if(std::make_reverse_iterator(max_element_it),
282  std::make_reverse_iterator(waveform.begin()),
283  [threshold](const double& x)->bool
284  {return x < threshold;} ).base();
285 
286  // integrate the area below the peak
287  // note that fSampling is in MHz and
288  // we convert it to GHz here so as to
289  // have an area in ADC*ns.
290  Area = std::accumulate(it_s, it_e, 0.0);
291  if (electronicsType == "daphne"){
292  Area = Area / (fSampling_Daphne / 1000.);}
293  else{
294  Area = Area / (fSampling / 1000.);};
295 
296  // TODO: try to just remove this
297  // TODO: better even, return iterator to last position
298  std::fill(it_s, it_e, 0.0); // zeroes out that peak
299  return true;
300  } // bool opHitFinderSBND::findAndSuppressPeak()
x
process_name opflash particleana ie x
Definition: run_opflash_electron.fcl:90
icarus::distance
double distance(geo::Point_t const &point, CathodeDesc_t const &cathode)
Returns the distance of a point from the cathode.
Definition: CathodeCrossingUtils.cxx:35
trkf::fill
void fill(const art::PtrVector< recob::Hit > &hits, int only_plane)
Definition: KHitContainer.cxx:18
bool opdet::opHitFinderSBND::findPeak ( TH1D *  h,
size_t &  time,
double &  Area,
double  rms,
double &  amplitude,
std::string  type 
)
private

Definition at line 91 of file opHitFinderSBND.cc.

91  {
92 
93  //Gets info from highest peak and suppress it in histogram
94  double aux = h->GetMaximum();
95  double max;
96  size_t time_end, bin, binmax = h->GetMaximumBin();
97  int threshold;
98 
99  if(type=="pmt" || type == "barepmt"){
100  threshold=fThresholdPMT;
101  }else{
102  threshold=fThresholdArapuca;
103  }
104 
105  bin=binmax;
106  amplitude=aux;
107  max=aux;
108 
109  if(aux<threshold)return false;
110 
111  while(aux>=rms){
112  bin++;
113  aux = h->GetBinContent(bin);
114  }
115  time_end=bin-1; //looking for the length of the peak
116 
117  aux=max;
118  while(aux>=rms){
119  bin--;
120  aux = h->GetBinContent(bin);
121  }
122  time = bin+1; //for rise time
123 
124  Area=(h->Integral(time,time_end))/fSampling;
125 
126  bin = time;
127  aux = h->GetBinContent(time);
128 
129  while(aux>=rms){
130  h->SetBinContent(bin,0.0);
131  bin++;
132  aux = h->GetBinContent(bin);
133  }
134  //std::cout << time << " " << time_end << " " << (time_end - time);
135  time=binmax; //returning the peak time
136  return true;
137  }
icarus::ns::util::bin
constexpr details::BinObj< T > bin(T value)
Returns a wrapper to print the specified data in binary format.
h
while getopts h
Definition: getStatistics.sh:13
std::vector< recob::OpHit > opdet::opHitFinderSBND::MakeHits ( const std::vector< raw::OpDetWaveform > &  waveforms)

Definition at line 30 of file opHitFinderSBND.cc.

30  {
31  std::vector< recob::OpHit > ret;
32 
33  size_t timebin=0;
34  double FWHM=1, Area=1, phelec, fasttotal=3./4., rms=0, amplitude=0, time=0;
35  unsigned short frame=1;
36  int histogram_number = 0;
37 
38  for(auto const& wvf : waveforms) {
39  fChNumber = wvf.ChannelNumber();
40  histname.str(std::string());
41  histname << "opchannel_" << fChNumber << "_histo_" << histogram_number;
42  wvfHist = new TH1D(histname.str().c_str(), "Histogram", wvf.size(),0, double(wvf.size()));
43  for(unsigned int i=0;i<wvf.size();i++){
44  wvfHist->SetBinContent(i,wvf[i]);
45  }
46  subtractBaseline(wvfHist, map.pdType(fChNumber), rms);
47  if((map.pdType(fChNumber)=="pmt") || (map.pdType(fChNumber)== "barepmt")){
48  }else{
49  if(fUseDenoising==1) denoise(wvfHist);
50  }
51  int i=1;
52  while(findPeak(wvfHist,timebin,Area,rms,amplitude,map.pdType(fChNumber))){
53  time = wvf.TimeStamp() + (double)timebin/fSampling;
54  if(map.pdType(fChNumber)=="pmt" || map.pdType(fChNumber) == "barepmt"){
55  phelec=Area/fArea1pePMT;
56  // std::cout << 0 << " " << time << " " << Area << " " << phelec << std::endl;
57  }else{
58  phelec=Area/fArea1peSiPM;
59  // std::cout << 1 << " " << time << " " << Area << " " << phelec << std::endl;
60  }
61  i++;
62  recob::OpHit opHit(fChNumber, time, time, frame, FWHM, Area, amplitude, phelec, fasttotal);//including hit info: OpChannel, PeakTime, PeakTimeAbs, Frame, Width, Area, PeakHeight, PE, FastToTotal
63  ret.push_back(opHit);
64  }
65  histogram_number += 1;
66  delete wvfHist;
67  }
68  return ret;
69  }
opdet::opHitFinderSBND::findPeak
bool findPeak(TH1D *h, size_t &time, double &Area, double rms, double &amplitude, std::string type)
Definition: opHitFinderSBND.cc:91
opdet::opHitFinderSBND::map
opdet::sbndPDMapAlg map
Definition: opHitFinderSBND.hh:48
opdet::opHitFinderSBND::subtractBaseline
void subtractBaseline(TH1D *hist, std::string pdtype, double &rms)
Definition: opHitFinderSBND.cc:71
opdet::opHitFinderSBND::histname
std::stringstream histname
Definition: opHitFinderSBND.hh:76
opdet::sbndPDMapAlg::pdType
std::string pdType(size_t ch) const override
Definition: sbndPDMapAlg_tool.cc:37
opHitFinderSBND& opdet::opHitFinderSBND::operator= ( opHitFinderSBND const &  )
delete
opHitFinderSBND& opdet::opHitFinderSBND::operator= ( opHitFinderSBND &&  )
delete
void opdet::opHitFinderSBND::produce ( art::Event &  e)
override

Definition at line 132 of file opHitFinderSBND_module.cc.

133  {
134  // Implementation of required member function here.
135  fEvNumber = e.id().event();
136  mf::LogInfo("opHitFinder") << "Event #" << fEvNumber;
137 
138  std::unique_ptr< std::vector< recob::OpHit > > pulseVecPtr(std::make_unique< std::vector< recob::OpHit > > ());
139  fwaveform.reserve(30000); // TODO: no hardcoded value
140  outwvform.reserve(30000); // TODO: no hardcoded value
141 
142  art::ServiceHandle<art::TFileService> tfs;
143  art::Handle< std::vector< raw::OpDetWaveform > > wvfHandle;
144  std::vector<art::Ptr<raw::OpDetWaveform>> wvfList;
145  if(e.getByLabel(fInputModuleName, wvfHandle))
146  art::fill_ptr_vector(wvfList, wvfHandle);
147 
148  if(!wvfHandle.isValid()) {
149  mf::LogWarning("opHitFinder") << Form("Did not find any waveform");
150  }
151 
152  size_t timebin = 0;
153  double FWHM = 1, Area = 0, phelec, fasttotal = 3./4., rms = 0, amplitude = 0, time = 0;
154  unsigned short frame = 1;
155  //int histogram_number = 0;
156  for(auto const& wvf_P : wvfList) {
157  auto const& wvf = *wvf_P;
158  if (wvf.size() == 0 ) {
159  mf::LogInfo("opHitFinder") << "Empty waveform, continue.";
160  continue;
161  }
162 
163  fChNumber = wvf.ChannelNumber();
164  opdetType = map.pdType(fChNumber);
165  electronicsType = map.electronicsType(fChNumber);
166  if(opdetType == "pmt_coated" || opdetType == "pmt_uncoated") {
167  threshold = fThresholdPMT;
168  }
169  else if((opdetType == "xarapuca_vuv") || (opdetType == "xarapuca_vis")) {
170  threshold = fThresholdArapuca;
171  }
172  else {
173  mf::LogWarning("opHitFinder") << "Unexpected OpChannel: " << opdetType;
174  continue;
175  }
176 
177  fwaveform.resize(wvf.size());
178  for(unsigned int i = 0; i < wvf.size(); i++) {
179  fwaveform[i] = wvf[i];
180  }
181 
182  subtractBaseline(fwaveform, opdetType, electronicsType, rms);
183 
184  if(fUseDenoising) {
185  if((opdetType == "pmt_coated") || (opdetType == "pmt_uncoated")) {
186  }
187  else if((opdetType == "xarapuca_vuv") || (opdetType == "xarapuca_vis")) {
188  denoise(fwaveform, outwvform);
189  }
190  else {
191  mf::LogInfo("opHitFinder") << "Unexpected OpChannel: " << opdetType
192  << ", continue." << std::endl;
193  std::terminate();
194  }
195  }
196 
197  // TODO: pass rms to this function once that's sorted. ~icaza
198  while(findAndSuppressPeak(fwaveform, timebin, Area, amplitude, threshold, opdetType, electronicsType)){
199  if(electronicsType == "daphne") time = wvf.TimeStamp() + (double)timebin / fSampling_Daphne;
200  else time = wvf.TimeStamp() + (double)timebin / fSampling;
201 
202  if(opdetType == "pmt_coated" || opdetType == "pmt_uncoated") {
203  phelec = Area / fArea1pePMT;
204  }
205  else if((opdetType == "xarapuca_vuv") || (opdetType == "xarapuca_vis")) {
206  phelec = Area / fArea1peSiPM;
207  }
208  else {
209  mf::LogWarning("opHitFinder") << "Unexpected OpChannel: " << opdetType
210  << ", continue.";
211  continue;
212  }
213 
214  //including hit info: OpChannel, PeakTime, PeakTimeAbs, Frame, Width, Area, PeakHeight, PE, FastToTotal
215  recob::OpHit opHit(fChNumber, time, time, frame, FWHM, Area, amplitude, phelec, fasttotal);
216  pulseVecPtr->emplace_back(opHit);
217  } // while findAndSuppressPeak()
218  } // for(auto const& wvf : (*wvfHandle)){
219  e.put(std::move(pulseVecPtr));
220  std::vector<double>().swap(fwaveform); // clear and release the memory of fwaveform
221  std::vector<double>().swap(outwvform); // clear and release the memory of outwvform
222  } // void opHitFinderSBND::produce(art::Event & e)
opdet::opHitFinderSBND::map
opdet::sbndPDMapAlg map
Definition: opHitFinderSBND.hh:48
opdet::opHitFinderSBND::outwvform
std::vector< double > outwvform
Definition: opHitFinderSBND_module.cc:92
opdet::opHitFinderSBND::subtractBaseline
void subtractBaseline(TH1D *hist, std::string pdtype, double &rms)
Definition: opHitFinderSBND.cc:71
opdet::opHitFinderSBND::fwaveform
std::vector< double > fwaveform
Definition: opHitFinderSBND_module.cc:91
opdet::sbndPDMapAlg::electronicsType
std::string electronicsType(size_t ch) const
Definition: sbndPDMapAlg_tool.cc:42
e
do i e
Definition: updateVersionList.sh:68
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
opdet::sbndPDMapAlg::pdType
std::string pdType(size_t ch) const override
Definition: sbndPDMapAlg_tool.cc:37
opdet::opHitFinderSBND::findAndSuppressPeak
bool findAndSuppressPeak(std::vector< double > &waveform, size_t &timebin, double &Area, double &amplitude, const int &threshold, const std::string &opdetType, const std::string &electronicsType)
Definition: opHitFinderSBND_module.cc:261
void opdet::opHitFinderSBND::subtractBaseline ( TH1D *  hist,
std::string  pdtype,
double &  rms 
)
private

Definition at line 71 of file opHitFinderSBND.cc.

71  {
72  double baseline = 0.0;
73  rms=0.0;
74  int cnt = 0;
75  for(int i=0; i<fBaselineSample; i++){
76  baseline+=h->GetBinContent(i);
77  rms+=pow((h->GetBinContent(i)),2.0);
78  cnt++;
79  }
80  baseline=baseline/cnt;
81  rms=sqrt(rms/cnt-baseline*baseline);
82  rms=rms/sqrt(cnt-1);
83 
84  if(pdtype=="pmt" || pdtype == "barepmt"){
85  for(int i=0; i<h->GetNbinsX(); i++)h->SetBinContent(i, (fPulsePolarityPMT*(h->GetBinContent(i)-baseline)));
86  }else{
87  for(int i=0; i<h->GetNbinsX(); i++)h->SetBinContent(i, (fPulsePolarityArapuca*(h->GetBinContent(i)-baseline)));
88  }
89  }
h
while getopts h
Definition: getStatistics.sh:13
baseline
BEGIN_PROLOG baseline
Definition: OnlineDecoder.fcl:14
void opdet::opHitFinderSBND::subtractBaseline ( std::vector< double > &  waveform,
std::string  pdtype,
std::string  electronicsType,
double &  rms 
)
private

Definition at line 226 of file opHitFinderSBND_module.cc.

228  {
229  double baseline = 0.0;
230  rms = 0.0;
231  int cnt = 0;
232  double NBins=fBaselineSample;
233  if (electronicsType=="daphne") NBins/=(fSampling/fSampling_Daphne);//correct the number of bins to the sampling frecuency. TODO: use a fixed time interval instead, then use the channel frequency to get the number of bins ~rodrigoa
234  // TODO: this is broken it assumes that the beginning of the
235  // waveform is only noise, which is not always the case. ~icaza.
236  // TODO: use std::accumulate instead of this loop. ~icaza.
237  for(int i = 0; i < NBins; i++) {
238  baseline += waveform[i];
239  rms += std::pow(waveform[i], 2);
240  cnt++;
241  }
242 
243  baseline = baseline / cnt;
244  rms = sqrt(rms / cnt - baseline * baseline);
245  rms = rms / sqrt(cnt - 1);
246 
247  if(pdtype == "pmt_coated" || pdtype == "pmt_uncoated") {
248  for(unsigned int i = 0; i < waveform.size(); i++) waveform[i] = fPulsePolarityPMT * (waveform[i] - baseline);
249  }
250  else if((opdetType == "xarapuca_vuv") || (opdetType == "xarapuca_vis")) {
251  for(unsigned int i = 0; i < waveform.size(); i++) waveform[i] = fPulsePolarityArapuca * (waveform[i] - baseline);
252  }
253  else {
254  mf::LogWarning("opHitFinder") << "Unexpected OpChannel: " << opdetType;
255  return;
256  }
257  }
baseline
BEGIN_PROLOG baseline
Definition: OnlineDecoder.fcl:14
void opdet::opHitFinderSBND::TV1D_denoise ( float *  input,
float *&  output,
const int  width,
const float  lambda 
)
private

Definition at line 168 of file opHitFinderSBND.cc.

168  {
169  if (width>0) { /*to avoid invalid memory access to input[0]*/
170  int k=0, k0=0; /*k: current sample location, k0: beginning of current segment*/
171  float umin=lambda, umax=-lambda; /*u is the dual variable*/
172  float vmin=input[0]-lambda, vmax=input[0]+lambda; /*bounds for the segment's value*/
173  int kplus=0, kminus=0; /*last positions where umax=-lambda, umin=lambda, respectively*/
174  const float twolambda=2.0*lambda; /*auxiliary variable*/
175  const float minlambda=-lambda; /*auxiliary variable*/
176  for (;;) { /*simple loop, the exit test is inside*/
177  while (k==width-1) { /*we use the right boundary condition*/
178  if (umin<0.0) { /*vmin is too high -> negative jump necessary*/
179  do output[k0++]=vmin; while (k0<=kminus);
180  umax=(vmin=input[kminus=k=k0])+(umin=lambda)-vmax;
181  } else if (umax>0.0) { /*vmax is too low -> positive jump necessary*/
182  do output[k0++]=vmax; while (k0<=kplus);
183  umin=(vmax=input[kplus=k=k0])+(umax=minlambda)-vmin;
184  } else {
185  vmin+=umin/(k-k0+1);
186  do output[k0++]=vmin; while(k0<=k);
187  return;
188  }
189  }
190  if ((umin+=input[k+1]-vmin)<minlambda) { /*negative jump necessary*/
191  do output[k0++]=vmin; while (k0<=kminus);
192  vmax=(vmin=input[kplus=kminus=k=k0])+twolambda;
193  umin=lambda; umax=minlambda;
194  } else if ((umax+=input[k+1]-vmax)>lambda) { /*positive jump necessary*/
195  do output[k0++]=vmax; while (k0<=kplus);
196  vmin=(vmax=input[kplus=kminus=k=k0])-twolambda;
197  umin=lambda; umax=minlambda;
198  } else { /*no jump necessary, we continue*/
199  k++;
200  if (umin>=lambda) { /*update of vmin*/
201  vmin+=(umin-lambda)/((kminus=k)-k0+1);
202  umin=lambda;
203  }
204  if (umax<=minlambda) { /*update of vmax*/
205  vmax+=(umax+lambda)/((kplus=k)-k0+1);
206  umax=minlambda;
207  }
208  }
209  }
210  }
211  }
evgen::ldm::lambda
double lambda(double a, double b, double c)
Definition: HNLMakeDecay_tool.cc:134
kplus
process_name eventweight kplus
Definition: fluxreader_evwgh.fcl:61
output
BEGIN_PROLOG sequence::SlidingWindowTriggerPatternsOppositeWindows END_PROLOG simSlidingORM6O6 effSlidingORW output
Definition: triggersim_eastmodule_icarus.fcl:134
k
pdgs k
Definition: selectors.fcl:22
kminus
process_name eventweight kminus
Definition: fluxreader_evwgh.fcl:61
bool opdet::opHitFinderSBND::TV1D_denoise ( std::vector< double > &  waveform,
std::vector< double > &  outwaveform,
const double  lambda 
)
private

Definition at line 331 of file opHitFinderSBND_module.cc.

334  {
335  int width = waveform.size();
336  int k = 0, k0 = 0; // k: current sample location, k0: beginning of current segment
337  double umin = lambda, umax = -lambda; // u is the dual variable
338  double vmin = waveform[0] - lambda, vmax = waveform[0] + lambda; // bounds for the segment's value
339  int kplus = 0, kminus = 0; // last positions where umax=-lambda, umin=lambda, respectively
340  const double twolambda = 2.0 * lambda; // auxiliary variable
341  const double minlambda = -lambda; // auxiliary variable
342  for (;;) { // simple loop, the exit test is inside
343  while (k == width - 1) { // we use the right boundary condition
344  if (umin < 0.0) { // vmin is too high -> negative jump necessary
345  do outwaveform[k0++] = vmin; while (k0 <= kminus);
346  umax = (vmin = waveform[kminus = k = k0]) + (umin = lambda) - vmax;
347  }
348  else if (umax > 0.0) { // vmax is too low -> positive jump necessary
349  do outwaveform[k0++] = vmax; while (k0 <= kplus);
350  umin = (vmax = waveform[kplus = k = k0]) + (umax = minlambda) - vmin;
351  }
352  else {
353  vmin += umin / (k - k0 + 1);
354  do outwaveform[k0++] = vmin; while(k0 <= k);
355  return true;
356  }
357  } // while (k == width - 1)
358  if ((umin += waveform[k + 1] - vmin) < minlambda) { // negative jump necessary
359  if (k0 > width) return false;
360  do outwaveform[k0++] = vmin; while (k0 <= kminus);
361  vmax = (vmin = waveform[kplus = kminus = k = k0]) + twolambda;
362  umin = lambda; umax = minlambda;
363  }
364  else if ((umax += waveform[k + 1] - vmax) > lambda) { // positive jump necessary
365  if (k0 > width) return false;
366  do outwaveform[k0++] = vmax; while (k0 <= kplus);
367  vmin = (vmax = waveform[kplus = kminus = k = k0]) - twolambda;
368  umin = lambda; umax = minlambda;
369  }
370  else { //no jump necessary, we continue
371  k++;
372  if (k > width) return false;
373  if (umin >= lambda) { // update of vmin
374  vmin += (umin - lambda) / ((kminus = k) - k0 + 1);
375  umin = lambda;
376  }
377  if (umax <= minlambda) { // update of vmax
378  vmax += (umax + lambda) / ((kplus = k) - k0 + 1);
379  umax = minlambda;
380  }
381  }
382  } // for (;;)
383  } // bool opHitFinderSBND::TV1D_denoise()
evgen::ldm::lambda
double lambda(double a, double b, double c)
Definition: HNLMakeDecay_tool.cc:134
kplus
process_name eventweight kplus
Definition: fluxreader_evwgh.fcl:61
k
pdgs k
Definition: selectors.fcl:22
kminus
process_name eventweight kminus
Definition: fluxreader_evwgh.fcl:61
void opdet::opHitFinderSBND::TV1D_denoise_v2 ( std::vector< double > &  input,
std::vector< double > &  output,
unsigned int  width,
const double  lambda 
)
private

Definition at line 386 of file opHitFinderSBND_module.cc.

388  {
389  // unsigned int* indstart_low = malloc(sizeof *indstart_low * width);
390  // unsigned int* indstart_up = malloc(sizeof *indstart_up * width);
391  std::vector<unsigned int> indstart_low(width);
392  std::vector<unsigned int> indstart_up(width);
393  unsigned int j_low = 0, j_up = 0, jseg = 0, indjseg = 0, i = 1, indjseg2, ind;
394  double output_low_first = input[0] - lambda;
395  double output_low_curr = output_low_first;
396  double output_up_first = input[0] + lambda;
397  double output_up_curr = output_up_first;
398  const double twolambda = 2.0 * lambda;
399  if (width == 1) {
400  output[0] = input[0];
401  return;
402  }
403  indstart_low[0] = 0;
404  indstart_up[0] = 0;
405  width--;
406  for (; i < width; i++) {
407  if (input[i] >= output_low_curr) {
408  if (input[i] <= output_up_curr) {
409  output_up_curr += (input[i] - output_up_curr) / (i - indstart_up[j_up] + 1);
410  output[indjseg] = output_up_first;
411  while ((j_up > jseg) && (output_up_curr <= output[ind = indstart_up[j_up - 1]]))
412  output_up_curr += (output[ind] - output_up_curr) *
413  ((double)(indstart_up[j_up--] - ind) / (i - ind + 1));
414  if (j_up == jseg) {
415  while ((output_up_curr <= output_low_first) && (jseg < j_low)) {
416  indjseg2 = indstart_low[++jseg];
417  output_up_curr += (output_up_curr - output_low_first) *
418  ((double)(indjseg2 - indjseg) / (i - indjseg2 + 1));
419  while (indjseg < indjseg2) output[indjseg++] = output_low_first;
420  output_low_first = output[indjseg];
421  }
422  output_up_first = output_up_curr;
423  indstart_up[j_up = jseg] = indjseg;
424  }
425  else output[indstart_up[j_up]] = output_up_curr;
426  }
427  else
428  output_up_curr = output[i] = input[indstart_up[++j_up] = i];
429  output_low_curr += (input[i] - output_low_curr) / (i - indstart_low[j_low] + 1);
430  output[indjseg] = output_low_first;
431  while ((j_low > jseg) && (output_low_curr >= output[ind = indstart_low[j_low - 1]]))
432  output_low_curr += (output[ind] - output_low_curr) *
433  ((double)(indstart_low[j_low--] - ind) / (i - ind + 1));
434  if (j_low == jseg) {
435  while ((output_low_curr >= output_up_first) && (jseg < j_up)) {
436  indjseg2 = indstart_up[++jseg];
437  output_low_curr += (output_low_curr - output_up_first) *
438  ((double)(indjseg2 - indjseg) / (i - indjseg2 + 1));
439  while (indjseg < indjseg2) output[indjseg++] = output_up_first;
440  output_up_first = output[indjseg];
441  }
442  if ((indstart_low[j_low = jseg] = indjseg) == i) output_low_first = output_up_first - twolambda;
443  else output_low_first = output_low_curr;
444  }
445  else output[indstart_low[j_low]] = output_low_curr;
446  }
447  else {
448  output_up_curr += ((output_low_curr = output[i] = input[indstart_low[++j_low] = i]) -
449  output_up_curr) / (i - indstart_up[j_up] + 1);
450  output[indjseg] = output_up_first;
451  while ((j_up > jseg) && (output_up_curr <= output[ind = indstart_up[j_up - 1]]))
452  output_up_curr += (output[ind] - output_up_curr) *
453  ((double)(indstart_up[j_up--] - ind) / (i - ind + 1));
454  if (j_up == jseg) {
455  while ((output_up_curr <= output_low_first) && (jseg < j_low)) {
456  indjseg2 = indstart_low[++jseg];
457  output_up_curr += (output_up_curr - output_low_first) *
458  ((double)(indjseg2 - indjseg) / (i - indjseg2 + 1));
459  while (indjseg < indjseg2) output[indjseg++] = output_low_first;
460  output_low_first = output[indjseg];
461  }
462  if ((indstart_up[j_up = jseg] = indjseg) == i) output_up_first = output_low_first + twolambda;
463  else output_up_first = output_up_curr;
464  }
465  else output[indstart_up[j_up]] = output_up_curr;
466  }
467  }
468  /* here i==width (with value the actual width minus one) */
469  if (input[i] + lambda <= output_low_curr) {
470  while (jseg < j_low) {
471  indjseg2 = indstart_low[++jseg];
472  while (indjseg < indjseg2) output[indjseg++] = output_low_first;
473  output_low_first = output[indjseg];
474  }
475  while (indjseg < i) output[indjseg++] = output_low_first;
476  output[indjseg] = input[i] + lambda;
477  }
478  else if (input[i] - lambda >= output_up_curr) {
479  while (jseg < j_up) {
480  indjseg2 = indstart_up[++jseg];
481  while (indjseg < indjseg2) output[indjseg++] = output_up_first;
482  output_up_first = output[indjseg];
483  }
484  while (indjseg < i) output[indjseg++] = output_up_first;
485  output[indjseg] = input[i] - lambda;
486  }
487  else {
488  output_low_curr += (input[i] + lambda - output_low_curr) / (i - indstart_low[j_low] + 1);
489  output[indjseg] = output_low_first;
490  while ((j_low > jseg) && (output_low_curr >= output[ind = indstart_low[j_low - 1]]))
491  output_low_curr += (output[ind] - output_low_curr) *
492  ((double)(indstart_low[j_low--] - ind) / (i - ind + 1));
493  if (j_low == jseg) {
494  if (output_up_first >= output_low_curr)
495  while (indjseg <= i) output[indjseg++] = output_low_curr;
496  else {
497  output_up_curr += (input[i] - lambda - output_up_curr) / (i - indstart_up[j_up] + 1);
498  output[indjseg] = output_up_first;
499  while ((j_up > jseg) && (output_up_curr <= output[ind = indstart_up[j_up - 1]]))
500  output_up_curr += (output[ind] - output_up_curr) *
501  ((double)(indstart_up[j_up--] - ind) / (i - ind + 1));
502  while (jseg < j_up) {
503  indjseg2 = indstart_up[++jseg];
504  while (indjseg < indjseg2) output[indjseg++] = output_up_first;
505  output_up_first = output[indjseg];
506  }
507  indjseg = indstart_up[j_up];
508  while (indjseg <= i) output[indjseg++] = output_up_curr;
509  }
510  }
511  else {
512  while (jseg < j_low) {
513  indjseg2 = indstart_low[++jseg];
514  while (indjseg < indjseg2) output[indjseg++] = output_low_first;
515  output_low_first = output[indjseg];
516  }
517  indjseg = indstart_low[j_low];
518  while (indjseg <= i) output[indjseg++] = output_low_curr;
519  }
520  }
521  // free(indstart_low);
522  // free(indstart_up);
523  }// void opHitFinderSBND::TV1D_denoise_v2()
evgen::ldm::lambda
double lambda(double a, double b, double c)
Definition: HNLMakeDecay_tool.cc:134

Member Data Documentation

std::string opdet::opHitFinderSBND::electronicsType
private

Definition at line 89 of file opHitFinderSBND_module.cc.

double opdet::opHitFinderSBND::fArea1pePMT
private

Definition at line 61 of file opHitFinderSBND.hh.

double opdet::opHitFinderSBND::fArea1peSiPM
private

Definition at line 62 of file opHitFinderSBND.hh.

double opdet::opHitFinderSBND::fBaselineSample
private

Definition at line 56 of file opHitFinderSBND.hh.

int opdet::opHitFinderSBND::fChNumber
private

Definition at line 66 of file opHitFinderSBND.hh.

int opdet::opHitFinderSBND::fEvNumber
private

Definition at line 65 of file opHitFinderSBND.hh.

std::string opdet::opHitFinderSBND::fInputModuleName
private

Definition at line 53 of file opHitFinderSBND.hh.

double opdet::opHitFinderSBND::fPulsePolarityArapuca
private

Definition at line 59 of file opHitFinderSBND.hh.

double opdet::opHitFinderSBND::fPulsePolarityPMT
private

Definition at line 58 of file opHitFinderSBND.hh.

double opdet::opHitFinderSBND::fSampling
private

Definition at line 55 of file opHitFinderSBND.hh.

double opdet::opHitFinderSBND::fSampling_Daphne
private

Definition at line 76 of file opHitFinderSBND_module.cc.

double opdet::opHitFinderSBND::fSaturation
private

Definition at line 60 of file opHitFinderSBND.hh.

int opdet::opHitFinderSBND::fThresholdArapuca
private

Definition at line 64 of file opHitFinderSBND.hh.

int opdet::opHitFinderSBND::fThresholdPMT
private

Definition at line 63 of file opHitFinderSBND.hh.

double opdet::opHitFinderSBND::fUseDenoising
private

Definition at line 57 of file opHitFinderSBND.hh.

bool opdet::opHitFinderSBND::fUseDenoising
private

Definition at line 83 of file opHitFinderSBND_module.cc.

std::vector<double> opdet::opHitFinderSBND::fwaveform
private

Definition at line 91 of file opHitFinderSBND_module.cc.

std::stringstream opdet::opHitFinderSBND::histname
private

Definition at line 76 of file opHitFinderSBND.hh.

opdet::sbndPDMapAlg opdet::opHitFinderSBND::map

Definition at line 48 of file opHitFinderSBND.hh.

std::string opdet::opHitFinderSBND::opdetType
private

Definition at line 88 of file opHitFinderSBND_module.cc.

std::vector<double> opdet::opHitFinderSBND::outwvform
private

Definition at line 92 of file opHitFinderSBND_module.cc.

int opdet::opHitFinderSBND::threshold
private

Definition at line 90 of file opHitFinderSBND_module.cc.

TH1D* opdet::opHitFinderSBND::wvfHist
private

Definition at line 70 of file opHitFinderSBND.hh.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.5