Adds more sophisticated plotting machinery (from filval), and adds new

plots for examining second hit resolutions

analysis/plots.py

@@ -1,163 +1,160 @@
 #!/usr/bin/env python
-
-import os
 import sys
-import matplotlib as mpl
-import matplotlib.pyplot as plt
 
-fv_path = os.path.dirname(os.path.abspath(__file__))+"/../filval/python"
-sys.path.append(fv_path)
-
-from result_set import ResultSet
-from plotter import histogram, plot_histogram, plot_histogram2d, histogram_slice
-mpl.rc('text', usetex=True)
-mpl.rc('savefig', dpi=180)
-mpl.rc('savefig', transparent=False)
+import matplotlib.pyplot as plt
 
+from filval.result_set import ResultSet
+from filval.histogram_utils import hist, hist2d, hist_slice
+from filval.plotter import make_plot, plot_registry, hist_plot, hist2d_plot
 
 
+@make_plot(scale=0.85)
 def first_hits_v_eta(rs):
-    scale = 0.85
-    plt.figure(figsize=(scale*10, scale*10))
+    r'''
+    Plots of $\Delta \phi$, $\Delta z$, and $\Delta r$ vs $\eta$ between RecHit and SimHit for
+    the first matched hit in the seed.
+    '''
 
     plt.subplot(321)
-    plot_histogram2d(rs.dphi_v_eta_first_hits_in_B1,
-                     ylabel="$\\Delta \\phi_1$(rad)",
-                     title="BPIX - Layer 1")
+    hist2d_plot(hist2d(rs.dphi_v_eta_first_hits_in_B1),
+                ylabel=r"$\Delta \phi_1$(rad)",
+                title="BPIX - Layer 1")
 
     plt.subplot(322)
-    plot_histogram2d(rs.dphi_v_eta_first_hits_in_B2,
-                     title="BPIX - Layer 2")
+    hist2d_plot(hist2d(rs.dphi_v_eta_first_hits_in_B2),
+                title="BPIX - Layer 2")
 
     plt.subplot(323)
-    plot_histogram2d(rs.dz_v_eta_first_hits_in_B1,
-                     ylabel="$\\Delta \\textrm{z}_1$(cm)")
+    hist2d_plot(hist2d(rs.dz_v_eta_first_hits_in_B1),
+                ylabel=r"$\Delta \textrm{z}_1$(cm)",
+                )
 
     plt.subplot(324)
-    plot_histogram2d(rs.dz_v_eta_first_hits_in_B2)
+    hist2d_plot(hist2d(rs.dz_v_eta_first_hits_in_B2))
 
     plt.subplot(325)
-    plot_histogram2d(rs.dr_v_eta_first_hits_in_B1,
-                     ylabel="$\\Delta r_1$(cm)",
-                     xlabel="$\\eta$")
+    hist2d_plot(hist2d(rs.dr_v_eta_first_hits_in_B1),
+                ylabel=r"$\Delta r_1$(cm)",
+                xlabel=r"$\eta$"
+                )
 
     plt.subplot(326)
-    plot_histogram2d(rs.dr_v_eta_first_hits_in_B2,
-                     xlabel="$\\eta$")
+    hist2d_plot(hist2d(rs.dr_v_eta_first_hits_in_B2),
+                xlabel=r"$\eta$"
+                )
 
-    plt.tight_layout()
-    plt.savefig("figures/first_hits_v_eta.png")
 
-
-plt.clf()
+@make_plot(scale=0.85)
 def second_hits_v_eta(rs):
     plt.subplot(321)
-    plot_histogram2d(rs.dphi_v_eta_second_hits_in_B2,
-                     ylabel="$\\Delta \\phi_2$(rad)",
-                     title="BPIX - Layer 2")
+    hist2d_plot(hist2d(rs.dphi_v_eta_second_hits_in_B2),
+                ylabel=r"$\Delta \phi_2$(rad)",
+                title="BPIX - Layer 2")
 
     plt.subplot(322)
-    plot_histogram2d(rs.dphi_v_eta_second_hits_in_B3,
-                     title="BPIX - Layer 3")
+    hist2d_plot(hist2d(rs.dphi_v_eta_second_hits_in_B3),
+                title="BPIX - Layer 3")
 
     plt.subplot(323)
-    plot_histogram2d(rs.dz_v_eta_second_hits_in_B2,
-                     ylabel="$\\Delta \\textrm{z}_2$(cm)")
+    hist2d_plot(hist2d(rs.dz_v_eta_second_hits_in_B2),
+                ylabel=r"$\Delta \textrm{z}_2$(cm)")
 
     plt.subplot(324)
-    plot_histogram2d(rs.dz_v_eta_second_hits_in_B3)
+    hist2d_plot(hist2d(rs.dz_v_eta_second_hits_in_B3))
 
     plt.subplot(325)
-    plot_histogram2d(rs.dr_v_eta_second_hits_in_B2,
-                     ylabel="$\\Delta r_2$(cm)",
-                     xlabel="$\\eta$")
+    hist2d_plot(hist2d(rs.dr_v_eta_second_hits_in_B2),
+                ylabel=r"$\Delta r_2$(cm)",
+                xlabel=r"$\eta$")
 
     plt.subplot(326)
-    plot_histogram2d(rs.dr_v_eta_second_hits_in_B3,
-                     xlabel="$\\eta$")
-
-    plt.tight_layout()
-    plt.savefig("figures/second_hits_v_eta.png")
+    hist2d_plot(hist2d(rs.dr_v_eta_second_hits_in_B3),
+                xlabel=r"$\eta$")
 
 
-plt.clf()
+@make_plot(scale=0.85)
 def first_hits(rs):
     plt.subplot(321)
-    plot_histogram(rs.dphi_v_eta_first_hits_in_B1.ProjectionY(),
-                   include_errors=True, xlabel="$\\Delta \\phi_1$(rad)",
-                   title="BPIX - Layer 1")
+    hist_plot(hist(rs.dphi_v_eta_first_hits_in_B1.ProjectionY()),
+              include_errors=True, xlabel=r"$\Delta \phi_1$(rad)",
+              title="BPIX - Layer 1")
 
     plt.subplot(322)
-    plot_histogram(rs.dphi_v_eta_first_hits_in_B2.ProjectionY(), include_errors=True,
-                   xlabel="$\\Delta \\phi_1$(rad)",
-                   title="BPIX - Layer 2")
+    hist_plot(hist(rs.dphi_v_eta_first_hits_in_B2.ProjectionY()),
+              include_errors=True,
+              xlabel=r"$\Delta \phi_1$(rad)",
+              title="BPIX - Layer 2")
 
     plt.subplot(323)
-    plot_histogram(rs.dz_v_eta_first_hits_in_B1.ProjectionY(), include_errors=True,
-                   xlabel="$\\Delta z_1$(cm)")
+    hist_plot(hist(rs.dz_v_eta_first_hits_in_B1.ProjectionY()),
+              include_errors=True,
+              xlabel=r"$\Delta z_1$(cm)")
 
     plt.subplot(324)
-    plot_histogram(rs.dz_v_eta_first_hits_in_B2.ProjectionY(), include_errors=True,
-                   xlabel="$\\Delta z_1$(cm)")
+    hist_plot(hist(rs.dz_v_eta_first_hits_in_B2.ProjectionY()),
+              include_errors=True,
+              xlabel=r"$\Delta z_1$(cm)")
 
     plt.subplot(325)
-    plot_histogram(rs.dr_v_eta_first_hits_in_B1.ProjectionY(), include_errors=True,
-                   xlabel="$\\Delta r_1$(cm)")
+    hist_plot(hist(rs.dr_v_eta_first_hits_in_B1.ProjectionY()),
+              include_errors=True,
+              xlabel=r"$\Delta r_1$(cm)")
 
     plt.subplot(326)
-    plot_histogram(rs.dr_v_eta_first_hits_in_B2.ProjectionY(), include_errors=True,
-                   xlabel="$\\Delta r_1$(cm)")
-
-    plt.tight_layout()
-    plt.savefig("figures/first_hits.png")
+    hist_plot(hist(rs.dr_v_eta_first_hits_in_B2.ProjectionY()),
+              include_errors=True,
+              xlabel="r$\Delta r_1$(cm)")
 
 
-plt.clf()
+@make_plot(scale=0.85)
 def second_hits(rs):
     plt.subplot(321)
-    plot_histogram(rs.dphi_v_eta_second_hits_in_B2.ProjectionY(), include_errors=True,
-                   xlabel="$\\Delta \\phi_2$(rad)",
-                   title="BPIX - Layer 2")
+    hist_plot(hist(rs.dphi_v_eta_second_hits_in_B2.ProjectionY()),
+              include_errors=True,
+              xlabel=r"$\Delta \phi_2$(rad)",
+              title="BPIX - Layer 2")
 
     plt.subplot(322)
-    plot_histogram(rs.dphi_v_eta_second_hits_in_B3.ProjectionY(), include_errors=True,
-                   xlabel="$\\Delta \\phi_2$(rad)",
-                   title="BPIX - Layer 3")
+    hist_plot(hist(rs.dphi_v_eta_second_hits_in_B3.ProjectionY()),
+              include_errors=True,
+              xlabel=r"$\Delta \phi_2$(rad)",
+              title="BPIX - Layer 3")
 
     plt.subplot(323)
-    plot_histogram(rs.dz_v_eta_second_hits_in_B2.ProjectionY(), include_errors=True,
-                   xlabel="$\\Delta z_2$(cm)")
+    hist_plot(hist(rs.dz_v_eta_second_hits_in_B2.ProjectionY()),
+              include_errors=True,
+              xlabel=r"$\Delta z_2$(cm)")
 
     plt.subplot(324)
-    plot_histogram(rs.dz_v_eta_second_hits_in_B3.ProjectionY(), include_errors=True,
-                   xlabel="$\\Delta z_2$(cm)")
+    hist_plot(hist(rs.dz_v_eta_second_hits_in_B3.ProjectionY()),
+              include_errors=True,
+              xlabel=r"$\Delta z_2$(cm)")
 
     plt.subplot(325)
-    plot_histogram(rs.dr_v_eta_second_hits_in_B2.ProjectionY(), include_errors=True,
-                   xlabel="$\\Delta r_2$(cm)")
+    hist_plot(hist(rs.dr_v_eta_second_hits_in_B2.ProjectionY()),
+              include_errors=True,
+              xlabel=r"$\Delta r_2$(cm)")
 
     plt.subplot(326)
-    plot_histogram(rs.dr_v_eta_second_hits_in_B3.ProjectionY(), include_errors=True,
-                   xlabel="$\\Delta r_2$(cm)")
+    hist_plot(hist(rs.dr_v_eta_second_hits_in_B3.ProjectionY()),
+              include_errors=True,
+              xlabel=r"$\Delta r_2$(cm)")
 
-    plt.tight_layout()
-    plt.savefig("figures/second_hits.png")
 
-
-plt.clf()
+@make_plot(scale=0.85)
 def delta_phi_z_v_ladder(rs):
     def even_odd_plot(even, odd, var, scale, unit, **kwargs):
         even_dist = even.ProjectionY()
         odd_dist = odd.ProjectionY()
-        plot_histogram(even_dist, include_errors=True, label="even ladders")
-        plot_histogram(odd_dist, include_errors=True, color='r', label="odd ladders",
-                       **kwargs)
+        hist_plot(hist(even_dist), include_errors=True, label="even ladders")
+        hist_plot(hist(odd_dist), include_errors=True, color='r', label="odd ladders",
+                  **kwargs)
 
         even_mean = even_dist.GetMean()*scale
         odd_mean = odd_dist.GetMean()*scale
         axes = plt.gca()
         txt = r"$ \hat{{\Delta {0} }}_{{1,\textrm{{ {1} }} }}={2:4.2g}${3}"
-        axes.text(0.05, .7, txt.format(var, "odd",  odd_mean, unit), transform=axes.transAxes)
+        axes.text(0.05, .7, txt.format(var, "odd", odd_mean, unit), transform=axes.transAxes)
         axes.text(0.05, .6, txt.format(var, "even", even_mean, unit), transform=axes.transAxes)
     plt.subplot(221)
     even_odd_plot(rs.dphi_v_eta_first_hits_in_B1_even_ladder, rs.dphi_v_eta_first_hits_in_B1_odd_ladder,
@@ -176,158 +173,147 @@ def delta_phi_z_v_ladder(rs):
     even_odd_plot(rs.dz_v_eta_first_hits_in_B2_even_ladder, rs.dz_v_eta_first_hits_in_B2_odd_ladder,
                   "z", 10**4, "um", xlabel=r"$\Delta z_1$(cm)")
 
-    plt.tight_layout()
-    plt.savefig("figures/delta_phi_z_v_ladder.png")
 
+@make_plot(scale=0.85)
 def sc_extrapolation_first(rs):
-    # Raphael's plots
+    ''' Raphael's plots '''
     norm = 1
     errors = True
+
     def preproc(h):
-        return histogram_slice(histogram(h.ProjectionY()), (-0.07, 0.07))
+        return hist_slice(hist(h.ProjectionY()), (-0.07, 0.07))
     plt.subplot(221)
-    plot_histogram(preproc(rs.sc_first_hits_in_B1_dz),
-                   include_errors=errors,
-                   norm=norm,
-                   log=True,
-                   xlabel="$\\Delta z_1$(cm)",
-                   title="BPIX - Layer 1")
+    hist_plot(hist(rs.sc_first_hits_in_B1_dz.ProjectionY()),
+              include_errors=errors,
+              norm=norm,
+              log=False,
+              xlabel=r"$\Delta z_1$(cm)",
+              title="BPIX - Layer 1")
 
     plt.subplot(222)
-    plot_histogram(preproc(rs.sc_first_hits_in_B2_dz),
-                   include_errors=errors,
-                   norm=norm,
-                   log=True,
-                   xlabel="$\\Delta z_1$(cm)",
-                   title="BPIX - Layer 2")
+    hist_plot(hist(rs.sc_first_hits_in_B2_dz.ProjectionY()),
+              include_errors=errors,
+              norm=norm,
+              log=False,
+              xlabel=r"$\Delta z_1$(cm)",
+              title="BPIX - Layer 2")
 
     plt.subplot(223)
-    plot_histogram(preproc(rs.sc_first_hits_in_B1_dphi),
-                   include_errors=errors,
-                   norm=norm,
-                   log=True,
-                   # ylim=(0.5E-3, 5),
-                   xlabel="$\\Delta \\phi_1$(rad)",
-                   # xlim=(-0.06, 0.06)
-                   )
+    hist_plot(preproc(rs.sc_first_hits_in_B1_dphi),
+              include_errors=errors,
+              norm=norm,
+              log=False,
+              # ylim=(0.5E-3, 5),
+              xlabel=r"$\Delta \phi_1$(rad)",
+              # xlim=(-0.06, 0.06)
+              )
 
     plt.subplot(224)
-    plot_histogram(preproc(rs.sc_first_hits_in_B2_dphi),
-                   include_errors=errors,
-                   norm=norm,
-                   log=True,
-                   # ylim=(0.5E-3, 5),
-                   xlabel="$\\Delta \\phi_1$(rad)",
-                   # xlim=(-0.06, 0.06)
-                   )
-
-    plt.tight_layout()
-    plt.savefig("figures/sc_extrapolation_first.png")
+    hist_plot(preproc(rs.sc_first_hits_in_B2_dphi),
+              include_errors=errors,
+              norm=norm,
+              log=False,
+              # ylim=(0.5E-3, 5),
+              xlabel=r"$\Delta \phi_1$(rad)",
+              # xlim=(-0.06, 0.06)
+              )
 
 
+@make_plot(scale=0.85)
 def sc_extrapolation_second(rs):
-    norm = 1
+    from scipy.stats import norm
+
+    def gauss(x, A, mu, sigma):
+        return A*norm.pdf(x, mu, sigma)
+
+    def gauss2(x, A1, mu1, sigma1, A2, mu2, sigma2):
+        return (A1*norm.pdf(x, mu1, sigma1) +
+                A2*norm.pdf(x, mu2, sigma2))
+
+    integral = 1
     errors = True
-    def preproc(h):
-        return histogram(h.ProjectionY()), (-0.07, 0.07)
+
+    def preproc(h, slice_=None):
+        h = hist(h.ProjectionY())
+        if slice_:
+            h = hist_slice(h, slice_)
+        return h
 
     plt.subplot(221)
-    plot_histogram(preproc(rs.sc_second_hits_in_B2_dz),
-                   include_errors=errors,
-                   norm=norm,
-                   log=True,
-                   xlabel="$\\Delta z_2$(cm)",
-                   title="BPIX - Layer 2")
+    hist_plot(preproc(rs.sc_second_hits_in_B2_dz),
+              include_errors=errors,
+              norm=integral,
+              log=False,
+              fit=(gauss2, (1, 0, 0.025, 1, 0, 0.005)),
+              xlabel=r"$\Delta z_2$(cm)",
+              title="BPIX - Layer 2")
 
     plt.subplot(222)
-    plot_histogram(preproc(rs.sc_second_hits_in_B3_dz),
-                   include_errors=errors,
-                   norm=norm,
-                   log=True,
-                   xlabel="$\\Delta z_2$(cm)",
-                   title="BPIX - Layer 3")
+    hist_plot(preproc(rs.sc_second_hits_in_B3_dz),
+              include_errors=errors,
+              norm=integral,
+              log=False,
+              fit=(gauss2, (1, 0, 0.025, 1, 0, 0.005)),
+              xlabel=r"$\Delta z_2$(cm)",
+              title="BPIX - Layer 3")
 
     plt.subplot(223)
-    plot_histogram(preproc(rs.sc_second_hits_in_B2_dphi),
-                   include_errors=errors,
-                   norm=norm,
-                   log=True,
-                   # ylim=(0.5E-3, 5),
-                   xlabel="$\\Delta \\phi_2$(rad)")
+    hist_plot(preproc(rs.sc_second_hits_in_B2_dphi, (-0.09, 0.09)),
+              include_errors=errors,
+              norm=integral,
+              log=False,
+              fit=(gauss2, (1, 0, 0.015, 1, 0, 0.005)),
+              xlabel=r"$\Delta \phi_2$(rad)")
 
     plt.subplot(224)
-    plot_histogram(preproc(rs.sc_second_hits_in_B3_dphi),
-                   include_errors=errors,
-                   norm=norm,
-                   log=True,
-                   # ylim=(0.5E-3, 5),
-                   xlabel="$\\Delta \\phi_2$(rad)")
-
-    plt.tight_layout()
-    plt.savefig("figures/sc_extrapolation_second.png")
-
-
-# First create a ResultSet object which loads all of the objects from output.root
-# into memory and makes them available as attributes
-if len(sys.argv) != 2:
-    raise ValueError("please supply root file")
-rs = ResultSet("DY2LL", sys.argv[1])
-
-try:
-    os.mkdir('figures')
-except FileExistsError:
-    pass
-
-plots = [
-    first_hits_v_eta,
-    second_hits_v_eta,
-    first_hits,
-    second_hits,
-    delta_phi_z_v_ladder,
-    sc_extrapolation_first,
-    sc_extrapolation_second,
-]
-
-for plot in plots:
-    plt.clf()
-    plot(rs)
-
-import glob
-from jinja2 import Environment, PackageLoader, select_autoescape
-from os.path import join
-from shutil import copytree, rmtree
-
-env = Environment(
-    loader=PackageLoader('plots', 'templates'),
-    autoescape=select_autoescape(['htm', 'html', 'xml'])
-)
-
-def render_to_file(template_name, **kwargs):
-    try:
-        os.mkdir('output')
-    except FileExistsError:
-        pass
-    with open(join('output', template_name), 'w') as tempout:
-        template = env.get_template(template_name)
-        tempout.write(template.render(**kwargs))
-
-
-try:
-    os.mkdir('output')
-except FileExistsError:
-    pass
-try:
-    rmtree('output/figures')
-except FileNotFoundError:
-    pass
-copytree('figures', 'output/figures')
-
-imgs = glob.glob("figures/*.png")
-
-
-def get_by_n(l, n=2):
-    while l:
-        yield l[:n]
-        l = l[n:]
-
-render_to_file('dashboard.htm', imgs=get_by_n(imgs, 6))
+    hist_plot(preproc(rs.sc_second_hits_in_B3_dphi, (-0.09, 0.09)),
+              include_errors=errors,
+              norm=integral,
+              log=False,
+              fit=(gauss2, (1, 0, 0.025, 1, 0, 0.005)),
+              xlabel=r"$\Delta \phi_2$(rad)")
+
+
+
+def generate_dashboard():
+    from jinja2 import Environment, PackageLoader, select_autoescape
+    from os.path import join
+    from urllib.parse import quote
+
+    env = Environment(
+        loader=PackageLoader('plots', 'templates'),
+        autoescape=select_autoescape(['htm', 'html', 'xml'])
+    )
+
+    def render_to_file(template_name, **kwargs):
+        with open(join('output', template_name), 'w') as tempout:
+            template = env.get_template(template_name)
+            tempout.write(template.render(**kwargs))
+
+
+    def get_by_n(l, n=2):
+        l = list(l)
+        while l:
+            yield l[:n]
+            l = l[n:]
+
+    render_to_file('dashboard.htm', plots=get_by_n(plot_registry.values(), 3),
+                   quote=quote)
+
+
+if __name__ == '__main__':
+    # First create a ResultSet object which loads all of the objects from output.root
+    # into memory and makes them available as attributes
+    if len(sys.argv) != 2:
+        raise ValueError("please supply root file")
+    rs = ResultSet("DY2LL", sys.argv[1])
+
+    first_hits_v_eta(rs)
+    second_hits_v_eta(rs)
+    first_hits(rs)
+    second_hits(rs)
+    delta_phi_z_v_ladder(rs)
+    sc_extrapolation_first(rs)
+    sc_extrapolation_second(rs)
+
+    generate_dashboard()

analysis/publish

@@ -0,0 +1,6 @@
+#!/usr/bin/env bash
+set -e
+echo "Will copy content of the 'output' directory to the 'public' directory on t3.unl.edu"
+scp -r output cfangmeier@t3.unl.edu:public
+echo "done"
+echo "files available at t3.unl.edu:8000/output/ (provided that the local server is running)"

analysis/templates/dashboard.htm

@@ -3,22 +3,50 @@
 <head>
   <title>EGamma Dashboard</title>
   <meta charset="utf-8">
-  <meta http-equiv="refresh" content="10">
   <meta name="viewport" content="width=device-width, initial-scale=2">
   <link rel="stylesheet" href="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.7/css/bootstrap.min.css">
   <script src="https://ajax.googleapis.com/ajax/libs/jquery/3.2.1/jquery.min.js"></script>
   <script src="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.7/js/bootstrap.min.js"></script>
+
+<script type="text/x-mathjax-config">
+MathJax.Hub.Config({
+  config: ["MMLorHTML.js"],
+  jax: ["input/TeX", "output/HTML-CSS", "output/NativeMML"],
+  extensions: ["MathMenu.js", "MathZoom.js"]
+});
+</script>
+<script type="text/javascript" src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS_HTML-full"> </script>
 </head>
 <body>
 <div class="container-fluid">
-  {% for img_row in imgs %}
+  {% for plot_row in plots %}
   <div class="row">
-    {% for img in img_row %}
-    <div class="col-md-2">
-      <div class="thumbnail">
-        <a href="{{img}}"><img src="{{img}}" style="width:100%"></a>
+    {% for plot in plot_row %}
+    <div class="col-md-4">
+      <div class="well">
+        {% if plot.title %}
+        <h3 class="text-center">{{ plot.title }}</h3>
+        {% endif %}
+        <a href="{{ plot.filename }}"><img src="{{ plot.filename }}" style="width:100%"></a>
         <div class="caption">
-          <p class="text-center"> {{img}} </p>
+          <p class="text-center"> {{ plot.name }} </p>
+          {% if plot.desc %}
+          <hr>
+          <p class="text-left">{{ plot.desc|safe }}</p>
+          <hr>
+          {% endif %}
+          {% if plot.args %}
+          <p class="text-left"><strong>Plot Arguments</strong></p>
+          <table class="table table-hover">
+            <tbody>
+              {% for key, val in plot.args.items() %}
+              <tr>
+                <td>{{ key }}</td> <td>{{ val }}</td>
+              </tr>
+              {% endfor %}
+            </tbody>
+          </table>
+          {% endif %}
         </div>
       </div>
     </div>

analysis/tracking_validation.cpp

@@ -132,7 +132,7 @@ void setup_skipped_layer_hit_pairs(TrackingDataSet& tds){
                          vector<SimHit>,
                          vector<Track>,
                          vector<SimTrack>,
-                         int, int, int)>("find_matched_nth_hit_in_layer_with_skip",
+                         int, int, int, bool)>("find_matched_nth_hit_in_layer_with_skip",
         FUNC(([](const vector<Seed>& seeds,
                  const vector<PixRecHit>& pixrec_hits,
                  const vector<SimHit>& sim_hits,
@@ -140,7 +140,8 @@ void setup_skipped_layer_hit_pairs(TrackingDataSet& tds){
                  const vector<SimTrack>& sim_tracks,
                  const unsigned int&& det,
                  const unsigned int&& pix_layer,
-                 const unsigned int&& skip){
+                 const unsigned int&& skip,
+                 const bool&& first){
             vector<HitPair> matched_hits;
             for(const Track &trk : tracks){ // loop over all tracks
                 const Seed &seed = seeds[trk.seedIdx];
@@ -154,14 +155,27 @@ void setup_skipped_layer_hit_pairs(TrackingDataSet& tds){
 
                 const PixRecHit &rec_hit1 = pixrec_hits[seed.hitIdx[0]];
                 const PixRecHit &rec_hit2 = pixrec_hits[seed.hitIdx[1]];
-                if(in_det(rec_hit1, det, pix_layer) && in_det(rec_hit2, det, pix_layer+1+skip)){
-                    // We have the RecHit we want to work with, now find a properly* matched SimHit
-                    if(rec_hit1.simHitIdx.size() == 0) continue; // if rechit is matched to no simhits, give up.
-                    if(trk.simTrkIdx.size() == 0) continue; // if track is matched to no simtracks, give up.
-                    const int &simTrkIdx = trk.simTrkIdx[0];  // get first matched track
-                    for(const int& simHitIdx : sim_tracks[simTrkIdx].simHitIdx){ // loop over SimHits from SimTrack
-                        if(simHitIdx == rec_hit1.simHitIdx[0]) // take first matched simhit (should be the closest one)
-                            matched_hits.push_back({rec_hit1, sim_hits[rec_hit1.simHitIdx[0]]});
+                if(first){ // Looking at first hit in the pair
+                    if(in_det(rec_hit1, det, pix_layer) && in_det(rec_hit2, det, pix_layer+1+skip)){
+                        // We have the RecHit we want to work with, now find a properly* matched SimHit
+                        if(rec_hit1.simHitIdx.size() == 0) continue; // if rechit is matched to no simhits, give up.
+                        if(trk.simTrkIdx.size() == 0) continue; // if track is matched to no simtracks, give up.
+                        const int &simTrkIdx = trk.simTrkIdx[0];  // get first matched track
+                        for(const int& simHitIdx : sim_tracks[simTrkIdx].simHitIdx){ // loop over SimHits from SimTrack
+                            if(simHitIdx == rec_hit1.simHitIdx[0]) // take first matched simhit (should be the closest one)
+                                matched_hits.push_back({rec_hit1, sim_hits[rec_hit1.simHitIdx[0]]});
+                        }
+                    }
+                }else{  // Looking at second hit in the pair
+                    if(in_det(rec_hit2, det, pix_layer) && in_det(rec_hit1, det, pix_layer-1-skip)){
+                        // We have the RecHit we want to work with, now find a properly* matched SimHit
+                        if(rec_hit2.simHitIdx.size() == 0) continue; // if rechit is matched to no simhits, give up.
+                        if(trk.simTrkIdx.size() == 0) continue; // if track is matched to no simtracks, give up.
+                        const int &simTrkIdx = trk.simTrkIdx[0];  // get first matched track
+                        for(const int& simHitIdx : sim_tracks[simTrkIdx].simHitIdx){ // loop over SimHits from SimTrack
+                            if(simHitIdx == rec_hit1.simHitIdx[0]) // take first matched simhit (should be the closest one)
+                                matched_hits.push_back({rec_hit2, sim_hits[rec_hit2.simHitIdx[0]]});
+                        }
                     }
                 }
             }
@@ -176,6 +190,9 @@ void setup_skipped_layer_hit_pairs(TrackingDataSet& tds){
     auto skip_zero = constant("skip_zero", 0);
     auto skip_one  = constant("skip_one", 1);
 
+    auto first = constant("first", true);
+    auto second = constant("second", false);
+
 
     // First hit on 1st/2nd bpix layers and second hit in 2nd/3rd
     auto first_hits_in_B1_skip_zero = fv::tup_apply(find_matched_nth_hit_in_layer_with_skip,

filval

@@ -1 +1 @@
-Subproject commit 1a85d61848cbde369c479b434883640cfd08e3f0
+Subproject commit 8fe4cbd0fe9636e6a421d18df688bc942455551b