EGamma_ElectronTrackingValidation/looper/analysis/tracking_eff.cpp

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

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

#include "analysis/TrackingNtupleObjs.hpp"

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

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

void register_objects(TrackingDataSet& tds){
    seeds.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;

template <typename T>
struct KinColl {
    T* pt;
    T* eta;
    T* phi;
};

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_fake(const SimTrack& sim_track) {
//    const auto& vertex = sim_vertices[sim_track.parentVtxIdx()];
//    return abs(sim_track.pdgId()) != 11 and sim_track.q() != 0 and in_lum_region(vertex);
//};

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 deltaR(float scl_eta, float scl_phi, const SimTrack& sim_track) {
    return sqrt(pow(sim_track.eta() - scl_eta, 2.0) +
                pow(sim_track.phi() - scl_phi, 2.0));
}

double deltaR(const SimTrack& sim_track, const Track& gsf_track) {
    return sqrt(pow(sim_track.eta() - gsf_track.eta(), 2.0) +
                pow(sim_track.phi() - gsf_track.phi(), 2.0));
}


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_tracks)
        if (is_good_sim(sim_track)) sim_matches[sim_track.idx] = {};
    for(const auto& gsf_track : gsf_tracks)
        gsf_matches[gsf_track.idx] = {};

    nSimTrack = sim_matches.size();
    nGSFTrack = gsf_matches.size();

    for (const auto& sim_track : sim_tracks) {
        if (!is_good_sim(sim_track)) continue;
        if (dRMatch) {
            for (const auto& gsf_track : gsf_tracks) {
                float dR = deltaR(sim_track, gsf_track);
                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];
                float dR = deltaR(sim_track, gsf_track);
                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++;
            }
        }
    }

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

    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)
            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) ? THParams::lookup("dRz_v_Et") : THParams::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) ? THParams::lookup("dPhi_v_Et") : THParams::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) ? THParams::lookup("dRz_v_Et") : THParams::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) ? THParams::lookup("dPhi_v_Et") : THParams::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) ? THParams::lookup("dPhi_v_eta") : THParams::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) ? THParams::lookup("dRz_v_eta") : THParams::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) ? THParams::lookup("dRz_v_Et") : THParams::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) ? THParams::lookup("dPhi_v_Et") : THParams::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<EfficiencyContainer<float>> seed_eff = {
            tds.register_container<EfficiencyContainer<float>>("seed_eff_v_pt",  THParams::lookup("pt")),
            tds.register_container<EfficiencyContainer<float>>("seed_eff_v_eta", THParams::lookup("eta")),
            tds.register_container<EfficiencyContainer<float>>("seed_eff_v_phi", THParams::lookup("phi"))};
    KinColl<EfficiencyContainer<float>> seed_pur = {
            tds.register_container<EfficiencyContainer<float>>("seed_pur_v_pt",  THParams::lookup("pt")),
            tds.register_container<EfficiencyContainer<float>>("seed_pur_v_eta", THParams::lookup("eta")),
            tds.register_container<EfficiencyContainer<float>>("seed_pur_v_phi", THParams::lookup("phi"))};
    KinColl<EfficiencyContainer<float>> tracking_eff = {
            tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_pt",  THParams::lookup("pt")),
            tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_eta", THParams::lookup("eta")),
            tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_phi", THParams::lookup("phi"))};
    KinColl<EfficiencyContainer<float>> tracking_pur = {
            tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_pt",  THParams::lookup("pt")),
            tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_eta", THParams::lookup("eta")),
            tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_phi", THParams::lookup("phi"))};

    KinColl<EfficiencyContainer<float>> tracking_eff_dR = {
            tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_pt_dR",  THParams::lookup("pt")),
            tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_eta_dR", THParams::lookup("eta")),
            tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_phi_dR", THParams::lookup("phi"))};
    KinColl<EfficiencyContainer<float>> tracking_pur_dR = {
            tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_pt_dR",  THParams::lookup("pt")),
            tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_eta_dR", THParams::lookup("eta")),
            tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_phi_dR", THParams::lookup("phi"))};

    KinColl<EfficiencyContainer<float>> partial_fake_rate = {
            tds.register_container<EfficiencyContainer<float>>("partial_fake_rate_v_pt",  THParams::lookup("pt")),
            tds.register_container<EfficiencyContainer<float>>("partial_fake_rate_v_eta", THParams::lookup("eta")),
            tds.register_container<EfficiencyContainer<float>>("partial_fake_rate_v_phi", THParams::lookup("phi"))};

    KinColl<EfficiencyContainer<float>> full_fake_rate = {
            tds.register_container<EfficiencyContainer<float>>("full_fake_rate_v_pt",  THParams::lookup("pt")),
            tds.register_container<EfficiencyContainer<float>>("full_fake_rate_v_eta", THParams::lookup("eta")),
            tds.register_container<EfficiencyContainer<float>>("full_fake_rate_v_phi", THParams::lookup("phi"))};

    KinColl<EfficiencyContainer<float>> clean_fake_rate = {
            tds.register_container<EfficiencyContainer<float>>("clean_fake_rate_v_pt",  THParams::lookup("pt")),
            tds.register_container<EfficiencyContainer<float>>("clean_fake_rate_v_eta", THParams::lookup("eta")),
            tds.register_container<EfficiencyContainer<float>>("clean_fake_rate_v_phi", THParams::lookup("phi"))};

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

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


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

    auto& n_seeds = *tds.register_container<ContainerTH1<size_t>>("n_seeds", THParams::lookup("n_seeds"));
    auto& n_good_sim = *tds.register_container<ContainerTH1<int>>("n_good_sim", THParams::lookup("n_seeds"));
    auto& n_gsf_tracks = *tds.register_container<ContainerTH1<int>>("n_gsf_track", THParams::lookup("n_seeds"));

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

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

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

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

    auto& sim_pt = *tds.register_container<ContainerTH1<float>>("sim_pt", THParams::lookup("pt_fine"));
    auto& sim_eta = *tds.register_container<ContainerTH1<float>>("sim_eta", THParams::lookup("eta_fine"));
    auto& sim_phi = *tds.register_container<ContainerTH1<float>>("sim_phi", THParams::lookup("phi_fine"));

    auto& reco_pt = *tds.register_container<ContainerTH1<float>>("reco_pt", THParams::lookup("pt_fine"));
    auto& reco_eta = *tds.register_container<ContainerTH1<float>>("reco_eta", THParams::lookup("eta_fine"));
    auto& reco_phi = *tds.register_container<ContainerTH1<float>>("reco_phi", THParams::lookup("phi_fine"));

    while (tds.next()) {

        n_seeds.fill(seeds.size());
        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_tracks) {
            if (!is_good_sim(sim_track)) continue;
            sim_pt.fill(sim_track.pt());
            sim_eta.fill(sim_track.eta());
            sim_phi.fill(sim_track.phi());
        }

        for (const auto& reco_track : gsf_tracks) {
            reco_pt.fill(reco_track.pt());
            reco_eta.fill(reco_track.eta());
            reco_phi.fill(reco_track.phi());
        }

        for (const auto& sim_track : sim_tracks) {
            if (!is_good_sim(sim_track)) continue;
            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_tracks) {
            if (!is_good_sim(sim_track)) continue;
            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.4)
                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.4 and sim_track.pt() > 20)
                seed_eff.phi->fill(sim_track.phi(), matched);
        }

        // Tracking Efficiency
        for (const auto& sim_track : sim_tracks) {
            if (!is_good_sim(sim_track)) continue;
            if (gsf_tracks.size() == 0) continue;
            bool matched = 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;
                    break;
                }
            }
            if (abs(sim_track.eta()) < 2.4)
                tracking_eff.pt->fill(sim_track.pt(), matched);
            if (sim_track.pt() > 20.0)
                tracking_eff.eta->fill(sim_track.eta(), matched);
            if (abs(sim_track.eta()) < 2.4 and sim_track.pt() > 20)
                tracking_eff.phi->fill(sim_track.phi(), matched);

            // dR-matching
            matched = false;
            for (const auto& gsf_track : gsf_tracks) {
                const Seed& seed = seeds[gsf_track.seedIdx()];
                double dR = deltaR(sim_track, gsf_track);
                if (is_good_seed(seed, hoe_cut) and dR < 0.2) {
                    matched = true;
                    break;
                }
            }
            if (abs(sim_track.eta()) < 2.4)
                tracking_eff_dR.pt->fill(sim_track.pt(), matched);
            if (sim_track.pt() > 20.0)
                tracking_eff_dR.eta->fill(sim_track.eta(), matched);
            if (abs(sim_track.eta()) < 2.4 and sim_track.pt() > 20)
                tracking_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_tracks[sim_track_idx])) {
                    match = true;
                    break;
                }
            }
            if (abs(seed.eta()) < 2.4)
                seed_pur.pt->fill(seed.pt(), match);
            if (seed.pt() > 20)
                seed_pur.eta->fill(seed.eta(), match);
            if (abs(seed.eta()) < 2.4 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 matched = false;
            for (const auto& sim_track_idx : gsf_track.simTrkIdx()) {
                if (is_good_sim(sim_tracks[sim_track_idx])) {
                    matched = true;
                    break;
                }
            }
            if (abs(gsf_track.eta()) < 2.4)
                tracking_pur.pt->fill(gsf_track.pt(), matched);
            if (gsf_track.pt() > 20)
                tracking_pur.eta->fill(gsf_track.eta(), matched);
            if (abs(gsf_track.eta()) < 2.4 and gsf_track.pt() > 20)
                tracking_pur.phi->fill(gsf_track.phi(), matched);

            // dR-matching
            matched = false;
            for (const auto& sim_track : sim_tracks) {
                double dR = deltaR(sim_track, gsf_track);
                if (is_good_sim(sim_track) and dR < 0.2) {
                    matched = true;
                    break;
                }
            }
            if (abs(gsf_track.eta()) < 2.4)
                tracking_pur_dR.pt->fill(gsf_track.pt(), matched);
            if (gsf_track.pt() > 20.0)
                tracking_pur_dR.eta->fill(gsf_track.eta(), matched);
            if (abs(gsf_track.eta()) < 2.4 and gsf_track.pt() > 20)
                tracking_pur_dR.phi->fill(gsf_track.phi(), 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_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
        for (const auto &scl : scls) {
            Vec4 p4 = scl_p4(scl);
            for(const auto& sim_track : sim_tracks) {
                if (!is_good_sim(sim_track)) continue;
                if (deltaR(p4.eta(), p4.phi(), sim_track) > 0.2) continue;
                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.py(), scl.px());

            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.pt->fill(scl_pt, partial_fake);
                partial_fake_rate.eta->fill(scl_eta, partial_fake);
                partial_fake_rate.phi->fill(scl_phi, partial_fake);

                full_fake_rate.pt->fill(scl_pt, full_fake);
                full_fake_rate.eta->fill(scl_eta, full_fake);
                full_fake_rate.phi->fill(scl_phi, full_fake);

                if (tm_scls.count(scl.idx) == 0) {
                    /* cout << "Filling clean_fake_rate" << endl; */
                    clean_fake_rate.pt->fill(scl_pt, full_fake);
                    clean_fake_rate.eta->fill(scl_eta, full_fake);
                    clean_fake_rate.phi->fill(scl_phi, full_fake);
                }
            }
            /* 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_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_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};
            }

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