#include #include #include #include #include #include #include #include #include "filval/filval.hpp" #include "filval/root/filval.hpp" #include #include #include "TrackingNtuple.h" #include "analysis/obj_types.cpp" using namespace std; using namespace fv; using namespace fv::root; typedef std::tuple MatchedTrack; typedef std::pair,std::vector> 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 subdet_names = {{1, "BPIX"}, {2, "FPIX"}}; template 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 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(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 float pseudorapidity(const T &t){ return pseudorapidity(t.x, t.y, t.z); } vector seedTypes = {"initialStepSeeds", "highPtTripletStepSeeds", "mixedTripletStepSeeds", "pixelLessStepSeeds", "tripletElectronSeeds", "pixelPairElectronSeeds", "stripPairElectronSeeds"}; Value>* seeds; Value>* pixrec_hits; Value>* tracks; Value>* sim_hits; Value>* 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 HitPair; auto& find_matched_nth_hit_in_layer = func(vector, vector, vector, vector, vector, int, int, int)>("find_matched_nth_hit_in_layer", FUNC(([](const vector& seeds, const vector& pixrec_hits, const vector& sim_hits, const vector& tracks, const vector& sim_tracks, const int&& det, const int&& pix_layer, const int&& hit_number){ vector 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, bool)>("select_even_odd_ladder_blade_hit_pairs", FUNC(([](const vector& hit_pairs, const bool &&odd){ vector even_pairs; for(const HitPair &hit_pair : hit_pairs){ // loop over all seeds if(std::get(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)>("calc_dphi_v_eta", FUNC(([](const vector& hit_pairs){ vector dphis; vector etas; for(auto hit_pair : hit_pairs){ auto& pixrec_hit = std::get(hit_pair); auto& sim_hit = std::get(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>("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>("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>("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>("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>("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>("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>("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>("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>("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>("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>("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>("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>("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)>("calc_dr_v_eta", FUNC(([](const vector& hit_pairs){ vector drs; vector etas; for(auto hit_pair : hit_pairs){ auto& pixrec_hit = std::get(hit_pair); auto& sim_hit = std::get(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>("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>("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>("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>("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>("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>("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>("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>("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>("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>("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>("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>("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>("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)>("calc_dz_v_eta", FUNC(([](const vector& hit_pairs){ vector dzs; vector etas; for(auto hit_pair : hit_pairs){ auto& pixrec_hit = std::get(hit_pair); auto& sim_hit = std::get(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>("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>("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>("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>("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>("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>("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>("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>("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>("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)>("calc_drho_v_eta", FUNC(([](const vector& hit_pairs){ vector drhos; vector etas; for(auto hit_pair : hit_pairs){ auto& pixrec_hit = std::get(hit_pair); auto& sim_hit = std::get(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>("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>("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>("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>("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& 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({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; }