#!/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) )