EGamma_ElectronTrackingValidation/analysis/tracking_validation.cpp

#include <iostream>
#include <vector>
#include <map>
#include <utility>
#include <numeric>
#include <limits>
#include <TSystem.h>

#include "filval/filval.hpp"
#include "filval/root/filval.hpp"

#include <boost/range/combine.hpp>
#include <boost/format.hpp>

#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::pair<std::vector<float>,std::vector<float>> pair_vec;

#define HIT_TYPE_PIXEL 0
#define HIT_TYPE_GLUED 1
#define HIT_TYPE_STRIP 2

enum PixSubDet{
    Void = 0,
    Barrel = 1,
    EndCap = 2
};

struct PixDetElem{
    PixSubDet subdet;
    unsigned char lay;

    PixDetElem():subdet(PixSubDet::Void), lay(0) {}

    PixDetElem(const PixSubDet& subdet, const unsigned char& lay){
        this->subdet = subdet;
        this->lay = lay;
    }

    template<typename T>
    bool contains(const T& obj) const{
        return obj.det == subdet && obj.lay == lay;
    }

    bool contains(const PixSubDet& subdet, const unsigned char& lay) const{
        return subdet == this->subdet && lay == this->lay;
    }
};

std::map<PixSubDet,string> subdet_names = {{PixSubDet::Barrel, "BPIX"}, {PixSubDet::EndCap, "FPIX"}};


typedef std::tuple<PixRecHit, SimHit> HitPair;

vector<string> seedTypes =
    {"initialStepSeeds",
     "highPtTripletStepSeeds",
     "mixedTripletStepSeeds",
     "pixelLessStepSeeds",
     "tripletElectronSeeds",
     "pixelPairElectronSeeds",
     "stripPairElectronSeeds"};

Value<vector<SuperCluster>>* super_clusters;
Value<vector<Seed>>* seeds;
Value<vector<PixRecHit>>* pixrec_hits;
Value<vector<Track>>* tracks;
Value<vector<SimHit>>* sim_hits;
Value<vector<SimTrack>>* sim_tracks;

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);
        })));
}

void register_objects(TrackingDataSet& tds){
    super_clusters = register_super_clusters(tds);
    seeds = register_seeds(tds);

    pixrec_hits = register_pixrec_hits(tds);
    tracks = register_tracks(tds);

    sim_hits = register_sim_hits(tds);
    sim_tracks = register_sim_tracks(tds);

}

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>,
                         PixDetElem, int, bool)>("find_matched_nth_hit_in_layer_with_skip",
        FUNC(([](const vector<Seed>& seeds,
                 const vector<PixRecHit>& pixrec_hits,
                 const vector<SimHit>& sim_hits,
                 const vector<Track>& tracks,
                 const vector<SimTrack>& sim_tracks,
                 const PixDetElem&& pixdet,
                 const unsigned int&& skip,
                 const bool&& first){
            vector<HitPair> matched_hits;
            for(const Track &trk : tracks){ // loop over all tracks
                const Seed &seed = seeds[trk.seedIdx];

                // 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(first){ // Looking at first hit in the pair
                    if(pixdet.contains(rec_hit1) && pixdet.contains((PixSubDet)rec_hit2.det, rec_hit2.lay-1-skip)){
                        // We have the RecHit we want to work with, now find a properly* matched SimHit
                        if(rec_hit1.simHitIdx.size() == 0) continue; // if rechit is matched to no simhits, give up.
                        if(trk.simTrkIdx.size() == 0) continue; // if track is matched to no simtracks, give up.
                        const int &simTrkIdx = trk.simTrkIdx[0];  // get first matched track
                        for(const int& simHitIdx : sim_tracks[simTrkIdx].simHitIdx){ // loop over SimHits from SimTrack
                            if(simHitIdx == rec_hit1.simHitIdx[0]) // take first matched simhit (should be the closest one)
                                matched_hits.push_back({rec_hit1, sim_hits[rec_hit1.simHitIdx[0]]});
                        }
                    }
                }else{  // Looking at second hit in the pair
                    if(pixdet.contains(rec_hit2) && pixdet.contains((PixSubDet)rec_hit1.det, rec_hit2.lay+1+skip)){
                        // We have the RecHit we want to work with, now find a properly* matched SimHit
                        if(rec_hit2.simHitIdx.size() == 0) continue; // if rechit is matched to no simhits, give up.
                        if(trk.simTrkIdx.size() == 0) continue; // if track is matched to no simtracks, give up.
                        const int &simTrkIdx = trk.simTrkIdx[0];  // get first matched track
                        for(const int& simHitIdx : sim_tracks[simTrkIdx].simHitIdx){ // loop over SimHits from SimTrack
                            if(simHitIdx == rec_hit1.simHitIdx[0]) // take first matched simhit (should be the closest one)
                                matched_hits.push_back({rec_hit2, sim_hits[rec_hit2.simHitIdx[0]]});
                        }
                    }
                }
            }
            return matched_hits;
        })));

    auto BPIX_L1 = constant("BPIX_L1", PixDetElem(PixSubDet::Barrel, 1));
    auto BPIX_L2 = constant("BPIX_L2", PixDetElem(PixSubDet::Barrel, 2));
    auto FPIX_L1 = constant("FPIX_L1", PixDetElem(PixSubDet::EndCap, 1));
    auto FPIX_L2 = constant("FPIX_L2", PixDetElem(PixSubDet::EndCap, 2));

    auto skip_zero = constant("skip_zero", 0);
    auto skip_one  = constant("skip_one", 1);

    auto first = constant("first", true);
    auto second = constant("second", false);


    // First hit on 1st/2nd bpix layers and second hit in 2nd/3rd
    auto first_hits_in_B1_skip_zero = fv::tup_apply(find_matched_nth_hit_in_layer_with_skip,
            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, BPIX_L1, skip_zero, first), "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, BPIX_L2, skip_zero, first), "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, FPIX_L1, skip_zero, first), "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, FPIX_L2, skip_zero, first), "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, BPIX_L1, skip_one, first), "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, BPIX_L2, skip_one, first), "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, FPIX_L1, skip_one, first), "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){

    auto find_matched_nth_hit_in_layer =
    func<vector<HitPair>(vector<Seed>,
                         vector<PixRecHit>,
                         vector<SimHit>,
                         vector<Track>,
                         vector<SimTrack>,
                         PixDetElem, int)>("find_matched_nth_hit_in_layer",
        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 PixDetElem&& pixdet,
                 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
                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(pixdet.contains(rec_hit)){
                    // 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.
                    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]]});
                    }
                }
            }
            return matched_hits;
        })));

    auto BPIX_L1 = constant("BPIX_L1", PixDetElem(PixSubDet::Barrel, 1));
    auto BPIX_L2 = constant("BPIX_L2", PixDetElem(PixSubDet::Barrel, 2));
    auto BPIX_L3 = constant("BPIX_L3", PixDetElem(PixSubDet::Barrel, 3));
    auto BPIX_L4 = constant("BPIX_L4", PixDetElem(PixSubDet::Barrel, 4));
    auto FPIX_L1 = constant("FPIX_L1", PixDetElem(PixSubDet::EndCap, 1));
    auto FPIX_L2 = constant("FPIX_L2", PixDetElem(PixSubDet::EndCap, 2));
    auto FPIX_L3 = constant("FPIX_L3", PixDetElem(PixSubDet::EndCap, 3));

    auto first_hit  = constant("1st", 0);
    auto second_hit = constant("2nd", 1);

    // First hits on inner two bpix layers
    auto first_hits_in_B1 = fv::tup_apply(find_matched_nth_hit_in_layer,
            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, BPIX_L1, first_hit), "first_hits_in_B1");
    auto first_hits_in_B2 = fv::tup_apply(find_matched_nth_hit_in_layer,
            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, BPIX_L2, first_hit), "first_hits_in_B2");

    // Second hits on outer three bpix layers
    auto second_hits_in_B2 = fv::tup_apply(find_matched_nth_hit_in_layer,
            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, BPIX_L2, second_hit), "second_hits_in_B2");
    auto second_hits_in_B3 = fv::tup_apply(find_matched_nth_hit_in_layer,
            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, BPIX_L3, second_hit), "second_hits_in_B3");
    auto second_hits_in_B4 = fv::tup_apply(find_matched_nth_hit_in_layer,
            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, BPIX_L4, second_hit), "second_hits_in_B4");

    auto first_hits_in_F1 = fv::tup_apply(find_matched_nth_hit_in_layer,
            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, FPIX_L1, first_hit), "first_hits_in_F1");
    auto first_hits_in_F2 = fv::tup_apply(find_matched_nth_hit_in_layer,
            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, FPIX_L2, first_hit), "first_hits_in_F2");

    auto second_hits_in_F2 = fv::tup_apply(find_matched_nth_hit_in_layer,
            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, FPIX_L2, second_hit), "second_hits_in_F2");
    auto second_hits_in_F3 = fv::tup_apply(find_matched_nth_hit_in_layer,
            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, FPIX_L3, 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 =
    func<vector<HitPair>(vector<HitPair>, bool)>("select_even_odd_ladder_blade_hit_pairs",
        FUNC(([](const vector<HitPair>& hit_pairs,
                 const bool &&odd){
            vector<HitPair> even_pairs;
            for(const HitPair &hit_pair : hit_pairs){ // loop over all seeds
                if(std::get<PixRecHit>(hit_pair).ladder_blade % 2 == odd){
                    even_pairs.push_back(hit_pair);
                }
            }
            return even_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::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::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::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 = lookup_function<pair_vec(vector<HitPair>)>("calc_dphi_v_eta");
    TH2Params params_dphi = {"$\\eta$",       100, -4,   4,
                             "$\\Delta \\phi$(rad)", 150,  -.0015, .0015};
    tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, second_hits_in_F3),
                                                    "Second Hit in FPIX-L3", params_dphi);

    //Plots for dr of collections defined above
    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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, second_hits_in_F3),
                                                    "Second Hit in FPIX-L3", params_dr);


    //Plots for dz of collections defined above
    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, 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, 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, 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, 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, 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, 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, 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, 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, second_hits_in_B4),
                                                    "Second Hit in BPIX-L4", params_dz);

    //Plots for drho of collections defined above
    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, 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, 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, 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, second_hits_in_F3),
                                                    "Second Hit in FPIX-L3", params_drho);
}


void sc_value_hists_for_lb(TrackingDataSet& tds, std::function<bool(int)>& parity_check,
                           const std::string& label){

    auto relabel = [label](const std::string&& base){
        return base + "_" + label;
    };
    auto BPIX_L1 = constant("BPIX_L1", PixDetElem(PixSubDet::Barrel, 1));
    auto FPIX_L1 = constant("FPIX_L1", PixDetElem(PixSubDet::EndCap, 1));
    auto BPIX_L2 = constant("BPIX_L2", PixDetElem(PixSubDet::Barrel, 2));
    auto FPIX_L2 = constant("FPIX_L2", PixDetElem(PixSubDet::EndCap, 2));
    auto BPIX_L3 = constant("BPIX_L3", PixDetElem(PixSubDet::Barrel, 3));
    auto FPIX_L3 = constant("FPIX_L2", PixDetElem(PixSubDet::EndCap, 3));
    auto BPIX_L4 = constant("BPIX_L4", PixDetElem(PixSubDet::Barrel, 4));


    auto first_hit  = constant("1st", 0);
    auto second_hit = constant("2nd", 1);

    auto dRz = constant("dRz", true);
    auto dPhi = constant("dPhi", false);

    auto pcheck = constant(label, parity_check);
    auto sc_hit_residuals = lookup_function<
        pair_vec(vector<SuperCluster>, PixDetElem, int, bool, std::function<bool(int)>)>("sc_hit_residuals");

    // First hits on inner three bpix layers
    auto sc_first_hits_in_B1_dz = fv::tup_apply(sc_hit_residuals,
            fv::tuple(super_clusters, BPIX_L1, first_hit, dRz, pcheck), relabel("sc_first_hits_in_B1_dz"));
    auto sc_first_hits_in_B2_dz = fv::tup_apply(sc_hit_residuals,
            fv::tuple(super_clusters, BPIX_L2, first_hit, dRz, pcheck), relabel("sc_first_hits_in_B2_dz"));
    auto sc_first_hits_in_B3_dz = fv::tup_apply(sc_hit_residuals,
            fv::tuple(super_clusters, BPIX_L3, first_hit, dRz, pcheck), relabel("sc_first_hits_in_B3_dz"));

    auto sc_first_hits_in_B1_dphi = fv::tup_apply(sc_hit_residuals,
            fv::tuple(super_clusters, BPIX_L1, first_hit, dPhi, pcheck), relabel("sc_first_hits_in_B1_dphi"));
    auto sc_first_hits_in_B2_dphi = fv::tup_apply(sc_hit_residuals,
            fv::tuple(super_clusters, BPIX_L2, first_hit, dPhi, pcheck), relabel("sc_first_hits_in_B2_dphi"));
    auto sc_first_hits_in_B3_dphi = fv::tup_apply(sc_hit_residuals,
            fv::tuple(super_clusters, BPIX_L3, first_hit, dPhi, pcheck), relabel("sc_first_hits_in_B3_dphi"));

    // Second hits on outer three bpix layers
    auto sc_second_hits_in_B2_dz = fv::tup_apply(sc_hit_residuals,
            fv::tuple(super_clusters, BPIX_L2, second_hit, dRz, pcheck), relabel("sc_second_hits_in_B2_dz"));
    auto sc_second_hits_in_B3_dz = fv::tup_apply(sc_hit_residuals,
            fv::tuple(super_clusters, BPIX_L3, second_hit, dRz, pcheck), relabel("sc_second_hits_in_B3_dz"));
    auto sc_second_hits_in_B4_dz = fv::tup_apply(sc_hit_residuals,
            fv::tuple(super_clusters, BPIX_L4, second_hit, dRz, pcheck), relabel("sc_second_hits_in_B4_dz"));

    auto sc_second_hits_in_B2_dphi = fv::tup_apply(sc_hit_residuals,
            fv::tuple(super_clusters, BPIX_L2, second_hit, dPhi, pcheck), relabel("sc_second_hits_in_B2_dphi"));
    auto sc_second_hits_in_B3_dphi = fv::tup_apply(sc_hit_residuals,
            fv::tuple(super_clusters, BPIX_L3, second_hit, dPhi, pcheck), relabel("sc_second_hits_in_B3_dphi"));
    auto sc_second_hits_in_B4_dphi = fv::tup_apply(sc_hit_residuals,
            fv::tuple(super_clusters, BPIX_L4, second_hit, dPhi, pcheck), relabel("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>>(relabel("sc_first_hits_in_B1_dz"),
                                                    sc_first_hits_in_B1_dz,
                                                    "First Hit in BPIX-L1", params_dz);
    tds.register_container<ContainerTH2Many<float>>(relabel("sc_first_hits_in_B2_dz"),
                                                    sc_first_hits_in_B2_dz,
                                                    "First Hit in BPIX-L2", params_dz);
    tds.register_container<ContainerTH2Many<float>>(relabel("sc_first_hits_in_B3_dz"),
                                                    sc_first_hits_in_B3_dz,
                                                    "First Hit in BPIX-L3", params_dz);

    tds.register_container<ContainerTH2Many<float>>(relabel("sc_second_hits_in_B2_dz"),
                                                    sc_second_hits_in_B2_dz,
                                                    "Second Hit in BPIX-L2", params_dz);
    tds.register_container<ContainerTH2Many<float>>(relabel("sc_second_hits_in_B3_dz"),
                                                    sc_second_hits_in_B3_dz,
                                                    "Second Hit in BPIX-L3", params_dz);
    tds.register_container<ContainerTH2Many<float>>(relabel("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)", 1500, -0.7,  0.7};
    tds.register_container<ContainerTH2Many<float>>(relabel("sc_first_hits_in_B1_dphi"),
                                                    sc_first_hits_in_B1_dphi,
                                                    "First Hit in BPIX-L1", params_dphi);
    tds.register_container<ContainerTH2Many<float>>(relabel("sc_first_hits_in_B2_dphi"),
                                                    sc_first_hits_in_B2_dphi,
                                                    "First Hit in BPIX-L2", params_dphi);
    tds.register_container<ContainerTH2Many<float>>(relabel("sc_first_hits_in_B3_dphi"),
                                                    sc_first_hits_in_B3_dphi,
                                                    "First Hit in BPIX-L3", params_dphi);

    tds.register_container<ContainerTH2Many<float>>(relabel("sc_second_hits_in_B2_dphi"),
                                                    sc_second_hits_in_B2_dphi,
                                                    "Second Hit in BPIX-L2", params_dphi);
    tds.register_container<ContainerTH2Many<float>>(relabel("sc_second_hits_in_B3_dphi"),
                                                    sc_second_hits_in_B3_dphi,
                                                    "Second Hit in BPIX-L3", params_dphi);
    tds.register_container<ContainerTH2Many<float>>(relabel("sc_second_hits_in_B4_dphi"),
                                                    sc_second_hits_in_B4_dphi,
                                                    "Second Hit in BPIX-L4", params_dphi);
}

void setup_sc_residuals(TrackingDataSet& tds){

    /* sc_hit_residuals looks through all the super clusters in the event and
     * examins the matched pixel hits. It calculates a residual based on the
     * extrapolated position of the SC track and the position of the matched
     * hit.
     *   det: 1 for barrel, 2 for endcap
     *   pix_layer: 1-4 for barrel, 1-3 for endcap
     *   hit_number: examin either the first or second matched hit
     *   sel_dRz: True: calculate dR, False: calculate dz
     *   pcheck: a boolean function that filters residuals based on
     *   ladder/blade. Normally used to select either even or odd ladders.
     */
    func<pair_vec(vector<SuperCluster>,
                  PixDetElem, int, bool, std::function<bool(int)>)>("sc_hit_residuals",
        FUNC(([](const vector<SuperCluster>& super_clusters,
                 const PixDetElem&& pixdet,
                 const int&& hit_number,
                 const bool&& sel_dRz,
                 const std::function<bool(int)>& pcheck){
            vector<float> residuals;
            vector<float> etas;
            for(const SuperCluster& super_cluster : super_clusters){ // loop over all super-clusters
                float dRz, dPhi;
                int ladder_blade;
                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(!pixdet.contains((PixSubDet)super_cluster.subDet1[seedIdx], super_cluster.lay1[seedIdx])) continue;
                        dRz = super_cluster.dRz1[seedIdx];
                        dPhi = super_cluster.dPhi1[seedIdx];
                        ladder_blade = super_cluster.ladder_blade1[seedIdx];
                    }else{
                        if(!pixdet.contains((PixSubDet)super_cluster.subDet2[seedIdx], super_cluster.lay2[seedIdx])) continue;
                        dRz = super_cluster.dRz2[seedIdx];
                        dPhi = super_cluster.dPhi2[seedIdx];
                        ladder_blade = super_cluster.ladder_blade2[seedIdx];
                    }
                    DEBUG("New Seed, Idx: " << seedIdx << endl
                          << "\tdRz: " << dRz << endl
                          << "\tdPhi: " << dPhi << endl);
                    if(!pcheck(ladder_blade)) continue;
                    if(sel_dRz)
                        residuals.push_back(dRz);
                    else
                        residuals.push_back(dPhi);
                    etas.push_back(pseudorapidityP(super_cluster));
                }
            }
            return std::make_pair(etas,residuals);
        })));

    std::function<bool(int)> even = [](const int& t){ return !(t%2); };
    std::function<bool(int)> odd = [](const int& t){ return (t%2); };

    sc_value_hists_for_lb(tds, even, "even");
    sc_value_hists_for_lb(tds, odd, "odd");

}


void setup_residual_cross_corrolations(TrackingDataSet& tds){
    /*
     * Interesting cross-corrolations
     * for hit1 in L1 or L2
     * dphi2 v dz1
     * dz2 v dz1
     * dphi2 v dphi1
     * dz2 v dphi1
     */
    enum Res{
        dz1,
        dz2,
        dPhi1,
        dPhi2
    };


    auto sc_cross_correlations = func<pair_vec(vector<SuperCluster>,
                  PixDetElem, PixDetElem, Res, Res)>("sc_cross_correlations",
        FUNC(([](const vector<SuperCluster>& super_clusters,
                 const PixDetElem&& subdet1,
                 const PixDetElem&& subdet2,
                 const Res&& res1,
                 const Res&& res2){
            vector<float> residuals1;
            vector<float> residuals2;
            for(const SuperCluster& super_cluster : super_clusters){ // loop over all supser-clusters
                unsigned int nSeeds = super_cluster.charge.size();
                for(unsigned int seedIdx=0; seedIdx<nSeeds; seedIdx++){ // loop over all seeds associated w/ sc
                    if (!subdet1.contains((PixSubDet)super_cluster.subDet1[seedIdx],super_cluster.lay1[seedIdx]) ||
                        !subdet2.contains((PixSubDet)super_cluster.subDet2[seedIdx],super_cluster.lay2[seedIdx]))
                        continue;
                    switch(res1){
                        case Res::dz1:
                            residuals1.push_back(super_cluster.dRz1[seedIdx]);
                            break;
                        case Res::dz2:
                            residuals1.push_back(super_cluster.dRz2[seedIdx]);
                            break;
                        case Res::dPhi1:
                            residuals1.push_back(super_cluster.dPhi1[seedIdx]);
                            break;
                        case Res::dPhi2:
                            residuals1.push_back(super_cluster.dPhi2[seedIdx]);
                            break;
                    }
                    switch(res2){
                        case Res::dz1:
                            residuals2.push_back(super_cluster.dRz1[seedIdx]);
                            break;
                        case Res::dz2:
                            residuals2.push_back(super_cluster.dRz2[seedIdx]);
                            break;
                        case Res::dPhi1:
                            residuals2.push_back(super_cluster.dPhi1[seedIdx]);
                            break;
                        case Res::dPhi2:
                            residuals2.push_back(super_cluster.dPhi2[seedIdx]);
                            break;
                    }
                }
            }
            return std::make_pair(residuals1,residuals2);
        })));

    auto BPIX_L1 = constant("BPIX_L1", PixDetElem(PixSubDet::Barrel, 1));
    auto BPIX_L2 = constant("BPIX_L2", PixDetElem(PixSubDet::Barrel, 2));
    auto BPIX_L3 = constant("BPIX_L3", PixDetElem(PixSubDet::Barrel, 3));
    auto BPIX_L4 = constant("BPIX_L4", PixDetElem(PixSubDet::Barrel, 4));

    auto v_dz1 = constant("dz1", Res::dz1);
    auto v_dz2 = constant("dz2", Res::dz2);
    auto v_dPhi1 = constant("dPhi1", Res::dPhi1);
    auto v_dPhi2 = constant("dPhi2", Res::dPhi2);

    auto sc_dphi2_v_dz1_L1_L2 = fv::tup_apply(sc_cross_correlations,
            fv::tuple(super_clusters, BPIX_L1, BPIX_L2, v_dPhi2, v_dz1), "sc_dphi1_v_dz1_L1_L2");

    auto sc_dz2_v_dz1_L1_L2 = fv::tup_apply(sc_cross_correlations,
            fv::tuple(super_clusters, BPIX_L1, BPIX_L2, v_dz2, v_dz1), "sc_dz2_v_dz1_L1_L2");

    auto sc_dphi2_v_dphi1_L1_L2 = fv::tup_apply(sc_cross_correlations,
            fv::tuple(super_clusters, BPIX_L1, BPIX_L2, v_dPhi2, v_dPhi1), "sc_dphi1_v_dphi1_L1_L2");

    auto sc_dz2_v_dphi1_L1_L2 = fv::tup_apply(sc_cross_correlations,
            fv::tuple(super_clusters, BPIX_L1, BPIX_L2, v_dz2, v_dPhi1), "sc_dz2_v_dphi1_L1_L2");

    TH2Params hist_params = TH2Params::lookup("sc_dphi2_v_dz1_L1_L2");
    tds.register_container<ContainerTH2Many<float>>("sc_dphi2_v_dz1_L1_L2",
                                                    sc_dphi2_v_dz1_L1_L2,
                                                    "sc_dphi2_v_dz1_L1_L2", hist_params);
    hist_params = TH2Params::lookup("sc_dz2_v_dz1_L1_L2");
    tds.register_container<ContainerTH2Many<float>>("sc_dz2_v_dz1_L1_L2",
                                                    sc_dz2_v_dz1_L1_L2,
                                                    "sc_dz2_v_dz1_L1_L2", hist_params);
    hist_params = TH2Params::lookup("sc_dphi2_v_dphi1_L1_L2");
    tds.register_container<ContainerTH2Many<float>>("sc_dphi2_v_dphi1_L1_L2",
                                                    sc_dphi2_v_dphi1_L1_L2,
                                                    "sc_dphi2_v_dphi1_L1_L2", hist_params);
    hist_params = TH2Params::lookup("sc_dz2_v_dphi1_L1_L2");
    tds.register_container<ContainerTH2Many<float>>("sc_dz2_v_dphi1_L1_L2",
                                                    sc_dz2_v_dphi1_L1_L2,
                                                    "sc_dz2_v_dphi1_L1_L2", hist_params);

    hist_params = TH2Params::lookup("sc_dz2_v_dz1_L1_L2_zoom");
    tds.register_container<ContainerTH2Many<float>>("sc_dz2_v_dz1_L1_L2_zoom",
                                                    sc_dz2_v_dz1_L1_L2,
                                                    "sc_dz2_v_dz1_L1_L2_zoom", hist_params);
    hist_params = TH2Params::lookup("sc_dz2_v_dphi1_L1_L2_zoom");
    tds.register_container<ContainerTH2Many<float>>("sc_dz2_v_dphi1_L1_L2_zoom",
                                                    sc_dz2_v_dphi1_L1_L2,
                                                    "sc_dz2_v_dphi1_L1_L2_zoom", hist_params);
    hist_params = TH2Params::lookup("sc_dphi2_v_dphi1_L1_L2_zoom");
    tds.register_container<ContainerTH2Many<float>>("sc_dphi2_v_dphi1_L1_L2_zoom",
                                                    sc_dphi2_v_dphi1_L1_L2,
                                                    "sc_dphi2_v_dphi1_L1_L2_zoom", hist_params);

    auto sc_dphi2_v_dz1_L1_L3 = fv::tup_apply(sc_cross_correlations,
            fv::tuple(super_clusters, BPIX_L1, BPIX_L3, v_dPhi2, v_dz1), "sc_dphi1_v_dz1_L1_L3");

    auto sc_dz2_v_dz1_L1_L3 = fv::tup_apply(sc_cross_correlations,
            fv::tuple(super_clusters, BPIX_L1, BPIX_L3, v_dz2, v_dz1), "sc_dz2_v_dz1_L1_L3");

    auto sc_dphi2_v_dphi1_L1_L3 = fv::tup_apply(sc_cross_correlations,
            fv::tuple(super_clusters, BPIX_L1, BPIX_L3, v_dPhi2, v_dPhi1), "sc_dphi1_v_dphi1_L1_L3");

    auto sc_dz2_v_dphi1_L1_L3 = fv::tup_apply(sc_cross_correlations,
            fv::tuple(super_clusters, BPIX_L1, BPIX_L3, v_dz2, v_dPhi1), "sc_dz2_v_dphi1_L1_L3");

    hist_params = TH2Params::lookup("sc_dphi2_v_dz1_L1_L2");
    tds.register_container<ContainerTH2Many<float>>("sc_dphi2_v_dz1_L1_L3",
                                                    sc_dphi2_v_dz1_L1_L3,
                                                    "sc_dphi2_v_dz1_L1_L3", hist_params);
    hist_params = TH2Params::lookup("sc_dz2_v_dz1_L1_L2");
    tds.register_container<ContainerTH2Many<float>>("sc_dz2_v_dz1_L1_L3",
                                                    sc_dz2_v_dz1_L1_L3,
                                                    "sc_dz2_v_dz1_L1_L3", hist_params);
    hist_params = TH2Params::lookup("sc_dphi2_v_dphi1_L1_L2");
    tds.register_container<ContainerTH2Many<float>>("sc_dphi2_v_dphi1_L1_L3",
                                                    sc_dphi2_v_dphi1_L1_L3,
                                                    "sc_dphi2_v_dphi1_L1_L3", hist_params);
    hist_params = TH2Params::lookup("sc_dz2_v_dphi1_L1_L2");
    tds.register_container<ContainerTH2Many<float>>("sc_dz2_v_dphi1_L1_L3",
                                                    sc_dz2_v_dphi1_L1_L3,
                                                    "sc_dz2_v_dphi1_L1_L3", hist_params);

    hist_params = TH2Params::lookup("sc_dz2_v_dz1_L1_L2_zoom");
    tds.register_container<ContainerTH2Many<float>>("sc_dz2_v_dz1_L1_L3_zoom",
                                                    sc_dz2_v_dz1_L1_L3,
                                                    "sc_dz2_v_dz1_L1_L3_zoom", hist_params);
    hist_params = TH2Params::lookup("sc_dz2_v_dphi1_L1_L2_zoom");
    tds.register_container<ContainerTH2Many<float>>("sc_dz2_v_dphi1_L1_L3_zoom",
                                                    sc_dz2_v_dphi1_L1_L3,
                                                    "sc_dz2_v_dphi1_L1_L3_zoom", hist_params);
    hist_params = TH2Params::lookup("sc_dphi2_v_dphi1_L1_L2_zoom");
    tds.register_container<ContainerTH2Many<float>>("sc_dphi2_v_dphi1_L1_L3_zoom",
                                                    sc_dphi2_v_dphi1_L1_L3,
                                                    "sc_dphi2_v_dphi1_L1_L3_zoom", hist_params);


}

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){
        return input.substr(0, input.find_last_of(".")) + new_suffix;
    };

    string log_filename = replace_suffix(output_filename, ".log");
    fv::util::Log::init_logger(log_filename, fv::util::LogPriority::kLogInfo);

    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);
    setup_residual_cross_corrolations(tds);


    tds.process(silent);
    tds.save_all();
}

int main(int argc, char * argv[]){
    fv::util::ArgParser args(argc, argv);
    bool silent = args.cmdOptionExists("-s");
    string output_filename = args.cmdOptionExists("-o") ? args.getCmdOption("-o") : "output.root";
    if(args.cmdOptionExists("-c")){
        fv::util::init_config(args.getCmdOption("-c"));
    }
    if(args.cmdOptionExists("-F")){
        auto file_list = fv::util::read_input_list(args.getCmdOption("-F"));
        run_analysis(file_list, output_filename, silent);
    }else if(args.cmdOptionExists("-f")){
        string input_filename = args.getCmdOption("-f");
        string data_label = args.cmdOptionExists("-l") ? args.getCmdOption("-l") : "";
        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} {-c DATA_CATEGORY} {-o output_filename} -f treefile.root" << endl;
    }
    return 0;
}

/* int main(int argc, char * argv[]){ */
/*     using namespace fv::util; */
/*     init_config(argv[1]); */
/*     int max_events = the_config.get("max-events").as<int>(); */
/*     cout << max_events << endl; */
/* } */