diff --git a/include/Fadc250Decoder.h b/include/Fadc250Decoder.h
index aa9479202467353fc1598736ad1a6ebfd20fea98..dae10b857ff1fc5e7f1c6f1759bdd48b92a438b6 100644
--- a/include/Fadc250Decoder.h
+++ b/include/Fadc250Decoder.h
@@ -36,7 +36,7 @@ struct Fadc250Chan
{
bool overflow;
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) {}
};
diff --git a/include/utils.h b/include/utils.h
index c5e3bf92306630b38151dc032931875bf650c27f..7f8fbfad76719d416a5065a7c24d4938320debe7 100644
--- a/include/utils.h
+++ b/include/utils.h
@@ -6,6 +6,7 @@
#include <functional>
#include <unordered_map>
#include "ConfigParser.h"
+#include "ConfigObject.h"
enum ModuleType
@@ -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
diff --git a/scripts/analyze_waveform.py b/scripts/analyze_waveform.py
index e1a509361d439debd84db7c5160c4a31d2986dee..d2e7b2e346715342411dd9fed4d6dead83369a67 100644
--- a/scripts/analyze_waveform.py
+++ b/scripts/analyze_waveform.py
@@ -4,6 +4,7 @@ import json
import numpy as np
import pandas as pd
import argparse
+from collections import OrderedDict
from matplotlib import pyplot as plt
from scipy import signal
from collections.abc import Iterable
@@ -13,14 +14,25 @@ def branch_to_array1d(br, t):
return np.ndarray((len(br), ), dtype=t, buffer=br)
-def get_channels(json_path):
- dbf = open(json_path)
- db = json.load(dbf)
- channels = dict()
- for mod, mod_prop in db.items():
- for ch in mod_prop['channels']:
- channels[ch['name']] = ch['type']
- return channels
+def get_maximum_peaks(tree, chan_pos, ref_pos=0.):
+ # init
+ peaks = np.zeros(len(chan_pos), dtype=np.float32)
+ poses = np.zeros(len(chan_pos), dtype=np.float32)
+
+ for i, (c, twin) in enumerate(chan_pos.items()):
+ cpeaks = branch_to_array1d(tree.__getattr__(c + '_Ppeak'), np.float32)
+ 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')
@@ -41,69 +53,59 @@ dbf = open(args.tconfig)
ch_pos = json.load(dbf)
# 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']
-# tr_time = (17, 23)
-tr = ['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']
-tr_time = (28, 32)
+# scintillator
+sc_ch = ['S2_1', 'S2_2', 'S2_3', 'S2_4', 'S2_5', 'S2_6', 'S2_7', 'S2_8', 'S2_9', 'S2_10', 'S2_11']
+sc_pos = OrderedDict([(ch, (17, 23)) for ch in sc_ch])
+
+# 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')
tree = f.EvTree
-
-trg_ch = np.ndarray(shape=(tree.GetEntries(), ), dtype=object)
-trg_val = np.ndarray(shape=(tree.GetEntries(), 3), dtype='float64')
-ch_val = np.ndarray(shape=(tree.GetEntries(), len(ch_pos)*2), dtype='float64')
-
-for iev in np.arange(0, tree.GetEntries()):
+nev = tree.GetEntries()
+if args.nev > 0 and args.nev <= nev:
+ nev = args.nev
+
+# init buffers
+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)
if iev % 1000 == 0:
print('processed {}'.format(iev), end='\r')
- # triggers
- ev_trg_peaks, ev_trg_poses = [], []
- for c in tr:
- peaks = branch_to_array1d(tree.__getattr__(c + '_Ppeak'), np.float32)
- poses = branch_to_array1d(tree.__getattr__(c + '_Ptime'), np.float32)
- tpeak, tpos = 0, 0
- for peak, pos in zip(peaks, poses):
- if pos <= tr_time[1] and pos >= tr_time[0] and peak > tpeak:
- tpeak = peak
- tpos = pos
- ev_trg_peaks.append(tpeak)
- ev_trg_poses.append(tpos)
-
- # no triggers
- if not len(ev_trg_peaks):
- continue
-
- # 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
+ sc_peaks, sc_poses = get_maximum_peaks(tree, sc_pos)
+ ec_peaks, ec_poses = get_maximum_peaks(tree, ec_pos)
+
+ # save result to buffer, EC
+ iec = np.argmax(ec_peaks)
+ ec_names[iev] = list(ec_pos)[iec]
+ ec_vals[iev] = (ec_peaks[iec], ec_poses[iec], np.sum(ec_peaks > 0), np.sum(ec_peaks))
+ ref = ec_poses[iec]
+
+ # SC
+ isc = np.argmax(sc_peaks)
+ 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))
+
+ # CHER
+ ch_peaks, ch_poses = get_maximum_peaks(tree, ch_pos, ref)
+ ch_vals[iev] = np.concatenate([ch_peaks, ch_poses])
print('processed {}'.format(iev))
-cols = ['trg_peak', 'trg_pos', 'trg_nhits'] \
- + [c + '_peak' for c, _ in ch_pos.items()] \
- + [c + '_pos' for c, _ in ch_pos.items()]
-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))
-result.loc[:, 'trg_ch'] = trg_ch[0:iev+1]
+cols = ['ec_' + col for col in trg_cols] + ['sc_' + col for col in trg_cols] \
+ + [c + '_peak' for c, _ in ch_pos.items()] + [c + '_pos' for c, _ in ch_pos.items()]
+result = pd.DataFrame(index=np.arange(nev), columns=cols,
+ data=np.concatenate((ec_vals, sc_vals, ch_vals), axis=1))
+result.loc[:, 'ec_ch'] = ec_names
+result.loc[:, 'sc_ch'] = sc_names
result.to_csv(args.output)
diff --git a/src/analyze.cpp b/src/analyze.cpp
index db816174b29664fa5e20712ef85db2fc76f579af..2b6c0af9beaef8db2d440ece00b59307a6e577ed 100644
--- a/src/analyze.cpp
+++ b/src/analyze.cpp
@@ -1,3 +1,5 @@
+#include <iostream>
+#include <iomanip>
#include <fstream>
#include "utils.h"
#include "ConfigParser.h"
@@ -6,21 +8,50 @@
#include "THaCodaFile.h"
#include "THaEvData.h"
#include "CodaDecoder.h"
-#include "Fadc250Module.h"
#include "evio.h"
+#include "TSpectrum.h"
#include "TTree.h"
#include "TFile.h"
#include "TH1.h"
+#include "simpleLib.h"
+#include "Fadc250Module.h"
+#include "Fadc250Decoder.h"
+#define FADC_BANK 3
#define PROGRESS_COUNT 1000
-std::vector<Module> read_modules(const std::string &path);
-std::unordered_map<std::string, std::vector<float>> read_timing(const std::string &path);
void write_raw_data(const std::string &dpath, const std::string &opath, int nev, const std::vector<Module> &modules);
-void process_data(const std::string &dpath, const std::string &opath, int nev, const std::vector<Module> &modules,
- std::unordered_map<std::string, std::vector<float>> timing, float nsig, float thres);
+
+// get fadc250 data from the module
+Fadc250Data pack_fadc250_data(Decoder::Fadc250Module* fadc)
+{
+ Fadc250Data res(0, 0, 1);
+ res.mode = fadc->GetMode();
+
+ res.events.resize(1);
+ auto &event = res.events[0];
+
+ for (int i = 0; i < FADC_NCHAN; ++i) {
+ for (int j = 0; j < fadc->GetNumFadcEvents(i); ++j) {
+ auto integral = fadc->GetPulseIntegralData(i, j);
+ auto peak = fadc->GetPulsePeakData(i, j);
+ auto pedestal = fadc->GetPulsePedestalData(i, j);
+ auto time = fadc->GetPulseTimeData(i, j);
+
+ event.channels[i].integral.push_back(integral);
+ event.channels[i].peak.push_back(peak);
+ event.channels[i].pedestal.push_back(pedestal);
+ event.channels[i].time.push_back(time);
+ }
+
+ if (fadc->GetNumSamples(i) > 0) {
+ event.channels[i].raw = fadc->GetPulseSamplesVector(i);
+ }
+ }
+ return res;
+}
int main(int argc, char* argv[])
@@ -29,18 +60,12 @@ int main(int argc, char* argv[])
ConfigOption conf_opt;
conf_opt.AddOpts(ConfigOption::help_message, 'h', "help");
conf_opt.AddOpt(ConfigOption::arg_require, 'n');
- conf_opt.AddOpt(ConfigOption::arg_require, 'c');
- conf_opt.AddOpt(ConfigOption::arg_require, 'f');
conf_opt.AddLongOpt(ConfigOption::arg_require, "config-module", 'm');
- conf_opt.AddLongOpt(ConfigOption::arg_require, "config-timing", 't');
- conf_opt.AddLongOpt(ConfigOption::arg_none, "raw-output", 'r');
+ conf_opt.AddLongOpt(ConfigOption::arg_require, "waveform-output", 'w');
conf_opt.SetDesc("usage: %0 <data_file> <out_file>");
conf_opt.SetDesc('n', "number of events to process (< 0 means all).");
- // conf_opt.SetDesc('c', "timing cut factor, default 3 (3 sigmas).");
- // conf_opt.SetDesc('f', "fired module threshold on ADC value (ped subtracted, default 50.0");
conf_opt.SetDesc('m', "configuration file for modules to be read-in, default \"config/mapmt_module.conf\"");
- // conf_opt.SetDesc('t', "configuration file for signal timing to be read-in, default \"config/mapmt_timing.conf\"");
if (!conf_opt.ParseArgs(argc, argv) || conf_opt.NbofArgs() != 2) {
std::cout << conf_opt.GetInstruction() << std::endl;
@@ -48,28 +73,16 @@ int main(int argc, char* argv[])
}
int nev = -1;
- float nsig = 3.;
- float thres = 50.;
std::string mconf = "config/esb_module.conf";
- std::string tconf = "config/mapmt_timing.conf";
for (auto &opt : conf_opt.GetOptions()) {
switch (opt.mark) {
case 'n':
nev = opt.var.Int();
break;
- case 'c':
- nsig = opt.var.Float();
- break;
- case 'f':
- thres = opt.var.Float();
- break;
case 'm':
mconf = opt.var.String();
break;
- case 't':
- tconf = opt.var.String();
- break;
default :
std::cout << conf_opt.GetInstruction() << std::endl;
return -1;
@@ -77,339 +90,224 @@ int main(int argc, char* argv[])
}
auto modules = read_modules(mconf);
- // auto timing = read_timing(tconf);
write_raw_data(conf_opt.GetArgument(0), conf_opt.GetArgument(1), nev, modules);
}
-std::vector<Module> read_modules(const std::string &path)
+#define MAX_NPEAKS 20
+#define MAX_RAW 300
+struct BranchData
{
-
- // 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())
- });
- }
+ int npul, nraw;
+ float integral[MAX_NPEAKS], peak[MAX_NPEAKS], time[MAX_NPEAKS];
+ int raw[MAX_RAW];
+ float ped_mean, ped_err;
+};
+
+#define BUF_SIZE 1000
+static double buffer[BUF_SIZE], wfbuf[BUF_SIZE], bkbuf[BUF_SIZE];
+void refine_pedestal(const std::vector<uint32_t> &raw, float &ped_mean, float &ped_err)
+{
+ int count = 0;
+ float mean = 0.;
+ for (auto &val : raw) {
+ if (std::abs(val - ped_mean) < 2.0 * ped_err) {
+ buffer[count] = val;
+ mean += val;
+ count ++;
}
- res.emplace_back(mod);
+ }
+ if (count == 0) {
+ return;
}
- // 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;
- }
+ mean /= count;
+ float change = std::abs(ped_mean - mean);
+
+ // no outliers
+ if (change < 0.05 * ped_err) {
+ return;
}
- return res;
+ // recursively refine
+ float err = 0.;
+ for (int i = 0; i < count; ++i) {
+ err += (buffer[i] - mean)*(buffer[i] - mean);
+ }
+ ped_mean = mean;
+ ped_err = std::sqrt(err/count);
+ refine_pedestal(raw, ped_mean, ped_err);
}
-
-std::unordered_map<std::string, std::vector<float>> read_timing(const std::string &path)
+// a function for a simple constant baseline fit
+void find_pedestal(const std::vector<uint32_t> &raw, float &ped_mean, float &ped_err)
{
- std::unordered_map<std::string, std::vector<float>> res;
- ConfigParser parser;
-
- std::cout << "Reading channel timing info from \"" << path << "\"." << std::endl;
- parser.ReadFile(path);
- std::string name;
- float center, sigma;
- while (parser.ParseLine()) {
- parser >> name >> center >> sigma;
- if (res.find(name) != res.end()) {
- std::cout << "WARNING: Duplicated timing information for channel " << name << std::endl;
- }
- res[name] = {center, sigma};
- std::cout << name << ": (" << center << ", " << sigma << ")" << std::endl;
+ ped_mean = 0;
+ ped_err = 0;
+
+ for (auto &val : raw) {
+ ped_mean += val;
}
+ ped_mean /= raw.size();
- return res;
-}
+ for (auto &val : raw) {
+ ped_err += (val - ped_mean)*(val - ped_mean);
+ }
+ ped_err = std::sqrt(ped_err/raw.size());
+ refine_pedestal(raw, ped_mean, ped_err);
+}
-void write_raw_data(const std::string &dpath, const std::string &opath, int nev, const std::vector<Module> &modules)
+// analyze the waveform data and fill the result in a branch data
+void waveform_analysis(const std::vector<uint32_t> &raw, BranchData &res)
{
- Decoder::THaCodaFile datafile(dpath.c_str());
- THaEvData *evdata = new Decoder::CodaDecoder();
-
-#define BUFFER_SIZE 100
- struct BranchData
- {
- int num;
- float integral[BUFFER_SIZE], peak[BUFFER_SIZE], time[BUFFER_SIZE];
- };
+ // no need to analyze
+ if (raw.empty()) {
+ return;
+ }
- auto *hfile = new TFile(opath.c_str(), "RECREATE", "MAPMT test results");
- auto tree = new TTree("EvTree", "Cherenkov Test Events");
- int evtype;
- tree->Branch("Event Type", &evtype, "EventType/I");
- // set up branches and a container to be combined with the branches
- std::unordered_map<std::string, BranchData> brdata;
- for (auto &mod : modules) {
- for (auto &ch : mod.channels) {
- auto n = ch.name;
- if (brdata.find(n) != brdata.end())
- std::cout << "WARNING: Duplicated channel names found! -- " << n << std::endl;
- brdata[n] = BranchData();
- tree->Branch((ch.name + "_Npulse").c_str(), &brdata[n].num, (ch.name + "_N/I").c_str());
- tree->Branch((ch.name + "_Pint").c_str(), &brdata[n].integral[0], (ch.name + "_Pint[" + ch.name + "_N]/F").c_str());
- tree->Branch((ch.name + "_Ppeak").c_str(), &brdata[n].peak[0], (ch.name + "_Ppeak[" + ch.name + "_N]/F").c_str());
- tree->Branch((ch.name + "_Ptime").c_str(), &brdata[n].time[0], (ch.name + "_Ptime[" + ch.name + "_N]/F").c_str());
- }
+ // fill in the raw data
+ res.nraw = raw.size();
+ for (size_t i = 0; i < raw.size(); ++i) {
+ res.raw[i] = raw[i];
}
- int count = 0;
- while ((datafile.codaRead() == S_SUCCESS) && (nev-- != 0)) {
- // read-in data
- evdata->LoadEvent(datafile.getEvBuffer());
- evtype = evdata->GetEvType();
+ // get pedestals
+ find_pedestal(raw, res.ped_mean, res.ped_err);
- if((++count % PROGRESS_COUNT) == 0) {
- std::cout << "Processed events - " << count << "\r" << std::flush;
- }
+ // fill in spectrum buffer
+ for (size_t i = 0; i < raw.size(); ++i) {
+ wfbuf[i] = raw[i] - res.ped_mean;
+ }
- // clear data buffer
- for (auto &br : brdata) {
- br.second.num = 0;
+ // find peaks
+ TSpectrum s;
+ s.SetResolution(0.5);
+ int npeaks = s.SearchHighRes(wfbuf, bkbuf, res.nraw, 3.0, 20, false, 5, true, 3);
+
+ // fill branch data
+ double *pos = s.GetPositionX();
+ res.npul = 0;
+ for (int i = 0; i < npeaks; ++i) {
+ int j = pos[i];
+ if (wfbuf[j] < 5.*res.ped_err) { continue; }
+ res.time[res.npul] = j;
+ res.peak[res.npul] = wfbuf[j];
+ /*
+ res.integral[res.npul] = wfbuf[j];
+ j = pos[i] - 1;
+ while ( (j > 0) && (wfbuf[j] - wfbuf[j - 1])*wfbuf[j] > 0. ) {
+ res.integral[res.npul] += wfbuf[j--];
+ }
+ j = pos[i] + 1;
+ while ( (j < res.nraw - 1) && (wfbuf[j] - wfbuf[j + 1])*wfbuf[j] > 0. ) {
+ res.integral[res.npul] += wfbuf[j++];
}
+ */
+ res.npul ++;
+ }
+}
- // read data into event buffer
- for (auto &mod : modules) {
- auto *fadc = dynamic_cast<Decoder::Fadc250Module*>(evdata->GetModule(mod.crate, mod.slot));
- if ((evdata->GetNumRaw(mod.crate, mod.slot) == 0) || (fadc == nullptr))
- std::cout << " lll " << evdata->GetNumRaw(mod.crate, mod.slot) << std::endl;
- continue;
-
- for (auto &ch : mod.channels) {
- auto pulses = get_pulses(fadc, ch.id, mod.type);
- auto &data = brdata[ch.name];
- data.num = pulses.size();
- for (int i = 0; i < data.num; ++i) {
- data.integral[i] = pulses[i].integral;
- data.peak[i] = pulses[i].peak;
- std::cout << " i " << i << " int " << pulses[i].integral << std::endl;
- data.time[i] = pulses[i].time;
- }
+
+// fill decoded FADC250 data into the container (branch data)
+void fill_branch(const Fadc250Data &slot_data, const Module &mod, std::unordered_map<std::string, BranchData> &brdata)
+{
+ for (auto &event : slot_data.events) {
+ for (auto &ch : mod.channels) {
+ auto &channel = event.channels[ch.id];
+ auto &bd = brdata[ch.name];
+ bd.npul = channel.integral.size();
+ for (uint32_t i = 0; i < channel.integral.size(); ++i) {
+ bd.integral[i] = channel.integral[i] - 23./4.*channel.pedestal[i];
+ bd.peak[i] = channel.peak[i] - 1./4.*channel.pedestal[i];
+ bd.time[i] = channel.time[i] * 0.0625;
}
+ waveform_analysis(channel.raw, bd);
}
+ }
+}
- // take reference time difference
- if (brdata["Ref"].num == 0)
- continue;
- auto ref_time = brdata["Ref"].time[0];
- for (auto &br: brdata) {
- for (int i = 0; i < br.second.num; ++i) {
- br.second.time[i] -= ref_time;
+// fill branch data into the root tree
+void fill_tree(TTree *tree, std::unordered_map<std::string, BranchData> &brdata, bool &init, int mode)
+{
+ if ((mode != 1) && (mode != 3)) {
+ std::cout << "Warning: unsupported mode " << mode << ", data won't be recorded." << std::endl;
+ return;
+ }
+
+ // initialize
+ if (!init) {
+ for (auto &it : brdata) {
+ auto n = it.first;
+ tree->Branch((n + "_Npulse").c_str(), &brdata[n].npul, (n + "_N/I").c_str());
+ tree->Branch((n + "_Pint").c_str(), &brdata[n].integral[0], (n + "_Pint[" + n + "_N]/F").c_str());
+ tree->Branch((n + "_Ptime").c_str(), &brdata[n].time[0], (n + "_Ptime[" + n + "_N]/F").c_str());
+ // raw waveform provides more information
+ if (mode == 1) {
+ tree->Branch((n + "_Ppeak").c_str(), &brdata[n].peak[0], (n + "_Ppeak[" + n + "_N]/F").c_str());
+ tree->Branch((n + "_Nraw").c_str(), &brdata[n].nraw, (n + "_Nraw/I").c_str());
+ tree->Branch((n + "_raw").c_str(), &brdata[n].raw[0], (n + "_raw[" + n + "_Nraw]/I").c_str());
+ tree->Branch((n + "_ped_mean").c_str(), &brdata[n].ped_mean, (n + "_ped_mean/F").c_str());
+ tree->Branch((n + "_ped_err").c_str(), &brdata[n].ped_err, (n + "_ped_err/F").c_str());
}
}
-
- tree->Fill();
+ init = true;
+ std::cout << "Initialized root file for mode " << mode << " data." << std::endl;
}
- std::cout << "Processed events - " << count << std::endl;
- hfile->Write();
- hfile->Close();
+ tree->Fill();
}
-void process_data(const std::string &dpath, const std::string &opath, int nev, const std::vector<Module> &modules,
- std::unordered_map<std::string, std::vector<float>> timing, float nsig, float thres)
+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());
THaEvData *evdata = new Decoder::CodaDecoder();
- auto *hfile = new TFile(opath.c_str(), "RECREATE", "MAPMT test results");
-
- std::unordered_map<std::string, std::vector<PulseData>> event;
- std::unordered_map<std::string, std::vector<TH1*>> hists;
- std::vector<std::string> signals;
- std::unordered_map<std::string, PulseData> pulses;
- // initialize histograms and timing
+ auto *hfile = new TFile(opath.c_str(), "RECREATE", "MAPMT test results");
+ auto tree = new TTree("EvTree", "Cherenkov Test Events");
+ int evtype;
+ tree->Branch("Event Type", &evtype, "EventType/I");
+ // set up branches and a container to be combined with the branches
+ std::unordered_map<std::string, BranchData> brdata;
for (auto &mod : modules) {
for (auto &ch : mod.channels) {
- if (ch.type != kSignal)
- continue;
- signals.push_back(ch.name);
- pulses[ch.name] = PulseData();
- if (hists.find(ch.name) != hists.end())
- std::cout << "WARNING: Duplicated channel names found! -- " << ch.name << std::endl;
-
- hists[ch.name] = {
- new TH1F((ch.name + "_tdiff").c_str(), (ch.name + "_tdiff").c_str(), 400, 0, 0),
- new TH1F((ch.name + "_tdiff_cut").c_str(), (ch.name + "_tdiff_cut").c_str(), 100, 0, 0),
- new TH1F((ch.name + "_Ppk").c_str(), (ch.name + "_Ppk").c_str(), 100, 0, 0)
- };
-
- if (timing.find(ch.name) == timing.end()) {
- std::cout << "WARNING: Cannot find timing information for channel " << ch.name
- << ", using default value (20, 1)" << std::endl;
- timing[ch.name] = {20., 1.};
+ auto n = ch.name;
+ if (brdata.find(n) != brdata.end()) {
+ std::cout << "WARNING: Duplicated channel names found! -- " << n << std::endl;
}
+ brdata[n] = BranchData();
}
}
- auto trghist = new TH1F("Trigger Timing", "CaloSum Trigger", 100, -200, 200);
- auto sumhist = new TH1F("Sum of Cerenkov Signals", "Signal Sum", 200, 0, 0);
- std::map<int, TH1*> sumhists {
- {2, new TH1F("Sum of Cerenkov Signals Nhit = 2", "Signal Sum Nhit = 2", 200, 0, 0)},
- {3, new TH1F("Sum of Cerenkov Signals Nhit = 3", "Signal Sum Nhit = 3", 200, 0, 0)},
- {4, new TH1F("Sum of Cerenkov Signals Nhit = 4", "Signal Sum Nhit = 4", 200, 0, 0)},
- {5, new TH1F("Sum of Cerenkov Signals Nhit = 5", "Signal Sum Nhit = 5", 200, 0, 0)},
- {6, new TH1F("Sum of Cerenkov Signals Nhit = 6", "Signal Sum Nhit = 6", 200, 0, 0)},
- {7, new TH1F("Sum of Cerenkov Signals Nhit = 7", "Signal Sum Nhit = 7", 200, 0, 0)},
- {8, new TH1F("Sum of Cerenkov Signals Nhit = 8", "Signal Sum Nhit = 8", 200, 0, 0)},
- };
- int sum_thres = 6;
- auto sumhist2 = new TH1F(Form("Sum of Cerenkov Signals (Mhit > %i)", sum_thres),
- Form("Signal Sum (Nhit > %i)", sum_thres),
- 200, 0, 0);
- auto cnthist = new TH1I("Number of Fired Modules", "Fire Count", 16, 1, 16);
-
- auto ocomp = ConfigParser::decompose_path(opath);
- ocomp.ext = "dat";
- std::ofstream output(ConfigParser::compose_path(ocomp));
- char csv_sep = ',';
-
- // csv file header
- output << "Event Number" << csv_sep << "Event Type" << csv_sep << "Fired Counts" << csv_sep <<"Signal Sum" << csv_sep << "LED";
- for (auto &sig : signals) {
- output << csv_sep << sig;
- }
- output << std::endl;
int count = 0;
+ bool init_tree = false;
while ((datafile.codaRead() == S_SUCCESS) && (nev-- != 0)) {
// read-in data
evdata->LoadEvent(datafile.getEvBuffer());
- if((++count % PROGRESS_COUNT) == 0) {
- std::cout << "Processed events - " << count << "\r" << std::flush;
- }
+ evtype = evdata->GetEvType();
- for (auto &pul : pulses) {
- pul.second = PulseData();
+ if((count % PROGRESS_COUNT) == 0) {
+ std::cout << "Processed events - " << count << "\r" << std::flush;
}
- // read data into event buffer
- bool empty_slot = false;
+ int data_mode = -1;
for (auto &mod : modules) {
auto *fadc = dynamic_cast<Decoder::Fadc250Module*>(evdata->GetModule(mod.crate, mod.slot));
- if ((evdata->GetNumRaw(mod.crate, mod.slot) == 0) || (fadc == nullptr))
- empty_slot = true;
-
- for (auto &ch : mod.channels) {
- if (empty_slot) {
- event[ch.name] = std::vector<PulseData>();
- } else {
- event[ch.name] = get_pulses(fadc, ch.id, mod.type);
- }
+ // no data
+ if ((evdata->GetNumRaw(mod.crate, mod.slot) == 0) || (fadc == nullptr)) { continue; }
+ auto slot_data = pack_fadc250_data(fadc);
+ // check data mode
+ if (slot_data.GetMode() > 0 && data_mode < 0) {
+ data_mode = slot_data.GetMode();
}
+ // fill branch data
+ fill_branch(slot_data, mod, brdata);
}
-
- // further process event
- // no reference time or triggered by LED
- if ((event["Ref"].size() == 0) || (event["LED"].size() > 0))
- continue;
- auto ref_time = event["Ref"].front().time;
- // Calosum timing
- float trg_pulse = 0.;
- float trg_time = 0.;
- for (auto &cal : event["CaloSum"]) {
- if (trg_pulse > cal.integral)
- continue;
- trg_time = cal.time;
- trg_pulse = cal.integral;
- }
-
- if (trg_time == 0.)
- continue;
- trghist->Fill(trg_time);
-
- int fire = 0;
- float sum = 0;
- for (auto &it : hists) {
- auto hist = it.second;
- float integral = thres;
- for (auto &ev : event[it.first]) {
- float tdiff = ev.time - trg_time;
- // tdiff
- hist[0]->Fill(tdiff);
- if ((std::abs(tdiff - timing[it.first][0]) < nsig*timing[it.first][1]) &&
- (ev.integral > integral)) {
- pulses[it.first] = ev;
- integral = ev.integral;
- // tdiff_cut
- hist[1]->Fill(tdiff);
- }
- }
- if (integral > thres) {
- // pint
- hist[2]->Fill(integral);
- fire ++;
- sum += integral;
- }
- }
-
- if (fire > 0) {
- sumhist->Fill(sum);
- cnthist->Fill(fire);
- auto it = sumhists.find(fire);
- if (it != sumhists.end()) {
- it->second->Fill(sum);
- } else if (fire > sum_thres) {
- sumhist2->Fill(sum);
- }
- }
-
- if (fire >= 3) {
- output << count << csv_sep << evdata->GetEvType() << csv_sep << fire << csv_sep << sum << csv_sep << 0;
- for (auto &sig : signals) {
- output << csv_sep << pulses[sig].integral;
- }
- output << std::endl;
- }
-
+ // fill event
+ fill_tree(tree, brdata, init_tree, data_mode);
+ count ++;
}
std::cout << "Processed events - " << count << std::endl;
@@ -417,4 +315,3 @@ void process_data(const std::string &dpath, const std::string &opath, int nev, c
hfile->Close();
}
-
diff --git a/src/esb_analyze.cpp b/src/esb_analyze.cpp
index 5a34c6eda20791e0dd493c3132f0f359b25994d1..146b3bf759967701749cfda2a98f6b74cbe1697f 100644
--- a/src/esb_analyze.cpp
+++ b/src/esb_analyze.cpp
@@ -21,8 +21,7 @@
#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, const std::string &wpath);
+void write_raw_data(const std::string &dpath, const std::string &opath, int nev, const std::vector<Module> &modules);
bool parseEvent(const uint32_t *buf, bool verbose);
uint32_t parseBlock(const uint32_t *buf);
@@ -135,7 +134,6 @@ int main(int argc, char* argv[])
int nev = -1;
std::string mconf = "config/esb_module.conf";
- std::string wpath = "";
for (auto &opt : conf_opt.GetOptions()) {
switch (opt.mark) {
@@ -145,9 +143,6 @@ int main(int argc, char* argv[])
case 'm':
mconf = opt.var.String();
break;
- case 'w':
- wpath = opt.var.String();
- break;
default :
std::cout << conf_opt.GetInstruction() << std::endl;
return -1;
@@ -156,73 +151,7 @@ int main(int argc, char* argv[])
auto modules = read_modules(mconf);
- write_raw_data(conf_opt.GetArgument(0), conf_opt.GetArgument(1), nev, modules, wpath);
-}
-
-
-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;
+ write_raw_data(conf_opt.GetArgument(0), conf_opt.GetArgument(1), nev, modules);
}
@@ -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());
THaEvData *evdata = new Decoder::CodaDecoder();
@@ -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;
bool init_tree = false;
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,
// fill event
fill_tree(tree, brdata, init_tree, data_mode);
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;
}
}