Guest User

Untitled

a guest
Nov 18th, 2018
87
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 2.75 KB | None | 0 0
  1. #!/usr/bin/env python
  2. from KDataPy.database import kdatadb
  3. import KDataPy.util
  4. import ROOT
  5. import scipy.signal as sig
  6. import numpy as np
  7. import matplotlib.pyplot as plt
  8. plt.ion()
  9.  
  10. #acquire data location from the database
  11. db = kdatadb()
  12. viewResults = db.view('proc/raw', key = 'mf10a002') #returns the location of all RAW .root files for this run
  13.  
  14. #create a KPulseAnalysisChain to remove the baseline and apply a highpass filter
  15. bas = ROOT.KBaselineRemoval() #create a new baseline removal object
  16. nyquist = 0.5/(2.016e-3) #the nyquist frequency for typical heat channels
  17. (b,a) = sig.iirfilter(2, 2.5/nyquist, btype='highpass') #b/a parameters for a second order butterworth highpass filter @ 2.5 Hz
  18. myFilter = ROOT.KIIRSecondOrder(a[1], a[2], b[0], b[1], b[2])
  19. (b,a) = sig.iirfilter(2, 40.0/nyquist, btype='lowpass') #b/a parameters for a second order butterworth lowpass filter @ 40.0 Hz
  20. myFilterLowPass = ROOT.KIIRSecondOrder(a[1], a[2], b[0], b[1], b[2])
  21. chain = ROOT.KPulseAnalysisChain()
  22. chain.AddProcessor(bas)
  23. chain.AddProcessor(myFilter)
  24. chain.AddProcessor(myFilterLowPass)
  25.  
  26. pulseWidths = [] #i will store my results in this array. i could use a ROOT Tree, if I want to save the results for later
  27. minVals = []
  28. minPos = []
  29. #my analysis function is defined here - which will be called by the 'looppulse' function below.
  30. def myAnalysis(pulse,pta,**kwargs):
  31. global pulseWidths # type this in order to access the results array initiliazed above
  32. global minVals
  33. outputPulse = KDataPy.util.get_out(pta) #this gives you the output of the KPulseAnalysisChain defined above
  34.  
  35. #assume a negative going pulse
  36. minval = np.amin(outputPulse)
  37. minpos = np.argmin(outputPulse)
  38. minPos.append(minpos)
  39.  
  40. #find the 90% position
  41. ninPos = 0
  42. for i in range(minpos+1):
  43. if outputPulse[i] <= 0.9 * minval: #remember, i'm assuming that min is a negative value
  44. ninPos = i
  45. break
  46.  
  47. #find the 10% position
  48. tenPos = ninPos
  49. for i in range(ninPos, len(outputPulse)):
  50. if outputPulse[i] >= 0.1 * minval:
  51. tenPos = i
  52. break
  53.  
  54. pulseWidths.append(tenPos - ninPos)
  55. minVals.append(minval)
  56. #print 'position', minpos, 'width', tenPos-ninPos, 'amplitude', minval
  57.  
  58.  
  59. #now loop through the files returned by the database
  60. for row in viewResults:
  61. doc = db[row['id']]
  62. filename = doc['proc1']['file'] # full path to each file of interest
  63.  
  64. #find each H1 ZM1 pulse in the file, apply the "chain" and then pass the chain to your myAnalysis function
  65. #KDataPy.util.plotpulse(filename, name="H1 ZM1", match=True, pta=chain, analysisFunction = myAnalysis)
  66. KDataPy.util.looppulse(filename, name="H1 ZM1", match=True, pta=chain, analysisFunction = myAnalysis)
  67.  
  68.  
  69. #analyze the results
  70. plt.scatter(np.array(minVals), np.array(pulseWidths))
  71.  
  72. raw_input('hit the enter key to quit...')
Add Comment
Please, Sign In to add comment