| 1 | #!/usr/bin/python -tti
|
|---|
| 2 |
|
|---|
| 3 | from pyfact import RawData
|
|---|
| 4 | from drs_spikes import DRSSpikes
|
|---|
| 5 | import sys
|
|---|
| 6 | import numpy as np
|
|---|
| 7 | from ctypes import *
|
|---|
| 8 | import os
|
|---|
| 9 | import os.path as path
|
|---|
| 10 | from ROOT import TFile, TCanvas, TH2F, TTree, TStyle, TObject
|
|---|
| 11 | import time
|
|---|
| 12 |
|
|---|
| 13 | dfn = sys.argv[1]
|
|---|
| 14 | cfn = sys.argv[2]
|
|---|
| 15 | run = RawData( dfn, cfn)
|
|---|
| 16 |
|
|---|
| 17 |
|
|---|
| 18 | file_base_name = path.splitext(path.basename(dfn))[0]
|
|---|
| 19 | root_filename = '/home_nfs/isdc/neise/' + file_base_name + '_' + time.strftime('%Y%m%d_%H%M%S') +'_spikeana.root'
|
|---|
| 20 |
|
|---|
| 21 | rootfile = TFile(root_filename, "RECREATE")
|
|---|
| 22 |
|
|---|
| 23 | baum = TTree('spiketree', 'spike ana tree')
|
|---|
| 24 |
|
|---|
| 25 | # prepare some vars for the tree
|
|---|
| 26 | chid = c_int(0)
|
|---|
| 27 | startcell = c_int(0)
|
|---|
| 28 | number_of_singles = c_int(0)
|
|---|
| 29 | position_of_spikes_in_logical_pipeline = np.zeros(100, np.int32)
|
|---|
| 30 | position_of_spikes_in_physical_pipeline = np.zeros(100, np.int32)
|
|---|
| 31 | time_to_previous = c_long(0)
|
|---|
| 32 |
|
|---|
| 33 |
|
|---|
| 34 | baum.Branch('chid',chid,'chid/I')
|
|---|
| 35 | baum.Branch('sc',startcell,'sc/I')
|
|---|
| 36 | baum.Branch('n',number_of_singles,'n/I')
|
|---|
| 37 | baum.Branch('logpos',position_of_spikes_in_logical_pipeline,'logpos[n]/I')
|
|---|
| 38 | baum.Branch('physpos',position_of_spikes_in_physical_pipeline,'physpos[n]/I')
|
|---|
| 39 | baum.Branch('time',time_to_previous,'time/I')
|
|---|
| 40 |
|
|---|
| 41 |
|
|---|
| 42 | def spikecallback(candidates, singles, doubles, data, ind):
|
|---|
| 43 | if len(singles) >0 :
|
|---|
| 44 | for s in singles:
|
|---|
| 45 | s = np.unravel_index(s, data.shape)
|
|---|
| 46 | hs.Fill( s[0], s[1])
|
|---|
| 47 | if len(doubles) >0 :
|
|---|
| 48 | for d in doubles:
|
|---|
| 49 | d = np.unravel_index(d, data.shape)
|
|---|
| 50 | hd.Fill( d[0], d[1])
|
|---|
| 51 |
|
|---|
| 52 |
|
|---|
| 53 |
|
|---|
| 54 | despike = DRSSpikes(user_action = spikecallback)
|
|---|
| 55 |
|
|---|
| 56 | def mars_spikes( data ):
|
|---|
| 57 | """
|
|---|
| 58 | should search for spikes, just as it is implemented in
|
|---|
| 59 | mcore/DrsCalib.h in DrsCalib::RemoveSpikes
|
|---|
| 60 | static void RemoveSpikes(float *vec, uint32_t roi)
|
|---|
| 61 | {
|
|---|
| 62 | if (roi<4)
|
|---|
| 63 | return;
|
|---|
| 64 |
|
|---|
| 65 | for (size_t ch=0; ch<1440; ch++)
|
|---|
| 66 | {
|
|---|
| 67 | float *p = vec + ch*roi;
|
|---|
| 68 |
|
|---|
| 69 | for (size_t i=1; i<roi-2; i++)
|
|---|
| 70 | {
|
|---|
| 71 | if (p[i]-p[i-1]>25 && p[i]-p[i+1]>25)
|
|---|
| 72 | {
|
|---|
| 73 | p[i] = (p[i-1]+p[i+1])/2;
|
|---|
| 74 | }
|
|---|
| 75 |
|
|---|
| 76 | if (p[i]-p[i-1]>22 && fabs(p[i]-p[i+1])<4 && p[i+1]-p[i+2]>22)
|
|---|
| 77 | {
|
|---|
| 78 | p[i] = (p[i-1]+p[i+2])/2;
|
|---|
| 79 | p[i+1] = p[i];
|
|---|
| 80 | }
|
|---|
| 81 | }
|
|---|
| 82 | }
|
|---|
| 83 | }
|
|---|
| 84 |
|
|---|
| 85 | """
|
|---|
| 86 |
|
|---|
| 87 | #: list, conaining the (chid,slice) tuple of the single spike positions
|
|---|
| 88 | singles = []
|
|---|
| 89 | #: list, conaining the (chid,slice) tuple of the 1st double spike positions
|
|---|
| 90 | doubles = []
|
|---|
| 91 |
|
|---|
| 92 | for chid, pdata in enumerate(data):
|
|---|
| 93 | single_cand = np.where( np.diff(pdata[:-1]) > 25)[0]
|
|---|
| 94 | for cand in single_cand:
|
|---|
| 95 | if -np.diff(pdata[1:])[cand] > 25:
|
|---|
| 96 | singles.append( (chid, cand) )
|
|---|
| 97 |
|
|---|
| 98 | double_cand = np.where( np.diff(pdata[:-1]) > 22 )[0]
|
|---|
| 99 | for cand in double_cand:
|
|---|
| 100 | if cand+1 < len(np.diff(pdata[1:])):
|
|---|
| 101 | if abs(-np.diff(pdata[1:])[cand])<4 and -np.diff(pdata[1:])[cand+1] > 22:
|
|---|
| 102 | doubles.append( (chid, cand) )
|
|---|
| 103 |
|
|---|
| 104 |
|
|---|
| 105 | return singles, doubles
|
|---|
| 106 |
|
|---|
| 107 |
|
|---|
| 108 |
|
|---|
| 109 |
|
|---|
| 110 |
|
|---|
| 111 | event = run.next()
|
|---|
| 112 | bt_old = event['board_times'].copy()
|
|---|
| 113 |
|
|---|
| 114 | for event in run:
|
|---|
| 115 | data = event['data']
|
|---|
| 116 | singles, doubles = mars_spikes(data)
|
|---|
| 117 | sc = event['start_cells']
|
|---|
| 118 | board_time = event['board_times'].copy()
|
|---|
| 119 |
|
|---|
| 120 | ss = []
|
|---|
| 121 | for i in range(1440):
|
|---|
| 122 | ss.append([])
|
|---|
| 123 | for s in singles:
|
|---|
| 124 | ss[s[0]].append(s[1])
|
|---|
| 125 |
|
|---|
| 126 | for pixel_id,pixel in enumerate(ss):
|
|---|
| 127 | if len(pixel) >0 :
|
|---|
| 128 | chid.value = pixel_id
|
|---|
| 129 | startcell.value = sc[ chid.value ]
|
|---|
| 130 | number_of_singles.value = len(pixel)
|
|---|
| 131 | time.value = board_time[pixel_id/36]-bt_old[pixel_id/36]
|
|---|
| 132 | for num,s in enumerate(pixel):
|
|---|
| 133 | log = s
|
|---|
| 134 | phys = (startcell.value+log)%1024
|
|---|
| 135 | position_of_spikes_in_logical_pipeline[num] = log
|
|---|
| 136 | position_of_spikes_in_physical_pipeline[num] = phys
|
|---|
| 137 | baum.Fill()
|
|---|
| 138 |
|
|---|
| 139 | bt_old = event['board_times'].copy()
|
|---|
| 140 |
|
|---|
| 141 | baum.Write()
|
|---|
| 142 |
|
|---|
| 143 | rootfile.Close()
|
|---|
