| 1 | #!/usr/bin/python -tt
|
|---|
| 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 | fbn = path.splitext(file_base_name)[0]
|
|---|
| 24 | fbn = fbn[0:8]+fbn[9:12]
|
|---|
| 25 |
|
|---|
| 26 |
|
|---|
| 27 | baum = TTree('spiketree', 'spike ana tree')
|
|---|
| 28 |
|
|---|
| 29 | # prepare some vars for the tree
|
|---|
| 30 | chid = c_int(0)
|
|---|
| 31 | startcell = c_int(0)
|
|---|
| 32 | number_of_singles = c_int(0)
|
|---|
| 33 | position_of_spikes_in_logical_pipeline = np.zeros(100, np.int32)
|
|---|
| 34 | position_of_spikes_in_physical_pipeline = np.zeros(100, np.int32)
|
|---|
| 35 | time_to_previous = c_int(0)
|
|---|
| 36 | event_id = c_int(0)
|
|---|
| 37 | file_id = c_long(long(fbn))
|
|---|
| 38 |
|
|---|
| 39 | baum.Branch('file',file_id,'file/L')
|
|---|
| 40 | baum.Branch('event',event_id,'event/I')
|
|---|
| 41 | baum.Branch('chid',chid,'chid/I')
|
|---|
| 42 | baum.Branch('sc',startcell,'sc/I')
|
|---|
| 43 | baum.Branch('n',number_of_singles,'n/I')
|
|---|
| 44 | baum.Branch('logpos',position_of_spikes_in_logical_pipeline,'logpos[n]/I')
|
|---|
| 45 | baum.Branch('physpos',position_of_spikes_in_physical_pipeline,'physpos[n]/I')
|
|---|
| 46 | baum.Branch('time',time_to_previous,'time/I')
|
|---|
| 47 |
|
|---|
| 48 |
|
|---|
| 49 | def spikecallback(candidates, singles, doubles, data, ind):
|
|---|
| 50 | if len(singles) >0 :
|
|---|
| 51 | for s in singles:
|
|---|
| 52 | s = np.unravel_index(s, data.shape)
|
|---|
| 53 | hs.Fill( s[0], s[1])
|
|---|
| 54 | if len(doubles) >0 :
|
|---|
| 55 | for d in doubles:
|
|---|
| 56 | d = np.unravel_index(d, data.shape)
|
|---|
| 57 | hd.Fill( d[0], d[1])
|
|---|
| 58 |
|
|---|
| 59 |
|
|---|
| 60 |
|
|---|
| 61 | despike = DRSSpikes(user_action = spikecallback)
|
|---|
| 62 |
|
|---|
| 63 | def mars_spikes( data ):
|
|---|
| 64 | """
|
|---|
| 65 | should search for spikes, just as it is implemented in
|
|---|
| 66 | mcore/DrsCalib.h in DrsCalib::RemoveSpikes
|
|---|
| 67 | static void RemoveSpikes(float *vec, uint32_t roi)
|
|---|
| 68 | {
|
|---|
| 69 | if (roi<4)
|
|---|
| 70 | return;
|
|---|
| 71 |
|
|---|
| 72 | for (size_t ch=0; ch<1440; ch++)
|
|---|
| 73 | {
|
|---|
| 74 | float *p = vec + ch*roi;
|
|---|
| 75 |
|
|---|
| 76 | for (size_t i=1; i<roi-2; i++)
|
|---|
| 77 | {
|
|---|
| 78 | if (p[i]-p[i-1]>25 && p[i]-p[i+1]>25)
|
|---|
| 79 | {
|
|---|
| 80 | p[i] = (p[i-1]+p[i+1])/2;
|
|---|
| 81 | }
|
|---|
| 82 |
|
|---|
| 83 | if (p[i]-p[i-1]>22 && fabs(p[i]-p[i+1])<4 && p[i+1]-p[i+2]>22)
|
|---|
| 84 | {
|
|---|
| 85 | p[i] = (p[i-1]+p[i+2])/2;
|
|---|
| 86 | p[i+1] = p[i];
|
|---|
| 87 | }
|
|---|
| 88 | }
|
|---|
| 89 | }
|
|---|
| 90 | }
|
|---|
| 91 |
|
|---|
| 92 | """
|
|---|
| 93 |
|
|---|
| 94 | #: list, conaining the (chid,slice) tuple of the single spike positions
|
|---|
| 95 | singles = []
|
|---|
| 96 | #: list, conaining the (chid,slice) tuple of the 1st double spike positions
|
|---|
| 97 | doubles = []
|
|---|
| 98 |
|
|---|
| 99 | for chid, pdata in enumerate(data):
|
|---|
| 100 | single_cand = np.where( np.diff(pdata[:-1]) > 25)[0]
|
|---|
| 101 | for cand in single_cand:
|
|---|
| 102 | if -np.diff(pdata[1:])[cand] > 25:
|
|---|
| 103 | singles.append( (chid, cand) )
|
|---|
| 104 |
|
|---|
| 105 | double_cand = np.where( np.diff(pdata[:-1]) > 22 )[0]
|
|---|
| 106 | for cand in double_cand:
|
|---|
| 107 | if cand+1 < len(np.diff(pdata[1:])):
|
|---|
| 108 | if abs(-np.diff(pdata[1:])[cand])<4 and -np.diff(pdata[1:])[cand+1] > 22:
|
|---|
| 109 | doubles.append( (chid, cand) )
|
|---|
| 110 |
|
|---|
| 111 |
|
|---|
| 112 | return singles, doubles
|
|---|
| 113 |
|
|---|
| 114 |
|
|---|
| 115 |
|
|---|
| 116 |
|
|---|
| 117 |
|
|---|
| 118 | event = run.next()
|
|---|
| 119 | board_time = event['board_times'].copy()
|
|---|
| 120 |
|
|---|
| 121 | for event in run:
|
|---|
| 122 | event_id.value = event['event_id']
|
|---|
| 123 |
|
|---|
| 124 |
|
|---|
| 125 | data = event['data']
|
|---|
| 126 | singles, doubles = mars_spikes(data)
|
|---|
| 127 | sc = event['start_cells']
|
|---|
| 128 | bt_old = board_time.copy()
|
|---|
| 129 | board_time = event['board_times'].copy()
|
|---|
| 130 | # print fbn, event['event_id']
|
|---|
| 131 |
|
|---|
| 132 | ss = []
|
|---|
| 133 | for i in range(1440):
|
|---|
| 134 | ss.append([])
|
|---|
| 135 | for s in singles:
|
|---|
| 136 | ss[s[0]].append(s[1])
|
|---|
| 137 |
|
|---|
| 138 | for pixel_id,pixel in enumerate(ss):
|
|---|
| 139 | if len(pixel) >0 :
|
|---|
| 140 | chid.value = pixel_id
|
|---|
| 141 | startcell.value = sc[ chid.value ]
|
|---|
| 142 | number_of_singles.value = len(pixel)
|
|---|
| 143 | time_to_previous.value = board_time[pixel_id/36]-bt_old[pixel_id/36]
|
|---|
| 144 | for num,s in enumerate(pixel):
|
|---|
| 145 | log = s
|
|---|
| 146 | phys = (startcell.value+log)%1024
|
|---|
| 147 | position_of_spikes_in_logical_pipeline[num] = log
|
|---|
| 148 | position_of_spikes_in_physical_pipeline[num] = phys
|
|---|
| 149 | baum.Fill()
|
|---|
| 150 |
|
|---|
| 151 |
|
|---|
| 152 | baum.Write()
|
|---|
| 153 |
|
|---|
| 154 | rootfile.Close()
|
|---|
