Adjust simulation analysis script

run.mac

@@ -4,7 +4,7 @@
 # beam starting point (x, y, z) in m
 /cher/gun/pos 0 0 -1
 # beam distribution (dx, dy, dz) in mm
-/cher/gun/sigma 0 0 0
+/cher/gun/sigma 80 80 0
 
 # beam polar angle in degree
 /cher/gun/theta 0

scripts/sim_analysis.py

@@ -4,6 +4,7 @@
 import pandas as pd
 import numpy as np
 import os
+import errno
 from matplotlib import pyplot as plt
 from scipy.optimize import curve_fit
 import csv_converter
@@ -32,14 +33,13 @@ parser = argparse.ArgumentParser('CSV converter')
 
 parser.add_argument('infile', help='simulation output in csv or root format')
 parser.add_argument('-o', '--output-dir', dest='outdir', default='.', help='path to the output directory')
-parser.add_argument('-n', dest='nh', type=int, default=6, help='number of hits to fit the signal sum')
 parser.add_argument('-t', dest='thres', type=float, default=50., help='threshold for fired PMT channel')
 parser.add_argument('-s', dest='spe', type=float, default=80., help='single photo electron peak height')
-parser.add_argument('--prefix', dest='prefix', default='sim', help='prefix for output files')
+parser.add_argument('--nhits-min', dest='nhmin', type=int, default=4, help='minimum number of hits of the signal sum')
+parser.add_argument('--nhits-max', dest='nhmax', type=int, default=7, help='maximum number of hits of the signal sum')
 
 args = parser.parse_args()
 
-
 df = pd.DataFrame()
 ftype = args.infile.split('.')[-1]
 if ftype == 'csv':
@@ -50,6 +50,13 @@ else:
     print('Unsupported file type {}'.format(ftype))
 
 
+try:
+    os.makedirs(args.outdir)
+except OSError as e:
+    if e.errno != errno.EEXIST:
+        raise
+
+
 channels = ['Cer_{}'.format(i + 1) for i in np.arange(16)]
 ncols, nrows = auto_split(len(channels))
 figsize = (ncols*6, nrows*4)
@@ -57,7 +64,7 @@ figsize = (ncols*6, nrows*4)
 cycle = plt.rcParams['axes.prop_cycle'].by_key()['color']
 bprops = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
 
-# peaks
+# signals (normalized by spe)
 prange = (0, 30)
 dbin = 0.2
 gridspec = {'left': 0.08, 'bottom': 0.08, 'right': 0.95, 'top': 0.95, 'hspace': 0., 'wspace': 0.}
@@ -81,11 +88,11 @@ for i, ax in enumerate(axs.flat):
 fig.legend(lines, labels, loc='upper center', ncol=5, fontsize=22)
 fig.text(0.5, 0.03, 'N.p.e', ha='center', fontsize=20)
 fig.text(0.03, 0.5, 'Counts', va='center', rotation=90, fontsize=20)
-fig.savefig(os.path.join(args.outdir, '{}_quad_peaks.png'.format(args.prefix)))
+fig.savefig(os.path.join(args.outdir, 'sim_quadsig.png'))
 
 # number of photo-electrons
-prange = (0, 10)
-dbin = 0.2
+prange = (0.5, 10.5)
+dbin = 1
 gridspec = {'left': 0.08, 'bottom': 0.08, 'right': 0.95, 'top': 0.95, 'hspace': 0., 'wspace': 0.}
 fig, axs = plt.subplots(nrows, ncols, figsize=figsize, sharex='all', sharey='all', gridspec_kw=gridspec)
 lines, labels = [], []
@@ -96,30 +103,51 @@ for i, ax in enumerate(axs.flat):
     ch = channels[i]
     ax.text(0.5, 0.90, ch, transform=ax.transAxes, fontsize=22, verticalalignment='top', bbox=bprops)
     mask = (df[ch + '.npe'] > 0)
-    data = df.loc[mask, ch + '.npe'].values/args.spe
+    if not np.sum(mask):
+        continue
+    data = df.loc[mask, ch + '.npe'].values
     vals, bins = np.histogram(data, np.arange(*prange, step=dbin))
     didx = (bins[1:] + bins[:-1])/2.
     ax.bar(didx, vals, width=dbin, label='Simulation')
-    ax.set_yscale('log')
     ax.tick_params(labelsize=18)
     lines, labels = ax.get_legend_handles_labels()
 
 fig.legend(lines, labels, loc='upper center', ncol=5, fontsize=22)
 fig.text(0.5, 0.03, 'N.p.e', ha='center', fontsize=20)
 fig.text(0.03, 0.5, 'Counts', va='center', rotation=90, fontsize=20)
-fig.savefig(os.path.join(args.outdir, '{}_quad_npe.png'.format(args.prefix)))
+fig.savefig(os.path.join(args.outdir, 'sim_quadnpe.png'))
 
 
 df.loc[:, 'nhits'] = (df[[ch + '.peak' for ch in channels]] > args.thres).sum(axis=1)
 df.loc[:, 'signal_sum'] = df[[ch + '.peak' for ch in channels]].sum(axis=1)/args.spe
 
+for nh in np.arange(args.nhmin, args.nhmax + 1):
+    sbins = np.linspace(0, 40, 100)
+    fig, ax = plt.subplots(1, 1, figsize=(16, 12))
+    ax.set_xlabel('N.p.e', fontsize=18)
+    ax.set_ylabel('Counts', fontsize=18)
+    ax.set_title('Signal Sum Distribution', fontsize=18)
+    hvals = ax.hist(df.loc[df['nhits'] == nh, 'signal_sum'].values, sbins, label=r'Simulation (Nhits = {})'.format(nh),
+                    color=cycle[0], edgecolor='black', alpha=0.75)
+    idx = (sbins[1:] + sbins[:-1])/2.
+    pars, _ = curve_fit(dgaus, idx, hvals[0], p0=[500, 15, 4, 30, 8])
+    ax.plot(idx, dgaus(idx, *pars), 'k-', lw=2, label='Fit, N.p.e = {:.2f}'.format(np.square(pars[1]/pars[2])))
+    ax.plot(idx, gaus(idx, *pars[:3]), 'r--', lw=2, label='Regular Signal Fit')
+    ax.plot(idx, gaus(idx, pars[3], pars[1]*2, pars[4]), 'g--', lw=2, label='Pair Signal Fit')
+    ax.tick_params('both', labelsize=16)
+    ax.legend(fontsize=18)
+    # ax.legend(fontsize=18, title='N.p.e = {:.2f}'.format(np.square(pars[1]/pars[2])), title_fontsize=18)
+    fig.tight_layout()
+    fig.savefig(os.path.join(args.outdir, 'sim_npe{}.png'.format(nh)))
+
+
 sbins = np.linspace(0, 40, 100)
 fig, ax = plt.subplots(1, 1, figsize=(16, 12))
 ax.set_xlabel('N.p.e', fontsize=18)
 ax.set_ylabel('Counts', fontsize=18)
 ax.set_title('Signal Sum Distribution', fontsize=18)
-hvals = ax.hist(df.loc[df['nhits'] == args.nh, 'signal_sum'].values, sbins, label=r'Simulation (Nhits = {})'.format(args.nh),
-                color=cycle[0], edgecolor='black', alpha=0.75)
+hvals = ax.hist(df.loc[df['nhits'].isin(np.arange(args.nhmin, args.nhmax + 1)), 'signal_sum'].values, sbins,
+                label=r'Simulation'.format(nh), color=cycle[0], edgecolor='black', alpha=0.75)
 idx = (sbins[1:] + sbins[:-1])/2.
 pars, _ = curve_fit(dgaus, idx, hvals[0], p0=[500, 15, 4, 30, 8])
 ax.plot(idx, dgaus(idx, *pars), 'k-', lw=2, label='Fit, N.p.e = {:.2f}'.format(np.square(pars[1]/pars[2])))
@@ -129,5 +157,5 @@ ax.tick_params('both', labelsize=16)
 ax.legend(fontsize=18)
 # ax.legend(fontsize=18, title='N.p.e = {:.2f}'.format(np.square(pars[1]/pars[2])), title_fontsize=18)
 fig.tight_layout()
-fig.savefig(os.path.join(args.outdir, '{}_sum_fits.png'.format(args.prefix)))
+fig.savefig(os.path.join(args.outdir, 'sim_npe.png'))
 

src/CherProtoPixelSD.cc

@@ -102,7 +102,7 @@ void CherProtoPixelSD::regAnalysis()
 
 void CherProtoPixelSD::Initialize(G4HCofThisEvent *HCE)
 {
-    hitsCollection = new CherProtoHitCollection(GetName(),collectionName[0]);
+    hitsCollection = new CherProtoHitCollection(GetName(), collectionName[0]);
 
     static G4int HCID = -1;
     if (HCID < 0) { HCID = GetCollectionID(0); }

src/EventAction.cc

@@ -73,4 +73,3 @@ void EventAction::EndOfEventAction(const G4Event * /*evt*/)
     // place holder
 }
 
-