adds studies using super-cluster information

.gitignore

@@ -3,6 +3,10 @@ legacy/
 data/
 build/
 env/
+analysis/output/
+analysis/figures/
 
 tags
 .ipynb_checkpoints/
+
+*/__pycache__/

analysis/TrackingNtuple.h

@@ -5,8 +5,8 @@
 // found on file: trackingNtuple.root
 //////////////////////////////////////////////////////////
 
-#ifndef tree_h
-#define tree_h
+#ifndef TrackingNtuple_h
+#define TrackingNtuple_h
 
 #include <TROOT.h>
 #include <TChain.h>
@@ -1181,4 +1181,4 @@ Int_t TrackingNtuple::Cut(Long64_t entry)
 // returns -1 otherwise.
    return 1;
 }
-#endif // #ifdef tree_cxx
+#endif // TrackingNtuple_h

analysis/common.hpp

@@ -0,0 +1,43 @@
+#ifndef COMMON_HPP
+#define COMMON_HPP
+
+#include <cmath>
+
+template<typename A, typename B>
+float displacement(const A& a, const B& b){
+    return std::sqrt(pow(a.x-b.x, 2) +
+                     pow(a.y-b.y, 2) +
+                     pow(a.z-b.z, 2));
+}
+
+template<typename T>
+float rho(const T& t){
+    return std::sqrt(pow(t.x,2)+pow(t.y,2));
+}
+
+float pseudorapidity(float z, float r){
+    float theta = atan2(r, z);
+    return -log(tan(theta/2.0));
+}
+
+float pseudorapidity(float x, float y, float z){
+    float r = sqrt(x*x + y*y);
+    return pseudorapidity(z, r);
+}
+
+template<typename T>
+float pseudorapidity(const T &t){
+    return pseudorapidity(t.x, t.y, t.z);
+}
+
+float pseudorapidityP(float px, float py, float pz){
+    float p_mag = sqrt(px*px + py*py + pz*pz);
+    return 0.5*log((p_mag+pz)/(p_mag-pz));
+}
+
+template<typename T>
+float pseudorapidityP(const T &t){
+    return pseudorapidityP(t.px, t.py, t.pz);
+}
+
+#endif // COMMON_HPP

analysis/hit_resolution_plots.py

@@ -1,184 +0,0 @@
-
-# coding: utf-8
-
-import os
-import sys
-import matplotlib as mpl
-import matplotlib.pyplot as plt
-
-sys.path.append("../filval/python/")
-from utils import ResultSet
-from plotter import plot_histogram, plot_histogram2d
-mpl.rc('text', usetex=True)
-mpl.rc('savefig', dpi=180)
-mpl.rc('savefig', transparent=True)
-
-# First create a ResultSet object which loads all of the objects from output.root
-# into memory and makes them available as attributes
-rs = ResultSet("DY2LL",  "/home/caleb/Sources/EGamma/build/output.root")
-
-try:
-    os.mkdir('figures')
-except FileExistsError:
-    pass
-
-scale = 0.85
-plt.figure(figsize=(scale*10, scale*10))
-
-plt.subplot(321)
-plot_histogram2d(rs.dphi_v_eta_first_hits_in_B1,
-                 ylabel="$\\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")
-
-plt.subplot(323)
-plot_histogram2d(rs.dz_v_eta_first_hits_in_B1,
-                 ylabel="$\\Delta \\textrm{z}_1$(cm)")
-
-plt.subplot(324)
-plot_histogram2d(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$")
-
-plt.subplot(326)
-plot_histogram2d(rs.dr_v_eta_first_hits_in_B2,
-                 xlabel="$\\eta$")
-
-plt.tight_layout()
-plt.savefig("figures/first_hits_v_eta.png")
-
-
-plt.clf()
-plt.subplot(321)
-plot_histogram2d(rs.dphi_v_eta_second_hits_in_B2,
-                 ylabel="$\\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")
-
-plt.subplot(323)
-plot_histogram2d(rs.dz_v_eta_second_hits_in_B2,
-                 ylabel="$\\Delta \\textrm{z}_2$(cm)")
-
-plt.subplot(324)
-plot_histogram2d(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$")
-
-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")
-
-
-plt.clf()
-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")
-
-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")
-
-plt.subplot(323)
-plot_histogram(rs.dz_v_eta_first_hits_in_B1.ProjectionY(), include_errors=True,
-               xlabel="$\\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)")
-
-plt.subplot(325)
-plot_histogram(rs.dr_v_eta_first_hits_in_B1.ProjectionY(), include_errors=True,
-               xlabel="$\\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")
-
-
-plt.clf()
-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")
-
-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")
-
-plt.subplot(323)
-plot_histogram(rs.dz_v_eta_second_hits_in_B2.ProjectionY(), include_errors=True,
-               xlabel="$\\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)")
-
-plt.subplot(325)
-plot_histogram(rs.dr_v_eta_second_hits_in_B2.ProjectionY(), include_errors=True,
-               xlabel="$\\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)")
-
-plt.tight_layout()
-plt.savefig("figures/second_hits.png")
-
-
-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)
-
-    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, .6, txt.format(var, "even", even_mean, unit), transform=axes.transAxes)
-
-
-plt.clf()
-
-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,
-              r"\phi", 10**6, "urad", xlabel=r"$\Delta \phi_1$(rad)", title="BPIX - Layer 1")
-
-plt.subplot(222)
-even_odd_plot(rs.dphi_v_eta_first_hits_in_B2_even_ladder, rs.dphi_v_eta_first_hits_in_B2_odd_ladder,
-              r"\phi", 10**6, "urad", xlabel=r"$\Delta \phi_1$(rad)", title="BPIX - Layer 2")
-plt.legend()
-
-plt.subplot(223)
-even_odd_plot(rs.dz_v_eta_first_hits_in_B1_even_ladder, rs.dz_v_eta_first_hits_in_B1_odd_ladder,
-              "z", 10**4, "um", xlabel=r"$\Delta z_1$(cm)")
-
-plt.subplot(224)
-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")

analysis/obj_types.cpp

@@ -9,8 +9,6 @@ using namespace std;
 using namespace fv;
 using namespace fv::root;
 
-using namespace std;
-
 typedef TreeDataSet<TrackingNtuple> TrackingDataSet;
 
 //#^#=========================================================================
@@ -46,26 +44,26 @@ struct PixRecHit {
 
 Value<vector<PixRecHit>>*
 register_pixrec_hits(TrackingDataSet &tds){
-    auto& builder = func<std::vector<PixRecHit>(vector<unsigned short>,
-                                                vector<unsigned short>,
-                                                vector<unsigned short>,
-                                                vector<unsigned int>,
-                                                vector<unsigned short>,
-                                                vector<float>,
-                                                vector<float>,
-                                                vector<float>,
-                                                vector<float>,
-                                                vector<float>,
-                                                vector<float>,
-                                                vector<float>,
-                                                vector<float>,
-                                                vector<float>,
-                                                vector<float>,
-                                                vector<float>,
-                                                vector<vector<int>>,
-                                                vector<vector<int>>,
-                                                vector<vector<int>>,
-                                                vector<vector<float>>)>("build_pixrec_hits",
+    auto builder = func<std::vector<PixRecHit>(vector<unsigned short>,
+                                               vector<unsigned short>,
+                                               vector<unsigned short>,
+                                               vector<unsigned int>,
+                                               vector<unsigned short>,
+                                               vector<float>,
+                                               vector<float>,
+                                               vector<float>,
+                                               vector<float>,
+                                               vector<float>,
+                                               vector<float>,
+                                               vector<float>,
+                                               vector<float>,
+                                               vector<float>,
+                                               vector<float>,
+                                               vector<float>,
+                                               vector<vector<int>>,
+                                               vector<vector<int>>,
+                                               vector<vector<int>>,
+                                               vector<vector<float>>)>("build_pixrec_hits",
         FUNC(([](const vector<unsigned short>& det,
                  const vector<unsigned short>& lay,
                  const vector<unsigned short>& ladder_blade,
@@ -87,7 +85,7 @@ register_pixrec_hits(TrackingDataSet &tds){
                  const vector<vector<int>>& simHitIdx,
                  const vector<vector<float>>& chargeFraction){
             std::vector<PixRecHit> pixrec_hits;
-            for(int i=0; i<det.size(); i++)
+            for(unsigned int i=0; i<det.size(); i++)
                 pixrec_hits.push_back({i,
                                        det[i],
                                        lay[i],
@@ -114,7 +112,7 @@ register_pixrec_hits(TrackingDataSet &tds){
         })));
 
 
-    return apply(builder,
+    return tup_apply(builder,
         fv::tuple(tds.track_branch_obj<vector<unsigned short>>("pix_det"),
                   tds.track_branch_obj<vector<unsigned short>>("pix_lay"),
                   tds.track_branch_obj<vector<unsigned short>>("pix_ladder_blade"),
@@ -160,7 +158,7 @@ struct SimHit {
 
 Value<vector<SimHit>>*
 register_sim_hits(TrackingDataSet &tds){
-    auto& builder =
+    auto builder =
         func<std::vector<SimHit>(vector<unsigned short>, //  det;
                                  vector<unsigned short>, //  lay;
                                  vector<unsigned int>,   //  detId;
@@ -190,7 +188,7 @@ register_sim_hits(TrackingDataSet &tds){
                  const vector<vector<int>>&    hitType
                  ){
             std::vector<SimHit> sim_hits;
-            for(int i=0; i<det.size(); i++)
+            for(unsigned int i=0; i<det.size(); i++)
                 sim_hits.push_back({i,
                                     det[i],
                                     lay[i],
@@ -209,7 +207,7 @@ register_sim_hits(TrackingDataSet &tds){
         })));
 
 
-    return apply(builder,
+    return tup_apply(builder,
         fv::tuple(tds.track_branch_obj<vector<unsigned short>>("simhit_det"),
                   tds.track_branch_obj<vector<unsigned short>>("simhit_lay"),
                   tds.track_branch_obj<vector<unsigned int>>("simhit_detId"),
@@ -284,7 +282,7 @@ struct Seed {
 
 Value<vector<Seed>>*
 register_seeds(TrackingDataSet &tds){
-    auto& builder = func<std::vector<Seed>(vector<short>,          // fitok;
+    auto builder = func<std::vector<Seed>(vector<short>,          // fitok;
                                                vector<float>,          // px;
                                                vector<float>,          // py;
                                                vector<float>,          // pz;
@@ -382,7 +380,7 @@ register_seeds(TrackingDataSet &tds){
                const vector<unsigned int>& ecalDriven,
                const vector<unsigned int>& trkDriven){
             std::vector<Seed> seeds;
-            for(int i=0; i<fitok.size(); i++)
+            for(unsigned int i=0; i<fitok.size(); i++)
                 seeds.push_back({i,
                                  fitok[i],
                                  px[i],
@@ -436,7 +434,7 @@ register_seeds(TrackingDataSet &tds){
         })));
 
 
-    return apply(builder,
+    return tup_apply(builder,
         fv::tuple(tds.track_branch_obj<vector<short>>("see_fitok"),
                  tds.track_branch_obj<vector<float>>("see_px"),
                  tds.track_branch_obj<vector<float>>("see_py"),
@@ -551,59 +549,60 @@ struct Track{
 
 Value<vector<Track>>*
 register_tracks(TrackingDataSet &tds){
-    auto& builder= func<std::vector<Track>(vector<float>,         //  px;
-                                                vector<float>,         //  py;
-                                                vector<float>,         //  pz;
-                                                vector<float>,         //  pt;
-                                                vector<float>,         //  inner_px;
-                                                vector<float>,         //  inner_py;
-                                                vector<float>,         //  inner_pz;
-                                                vector<float>,         //  inner_pt;
-                                                vector<float>,         //  outer_px;
-                                                vector<float>,         //  outer_py;
-                                                vector<float>,         //  outer_pz;
-                                                vector<float>,         //  outer_pt;
-                                                vector<float>,         //  eta;
-                                                vector<float>,         //  lambda;
-                                                vector<float>,         //  cotTheta;
-                                                vector<float>,         //  phi;
-                                                vector<float>,         //  dxy;
-                                                vector<float>,         //  dz;
-                                                vector<float>,         //  ptErr;
-                                                vector<float>,         //  etaErr;
-                                                vector<float>,         //  lambdaErr;
-                                                vector<float>,         //  phiErr;
-                                                vector<float>,         //  dxyErr;
-                                                vector<float>,         //  dzErr;
-                                                vector<float>,         //  refpoint_x;
-                                                vector<float>,         //  refpoint_y;
-                                                vector<float>,         //  refpoint_z;
-                                                vector<float>,         //  nChi2;
-                                                vector<int>,           //  q;
-                                                vector<unsigned int>,  //  nValid;
-                                                vector<unsigned int>,  //  nInvalid;
-                                                vector<unsigned int>,  //  nPixel;
-                                                vector<unsigned int>,  //  nStrip;
-                                                vector<unsigned int>,  //  nPixelLay;
-                                                vector<unsigned int>,  //  nStripLay;
-                                                vector<unsigned int>,  //  n3DLay;
-                                                vector<unsigned int>,  //  nOuterLost;
-                                                vector<unsigned int>,  //  nInnerLost;
-                                                vector<unsigned int>,  //  algo;
-                                                vector<unsigned int>,  //  originalAlgo;
-                                                vector<ULong64_t>,     //  algoMask;
-                                                vector<unsigned int>,  //  stopReason;
-                                                vector<short>,         //  isHP;
-                                                vector<int>,           //  seedIdx;
-                                                vector<float>,         //  vtxx;
-                                                vector<float>,         //  vtxy;
-                                                vector<float>,         //  vtxz;
-                                                vector<vector<float>>, //  shareFrac;
-                                                vector<vector<int>>,   //  simTrkIdx;
-                                                vector<vector<int>>,   //  hitIdx;
-                                                vector<vector<int>>,   //  hitType;
-                                                vector<int>,           //  genIdx;
-                                                vector<float>          //  genDR;
+    auto builder =
+        func<std::vector<Track>(vector<float>,         //  px;
+                                vector<float>,         //  py;
+                                vector<float>,         //  pz;
+                                vector<float>,         //  pt;
+                                vector<float>,         //  inner_px;
+                                vector<float>,         //  inner_py;
+                                vector<float>,         //  inner_pz;
+                                vector<float>,         //  inner_pt;
+                                vector<float>,         //  outer_px;
+                                vector<float>,         //  outer_py;
+                                vector<float>,         //  outer_pz;
+                                vector<float>,         //  outer_pt;
+                                vector<float>,         //  eta;
+                                vector<float>,         //  lambda;
+                                vector<float>,         //  cotTheta;
+                                vector<float>,         //  phi;
+                                vector<float>,         //  dxy;
+                                vector<float>,         //  dz;
+                                vector<float>,         //  ptErr;
+                                vector<float>,         //  etaErr;
+                                vector<float>,         //  lambdaErr;
+                                vector<float>,         //  phiErr;
+                                vector<float>,         //  dxyErr;
+                                vector<float>,         //  dzErr;
+                                vector<float>,         //  refpoint_x;
+                                vector<float>,         //  refpoint_y;
+                                vector<float>,         //  refpoint_z;
+                                vector<float>,         //  nChi2;
+                                vector<int>,           //  q;
+                                vector<unsigned int>,  //  nValid;
+                                vector<unsigned int>,  //  nInvalid;
+                                vector<unsigned int>,  //  nPixel;
+                                vector<unsigned int>,  //  nStrip;
+                                vector<unsigned int>,  //  nPixelLay;
+                                vector<unsigned int>,  //  nStripLay;
+                                vector<unsigned int>,  //  n3DLay;
+                                vector<unsigned int>,  //  nOuterLost;
+                                vector<unsigned int>,  //  nInnerLost;
+                                vector<unsigned int>,  //  algo;
+                                vector<unsigned int>,  //  originalAlgo;
+                                vector<ULong64_t>,     //  algoMask;
+                                vector<unsigned int>,  //  stopReason;
+                                vector<short>,         //  isHP;
+                                vector<int>,           //  seedIdx;
+                                vector<float>,         //  vtxx;
+                                vector<float>,         //  vtxy;
+                                vector<float>,         //  vtxz;
+                                vector<vector<float>>, //  shareFrac;
+                                vector<vector<int>>,   //  simTrkIdx;
+                                vector<vector<int>>,   //  hitIdx;
+                                vector<vector<int>>,   //  hitType;
+                                vector<int>,           //  genIdx;
+                                vector<float>          //  genDR;
         )>("build_tracks", FUNC(([](const vector<float>&         px,
                                     const vector<float>&         py,
                                     const vector<float>&         pz,
@@ -658,7 +657,7 @@ register_tracks(TrackingDataSet &tds){
                                     const vector<int>&           genIdx,
                                     const vector<float>&         genDR){
             std::vector<Track> tracks;
-            for(int i=0; i<px.size(); i++)
+            for(unsigned int i=0; i<px.size(); i++)
                 tracks.push_back({i,
                                   px[i],
                                   py[i],
@@ -715,7 +714,7 @@ register_tracks(TrackingDataSet &tds){
                                   genDR[i]});
             return tracks;
         })));
-    return apply(builder,
+    return tup_apply(builder,
         fv::tuple(tds.track_branch_obj<vector<float>>("trk_px"),
                   tds.track_branch_obj<vector<float>>("trk_py"),
                   tds.track_branch_obj<vector<float>>("trk_pz"),
@@ -811,34 +810,35 @@ struct SimTrack{
 
 Value<vector<SimTrack>>*
 register_sim_tracks(TrackingDataSet &tds){
-    auto& builder = func<std::vector<SimTrack>(vector<int>,           //  event;
-                                               vector<int>,           //  bunchCrossing;
-                                               vector<int>,           //  pdgId;
-                                               vector<float>,         //  px;
-                                               vector<float>,         //  py;
-                                               vector<float>,         //  pz;
-                                               vector<float>,         //  pt;
-                                               vector<float>,         //  eta;
-                                               vector<float>,         //  phi;
-                                               vector<float>,         //  pca_pt;
-                                               vector<float>,         //  pca_eta;
-                                               vector<float>,         //  pca_lambda;
-                                               vector<float>,         //  pca_cotTheta;
-                                               vector<float>,         //  pca_phi;
-                                               vector<float>,         //  pca_dxy;
-                                               vector<float>,         //  pca_dz;
-                                               vector<int>,           //  q;
-                                               vector<unsigned int>,  //  nValid;
-                                               vector<unsigned int>,  //  nPixel;
-                                               vector<unsigned int>,  //  nStrip;
-                                               vector<unsigned int>,  //  nLay;
-                                               vector<unsigned int>,  //  nPixelLay;
-                                               vector<unsigned int>,  //  n3DLay;
-                                               vector<vector<int>>,   //  trkIdx;
-                                               vector<vector<float>>, //  shareFrac;
-                                               vector<int>,           //  parentVtxIdx;
-                                               vector<vector<int>>,   //  decayVtxIdx;
-                                               vector<vector<int>>    //  simHitIdx;
+    auto builder =
+        func<std::vector<SimTrack>(vector<int>,           //  event;
+                                   vector<int>,           //  bunchCrossing;
+                                   vector<int>,           //  pdgId;
+                                   vector<float>,         //  px;
+                                   vector<float>,         //  py;
+                                   vector<float>,         //  pz;
+                                   vector<float>,         //  pt;
+                                   vector<float>,         //  eta;
+                                   vector<float>,         //  phi;
+                                   vector<float>,         //  pca_pt;
+                                   vector<float>,         //  pca_eta;
+                                   vector<float>,         //  pca_lambda;
+                                   vector<float>,         //  pca_cotTheta;
+                                   vector<float>,         //  pca_phi;
+                                   vector<float>,         //  pca_dxy;
+                                   vector<float>,         //  pca_dz;
+                                   vector<int>,           //  q;
+                                   vector<unsigned int>,  //  nValid;
+                                   vector<unsigned int>,  //  nPixel;
+                                   vector<unsigned int>,  //  nStrip;
+                                   vector<unsigned int>,  //  nLay;
+                                   vector<unsigned int>,  //  nPixelLay;
+                                   vector<unsigned int>,  //  n3DLay;
+                                   vector<vector<int>>,   //  trkIdx;
+                                   vector<vector<float>>, //  shareFrac;
+                                   vector<int>,           //  parentVtxIdx;
+                                   vector<vector<int>>,   //  decayVtxIdx;
+                                   vector<vector<int>>    //  simHitIdx;
             )>("build_sim_tracks", FUNC(([](
                const vector<int>&           event,
                const vector<int>&           bunchCrossing,
@@ -869,7 +869,7 @@ register_sim_tracks(TrackingDataSet &tds){
                const vector<vector<int>>&   decayVtxIdx,
                const vector<vector<int>>&   simHitIdx){
             std::vector<SimTrack> sim_tracks;
-            for(int i=0; i<event.size(); i++)
+            for(unsigned int i=0; i<event.size(); i++)
                 sim_tracks.push_back({i,
                                       event[i],
                                       bunchCrossing[i],
@@ -901,7 +901,7 @@ register_sim_tracks(TrackingDataSet &tds){
                                       simHitIdx[i]});
             return sim_tracks;
         })));
-    return apply(builder,
+    return tup_apply(builder,
         fv::tuple(tds.track_branch_obj<vector<int>>("sim_event"),
                   tds.track_branch_obj<vector<int>>("sim_bunchCrossing"),
                   tds.track_branch_obj<vector<int>>("sim_pdgId"),
@@ -932,3 +932,117 @@ register_sim_tracks(TrackingDataSet &tds){
                   tds.track_branch_obj<vector<vector<int>>>("sim_simHitIdx")),
         "sim_tracks");
 }
+
+//#^#=========================================================================
+//==============SUPER CLUSTERS================================================
+//============================================================================
+struct SuperCluster {
+    unsigned int          idx;
+    float                 e;
+    float                 px;
+    float                 py;
+    float                 pz;
+    float                 x;
+    float                 y;
+    float                 z;
+    vector<int>           charge;
+    vector<int>           lay1;
+    vector<int>           lay2;
+    vector<int>           subDet1;
+    vector<int>           subDet2;
+    vector<float>         dRz1;
+    vector<float>         dPhi1;
+    vector<float>         dRz2;
+    vector<float>         dPhi2;
+    vector<int>           seedType;
+    vector<unsigned char> hitsMask;
+};
+
+Value<vector<SuperCluster>>*
+register_super_clusters(TrackingDataSet &tds){
+    auto builder =
+        func<std::vector<SuperCluster>(vector<float>,                 // e,
+                                       vector<float>,                 // px,
+                                       vector<float>,                 // py,
+                                       vector<float>,                 // pz,
+                                       vector<float>,                 // x,
+                                       vector<float>,                 // y,
+                                       vector<float>,                 // z,
+                                       vector<vector<int>>,           // charge;
+                                       vector<vector<int>>,           // lay1;
+                                       vector<vector<int>>,           // lay2;
+                                       vector<vector<int>>,           // subDet1;
+                                       vector<vector<int>>,           // subDet2;
+                                       vector<vector<float>>,         // dRz1;
+                                       vector<vector<float>>,         // dPhi1;
+                                       vector<vector<float>>,         // dRz2;
+                                       vector<vector<float>>,         // dPhi2;
+                                       vector<vector<int>>,           // seedType;
+                                       vector<vector<unsigned char>>  // hitsMask;
+                                        )>("build_super_clusters",
+        FUNC(([](const vector<float>&                  e,
+                 const vector<float>&                  px,
+                 const vector<float>&                  py,
+                 const vector<float>&                  pz,
+                 const vector<float>&                  x,
+                 const vector<float>&                  y,
+                 const vector<float>&                  z,
+                 const vector<vector<int>>&            charge,
+                 const vector<vector<int>>&            lay1,
+                 const vector<vector<int>>&            lay2,
+                 const vector<vector<int>>&            subDet1,
+                 const vector<vector<int>>&            subDet2,
+                 const vector<vector<float>>&          dRz1,
+                 const vector<vector<float>>&          dPhi1,
+                 const vector<vector<float>>&          dRz2,
+                 const vector<vector<float>>&          dPhi2,
+                 const vector<vector<int>>&            seedType,
+                 const vector<vector<unsigned char>>& hitsMask){
+            std::vector<SuperCluster> super_clusters;
+            for(unsigned int i=0; i<e.size(); i++)
+                super_clusters.push_back({i,
+                                          e[i],
+                                          px[i],
+                                          py[i],
+                                          pz[i],
+                                          x[i],
+                                          y[i],
+                                          z[i],
+                                          charge[i],
+                                          lay1[i],
+                                          lay2[i],
+                                          subDet1[i],
+                                          subDet2[i],
+                                          dRz1[i],
+                                          dPhi1[i],
+                                          dRz2[i],
+                                          dPhi2[i],
+                                          seedType[i],
+                                          hitsMask[i]});
+            return super_clusters;
+        })));
+
+    return tup_apply(builder,
+        fv::tuple(tds.track_branch_obj<vector<float>>("scl_e"),
+                  tds.track_branch_obj<vector<float>>("scl_px"),
+                  tds.track_branch_obj<vector<float>>("scl_py"),
+                  tds.track_branch_obj<vector<float>>("scl_pz"),
+                  tds.track_branch_obj<vector<float>>("scl_x"),
+                  tds.track_branch_obj<vector<float>>("scl_y"),
+                  tds.track_branch_obj<vector<float>>("scl_z"),
+                  tds.track_branch_obj<vector<vector<int>>>("scl_charge"),
+                  tds.track_branch_obj<vector<vector<int>>>("scl_lay1"),
+                  tds.track_branch_obj<vector<vector<int>>>("scl_lay2"),
+                  tds.track_branch_obj<vector<vector<int>>>("scl_subDet1"),
+                  tds.track_branch_obj<vector<vector<int>>>("scl_subDet2"),
+                  tds.track_branch_obj<vector<vector<float>>>("scl_dRz1"),
+                  tds.track_branch_obj<vector<vector<float>>>("scl_dPhi1"),
+                  tds.track_branch_obj<vector<vector<float>>>("scl_dRz2"),
+                  tds.track_branch_obj<vector<vector<float>>>("scl_dPhi2"),
+                  tds.track_branch_obj<vector<vector<int>>>("scl_seedType"),
+                  tds.track_branch_obj<vector<vector<unsigned char>>>("scl_hitsMask")),
+        "super_clusters");
+}
+
+
+

analysis/plots.py

@@ -0,0 +1,333 @@
+#!/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)
+
+
+
+def first_hits_v_eta(rs):
+    scale = 0.85
+    plt.figure(figsize=(scale*10, scale*10))
+
+    plt.subplot(321)
+    plot_histogram2d(rs.dphi_v_eta_first_hits_in_B1,
+                     ylabel="$\\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")
+
+    plt.subplot(323)
+    plot_histogram2d(rs.dz_v_eta_first_hits_in_B1,
+                     ylabel="$\\Delta \\textrm{z}_1$(cm)")
+
+    plt.subplot(324)
+    plot_histogram2d(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$")
+
+    plt.subplot(326)
+    plot_histogram2d(rs.dr_v_eta_first_hits_in_B2,
+                     xlabel="$\\eta$")
+
+    plt.tight_layout()
+    plt.savefig("figures/first_hits_v_eta.png")
+
+
+plt.clf()
+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")
+
+    plt.subplot(322)
+    plot_histogram2d(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)")
+
+    plt.subplot(324)
+    plot_histogram2d(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$")
+
+    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")
+
+
+plt.clf()
+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")
+
+    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")
+
+    plt.subplot(323)
+    plot_histogram(rs.dz_v_eta_first_hits_in_B1.ProjectionY(), include_errors=True,
+                   xlabel="$\\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)")
+
+    plt.subplot(325)
+    plot_histogram(rs.dr_v_eta_first_hits_in_B1.ProjectionY(), include_errors=True,
+                   xlabel="$\\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")
+
+
+plt.clf()
+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")
+
+    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")
+
+    plt.subplot(323)
+    plot_histogram(rs.dz_v_eta_second_hits_in_B2.ProjectionY(), include_errors=True,
+                   xlabel="$\\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)")
+
+    plt.subplot(325)
+    plot_histogram(rs.dr_v_eta_second_hits_in_B2.ProjectionY(), include_errors=True,
+                   xlabel="$\\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)")
+
+    plt.tight_layout()
+    plt.savefig("figures/second_hits.png")
+
+
+plt.clf()
+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)
+
+        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, .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,
+                  r"\phi", 10**6, "urad", xlabel=r"$\Delta \phi_1$(rad)", title="BPIX - Layer 1")
+
+    plt.subplot(222)
+    even_odd_plot(rs.dphi_v_eta_first_hits_in_B2_even_ladder, rs.dphi_v_eta_first_hits_in_B2_odd_ladder,
+                  r"\phi", 10**6, "urad", xlabel=r"$\Delta \phi_1$(rad)", title="BPIX - Layer 2")
+    plt.legend()
+
+    plt.subplot(223)
+    even_odd_plot(rs.dz_v_eta_first_hits_in_B1_even_ladder, rs.dz_v_eta_first_hits_in_B1_odd_ladder,
+                  "z", 10**4, "um", xlabel=r"$\Delta z_1$(cm)")
+
+    plt.subplot(224)
+    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")
+
+def sc_extrapolation_first(rs):
+    # Raphael's plots
+    norm = 1
+    errors = True
+    def preproc(h):
+        return histogram_slice(histogram(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")
+
+    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")
+
+    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)
+                   )
+
+    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")
+
+
+def sc_extrapolation_second(rs):
+    norm = 1
+    errors = True
+    def preproc(h):
+        return histogram(h.ProjectionY()), (-0.07, 0.07)
+
+    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")
+
+    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")
+
+    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)")
+
+    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))

analysis/templates/dashboard.htm

@@ -0,0 +1,30 @@
+<!DOCTYPE html>
+<html lang="en">
+<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>
+</head>
+<body>
+<div class="container-fluid">
+  {% for img_row in imgs %}
+  <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>
+        <div class="caption">
+          <p class="text-center"> {{img}} </p>
+        </div>
+      </div>
+    </div>
+    {% endfor %}
+  </div>
+  {% endfor %}
+</div>
+</body>
+</html>

analysis/tracking_validation.cpp

@@ -4,7 +4,6 @@
 #include <utility>
 #include <numeric>
 #include <limits>
-#include <cmath>
 #include <TSystem.h>
 
 #include "filval/filval.hpp"
@@ -13,14 +12,14 @@
 #include <boost/range/combine.hpp>
 #include <boost/format.hpp>
 
-#include "TrackingNtuple.h"
+#include "analysis/TrackingNtuple.h"
 #include "analysis/obj_types.cpp"
+#include "analysis/common.hpp"
 
 using namespace std;
 using namespace fv;
 using namespace fv::root;
 
-typedef std::tuple<Seed, PixRecHit, Track, SimHit, SimTrack> MatchedTrack;
 typedef std::pair<std::vector<float>,std::vector<float>> pair_vec;
 
 #define HIT_TYPE_PIXEL 0
@@ -31,36 +30,64 @@ typedef std::pair<std::vector<float>,std::vector<float>> pair_vec;
 #define PIXEL_ENDCAP 2
 std::map<int,string> subdet_names = {{1, "BPIX"}, {2, "FPIX"}};
 
-template<typename A, typename B>
-float displacement(const A& a, const B& b){
-    return std::sqrt(pow(a.x-b.x, 2) +
-                     pow(a.y-b.y, 2) +
-                     pow(a.z-b.z, 2));
-}
 
-template<typename T>
-float rho(const T& t){
-    return std::sqrt(pow(t.x,2)+pow(t.y,2));
-}
+typedef std::tuple<PixRecHit, SimHit> HitPair;
 
-bool in_det(const MatchedTrack &mt, int &&det, int &&layer){
-    auto& hit = std::get<PixRecHit>(mt);
-    return hit.det == det && hit.lay == layer;
-};
 
-float pseudorapidity(float z, float r){
-    float theta = atan2(r, z);
-    return -log(tan(theta/2.0));
-}
-float pseudorapidity(float x, float y, float z){
-    float r = sqrt(x*x + y*y);
-    return pseudorapidity(z, r);
-}
-template<typename T>
-float pseudorapidity(const T &t){
-    return pseudorapidity(t.x, t.y, t.z);
-}
+void setup_hit_pair_functions(){
+
+    func<pair_vec(vector<HitPair>)>("calc_dphi_v_eta",
+        FUNC(([](const vector<HitPair>& hit_pairs){
+            vector<float> dphis;
+            vector<float> etas;
+            for(auto hit_pair : hit_pairs){
+                auto& pixrec_hit = std::get<PixRecHit>(hit_pair);
+                auto& sim_hit = std::get<SimHit>(hit_pair);
+                dphis.push_back(atan2(sim_hit.x, sim_hit.y) - atan2(pixrec_hit.x, pixrec_hit.y));
+                etas.push_back(pseudorapidity(pixrec_hit));
+            }
+            return std::make_pair(etas, dphis);
+        })));
 
+    func<pair_vec(vector<HitPair>)>("calc_dr_v_eta",
+        FUNC(([](const vector<HitPair>& hit_pairs){
+            vector<float> drs;
+            vector<float> etas;
+            for(auto hit_pair : hit_pairs){
+                auto& pixrec_hit = std::get<PixRecHit>(hit_pair);
+                auto& sim_hit = std::get<SimHit>(hit_pair);
+                drs.push_back(displacement(pixrec_hit, sim_hit));
+                etas.push_back(pseudorapidity(pixrec_hit));
+            }
+            return std::make_pair(etas, drs);
+        })));
+
+    func<pair_vec(vector<HitPair>)>("calc_dz_v_eta",
+        FUNC(([](const vector<HitPair>& hit_pairs){
+            vector<float> dzs;
+            vector<float> etas;
+            for(auto hit_pair : hit_pairs){
+                auto& pixrec_hit = std::get<PixRecHit>(hit_pair);
+                auto& sim_hit = std::get<SimHit>(hit_pair);
+                dzs.push_back(sim_hit.z - pixrec_hit.z);
+                etas.push_back(pseudorapidity(pixrec_hit));
+            }
+            return std::make_pair(etas, dzs);
+        })));
+
+    func<pair_vec(vector<HitPair>)>("calc_drho_v_eta",
+        FUNC(([](const vector<HitPair>& hit_pairs){
+            vector<float> drhos;
+            vector<float> etas;
+            for(auto hit_pair : hit_pairs){
+                auto& pixrec_hit = std::get<PixRecHit>(hit_pair);
+                auto& sim_hit = std::get<SimHit>(hit_pair);
+                drhos.push_back(rho(sim_hit) - rho(pixrec_hit));
+                etas.push_back(pseudorapidity(pixrec_hit));
+            }
+            return std::make_pair(etas, drhos);
+        })));
+}
 
 vector<string> seedTypes =
     {"initialStepSeeds",
@@ -71,6 +98,7 @@ vector<string> seedTypes =
      "pixelPairElectronSeeds",
      "stripPairElectronSeeds"};
 
+Value<vector<SuperCluster>>* super_clusters;
 Value<vector<Seed>>* seeds;
 Value<vector<PixRecHit>>* pixrec_hits;
 Value<vector<Track>>* tracks;
@@ -78,6 +106,7 @@ Value<vector<SimHit>>* sim_hits;
 Value<vector<SimTrack>>* sim_tracks;
 
 void register_objects(TrackingDataSet& tds){
+    super_clusters = register_super_clusters(tds);
     seeds = register_seeds(tds);
 
     pixrec_hits = register_pixrec_hits(tds);
@@ -85,12 +114,127 @@ void register_objects(TrackingDataSet& tds){
 
     sim_hits = register_sim_hits(tds);
     sim_tracks = register_sim_tracks(tds);
+
+}
+
+template<typename A>
+bool in_det(const A& hit, const unsigned int& det, const unsigned int& lay){
+    return hit.det == det && hit.lay == lay;
+}
+
+void setup_skipped_layer_hit_pairs(TrackingDataSet& tds){
+
+    /* Finds SimHit/RecHit pairs
+     */
+    auto find_matched_nth_hit_in_layer_with_skip =
+    func<vector<HitPair>(vector<Seed>,
+                         vector<PixRecHit>,
+                         vector<SimHit>,
+                         vector<Track>,
+                         vector<SimTrack>,
+                         int, int, int)>("find_matched_nth_hit_in_layer_with_skip",
+        FUNC(([](const vector<Seed>& seeds,
+                 const vector<PixRecHit>& pixrec_hits,
+                 const vector<SimHit>& sim_hits,
+                 const vector<Track>& tracks,
+                 const vector<SimTrack>& sim_tracks,
+                 const unsigned int&& det,
+                 const unsigned int&& pix_layer,
+                 const unsigned int&& skip){
+            vector<HitPair> matched_hits;
+            for(const Track &trk : tracks){ // loop over all tracks
+                const Seed &seed = seeds[trk.seedIdx];
+
+                // Require seed's source algorithm is one of those in seedTypes
+                if(seed.algoOriginal < 0 || seed.algoOriginal >= seedTypes.size()) continue;
+                // Require seed w/ at least 2 hits
+                if(seed.hitIdx.size() < 2) continue;
+                // take only pixel hits for now
+                if(seed.hitType[0] != HIT_TYPE_PIXEL || seed.hitType[1] != HIT_TYPE_PIXEL)   continue;
+
+                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]]});
+                    }
+                }
+            }
+            return matched_hits;
+        })));
+
+    auto pix_barrel = constant("PIXEL_BARREL", PIXEL_BARREL);
+    auto pix_endcap = constant("PIXEL_ENDCAP", PIXEL_ENDCAP);
+    auto L1 = constant("L1", 1);
+    auto L2 = constant("L2", 2);
+
+    auto skip_zero = constant("skip_zero", 0);
+    auto skip_one  = constant("skip_one", 1);
+
+
+    // 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,
+            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, pix_barrel, L1, skip_zero), "first_hits_in_B1_skip_zero");
+    auto first_hits_in_B2_skip_zero = fv::tup_apply(find_matched_nth_hit_in_layer_with_skip,
+            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, pix_barrel, L2, skip_zero), "first_hits_in_B2_skip_zero");
+
+
+    // First hit on 1st/2nd fpix layers and second hit in 2nd/3rd
+    auto first_hits_in_F1_skip_zero = fv::tup_apply(find_matched_nth_hit_in_layer_with_skip,
+            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, pix_endcap, L1, skip_zero), "first_hits_in_F1_skip_zero");
+    auto first_hits_in_F2_skip_zero = fv::tup_apply(find_matched_nth_hit_in_layer_with_skip,
+            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, pix_endcap, L2, skip_zero), "first_hits_in_F2_skip_zero");
+
+    // First hit on 1st/2nd fpix layers and second hit in 3rd/4th
+    auto first_hits_in_B1_skip_one = fv::tup_apply(find_matched_nth_hit_in_layer_with_skip,
+            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, pix_barrel, L1, skip_one), "first_hits_in_B1_skip_one");
+    auto first_hits_in_B2_skip_one = fv::tup_apply(find_matched_nth_hit_in_layer_with_skip,
+            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, pix_barrel, L2, skip_one), "first_hits_in_B2_skip_one");
+
+
+    // First hit on 1st fpix layer and second hit in 3rd layer
+    auto first_hits_in_F1_skip_one = fv::tup_apply(find_matched_nth_hit_in_layer_with_skip,
+            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, pix_endcap, L1, skip_one), "first_hits_in_F1_skip_one");
+
+
+    TH2Params params_dz = {"$\\eta$",       100, -4,   4,
+                           "$\\Delta z$(rad)", 50,  -.01, .01};
+    auto calc_dz_v_eta = lookup_function<pair_vec(vector<HitPair>)>("calc_dz_v_eta");
+    tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B1_skip_zero",
+                                                    fv::apply(calc_dz_v_eta, first_hits_in_B1_skip_zero),
+                                                    "First Hit in BPIX-L1 - Skip 0", params_dz);
+    tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B2_skip_zero",
+                                                    fv::apply(calc_dz_v_eta, first_hits_in_B2_skip_zero),
+                                                    "First Hit in BPIX-L2 - Skip 0", params_dz);
+
+    tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_F1_skip_zero",
+                                                    fv::apply(calc_dz_v_eta, first_hits_in_F1_skip_zero),
+                                                    "First Hit in FPIX-L1 - Skip 0", params_dz);
+    tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_F2_skip_zero",
+                                                    fv::apply(calc_dz_v_eta, first_hits_in_F2_skip_zero),
+                                                    "First Hit in FPIX-L2 - Skip 0", params_dz);
+
+    tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B1_skip_one",
+                                                    fv::apply(calc_dz_v_eta, first_hits_in_B1_skip_one),
+                                                    "First Hit in BPIX-L1 - Skip 1", params_dz);
+    tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B2_skip_one",
+                                                    fv::apply(calc_dz_v_eta, first_hits_in_B2_skip_one),
+                                                    "First Hit in BPIX-L2 - Skip 1", params_dz);
+
+    tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_F1_skip_one",
+                                                    fv::apply(calc_dz_v_eta, first_hits_in_F1_skip_one),
+                                                    "First Hit in FPIX-L1 - Skip 1", params_dz);
+
 }
 
 void setup_first_hit_pairs(TrackingDataSet& tds){
-    typedef std::tuple<PixRecHit, SimHit> HitPair;
 
-    auto& find_matched_nth_hit_in_layer =
+    auto find_matched_nth_hit_in_layer =
     func<vector<HitPair>(vector<Seed>,
                          vector<PixRecHit>,
                          vector<SimHit>,
@@ -102,24 +246,24 @@ void setup_first_hit_pairs(TrackingDataSet& tds){
                  const vector<SimHit>& sim_hits,
                  const vector<Track>& tracks,
                  const vector<SimTrack>& sim_tracks,
-                 const int&& det,
-                 const int&& pix_layer,
-                 const int&& hit_number){
+                 const unsigned int&& det,
+                 const unsigned int&& pix_layer,
+                 const unsigned int&& hit_number){
             vector<HitPair> matched_hits;
             for(const Track &trk : tracks){ // loop over all tracks
                 const Seed &seed = seeds[trk.seedIdx];
-                if(seed.hitIdx.size() <= hit_number) continue; // looking for hit_number'th hit, which this seed doesn't have
+                if(seed.hitIdx.size() <= hit_number) continue; // looking for hit_number'th hit
                 if(seed.algoOriginal < 0 || seed.algoOriginal >= seedTypes.size()) continue;
                 if(seed.hitType[hit_number] != HIT_TYPE_PIXEL) continue; // take only pixel hits for now
                 const PixRecHit &rec_hit = pixrec_hits[seed.hitIdx[hit_number]];
                 if(rec_hit.det == det && rec_hit.lay == pix_layer){
                     // We have the RecHit we want to work with, now find a properly* matched SimHit
                     if(rec_hit.simHitIdx.size() == 0) continue; // if rechit is matched to no simhits, give up.
-                    for(const int &simTrkIdx : trk.simTrkIdx){  // loop over SimTracks matched to Track
-                        for(const int& simHitIdx : sim_tracks[simTrkIdx].simHitIdx){ // loop over SimHits from SimTrack
-                            if(simHitIdx == rec_hit.simHitIdx[0]) // take first matched simhit (should be the closest one)
-                                matched_hits.push_back({rec_hit, sim_hits[rec_hit.simHitIdx[0]]});
-                        }
+                    if(trk.simTrkIdx.size() == 0) continue; // if rechit is matched to no simhits, 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_hit.simHitIdx[0]) // take first matched simhit (should be the closest one)
+                            matched_hits.push_back({rec_hit, sim_hits[rec_hit.simHitIdx[0]]});
                     }
                 }
             }
@@ -132,34 +276,34 @@ void setup_first_hit_pairs(TrackingDataSet& tds){
     auto second_hit = constant("2nd", 1);
 
     // First hits on inner three bpix layers
-    auto first_hits_in_B1 = fv::apply(find_matched_nth_hit_in_layer,
+    auto first_hits_in_B1 = fv::tup_apply(find_matched_nth_hit_in_layer,
             fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, barrel_val, constant("L1", 1), first_hit), "first_hits_in_B1");
-    auto first_hits_in_B2 = fv::apply(find_matched_nth_hit_in_layer,
+    auto first_hits_in_B2 = fv::tup_apply(find_matched_nth_hit_in_layer,
             fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, barrel_val, constant("L2", 2), first_hit), "first_hits_in_B2");
 
     // Second hits on outer three bpix layers
-    auto second_hits_in_B2 = fv::apply(find_matched_nth_hit_in_layer,
+    auto second_hits_in_B2 = fv::tup_apply(find_matched_nth_hit_in_layer,
             fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, barrel_val, constant("L2", 2), second_hit), "second_hits_in_B2");
-    auto second_hits_in_B3 = fv::apply(find_matched_nth_hit_in_layer,
+    auto second_hits_in_B3 = fv::tup_apply(find_matched_nth_hit_in_layer,
             fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, barrel_val, constant("L3", 3), second_hit), "second_hits_in_B3");
-    auto second_hits_in_B4 = fv::apply(find_matched_nth_hit_in_layer,
+    auto second_hits_in_B4 = fv::tup_apply(find_matched_nth_hit_in_layer,
             fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, barrel_val, constant("L4", 4), second_hit), "second_hits_in_B4");
 
-    auto first_hits_in_F1 = fv::apply(find_matched_nth_hit_in_layer,
+    auto first_hits_in_F1 = fv::tup_apply(find_matched_nth_hit_in_layer,
             fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, endcap_val, constant("L1", 1), first_hit), "first_hits_in_F1");
-    auto first_hits_in_F2 = fv::apply(find_matched_nth_hit_in_layer,
+    auto first_hits_in_F2 = fv::tup_apply(find_matched_nth_hit_in_layer,
             fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, endcap_val, constant("L2", 2), first_hit), "first_hits_in_F2");
 
-    auto second_hits_in_F2 = fv::apply(find_matched_nth_hit_in_layer,
+    auto second_hits_in_F2 = fv::tup_apply(find_matched_nth_hit_in_layer,
             fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, endcap_val, constant("L2", 2), second_hit), "second_hits_in_F2");
-    auto second_hits_in_F3 = fv::apply(find_matched_nth_hit_in_layer,
+    auto second_hits_in_F3 = fv::tup_apply(find_matched_nth_hit_in_layer,
             fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, endcap_val, constant("L3", 3), second_hit), "second_hits_in_F3");
 
     // Even vs Odd Ladders
     auto even = constant("even", false);
     auto odd = constant("odd", true);
 
-    auto& select_even_odd_ladder_blade_hit_pairs =
+    auto select_even_odd_ladder_blade_hit_pairs =
     func<vector<HitPair>(vector<HitPair>, bool)>("select_even_odd_ladder_blade_hit_pairs",
         FUNC(([](const vector<HitPair>& hit_pairs,
                  const bool &&odd){
@@ -172,210 +316,292 @@ void setup_first_hit_pairs(TrackingDataSet& tds){
             return even_pairs;
         })));
 
-    auto first_hits_in_B1_even_ladder = fv::apply(select_even_odd_ladder_blade_hit_pairs,
+    auto first_hits_in_B1_even_ladder = fv::tup_apply(select_even_odd_ladder_blade_hit_pairs,
             fv::tuple(first_hits_in_B1, even), "first_hits_in_B1_even_ladder");
-    auto first_hits_in_B1_odd_ladder = fv::apply(select_even_odd_ladder_blade_hit_pairs,
+    auto first_hits_in_B1_odd_ladder = fv::tup_apply(select_even_odd_ladder_blade_hit_pairs,
             fv::tuple(first_hits_in_B1, odd), "first_hits_in_B1_odd_ladder");
 
-    auto first_hits_in_B2_even_ladder = fv::apply(select_even_odd_ladder_blade_hit_pairs,
+    auto first_hits_in_B2_even_ladder = fv::tup_apply(select_even_odd_ladder_blade_hit_pairs,
             fv::tuple(first_hits_in_B2, even), "first_hits_in_B2_even_ladder");
-    auto first_hits_in_B2_odd_ladder = fv::apply(select_even_odd_ladder_blade_hit_pairs,
+    auto first_hits_in_B2_odd_ladder = fv::tup_apply(select_even_odd_ladder_blade_hit_pairs,
             fv::tuple(first_hits_in_B2, odd), "first_hits_in_B2_odd_ladder");
 
     //Plots for dPhi of collections defined above
-    auto& calc_dphi_v_eta = func<pair_vec(vector<HitPair>)>("calc_dphi_v_eta",
-        FUNC(([](const vector<HitPair>& hit_pairs){
-            vector<float> dphis;
-            vector<float> etas;
-            for(auto hit_pair : hit_pairs){
-                auto& pixrec_hit = std::get<PixRecHit>(hit_pair);
-                auto& sim_hit = std::get<SimHit>(hit_pair);
-                dphis.push_back(atan2(sim_hit.x, sim_hit.y) - atan2(pixrec_hit.x, pixrec_hit.y));
-                etas.push_back(pseudorapidity(pixrec_hit));
-            }
-            return std::make_pair(etas, dphis);
-        })));
+    auto calc_dphi_v_eta = lookup_function<pair_vec(vector<HitPair>)>("calc_dphi_v_eta");
     TH2Params params_dphi = {"$\\eta$",       100, -4,   4,
                              "$\\Delta \\phi$(rad)", 50,  -.0015, .0015};
     tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_B1",
-                                                    fv::apply(calc_dphi_v_eta, fv::tuple(first_hits_in_B1)),
+                                                    fv::apply(calc_dphi_v_eta, first_hits_in_B1),
                                                     "First Hit in BPIX-L1", params_dphi);
     tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_B2",
-                                                    fv::apply(calc_dphi_v_eta, fv::tuple(first_hits_in_B2)),
+                                                    fv::apply(calc_dphi_v_eta, first_hits_in_B2),
                                                     "First Hit in BPIX-L2", params_dphi);
 
     tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_B1_even_ladder",
-                                                    fv::apply(calc_dphi_v_eta, fv::tuple(first_hits_in_B1_even_ladder)),
+                                                    fv::apply(calc_dphi_v_eta, first_hits_in_B1_even_ladder),
                                                     "First Hit in BPIX-L1 - Even Ladders", params_dphi);
     tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_B1_odd_ladder",
-                                                    fv::apply(calc_dphi_v_eta, fv::tuple(first_hits_in_B1_odd_ladder)),
+                                                    fv::apply(calc_dphi_v_eta, first_hits_in_B1_odd_ladder),
                                                     "First Hit in BPIX-L1 - Odd Ladders", params_dphi);
     tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_B2_even_ladder",
-                                                    fv::apply(calc_dphi_v_eta, fv::tuple(first_hits_in_B2_even_ladder)),
+                                                    fv::apply(calc_dphi_v_eta, first_hits_in_B2_even_ladder),
                                                     "First Hit in BPIX-L2 - Even Ladders", params_dphi);
     tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_B2_odd_ladder",
-                                                    fv::apply(calc_dphi_v_eta, fv::tuple(first_hits_in_B2_odd_ladder)),
+                                                    fv::apply(calc_dphi_v_eta, first_hits_in_B2_odd_ladder),
                                                     "First Hit in BPIX-L2 - Odd Ladders", params_dphi);
 
     tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_second_hits_in_B2",
-                                                    fv::apply(calc_dphi_v_eta, fv::tuple(second_hits_in_B2)),
+                                                    fv::apply(calc_dphi_v_eta, second_hits_in_B2),
                                                     "Second Hit in BPIX-L2", params_dphi);
     tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_second_hits_in_B3",
-                                                    fv::apply(calc_dphi_v_eta, fv::tuple(second_hits_in_B3)),
+                                                    fv::apply(calc_dphi_v_eta, second_hits_in_B3),
                                                     "Second Hit in BPIX-L3", params_dphi);
     tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_second_hits_in_B4",
-                                                    fv::apply(calc_dphi_v_eta, fv::tuple(second_hits_in_B4)),
+                                                    fv::apply(calc_dphi_v_eta, second_hits_in_B4),
                                                     "Second Hit in BPIX-L4", params_dphi);
 
     tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_F1",
-                                                    fv::apply(calc_dphi_v_eta, fv::tuple(first_hits_in_F1)),
+                                                    fv::apply(calc_dphi_v_eta, first_hits_in_F1),
                                                     "First Hit in FPIX-L1", params_dphi);
     tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_F2",
-                                                    fv::apply(calc_dphi_v_eta, fv::tuple(first_hits_in_F2)),
+                                                    fv::apply(calc_dphi_v_eta, first_hits_in_F2),
                                                     "First Hit in FPIX-L2", params_dphi);
 
     tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_second_hits_in_F2",
-                                                    fv::apply(calc_dphi_v_eta, fv::tuple(second_hits_in_F2)),
+                                                    fv::apply(calc_dphi_v_eta, second_hits_in_F2),
                                                     "Second Hit in FPIX-L2", params_dphi);
     tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_second_hits_in_F3",
-                                                    fv::apply(calc_dphi_v_eta, fv::tuple(second_hits_in_F3)),
+                                                    fv::apply(calc_dphi_v_eta, second_hits_in_F3),
                                                     "Second Hit in FPIX-L3", params_dphi);
 
     //Plots for dr of collections defined above
-    auto& calc_dr_v_eta = func<pair_vec(vector<HitPair>)>("calc_dr_v_eta",
-        FUNC(([](const vector<HitPair>& hit_pairs){
-            vector<float> drs;
-            vector<float> etas;
-            for(auto hit_pair : hit_pairs){
-                auto& pixrec_hit = std::get<PixRecHit>(hit_pair);
-                auto& sim_hit = std::get<SimHit>(hit_pair);
-                drs.push_back(displacement(pixrec_hit, sim_hit));
-                etas.push_back(pseudorapidity(pixrec_hit));
-            }
-            return std::make_pair(etas, drs);
-        })));
+    auto calc_dr_v_eta = lookup_function<pair_vec(vector<HitPair>)>("calc_dr_v_eta");
     TH2Params params_dr = {"$\\eta$",       100, -4,   4,
                            "$\\Delta \\r$(rad)", 50, 0, .006};
     tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_B1",
-                                                    fv::apply(calc_dr_v_eta, fv::tuple(first_hits_in_B1)),
+                                                    fv::apply(calc_dr_v_eta, first_hits_in_B1),
                                                     "First Hit in BPIX-L1", params_dr);
     tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_B2",
-                                                    fv::apply(calc_dr_v_eta, fv::tuple(first_hits_in_B2)),
+                                                    fv::apply(calc_dr_v_eta, first_hits_in_B2),
                                                     "First Hit in BPIX-L2", params_dr);
 
     tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_B1_even_ladder",
-                                                    fv::apply(calc_dr_v_eta, fv::tuple(first_hits_in_B1_even_ladder)),
+                                                    fv::apply(calc_dr_v_eta, first_hits_in_B1_even_ladder),
                                                     "First Hit in BPIX-L1 - Even Ladders", params_dr);
     tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_B1_odd_ladder",
-                                                    fv::apply(calc_dr_v_eta, fv::tuple(first_hits_in_B1_odd_ladder)),
+                                                    fv::apply(calc_dr_v_eta, first_hits_in_B1_odd_ladder),
                                                     "First Hit in BPIX-L1 - Odd Ladders", params_dr);
     tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_B2_even_ladder",
-                                                    fv::apply(calc_dr_v_eta, fv::tuple(first_hits_in_B2_even_ladder)),
+                                                    fv::apply(calc_dr_v_eta, first_hits_in_B2_even_ladder),
                                                     "First Hit in BPIX-L2 - Even Ladders", params_dr);
     tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_B2_odd_ladder",
-                                                    fv::apply(calc_dr_v_eta, fv::tuple(first_hits_in_B2_odd_ladder)),
+                                                    fv::apply(calc_dr_v_eta, first_hits_in_B2_odd_ladder),
                                                     "First Hit in BPIX-L2 - Odd Ladders", params_dr);
 
     tds.register_container<ContainerTH2Many<float>>("dr_v_eta_second_hits_in_B2",
-                                                    fv::apply(calc_dr_v_eta, fv::tuple(second_hits_in_B2)),
+                                                    fv::apply(calc_dr_v_eta, second_hits_in_B2),
                                                     "Second Hit in BPIX-L2", params_dr);
     tds.register_container<ContainerTH2Many<float>>("dr_v_eta_second_hits_in_B3",
-                                                    fv::apply(calc_dr_v_eta, fv::tuple(second_hits_in_B3)),
+                                                    fv::apply(calc_dr_v_eta, second_hits_in_B3),
                                                     "Second Hit in BPIX-L3", params_dr);
     tds.register_container<ContainerTH2Many<float>>("dr_v_eta_second_hits_in_B4",
-                                                    fv::apply(calc_dr_v_eta, fv::tuple(second_hits_in_B4)),
+                                                    fv::apply(calc_dr_v_eta, second_hits_in_B4),
                                                     "Second Hit in BPIX-L4", params_dr);
 
     tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_F1",
-                                                    fv::apply(calc_dr_v_eta, fv::tuple(first_hits_in_F1)),
+                                                    fv::apply(calc_dr_v_eta, first_hits_in_F1),
                                                     "First Hit in FPIX-L1", params_dr);
     tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_F2",
-                                                    fv::apply(calc_dr_v_eta, fv::tuple(first_hits_in_F2)),
+                                                    fv::apply(calc_dr_v_eta, first_hits_in_F2),
                                                     "First Hit in FPIX-L2", params_dr);
 
     tds.register_container<ContainerTH2Many<float>>("dr_v_eta_second_hits_in_F2",
-                                                    fv::apply(calc_dr_v_eta, fv::tuple(second_hits_in_F2)),
+                                                    fv::apply(calc_dr_v_eta, second_hits_in_F2),
                                                     "Second Hit in FPIX-L2", params_dr);
     tds.register_container<ContainerTH2Many<float>>("dr_v_eta_second_hits_in_F3",
-                                                    fv::apply(calc_dr_v_eta, fv::tuple(second_hits_in_F3)),
+                                                    fv::apply(calc_dr_v_eta, second_hits_in_F3),
                                                     "Second Hit in FPIX-L3", params_dr);
 
 
     //Plots for dz of collections defined above
-    auto& calc_dz_v_eta = func<pair_vec(vector<HitPair>)>("calc_dz_v_eta",
-        FUNC(([](const vector<HitPair>& hit_pairs){
-            vector<float> dzs;
-            vector<float> etas;
-            for(auto hit_pair : hit_pairs){
-                auto& pixrec_hit = std::get<PixRecHit>(hit_pair);
-                auto& sim_hit = std::get<SimHit>(hit_pair);
-                dzs.push_back(sim_hit.z - pixrec_hit.z);
-                etas.push_back(pseudorapidity(pixrec_hit));
-            }
-            return std::make_pair(etas, dzs);
-        })));
+    auto calc_dz_v_eta = lookup_function<pair_vec(vector<HitPair>)>("calc_dz_v_eta");
     TH2Params params_dz = {"$\\eta$",       100, -4,   4,
                            "$\\Delta z$(rad)", 50,  -.01, .01};
     tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B1",
-                                                    fv::apply(calc_dz_v_eta, fv::tuple(first_hits_in_B1)),
+                                                    fv::apply(calc_dz_v_eta, first_hits_in_B1),
                                                     "First Hit in BPIX-L1", params_dz);
     tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B2",
-                                                    fv::apply(calc_dz_v_eta, fv::tuple(first_hits_in_B2)),
+                                                    fv::apply(calc_dz_v_eta, first_hits_in_B2),
                                                     "First Hit in BPIX-L2", params_dz);
 
     tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B1_even_ladder",
-                                                    fv::apply(calc_dz_v_eta, fv::tuple(first_hits_in_B1_even_ladder)),
+                                                    fv::apply(calc_dz_v_eta, first_hits_in_B1_even_ladder),
                                                     "First Hit in BPIX-L1 - Even Ladders", params_dz);
     tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B1_odd_ladder",
-                                                    fv::apply(calc_dz_v_eta, fv::tuple(first_hits_in_B1_odd_ladder)),
+                                                    fv::apply(calc_dz_v_eta, first_hits_in_B1_odd_ladder),
                                                     "First Hit in BPIX-L1 - Odd Ladders", params_dz);
     tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B2_even_ladder",
-                                                    fv::apply(calc_dz_v_eta, fv::tuple(first_hits_in_B2_even_ladder)),
+                                                    fv::apply(calc_dz_v_eta, first_hits_in_B2_even_ladder),
                                                     "First Hit in BPIX-L2 - Even Ladders", params_dz);
     tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B2_odd_ladder",
-                                                    fv::apply(calc_dz_v_eta, fv::tuple(first_hits_in_B2_odd_ladder)),
+                                                    fv::apply(calc_dz_v_eta, first_hits_in_B2_odd_ladder),
                                                     "First Hit in BPIX-L2 - Odd Ladders", params_dz);
 
     tds.register_container<ContainerTH2Many<float>>("dz_v_eta_second_hits_in_B2",
-                                                    fv::apply(calc_dz_v_eta, fv::tuple(second_hits_in_B2)),
+                                                    fv::apply(calc_dz_v_eta, second_hits_in_B2),
                                                     "Second Hit in BPIX-L2", params_dz);
     tds.register_container<ContainerTH2Many<float>>("dz_v_eta_second_hits_in_B3",
-                                                    fv::apply(calc_dz_v_eta, fv::tuple(second_hits_in_B3)),
+                                                    fv::apply(calc_dz_v_eta, second_hits_in_B3),
                                                     "Second Hit in BPIX-L3", params_dz);
     tds.register_container<ContainerTH2Many<float>>("dz_v_eta_second_hits_in_B4",
-                                                    fv::apply(calc_dz_v_eta, fv::tuple(second_hits_in_B4)),
+                                                    fv::apply(calc_dz_v_eta, second_hits_in_B4),
                                                     "Second Hit in BPIX-L4", params_dz);
 
     //Plots for drho of collections defined above
-    auto& calc_drho_v_eta = func<pair_vec(vector<HitPair>)>("calc_drho_v_eta",
-        FUNC(([](const vector<HitPair>& hit_pairs){
-            vector<float> drhos;
-            vector<float> etas;
-            for(auto hit_pair : hit_pairs){
-                auto& pixrec_hit = std::get<PixRecHit>(hit_pair);
-                auto& sim_hit = std::get<SimHit>(hit_pair);
-                drhos.push_back(rho(sim_hit) - rho(pixrec_hit));
-                etas.push_back(pseudorapidity(pixrec_hit));
-            }
-            return std::make_pair(etas, drhos);
-        })));
+    auto calc_drho_v_eta = lookup_function<pair_vec(vector<HitPair>)>("calc_drho_v_eta");
     TH2Params params_drho = {"$\\eta$",       100, -4,   4,
                              "$\\Delta \\rho$(cm)", 50,  -.01, .01};
     tds.register_container<ContainerTH2Many<float>>("drho_v_eta_first_hits_in_F1",
-                                                    fv::apply(calc_drho_v_eta, fv::tuple(first_hits_in_F1)),
+                                                    fv::apply(calc_drho_v_eta, first_hits_in_F1),
                                                     "First Hit in FPIX-L1", params_drho);
     tds.register_container<ContainerTH2Many<float>>("drho_v_eta_first_hits_in_F2",
-                                                    fv::apply(calc_drho_v_eta, fv::tuple(first_hits_in_F2)),
+                                                    fv::apply(calc_drho_v_eta, first_hits_in_F2),
                                                     "First Hit in FPIX-L2", params_drho);
 
     tds.register_container<ContainerTH2Many<float>>("drho_v_eta_second_hits_in_F2",
-                                                    fv::apply(calc_drho_v_eta, fv::tuple(second_hits_in_F2)),
+                                                    fv::apply(calc_drho_v_eta, second_hits_in_F2),
                                                     "Second Hit in FPIX-L2", params_drho);
     tds.register_container<ContainerTH2Many<float>>("drho_v_eta_second_hits_in_F3",
-                                                    fv::apply(calc_drho_v_eta, fv::tuple(second_hits_in_F3)),
+                                                    fv::apply(calc_drho_v_eta, second_hits_in_F3),
                                                     "Second Hit in FPIX-L3", params_drho);
 }
 
+
+void setup_sc_residuals(TrackingDataSet& tds){
+
+    auto sc_hit_residuals =
+    func<pair_vec(vector<SuperCluster>,
+                       int, int, int, bool)>("sc_hit_residuals",
+        FUNC(([](const vector<SuperCluster>& super_clusters,
+                 const int&& det,
+                 const int&& pix_layer,
+                 const int&& hit_number,
+                 const bool&& sel_dRz){
+            vector<float> residuals;
+            vector<float> etas;
+            for(const SuperCluster& super_cluster : super_clusters){ // loop over all supser-clusters
+                float dRz, dPhi;
+                unsigned int nSeeds = super_cluster.charge.size();
+                for(unsigned int seedIdx=0; seedIdx<nSeeds; seedIdx++){ // loop over all seeds associated w/ sc
+                    if(hit_number == 0){
+                        if(super_cluster.lay1[seedIdx] != pix_layer || super_cluster.subDet1[seedIdx] != det) continue;
+                        dRz = super_cluster.dRz1[seedIdx];
+                        dPhi = super_cluster.dPhi1[seedIdx];
+                    }else{
+                        if(super_cluster.lay2[seedIdx] != pix_layer || super_cluster.subDet2[seedIdx] != det) continue;
+                        dRz = super_cluster.dRz2[seedIdx];
+                        dPhi = super_cluster.dPhi2[seedIdx];
+                    }
+                    DEBUG("New Seed, Idx: " << seedIdx << endl
+                          << "\tdRz: " << dRz << endl
+                          << "\tdPhi: " << dPhi << endl);
+                    if(sel_dRz)
+                        residuals.push_back(dRz);
+                    else
+                        residuals.push_back(dPhi);
+                    etas.push_back(pseudorapidityP(super_cluster));
+                }
+            }
+            return std::make_pair(etas,residuals);
+        })));
+
+
+    auto pix_barrel = constant("PIXEL_BARREL", PIXEL_BARREL);
+    /* auto pix_endcap = constant<unsigned int>("PIXEL_ENDCAP", PIXEL_ENDCAP); */
+    auto first_hit  = constant("1st", 0);
+    auto second_hit = constant("2nd", 1);
+    auto L1 = constant("L1", 1);
+    auto L2 = constant("L2", 2);
+    auto L3 = constant("L3", 3);
+    auto L4 = constant("L4", 4);
+
+    auto dRz = constant("dRz", true);
+    auto dPhi = constant("dPhi", false);
+
+    // First hits on inner two bpix layers
+    auto sc_first_hits_in_B1_dz = fv::tup_apply(sc_hit_residuals,
+            fv::tuple(super_clusters, pix_barrel, L1, first_hit, dRz), "sc_first_hits_in_B1_dz");
+    auto sc_first_hits_in_B2_dz = fv::tup_apply(sc_hit_residuals,
+            fv::tuple(super_clusters, pix_barrel, L2, first_hit, dRz), "sc_first_hits_in_B2_dz");
+    auto sc_first_hits_in_B3_dz = fv::tup_apply(sc_hit_residuals,
+            fv::tuple(super_clusters, pix_barrel, L3, first_hit, dRz), "sc_first_hits_in_B3_dz");
+
+    auto sc_second_hits_in_B2_dz = fv::tup_apply(sc_hit_residuals,
+            fv::tuple(super_clusters, pix_barrel, L2, second_hit, dRz), "sc_second_hits_in_B2_dz");
+    auto sc_second_hits_in_B3_dz = fv::tup_apply(sc_hit_residuals,
+            fv::tuple(super_clusters, pix_barrel, L3, second_hit, dRz), "sc_second_hits_in_B3_dz");
+    auto sc_second_hits_in_B4_dz = fv::tup_apply(sc_hit_residuals,
+            fv::tuple(super_clusters, pix_barrel, L4, second_hit, dRz), "sc_second_hits_in_B4_dz");
+
+    auto sc_first_hits_in_B1_dphi = fv::tup_apply(sc_hit_residuals,
+            fv::tuple(super_clusters, pix_barrel, L1, first_hit, dPhi), "sc_first_hits_in_B1_dphi");
+    auto sc_first_hits_in_B2_dphi = fv::tup_apply(sc_hit_residuals,
+            fv::tuple(super_clusters, pix_barrel, L2, first_hit, dPhi), "sc_first_hits_in_B2_dphi");
+    auto sc_first_hits_in_B3_dphi = fv::tup_apply(sc_hit_residuals,
+            fv::tuple(super_clusters, pix_barrel, L3, first_hit, dPhi), "sc_first_hits_in_B3_dphi");
+
+    auto sc_second_hits_in_B2_dphi = fv::tup_apply(sc_hit_residuals,
+            fv::tuple(super_clusters, pix_barrel, L2, second_hit, dPhi), "sc_second_hits_in_B2_dphi");
+    auto sc_second_hits_in_B3_dphi = fv::tup_apply(sc_hit_residuals,
+            fv::tuple(super_clusters, pix_barrel, L3, second_hit, dPhi), "sc_second_hits_in_B3_dphi");
+    auto sc_second_hits_in_B4_dphi = fv::tup_apply(sc_hit_residuals,
+            fv::tuple(super_clusters, pix_barrel, L4, second_hit, dPhi), "sc_second_hits_in_B4_dphi");
+
+    TH2Params params_dz = {"$\\eta$",       100, -4,   4,
+                           "$\\delta z$(cm)", 100, -0.7,  0.7};
+    tds.register_container<ContainerTH2Many<float>>("sc_first_hits_in_B1_dz",
+                                                    sc_first_hits_in_B1_dz,
+                                                    "First Hit in BPIX-L1", params_dz);
+    tds.register_container<ContainerTH2Many<float>>("sc_first_hits_in_B2_dz",
+                                                    sc_first_hits_in_B2_dz,
+                                                    "First Hit in BPIX-L2", params_dz);
+    tds.register_container<ContainerTH2Many<float>>("sc_first_hits_in_B3_dz",
+                                                    sc_first_hits_in_B3_dz,
+                                                    "First Hit in BPIX-L3", params_dz);
+
+    tds.register_container<ContainerTH2Many<float>>("sc_second_hits_in_B2_dz",
+                                                    sc_second_hits_in_B2_dz,
+                                                    "Second Hit in BPIX-L2", params_dz);
+    tds.register_container<ContainerTH2Many<float>>("sc_second_hits_in_B3_dz",
+                                                    sc_second_hits_in_B3_dz,
+                                                    "Second Hit in BPIX-L3", params_dz);
+    tds.register_container<ContainerTH2Many<float>>("sc_second_hits_in_B4_dz",
+                                                    sc_second_hits_in_B4_dz,
+                                                    "Second Hit in BPIX-L4", params_dz);
+
+    TH2Params params_dphi = {"$\\eta$",       100, -4,   4,
+                             "$\\delta \\phi$(rad)", 500, -0.7,  0.7};
+    tds.register_container<ContainerTH2Many<float>>("sc_first_hits_in_B1_dphi",
+                                                    sc_first_hits_in_B1_dphi,
+                                                    "First Hit in BPIX-L1", params_dphi);
+    tds.register_container<ContainerTH2Many<float>>("sc_first_hits_in_B2_dphi",
+                                                    sc_first_hits_in_B2_dphi,
+                                                    "First Hit in BPIX-L2", params_dphi);
+    tds.register_container<ContainerTH2Many<float>>("sc_first_hits_in_B3_dphi",
+                                                    sc_first_hits_in_B3_dphi,
+                                                    "First Hit in BPIX-L3", params_dphi);
+
+    tds.register_container<ContainerTH2Many<float>>("sc_second_hits_in_B2_dphi",
+                                                    sc_second_hits_in_B2_dphi,
+                                                    "Second Hit in BPIX-L2", params_dphi);
+    tds.register_container<ContainerTH2Many<float>>("sc_second_hits_in_B3_dphi",
+                                                    sc_second_hits_in_B3_dphi,
+                                                    "Second Hit in BPIX-L3", params_dphi);
+    tds.register_container<ContainerTH2Many<float>>("sc_second_hits_in_B4_dphi",
+                                                    sc_second_hits_in_B4_dphi,
+                                                    "Second Hit in BPIX-L4", params_dphi);
+}
+
 void run_analysis(const vector<fv::util::DataFileDescriptor>& dfds, const string output_filename, bool silent){
     gSystem->Load("libfilval.so");
     auto replace_suffix = [](const std::string& input, const std::string& new_suffix){
@@ -388,8 +614,11 @@ void run_analysis(const vector<fv::util::DataFileDescriptor>& dfds, const string
     TrackingDataSet tds(output_filename, dfds, "trackingNtuple/tree");
     /* tds.set_max_events(10); */
     register_objects(tds);
+    setup_hit_pair_functions();
 
     setup_first_hit_pairs(tds);
+    setup_skipped_layer_hit_pairs(tds);
+    setup_sc_residuals(tds);
 
 
     tds.process(silent);
@@ -406,11 +635,12 @@ int main(int argc, char * argv[]){
     }else if(args.cmdOptionExists("-f")){
         string input_filename = args.getCmdOption("-f");
         string data_label = args.cmdOptionExists("-l") ? args.getCmdOption("-l") : "";
-        fv::util::DataFileDescriptor dfd(input_filename, data_label);
+        string data_category = args.cmdOptionExists("-c") ? args.getCmdOption("-c") : "";
+        fv::util::DataFileDescriptor dfd(input_filename, data_label, data_category);
         run_analysis(std::vector<fv::util::DataFileDescriptor>({dfd}), output_filename, silent);
     }else {
         cout << "Usage: ./tracking_validation (-s) {-o output_filename} -F datafiles.txt" << endl;
-        cout << "       ./tracking_validation (-s) {-l DATA_LABEL} {-o output_filename} -f treefile.root" << endl;
+        cout << "       ./tracking_validation (-s) {-l DATA_LABEL} {-c DATA_CATEGORY} {-o output_filename} -f treefile.root" << endl;
     }
     return 0;
 }

filval

@@ -1 +1 @@
-Subproject commit 9f48da0b66d494cdec8bfadc49359f349e99c36b
+Subproject commit 1a85d61848cbde369c479b434883640cfd08e3f0