EGamma_ElectronTrackingValidation/looper/analysis/tracking_eff.cpp

#include <iostream>
#include <vector>
#include <cmath>
#include <numeric>
#include <utility>
#include <TSystem.h>
#include <Math/VectorUtil.h>

#include "filval.hpp"
#include "root_filval.hpp"
#include "common.hpp"

#include "analysis/TrackingNtupleObjs.hpp"

const double pi = 3.1415926535897932384626;
const double pi2 = 2*pi;

typedef ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<float>> Vec4;

SeedCollection seeds;
GenCollection gens;
SimTrackCollection sim_tracks;
SimVertexCollection sim_vertices;
TrackCollection gsf_tracks;
SuperClusterCollection scls;

std::vector<SimTrack> sim_els;

Value<float>* truePileup;
Value<int>* inTimePileup;
/* Value<vector<int>>* outOfTimePileup; */

/* Value<vector<unsigned int>>* totalSeedsAlgo; */
Value<unsigned int>* totalSeeds;

void register_objects(TrackingDataSet& tds){
    seeds.init(&tds);
    gens.init(&tds);
    sim_tracks.init(&tds);
    sim_vertices.init(&tds);
    gsf_tracks.init(&tds);
    scls.init(&tds);
}

struct {
    Value<float>* x;
    Value<float>* y;
    Value<float>* z;
    Value<float>* sigma_x;
    Value<float>* sigma_y;
    Value<float>* sigma_z;
} bs;

struct EffColl {
    EfficiencyContainer<float>* pt;
    EfficiencyContainer<float>* eta;
    EfficiencyContainer<float>* phi;
    EffColl(TrackingDataSet& tds, const std::string& pfx,
            THParams pt_params, THParams eta_params, THParams phi_params,
            const std::string& sfx="") {
        std::string _sfx = "";
        if (sfx != "") _sfx = "_"+sfx;

        pt = tds.register_container<EfficiencyContainer<float>>(pfx+"_v_pt"+_sfx, pt_params);
        eta = tds.register_container<EfficiencyContainer<float>>(pfx+"_v_eta"+_sfx, eta_params);
        phi = tds.register_container<EfficiencyContainer<float>>(pfx+"_v_phi"+_sfx, phi_params);
    }
};

struct KinColl {
    ContainerTH1<float>* pt;
    ContainerTH1<float>* eta;
    ContainerTH1<float>* phi;
    ContainerTH2<float>* eta_phi;

    KinColl(TrackingDataSet& tds, const std::string& pfx,
            THParams pt_params, THParams eta_params, THParams phi_params,
            const std::string& sfx="") {
        std::string _sfx = "";
        if (sfx != "") _sfx = "_"+sfx;

        pt = tds.register_container<ContainerTH1<float>>(pfx+"_v_pt"+_sfx, pt_params);
        eta = tds.register_container<ContainerTH1<float>>(pfx+"_v_eta"+_sfx, eta_params);
        phi = tds.register_container<ContainerTH1<float>>(pfx+"_v_phi"+_sfx, phi_params);

        THParams eta_phi_params = eta_params;
        eta_phi_params.low_y = phi_params.low_x;
        eta_phi_params.high_y = phi_params.high_x;
        eta_phi_params.label_y = phi_params.label_x;

        eta_phi = tds.register_container<ContainerTH2<float>>(pfx+"_v_eta_phi"+_sfx, eta_phi_params);
    }
};

template <typename T1, typename T2>
struct AugKinColl {
    T1* pt;
    T1* eta;
    T1* phi;
    T2* eta_phi;
};

int get_n_pv() {
    /* std::cout << "------------------------------------------------" << std::endl; */
    int n_pv = 0;
    for (const auto& sim_vertex : sim_vertices) {
        if (sim_vertex.sourceSimIdx().size() == 0) {
            /* std::cout << "(" */
            /*     << sim_vertex.x() << ", " */
            /*     << sim_vertex.y() << ", " */
            /*     << sim_vertex.z() << ") " */
            /*     << sim_vertex.processType() << std::endl; */
            n_pv++;
        }
    }
    /* std::cout << "::::: " << pileup << std::endl; */
    return n_pv / 2; // for some reason, all vertices are repeated twice...
}

bool in_lum_region(const SimVertex& vertex) {
    /* Here approximate the luminous region as a cylinder with length 5*\sigma_z and radius
     * sqrt(\simga_x^2 + \sigma_y^2) centered as beamspot position.
     */
    float dx = vertex.x() - bs.x->get();
    float dy = vertex.y() - bs.y->get();
    float dz = vertex.z() - bs.z->get();
    auto radius = static_cast<float>(5 * sqrt(pow(bs.sigma_x->get(), 2) + pow(bs.sigma_y->get(), 2)));
    float half_len = 5*bs.sigma_z->get();
    return sqrt(dx*dx + dy*dy) < radius and abs(dz) < half_len;
};

bool is_good_sim(const SimTrack& sim_track) {
    const auto& vertex = sim_vertices[sim_track.parentVtxIdx()];
    return abs(sim_track.pdgId()) == 11 and  // electrons
           in_lum_region(vertex) and         // from luminous region
           abs(sim_track.eta())<3.0;         // inside ECAL acceptance
};

bool is_good_sim(const std::vector<SimTrack> prompt_sims, const SimTrack& sim_track) {
    return find(prompt_sims.begin(), prompt_sims.end(), sim_track) != prompt_sims.end();
};

bool is_good_seed(const Seed& seed, float hoe_cut) {
    return seed.isECALDriven() and seed.hoe() < hoe_cut;
}

float reco_energy_rel_err(const SimTrack& sim_track, const Seed& seed) {
    return (sim_track.pt() - seed.Et()) / sim_track.pt() ;
}

float reco_energy_rel_err(const SimTrack& sim_track, const Track& track) {
    return (sim_track.pt() - track.pt()) / sim_track.pt() ;
}

template<typename TrackOrSeed>
bool reco_energy_consistent(const SimTrack& sim_track, const TrackOrSeed& track, float consistency_cut=0.1) {
    return fabs(reco_energy_rel_err(sim_track, track)) < consistency_cut;
}

Vec4 scl_p4(const SuperCluster& scl) {
    return Vec4(scl.px(), scl.py(), scl.pz(), scl.e());
}

double dphi(double phi1, double phi2) {
    auto spec_mod = [](double a) {
        return a - floor(a/pi2)*pi2;
    };
    return spec_mod(phi2 - phi1 + pi) - pi;
}


double deltaR(double eta1, double phi1, double eta2, double phi2) {
    return sqrt(pow(eta1 - eta2, 2.0) +
                pow(dphi(phi1, phi2), 2.0));
}

void update_sim_els() {
    sim_els.clear();
    for (const SimTrack sim_track: sim_tracks) {
        if (is_good_sim(sim_track)) {
            sim_els.push_back(sim_track);
        }
    }
}

std::tuple<std::vector<SimTrack>, std::vector<SimTrack>> get_prompt_sims() {
    /*
     * Returns sim tracks that pass is_good_sim as well as match well with a
     * prompt gen electron
     */
    std::vector<Gen> prompt_gen;
    for (const auto& gen : gens) {
        if (abs(gen.pdgId()) == 11 and gen.isPrompt()) {
            prompt_gen.push_back(gen);
        }
    }

    std::vector<SimTrack> prompt_sim;
    for (const auto& gen : prompt_gen) {
        float gen_phi = atan2(gen.py(), gen.px());
        float gen_eta = pseudorapidityP(gen.px(), gen.py(), gen.pz());
        int gen_id = gen.pdgId();
        for (const auto& sim_el : sim_els) {
            float sim_phi = sim_el.phi();
            float sim_eta = sim_el.eta();
            int sim_id = sim_el.pdgId();
            if (gen_id == sim_id and deltaR(gen_eta, gen_phi, sim_eta, sim_phi) < 0.05) {
                prompt_sim.push_back(sim_el);
                break;
            }
        }
    }
    std::vector<SimTrack> nonprompt_sim;
    for (const auto& sim_el : sim_els) {
        if (find(prompt_sim.begin(), prompt_sim.end(), sim_el) == prompt_sim.end()) {
            nonprompt_sim.push_back(sim_el);
        }
    }
    return make_tuple(prompt_sim, nonprompt_sim);
}


auto get_match_stats(bool dRMatch) {
    int nSimTrack = 0;
    int nGSFTrack = 0;
    int nMatched = 0;  // 1-to-1
    int nMerged = 0; // n-to-1
    int nLost = 0;   // 1-to-0
    int nSplit = 0;  // 1-to-n
    int nFaked = 0;   // 0-to-1
    int nFlipped = 0;
    vector<float> matched_dR;
    vector<float> matched_dpT;
    std::map<int, std::vector<int>> sim_matches;
    std::map<int, std::vector<int>> gsf_matches;
    for(const auto& sim_track : sim_els)
        sim_matches[sim_track.idx] = {};
    for(const auto& gsf_track : gsf_tracks)
        gsf_matches[gsf_track.idx] = {};

    nSimTrack = static_cast<int>(sim_matches.size());
    nGSFTrack = static_cast<int>(gsf_matches.size());

    for (const auto& sim_track : sim_els) {
        if (dRMatch) {
            for (const auto& gsf_track : gsf_tracks) {
                double dR = deltaR(sim_track.eta(), sim_track.phi(), gsf_track.eta(), gsf_track.phi());
                if (dR < 0.2) {
                    sim_matches[sim_track.idx].push_back(gsf_track.idx);
                    gsf_matches[gsf_track.idx].push_back(sim_track.idx);
                    matched_dR.push_back(dR);
                    matched_dpT.push_back((sim_track.pt() - gsf_track.pt())/sim_track.pt());
                    if (gsf_track.q() != sim_track.q()) {
                        nFlipped++;
                    }
                }
            }
        } else {
            for (const auto& gsfIdx : sim_track.trkIdx()) {
                const Track& gsf_track = gsf_tracks[gsfIdx];
                double dR = deltaR(sim_track.eta(), sim_track.phi(), gsf_track.eta(), gsf_track.phi());
                sim_matches[sim_track.idx].push_back(gsf_track.idx);
                gsf_matches[gsf_track.idx].push_back(sim_track.idx);
                matched_dR.push_back(static_cast<float &&>(dR));
                matched_dpT.push_back((sim_track.pt() - gsf_track.pt())/sim_track.pt());
                if (gsf_track.q() != sim_track.q()) {
                    nFlipped++;
                }
            }
        }
    }

    for (const auto& p : sim_matches) {
        const auto& matches = p.second;
        if (matches.size() == 0) nLost++;
        if (matches.size() == 1 and gsf_matches[matches[0]].size() == 1) nMatched++;
        if (matches.size() > 1) nSplit++;
    }
    for (const auto& p : gsf_matches) {
        const auto& matches = p.second;
        if (matches.size() == 0) nFaked++;
        if (matches.size() > 1) nMerged++;
    }
    return std::make_tuple(nSimTrack, nGSFTrack, nMatched, nMerged, nLost, nSplit, nFaked, nFlipped, matched_dR, matched_dpT);
}

void run(){
    using namespace std;
    using namespace fv;
    using namespace fv_root;
    using namespace fv_util;

    auto lookup = [](const std::string&& key) {
        return THParams::lookup(std::forward<const std::string>(key));
    };

    auto file_list = the_config->get_source_files();
    string output_filename = the_config->get_output_filename();
    TrackingDataSet tds(output_filename, file_list, "trackingNtuple/tree");
    register_objects(tds);

    auto hoe_cut = the_config->get("hoe_cut").as<float>(999);

    bs = {
        tds.track_branch<float>("bsp_x"),
        tds.track_branch<float>("bsp_y"),
        tds.track_branch<float>("bsp_z"),
        tds.track_branch<float>("bsp_sigmax"),
        tds.track_branch<float>("bsp_sigmay"),
        tds.track_branch<float>("bsp_sigmaz")
    };

    inTimePileup = tds.track_branch<int>("inTimePileup");
    truePileup = tds.track_branch<float>("truePileup");


    totalSeeds = tds.track_branch<unsigned int>("initsee_tot");

    enum TMType {
        NoMatch = 0,
        SeedMatched = 1,
        TrackMatched = 2,
        SeedAndTrackMatched = 3,
    };

    std::map<std::tuple<int, int>, ContainerTH2<float>*> residuals_dRz;
    std::map<std::tuple<int, int>, ContainerTH2<float>*> residuals_dPhi;

    std::map<std::tuple<int, int, int>, ContainerTH2<float>*> BPIX_residuals_dRz;
    std::map<std::tuple<int, int, int>, ContainerTH2<float>*> BPIX_residuals_dPhi;
    std::map<std::tuple<int, int, int>, ContainerTH2<float>*> FPIX_residuals_dRz;
    std::map<std::tuple<int, int, int>, ContainerTH2<float>*> FPIX_residuals_dPhi;

    stringstream name;
    auto set_name = [&name](const std::string& var, const std::string& region,
                            const int& layer, const int& hit, const int& tm_type) {
        name.str("");
        name << var << "_" << region;
        if (layer>0)
            name << "_L" << layer;
        name << "_H" << hit << "_v_Et";
        switch (tm_type) {
            case NoMatch:
                name << "_NoMatch";
                break;
            case SeedMatched:
                name << "_SeedMatched";
                break;
            case TrackMatched:
                name << "_TrackMatched";
                break;
            case SeedAndTrackMatched:
                name << "_SeedAndTrackMatched";
                break;
            default:
                break;
        }
    };

    THParams hist_params;
    for (int hit=1; hit<=3; hit++) {
        for (int tm_type = NoMatch; tm_type <= SeedAndTrackMatched; tm_type++) {
            for (int layer=1; layer<=4; layer++) {
                hist_params = (hit == 1) ? lookup("dRz_v_Et") : lookup("dRz_v_Et_outer_hits");
                set_name("dRz", "BPIX", layer, hit, tm_type);
                BPIX_residuals_dRz[make_tuple(layer, hit, tm_type)] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
                set_name("dRz", "FPIX", layer, hit, tm_type);
                FPIX_residuals_dRz[make_tuple(layer, hit, tm_type)] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);

                hist_params = (hit == 1) ? lookup("dPhi_v_Et") : lookup("dPhi_v_Et_outer_hits");
                set_name("dPhi", "BPIX", layer, hit, tm_type);
                BPIX_residuals_dPhi[make_tuple(layer, hit, tm_type)] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
                set_name("dPhi", "FPIX", layer, hit, tm_type);
                FPIX_residuals_dPhi[make_tuple(layer, hit, tm_type)] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
            }
            hist_params = (hit == 1) ? lookup("dRz_v_Et") : lookup("dRz_v_Et_outer_hits");
            set_name("dRz", "ALL", 0, hit, tm_type);
            residuals_dRz[std::make_tuple(hit, tm_type)] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);

            hist_params = (hit == 1) ? lookup("dPhi_v_Et") : lookup("dPhi_v_Et_outer_hits");
            set_name("dPhi", "ALL", 0, hit, tm_type);
            residuals_dPhi[make_tuple(hit, tm_type)] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
        }
    }

    std::map<int,vector<float>> eta_regions = {{1, {0.0, 1.1, 1.8, 3.0}},
                                               {2, {0.0, 1.4, 2.3, 3.0}},
                                               {3, {0.0, 1.0, 2.0, 3.0}}};

    auto get_region = [&eta_regions](int hit, float eta) {
        auto regions = eta_regions[hit];
        if (regions.size() <= 1) return 1;
        float abseta = abs(eta);
        for (int r_idx=1; r_idx < regions.size(); r_idx++) {
            if (abseta > regions[r_idx-1] and abseta < regions[r_idx]) {
                return r_idx;
            }
        }
        return static_cast<int>(regions.size());
    };
    std::map<int, ContainerTH2<float>*> dPhi_residuals_v_eta;
    std::map<int, ContainerTH2<float>*> dRz_residuals_v_eta;

    std::map<std::pair<int, int>, ContainerTH2<float>*> dPhi_residuals_v_region;
    std::map<std::pair<int, int>, ContainerTH2<float>*> dRz_residuals_v_region;
    for (int hit=1; hit<=3; hit++) {
        name.str("");
        name << "dPhi_residuals_v_eta_H" << hit;
        hist_params = (hit==1) ? lookup("dPhi_v_eta") : lookup("dPhi_v_eta_outer_hits");
        dPhi_residuals_v_eta[hit] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);

        name.str("");
        name << "dRz_residuals_v_eta_H" << hit;
        hist_params = (hit==1) ? lookup("dRz_v_eta") : lookup("dRz_v_eta_outer_hits");
        dRz_residuals_v_eta[hit] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
        for (int region=1; region<=eta_regions[hit].size()-1; region++){
            hist_params = (hit==1) ? lookup("dRz_v_Et") : lookup("dRz_v_Et_outer_hits");
            name.str("");
            name << "dRz_residuals_H"  << hit << "_R" << region;
            dRz_residuals_v_region[{hit, region}] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);

            hist_params = (hit==1) ? lookup("dPhi_v_Et") : lookup("dPhi_v_Et_outer_hits");
            name.str("");
            name << "dPhi_residuals_H"  << hit << "_R" << region;
            dPhi_residuals_v_region[{hit, region}] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
        }
    }


    KinColl good_sim_kinems    (tds, "good_sim",     lookup("pt_fine"), lookup("eta_fine"), lookup("phi_fine"));
    KinColl gsf_kinems         (tds, "gsf",          lookup("pt_fine"), lookup("eta_fine"), lookup("phi_fine"));
    KinColl gsf_high_pt1_kinems(tds, "gsf_high_pt1", lookup("pt_fine"), lookup("eta_fine"), lookup("phi_fine"));
    KinColl gsf_high_pt2_kinems(tds, "gsf_high_pt2", lookup("pt_fine"), lookup("eta_fine"), lookup("phi_fine"));
    KinColl seed_kinems        (tds, "seed",         lookup("pt_fine"), lookup("eta_fine"), lookup("phi_fine"));
    KinColl scl_kinems         (tds, "scl",          lookup("pt_fine"), lookup("eta_fine"), lookup("phi_fine"));
    KinColl prompt_kinems      (tds, "prompt",       lookup("pt_fine"), lookup("eta_fine"), lookup("phi_fine"));
    KinColl nonprompt_kinems   (tds, "nonprompt",    lookup("pt_fine"), lookup("eta_fine"), lookup("phi_fine"));


    EffColl seed_eff(tds, "seed_eff", lookup("pt"), lookup("eta"), lookup("phi"));
    EffColl seed_pur(tds, "seed_pur", lookup("pt"), lookup("eta"), lookup("phi"));

    EffColl tracking_eff (tds, "tracking_eff",  lookup("pt"), lookup("eta"), lookup("phi"));
    EffColl tracking_pur (tds, "tracking_pur",  lookup("pt"), lookup("eta"), lookup("phi"));
    EffColl prompt_eff   (tds, "prompt_eff",    lookup("pt"), lookup("eta"), lookup("phi"));
    EffColl prompt_pur   (tds, "prompt_pur",    lookup("pt"), lookup("eta"), lookup("phi"));
    EffColl nonprompt_eff(tds, "nonprompt_eff", lookup("pt"), lookup("eta"), lookup("phi"));
    EffColl nonprompt_pur(tds, "nonprompt_pur", lookup("pt"), lookup("eta"), lookup("phi"));

    EffColl tracking_eff_dR (tds, "tracking_eff",  lookup("pt"), lookup("eta"), lookup("phi"), "dR");
    EffColl tracking_pur_dR (tds, "tracking_pur",  lookup("pt"), lookup("eta"), lookup("phi"), "dR");
    EffColl prompt_eff_dR   (tds, "prompt_eff",    lookup("pt"), lookup("eta"), lookup("phi"), "dR");
    EffColl prompt_pur_dR   (tds, "prompt_pur",    lookup("pt"), lookup("eta"), lookup("phi"), "dR");
    EffColl nonprompt_eff_dR(tds, "nonprompt_eff", lookup("pt"), lookup("eta"), lookup("phi"), "dR");
    EffColl nonprompt_pur_dR(tds, "nonprompt_pur", lookup("pt"), lookup("eta"), lookup("phi"), "dR");

    auto& tracking_eff_v_PU_dR     = *tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_PU_dR",     lookup("PU"));
    auto& tracking_eff_v_PU_dR_inc = *tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_PU_dR_inc", lookup("PU"));
    auto& tracking_pur_v_PU_dR     = *tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_PU_dR",     lookup("PU"));
    auto& tracking_pur_v_PU_dR_inc = *tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_PU_dR_inc", lookup("PU"));

    // Full Efficiency/Purity includes tracker-driven seeds
    EffColl tracking_eff_full   (tds, "tracking_eff_full", lookup("pt"), lookup("eta"), lookup("phi"));
    EffColl tracking_eff_full_dR(tds, "tracking_eff_full", lookup("pt"), lookup("eta"), lookup("phi"), "dR");
    EffColl tracking_pur_full   (tds, "tracking_pur_full", lookup("pt"), lookup("eta"), lookup("phi"));
    EffColl tracking_pur_full_dR(tds, "tracking_pur_full", lookup("pt"), lookup("eta"), lookup("phi"), "dR");

    EffColl fake_rate           (tds, "fake_rate",            lookup("pt"), lookup("eta"), lookup("phi"));
    EffColl fake_rate_no_e_match(tds, "fake_rate_no_e_match", lookup("pt"), lookup("eta"), lookup("phi"));

    auto& fake_rate_no_e_match_v_PU     = *tds.register_container<EfficiencyContainer<float>>("fake_rate_no_e_match_v_PU",     lookup("PU"));
    auto& fake_rate_no_e_match_v_PU_inc = *tds.register_container<EfficiencyContainer<float>>("fake_rate_no_e_match_v_PU_inc", lookup("PU"));

    EffColl partial_fake_rate        (tds, "partial_fake_rate",         lookup("pt"), lookup("eta"), lookup("phi"));
    EffColl full_fake_rate           (tds, "full_fake_rate",            lookup("pt"), lookup("eta"), lookup("phi"));
    EffColl clean_fake_rate          (tds, "clean_fake_rate",           lookup("pt"), lookup("eta"), lookup("phi"));
    EffColl fake_rate_incl           (tds, "fake_rate_incl",            lookup("pt"), lookup("eta"), lookup("phi"));
    EffColl fake_rate_no_e_match_incl(tds, "fake_rate_no_e_match_incl", lookup("pt"), lookup("eta"), lookup("phi"));
    EffColl partial_fake_rate_incl   (tds, "partial_fake_rate_incl",    lookup("pt"), lookup("eta"), lookup("phi"));
    EffColl full_fake_rate_incl      (tds, "full_fake_rate_incl",       lookup("pt"), lookup("eta"), lookup("phi"));
    EffColl clean_fake_rate_incl     (tds, "clean_fake_rate_incl",      lookup("pt"), lookup("eta"), lookup("phi"));

    auto& hit_vs_layer_barrel  = *tds.register_container<ContainerTH2<int>>("hit_vs_layer_barrel", lookup("hit_vs_layer"));
    auto& hit_vs_layer_forward = *tds.register_container<ContainerTH2<int>>("hit_vs_layer_forward", lookup("hit_vs_layer"));

    auto& gsf_tracks_nmatch_sim_tracks = *tds.register_container<ContainerTH1<size_t>>("gsf_tracks_nmatch_sim_tracks ", lookup("gsf_tracks_nmatch_sim_tracks"));

    auto& ecal_energy_resolution = *tds.register_container<ContainerTH1<float>>("ecal_energy_resolution ", lookup("ecal_energy_resolution"));

    auto& n_seeds      = *tds.register_container<ContainerTH1<size_t>>("n_seeds",      lookup("n_seeds"));
    auto& n_good_seeds = *tds.register_container<ContainerTH1<size_t>>("n_good_seeds", lookup("n_seeds"));
    auto& n_good_sim   = *tds.register_container<ContainerTH1<size_t>>("n_good_sim",   lookup("n_seeds"));
    auto& n_prompt     = *tds.register_container<ContainerTH1<size_t>>("n_prompt",     lookup("n_seeds"));
    auto& n_nonprompt  = *tds.register_container<ContainerTH1<size_t>>("n_nonprompt",  lookup("n_seeds"));
    auto& n_gsf_tracks = *tds.register_container<ContainerTH1<size_t>>("n_gsf_track",  lookup("n_seeds"));
    auto& n_scl        = *tds.register_container<ContainerTH1<size_t>>("n_scl",        lookup("n_seeds"));
    auto& n_good_scl   = *tds.register_container<ContainerTH1<size_t>>("n_good_scl",   lookup("n_seeds"));

    auto& n_matched = *tds.register_container<ContainerTH1<int>>("n_matched", lookup("n_seeds"));
    auto& n_merged  = *tds.register_container<ContainerTH1<int>>("n_merged",  lookup("n_seeds"));
    auto& n_lost    = *tds.register_container<ContainerTH1<int>>("n_lost",    lookup("n_seeds"));
    auto& n_split   = *tds.register_container<ContainerTH1<int>>("n_split",   lookup("n_seeds"));
    auto& n_faked   = *tds.register_container<ContainerTH1<int>>("n_faked",   lookup("n_seeds"));
    auto& n_flipped = *tds.register_container<ContainerTH1<int>>("n_flipped", lookup("n_seeds"));

    auto& matched_dR  = *tds.register_container<ContainerTH1<float>>("matched_dR", lookup("matched_dR"));
    auto& matched_dpT = *tds.register_container<ContainerTH1<float>>("matched_dpT", lookup("matched_dpT"));

    auto& n_matched_dR = *tds.register_container<ContainerTH1<int>>("n_matched_dR", lookup("n_seeds"));
    auto& n_merged_dR  = *tds.register_container<ContainerTH1<int>>("n_merged_dR",  lookup("n_seeds"));
    auto& n_lost_dR    = *tds.register_container<ContainerTH1<int>>("n_lost_dR",    lookup("n_seeds"));
    auto& n_split_dR   = *tds.register_container<ContainerTH1<int>>("n_split_dR",   lookup("n_seeds"));
    auto& n_faked_dR   = *tds.register_container<ContainerTH1<int>>("n_faked_dR",   lookup("n_seeds"));
    auto& n_flipped_dR = *tds.register_container<ContainerTH1<int>>("n_flipped_dR", lookup("n_seeds"));

    auto& matched_dR_dR = *tds.register_container<ContainerTH1<float>>("matched_dR_dR", lookup("matched_dR"));
    auto& matched_dpT_dR = *tds.register_container<ContainerTH1<float>>("matched_dpT_dR", lookup("matched_dpT"));

    auto& tm_corr = *tds.register_container<ContainerTH2<int>>("tm_corr", THParams(2, -0.5, 1.5, 2, -0.5, 1.5));
    auto& n_pileup = *tds.register_container<ContainerTH1<int>>("n_pileup", lookup("n_seeds"));


    auto& n_seeds_v_PU = *tds.register_container<ContainerTH2<int>>("n_seeds_v_PU",           lookup("n_seeds_v_PU"));
    auto& n_good_seeds_v_PU = *tds.register_container<ContainerTH2<int>>("n_good_seeds_v_PU", lookup("n_seeds_v_PU"));
    auto& n_tracks_v_PU = *tds.register_container<ContainerTH2<int>>("n_tracks_v_PU",         lookup("n_seeds_v_PU"));
    auto& n_pv_v_PU = *tds.register_container<ContainerTH2<int>>("n_pv_v_PU",                 lookup("n_seeds_v_PU"));

    auto& n_seeds_v_tPU = *tds.register_container<ContainerTH2<float>>("n_seeds_v_tPU",           lookup("n_seeds_v_PU"));
    auto& n_good_seeds_v_tPU = *tds.register_container<ContainerTH2<float>>("n_good_seeds_v_tPU", lookup("n_seeds_v_PU"));
    auto& n_tracks_v_tPU = *tds.register_container<ContainerTH2<float>>("n_tracks_v_tPU",         lookup("n_seeds_v_PU"));
    auto& n_pv_v_tPU = *tds.register_container<ContainerTH2<float>>("n_pv_v_tPU",                 lookup("n_seeds_v_PU"));


    auto& n_seeds_v_n_scl = *tds.register_container<ContainerTH2<int>>("n_seeds_v_n_scl",                     lookup("n_seeds_v_PU"));
    auto& n_good_seeds_v_n_scl = *tds.register_container<ContainerTH2<int>>("n_good_seeds_v_n_scl",           lookup("n_seeds_v_PU"));
    auto& n_good_seeds_v_n_good_scl = *tds.register_container<ContainerTH2<int>>("n_good_seeds_v_n_good_scl", lookup("n_seeds_v_PU"));

    auto& n_initseeds_v_PU = *tds.register_container<ContainerTH2<int>>("n_initseeds_v_PU",     lookup("n_initseeds_v_PU"));
    auto& n_initseeds_v_tPU = *tds.register_container<ContainerTH2<float>>("n_initseeds_v_tPU", lookup("n_initseeds_v_PU"));

    while (tds.next()) {

        update_sim_els();
        vector<SimTrack> prompt_sims, nonprompt_sims;
        std::tie(prompt_sims, nonprompt_sims) = get_prompt_sims();

        int n_pv = get_n_pv();
        n_pileup.fill(n_pv);
        n_pv_v_PU.fill(n_pv, inTimePileup->get());
        n_pv_v_tPU.fill(n_pv, truePileup->get());

        size_t _n_good_seeds = 0;
        /* cout << "Event: " << tds.get_current_event() << endl; */
        int idx = 0;
        for (const auto& seed : seeds) {
            /* cout << idx << ") [" << seed.pt() << ", " << seed.eta() << ", " << seed.phi() << "] " << seed.sclIdx() << " "; */

            if (is_good_seed(seed, hoe_cut)) {
                _n_good_seeds++;
                seed_kinems.pt->fill(seed.pt());
                seed_kinems.eta->fill(seed.eta());
                seed_kinems.phi->fill(seed.phi());
                seed_kinems.eta_phi->fill(seed.eta(), seed.phi());
            }
            idx++;
        }
        size_t _n_good_scl = 0;
        for (const auto& scl : scls) {
            if (scl.hoe() < hoe_cut) {
                _n_good_scl++;
                scl_kinems.pt->fill(hypot(scl.px(), scl.py()));
                float eta = pseudorapidity(scl.x(), scl.y(), scl.z());
                float phi = atan2(scl.y(), scl.x());
                scl_kinems.eta->fill(eta);
                scl_kinems.phi->fill(phi);
                scl_kinems.eta_phi->fill(eta, phi);
            }
        }
        /* idx = 0; */
        /* for (const auto& track : gsf_tracks) { */
        /*     cout << idx << ") [" << track.pt() << ", " << track.eta() << ", " << track.phi() << "] " << track.seedIdx() << endl; */
        /*     idx++; */
        /* } */
        /* string tmp; */
        /* cin >> tmp; */
        /* cout << "total seeds: " << totalSeeds->get() << endl; */
        n_initseeds_v_PU.fill(totalSeeds->get(), inTimePileup->get());
        n_initseeds_v_tPU.fill(totalSeeds->get(), truePileup->get());

        n_seeds_v_n_scl.fill(seeds.size(), scls.size());
        n_good_seeds_v_n_scl.fill(_n_good_seeds, scls.size());
        n_good_seeds_v_n_good_scl.fill(_n_good_seeds, _n_good_scl);

        n_seeds.fill(seeds.size());
        n_good_seeds.fill(_n_good_seeds);

        n_seeds_v_PU.fill(seeds.size(), inTimePileup->get());
        n_good_seeds_v_PU.fill(_n_good_seeds, inTimePileup->get());

        n_seeds_v_tPU.fill(seeds.size(), truePileup->get());
        n_good_seeds_v_tPU.fill(_n_good_seeds, truePileup->get());

        for (const auto& sim_track : sim_els) {
            good_sim_kinems.pt->fill(sim_track.pt());
            good_sim_kinems.eta->fill(sim_track.eta());
            good_sim_kinems.phi->fill(sim_track.phi());
            good_sim_kinems.eta_phi->fill(sim_track.eta(), sim_track.phi());
        }

        n_prompt.fill(prompt_sims.size());
        for (const auto& prompt : prompt_sims) {
            prompt_kinems.pt->fill(prompt.pt());
            prompt_kinems.eta->fill(prompt.eta());
            prompt_kinems.phi->fill(prompt.phi());
            prompt_kinems.eta_phi->fill(prompt.eta(), prompt.phi());
        }

        n_nonprompt.fill(nonprompt_sims.size());
        for (const auto& nonprompt : nonprompt_sims) {
            nonprompt_kinems.pt->fill(nonprompt.pt());
            nonprompt_kinems.eta->fill(nonprompt.eta());
            nonprompt_kinems.phi->fill(nonprompt.phi());
            nonprompt_kinems.eta_phi->fill(nonprompt.eta(), nonprompt.phi());
        }

        for (const auto& gsf_track : gsf_tracks) {
            if (!is_good_seed(seeds[gsf_track.seedIdx()], hoe_cut)) continue;
            gsf_kinems.pt->fill(gsf_track.pt());
            gsf_kinems.eta->fill(gsf_track.eta());
            gsf_kinems.phi->fill(gsf_track.phi());
            gsf_kinems.eta_phi->fill(gsf_track.eta(), gsf_track.phi());
            if (gsf_track.pt() > 10) {
                gsf_high_pt1_kinems.pt->fill(gsf_track.pt());
                gsf_high_pt1_kinems.eta->fill(gsf_track.eta());
                gsf_high_pt1_kinems.phi->fill(gsf_track.phi());
                gsf_high_pt1_kinems.eta_phi->fill(gsf_track.eta(), gsf_track.phi());
            }
            if (gsf_track.pt() > 20) {
                gsf_high_pt2_kinems.pt->fill(gsf_track.pt());
                gsf_high_pt2_kinems.eta->fill(gsf_track.eta());
                gsf_high_pt2_kinems.phi->fill(gsf_track.phi());
                gsf_high_pt2_kinems.eta_phi->fill(gsf_track.eta(), gsf_track.phi());
            }

        }
        n_tracks_v_PU.fill(gsf_tracks.size(), inTimePileup->get());
        n_tracks_v_tPU.fill(gsf_tracks.size(), truePileup->get());

        n_scl.fill(scls.size());
        n_good_scl.fill(_n_good_scl);

        auto match_stats =  get_match_stats(false);
        n_good_sim.fill(std::get<0>(match_stats));
        n_gsf_tracks.fill(std::get<1>(match_stats));
        n_matched.fill(std::get<2>(match_stats));
        n_merged.fill(std::get<3>(match_stats));
        n_lost.fill(std::get<4>(match_stats));
        n_split.fill(std::get<5>(match_stats));
        n_faked.fill(std::get<6>(match_stats));
        n_flipped.fill(std::get<7>(match_stats));
        for(float dR : std::get<8>(match_stats)) matched_dR.fill(dR);
        for(float dpT : std::get<9>(match_stats)) matched_dpT.fill(dpT);

        match_stats =  get_match_stats(true);
        n_matched_dR.fill(std::get<2>(match_stats));
        n_merged_dR.fill(std::get<3>(match_stats));
        n_lost_dR.fill(std::get<4>(match_stats));
        n_split_dR.fill(std::get<5>(match_stats));
        n_faked_dR.fill(std::get<6>(match_stats));
        n_flipped_dR.fill(std::get<7>(match_stats));
        for(float dR : std::get<8>(match_stats)) matched_dR_dR.fill(dR);
        for(float dpT : std::get<9>(match_stats)) matched_dpT_dR.fill(dpT);

        for (const auto& sim_track : sim_els) {
            if (seeds.size() == 0) continue;
            for (const auto &seed_idx : sim_track.seedIdx()) {
                const Seed& seed = seeds[seed_idx];
                if (!is_good_seed(seed, hoe_cut)) continue;
                ecal_energy_resolution.fill(reco_energy_rel_err(sim_track, seed));
            }
        }

        // Seeding Efficiency
        for (const auto& sim_track : sim_els) {
            if (seeds.size() == 0) continue;

            bool matched = false;
            for (const auto& seed_idx : sim_track.seedIdx()) {
                const Seed& seed = seeds[seed_idx];
                if (is_good_seed(seed, hoe_cut)) {
                    matched=true;
                    break;
                }
            }
            if (abs(sim_track.eta()) < 2.5)
                seed_eff.pt->fill(sim_track.pt(), matched);
            if (sim_track.pt() > 20.0)
                seed_eff.eta->fill(sim_track.eta(), matched);
            if (abs(sim_track.eta()) < 2.5 and sim_track.pt() > 20)
                seed_eff.phi->fill(sim_track.phi(), matched);
        }

        // Tracking Efficiency
        for (const auto& sim_track : sim_els) {
            if (gsf_tracks.size() == 0) continue;
            bool matched = false;
            bool matched_any = false;
            for (const auto& track_idx : sim_track.trkIdx()) {
                const Seed& seed = seeds[gsf_tracks[track_idx].seedIdx()];
                if (is_good_seed(seed, hoe_cut)) {
                    matched = true;
                    matched_any = true;
                } else if (!seed.isECALDriven()) {
                    matched_any = true;
                }
            }
            if (abs(sim_track.eta()) < 2.5) {
                tracking_eff.pt->fill(sim_track.pt(), matched);
                tracking_eff_full.pt->fill(sim_track.pt(), matched_any);
            }
            if (sim_track.pt() > 20.0) {
                tracking_eff.eta->fill(sim_track.eta(), matched);
                tracking_eff_full.eta->fill(sim_track.eta(), matched_any);
            }
            if (abs(sim_track.eta()) < 2.5 and sim_track.pt() > 20) {
                tracking_eff.phi->fill(sim_track.phi(), matched);
                tracking_eff_full.phi->fill(sim_track.phi(), matched_any);
            }

            if (std::find(prompt_sims.begin(), prompt_sims.end(), sim_track) != prompt_sims.end()) {
                if (abs(sim_track.eta()) < 2.5)
                    prompt_eff.pt->fill(sim_track.pt(), matched);
                if (sim_track.pt() > 20.0)
                    prompt_eff.eta->fill(sim_track.eta(), matched);
                if (abs(sim_track.eta()) < 2.5 and sim_track.pt() > 20)
                    prompt_eff.phi->fill(sim_track.phi(), matched);
            }
            if (std::find(nonprompt_sims.begin(), nonprompt_sims.end(), sim_track) != nonprompt_sims.end()) {
                if (abs(sim_track.eta()) < 2.5)
                    nonprompt_eff.pt->fill(sim_track.pt(), matched);
                if (sim_track.pt() > 20.0)
                    nonprompt_eff.eta->fill(sim_track.eta(), matched);
                if (abs(sim_track.eta()) < 2.5 and sim_track.pt() > 20)
                    nonprompt_eff.phi->fill(sim_track.phi(), matched);
            }

            // dR-matching
            matched = false;
            matched_any = false;
            for (const auto& gsf_track : gsf_tracks) {
                const Seed& seed = seeds[gsf_track.seedIdx()];
                double dR = deltaR(sim_track.eta(), sim_track.phi(), gsf_track.eta(), gsf_track.phi());
                if (is_good_seed(seed, hoe_cut) and dR < 0.2) {
                    matched = true;
                    matched_any = true;
                } else if (!seed.isECALDriven()) {
                    matched_any = true;
                }
            }
            if (abs(sim_track.eta()) < 2.5) {
                tracking_eff_dR.pt->fill(sim_track.pt(), matched);
                tracking_eff_full_dR.pt->fill(sim_track.pt(), matched_any);
            }
            if (sim_track.pt() > 20.0) {
                tracking_eff_dR.eta->fill(sim_track.eta(), matched);
                tracking_eff_full_dR.eta->fill(sim_track.eta(), matched_any);
            }
            if (abs(sim_track.eta()) < 2.5 and sim_track.pt() > 20) {
                tracking_eff_dR.phi->fill(sim_track.phi(), matched);
                tracking_eff_full_dR.phi->fill(sim_track.phi(), matched_any);
                tracking_eff_v_PU_dR.fill(inTimePileup->get(), matched);
            }

            tracking_eff_v_PU_dR_inc.fill(inTimePileup->get(), matched);

            if (find(prompt_sims.begin(), prompt_sims.end(), sim_track) != prompt_sims.end()) {
                if (abs(sim_track.eta()) < 2.5)
                    prompt_eff_dR.pt->fill(sim_track.pt(), matched);
                if (sim_track.pt() > 20.0)
                    prompt_eff_dR.eta->fill(sim_track.eta(), matched);
                if (abs(sim_track.eta()) < 2.5 and sim_track.pt() > 20)
                    prompt_eff_dR.phi->fill(sim_track.phi(), matched);
            }
            if (find(nonprompt_sims.begin(), nonprompt_sims.end(), sim_track) != nonprompt_sims.end()) {
                if (abs(sim_track.eta()) < 2.5)
                    nonprompt_eff_dR.pt->fill(sim_track.pt(), matched);
                if (sim_track.pt() > 20.0)
                    nonprompt_eff_dR.eta->fill(sim_track.eta(), matched);
                if (abs(sim_track.eta()) < 2.5 and sim_track.pt() > 20)
                    nonprompt_eff_dR.phi->fill(sim_track.phi(), matched);
            }
        }

        // Seeding Purity
        for (const auto& seed : seeds) {
            if (!is_good_seed(seed, hoe_cut)) continue;
            bool match = false;
            for (const auto& sim_track_idx : seed.simTrkIdx()) {
                if (is_good_sim(sim_els, sim_tracks[sim_track_idx])) {
                    match = true;
                    break;
                }
            }
            if (abs(seed.eta()) < 2.5)
                seed_pur.pt->fill(seed.pt(), match);
            if (seed.pt() > 20)
                seed_pur.eta->fill(seed.eta(), match);
            if (abs(seed.eta()) < 2.5 and seed.pt() > 20)
                seed_pur.phi->fill(seed.phi(), match);
        }

        // Tracking Purity
        for (const auto& gsf_track : gsf_tracks) {
            gsf_tracks_nmatch_sim_tracks.fill(gsf_track.simTrkIdx().size());
            const Seed& seed = seeds[gsf_track.seedIdx()];
            /* if (!is_good_seed(seed, hoe_cut)) continue; */
            bool good_ecal_driven_seed = is_good_seed(seed, hoe_cut);
            bool good_seed = good_ecal_driven_seed or !seed.isECALDriven();
            bool matched = false;
            bool prompt_matched = false;
            bool nonprompt_matched = false;
            for (const auto& sim_track_idx : gsf_track.simTrkIdx()) {
                const auto& sim_track = sim_tracks[sim_track_idx];
                if (!is_good_sim(sim_els, sim_track)) continue;
                matched = true;
                if (find(prompt_sims.begin(), prompt_sims.end(), sim_track) != prompt_sims.end()) {
                    prompt_matched = true;
                }
                if (find(nonprompt_sims.begin(), nonprompt_sims.end(), sim_track) != nonprompt_sims.end()) {
                    nonprompt_matched = true;
                }
            }

            if (abs(gsf_track.eta()) < 2.5) {
                if (good_ecal_driven_seed) tracking_pur.pt->fill(gsf_track.pt(), matched);
                if (good_seed) tracking_pur_full.pt->fill(gsf_track.pt(), matched);
            }
            if (gsf_track.pt() > 20) {
                if (good_ecal_driven_seed) tracking_pur.eta->fill(gsf_track.eta(), matched);
                if (good_seed) tracking_pur_full.eta->fill(gsf_track.eta(), matched);
            }
            if (abs(gsf_track.eta()) < 2.5 and gsf_track.pt() > 20) {
                if (good_ecal_driven_seed) tracking_pur.phi->fill(gsf_track.phi(), matched);
                if (good_seed) tracking_pur_full.phi->fill(gsf_track.phi(), matched);
            }

            if (abs(gsf_track.eta()) < 2.5)
                if (good_ecal_driven_seed) prompt_pur.pt->fill(gsf_track.pt(), prompt_matched);
            if (gsf_track.pt() > 20)
                if (good_ecal_driven_seed) prompt_pur.eta->fill(gsf_track.eta(), prompt_matched);
            if (abs(gsf_track.eta()) < 2.5 and gsf_track.pt() > 20)
                if (good_ecal_driven_seed) prompt_pur.phi->fill(gsf_track.phi(), prompt_matched);

            if (abs(gsf_track.eta()) < 2.5)
                if (good_ecal_driven_seed) nonprompt_pur.pt->fill(gsf_track.pt(), nonprompt_matched);
            if (gsf_track.pt() > 20)
                if (good_ecal_driven_seed) nonprompt_pur.eta->fill(gsf_track.eta(), nonprompt_matched);
            if (abs(gsf_track.eta()) < 2.5 and gsf_track.pt() > 20)
                if (good_ecal_driven_seed) nonprompt_pur.phi->fill(gsf_track.phi(), nonprompt_matched);

            // dR-matching
            matched = false;
            prompt_matched = false;
            nonprompt_matched = false;
            for (const auto& sim_track : sim_els) {
                double dR = deltaR(sim_track.eta(), sim_track.phi(), gsf_track.eta(), gsf_track.phi());
                if (dR < 0.2) {
                    matched = true;
                    if (find(prompt_sims.begin(), prompt_sims.end(), sim_track) != prompt_sims.end()) {
                        prompt_matched = true;
                    }
                    if (find(nonprompt_sims.begin(), nonprompt_sims.end(), sim_track) != nonprompt_sims.end()) {
                        nonprompt_matched = true;
                    }
                }
            }

            if (abs(gsf_track.eta()) < 2.5) {
                if (good_ecal_driven_seed) tracking_pur_dR.pt->fill(gsf_track.pt(), matched);
                if (good_seed) tracking_pur_full_dR.pt->fill(gsf_track.pt(), matched);
            }
            if (gsf_track.pt() > 20.0) {
                if (good_ecal_driven_seed) tracking_pur_dR.eta->fill(gsf_track.eta(), matched);
                if (good_seed) tracking_pur_full_dR.eta->fill(gsf_track.eta(), matched);
            }
            if (abs(gsf_track.eta()) < 2.5 and gsf_track.pt() > 20) {
                if (good_ecal_driven_seed) tracking_pur_dR.phi->fill(gsf_track.phi(), matched);
                if (good_seed) tracking_pur_full_dR.phi->fill(gsf_track.phi(), matched);
                if (good_ecal_driven_seed) tracking_pur_v_PU_dR.fill(inTimePileup->get(), matched);
            }
            tracking_pur_v_PU_dR_inc.fill(inTimePileup->get(), matched);

            if (abs(gsf_track.eta()) < 2.5)
                if (good_ecal_driven_seed) prompt_pur_dR.pt->fill(gsf_track.pt(), prompt_matched);
            if (gsf_track.pt() > 20)
                if (good_ecal_driven_seed) prompt_pur_dR.eta->fill(gsf_track.eta(), prompt_matched);
            if (abs(gsf_track.eta()) < 2.5 and gsf_track.pt() > 20)
                if (good_ecal_driven_seed) prompt_pur_dR.phi->fill(gsf_track.phi(), prompt_matched);

            if (abs(gsf_track.eta()) < 2.5)
                if (good_ecal_driven_seed) nonprompt_pur_dR.pt->fill(gsf_track.pt(), nonprompt_matched);
            if (gsf_track.pt() > 20)
                if (good_ecal_driven_seed) nonprompt_pur_dR.eta->fill(gsf_track.eta(), nonprompt_matched);
            if (abs(gsf_track.eta()) < 2.5 and gsf_track.pt() > 20)
                if (good_ecal_driven_seed) nonprompt_pur_dR.phi->fill(gsf_track.phi(), nonprompt_matched);
        }

        // Fake-Rate
        std::map<size_t, std::vector<int>> scl2trk; // Maps superclusters to tm-status of matched gsf-tracks
        for (const auto &gsf_track : gsf_tracks) {
            bool gsf_truth_matched = false;
            for(const auto& sim_track_idx : gsf_track.simTrkIdx()) {
                if (is_good_sim(sim_els, sim_tracks[sim_track_idx])) {
                    gsf_truth_matched = true;
                    break;
                }
            }
            int scl_idx = seeds[gsf_track.seedIdx()].sclIdx();
            scl2trk[scl_idx].push_back(gsf_truth_matched);
        }
        std::set<size_t> tm_scls; // Set of super-clusters with well matched sim-electrons
        std::set<size_t> tm_scls_no_e_match; // Set of super-clusters with well matched sim-electrons w/o energy matching
        for (const auto &scl : scls) {
            Vec4 p4 = scl_p4(scl);
            for(const auto& sim_track : sim_els) {
                if (deltaR(p4.eta(), p4.phi(), sim_track.eta(), sim_track.phi()) > 0.2) continue;
                tm_scls_no_e_match.insert(scl.idx);
                if (((p4.Et() - sim_track.pt()) / p4.Et()) > 0.1) continue;
                tm_scls.insert(scl.idx);
                break;
            }
        }
        /* cout << "scls: "; */
        /* for (const auto& scl: scls) cout << scl.idx << ", "; */
        /* cout << endl << "tm_scls: "; */
        /* for (const auto& idx: tm_scls) cout << idx << ", "; */
        /* cout << endl; */
        for (const auto &scl : scls) {
            if (scl.hoe() > hoe_cut) continue;
            float scl_pt = hypot(scl.px(), scl.py());
            float scl_eta = pseudorapidity(scl.x(), scl.y(), scl.z());
            float scl_phi = atan2(scl.y(), scl.x());
            float rad = hypot(scl.x(), scl.y());
            /* cout << "weird phi: " << scl_phi << " (" << scl.x() << ", " << scl.y() << ") - " << rad << "\n"; */
            /* if (scl_phi > 3.141592653589793238 or scl_phi < -3.141592653589793238) { */
            /*     cout << "weird phi: " << scl_phi << " (" << scl.x() << ", " << scl.y() << ")\n"; */
            /* } */

            int ntracks = scl2trk[scl.idx].size();
            int ntmtracks = std::accumulate(scl2trk[scl.idx].begin(), scl2trk[scl.idx].end(), 0);
            bool partial_fake = ntmtracks > 0 and ntracks > ntmtracks;
            bool full_fake = ntmtracks == 0 and ntracks > 0;
            /* cout << "ntracks: " << ntracks << " "; */
            /* cout << "count: " << tm_scls.count(scl.idx); */
            if (ntracks > 0) {
                partial_fake_rate_incl.pt->fill(scl_pt, partial_fake);
                partial_fake_rate_incl.eta->fill(scl_eta, partial_fake);
                partial_fake_rate_incl.phi->fill(scl_phi, partial_fake);
                if (abs(scl_eta) < 2.5)                 partial_fake_rate.pt->fill(scl_pt, partial_fake);
                if (scl_pt > 20.0)                      partial_fake_rate.eta->fill(scl_eta, partial_fake);
                if (abs(scl_eta) < 2.5 and scl_pt > 20) partial_fake_rate.phi->fill(scl_phi, partial_fake);

                full_fake_rate_incl.pt->fill(scl_pt, full_fake);
                full_fake_rate_incl.eta->fill(scl_eta, full_fake);
                full_fake_rate_incl.phi->fill(scl_phi, full_fake);
                if (abs(scl_eta) < 2.5)                 full_fake_rate.pt->fill(scl_pt, full_fake);
                if (scl_pt > 20.0)                      full_fake_rate.eta->fill(scl_eta, full_fake);
                if (abs(scl_eta) < 2.5 and scl_pt > 20) full_fake_rate.phi->fill(scl_phi, full_fake);

                if (tm_scls.count(scl.idx) == 0) {
                    clean_fake_rate_incl.pt->fill(scl_pt, full_fake);
                    clean_fake_rate_incl.eta->fill(scl_eta, full_fake);
                    clean_fake_rate_incl.phi->fill(scl_phi, full_fake);
                    if (abs(scl_eta) < 2.5)                 clean_fake_rate.pt->fill(scl_pt, full_fake);
                    if (scl_pt > 20.0)                      clean_fake_rate.eta->fill(scl_eta, full_fake);
                    if (abs(scl_eta) < 2.5 and scl_pt > 20) clean_fake_rate.phi->fill(scl_phi, full_fake);
                }
            }
            if (tm_scls.count(scl.idx) == 0) {
                fake_rate_incl.pt->fill(scl_pt, ntracks>0);
                fake_rate_incl.eta->fill(scl_eta, ntracks>0);
                fake_rate_incl.phi->fill(scl_phi, ntracks>0);
                if (abs(scl_eta) < 2.5)                 fake_rate.pt->fill(scl_pt, ntracks>0);
                if (scl_pt > 20.0)                      fake_rate.eta->fill(scl_eta, ntracks>0);
                if (abs(scl_eta) < 2.5 and scl_pt > 20) fake_rate.phi->fill(scl_phi, ntracks>0);
            }
            if (tm_scls_no_e_match.count(scl.idx) == 0) {
                fake_rate_no_e_match_incl.pt->fill(scl_pt, ntracks>0);
                fake_rate_no_e_match_incl.eta->fill(scl_eta, ntracks>0);
                fake_rate_no_e_match_incl.phi->fill(scl_phi, ntracks>0);
                if (abs(scl_eta) < 2.5)                 fake_rate_no_e_match.pt->fill(scl_pt, ntracks>0);
                if (scl_pt > 20.0)                      fake_rate_no_e_match.eta->fill(scl_eta, ntracks>0);
                if (abs(scl_eta) < 2.5 and scl_pt > 20) {
                    fake_rate_no_e_match.phi->fill(scl_phi, ntracks>0);
                    fake_rate_no_e_match_v_PU.fill(inTimePileup->get(), ntracks>0);
                }
                fake_rate_no_e_match_v_PU_inc.fill(inTimePileup->get(), ntracks>0);
            }
            /* cout << endl; */
        }

        // Hit Residuals
        for (const auto& seed : seeds) {
            if (seed.trkIdx() < 0) continue;  // require that seed produced gsf-track
            if (!is_good_seed(seed, hoe_cut)) continue;
            bool is_seed_sim_matched = false;
            for (const auto& sim_track_idx : seed.simTrkIdx()) {
                const SimTrack& sim_track = sim_tracks[sim_track_idx];
                if (is_good_sim(sim_els, sim_track) and reco_energy_consistent(sim_track, seed)) {
                    is_seed_sim_matched = true;
                    break;
                }
            }

            bool is_track_sim_matched = false;
            const auto the_track = gsf_tracks[seed.trkIdx()];
            for (const auto& sim_track_idx : the_track.simTrkIdx()) {
                const SimTrack& sim_track = sim_tracks[sim_track_idx];
                if (is_good_sim(sim_els, sim_track) and reco_energy_consistent(sim_track, the_track)) {
                    is_track_sim_matched = true;
                    break;
                }
            }

            std::vector<TMType> tm_types;
            if (is_seed_sim_matched and is_track_sim_matched) {
                tm_types = {SeedAndTrackMatched, SeedMatched, TrackMatched};
            } else if (is_seed_sim_matched) {
                tm_types = {SeedMatched};
            } else if (is_track_sim_matched) {
                tm_types = {TrackMatched};
            } else {
                tm_types = {NoMatch};
            }
            tm_corr.fill(is_seed_sim_matched, is_track_sim_matched);

            vector<int>   hits_valid;
            vector<float> hits_dRz;
            vector<float> hits_dPhi;
            if (the_track.q() == 1) {
                hits_valid = seed.isValidPos();
                hits_dRz   = seed.dRZPos();
                hits_dPhi  = seed.dPhiPos();
            } else {
                hits_valid = seed.isValidNeg();
                hits_dRz   = seed.dRZNeg();
                hits_dPhi  = seed.dPhiNeg();
            }
            const vector<int>& hits_isBarrel = seed.isBarrel();
            const vector<int>& hits_layerOrDiskNr = seed.layerOrDiskNr();

            for (const auto& tm_type : tm_types) {
                size_t nHits = hits_valid.size();
                for (size_t hit_idx = 0; hit_idx < nHits; hit_idx++) {
                    if (!hits_valid[hit_idx]) continue;
                    bool isBarrel = hits_isBarrel[hit_idx] == 1;
                    int layerOrDiskNr = hits_layerOrDiskNr[hit_idx];
                    float dRz = abs(hits_dRz[hit_idx]);
                    float dPhi = abs(hits_dPhi[hit_idx]);
                    int eta_region = get_region(static_cast<int>(hit_idx + 1), seed.eta());


                    if (is_seed_sim_matched) {
                        dRz_residuals_v_eta[hit_idx + 1]->fill(dRz, abs(seed.eta()));
                        dPhi_residuals_v_eta[hit_idx + 1]->fill(dPhi, abs(seed.eta()));

                        dRz_residuals_v_region[make_pair(hit_idx + 1, eta_region)]->fill(dRz, seed.Et());
                        dPhi_residuals_v_region[make_pair(hit_idx + 1, eta_region)]->fill(dPhi, seed.Et());
                    }

                    residuals_dRz[make_tuple(hit_idx + 1, tm_type)]->fill(dRz, seed.Et());
                    residuals_dPhi[make_tuple(hit_idx + 1, tm_type)]->fill(dPhi, seed.Et());
                    if (layerOrDiskNr == -1) continue;  // subDet is not set w/ old seeding

                    if (isBarrel)
                        hit_vs_layer_barrel.fill(layerOrDiskNr, static_cast<const int &>(hit_idx + 1));
                    else
                        hit_vs_layer_forward.fill(layerOrDiskNr, static_cast<const int &>(hit_idx + 1));

                    if (isBarrel) {
                        BPIX_residuals_dRz[make_tuple(layerOrDiskNr, hit_idx + 1, tm_type)]->fill(dRz, seed.Et());
                        BPIX_residuals_dPhi[make_tuple(layerOrDiskNr, hit_idx + 1, tm_type)]->fill(dPhi, seed.Et());
                    } else {
                        FPIX_residuals_dRz[make_tuple(layerOrDiskNr, hit_idx + 1, tm_type)]->fill(dRz, seed.Et());
                        FPIX_residuals_dPhi[make_tuple(layerOrDiskNr, hit_idx + 1, tm_type)]->fill(dPhi, seed.Et());
                    }
                }
            }
        }
    }

    tds.save_all();
}

int main(int argc, char * argv[]){
    using namespace fv_util;
    ArgParser args(argc, argv);
    if(args.cmd_option_exists("-c")) {
        init_config(args.get_cmd_option("-c"));
        args.update_config();
        if (args.cmd_option_exists("-s")) {
            the_config->update_key("silent", true);
        }
        if (args.cmd_option_exists("-b")) {
            the_config->update_key("batch", true);
        }
        init_log(LogPriority::kLogInfo);
//        gSystem->Load("libfilval.so");
        run();
    } else {
        cout << "Usage: ./" << argv[0] << " (-s) (-b) -c config_file.yaml" << endl;
    }
    return 0;
}