Index: /fact/tools/pyscripts/sandbox/dneise/spikes/ana.py =================================================================== --- /fact/tools/pyscripts/sandbox/dneise/spikes/ana.py (revision 14170) +++ /fact/tools/pyscripts/sandbox/dneise/spikes/ana.py (revision 14171) @@ -5,4 +5,5 @@ import sys import numpy as np +from ctypes import * import os import os.path as path @@ -12,94 +13,29 @@ dfn = sys.argv[1] cfn = sys.argv[2] +run = RawData( dfn, cfn) + file_base_name = path.splitext(path.basename(dfn))[0] -root_filename = '/home_nfs/isdc/neise/' + file_base_name + time.strftime('%Y%m%d_%H%M%S') +'.root' - -run = RawData( dfn, cfn) +root_filename = '/home_nfs/isdc/neise/' + file_base_name + '_' + time.strftime('%Y%m%d_%H%M%S') +'_spikeana.root' rootfile = TFile(root_filename, "RECREATE") -canv = TCanvas('canv', 'test canvas',10,10,400,400) -canv.Divide(2,4) +baum = TTree('spiketree', 'spike ana tree') + +# prepare some vars for the tree +chid = c_int(0) +startcell = c_int(0) +number_of_singles = c_int(0) +position_of_spikes_in_logical_pipeline = np.zeros(100, np.int32) +position_of_spikes_in_physical_pipeline = np.zeros(100, np.int32) +time_to_previous = c_long(0) - -hs_log_chid = TH2F("hs_log_chid", - "singles - logical pipe vs chid", - 1440, -0.5, 1439.5, - 300, -0.5, 299.5) -hd_log_chid = TH2F("hd_log_chid", - "doubles - logical pipe vs chid", - 1440, -0.5, 1439.5, - 300, -0.5, 299.5) - -hs_phys_log = TH2F("hs_phys_log", - "singles logical pipe vs physical pipe", - 1024, -0.5, 1023.5, - 300, -0.5, 299.5) - -hd_phys_log = TH2F("hd_phys_log", - "doubles logical pipe vs physical pipe", - 1024, -0.5, 1023.5, - 300, -0.5, 299.5) - -hs_phys_sc = TH2F("hs_phys_sc", - "singles physical pipe vs startcell", - 1024, -0.5, 1023.5, - 1024, -0.5, 1023.5) - -hd_phys_sc = TH2F("hd_phys_sc", - "doubles physical pipe vs startcell", - 1024, -0.5, 1023.5, - 1024, -0.5, 1023.5) - -hs_log_sc = TH2F("hs_log_sc", - "singles logical pipe vs startcell", - 1024, -0.5, 1023.5, - 300, -0.5, 299.5 - ) - -hd_log_sc = TH2F("hd_log_sc", - "doubles logical pipe vs startcell", - 1024, -0.5, 1023.5, - 300, -0.5, 299.5, - ) - -rootfile.mkdir('per_chid') -rootfile.cd('per_chid') - -hd_log_sc_pp = [] -for chid in range(1440): - hd_log_sc_pp.append( - TH2F("chid_"+str(chid), - "doubles - log pipe vs startcell - chid:" +str(chid), - 1024, -0.5, 1023.5, - 300, -0.5, 299.5) - ) - -rootfile.cd('..') - - -singlebaum = TTree('singlebaum', 'spike ana tree') -doublebaum = TTree('doublebaum', 'spike ana tree') - - -# prepare some vars for the tree -chid = 0 -startcell = 0 -number_of_singles = 0 -number_of_doubles = 0 -single_positions = np.zeros(10) -double_positions = np.zeros(10) -singlebaum.Branch('chid',chid,'chid/I') -singlebaum.Branch('sc',startcell,'sc/I') -singlebaum.Branch('n',number_of_singles,'n/I') -singlebaum.Branch('spikes',single_positions,'spikes/I') - -doublebaum.Branch('chid',chid,'chid/I') -doublebaum.Branch('sc',startcell,'sc/I') -doublebaum.Branch('n',number_of_doubles,'n/I') -doublebaum.Branch('spikes',double_positions,'spikes/I') - +baum.Branch('chid',chid,'chid/I') +baum.Branch('sc',startcell,'sc/I') +baum.Branch('n',number_of_singles,'n/I') +baum.Branch('logpos',position_of_spikes_in_logical_pipeline,'logpos[n]/I') +baum.Branch('physpos',position_of_spikes_in_physical_pipeline,'physpos[n]/I') +baum.Branch('time',time_to_previous,'time/I') @@ -116,4 +52,6 @@ +despike = DRSSpikes(user_action = spikecallback) + def mars_spikes( data ): """ @@ -167,5 +105,10 @@ return singles, doubles -despike = DRSSpikes(user_action = spikecallback) + + + + +event = run.next() +bt_old = event['board_times'].copy for event in run: @@ -173,85 +116,21 @@ s, d = mars_spikes(data) sc = event['start_cells'] + bt = event['board_times'].copy() + + if len(s) >0 : + chid.value = s[0][0] + startcell.value = sc[ chid.value ] + number_of_singles.value = len(s) + time.value = bt[s[0][0]/9]-bt_old[s[0][0]/9] + for i in s: + log = i[1] + phys = (startcell.value+log)%1024 + position_of_spikes_in_logical_pipeline[i] = log + position_of_spikes_in_physical_pipeline[i] = phys + baum.Fill() + + bt_old = event['board_times'].copy() - if len(s) >0 : - for i in s: - log = i[1] - chid = i[0] - stace = sc[ chid ] - phys = (stace+log)%1024 - - hs_log_chid.Fill( chid, log) - hs_phys_log.Fill( phys, log) - hs_phys_sc.Fill( stace, phys) - hs_log_sc.Fill( stace, log) - - - if len(d) >0 : - for i in d: - log = i[1] - chid = i[0] - stace = sc[ chid ] - phys = (stace+log)%1024 +baum.Write() - hd_log_chid.Fill( chid, log) - hd_phys_log.Fill( phys, log) - hd_phys_sc.Fill( stace, phys) - hd_log_sc.Fill( stace, log) - hd_log_sc_pp[chid].Fill(stace, log) - - - - if event['event_id'] % 100 == 0: - print event['event_id'] - canv.cd(1) - hs_log_chid.Draw("COLZ") - canv.cd(2) - hd_log_chid.Draw("COLZ") - canv.cd(3) - hs_phys_log.Draw("COLZ") - canv.cd(4) - hd_phys_log.Draw("COLZ") - - canv.cd(5) - hs_phys_sc.Draw("COLZ") - canv.cd(6) - hd_phys_sc.Draw("COLZ") - canv.cd(7) - hs_log_sc.Draw("COLZ") - canv.cd(8) - hd_log_sc.Draw("COLZ") - canv.Modified() - canv.Update() - - - - if event['event_id'] % 1000 == 0: - print event['event_id'] - hs_log_chid.Write('',TObject.kOverwrite) - hd_log_chid.Write('',TObject.kOverwrite) - hs_phys_log.Write('',TObject.kOverwrite) - hd_phys_log.Write('',TObject.kOverwrite) - hs_phys_sc.Write('',TObject.kOverwrite) - hd_phys_sc.Write('',TObject.kOverwrite) - hs_log_sc.Write('',TObject.kOverwrite) - hd_log_sc.Write('',TObject.kOverwrite) - - rootfile.cd('per_chid') - for h in hd_log_sc_pp: - h.Write('',TObject.kOverwrite) - rootfile.cd('..') - -hs_log_chid.Write('',TObject.kOverwrite) -hd_log_chid.Write('',TObject.kOverwrite) -hs_phys_log.Write('',TObject.kOverwrite) -hd_phys_log.Write('',TObject.kOverwrite) -hs_phys_sc.Write('',TObject.kOverwrite) -hd_phys_sc.Write('',TObject.kOverwrite) -hs_log_sc.Write('',TObject.kOverwrite) -hd_log_sc.Write('',TObject.kOverwrite) - -rootfile.cd('per_chid') -for h in hd_log_sc_pp: - h.Write() -rootfile.cd('..') rootfile.Close() Index: /fact/tools/pyscripts/sandbox/dneise/spikes/ana_old.py =================================================================== --- /fact/tools/pyscripts/sandbox/dneise/spikes/ana_old.py (revision 14171) +++ /fact/tools/pyscripts/sandbox/dneise/spikes/ana_old.py (revision 14171) @@ -0,0 +1,279 @@ +#!/usr/bin/python -tti + +from pyfact import RawData +from drs_spikes import DRSSpikes +import sys +import numpy as np +from ctypes import * +import os +import os.path as path +from ROOT import TFile, TCanvas, TH2F, TTree, TStyle, TObject +import time + +dfn = sys.argv[1] +cfn = sys.argv[2] +run = RawData( dfn, cfn) + + +file_base_name = path.splitext(path.basename(dfn))[0] +root_filename = '/home_nfs/isdc/neise/' + file_base_name + '_' + time.strftime('%Y%m%d_%H%M%S') +'_spikeana.root' + +rootfile = TFile(root_filename, "RECREATE") + +singlebaum = TTree('singlebaum', 'spike ana tree') + +# prepare some vars for the tree +chid = c_int(0) +startcell = c_int(0) +number_of_singles = c_int(0) +position_of_spikes_in_logical_pipeline = np.zeros(100, np.int32) +position_of_spikes_in_physical_pipeline = np.zeros(100, np.int32) +time_to_previous = c_long(0) + +singlebaum.Branch('chid',chid,'chid/I') +singlebaum.Branch('sc',startcell,'sc/I') +singlebaum.Branch('n',number_of_singles,'n/I') +singlebaum.Branch('logpos',position_of_spikes_in_logical_pipeline,'logpos[n]/I') +singlebaum.Branch('physpos',position_of_spikes_in_physical_pipeline,'physpos[n]/I') +singlebaum.Branch('time',time_to_previous,'time/I') + + +def spikecallback(candidates, singles, doubles, data, ind): + if len(singles) >0 : + for s in singles: + s = np.unravel_index(s, data.shape) + hs.Fill( s[0], s[1]) + if len(doubles) >0 : + for d in doubles: + d = np.unravel_index(d, data.shape) + hd.Fill( d[0], d[1]) + + + +despike = DRSSpikes(user_action = spikecallback) + +def mars_spikes( data ): + """ + should search for spikes, just as it is implemented in + mcore/DrsCalib.h in DrsCalib::RemoveSpikes + static void RemoveSpikes(float *vec, uint32_t roi) + { + if (roi<4) + return; + + for (size_t ch=0; ch<1440; ch++) + { + float *p = vec + ch*roi; + + for (size_t i=1; i25 && p[i]-p[i+1]>25) + { + p[i] = (p[i-1]+p[i+1])/2; + } + + if (p[i]-p[i-1]>22 && fabs(p[i]-p[i+1])<4 && p[i+1]-p[i+2]>22) + { + p[i] = (p[i-1]+p[i+2])/2; + p[i+1] = p[i]; + } + } + } + } + + """ + + #: list, conaining the (chid,slice) tuple of the single spike positions + singles = [] + #: list, conaining the (chid,slice) tuple of the 1st double spike positions + doubles = [] + + for chid, pdata in enumerate(data): + single_cand = np.where( np.diff(pdata[:-1]) > 25)[0] + for cand in single_cand: + if -np.diff(pdata[1:])[cand] > 25: + singles.append( (chid, cand) ) + + double_cand = np.where( np.diff(pdata[:-1]) > 22 )[0] + for cand in double_cand: + if cand+1 < len(np.diff(pdata[1:])): + if abs(-np.diff(pdata[1:])[cand])<4 and -np.diff(pdata[1:])[cand+1] > 22: + doubles.append( (chid, cand) ) + + + return singles, doubles + + + + +for event in run: + data = event['data'] + s, d = mars_spikes(data) + sc = event['start_cells'] + bt = event['board_times'] + + if len(s) >0 : + chid.value = s[0][0] + startcell.value = sc[ chid.value ] + number_of_singles.value = len(s) + for i in s: + log = i[1] + phys = (startcell.value+log)%1024 + + hs_log_chid.Fill( chid, log) + hs_phys_log.Fill( phys, log) + hs_phys_sc.Fill( stace, phys) + hs_log_sc.Fill( stace, log) + + + if len(d) >0 : + for i in d: + log = i[1] + chid = i[0] + stace = sc[ chid ] + phys = (stace+log)%1024 + + hd_log_chid.Fill( chid, log) + hd_phys_log.Fill( phys, log) + hd_phys_sc.Fill( stace, phys) + hd_log_sc.Fill( stace, log) + hd_log_sc_pp[chid].Fill(stace, log) + + + + if event['event_id'] % 100 == 0: + print event['event_id'] + canv.cd(1) + hs_log_chid.Draw("COLZ") + canv.cd(2) + hd_log_chid.Draw("COLZ") + canv.cd(3) + hs_phys_log.Draw("COLZ") + canv.cd(4) + hd_phys_log.Draw("COLZ") + + canv.cd(5) + hs_phys_sc.Draw("COLZ") + canv.cd(6) + hd_phys_sc.Draw("COLZ") + canv.cd(7) + hs_log_sc.Draw("COLZ") + canv.cd(8) + hd_log_sc.Draw("COLZ") + canv.Modified() + canv.Update() + + + + if event['event_id'] % 1000 == 0: + print event['event_id'] + hs_log_chid.Write('',TObject.kOverwrite) + hd_log_chid.Write('',TObject.kOverwrite) + hs_phys_log.Write('',TObject.kOverwrite) + hd_phys_log.Write('',TObject.kOverwrite) + hs_phys_sc.Write('',TObject.kOverwrite) + hd_phys_sc.Write('',TObject.kOverwrite) + hs_log_sc.Write('',TObject.kOverwrite) + hd_log_sc.Write('',TObject.kOverwrite) + + rootfile.cd('per_chid') + for h in hd_log_sc_pp: + h.Write('',TObject.kOverwrite) + rootfile.cd('..') + +hs_log_chid.Write('',TObject.kOverwrite) +hd_log_chid.Write('',TObject.kOverwrite) +hs_phys_log.Write('',TObject.kOverwrite) +hd_phys_log.Write('',TObject.kOverwrite) +hs_phys_sc.Write('',TObject.kOverwrite) +hd_phys_sc.Write('',TObject.kOverwrite) +hs_log_sc.Write('',TObject.kOverwrite) +hd_log_sc.Write('',TObject.kOverwrite) + +rootfile.cd('per_chid') +for h in hd_log_sc_pp: + h.Write() +rootfile.cd('..') +rootfile.Close() + + +def book_histos(): + + histograms = {} + histograms['overview'] = {} + ov = histograms['overview'] + histograms['per_pixel'] = {} + pp = histograms['per_pixel'] + + ov['hs_log_chid'] = TH2F("hs_log_chid", + "singles", + 1440, -0.5, 1439.5, + 300, -0.5, 299.5) + ov['hs_log_chid'].GetXaxis().SetTitle("cont hardware ID (CHID)") + ov['hs_log_chid'].GetYaxis().SetTitle("spike position in the logical pipeline") + + ov['hd_log_chid'] = TH2F("hd_log_chid", + "doubles", + 1440, -0.5, 1439.5, + 300, -0.5, 299.5) + ov['hd_log_chid'].GetXaxis().SetTitle("cont hardware ID (CHID)") + ov['hd_log_chid'].GetYaxis().SetTitle("spike position in the logical pipeline") + + + ov['hs_phys_log'] = TH2F("hs_phys_log", + "singles", + 1024, -0.5, 1023.5, + 300, -0.5, 299.5) + ov['hs_phys_log'].GetXaxis().SetTitle("spike position in the physical pipeline") + ov['hs_phys_log'].GetYaxis().SetTitle("spike position in the logical pipeline") + + ov['hd_phys_log'] = TH2F("hd_phys_log", + "doubles", + 1024, -0.5, 1023.5, + 300, -0.5, 299.5) + ov['hd_phys_log'].GetXaxis().SetTitle("spike position in the physical pipeline") + ov['hd_phys_log'].GetYaxis().SetTitle("spike position in the logical pipeline") + + + ov['hs_phys_sc'] = TH2F("hs_phys_sc", + "singles", + 1024, -0.5, 1023.5, + 1024, -0.5, 1023.5) + ov['hs_phys_sc'].GetXaxis().SetTitle("DRS startcell") + ov['hs_phys_sc'].GetYaxis().SetTitle("spike position in the physical pipeline") + + ov['hd_phys_sc'] = TH2F("hd_phys_sc", + "doubles", + 1024, -0.5, 1023.5, + 1024, -0.5, 1023.5) + ov['hd_phys_sc'].GetXaxis().SetTitle("DRS startcell") + ov['hd_phys_sc'].GetYaxis().SetTitle("spike position in the physical pipeline") + + + + ov['hs_log_sc'] = TH2F("hs_log_sc", + "singles", + 1024, -0.5, 1023.5, + 1024, -0.5, 1023.5) + ov['hs_log_sc'].GetXaxis().SetTitle("DRS startcell") + ov['hs_log_sc'].GetYaxis().SetTitle("spike position in the logical pipeline") + + ov['hd_log_sc'] = TH2F("hd_log_sc", + "doubles", + 1024, -0.5, 1023.5, + 1024, -0.5, 1023.5) + ov['hd_log_sc'].GetXaxis().SetTitle("DRS startcell") + ov['hd_log_sc'].GetYaxis().SetTitle("spike position in the logical pipeline") + + + hd_log_sc_pp = [] + for chid in range(1440): + hd_log_sc_pp.append( + TH2F("chid_"+str(chid), + "doubles - log pipe vs startcell - chid:" +str(chid), + 1024, -0.5, 1023.5, + 300, -0.5, 299.5) + ) + + +