example

example.py

#!/usr/bin/env python

import couchdbkit
import operator
import pywt
import numpy as np
import string, os, sys
import KDataPy.util as kutil
import ROOT
import matplotlib.pyplot as plt
plt.ion()


chan = "ionisA FID807"  #first try a BBv1 channel....

pulseSize = 8192
windowAlpha = 0.1

bbVersion = 1
drawOpt = None

#pulseSize = 512
#windowAlpha = 0.5

afile = '/Users/adam/analysis/edelweiss/data/kdata/raw/ma22a000_000.root'

db = couchdbkit.Server('http://127.0.0.1:5984')['pulsetemplates']
vr = db.view('analytical/bychandate',reduce=False, descending=True, startkey=[chan, "2012-01-22 00:00:00.0"], limit=1, include_docs=True)
doc = vr.first()['doc']
vp = ROOT.std.vector("double")()

exec(doc['formula']['python']) #defines 'template' function

#for simplicity - assume the pulse length here... will have to be fixed up later. 
timeFactor = 1.0
if pulseSize == 512:
    timeFactor = 2.016

for i in range( pulseSize ):
  #print i, template(i*2.016, doc['formula']['par'])
  floatPars = np.array(doc['formula']['par'])
  vp.push_back( template(float(i)*timeFactor, floatPars) )
  

# scaleFactor = 1./abs(min(np.array(vp)))
# print 'scaling by', scaleFactor
# print vp.size(), vp[i], vp[vp.size()-1], vp[i]*scaleFactor

# for i in range(vp.size()):
#   vp[i] = scaleFactor * vp[i]


#plt.plot(np.array(vp))
#plt.show()
#raw_input()
#plt.cla()

cham = ROOT.KChamonixKAmpSite()
#cham.CreateIonWindow(pulseSize)
#cham.CreateHeatWindow(512)
cham.SetTemplate(chan, vp, -1 * int(doc['formula']['par'][0]+ 0.5), pulseSize != 512)
cham.SetNeedScout(False)  #gonna do my own scouting. 


noisePower = None

#define the parameters use in the wavelet analysis
noiseFactor = 1.6
ignoreSpace = 10

hc2p = ROOT.KHalfComplexPower()
r2hc = ROOT.KRealToHalfComplexDFT()


chain = ROOT.KPulseAnalysisChain()

window = None

if pulseSize != 512:
  if bbVersion == 2:
    bas = ROOT.KLinearRemoval()
  else: 
    bas = ROOT.KBaselineRemoval()

  chain.AddProcessor(bas)
  pat = ROOT.KPatternRemoval()
  pat2 = ROOT.KPatternRemoval()
  if bbVersion == 1:
    pat.SetPatternLength(200)
    pat2.SetPatternLength(400)
  else:
    pat.SetPatternLength(10)
    pat2.SetPatternLength(20)
    
  chain.AddProcessor(pat)
  chain.AddProcessor(pat2)
  chain.AddProcessor(cham.GetIonWindow())
  window = cham.GetIonWindow()


else:
  bas = ROOT.KLinearRemoval()
  chain.AddProcessor(bas)
  chain.AddProcessor(cham.GetHeatWindow())
  window = cham.GetHeatWindow()

erapf = ROOT.KEraPeakFinder()
erapf.SetOrder(5)
erapf.SetNumRms(5.0)

polCalc = ROOT.KPulsePolarityCalculator()

def fillNoise(x):  #x is a numpy array
  global noisePower
  global hc2p
  global r2hc

  r2hc.SetInputPulse(x, len(x))
  r2hc.RunProcess()
  hc2p.SetInputPulse(r2hc)
  hc2p.RunProcess()

  if noisePower != None:
    noisePower = np.append(noisePower, [ kutil.get_out(hc2p) ], axis=0)
  else:
    noisePower = np.array([kutil.get_out(hc2p)])

def findNoiseWavelet(pulse, pta=None, **kwargs):


  print pulse.GetChannelName()

  if pta:
    x = kutil.get_out(pta)
  else:
    x = np.array(pulse.GetTrace())
  
  coeffs = pywt.wavedec(x, 'db2', level = kwargs["lvl"])
  
  isEvent = False

  for i in range(2,kwargs["lvl"]):

    pulseMax = np.amax(np.abs(coeffs[i] [ignoreSpace: np.size(coeffs[i]) - ignoreSpace ] ))  #ignore the last number of bins... 
    noiseMax =  np.amax(np.abs( coeffs[i+1][ignoreSpace: np.size(coeffs[i+1]) - ignoreSpace ] ))
    #print i, pulseMax, noiseMax

    #condition = 'np.abs(pulseMax - noiseMax) < %f*noiseMax' % noiseFactor

    if pulseMax > noiseFactor*noiseMax:   #use 2.0 for heat pulses and 1.03 for ionization pulses. 
      isEvent = True
  
  #print 'is event', isEvent
   
  if isEvent == False:
    fillNoise(x)

  if kwargs.has_key("draw") == False: return

  if kwargs["draw"] == isEvent:
    plt.figure(2)
    plt.subplot(kwargs["lvl"]+1,2,1)
    plt.plot(x)

    for i in range(2, kwargs["lvl"]):
      plt.subplot(kwargs["lvl"]+1,2,2*i+3)
      plt.plot(coeffs[i])
      plt.subplot(kwargs["lvl"]+1,2,2*i+4)
      plt.plot(coeffs[i+1])

  
    raw_input()
    plt.clf()  
    

def findNoiseEra(pulse, pta, **kwargs):

  global erapf

  erapf.SetInputPulse(pta)
  #erapf.SetPolarity( polCalc.GetExpectedPolarity(pulse) )
  erapf.SetPolarity( 0 )
  erapf.RunProcess()

  isEvent = True
  if erapf.GetPeakBins().size() == 0:
    isEvent = False;

  #print 'is event', isEvent
   
  if isEvent == False:
    fillNoise( kutil.get_out(pta) )
   

  if kwargs.has_key("draw") == False: return

  if kwargs["draw"] == isEvent:
    plt.figure(2)
    plt.subplot(3,1,1)
    plt.plot( np.array(pulse.GetTrace()) )
    plt.subplot(3,1,2)
    plt.plot(kutil.get_out(pta))
    plt.subplot(3,1,3)
    plt.plot(kutil.get_out(erapf))

    raw_input()
    plt.clf()


#kutil.looppulse(afile, chan, match=True, pta = chain, analysisFunction = findNoiseWavelet, lvl = 6, draw=False)
kutil.looppulse(afile, chan, match=True, pta = chain, analysisFunction = findNoiseEra, lvl = 6, draw=drawOpt)


if noisePower.shape[0]:

  print noisePower.shape

  medianPower = np.zeros(noisePower.shape[1])
  modePower = np.zeros(noisePower.shape[1])
  meanPower = np.zeros(noisePower.shape[1])
  noisePowerTrans = noisePower.transpose()
  noisePowerTrans.shape

  noisePowerTrans.sort()
  print medianPower.size
  for i in range(medianPower.size):
    medianPower[i] = noisePowerTrans[i][ noisePowerTrans[i].size/2 + 1]
    meanPower[i] = np.mean(noisePowerTrans[i])
    (hist,binEdge) = np.histogram(noisePowerTrans[i], 100.0 * noisePower.shape[0], range=[0, medianPower[i]])
    modePower[i] = binEdge[np.argmax(hist)]

  plt.cla()
  print medianPower
  plt.loglog(medianPower)
  plt.loglog(meanPower)
  plt.loglog(modePower)
  window.SetInputPulse(vp)
  window.RunProcess()
  r2hc.SetInputPulse(window)
  r2hc.RunProcess()
  hc2p.SetInputPulse(r2hc)
  hc2p.RunProcess()
  plt.loglog( kutil.get_out(hc2p) )
  plt.show()
  raw_input()

  for i in range(medianPower.size):
    print 'frequency bin', i
    plt.cla()
    plt.hist(noisePowerTrans[i], noisePower.shape[0])
    print 'median value', medianPower[i]
    print 'mean value', meanPower[i]
    print 'mode value', modePower[i]

    plt.show()
    raw_input()

else:
  print 'ahhh'