Skip to content
Snippets Groups Projects
Commit 71d428f8 authored by Chao Peng's avatar Chao Peng
Browse files

clean up the code

parent 143171ac
Branches
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
include/Fadc250Decoder.h
+ 1
− 1
View file @ 71d428f8
...@@ -36,7 +36,7 @@ struct Fadc250Chan ...@@ -36,7 +36,7 @@ struct Fadc250Chan
{ {
bool overflow; bool overflow;
uint32_t win_sum; uint32_t win_sum;
std::vector<uint32_t> integral, amplitude, time, pulse_raw, raw; std::vector<uint32_t> integral, peak, time, pedestal, pulse_raw, raw;
Fadc250Chan() : overflow(false), win_sum(0) {} Fadc250Chan() : overflow(false), win_sum(0) {}
}; };
......
include/utils.h
+ 66
− 0
View file @ 71d428f8
...@@ -6,6 +6,7 @@ ...@@ -6,6 +6,7 @@
#include <functional> #include <functional>
#include <unordered_map> #include <unordered_map>
#include "ConfigParser.h" #include "ConfigParser.h"
#include "ConfigObject.h"
enum ModuleType enum ModuleType
...@@ -89,4 +90,69 @@ std::vector<PulseData> get_pulses(Module *mod, int ch, ModuleType t = kMaxModule ...@@ -89,4 +90,69 @@ std::vector<PulseData> get_pulses(Module *mod, int ch, ModuleType t = kMaxModule
} }
std::vector<Module> read_modules(const std::string &path, bool verbose=false)
{
// read in file
std::string buffer = ConfigParser::file_to_string(path);
// remove comments
ConfigParser::comment_line(buffer, "#", "\n");
// get content blocks
auto text = ConfigParser::break_into_blocks(buffer, "{", "}");
ConfigObject conf;
std::vector<Module> res;
for(auto &b : text.blocks)
{
// module
if (!ConfigParser::case_ins_equal("Module", b.label))
continue;
auto btext = ConfigParser::break_into_blocks(b.content, "{", "}");
conf.ReadConfigString(btext.residual);
// module attributes
Module mod;
mod.crate = conf.GetConfigValue("crate").Int();
mod.slot = conf.GetConfigValue("slot").Int();
mod.type = str2ModuleType(conf.GetConfigValue("type").c_str());
// channels
for (auto &sub : btext.blocks) {
if (!ConfigParser::case_ins_equal("Channels", sub.label))
continue;
ConfigParser parser;
parser.ReadBuffer(sub.content.c_str());
while(parser.ParseLine()) {
mod.channels.emplace_back(Channel{
parser.TakeFirst().Int(),
parser.TakeFirst().String(),
str2ChannelType(parser.TakeFirst().c_str())
});
}
}
res.emplace_back(mod);
}
if (verbose) {
// print out all channels
std::cout << "Read-in " << res.size() << " modules from \"" << path << "\"." << std::endl;
for (auto &mod : res) {
std::cout << "Module: crate = " << mod.crate
<< ", slot = " << mod.slot
<< ", type = " << ModuleType2str(mod.type)
<< std::endl;
for (auto &ch : mod.channels) {
std::cout << "\t Channel " << ch.id
<< ": name = " << ch.name
<< ", type = " << ChannelType2str(ch.type)
<< std::endl;
}
}
}
return res;
}
#endif // MY_UTILS_H #endif // MY_UTILS_H
scripts/analyze_waveform.py
+ 64
− 62
View file @ 71d428f8
...@@ -4,6 +4,7 @@ import json ...@@ -4,6 +4,7 @@ import json
import numpy as np import numpy as np
import pandas as pd import pandas as pd
import argparse import argparse
from collections import OrderedDict
from matplotlib import pyplot as plt from matplotlib import pyplot as plt
from scipy import signal from scipy import signal
from collections.abc import Iterable from collections.abc import Iterable
...@@ -13,14 +14,25 @@ def branch_to_array1d(br, t): ...@@ -13,14 +14,25 @@ def branch_to_array1d(br, t):
return np.ndarray((len(br), ), dtype=t, buffer=br) return np.ndarray((len(br), ), dtype=t, buffer=br)
def get_channels(json_path): def get_maximum_peaks(tree, chan_pos, ref_pos=0.):
dbf = open(json_path) # init
db = json.load(dbf) peaks = np.zeros(len(chan_pos), dtype=np.float32)
channels = dict() poses = np.zeros(len(chan_pos), dtype=np.float32)
for mod, mod_prop in db.items():
for ch in mod_prop['channels']: for i, (c, twin) in enumerate(chan_pos.items()):
channels[ch['name']] = ch['type'] cpeaks = branch_to_array1d(tree.__getattr__(c + '_Ppeak'), np.float32)
return channels cposes = branch_to_array1d(tree.__getattr__(c + '_Ptime'), np.float32) - ref_pos
# find maximum peak inside window
mask = (cposes <= twin[1]) & (cposes >= twin[0])
cpeaks = cpeaks[mask]
cposes = cposes[mask]
if not len(cpeaks):
continue
idx = np.argmax(cpeaks)
if cpeaks[idx] > 1e-5:
peaks[i] = cpeaks[idx]
poses[i] = cposes[idx]
return peaks, poses
parser = argparse.ArgumentParser('Raw waveform analysis') parser = argparse.ArgumentParser('Raw waveform analysis')
...@@ -41,69 +53,59 @@ dbf = open(args.tconfig) ...@@ -41,69 +53,59 @@ dbf = open(args.tconfig)
ch_pos = json.load(dbf) ch_pos = json.load(dbf)
# trigger channels # trigger channels
# tr = ['S2_1', 'S2_2', 'S2_3', 'S2_4', 'S2_5', 'S2_6', 'S2_7', 'S2_8', 'S2_9', 'S2_10', 'S2_11'] # scintillator
# tr_time = (17, 23) sc_ch = ['S2_1', 'S2_2', 'S2_3', 'S2_4', 'S2_5', 'S2_6', 'S2_7', 'S2_8', 'S2_9', 'S2_10', 'S2_11']
tr = ['C1', 'C2', 'C3', 'C4', 'C5_1', 'C5_2', 'C5_3', 'C5_4', 'C6_1', 'C6_2', 'C6_3', 'C6_4', sc_pos = OrderedDict([(ch, (17, 23)) for ch in sc_ch])
'C7_1', 'C7_2', 'C7_3', 'C7_4', 'C8_1', 'C8_2', 'C8_3', 'C8_4', 'C9_1', 'C9_2', 'C9_3', 'C9_4']
tr_time = (28, 32) # calorimeter
ec_ch = ['C1', 'C2', 'C3', 'C4', 'C5_1', 'C5_2', 'C5_3', 'C5_4', 'C6_1', 'C6_2', 'C6_3', 'C6_4',
'C7_1', 'C7_2', 'C7_3', 'C7_4', 'C8_1', 'C8_2', 'C8_3', 'C8_4', 'C9_1', 'C9_2', 'C9_3', 'C9_4']
ec_pos = OrderedDict([(ch, (28, 32)) for ch in ec_ch])
f = ROOT.TFile.Open(args.root_file, 'read') f = ROOT.TFile.Open(args.root_file, 'read')
tree = f.EvTree tree = f.EvTree
nev = tree.GetEntries()
trg_ch = np.ndarray(shape=(tree.GetEntries(), ), dtype=object) if args.nev > 0 and args.nev <= nev:
trg_val = np.ndarray(shape=(tree.GetEntries(), 3), dtype='float64') nev = args.nev
ch_val = np.ndarray(shape=(tree.GetEntries(), len(ch_pos)*2), dtype='float64')
# init buffers
for iev in np.arange(0, tree.GetEntries()): trg_cols = ['peak', 'pos', 'nhits', 'sum']
ec_names = np.ndarray(shape=(nev, ), dtype=object)
ec_vals = np.zeros(shape=(nev, len(trg_cols)), dtype=np.float32)
sc_names = np.ndarray(shape=(nev, ), dtype=object)
sc_vals = np.zeros(shape=(nev, len(trg_cols)), dtype=np.float32)
ch_vals = np.zeros(shape=(nev, len(ch_pos)*2), dtype=np.float32)
for iev in np.arange(0, nev):
tree.GetEntry(iev) tree.GetEntry(iev)
if iev % 1000 == 0: if iev % 1000 == 0:
print('processed {}'.format(iev), end='\r') print('processed {}'.format(iev), end='\r')
# triggers sc_peaks, sc_poses = get_maximum_peaks(tree, sc_pos)
ev_trg_peaks, ev_trg_poses = [], [] ec_peaks, ec_poses = get_maximum_peaks(tree, ec_pos)
for c in tr:
peaks = branch_to_array1d(tree.__getattr__(c + '_Ppeak'), np.float32) # save result to buffer, EC
poses = branch_to_array1d(tree.__getattr__(c + '_Ptime'), np.float32) iec = np.argmax(ec_peaks)
tpeak, tpos = 0, 0 ec_names[iev] = list(ec_pos)[iec]
for peak, pos in zip(peaks, poses): ec_vals[iev] = (ec_peaks[iec], ec_poses[iec], np.sum(ec_peaks > 0), np.sum(ec_peaks))
if pos <= tr_time[1] and pos >= tr_time[0] and peak > tpeak: ref = ec_poses[iec]
tpeak = peak
tpos = pos # SC
ev_trg_peaks.append(tpeak) isc = np.argmax(sc_peaks)
ev_trg_poses.append(tpos) sc_names[iev] = list(sc_pos)[isc]
sc_vals[iev] = (sc_peaks[isc], sc_poses[isc], np.sum(sc_peaks > 0), np.sum(sc_peaks))
# no triggers
if not len(ev_trg_peaks): # CHER
continue ch_peaks, ch_poses = get_maximum_peaks(tree, ch_pos, ref)
ch_vals[iev] = np.concatenate([ch_peaks, ch_poses])
# get the maximum peak from all trigger channels
ich = np.argmax(ev_trg_peaks)
trg_ch[iev] = tr[ich]
trg_val[iev] = (ev_trg_peaks[ich], ev_trg_poses[ich], len(ev_trg_peaks))
# channels
for i, (c, twin) in enumerate(ch_pos.items()):
peaks = branch_to_array1d(tree.__getattr__(c + '_Ppeak'), np.float32)
poses = branch_to_array1d(tree.__getattr__(c + '_Ptime'), np.float32)
# check timing, maximum peak inside the time window
cpeak, cpos = 0, 0
for peak, pos in zip(peaks, poses):
pos -= ev_trg_poses[ich]
if pos >= twin[0] and pos <= twin[1] and peak > cpeak:
cpeak = peak
cpos = pos
ch_val[iev][[i, len(ch_pos) + i]] = (cpeak, cpos)
if args.nev > 0 and iev >= args.nev:
break
print('processed {}'.format(iev)) print('processed {}'.format(iev))
cols = ['trg_peak', 'trg_pos', 'trg_nhits'] \ cols = ['ec_' + col for col in trg_cols] + ['sc_' + col for col in trg_cols] \
+ [c + '_peak' for c, _ in ch_pos.items()] \ + [c + '_peak' for c, _ in ch_pos.items()] + [c + '_pos' for c, _ in ch_pos.items()]
+ [c + '_pos' for c, _ in ch_pos.items()] result = pd.DataFrame(index=np.arange(nev), columns=cols,
result = pd.DataFrame(index=np.arange(iev+1), columns=cols, data=np.concatenate((trg_val[0:iev+1], ch_val[0:iev+1]), axis=1)) data=np.concatenate((ec_vals, sc_vals, ch_vals), axis=1))
result.loc[:, 'trg_ch'] = trg_ch[0:iev+1] result.loc[:, 'ec_ch'] = ec_names
result.loc[:, 'sc_ch'] = sc_names
result.to_csv(args.output) result.to_csv(args.output)
src/analyze.cpp
+ 198
− 301
View file @ 71d428f8
This diff is collapsed.
src/esb_analyze.cpp
+ 3
− 103
View file @ 71d428f8
...@@ -21,8 +21,7 @@ ...@@ -21,8 +21,7 @@
#define PROGRESS_COUNT 1000 #define PROGRESS_COUNT 1000
std::vector<Module> read_modules(const std::string &path); void write_raw_data(const std::string &dpath, const std::string &opath, int nev, const std::vector<Module> &modules);
void write_raw_data(const std::string &dpath, const std::string &opath, int nev, const std::vector<Module> &modules, const std::string &wpath);
bool parseEvent(const uint32_t *buf, bool verbose); bool parseEvent(const uint32_t *buf, bool verbose);
uint32_t parseBlock(const uint32_t *buf); uint32_t parseBlock(const uint32_t *buf);
...@@ -135,7 +134,6 @@ int main(int argc, char* argv[]) ...@@ -135,7 +134,6 @@ int main(int argc, char* argv[])
int nev = -1; int nev = -1;
std::string mconf = "config/esb_module.conf"; std::string mconf = "config/esb_module.conf";
std::string wpath = "";
for (auto &opt : conf_opt.GetOptions()) { for (auto &opt : conf_opt.GetOptions()) {
switch (opt.mark) { switch (opt.mark) {
...@@ -145,9 +143,6 @@ int main(int argc, char* argv[]) ...@@ -145,9 +143,6 @@ int main(int argc, char* argv[])
case 'm': case 'm':
mconf = opt.var.String(); mconf = opt.var.String();
break; break;
case 'w':
wpath = opt.var.String();
break;
default : default :
std::cout << conf_opt.GetInstruction() << std::endl; std::cout << conf_opt.GetInstruction() << std::endl;
return -1; return -1;
...@@ -156,73 +151,7 @@ int main(int argc, char* argv[]) ...@@ -156,73 +151,7 @@ int main(int argc, char* argv[])
auto modules = read_modules(mconf); auto modules = read_modules(mconf);
write_raw_data(conf_opt.GetArgument(0), conf_opt.GetArgument(1), nev, modules, wpath); write_raw_data(conf_opt.GetArgument(0), conf_opt.GetArgument(1), nev, modules);
}
std::vector<Module> read_modules(const std::string &path)
{
// read in file
std::string buffer = ConfigParser::file_to_string(path);
// remove comments
ConfigParser::comment_line(buffer, "#", "\n");
// get content blocks
auto text = ConfigParser::break_into_blocks(buffer, "{", "}");
ConfigObject conf;
std::vector<Module> res;
for(auto &b : text.blocks)
{
// module
if (!ConfigParser::case_ins_equal("Module", b.label))
continue;
auto btext = ConfigParser::break_into_blocks(b.content, "{", "}");
conf.ReadConfigString(btext.residual);
// module attributes
Module mod;
mod.crate = conf.GetConfigValue("crate").Int();
mod.slot = conf.GetConfigValue("slot").Int();
mod.type = str2ModuleType(conf.GetConfigValue("type").c_str());
// channels
for (auto &sub : btext.blocks) {
if (!ConfigParser::case_ins_equal("Channels", sub.label))
continue;
ConfigParser parser;
parser.ReadBuffer(sub.content.c_str());
while(parser.ParseLine()) {
mod.channels.emplace_back(Channel{
parser.TakeFirst().Int(),
parser.TakeFirst().String(),
str2ChannelType(parser.TakeFirst().c_str())
});
}
}
res.emplace_back(mod);
}
/*
// print out all channels
std::cout << "Read-in " << res.size() << " modules from \"" << path << "\"." << std::endl;
for (auto &mod : res) {
std::cout << "Module: crate = " << mod.crate
<< ", slot = " << mod.slot
<< ", type = " << ModuleType2str(mod.type)
<< std::endl;
for (auto &ch : mod.channels) {
std::cout << "\t Channel " << ch.id
<< ": name = " << ch.name
<< ", type = " << ChannelType2str(ch.type)
<< std::endl;
}
}
*/
return res;
} }
...@@ -390,7 +319,7 @@ void fill_tree(TTree *tree, std::unordered_map<std::string, BranchData> &brdata, ...@@ -390,7 +319,7 @@ void fill_tree(TTree *tree, std::unordered_map<std::string, BranchData> &brdata,
} }
void write_raw_data(const std::string &dpath, const std::string &opath, int nev, const std::vector<Module> &modules, const std::string &wpath) void write_raw_data(const std::string &dpath, const std::string &opath, int nev, const std::vector<Module> &modules)
{ {
Decoder::THaCodaFile datafile(dpath.c_str()); Decoder::THaCodaFile datafile(dpath.c_str());
THaEvData *evdata = new Decoder::CodaDecoder(); THaEvData *evdata = new Decoder::CodaDecoder();
...@@ -411,20 +340,6 @@ void write_raw_data(const std::string &dpath, const std::string &opath, int nev, ...@@ -411,20 +340,6 @@ void write_raw_data(const std::string &dpath, const std::string &opath, int nev,
} }
} }
std::ofstream output;
char csv_sep = ',';
if (!wpath.empty()) {
output.open(wpath.c_str());
// csv file header
output << "Event Number";
for (auto &mod : modules) {
for (auto &ch : mod.channels) {
output << csv_sep << ch.name;
}
}
output << std::endl;
}
int count = 0; int count = 0;
bool init_tree = false; bool init_tree = false;
while ((datafile.codaRead() == S_SUCCESS) && (nev-- != 0)) { while ((datafile.codaRead() == S_SUCCESS) && (nev-- != 0)) {
...@@ -482,21 +397,6 @@ void write_raw_data(const std::string &dpath, const std::string &opath, int nev, ...@@ -482,21 +397,6 @@ void write_raw_data(const std::string &dpath, const std::string &opath, int nev,
// fill event // fill event
fill_tree(tree, brdata, init_tree, data_mode); fill_tree(tree, brdata, init_tree, data_mode);
count ++; count ++;
// output csv file
if (wpath.empty()) { continue; }
output << count;
for (auto &mod : modules) {
for (auto &ch : mod.channels) {
auto &bd = brdata[ch.name];
output << csv_sep << "\"[";
for (uint32_t i = 0; i < 192; ++i) {
output << bd.raw[i] << csv_sep;
}
output << "]\"";
}
}
output << std::endl;
} }
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment