gistfile1.txt

from ROOT import *
import numpy as np
from scipy import stats
from scipy import signal

def fill(h,v):
    h.SetBins(len(v), 0 , len(v))
    for i in range(len(v)):
        h.SetBinContent(i+1, v[i])


f = KDataReader('/sps/edelweis/kdata/data/run15/raw/lf16b009_008.root')
e = f.GetEvent()

h = TH1D()
hh = TH1D()
hhh = TH1D()

bas = KBaselineRemoval()
pat = KPatternRemoval()
ana = KPulseAnalysisChain()
ana.AddProcessor(bas)
ana.AddProcessor(pat)

b = e.GetBolo(0)
p = b.GetPulseRecord(4)  #this is an ion channel pulse.

if p.GetIsHeatPulse(): 
  pat.SetPatternLength(0)
else:
  pat.SetPatternLength( int(p.GetHeatPulseStampWidth()) )

ana.SetInputPulse( p.GetTrace() )

ana.RunProcess()

fill(h, bas.GetOutputPulse())

filterHz = 50.0	
fNy = p.GetNyquistFrequency()  #scipy filter design requires specifying the filter frequency relative to fsample/2
(b,a) = signal.butter(3, filterHz/fNy, 'low')
zi = signal.lfilter_zi(b, a)


v = np.array(ana.GetOutputPulse()) #this is the baseline / pattern removed pulse stored in a fancy numpy array

(vv,zo)  = signal.lfilter(b,a,v, zi=zi*v[0])

fill(hh,vv)

filterHz = 50.0	
(b,a) = signal.butter(3, filterHz/fNy, 'high')
zi = signal.lfilter_zi(b, a)

(vvv,zzo) = signal.lfilter(b,a,v, zi=zi*v[0])

fill(hhh,vvv)
h.Draw()
hh.SetLineColor(kRed)
hhh.SetLineColor(kBlue)
hh.Draw('same')
hhh.Draw('same')