diff --git a/database/lappd_raw_dist.json b/database/lappd_raw_dist.json index 148fc34d0642b9100b1f76e0f1fee539177ee0dc..29f8f9ce619e38ef9d7645fa373e635d973ce6c2 100644 --- a/database/lappd_raw_dist.json +++ b/database/lappd_raw_dist.json @@ -1,4 +1,44 @@ { + "Calorimeter": { + "plots": [ + { + "title": "Calorimeter Peak Time Distribution", + "func_code": "unpack_time({ch}, true)", + "xlabel": "Time (ns)", + "ylabel": "Counts", + "hist_range": [64, 0, 256], + "path": "ecal_peak_time.png" + }, + { + "title": "Calorimeter Peak Height Distribution", + "func_code": "unpack_height({ch}, true)", + "xlabel": "ADC Value", + "ylabel": "Counts", + "hist_range": [100, 0, 500], + "path": "ecal_peak_height.png" + } + ] + }, + "Scintillator2": { + "plots": [ + { + "title": "Scintillator2 Peak Time Distribution", + "func_code": "unpack_time({ch}, true)", + "xlabel": "Time (ns)", + "ylabel": "Counts", + "hist_range": [64, 0, 256], + "path": "scint2_peak_time.png" + }, + { + "title": "Scintillator2 Peak Height Distribution", + "func_code": "unpack_height({ch}, true)", + "xlabel": "ADC Value", + "ylabel": "Counts", + "hist_range": [100, 0, 500], + "path": "scint2_peak_height.png" + } + ] + }, "LAPPD": { "plots": [ { diff --git a/database/lappd_tcuts_dist.json b/database/lappd_tcuts_dist.json new file mode 100644 index 0000000000000000000000000000000000000000..c3dc6a91462d3e1bf30b5e5dffb81381519cb115 --- /dev/null +++ b/database/lappd_tcuts_dist.json @@ -0,0 +1,25 @@ +{ + "LAPPD": { + "plots": [ + { + "title": "LAPPD Peak Time Distribution", + "func_code": "{ch}_peak.time - ecal.time", + "xlabel": "Time (ns)", + "ylabel": "Counts", + "hist_range": [80, -160, 160], + "path": "tcuts_lappd_peak_time.png", + "filter_code": "{ch}_peak.height > 0" + }, + { + "title": "LAPPD Peak Height Distribution", + "func_code": "{ch}_peak.height", + "xlabel": "ADC Value", + "ylabel": "Counts", + "hist_range": [100, 0, 100], + "path": "tcuts_lappd_peak_height.png", + "filter_code": "{ch}_peak.height > 0" + } + ] + } +} + diff --git a/database/run443_timing.json b/database/run443_timing.json new file mode 100644 index 0000000000000000000000000000000000000000..74e7b089dc29c099e4d0226bd4448237af360b7b --- /dev/null +++ b/database/run443_timing.json @@ -0,0 +1,618 @@ +{ + "C4": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "C6_1": { + "cut": [ + 110.15, + 130.15 + ], + "ref": "0" + }, + "C6_4": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "C6_2": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "C6_3": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "C7_1": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "C7_4": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "C7_2": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "C7_3": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "C5_1": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "C5_4": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "C5_2": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "C5_3": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "C9_1": { + "cut": [ + 110.15, + 130.15 + ], + "ref": "0" + }, + "C9_2": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "C9_4": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "C9_3": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "C8_2": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "C8_3": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "C8_1": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "C8_4": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "C1": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "C2": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "C3": { + "cut": [ + 105.95, + 125.95 + ], + "ref": "0" + }, + "LAPPD_E7": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_E2": { + "cut": [ + -111.75, + -71.75 + ], + "ref": "ecal.time" + }, + "LAPPD_E8": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_E4": { + "cut": [ + -99.75, + -59.75 + ], + "ref": "ecal.time" + }, + "LAPPD_D4": { + "cut": [ + -99.75, + -59.75 + ], + "ref": "ecal.time" + }, + "LAPPD_D8": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_D2": { + "cut": [ + -99.75, + -59.75 + ], + "ref": "ecal.time" + }, + "LAPPD_D7": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_F4": { + "cut": [ + -95.75, + -55.75 + ], + "ref": "ecal.time" + }, + "LAPPD_E5": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_E1": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_E6": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_D6": { + "cut": [ + -99.75, + -59.75 + ], + "ref": "ecal.time" + }, + "LAPPD_D1": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_D5": { + "cut": [ + -95.75, + -55.75 + ], + "ref": "ecal.time" + }, + "LAPPD_C4": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_F8": { + "cut": [ + -95.75, + -55.75 + ], + "ref": "ecal.time" + }, + "LAPPD_F3": { + "cut": [ + -95.75, + -55.75 + ], + "ref": "ecal.time" + }, + "LAPPD_F7": { + "cut": [ + -91.75, + -51.75 + ], + "ref": "ecal.time" + }, + "LAPPD_E3": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_D3": { + "cut": [ + -95.75, + -55.75 + ], + "ref": "ecal.time" + }, + "LAPPD_C7": { + "cut": [ + -95.75, + -55.75 + ], + "ref": "ecal.time" + }, + "LAPPD_C3": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_C8": { + "cut": [ + -95.75, + -55.75 + ], + "ref": "ecal.time" + }, + "LAPPD_F5": { + "cut": [ + -99.75, + -59.75 + ], + "ref": "ecal.time" + }, + "LAPPD_F1": { + "cut": [ + -99.75, + -59.75 + ], + "ref": "ecal.time" + }, + "LAPPD_F6": { + "cut": [ + -91.75, + -51.75 + ], + "ref": "ecal.time" + }, + "LAPPD_F2": { + "cut": [ + -83.75, + -43.75 + ], + "ref": "ecal.time" + }, + "LAPPD_C2": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_C6": { + "cut": [ + -95.75, + -55.75 + ], + "ref": "ecal.time" + }, + "LAPPD_C1": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_C5": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_G4": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_G8": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_G3": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_G7": { + "cut": [ + -119.75, + -79.75 + ], + "ref": "ecal.time" + }, + "LAPPD_B7": { + "cut": [ + -119.75, + -79.75 + ], + "ref": "ecal.time" + }, + "LAPPD_B3": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_B8": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_B4": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_G1": { + "cut": [ + -95.75, + -55.75 + ], + "ref": "ecal.time" + }, + "LAPPD_G6": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_G2": { + "cut": [ + -87.75, + -47.75 + ], + "ref": "ecal.time" + }, + "LAPPD_H6": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_A6": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_B2": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_B6": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_B1": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_G5": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_H4": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_H8": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_H3": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_A3": { + "cut": [ + -119.75, + -79.75 + ], + "ref": "ecal.time" + }, + "LAPPD_A8": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_A4": { + "cut": [ + -75.75, + -35.75 + ], + "ref": "ecal.time" + }, + "LAPPD_B5": { + "cut": [ + -75.75, + -35.75 + ], + "ref": "ecal.time" + }, + "LAPPD_H2": { + "cut": [ + -71.75, + -31.75 + ], + "ref": "ecal.time" + }, + "LAPPD_H7": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_H1": { + "cut": [ + -79.75, + -39.75 + ], + "ref": "ecal.time" + }, + "LAPPD_H5": { + "cut": [ + -115.75, + -75.75 + ], + "ref": "ecal.time" + }, + "LAPPD_A1": { + "cut": [ + -87.75, + -47.75 + ], + "ref": "ecal.time" + }, + "LAPPD_A5": { + "cut": [ + -87.75, + -47.75 + ], + "ref": "ecal.time" + }, + "LAPPD_A7": { + "cut": [ + -75.75, + -35.75 + ], + "ref": "ecal.time" + }, + "LAPPD_A2": { + "cut": [ + -71.75, + -31.75 + ], + "ref": "ecal.time" + } +} \ No newline at end of file diff --git a/scripts/root/dataframe_analysis.C b/scripts/root/dataframe_analysis.C index 22c1828598d2e237494ce5d889a8b5c4d21c21af..e548f8e5b7dffad23dcba1209497743968b08ab0 100644 --- a/scripts/root/dataframe_analysis.C +++ b/scripts/root/dataframe_analysis.C @@ -69,3 +69,19 @@ Peak get_maximum_peak(const std::vector<Peak> &peaks, Filter &&f, Args&&... args return res; } +Peak peak_sum(const std::vector<Peak> &peaks, double thres = 0.) +{ + Peak res(0., 0., 0., 0); + + for (auto &peak : peaks) { + if (peak.height > thres) { + res.time += peak.time; + res.integral += peak.integral; + res.height += peak.height; + res.pos ++; + } + } + res.time /= (double)res.pos; + return res; +} + diff --git a/scripts/signal_timing_cuts.py b/scripts/signal_timing_cuts.py index 6c6d6ba149a0f7de520aa39e8897788bb42b0ffb..0e2e30b148faca5fcef6a09ad3d44425d046e243 100644 --- a/scripts/signal_timing_cuts.py +++ b/scripts/signal_timing_cuts.py @@ -49,7 +49,10 @@ gROOT.SetBatch(args.batch) # root dataframe rdf = ROOT.RDataFrame('EvTree', args.root_file) -save_list = ['ecal', 'ecal_ch', 'scint2', 'scint2_ch'] +save_list = [ + 'ecal', 'ecal_ch', + # 'scint2', 'scint2_ch' +] # read json json_data = dict() diff --git a/scripts/signal_timing_cuts_lappd.py b/scripts/signal_timing_cuts_lappd.py new file mode 100644 index 0000000000000000000000000000000000000000..ad5c3f78edd3ff36c05a759e358b49c07ec880bd --- /dev/null +++ b/scripts/signal_timing_cuts_lappd.py @@ -0,0 +1,173 @@ +import ROOT +import os +import numpy as np +import pandas as pd +import argparse +import json +from matplotlib import pyplot as plt +from ROOT import gROOT, gSystem, gInterpreter, gStyle, gPad, kRed, kAzure +from plot_utils import prepare_canvas, prepare_figure, my_style, get_hist_contents + + +# make sure the relative path for root scripts is correct +owd = os.getcwd() +script_dir = os.path.dirname(os.path.realpath(__file__)) +os.chdir(os.path.join(script_dir, '..')) +print('working directory is ' + os.getcwd()) + +# rootlogon macro +gROOT.Macro('rootlogon.C') + +# the functions that process data are defined in scripts/root/dataframe_analysis.C +func_code = open(os.path.join(script_dir, 'root', 'dataframe_analysis.C')).read() +gInterpreter.Declare(func_code) + +# global style +my_style.cd() +bprops = dict(boxstyle='round', facecolor='wheat', alpha=0.2) + +# argument parser +parser = argparse.ArgumentParser('timing cuts study') +parser.add_argument('root_file', help='root file output from analyzer') +parser.add_argument('json_file', help='output timing cut in json') +parser.add_argument('-b', dest='batch', action='store_true', help='set to batch mode') +parser.add_argument('--read-cuts', dest='read', action='store_true', help='directly use cuts defined in json file.') +parser.add_argument('--snapshot', dest='snap', default='', type=str, help='path to save the data with trigger cuts') +parser.add_argument('--plot-dir', dest='out_dir', default='.', type=str, help='directory to save the plot') +parser.add_argument('--layout', dest='layout', default=os.path.join(script_dir, '..', 'database/channels_layout.json'), + help='json data for channels layout') +args = parser.parse_args() + +# recover paths +for attr in ['root_file', 'json_file', 'snap', 'out_dir', 'layout']: + setattr(args, attr, os.path.join(owd, getattr(args, attr))) + +os.makedirs(args.out_dir, exist_ok=True) + +# batch mode +gROOT.SetBatch(args.batch) + +# root dataframe +rdf = ROOT.RDataFrame('EvTree', args.root_file) +save_list = [ + 'ecal', 'ecal_ch' +] + +# read json +json_data = dict() +if os.path.exists(args.json_file): + with open(args.json_file, 'r') as f: + json_data.update(json.load(f)) + +with open(args.layout, 'r') as f: + channels_layout = json.load(f) + + +# main function +def timing_cuts_search(df, cuts_dict, cut_width, hbins=100): + for ch, cuts in cuts_dict.items(): + search = cuts['cut'] + ref = cuts['ref'] + hist = df.Define('peak_time', 'unpack_time({}, true, {})'.format(ch, ref))\ + .Histo1D(ROOT.RDF.TH1DModel('time_hist1', '', hbins, *search), 'peak_time') + + pos = hist.GetBinCenter(hist.GetMaximumBin()) + height = hist.GetMaximum() + if height > 0: + cuts_dict[ch]['cut'] = (pos - cut_width, pos + cut_width) + return cuts_dict + + +# plot distributions with cuts +def plot_dist_tcuts(df, canvas, pads_dict, cuts_dict, attr, hrange=(128, -256, 256), use_ref=False, hthres=0): + stacks, hists = [], [] + for ch, pad in pads_dict.items(): + cut = cuts_dict[ch]['cut'] + ref = cuts_dict[ch]['ref'] if use_ref else '0' + stacks.append(ROOT.THStack('hs_' + ch, ch)) + attr_name = '{}_{}'.format(ch, attr) + hists.append( + df.Define(attr_name, 'unpack_{}({}, true, {})'.format(attr, ch, ref))\ + .Histo1D(ROOT.RDF.TH1DModel(attr_name + '_h', '', *hrange), attr_name) + ) + hists[-1].SetFillColorAlpha(38, 0.7) + hists[-1].SetFillStyle(3000) + hists[-1].SetLineStyle(2) + hists[-1].SetLineColor(kAzure - 6) + stacks[-1].Add(hists[-1].GetPtr()) + + # cut peaks + hists.append( + df.Filter('{}_peak.height > {}'.format(ch, hthres))\ + .Define(attr_name + '_cut', '{}_peak.{} - {}'.format(ch, attr, ref))\ + .Histo1D(ROOT.RDF.TH1DModel(attr_name + '_hcut', '', *hrange), attr_name + '_cut') + ) + hists[-1].SetFillColorAlpha(46, 0.7) + hists[-1].SetFillStyle(3000) + hists[-1].SetLineStyle(2) + hists[-1].SetLineColor(kRed + 3) + stacks[-1].Add(hists[-1].GetPtr()) + + pad.cd() + stacks[-1].Draw('nostack') + canvas.Update() + return stacks, hists + + + +# ============================================================================= +# MaPMT sum channel timing study +# ============================================================================= +print('Study timing cuts for MaPMT sum channels...') +# calorimeter (trigger) channels +sum_channels = channels_layout['LAPPD'] +# peak search range +if args.read: + sum_cuts = {c: json_data[c] for chs in sum_channels for c in chs} +else: + sum_search = {c: {'cut': (-100, -50), 'ref': 'ecal.time'} for chs in sum_channels for c in chs} + # change sum timing search range for special cases + sum_cuts = timing_cuts_search(rdf, sum_search, cut_width=20) + json_data.update(sum_cuts) + +sum_height_cut = 0 + +# apply cuts +for ch, cut in sum_cuts.items(): + rdf = rdf.Define('{}_peak'.format(ch), 'get_maximum_peak({}.peaks, time_cut_rel, {}, {}, {})'\ + .format(ch, *cut['cut'], cut['ref'])) + save_list.append('{}_peak'.format(ch)) + +# plot timing distributions with t-cuts +c1 = ROOT.TCanvas('sum_time', '', 1920, 1080) +pads1 = prepare_canvas(c1, sum_channels, pad_label='MaPMT Sum Peak Time', + x_label='Time Rel. to ECal. (ns)', y_label='Counts') +res1 = plot_dist_tcuts(rdf, c1, pads1, sum_cuts, 'time', hrange=(80, -160, 160), hthres=sum_height_cut, use_ref=True) +c1.SaveAs(os.path.join(args.out_dir, 'tcuts_sum_time.png')) + +# plot height distributions with t-cuts +c2 = ROOT.TCanvas('sum_height', '', 1920, 1080) +pads2 = prepare_canvas(c2, sum_channels, pad_label='MaPMT Sum Peak Height', x_label='ADC Value', y_label='Counts') +res2 = plot_dist_tcuts(rdf, c2, pads2, sum_cuts, 'height', hrange=(100, 0, 1000), hthres=sum_height_cut) +c2.SaveAs(os.path.join(args.out_dir, 'tcuts_sum_height.png')) + + +# ============================================================================= +# Ending +# ============================================================================= +rdf.Report().Print() +# save json +if not args.read: + with open(args.json_file, 'w') as f: + json.dump(json_data, f, indent=2) + +# save data +if args.snap: + print('Saving data snapshot to \"{}\", with columns:\n{}'.format(args.snap, save_list)) + rdf.Snapshot('EvTree', args.snap, save_list) + +# this keeps the canvas +if not args.batch: + print("press enter to quit") + input() + diff --git a/scripts/trigger_timing_cuts.py b/scripts/trigger_timing_cuts.py index 554bbe03f62b3bbadbabaa77363a52ae77289aab..71e606860fd45bc35ebf5456b4da469d10fee475 100644 --- a/scripts/trigger_timing_cuts.py +++ b/scripts/trigger_timing_cuts.py @@ -49,7 +49,8 @@ gROOT.SetBatch(args.batch) # root dataframe rdf = ROOT.RDataFrame('EvTree', args.root_file) -save_list = ['Cer{}{}_{}'.format(i, j, k) for i in np.arange(1, 5) for j in np.arange(1, 5) for k in np.arange(1, 6)] +# save_list = ['Cer{}{}_{}'.format(i, j, k) for i in np.arange(1, 5) for j in np.arange(1, 5) for k in np.arange(1, 6)] +save_list = ['LAPPD_{}{}'.format(l, i) for i in np.arange(1, 9) for l in 'ABCDEFGH'] # read json json_data = dict() @@ -125,7 +126,7 @@ else: ecal_search = {c: {'cut': (105, 135), 'ref': '0'} for chs in ecal_channels for c in chs} # change ecal timing search range for special cases # ecal_search['C4']['cut'] = (100, 130) - ecal_cuts = timing_cuts_search(rdf, ecal_search, cut_width=20) + ecal_cuts = timing_cuts_search(rdf, ecal_search, cut_width=10) json_data.update(ecal_cuts) ecal_height_cut = 0 @@ -159,6 +160,7 @@ save_list += ['ecal', 'ecal_ch'] +""" # ============================================================================= # Scint channel timing study # ============================================================================= @@ -199,6 +201,7 @@ rdf = rdf.Define('scint2', 'get_maximum_peak({})'.format(s2_vec))\ .Define('scint2_ch', 'std::vector<TString> ch_list{}; return ch_list[scint2.pos];'.format(s2_vec_str))\ .Filter('scint2.height > {}'.format(s2_height_cut), 'Scint2 Signal Cut') save_list += ['scint2', 'scint2_ch'] +"""