EGamma_ElectronTrackingValidation/analysis/tracking_validation.cpp

#include <iostream>
#include <vector>
#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;

using namespace std;

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
vector<string> subdet_names = {"", "BPIX", "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.x-b.x, 2) +
                     pow(a.x-b.x, 2));
}

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



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){
    // Finds pairs of (rechit, simhit) where the rechit is the innermost hit on
    // a seed in a particular barrel layer.
    typedef std::tuple<PixRecHit, SimHit> HitPair;
    auto& matched_hits_in_layer =
    func<vector<HitPair>(vector<Seed>,
                         vector<PixRecHit>,
                         vector<SimHit>, int)>("find_matched_innermost_hits_in_layer",
        FUNC(([](const vector<Seed>& seeds,
                 const vector<PixRecHit>& pixrec_hits,
                 const vector<SimHit>& sim_hits,
                 const int&& bpix_layer){
            vector<HitPair> matched_hits;
            for(const Seed &seed : seeds){
                if( seed.algoOriginal < 0 || seed.algoOriginal >= seedTypes.size()) continue;
                for(auto tup : boost::combine(seed.hitIdx, seed.hitType)){ //for hits in track's seed
                    int hitIdx, hitType;
                    boost::tie(hitIdx, hitType) = tup;
                    if(hitType != HIT_TYPE_PIXEL) continue; // take only pixel hits for now
                    const PixRecHit &rec_hit = pixrec_hits[hitIdx];
                    if(rec_hit.det == PIXEL_BARREL && rec_hit.lay == bpix_layer){
                        if(rec_hit.simHitIdx.size() > 0){
                            matched_hits.push_back({rec_hit, sim_hits[rec_hit.simHitIdx[0]]});
                        }
                    }
                    break;
                }
            }
            return matched_hits;
        })));

    auto first_hits_in_B1 = fv::apply(matched_hits_in_layer, fv::tuple(seeds, pixrec_hits, sim_hits, constant("1", 1)), "first_hits_in_B1");
    auto first_hits_in_B2 = fv::apply(matched_hits_in_layer, fv::tuple(seeds, pixrec_hits, sim_hits, constant("2", 2)), "first_hits_in_B2");

    auto& calc_dphi = func<float(HitPair)>("calc_dphi",
        FUNC(([](const HitPair& hit_pair){
            const auto &rec_hit = std::get<PixRecHit>(hit_pair);
            const auto &sim_hit = std::get<SimHit>(hit_pair);
            return atan2(rec_hit.x, rec_hit.y) - atan2(sim_hit.x, sim_hit.y);
        })));

    tds.register_container<ContainerTH1Many<float>>("matched_hits_from_B1", "Matched Hits dPhi - B1",
                                                    fv::map(calc_dphi, first_hits_in_B1),
                                                    100, -0.001, 0.001,"dPhi(rad)");
    tds.register_container<ContainerTH1Many<float>>("matched_hits_from_B2", "Matched Hits dPhi - B2",
                                                    fv::map(calc_dphi, first_hits_in_B2),
                                                    100, -0.001, 0.001,"dPhi(rad)");
}

void setup_matched_tracks(TrackingDataSet& tds){
    // Finds sets of (seed, rechit, track, simhit sim_track) where the cycle
    // track->seed->rec_hit->sim_hit->sim_track->track
    // is closed
    auto& find_matched_tracks =
    func<vector<MatchedTrack>(vector<Seed>,
                              vector<PixRecHit>,
                              vector<Track>,
                              vector<SimHit>,
                              vector<SimTrack>)>("find_matched_tracks",
        FUNC(([](const vector<Seed>& seeds,
                 const vector<PixRecHit> pixrec_hits,
                 const vector<Track>& tracks,
                 const vector<SimHit> sim_hits,
                 const vector<SimTrack> sim_tracks){
            INFO("New event");
            INFO(boost::format("Number of tracks is %d, number of seeds is %d, number of hits %d, number of simhits %d")
                    % tracks.size() % seeds.size() % pixrec_hits.size() % sim_hits.size());
            vector<MatchedTrack> matched_tracks;
            for(const Track &track : tracks){
                const Seed& seed = seeds[track.seedIdx];
                if( seed.algoOriginal < 0 || seed.algoOriginal >= seedTypes.size()) continue;
                INFO(boost::format("       seed indx= %d   algorithm= %d   type %s\n")
                        % track.seedIdx % seed.algoOriginal % seedTypes[seed.algoOriginal]);
                for(auto tup : boost::combine(seed.hitIdx, seed.hitType)){ //for hits in track's seed
                    int hitIdx, hitType;
                    boost::tie(hitIdx, hitType) = tup;
                    if(hitType != HIT_TYPE_PIXEL) continue; // take only pixel hits for now
                    const PixRecHit &rec_hit = pixrec_hits[hitIdx];

                    TVector3 recoHitPoint;
                    recoHitPoint.SetXYZ( rec_hit.x, rec_hit.y, rec_hit.z);
                    float Rreco = recoHitPoint.Perp();
                    INFO(boost::format("            hit= %d is in %s-L%d at radius= %f\n")
                                       % hitIdx % subdet_names[rec_hit.det] % rec_hit.lay % Rreco);

                    for(const int& simHitIdx : rec_hit.simHitIdx){ // for sim-hits matching rec-hit
                        const SimHit &sim_hit = sim_hits[simHitIdx];
                        const SimTrack &sim_track = sim_tracks[sim_hit.simTrkIdx];
                        for(const int &trkIdx : sim_track.trkIdx){ // for tracks matched to this sim_track
                            if(trkIdx == track.idx){  // sim_track matches with our original track
                                matched_tracks.push_back({seed, rec_hit, track, sim_hit, sim_track});
                                break;
                            }
                        }
                    }
                }
            }
            return matched_tracks;
        })));

    fv::apply(find_matched_tracks, fv::tuple(seeds, pixrec_hits, tracks, sim_hits, sim_tracks), "matched_tracks");
}

void setup_residuals(TrackingDataSet &tds){


    auto matched_tracks = lookup<vector<MatchedTrack>>("matched_tracks");

    // Break matched_tracks into catagories based on Barrel/Endcap and layer
    auto matched_tracks_B1 = filter(func<bool(MatchedTrack)>("matched_tracks_B1",
        FUNC(([](const MatchedTrack& matched_track){
            return in_det(matched_track, PIXEL_BARREL, 1);
        }))), matched_tracks);
    auto matched_tracks_B2 = filter(func<bool(MatchedTrack)>("matched_tracks_B2",
        FUNC(([](const MatchedTrack& matched_track){
            return in_det(matched_track, PIXEL_BARREL, 2);
        }))), matched_tracks);
    auto matched_tracks_B3 = filter(func<bool(MatchedTrack)>("matched_tracks_B3",
        FUNC(([](const MatchedTrack& matched_track){
            return in_det(matched_track, PIXEL_BARREL, 3);
        }))), matched_tracks);
    auto matched_tracks_B4 = filter(func<bool(MatchedTrack)>("matched_tracks_B4",
        FUNC(([](const MatchedTrack& matched_track){
            return in_det(matched_track, PIXEL_BARREL, 4);
        }))), matched_tracks);

    auto matched_tracks_F1 = filter(func<bool(MatchedTrack)>("matched_tracks_F1",
        FUNC(([](const MatchedTrack& matched_track){
            return in_det(matched_track, PIXEL_ENDCAP, 1);
        }))), matched_tracks);
    auto matched_tracks_F2 = filter(func<bool(MatchedTrack)>("matched_tracks_F2",
        FUNC(([](const MatchedTrack& matched_track){
            return in_det(matched_track, PIXEL_ENDCAP, 2);
        }))), matched_tracks);
    auto matched_tracks_F3 = filter(func<bool(MatchedTrack)>("matched_tracks_F3",
        FUNC(([](const MatchedTrack& matched_track){
            return in_det(matched_track, PIXEL_ENDCAP, 3);
        }))), matched_tracks);


    auto& calc_residuals_v_eta = func<pair_vec(vector<MatchedTrack>)>("matched_track_hit_residuals",
        FUNC(([](const vector<MatchedTrack>& matched_tracks){
            vector<float> residuals;
            vector<float> etas;
            for(auto matched_track : matched_tracks){
                auto& pixrec_hit = std::get<PixRecHit>(matched_track);
                auto& sim_hit = std::get<SimHit>(matched_track);
                auto& sim_track = std::get<SimTrack>(matched_track);
                residuals.push_back(displacement(pixrec_hit, sim_hit));
                etas.push_back(sim_track.eta);
            }
            return std::make_pair(etas, residuals);
        })));
    tds.register_container<ContainerTH2Many<float>>("matched_track_residuals_v_eta_all", "Matched Track Residuals - All",
                                                    fv::apply(calc_residuals_v_eta, fv::tuple(matched_tracks)),
                                                    100, -4, 4, 100, 0, .01,"Eta","Residuals(cm)");
    tds.register_container<ContainerTH2Many<float>>("matched_track_residuals_v_eta_B1", "Matched Track Residuals - B1",
                                                    fv::apply(calc_residuals_v_eta, fv::tuple(matched_tracks_B1)),
                                                    100, -4, 4, 100, 0, .01,"Eta","Residuals(cm)");
    tds.register_container<ContainerTH2Many<float>>("matched_track_residuals_v_eta_B2", "Matched Track Residuals - B2",
                                                    fv::apply(calc_residuals_v_eta, fv::tuple(matched_tracks_B2)),
                                                    100, -4, 4, 100, 0, .01,"Eta","Residuals(cm)");
    tds.register_container<ContainerTH2Many<float>>("matched_track_residuals_v_eta_B3", "Matched Track Residuals - B3",
                                                    fv::apply(calc_residuals_v_eta, fv::tuple(matched_tracks_B3)),
                                                    100, -4, 4, 100, 0, .01,"Eta","Residuals(cm)");
    tds.register_container<ContainerTH2Many<float>>("matched_track_residuals_v_eta_B4", "Matched Track Residuals - B4",
                                                    fv::apply(calc_residuals_v_eta, fv::tuple(matched_tracks_B4)),
                                                    100, -4, 4, 100, 0, .01,"Eta","Residuals(cm)");
    tds.register_container<ContainerTH2Many<float>>("matched_track_residuals_v_eta_F1", "Matched Track Residuals - F1",
                                                    fv::apply(calc_residuals_v_eta, fv::tuple(matched_tracks_F1)),
                                                    100, -4, 4, 100, 0, .01,"Eta","Residuals(cm)");
    tds.register_container<ContainerTH2Many<float>>("matched_track_residuals_v_eta_F2", "Matched Track Residuals - F2",
                                                    fv::apply(calc_residuals_v_eta, fv::tuple(matched_tracks_F2)),
                                                    100, -4, 4, 100, 0, .01,"Eta","Residuals(cm)");
    tds.register_container<ContainerTH2Many<float>>("matched_track_residuals_v_eta_F3", "Matched Track Residuals - F3",
                                                    fv::apply(calc_residuals_v_eta, fv::tuple(matched_tracks_F3)),
                                                    100, -4, 4, 100, 0, .01,"Eta","Residuals(cm)");

    auto& calc_dphi_v_eta = func<pair_vec(vector<MatchedTrack>)>("matched_track_hit_dphis",
        FUNC(([](const vector<MatchedTrack>& matched_tracks){
            vector<float> dphis;
            vector<float> etas;
            for(auto matched_track : matched_tracks){
                auto& pixrec_hit = std::get<PixRecHit>(matched_track);
                auto& sim_hit = std::get<SimHit>(matched_track);
                auto& sim_track = std::get<SimTrack>(matched_track);
                dphis.push_back(atan2(sim_hit.x, sim_hit.y) - atan2(pixrec_hit.x, pixrec_hit.y));
                etas.push_back(sim_track.eta);
            }
            return std::make_pair(etas, dphis);
        })));
    tds.register_container<ContainerTH2Many<float>>("matched_track_dphis_v_eta_all", "Matched Track dphis - All",
                                                    fv::apply(calc_dphi_v_eta, fv::tuple(matched_tracks)),
                                                    100, -4, 4, 100, -.002, .002,"Eta","Residuals(rad)");
    tds.register_container<ContainerTH2Many<float>>("matched_track_dphis_v_eta_B1",  "Matched Track dphis - B1",
                                                    fv::apply(calc_dphi_v_eta, fv::tuple(matched_tracks_B1)),
                                                    100, -4, 4, 100, -.002, .002,"Eta","dphis(rad)");
    tds.register_container<ContainerTH2Many<float>>("matched_track_dphis_v_eta_B2",  "Matched Track dphis - B2",
                                                    fv::apply(calc_dphi_v_eta, fv::tuple(matched_tracks_B2)),
                                                    100, -4, 4, 100, -.002, .002,"Eta","dphis(rad)");
    tds.register_container<ContainerTH2Many<float>>("matched_track_dphis_v_eta_B3",  "Matched Track dphis - B3",
                                                    fv::apply(calc_dphi_v_eta, fv::tuple(matched_tracks_B3)),
                                                    100, -4, 4, 100, -.002, .002,"Eta","dphis(rad)");
    tds.register_container<ContainerTH2Many<float>>("matched_track_dphis_v_eta_B4",  "Matched Track dphis - B4",
                                                    fv::apply(calc_dphi_v_eta, fv::tuple(matched_tracks_B4)),
                                                    100, -4, 4, 100, -.002, .002,"Eta","dphis(rad)");
    tds.register_container<ContainerTH2Many<float>>("matched_track_dphis_v_eta_F1",  "Matched Track dphis - F1",
                                                    fv::apply(calc_dphi_v_eta, fv::tuple(matched_tracks_F1)),
                                                    100, -4, 4, 100, -.002, .002,"Eta","dphis(rad)");
    tds.register_container<ContainerTH2Many<float>>("matched_track_dphis_v_eta_F2",  "Matched Track dphis - F2",
                                                    fv::apply(calc_dphi_v_eta, fv::tuple(matched_tracks_F2)),
                                                    100, -4, 4, 100, -.002, .002,"Eta","dphis(rad)");
    tds.register_container<ContainerTH2Many<float>>("matched_track_dphis_v_eta_F3",  "Matched Track dphis - F3",
                                                    fv::apply(calc_dphi_v_eta, fv::tuple(matched_tracks_F3)),
                                                    100, -4, 4, 100, -.002, .002,"Eta","dphis(rad)");


    auto& calc_dz_v_eta = func<pair_vec(vector<MatchedTrack>)>("matched_track_hit_dphis",
        FUNC(([](const vector<MatchedTrack>& matched_tracks){
            vector<float> dzs;
            vector<float> etas;
            for(auto matched_track : matched_tracks){
                auto& pixrec_hit = std::get<PixRecHit>(matched_track);
                auto& sim_hit = std::get<SimHit>(matched_track);
                auto& sim_track = std::get<SimTrack>(matched_track);
                dzs.push_back(sim_hit.z - pixrec_hit.z);
                etas.push_back(sim_track.eta);
            }
            return std::make_pair(etas, dzs);
        })));
    tds.register_container<ContainerTH2Many<float>>("matched_track_dzs_v_eta_all", "Matched Track dz - All",
                                                    fv::apply(calc_dz_v_eta, fv::tuple(matched_tracks)),
                                                    100, -4, 4, 100, -.002, .002,"Eta","dz(cm)");
    tds.register_container<ContainerTH2Many<float>>("matched_track_dzs_v_eta_B1", "Matched Track dz - B1",
                                                    fv::apply(calc_dz_v_eta, fv::tuple(matched_tracks_B1)),
                                                    100, -4, 4, 100, -.002, .002,"Eta","dz(cm)");
    tds.register_container<ContainerTH2Many<float>>("matched_track_dzs_v_eta_B2", "Matched Track dz - B2",
                                                    fv::apply(calc_dz_v_eta, fv::tuple(matched_tracks_B2)),
                                                    100, -4, 4, 100, -.002, .002,"Eta","dz(cm)");
    tds.register_container<ContainerTH2Many<float>>("matched_track_dzs_v_eta_B3", "Matched Track dz - B3",
                                                    fv::apply(calc_dz_v_eta, fv::tuple(matched_tracks_B3)),
                                                    100, -4, 4, 100, -.002, .002,"Eta","dz(cm)");
    tds.register_container<ContainerTH2Many<float>>("matched_track_dzs_v_eta_B4", "Matched Track dz - B4",
                                                    fv::apply(calc_dz_v_eta, fv::tuple(matched_tracks_B4)),
                                                    100, -4, 4, 100, -.002, .002,"Eta","dz(cm)");
    auto& calc_dr_v_eta = func<pair_vec(vector<MatchedTrack>)>("matched_track_hit_drs",
        FUNC(([](const vector<MatchedTrack>& matched_tracks){
            vector<float> drs;
            vector<float> etas;
            for(auto matched_track : matched_tracks){
                auto& rec_hit = std::get<PixRecHit>(matched_track);
                auto& sim_hit = std::get<SimHit>(matched_track);
                float sim_hit_r = sqrt(pow(sim_hit.x, 2) + pow(sim_hit.y, 2));
                float rec_hit_r = sqrt(pow(rec_hit.x, 2) + pow(rec_hit.y, 2));
                drs.push_back(sim_hit_r - rec_hit_r);
                etas.push_back(std::get<SimTrack>(matched_track).eta);
            }
            return std::make_pair(etas, drs);
        })));
    tds.register_container<ContainerTH2Many<float>>("matched_track_drs_v_eta_F1", "Matched Track dz - F1",
                                                    fv::apply(calc_dr_v_eta, fv::tuple(matched_tracks_F1)),
                                                    100, -4, 4, 100, -.002, .002,"Eta","dr(cm)");
    tds.register_container<ContainerTH2Many<float>>("matched_track_drs_v_eta_F2", "Matched Track dz - F2",
                                                    fv::apply(calc_dr_v_eta, fv::tuple(matched_tracks_F2)),
                                                    100, -4, 4, 100, -.002, .002,"Eta","dr(cm)");
    tds.register_container<ContainerTH2Many<float>>("matched_track_drs_v_eta_F3", "Matched Track dz - F3",
                                                    fv::apply(calc_dr_v_eta, fv::tuple(matched_tracks_F3)),
                                                    100, -4, 4, 100, -.002, .002,"Eta","dr(cm)");
}

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::kLogDebug);

    TrackingDataSet tds(output_filename, dfds, "trackingNtuple/tree");
    /* tds.set_max_events(10); */
    register_objects(tds);

    setup_matched_tracks(tds);
    setup_residuals(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;
    }
}