Skip to content

Commit 593fbbd

Browse files
crogancrogan
committed
added pulse processing (a la DRS merge exec) to scopehdf2root -- now doing same operations as Caltech code
1 parent f80f095 commit 593fbbd

File tree

7 files changed

+223
-10
lines changed

7 files changed

+223
-10
lines changed

‎DRS/Makefile‎

Lines changed: 1 addition & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -51,15 +51,12 @@ $(TARGET4) : $(OBJ4)
5151
@echo $@
5252
$(LD) $(CPPFLAGS) -o $(TARGET4) $(OBJ4) $(LDFLAGS)
5353

54-
scopehdf2root: src/NetScope/scopehdf2root.C src/NetScope/hdf5io.o
54+
scopehdf2root: src/NetScope/scopehdf2root.C src/NetScope/hdf5io.o src/Aux.o
5555
$(LD) $(CPPFLAGS) $^ $(LIBS) $(LDFLAGS) -o $@
5656

5757
convert_ots2hdf5: src/NetScope/convert_ots2hdf5.C src/NetScope/hdf5io.o src/NetScope/fifo.o
5858
$(LD) $(CPPFLAGS) $^ $(LIBS) $(LDFLAGS) -o $@
5959

60-
#src/NetScope/hdf5io.o: src/NetScope/hdf5io.cc
61-
# $(CC) -Wall -O2 -DH5_NO_DEPRECATED_SYMBOLS -I$(INC)/include/NetScope/ -o $@ -c $<
62-
6360
%.o : %.cc
6461
@echo $@
6562
$(CXX) $(CPPFLAGS) -o $@ -c $<

‎DRS/dat2rootPixels.cc‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -442,7 +442,8 @@ int main(int argc, char **argv) {
442442

443443
// Check if all channels were active (if 8 channels active return 3072)
444444
int nsample = (event_header & 0xfff) / 3;
445-
//cout << nsample << " nsample" << endl;
445+
446+
// hard-coded to protect against corruption
446447
nsample = 1024;
447448

448449
// Define time coordinate

‎DRS/include/Aux.hh‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,13 +20,15 @@ static float amplificationFactor[nPoints] = { 27.9/1.37*10 , 26.8/1.37*10, 31.6/
2020
double GetAmplificationFactor ( double measuredAmplitude );
2121
TGraphErrors* WeierstrassTransform( short* channel, float* time, TString pulseName, double sigma = 1.0, bool makePlot = false );
2222
TGraphErrors* GetTGraph( double* channel, float* time );
23+
TGraphErrors* GetTGraph( float* channel, float* time );
2324
TGraphErrors GetTGraph( short* channel, float* time );
2425
double GetGaussTime( TGraphErrors* pulse );
2526
void HighPassFilter( short* channel, double* filteredCurrent, float* time, double R = -1.0, double C = -1.0 );
2627
void NotchFilter( short* channel, double* filteredCurrent, float* time, double R = -1.0, double C = -1.0, double L = -1.0 );
2728
int FindMin( int n, short *a);
2829
int FindRealMin( int n, short *a);
2930
int FindMinAbsolute( int n, double *a);
31+
int FindMinAbsolute( int n, float *a);
3032
int FindMinAbsolute( int n, short *a);
3133
int FindMinFirstPeakAboveNoise( int n, short *a);
3234
float GausFit_MeanTime(TGraphErrors * pulse, const float index_first, const float index_last);
@@ -41,9 +43,12 @@ float GausFit_MeanTime(TGraphErrors* pulse, const float index_first, const float
4143
float GetBaseline( int peak, short *a );
4244
float GetBaseline(TGraphErrors * pulse, int i_low, int i_high, TString fname );
4345
float GetPulseIntegral(int peak, short *a, std::string option = "");
46+
float GetPulseIntegral(int peak, float *a, std::string option = "");
4447
float GetPulseIntegral(int peak, int nsamples, short *a, float *t);
48+
float GetPulseIntegral(int peak, int nsamples, float *a, float *t);
4549
float ConstantThresholdTime(TGraphErrors * pulse, const float threshold);
4650
bool isRinging( int peak, short *a );
51+
bool isRinging( int peak, float *a );
4752

4853

4954
#endif

‎DRS/processNetScope_ETL_122017.sh‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,12 @@ NC='\033[0m' # No Color
1010
runNum=$1
1111
echo "Processing NetScope data for run ${runNum}"
1212

13+
# ./envHDF5.sh
14+
export HDF5_ROOT_PATH=`pwd`/hdf5
15+
export HDF5_INCLUDE_PATH=${HDF5_ROOT_PATH}/include
16+
export HDF5_LIBRARY_PATH=${HDF5_ROOT_PATH}/lib
17+
export LD_LIBRARY_PATH=${HDF5_LIBRARY_PATH}:$LD_LIBRARY_PATH
18+
1319
echo ""
1420
echo "Converting raw NetScope data to HDF5 format"
1521
echo -e "./convert_ots2hdf5 /eos/uscms/store/user/cmstestbeam/ETL/MT6Section1Data/122017/OTSDAQ/NetScopeRaw/RawDataSaver0NetScope_Run${runNum}_0_Raw.dat /eos/uscms/store/user/cmstestbeam/ETL/MT6Section1Data/122017/OTSDAQ/NetScopeHDF5/Run${runNum}.hdf5"

‎DRS/src/Aux.cc‎

Lines changed: 78 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -96,6 +96,21 @@ void NotchFilter( short* channel, double* filteredCurrent, float* time, double R
9696
return;
9797
};
9898

99+
TGraphErrors* GetTGraph( float* channel, float* time )
100+
{
101+
//Setting Errors
102+
float errorX[1024], errorY[1024], channelFloat[1024];
103+
float _errorY = 0.00; //5%error on Y
104+
for ( int i = 0; i < 1024; i++ )
105+
{
106+
errorX[i] = .0;
107+
errorY[i] = _errorY*channel[i];
108+
channelFloat[i] = -channel[i];
109+
}
110+
TGraphErrors* tg = new TGraphErrors( 1024, time, channelFloat, errorX, errorY );
111+
return tg;
112+
};
113+
99114
TGraphErrors* GetTGraph( double* channel, float* time )
100115
{
101116
//Setting Errors
@@ -151,6 +166,21 @@ int FindMinAbsolute( int n, short *a) {
151166
return loc;
152167
}
153168

169+
int FindMinAbsolute( int n, float *a) {
170+
171+
if (n <= 0 || !a) return -1;
172+
float xmin = a[5];
173+
int loc = 0;
174+
for (int i = 5; i < n-10; i++) {
175+
if ( a[i] < xmin && a[i+1] < 0.5*a[i] && a[i] < -40. )
176+
{
177+
xmin = a[i];
178+
loc = i;
179+
}
180+
}
181+
182+
return loc;
183+
}
154184

155185
int FindMinAbsolute( int n, double *a) {
156186

@@ -707,6 +737,24 @@ float GetPulseIntegral(int peak, short *a, std::string option)
707737

708738
}
709739

740+
float GetPulseIntegral(int peak, float *a, std::string option)
741+
{
742+
float integral = 0.;
743+
744+
if (option == "full") {
745+
for (int i=5; i < 1020; i++) {
746+
integral += a[i] * 0.2 * 1e-9 * (1.0/4096.0) * (1.0/50.0) * 1e12; //in units of pC, for 50Ohm termination
747+
}
748+
}
749+
else {
750+
for (int i=peak-20; i < peak+25; i++) {
751+
integral += a[i] * 0.2 * 1e-9 * (1.0/4096.0) * (1.0/50.0) * 1e12; //in units of pC, for 50Ohm termination
752+
}
753+
}
754+
return -1.0 * integral;
755+
756+
}
757+
710758
float GetPulseIntegral(int peak, int nsamples, short *a, float *t) //returns charge in pC asssuming 50 Ohm termination
711759
{
712760
float integral = 0.;
@@ -720,6 +768,19 @@ float GetPulseIntegral(int peak, int nsamples, short *a, float *t) //returns cha
720768

721769
}
722770

771+
float GetPulseIntegral(int peak, int nsamples, float *a, float *t) //returns charge in pC asssuming 50 Ohm termination
772+
{
773+
float integral = 0.;
774+
for (int i=peak-nsamples; i < peak+nsamples; i++) {
775+
//Simpson's Rule for equaled space-->Cartwright correction for unequaled space, only worked for odd points
776+
integral += ( (t[i+2]-t[i]) / 6.0 ) * ( ( 2-(t[i+2]-t[i+1])/(t[i+1]-t[i]) )* a[i] + (t[i+2]-t[i])*(t[i+2]-t[i])/((t[i+2]-t[i+1])*(t[i+1]-t[i])) * a[i+1] + ( 2-(t[i+1]-t[i])/(t[i+2]-t[i+1]) ) * a[i+2] ) * 1e-9 * (1.0/4096.0) * (1.0/50.0) * 1e12; //in units of pC, for 50Ohm termination
777+
i++;
778+
}
779+
780+
return -1.0 * integral;
781+
782+
}
783+
723784
//----------------------------------------------
724785
//Gaussian Filter to reduce high frequency noise
725786
//----------------------------------------------
@@ -813,3 +874,20 @@ bool isRinging( int peak, short *a )
813874
if ( right_max > 0.5*fabs(a[peak]) ) return true;
814875
return false;
815876
};
877+
878+
bool isRinging( int peak, float *a )
879+
{
880+
short left_max, right_max;//Left and right maximum amplitude from the minimum
881+
left_max = a[peak];
882+
right_max = a[peak];
883+
884+
for ( int i = 1; i < 150; i++)
885+
{
886+
if ( a[peak-i] > left_max ) left_max = a[peak-i];
887+
if ( a[peak+i] > right_max ) right_max = a[peak+i];
888+
}
889+
890+
if ( left_max > 0.5*fabs(a[peak]) ) return true;
891+
if ( right_max > 0.5*fabs(a[peak]) ) return true;
892+
return false;
893+
};

‎DRS/src/NetScope/hdf5io.o‎

9.46 KB
Binary file not shown.

‎DRS/src/NetScope/scopehdf2root.C‎

Lines changed: 131 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,10 @@
1919
#include <TGraphErrors.h>
2020
#include <TCanvas.h>
2121

22+
//LOCAL INCLUDES
23+
#include "Aux.hh"
24+
25+
// NetScope INCLUDES
2226
#include "common.h"
2327
#include "hdf5io.h"
2428

@@ -44,6 +48,13 @@ std::string ParseCommandLine( int argc, char* argv[], std::string opt )
4448
int main(int argc, char **argv) {
4549
size_t i, j, iCh, iEvent=0, nEvents=0, frameSize, nEventsInFile;
4650
char *inFileName;
51+
52+
// hardcoded number of time samples for pulses
53+
int NSAMPLES = 1000;
54+
55+
// need to do a reverse ADC to use same
56+
// code as DRS, then reverse
57+
float ADC = (1.0 / 4096.0);
4758

4859
struct hdf5io_waveform_file *waveformFile;
4960
struct waveform_attribute waveformAttr;
@@ -106,8 +117,8 @@ int main(int argc, char **argv) {
106117

107118
// initialization of variables for TTree output
108119
int event;
109-
float time[2][1000]; // calibrated time
110-
float channel[4][1000]; // input (in V)
120+
float time[2][NSAMPLES]; // calibrated time
121+
float channel[4][NSAMPLES]; // input (in V)
111122
float xmin[4]; // location of peak
112123
float base[4]; // baseline voltage
113124
float amp[4]; // pulse amplitude
@@ -126,6 +137,24 @@ int main(int argc, char **argv) {
126137
float constantThresholdTime[4];
127138
bool _isRinging[4];
128139

140+
for(iCh=0; iCh<4; iCh++) {
141+
xmin[iCh] = 0.;
142+
amp [iCh] = 0.;
143+
base[iCh] = 0.;
144+
integral[iCh] = 0.;
145+
integralFull[iCh] = 0.;
146+
gauspeak[iCh] = 0.;
147+
sigmoidTime[iCh] = 0.;
148+
linearTime0[iCh] = 0.;
149+
linearTime15[iCh] = 0.;
150+
linearTime30[iCh] = 0.;
151+
linearTime45[iCh] = 0.;
152+
linearTime60[iCh] = 0.;
153+
risetime[iCh] = 0.;
154+
constantThresholdTime[iCh] = 0.;
155+
_isRinging[iCh] = false;
156+
}
157+
129158
tree->Branch("event", &event, "event/I");
130159
tree->Branch("channel", channel, "channel[4][1000]/F");
131160
tree->Branch("time", time, "time[2][1000]/F");
@@ -162,7 +191,7 @@ int main(int argc, char **argv) {
162191

163192
event = waveformEvent.eventId;
164193

165-
for(i = 0; i < 1000; i++) {
194+
for(i = 0; i < NSAMPLES; i++) {
166195
if(i < waveformFile->nPt)
167196
time[0][i] = waveformAttr.dt*(i%frameSize);
168197
else
@@ -172,10 +201,10 @@ int main(int argc, char **argv) {
172201
for(iCh=0; iCh<SCOPE_NCH; iCh++) {
173202
if((1<<iCh) & waveformAttr.chMask) {
174203
if(i < waveformFile->nPt){
175-
channel[iCh][i] = (waveformBuf[j * waveformFile->nPt + i]
204+
channel[iCh][i] = ((waveformBuf[j * waveformFile->nPt + i]
176205
- waveformAttr.yoff[iCh])
177206
* waveformAttr.ymult[iCh]
178-
+ waveformAttr.yzero[iCh];
207+
+ waveformAttr.yzero[iCh]) / ADC;
179208
} else {
180209
channel[iCh][i] = 0.;
181210
}
@@ -191,6 +220,103 @@ int main(int argc, char **argv) {
191220
// printf("\n");
192221
}
193222
// printf("\n");
223+
224+
// loop over channels for pulse processing
225+
for(iCh=0; iCh<SCOPE_NCH; iCh++) {
226+
if((1<<iCh) & waveformAttr.chMask) {
227+
xmin[iCh] = 0.;
228+
amp [iCh] = 0.;
229+
base[iCh] = 0.;
230+
integral[iCh] = 0.;
231+
integralFull[iCh] = 0.;
232+
gauspeak[iCh] = 0.;
233+
sigmoidTime[iCh] = 0.;
234+
linearTime0[iCh] = 0.;
235+
linearTime15[iCh] = 0.;
236+
linearTime30[iCh] = 0.;
237+
linearTime45[iCh] = 0.;
238+
linearTime60[iCh] = 0.;
239+
risetime[iCh] = 0.;
240+
constantThresholdTime[iCh] = 0.;
241+
242+
// Make pulse shape graph
243+
TString pulseName = Form("pulse_event%d_ch%d", event, iCh);
244+
TGraphErrors* pulse = GetTGraph( channel[iCh], time[0] );
245+
246+
// Estimate baseline
247+
float baseline;
248+
baseline = GetBaseline( pulse, 5 ,150, pulseName );
249+
base[iCh] = baseline;
250+
251+
// Correct pulse shape for baseline offset + amp/att
252+
for(int j = 0; j < NSAMPLES; j++) {
253+
channel[iCh][j] = channel[iCh][j] + baseline;
254+
}
255+
256+
int index_min = FindMinAbsolute(NSAMPLES, channel[iCh]);
257+
258+
// Recreate the pulse TGraph using baseline-subtracted channel data
259+
delete pulse;
260+
pulse = GetTGraph( channel[iCh], time[0] );
261+
262+
xmin[iCh] = index_min;
263+
264+
Double_t tmpAmp = 0.0;
265+
Double_t tmpMin = 0.0;
266+
pulse->GetPoint(index_min, tmpMin, tmpAmp);
267+
amp[iCh] = tmpAmp * ADC;
268+
269+
// Get pulse integral
270+
if ( xmin[iCh] != 0 ) {
271+
integral[iCh] = GetPulseIntegral( index_min, 20, channel[iCh], time[0] );
272+
integralFull[iCh] = GetPulseIntegral( index_min , channel[iCh], "full");
273+
}
274+
else {
275+
integral[iCh] = 0.0;
276+
integralFull[iCh] = 0.0;
277+
}
278+
279+
// Gaussian time stamp and constant-fraction fit
280+
Double_t min = 0.; Double_t low_edge = 0.; Double_t high_edge = 0.; Double_t y = 0.;
281+
pulse->GetPoint(index_min, min, y);
282+
pulse->GetPoint(index_min-4, low_edge, y); // time of the low edge of the fit range
283+
pulse->GetPoint(index_min+4, high_edge, y); // time of the upper edge of the fit range
284+
285+
float timepeak = 0;
286+
float fs[6]; // constant-fraction fit output
287+
float fs_falling[6]; // falling exp timestapms
288+
float cft_low_range = 0.03;
289+
float cft_high_range = 0.20;
290+
291+
if(xmin[iCh] != 0.0) {
292+
timepeak = GausFit_MeanTime(pulse, low_edge, high_edge);
293+
RisingEdgeFitTime( pulse, index_min, cft_low_range, cft_high_range, fs, event, "linearFit_" + pulseName, false );
294+
}
295+
296+
_isRinging[iCh] = isRinging( index_min, channel[iCh] );
297+
// for output tree
298+
gauspeak[iCh] = timepeak;
299+
risetime[iCh] = fs[0];
300+
linearTime0[iCh] = fs[1];
301+
linearTime15[iCh] = fs[2];
302+
linearTime30[iCh] = fs[3];
303+
linearTime45[iCh] = fs[4];
304+
linearTime60[iCh] = fs[5];
305+
fallingTime[iCh] = fs_falling[0];
306+
constantThresholdTime[iCh] = ConstantThresholdTime( pulse, 75);
307+
308+
// reconvert to "mV" with ADC factor
309+
for(int s = 0; s < NSAMPLES; s++)
310+
channel[iCh][s] *= ADC;
311+
312+
delete pulse;
313+
314+
315+
316+
}
317+
}
318+
319+
194320
tree->Fill();
195321
}
196322

0 commit comments

Comments
 (0)