Functions
Protodune.cxx File Reference
#include "WireCellSigProc/Microboone.h"
#include "WireCellSigProc/Protodune.h"
#include "WireCellSigProc/Derivations.h"
#include "WireCellUtil/NamedFactory.h"
#include <cmath>
#include <complex>
#include <iostream>
#include <set>

Go to the source code of this file.

Functions

 WIRECELL_FACTORY (pdStickyCodeMitig, WireCell::SigProc::Protodune::StickyCodeMitig, WireCell::IChannelFilter, WireCell::IConfigurable) WIRECELL_FACTORY(pdOneChannelNoise
 
WireCell::IConfigurable WIRECELL_FACTORY(pdRelGainCalib, WireCell::SigProc::Protodune::RelGainCalib, WireCell::IChannelFilter, WireCell::IConfigurable) using namespace WireCell int LedgeIdentify1 (WireCell::Waveform::realseq_t &signal, double baseline, int LedgeStart[3], int LedgeEnd[3])
 
bool LedgeIdentify (WireCell::Waveform::realseq_t &signal, double baseline, int &LedgeStart, int &LedgeEnd)
 

Function Documentation

bool LedgeIdentify ( WireCell::Waveform::realseq_t &  signal,
double  baseline,
int &  LedgeStart,
int &  LedgeEnd 
)

Definition at line 213 of file Protodune.cxx.

213  {
214  int ledge = false;
215  int UNIT = 10;//5; // rebin unit
216  int CONTIN = 10;//20; // length of the continuous region
217  int JUMPS = 2;//4; // how many bins can accidental jump
218  std::vector<int> averaged; // store the rebinned waveform
219  int up = signal.size()/UNIT;// h2->GetNbinsX()/UNIT;
220  int nticks = signal.size();// h2->GetNbinsX();
221  // rebin
222  for(int i=0;i<up;i++){
223  int temp = 0;
224  for(int j=0;j<UNIT;j++){
225  temp += signal.at(i*UNIT+j); //h2->GetBinContent(i*UNIT+1+j);
226  }
227  averaged.push_back(temp);
228  }
229  // refine the selection cuts if there is a large signal
230  auto imax = std::min_element(signal.begin(), signal.end());
231  double max_value = *imax;
232 
233  // if(h2->GetMaximum()-baseline>1000) { CONTIN = 16; JUMPS = 5; }
234  if(max_value-baseline>1000) { CONTIN = 16; JUMPS = 5; }
235  // start judging
236  int decreaseD = 0, tolerence=0;
237  int start = 0;
238  // int end = 0; // never used?
239  for(int i=1;i<up-1;i++){
240  if(averaged.at(i)<averaged.at(i-1)) {
241  if(decreaseD==0) start = i;
242  decreaseD +=1;
243  }
244  else {
245  if(averaged.at(i+1)<averaged.at(i-1)&&tolerence<JUMPS&&decreaseD>0){ // we can ignore several jumps in the decreasing region
246  decreaseD+=2;
247  tolerence++;
248  i = i+1;
249  }
250  else{
251  if(decreaseD>CONTIN){
252  ledge = true;
253  LedgeStart = start*UNIT;
254  break;
255  }
256  else{
257  decreaseD = 0;
258  tolerence=0;
259  start = 0;
260  // end = 0;// end never used?
261  }
262  }
263  }
264  }
265  // find the sharp start edge
266  if(ledge &&LedgeStart>30){
267  // int edge = 0;
268  // int i = LedgeStart/UNIT-1;
269  // if(averaged.at(i)>averaged.at(i-1)&&averaged.at(i-1)>averaged.at(i-2)){ // find a edge
270  // edge = 1;
271  // }
272  // if(edge == 0) ledge = false; // if no edge, this is not ledge
273  // if((averaged.at(i)-averaged.at(i-2)<10*UNIT)&&(averaged.at(i)-averaged.at(i-3)<10*UNIT)) // slope cut
274  // ledge = false;
275  if(averaged.at(LedgeStart/UNIT)-baseline*UNIT>150*UNIT) ledge = false; // ledge is close to the baseline
276  }
277  // test the decay time
278  if(ledge &&LedgeStart>20){
279  double height = 0;
280  if(LedgeStart<5750) { // calculate the height of edge
281  // double tempHeight = h2 ->GetBinContent(LedgeStart+1+200) + h2 ->GetBinContent(LedgeStart+1+220) + h2 ->GetBinContent(LedgeStart+1+180) + h2 ->GetBinContent(LedgeStart+1+240);
282  // height = h2 ->GetBinContent(LedgeStart+1) - tempHeight/4;
283  double tempHeight = signal.at(LedgeStart+200) + signal.at(LedgeStart+220) + signal.at(LedgeStart+180) + signal.at(LedgeStart+240);
284  height = signal.at(LedgeStart) - tempHeight/4;
285  height /= 0.7;
286  }
287  // else height = h2 ->GetBinContent(LedgeStart+1) - h2 ->GetBinContent(6000);
288  else height = signal.at(LedgeStart) - signal.back();
289  if(height<0) height = 80; // norminal value
290  if(height>30&&LedgeStart<5900){ // test the decay with a relatively large height
291  double height50 = 0, height100 = 0;
292  // height50 = h2 ->GetBinContent(LedgeStart+51);
293  // height100 = h2 ->GetBinContent(LedgeStart+101);
294  // double height50Pre = h2 ->GetBinContent(LedgeStart+1)- height*(1-exp(-50/100.)); // minimum 100 ticks decay time
295  // double height100Pre = h2 ->GetBinContent(LedgeStart+1) - height*(1-exp(-100./100)); // minimum 100 ticks decay time
296 
297  height50 = signal.at(LedgeStart+50);
298  height100 = signal.at(LedgeStart+100);
299  double height50Pre = signal.at(LedgeStart)- height*(1-std::exp(-50/100.)); // minimum 100 ticks decay time
300  double height100Pre = signal.at(LedgeStart) - height*(1-std::exp(-100./100)); // minimum 100 ticks decay time
301  // if the measured is much smaller than the predicted, this is not ledge
302  if(height50-height50Pre<-8) ledge = false;
303  if(height100-height100Pre<-8) ledge = false;
304  }
305  }
306 
307  // determine the end of ledge
308  // case 1: find a jump of 10 ADC in the rebinned waveform
309  // case 2: a continuous 20 ticks has an average close to baseline, and a RMS larger than 3 ADC
310  // case 3: reaching the tick 6000 but neither 1 nor 2 occurs
311  if(ledge){
312  LedgeEnd = 0;
313  for(int i = LedgeStart/UNIT; i<up-1; i++){ // case 1
314  if(averaged.at(i+1)-averaged.at(i)>50) {
315  LedgeEnd = i*UNIT+5;
316  break;
317  }
318  }
319  if(LedgeEnd == 0) { // not find a jump, case 2
320  // double tempA[20];
321  WireCell::Waveform::realseq_t tempA(20);
322  for(int i = LedgeStart+80;i<nticks-20;i+=20){
323  for(int j=i;j<20+i;j++){
324  // tempA[j-i] = h2->GetBinContent(j+1);
325  tempA.at(j-i) = signal.at(j);
326  }
327  auto wfinfo = WireCell::Waveform::mean_rms(tempA);
328  // if(TMath::Mean(20,tempA)-baseline<2&&TMath::RMS(20,tempA)>3){
329  if(wfinfo.first - baseline < 2 && wfinfo.second > 3){
330  LedgeEnd = i;
331  break;
332  }
333  }
334  }
335  if(LedgeEnd == 0) LedgeEnd = 6000;
336  }
337  // done, release the memory
338  // vector<int>(averaged).swap(averaged); // is it necessary?
339  return ledge;
340 
341 }
WireCell::Waveform::realseq_t
Sequence< real_t > realseq_t
Definition: Waveform.h:31
nticks
const int nticks
Definition: test_binneddiffusion.cxx:50
height
const double height
Definition: test_raytiling.cxx:34
WireCell::Waveform::mean_rms
std::pair< double, double > mean_rms(const realseq_t &wave)
Definition: Waveform.cxx:24
WireCell::IConfigurable WIRECELL_FACTORY (pdRelGainCalib, WireCell::SigProc::Protodune::RelGainCalib, WireCell::IChannelFilter, WireCell::IConfigurable) using namespace WireCell int LedgeIdentify1 ( WireCell::Waveform::realseq_t &  signal,
double  baseline,
int  LedgeStart[3],
int  LedgeEnd[3] 
)

Definition at line 35 of file Protodune.cxx.

35  {
36  // find a maximum of 3 ledges in one waveform
37  // the number of ledges is returned
38  // the start and end of ledge is stored in the array
39  int UNIT = 10;//5; // rebin unit
40  int CONTIN = 10;//20; // length of the continuous region
41  int JUMPS = 2;//4; // how many bins can accidental jump
42  // We recommend UNIT*CONTIN = 100, UNIT*JUMPS = 20
43  std::vector<int> averaged; // store the rebinned waveform
44  int up = signal.size()/UNIT;
45  int nticks = signal.size();
46  // rebin
47  for(int i=0;i<up;i++){
48  int temp = 0;
49  for(int j=0;j<UNIT;j++){
50  temp += signal.at(i*UNIT+j);
51  }
52  averaged.push_back(temp);
53  }
54  // relax the selection cuts if there is a large signal
55  auto imax = std::min_element(signal.begin(), signal.end());
56  double max_value = *imax;
57  if(max_value-baseline>1000) { CONTIN /=1.25; JUMPS *= 1.2; }
58  // start judging
59  int NumberOfLedge = 0 ; // to be returned
60  int StartOfLastLedgeCandidate = 0; // store the last ledge candidate
61  for(int LE = 0; LE < 3; LE++){ // find three ledges
62  if(LE>0 && StartOfLastLedgeCandidate == 0 ) break; // no ledge candidate in the last round search
63  if(StartOfLastLedgeCandidate>nticks-200) break;// the last candidate is too late in the window
64  if(NumberOfLedge>0 && (LedgeEnd[NumberOfLedge-1]+200)>nticks) break; // the end last ledge reaches the end of readout window
65  int ledge = 0;
66  int decreaseD = 0, tolerence=0;
67  int start = 0;// end = 0; // temporary start and end
68  int StartWindow = 0; // where to start
69  if(StartOfLastLedgeCandidate==0) {
70  StartWindow = 0;
71  }
72  else{
73  if(NumberOfLedge == 0) // no ledge found. start from 200 ticks after the last candidate
74  StartWindow = (StartOfLastLedgeCandidate+200)/UNIT;
75  else{
76  if(StartOfLastLedgeCandidate>LedgeEnd[NumberOfLedge-1]) // the last candidate is not a real ledge
77  // StartWindow = (StartOfLastLedgeCandidate+200)/UNIT;
78  StartWindow = (StartOfLastLedgeCandidate+50)/UNIT;
79  else // the last candidate is a real ledge
80  StartWindow = (LedgeEnd[NumberOfLedge-1]+30)/UNIT;
81  }
82  }
83  for(int i=StartWindow+1;i<up-1;i++){
84  if(averaged.at(i)<averaged.at(i-1)) {
85  if(decreaseD==0) start = i;
86  decreaseD +=1;
87  }
88  else {
89  if(averaged.at(i+1)<averaged.at(i-1)&&tolerence<JUMPS&&decreaseD>0){ // we can ignore several jumps in the decreasing region
90  decreaseD+=2;
91  tolerence++;
92  i = i+1;
93  }
94  else{
95  if(decreaseD>CONTIN){
96  ledge = 1;
97  StartOfLastLedgeCandidate = start*UNIT;
98  break;
99  }
100  else{
101  decreaseD = 0;
102  tolerence=0;
103  start = 0;
104  // end = 0;
105  }
106  }
107  }
108  }
109  // find the sharp start edge
110  // if(ledge == 1&&LedgeStart>30){
111  // int edge = 0;
112  // int i = LedgeStart/UNIT-1;
113  // if(averaged.at(i)>averaged.at(i-1)&&averaged.at(i-1)>averaged.at(i-2)){ // find a edge
114  // edge = 1;
115  // }
116  // if(edge == 0) ledge = 0; // if no edge, this is not ledge
117  // if((averaged.at(i)-averaged.at(i-2)<10*UNIT)&&(averaged.at(i)-averaged.at(i-3)<10*UNIT)) // slope cut
118  // ledge = 0;
119  // if(averaged.at(LedgeStart/UNIT)-baseline*UNIT>300*UNIT) ledge = 0; // ledge is close to the baseline
120  // }
121  // determine the end of ledge
122  // case 1: find a jump of 5 ADC in the rebinned waveform
123  // case 2: a continuous 20 ticks has an average close to baseline, and a RMS larger than 3 ADC
124  // case 3: reaching the tick 6000 but neither 1 nor 2 occurs
125  int tempLedgeEnd = 0;
126  if(ledge ==1){
127  for(int i = StartOfLastLedgeCandidate/UNIT; i<up-1; i++){ // case 1
128  if(averaged.at(i+1)-averaged.at(i)>5*UNIT) {
129  tempLedgeEnd = i*UNIT+5;
130  if(tempLedgeEnd>nticks) tempLedgeEnd = nticks-1;
131  break;
132  }
133  }
134 
135  if(tempLedgeEnd == 0) { // not find a jump, case 2
136  WireCell::Waveform::realseq_t tempA(20);
137  for(int i = StartOfLastLedgeCandidate+80;i<nticks-20;i+=20){
138  for(int j=i;j<20+i;j++){
139  tempA.at(j-i) = signal.at(j);
140  }
141  auto stat = WireCell::Waveform::mean_rms(tempA);
142  if(stat.first-baseline<2&&stat.second>3){
143  tempLedgeEnd = i;
144  break;
145  }
146  }
147  }
148  if(tempLedgeEnd == 0) tempLedgeEnd = nticks-1;
149  }
150  // test the decay time
151  // if(ledge == 1&&StartOfLastLedgeCandidate>20){
152  if(ledge==1){
153  double height = signal.at(StartOfLastLedgeCandidate+1)- signal.at(tempLedgeEnd);
154  if(height<0) ledge = 0; // not a ledge
155  if((tempLedgeEnd-StartOfLastLedgeCandidate) > 80){ // test the decay if the ledge length is longenough
156  double height50 = 0;
157  height50 = signal.at(StartOfLastLedgeCandidate+51);
158  double height50Pre = signal.at(StartOfLastLedgeCandidate+1)- height*(1-exp(-50/100.)); // minimum 100 ticks decay time
159  // if the measured is much smaller than the predicted, this is not ledge
160  if(height50-height50Pre<-12/*-8*/) ledge = 0;
161  // cout << "height50-height50Pre: " << height50-height50Pre << endl;
162  }
163  }
164 
165  // // // find the sharp start edge
166  // if(ledge == 1&&StartOfLastLedgeCandidate>30){
167  // // int edge = 0;
168  // // int i = StartOfLastLedgeCandidate/UNIT-1;
169  // // if(averaged.at(i)>averaged.at(i-1)&&averaged.at(i-1)>averaged.at(i-2)){ // find a edge
170  // // edge = 1;
171  // // }
172  // // if(edge == 0) ledge = 0; // if no edge, this is not ledge
173  // // if((averaged.at(i)-averaged.at(i-2)<10*UNIT)&&(averaged.at(i)-averaged.at(i-3)<10*UNIT)) // slope cut
174  // // ledge = 0;
175  // // if(averaged.at(StartOfLastLedgeCandidate/UNIT)-baseline*UNIT>150*UNIT) ledge = 0; // ledge is close to the baseline
176 
177  // // if(signal.at(tempLedgeEnd) - baseline > 100) ledge=0; // [wgu] ledge end is close to the baseline
178  // if(averaged.at(tempLedgeEnd/UNIT)-baseline*UNIT>5.*UNIT) ledge = 0;
179  // // cout << "averaged.at(StartOfLastLedgeCandidate/UNIT) - baseline*UNIT = " << averaged.at(StartOfLastLedgeCandidate/UNIT)-baseline*UNIT << std::endl;
180  // }
181 
182  if(ledge==1){ // ledge is close to the baseline
183  if(averaged.at(tempLedgeEnd/UNIT)-baseline*UNIT>5.*UNIT) ledge = 0;
184 
185  if(averaged.at(StartOfLastLedgeCandidate/UNIT)-baseline*UNIT>100*UNIT) ledge = 0;
186  }
187 
188  if(ledge==1 && StartOfLastLedgeCandidate>1000){ // ledge always follows a pulse
189  // std::cerr << "[wgu] ch: " << m_ch << " StartOfLastLedgeCandidate: " << StartOfLastLedgeCandidate << std::endl;
190  float hmax=0;
191  for(int i=StartOfLastLedgeCandidate-1000; i<StartOfLastLedgeCandidate-100; i++){
192  if(signal.at(i)>hmax) hmax= signal.at(i);
193  }
194  if(hmax-baseline<200) ledge=0; // no large pre-signal
195  // std::cerr << "[wgu] hmax: " << hmax << std::endl;
196  }
197 
198 
199  if(ledge == 1){
200  LedgeStart[NumberOfLedge] = std::max(StartOfLastLedgeCandidate-20, 0);
201  LedgeEnd[NumberOfLedge] = std::min(tempLedgeEnd+10, (int)signal.size()); // FIXME: ends at a threshold
202  NumberOfLedge++;
203  }
204  } // LE
205 
206  return NumberOfLedge;
207 
208 }
WireCell::Waveform::realseq_t
Sequence< real_t > realseq_t
Definition: Waveform.h:31
python.larbatch_posix.stat
def stat(path)
Definition: larbatch_posix.py:483
nticks
const int nticks
Definition: test_binneddiffusion.cxx:50
height
const double height
Definition: test_raytiling.cxx:34
max
static int max(int a, int b)
Definition: fortranscanner.cpp:60366
cet::sqlite::min
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
WireCell::Waveform::mean_rms
std::pair< double, double > mean_rms(const realseq_t &wave)
Definition: Waveform.cxx:24
WIRECELL_FACTORY ( pdStickyCodeMitig  ,
WireCell::SigProc::Protodune::StickyCodeMitig  ,
WireCell::IChannelFilter  ,
WireCell::IConfigurable   
)
Generated by   doxygen 1.8.11