Adds plots of residuals vs eta for all layers/sections of pixel

analysis/tracking_validation.cpp

@@ -4,7 +4,6 @@
 #include <numeric>
 #include <limits>
 #include <cmath>
-
 #include <TSystem.h>
 
 #include "filval/filval.hpp"
@@ -26,6 +25,19 @@ using namespace std;
 #define HIT_TYPE_GLUED 1
 #define HIT_TYPE_STRIP 2
 
+#define PIXEL_BARREL 1
+#define PIXEL_ENDCAP 2
+
+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));
+}
+
+typedef std::tuple<Seed, PixRecHit, Track, SimHit, SimTrack> MatchedTrack;
+typedef std::pair<std::vector<float>,std::vector<float>> pair_vec;
+
 vector<string> seedTypes =
     {"initialStepSeeds",
      "highPtTripletStepSeeds",
@@ -56,8 +68,6 @@ void setup_matched_tracks(TrackingDataSet& tds){
     // Finds pairs of (track,sim_track) where the cycle
     // track->seed->rec_hit->sim_hit->sim_track->track
     // is closed
-    typedef std::tuple<Seed, PixRecHit, Track, SimHit, SimTrack> MatchedTrack;
-
     auto& find_matched_tracks =
     func<vector<MatchedTrack>(vector<Seed>,
                               vector<PixRecHit>,
@@ -104,96 +114,166 @@ void setup_matched_tracks(TrackingDataSet& tds){
             return matched_tracks;
         })));
 
-    auto matched_tracks =
-        fv::apply(find_matched_tracks, fv::tuple(seeds, pixrec_hits, tracks, sim_hits, sim_tracks), "matched_tracks");
+    fv::apply(find_matched_tracks, fv::tuple(seeds, pixrec_hits, tracks, sim_hits, sim_tracks), "matched_tracks");
+}
 
-    auto delta_pt = fv::map(func<float(MatchedTrack)>("matched_track_delta_pt",
-        FUNC(([](const MatchedTrack& matched_track){
-            const Track& track = std::get<Track>(matched_track);
-            const SimTrack& sim_track = std::get<SimTrack>(matched_track);
-            return track.pt - sim_track.pt;
-        }))), matched_tracks, "matched_track_delta_pt");
-
-    tds.register_container<ContainerTH1Many<float>>("matched_track_delta_pt", "Matched Track delta P_t",
-                                                      delta_pt, 20, -10, 10, "P_t");
-    /* auto dphi_vs_eta_B1 = */
-    /*     fv::apply(calc_dphi_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("1", 1)), "dphi_vs_eta_B1"); */
-    /* auto dphi_vs_eta_B2 = */
-    /*     fv::apply(calc_dphi_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("2", 2)), "dphi_vs_eta_B2"); */
-    /* auto dphi_vs_eta_B3 = */
-    /*     fv::apply(calc_dphi_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("3", 3)), "dphi_vs_eta_B3"); */
-    /* auto dphi_vs_eta_B4 = */
-    /*     fv::apply(calc_dphi_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("4", 4)), "dphi_vs_eta_B4"); */
-
-
-    /* tds.register_container<ContainerTH2Many<float>>("dphi_vs_eta_B1", "dphi vs eta (BPIX L1)", */
-    /*                                                   dphi_vs_eta_B1, 20, -3, 3, 100, -0.002, 0.002, */
-    /*                                                   "eta",  "dphi"); */
-    /* tds.register_container<ContainerTH2Many<float>>("dphi_vs_eta_B2", "dphi vs eta (BPIX L2)", */
-    /*                                                   dphi_vs_eta_B2, 20, -3, 3, 100, -0.002, 0.002, */
-    /*                                                   "eta",  "dphi"); */
-    /* tds.register_container<ContainerTH2Many<float>>("dphi_vs_eta_B3", "dphi vs eta (BPIX L3)", */
-    /*                                                   dphi_vs_eta_B3, 20, -3, 3, 100, -0.002, 0.002, */
-    /*                                                   "eta",  "dphi"); */
-    /* tds.register_container<ContainerTH2Many<float>>("dphi_vs_eta_B4", "dphi vs eta (BPIX L4)", */
-    /*                                                   dphi_vs_eta_B4, 20, -3, 3, 100, -0.002, 0.002, */
-    /*                                                   "eta",  "dphi"); */
+bool in_det(const MatchedTrack &mt, int &&det, int &&layer){
+    auto& hit = std::get<PixRecHit>(mt);
+    return hit.det == det && hit.lay == layer;
+};
+
+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);
 
-/* void setup_dz_vs_eta(TrackingDataSet& tds){ */
-/*     auto& calc_dz_vs_eta = */
-/*     func<std::pair<vector<float>,vector<float>>(vector<Seed>, */
-/*                                                 vector<PixRecHit>, */
-/*                                                 vector<SimHit>, int)>("calc_dz_vs_eta", */
-/*             FUNC(([](const vector<Seed>& seeds, */
-/*                      const vector<PixRecHit> pixrec_hits, */
-/*                      const vector<SimHit> sim_hits, */
-/*                      int layer){ */
-/*                 vector<float> dz; */
-/*                 vector<float> eta; */
-/*                 for(const Seed &seed : seeds){ */
-/*                     int hitIdx, hitType; */
-/*                     for(auto tup : boost::combine(seed.hitIdx, seed.hitType)){ */
-/*                         boost::tie(hitIdx, hitType) = tup; */
-/*                         if(hitType == HIT_TYPE_PIXEL){ */
-/*                             const PixRecHit &hit = pixrec_hits[hitIdx]; */
-/*                             if(layer && hit.lay != layer) continue; */
-/*                             for(const int& simHitIdx : hit.simHitIdx){ */
-/*                                 dz.push_back(hit.z-sim_hits[simHitIdx].z); */
-/*                                 eta.push_back(seed.eta); */
-/*                             } */
-/*                         } */
-/*                     } */
-/*                 } */
-/*                 return std::make_pair(eta,dz); */
-/*             }))); */
-
-/*     auto dz_vs_eta_B1 = */
-/*     fv::apply(calc_dz_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("1", 1)), "dz_vs_eta_B1"); */
-/*     auto dz_vs_eta_B2 = */
-/*     fv::apply(calc_dz_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("2", 2)), "dz_vs_eta_B2"); */
-/*     auto dz_vs_eta_B3 = */
-/*     fv::apply(calc_dz_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("3", 3)), "dz_vs_eta_B3"); */
-/*     auto dz_vs_eta_B4 = */
-/*     fv::apply(calc_dz_vs_eta, fv::tuple(seeds, pixrec_hits, sim_hits, constant<int>("4", 4)), "dz_vs_eta_B4"); */
-
-
-/*     tds.register_container<ContainerTH2Many<float>>("dz_vs_eta_B1", "dz vs eta (BPIX L1)", */
-/*                                                       dz_vs_eta_B1, 20, -3, 3, 20, -.01, .01, */
-/*                                                       "eta",  "dz"); */
-/*     tds.register_container<ContainerTH2Many<float>>("dz_vs_eta_B2", "dz vs eta (BPIX L2)", */
-/*                                                       dz_vs_eta_B2, 20, -3, 3, 20, -.01, .01, */
-/*                                                       "eta",  "dz"); */
-/*     tds.register_container<ContainerTH2Many<float>>("dz_vs_eta_B3", "dz vs eta (BPIX L3)", */
-/*                                                       dz_vs_eta_B3, 20, -3, 3, 20, -.01, .01, */
-/*                                                       "eta",  "dz"); */
-/*     tds.register_container<ContainerTH2Many<float>>("dz_vs_eta_B4", "dz vs eta (BPIX L4)", */
-/*                                                       dz_vs_eta_B4, 20, -3, 3, 20, -.01, .01, */
-/*                                                       "eta",  "dz"); */
-
-/* } */
+    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)");
+    tds.register_container<ContainerTH2Many<float>>("matched_track_dzs_v_eta_F1", "Matched Track dz - F1",
+                                                    fv::apply(calc_dz_v_eta, fv::tuple(matched_tracks_F1)),
+                                                    100, -4, 4, 100, -.002, .002,"Eta","dz(cm)");
+    tds.register_container<ContainerTH2Many<float>>("matched_track_dzs_v_eta_F2", "Matched Track dz - F2",
+                                                    fv::apply(calc_dz_v_eta, fv::tuple(matched_tracks_F2)),
+                                                    100, -4, 4, 100, -.002, .002,"Eta","dz(cm)");
+    tds.register_container<ContainerTH2Many<float>>("matched_track_dzs_v_eta_F3", "Matched Track dz - F3",
+                                                    fv::apply(calc_dz_v_eta, fv::tuple(matched_tracks_F3)),
+                                                    100, -4, 4, 100, -.002, .002,"Eta","dz(cm)");
+}
 
 void run_analysis(const vector<fv::util::DataFileDescriptor>& dfds, const string output_filename, bool silent){
     gSystem->Load("libfilval.so");
@@ -208,10 +288,8 @@ void run_analysis(const vector<fv::util::DataFileDescriptor>& dfds, const string
     /* tds.set_max_events(10); */
     register_objects(tds);
 
-    /* setup_dz_vs_eta(tds); */
-    /* setup_dphi_vs_eta(tds); */
     setup_matched_tracks(tds);
-
+    setup_residuals(tds);
 
     tds.process(silent);
     tds.save_all();
@@ -225,7 +303,7 @@ int main(int argc, char * argv[])
     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")) {
+    }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);

filval

@@ -1 +1 @@
-Subproject commit 7f7e4233bec36e5ba9794ea2ada4dc74bef8eb38
+Subproject commit d7237e5c01868713d6b8d7491a28ffbebe2bbd8d