Index: /fact/tools/pyscripts/sandbox/dneise/spikes/ana.py =================================================================== --- /fact/tools/pyscripts/sandbox/dneise/spikes/ana.py (revision 14157) +++ /fact/tools/pyscripts/sandbox/dneise/spikes/ana.py (revision 14157) @@ -0,0 +1,257 @@ +#!/usr/bin/python -tti + +from pyfact import RawData +from drs_spikes import DRSSpikes +import sys +import numpy as np +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] + +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) + +rootfile = TFile(root_filename, "RECREATE") + +canv = TCanvas('canv', 'test canvas',10,10,400,400) +canv.Divide(2,4) + + + +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') + + + +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]) + + + +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 + +despike = DRSSpikes(user_action = spikecallback) + +for event in run: + data = event['data'] + s, d = mars_spikes(data) + sc = event['start_cells'] + + 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 + + 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()