EGamma_ElectronTrackingValidation/analysis/tracking_validation.cpp

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

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

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

#include "TrackingNtuple.h"
#include "analysis/obj_types.cpp"

using namespace std;
using namespace fv;
using namespace fv::root;

typedef std::tuple<Seed, PixRecHit, Track, SimHit, SimTrack> MatchedTrack;
typedef std::pair<std::vector<float>,std::vector<float>> pair_vec;

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

#define PIXEL_BARREL 1
#define PIXEL_ENDCAP 2
std::map<int,string> subdet_names = {{1, "BPIX"}, {2, "FPIX"}};

template<typename A, typename B>
float displacement(const A& a, const B& b){
    return std::sqrt(pow(a.x-b.x, 2) +
                     pow(a.y-b.y, 2) +
                     pow(a.z-b.z, 2));
}

template<typename T>
float rho(const T& t){
    return std::sqrt(pow(t.x,2)+pow(t.y,2));
}

bool in_det(const MatchedTrack &mt, int &&det, int &&layer){
    auto& hit = std::get<PixRecHit>(mt);
    return hit.det == det && hit.lay == layer;
};

float pseudorapidity(float z, float r){
    float theta = atan2(r, z);
    return -log(tan(theta/2.0));
}
float pseudorapidity(float x, float y, float z){
    float r = sqrt(x*x + y*y);
    return pseudorapidity(z, r);
}
template<typename T>
float pseudorapidity(const T &t){
    return pseudorapidity(t.x, t.y, t.z);
}


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

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

void register_objects(TrackingDataSet& 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_first_hit_pairs(TrackingDataSet& tds){
    typedef std::tuple<PixRecHit, SimHit> HitPair;

    auto& find_matched_nth_hit_in_layer =
    func<vector<HitPair>(vector<Seed>,
                         vector<PixRecHit>,
                         vector<SimHit>,
                         vector<Track>,
                         vector<SimTrack>,
                         int, int, 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 int&& det,
                 const int&& pix_layer,
                 const int&& hit_number){
            vector<HitPair> matched_hits;
            for(const Track &trk : tracks){ // loop over all tracks
                const Seed &seed = seeds[trk.seedIdx];
                if(seed.hitIdx.size() <= hit_number) continue; // looking for hit_number'th hit, which this seed doesn't have
                if(seed.algoOriginal < 0 || seed.algoOriginal >= seedTypes.size()) continue;
                if(seed.hitType[hit_number] != HIT_TYPE_PIXEL) continue; // take only pixel hits for now
                const PixRecHit &rec_hit = pixrec_hits[seed.hitIdx[hit_number]];
                if(rec_hit.det == det && rec_hit.lay == pix_layer){
                    // We have the RecHit we want to work with, now find a properly* matched SimHit
                    if(rec_hit.simHitIdx.size() == 0) continue; // if rechit is matched to no simhits, give up.
                    for(const int &simTrkIdx : trk.simTrkIdx){  // loop over SimTracks matched to Track
                        for(const int& simHitIdx : sim_tracks[simTrkIdx].simHitIdx){ // loop over SimHits from SimTrack
                            if(simHitIdx == rec_hit.simHitIdx[0]) // take first matched simhit (should be the closest one)
                                matched_hits.push_back({rec_hit, sim_hits[rec_hit.simHitIdx[0]]});
                        }
                    }
                }
            }
            return matched_hits;
        })));

    auto barrel_val = constant("PIXEL_BARREL", PIXEL_BARREL);
    auto endcap_val = constant("PIXEL_ENDCAP", PIXEL_ENDCAP);
    auto first_hit  = constant("1st", 0);
    auto second_hit = constant("2nd", 1);

    // First hits on inner three bpix layers
    auto first_hits_in_B1 = fv::apply(find_matched_nth_hit_in_layer,
            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, barrel_val, constant("L1", 1), first_hit), "first_hits_in_B1");
    auto first_hits_in_B2 = fv::apply(find_matched_nth_hit_in_layer,
            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, barrel_val, constant("L2", 2), first_hit), "first_hits_in_B2");

    // Second hits on outer three bpix layers
    auto second_hits_in_B2 = fv::apply(find_matched_nth_hit_in_layer,
            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, barrel_val, constant("L2", 2), second_hit), "second_hits_in_B2");
    auto second_hits_in_B3 = fv::apply(find_matched_nth_hit_in_layer,
            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, barrel_val, constant("L3", 3), second_hit), "second_hits_in_B3");
    auto second_hits_in_B4 = fv::apply(find_matched_nth_hit_in_layer,
            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, barrel_val, constant("L4", 4), second_hit), "second_hits_in_B4");

    auto first_hits_in_F1 = fv::apply(find_matched_nth_hit_in_layer,
            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, endcap_val, constant("L1", 1), first_hit), "first_hits_in_F1");
    auto first_hits_in_F2 = fv::apply(find_matched_nth_hit_in_layer,
            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, endcap_val, constant("L2", 2), first_hit), "first_hits_in_F2");

    auto second_hits_in_F2 = fv::apply(find_matched_nth_hit_in_layer,
            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, endcap_val, constant("L2", 2), second_hit), "second_hits_in_F2");
    auto second_hits_in_F3 = fv::apply(find_matched_nth_hit_in_layer,
            fv::tuple(seeds, pixrec_hits, sim_hits, tracks, sim_tracks, endcap_val, constant("L3", 3), second_hit), "second_hits_in_F3");

    // Even vs Odd Ladders
    auto even = constant("even", false);
    auto odd = constant("odd", true);

    auto& select_even_odd_ladder_blade_hit_pairs =
    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::apply(select_even_odd_ladder_blade_hit_pairs,
            fv::tuple(first_hits_in_B1, even), "first_hits_in_B1_even_ladder");
    auto first_hits_in_B1_odd_ladder = fv::apply(select_even_odd_ladder_blade_hit_pairs,
            fv::tuple(first_hits_in_B1, odd), "first_hits_in_B1_odd_ladder");

    auto first_hits_in_B2_even_ladder = fv::apply(select_even_odd_ladder_blade_hit_pairs,
            fv::tuple(first_hits_in_B2, even), "first_hits_in_B2_even_ladder");
    auto first_hits_in_B2_odd_ladder = fv::apply(select_even_odd_ladder_blade_hit_pairs,
            fv::tuple(first_hits_in_B2, odd), "first_hits_in_B2_odd_ladder");

    //Plots for dPhi of collections defined above
    auto& calc_dphi_v_eta = func<pair_vec(vector<HitPair>)>("calc_dphi_v_eta",
        FUNC(([](const vector<HitPair>& hit_pairs){
            vector<float> dphis;
            vector<float> etas;
            for(auto hit_pair : hit_pairs){
                auto& pixrec_hit = std::get<PixRecHit>(hit_pair);
                auto& sim_hit = std::get<SimHit>(hit_pair);
                dphis.push_back(atan2(sim_hit.x, sim_hit.y) - atan2(pixrec_hit.x, pixrec_hit.y));
                etas.push_back(pseudorapidity(pixrec_hit));
            }
            return std::make_pair(etas, dphis);
        })));
    TH2Params params_dphi = {"$\\eta$",       100, -4,   4,
                             "$\\Delta \\phi$(rad)", 50,  -.0015, .0015};
    tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_B1",
                                                    fv::apply(calc_dphi_v_eta, fv::tuple(first_hits_in_B1)),
                                                    "First Hit in BPIX-L1", params_dphi);
    tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_B2",
                                                    fv::apply(calc_dphi_v_eta, fv::tuple(first_hits_in_B2)),
                                                    "First Hit in BPIX-L2", params_dphi);

    tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_B1_even_ladder",
                                                    fv::apply(calc_dphi_v_eta, fv::tuple(first_hits_in_B1_even_ladder)),
                                                    "First Hit in BPIX-L1 - Even Ladders", params_dphi);
    tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_B1_odd_ladder",
                                                    fv::apply(calc_dphi_v_eta, fv::tuple(first_hits_in_B1_odd_ladder)),
                                                    "First Hit in BPIX-L1 - Odd Ladders", params_dphi);
    tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_B2_even_ladder",
                                                    fv::apply(calc_dphi_v_eta, fv::tuple(first_hits_in_B2_even_ladder)),
                                                    "First Hit in BPIX-L2 - Even Ladders", params_dphi);
    tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_B2_odd_ladder",
                                                    fv::apply(calc_dphi_v_eta, fv::tuple(first_hits_in_B2_odd_ladder)),
                                                    "First Hit in BPIX-L2 - Odd Ladders", params_dphi);

    tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_second_hits_in_B2",
                                                    fv::apply(calc_dphi_v_eta, fv::tuple(second_hits_in_B2)),
                                                    "Second Hit in BPIX-L2", params_dphi);
    tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_second_hits_in_B3",
                                                    fv::apply(calc_dphi_v_eta, fv::tuple(second_hits_in_B3)),
                                                    "Second Hit in BPIX-L3", params_dphi);
    tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_second_hits_in_B4",
                                                    fv::apply(calc_dphi_v_eta, fv::tuple(second_hits_in_B4)),
                                                    "Second Hit in BPIX-L4", params_dphi);

    tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_F1",
                                                    fv::apply(calc_dphi_v_eta, fv::tuple(first_hits_in_F1)),
                                                    "First Hit in FPIX-L1", params_dphi);
    tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_first_hits_in_F2",
                                                    fv::apply(calc_dphi_v_eta, fv::tuple(first_hits_in_F2)),
                                                    "First Hit in FPIX-L2", params_dphi);

    tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_second_hits_in_F2",
                                                    fv::apply(calc_dphi_v_eta, fv::tuple(second_hits_in_F2)),
                                                    "Second Hit in FPIX-L2", params_dphi);
    tds.register_container<ContainerTH2Many<float>>("dphi_v_eta_second_hits_in_F3",
                                                    fv::apply(calc_dphi_v_eta, fv::tuple(second_hits_in_F3)),
                                                    "Second Hit in FPIX-L3", params_dphi);

    //Plots for dr of collections defined above
    auto& calc_dr_v_eta = func<pair_vec(vector<HitPair>)>("calc_dr_v_eta",
        FUNC(([](const vector<HitPair>& hit_pairs){
            vector<float> drs;
            vector<float> etas;
            for(auto hit_pair : hit_pairs){
                auto& pixrec_hit = std::get<PixRecHit>(hit_pair);
                auto& sim_hit = std::get<SimHit>(hit_pair);
                drs.push_back(displacement(pixrec_hit, sim_hit));
                etas.push_back(pseudorapidity(pixrec_hit));
            }
            return std::make_pair(etas, drs);
        })));
    TH2Params params_dr = {"$\\eta$",       100, -4,   4,
                           "$\\Delta \\r$(rad)", 50, 0, .006};
    tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_B1",
                                                    fv::apply(calc_dr_v_eta, fv::tuple(first_hits_in_B1)),
                                                    "First Hit in BPIX-L1", params_dr);
    tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_B2",
                                                    fv::apply(calc_dr_v_eta, fv::tuple(first_hits_in_B2)),
                                                    "First Hit in BPIX-L2", params_dr);

    tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_B1_even_ladder",
                                                    fv::apply(calc_dr_v_eta, fv::tuple(first_hits_in_B1_even_ladder)),
                                                    "First Hit in BPIX-L1 - Even Ladders", params_dr);
    tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_B1_odd_ladder",
                                                    fv::apply(calc_dr_v_eta, fv::tuple(first_hits_in_B1_odd_ladder)),
                                                    "First Hit in BPIX-L1 - Odd Ladders", params_dr);
    tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_B2_even_ladder",
                                                    fv::apply(calc_dr_v_eta, fv::tuple(first_hits_in_B2_even_ladder)),
                                                    "First Hit in BPIX-L2 - Even Ladders", params_dr);
    tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_B2_odd_ladder",
                                                    fv::apply(calc_dr_v_eta, fv::tuple(first_hits_in_B2_odd_ladder)),
                                                    "First Hit in BPIX-L2 - Odd Ladders", params_dr);

    tds.register_container<ContainerTH2Many<float>>("dr_v_eta_second_hits_in_B2",
                                                    fv::apply(calc_dr_v_eta, fv::tuple(second_hits_in_B2)),
                                                    "Second Hit in BPIX-L2", params_dr);
    tds.register_container<ContainerTH2Many<float>>("dr_v_eta_second_hits_in_B3",
                                                    fv::apply(calc_dr_v_eta, fv::tuple(second_hits_in_B3)),
                                                    "Second Hit in BPIX-L3", params_dr);
    tds.register_container<ContainerTH2Many<float>>("dr_v_eta_second_hits_in_B4",
                                                    fv::apply(calc_dr_v_eta, fv::tuple(second_hits_in_B4)),
                                                    "Second Hit in BPIX-L4", params_dr);

    tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_F1",
                                                    fv::apply(calc_dr_v_eta, fv::tuple(first_hits_in_F1)),
                                                    "First Hit in FPIX-L1", params_dr);
    tds.register_container<ContainerTH2Many<float>>("dr_v_eta_first_hits_in_F2",
                                                    fv::apply(calc_dr_v_eta, fv::tuple(first_hits_in_F2)),
                                                    "First Hit in FPIX-L2", params_dr);

    tds.register_container<ContainerTH2Many<float>>("dr_v_eta_second_hits_in_F2",
                                                    fv::apply(calc_dr_v_eta, fv::tuple(second_hits_in_F2)),
                                                    "Second Hit in FPIX-L2", params_dr);
    tds.register_container<ContainerTH2Many<float>>("dr_v_eta_second_hits_in_F3",
                                                    fv::apply(calc_dr_v_eta, fv::tuple(second_hits_in_F3)),
                                                    "Second Hit in FPIX-L3", params_dr);


    //Plots for dz of collections defined above
    auto& calc_dz_v_eta = func<pair_vec(vector<HitPair>)>("calc_dz_v_eta",
        FUNC(([](const vector<HitPair>& hit_pairs){
            vector<float> dzs;
            vector<float> etas;
            for(auto hit_pair : hit_pairs){
                auto& pixrec_hit = std::get<PixRecHit>(hit_pair);
                auto& sim_hit = std::get<SimHit>(hit_pair);
                dzs.push_back(sim_hit.z - pixrec_hit.z);
                etas.push_back(pseudorapidity(pixrec_hit));
            }
            return std::make_pair(etas, dzs);
        })));
    TH2Params params_dz = {"$\\eta$",       100, -4,   4,
                           "$\\Delta z$(rad)", 50,  -.01, .01};
    tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B1",
                                                    fv::apply(calc_dz_v_eta, fv::tuple(first_hits_in_B1)),
                                                    "First Hit in BPIX-L1", params_dz);
    tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B2",
                                                    fv::apply(calc_dz_v_eta, fv::tuple(first_hits_in_B2)),
                                                    "First Hit in BPIX-L2", params_dz);

    tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B1_even_ladder",
                                                    fv::apply(calc_dz_v_eta, fv::tuple(first_hits_in_B1_even_ladder)),
                                                    "First Hit in BPIX-L1 - Even Ladders", params_dz);
    tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B1_odd_ladder",
                                                    fv::apply(calc_dz_v_eta, fv::tuple(first_hits_in_B1_odd_ladder)),
                                                    "First Hit in BPIX-L1 - Odd Ladders", params_dz);
    tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B2_even_ladder",
                                                    fv::apply(calc_dz_v_eta, fv::tuple(first_hits_in_B2_even_ladder)),
                                                    "First Hit in BPIX-L2 - Even Ladders", params_dz);
    tds.register_container<ContainerTH2Many<float>>("dz_v_eta_first_hits_in_B2_odd_ladder",
                                                    fv::apply(calc_dz_v_eta, fv::tuple(first_hits_in_B2_odd_ladder)),
                                                    "First Hit in BPIX-L2 - Odd Ladders", params_dz);

    tds.register_container<ContainerTH2Many<float>>("dz_v_eta_second_hits_in_B2",
                                                    fv::apply(calc_dz_v_eta, fv::tuple(second_hits_in_B2)),
                                                    "Second Hit in BPIX-L2", params_dz);
    tds.register_container<ContainerTH2Many<float>>("dz_v_eta_second_hits_in_B3",
                                                    fv::apply(calc_dz_v_eta, fv::tuple(second_hits_in_B3)),
                                                    "Second Hit in BPIX-L3", params_dz);
    tds.register_container<ContainerTH2Many<float>>("dz_v_eta_second_hits_in_B4",
                                                    fv::apply(calc_dz_v_eta, fv::tuple(second_hits_in_B4)),
                                                    "Second Hit in BPIX-L4", params_dz);

    //Plots for drho of collections defined above
    auto& calc_drho_v_eta = func<pair_vec(vector<HitPair>)>("calc_drho_v_eta",
        FUNC(([](const vector<HitPair>& hit_pairs){
            vector<float> drhos;
            vector<float> etas;
            for(auto hit_pair : hit_pairs){
                auto& pixrec_hit = std::get<PixRecHit>(hit_pair);
                auto& sim_hit = std::get<SimHit>(hit_pair);
                drhos.push_back(rho(sim_hit) - rho(pixrec_hit));
                etas.push_back(pseudorapidity(pixrec_hit));
            }
            return std::make_pair(etas, drhos);
        })));
    TH2Params params_drho = {"$\\eta$",       100, -4,   4,
                             "$\\Delta \\rho$(cm)", 50,  -.01, .01};
    tds.register_container<ContainerTH2Many<float>>("drho_v_eta_first_hits_in_F1",
                                                    fv::apply(calc_drho_v_eta, fv::tuple(first_hits_in_F1)),
                                                    "First Hit in FPIX-L1", params_drho);
    tds.register_container<ContainerTH2Many<float>>("drho_v_eta_first_hits_in_F2",
                                                    fv::apply(calc_drho_v_eta, fv::tuple(first_hits_in_F2)),
                                                    "First Hit in FPIX-L2", params_drho);

    tds.register_container<ContainerTH2Many<float>>("drho_v_eta_second_hits_in_F2",
                                                    fv::apply(calc_drho_v_eta, fv::tuple(second_hits_in_F2)),
                                                    "Second Hit in FPIX-L2", params_drho);
    tds.register_container<ContainerTH2Many<float>>("drho_v_eta_second_hits_in_F3",
                                                    fv::apply(calc_drho_v_eta, fv::tuple(second_hits_in_F3)),
                                                    "Second Hit in FPIX-L3", params_drho);
}

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_first_hit_pairs(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("-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") : "";
        fv::util::DataFileDescriptor dfd(input_filename, data_label);
        run_analysis(std::vector<fv::util::DataFileDescriptor>({dfd}), output_filename, silent);
    }else {
        cout << "Usage: ./tracking_validation (-s) {-o output_filename} -F datafiles.txt" << endl;
        cout << "       ./tracking_validation (-s) {-l DATA_LABEL} {-o output_filename} -f treefile.root" << endl;
    }
    return 0;
}