Adds eff/pur for dR matching

looper/CMakeLists.txt

@@ -3,7 +3,7 @@ CMAKE_MINIMUM_REQUIRED (VERSION 3.1)
 
 PROJECT (EGAMMA CXX)
 
-SET(CMAKE_CXX_STANDARD 14)
+SET(CMAKE_CXX_STANDARD 17)
 
 ADD_SUBDIRECTORY(filval)
 
@@ -15,3 +15,8 @@ ADD_EXECUTABLE(tracking_eff analysis/tracking_eff.cpp)
 # TARGET_LINK_LIBRARIES(tracking_eff filval yaml-cpp ${ROOT_LIBRARIES})
 TARGET_LINK_LIBRARIES(tracking_eff yaml-cpp ${ROOT_LIBRARIES})
 TARGET_INCLUDE_DIRECTORIES(tracking_eff PUBLIC ${ROOT_INCLUDE_DIR} ${CMAKE_SOURCE_DIR} filval/include/ filval/root/include/ filval/yaml-cpp/include/)
+
+
+ADD_EXECUTABLE(sim_track_viz analysis/sim_track_viz.cpp)
+TARGET_LINK_LIBRARIES(sim_track_viz yaml-cpp ${ROOT_LIBRARIES})
+TARGET_INCLUDE_DIRECTORIES(sim_track_viz PUBLIC ${ROOT_INCLUDE_DIR} ${CMAKE_SOURCE_DIR} filval/include/ filval/root/include/ filval/yaml-cpp/include/)

looper/analysis/config.yaml

@@ -1,8 +1,8 @@
 #max-events: 50000
 debug: false
 
-output-file: ../hists/ttbar-old-seeding.root
-source-file-key: ttbar-old-seeding-t3
+#output-file: ../hists/ttbar-old-seeding.root
+source-file-key: wide-window
 
 hoe_cut: 0.15
 

looper/analysis/sim_track_viz.cpp

@@ -0,0 +1,93 @@
+//
+// Created by caleb on 3/26/18.
+//
+
+#include <iostream>
+#include <vector>
+#include <cmath>
+#include <TSystem.h>
+
+#include "filval.hpp"
+#include "root_filval.hpp"
+
+#include "analysis/TrackingNtupleObjs.hpp"
+
+#include <memory>
+#include <iostream>
+#include <string>
+#include <cstdio>
+
+
+template<typename ... Args>
+string fmt( const std::string& format, Args ... args ) {
+    using namespace std;
+    size_t size = snprintf( nullptr, 0, format.c_str(), args ... ) + 1; // Extra space for '\0'
+    unique_ptr<char[]> buf( new char[ size ] );
+    snprintf( buf.get(), size, format.c_str(), args ... );
+    return string( buf.get(), buf.get() + size - 1 ); // We don't want the '\0' inside
+}
+
+
+SimTrackCollection sim_tracks;
+SimVertexCollection sim_vertices;
+
+void register_objects(TrackingDataSet& tds){
+    sim_tracks.init(tds);
+    sim_vertices.init(tds);
+}
+
+void run(bool silent) {
+    using namespace std;
+    using namespace fv;
+    using namespace fv_root;
+    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);
+
+
+    while (tds.next(!silent)) {
+
+
+        auto print_vtx_pos = [](int vtx_idx) {
+            const SimVertex vtx = sim_vertices[vtx_idx];
+            cout << fmt("(%.3f, %.3f, %.3f)", vtx.x(), vtx.y(), vtx.z());
+        };
+
+//        cout << tds.get_current_event();
+        for (const SimTrack sim_track : sim_tracks) {
+            if (sim_track.decayVtxIdx().size() <= 1) continue;
+            cout << fmt("  %d: %d \n", sim_track.pdgId(), sim_track.parentVtxIdx());
+            for (const auto vtx_idx : sim_track.decayVtxIdx()) {
+                cout << "    " << vtx_idx ;
+                print_vtx_pos(vtx_idx);
+                cout << ": ";
+                for (const auto simtrk_idx : sim_vertices[vtx_idx].sourceSimIdx())
+                    cout << sim_tracks[simtrk_idx].pdgId() << " ";
+                cout << " -> " << vtx_idx << ": ";
+                for (const auto simtrk_idx : sim_vertices[vtx_idx].daughterSimIdx())
+                    cout << sim_tracks[simtrk_idx].pdgId() << " ";
+                cout << endl;
+            }
+        }
+    }
+
+    tds.save_all();
+}
+
+int main(int argc, char * argv[]){
+    using namespace fv_util;
+    ArgParser args(argc, argv);
+    bool silent = args.cmd_option_exists("-s");
+    if(args.cmd_option_exists("-c")) {
+        init_config(args.get_cmd_option("-c"));
+        args.update_config();
+        init_log(LogPriority::kLogInfo);
+//        gSystem->Load("libfilval.so");
+        run(silent);
+    } else {
+        cout << "Usage: ./" << argv[0] << " (-s) -c config_file.yaml" << endl;
+    }
+    return 0;
+}
+

looper/analysis/tracking_eff.cpp

@@ -29,6 +29,14 @@ struct {
     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.
@@ -46,6 +54,11 @@ bool is_good_sim(const SimTrack& sim_track) {
     return abs(sim_track.pdgId()) == 11 and in_lum_region(vertex);
 };
 
+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;
 }
@@ -54,18 +67,20 @@ float reco_energy_rel_err(const SimTrack& sim_track, const Seed& seed) {
     return (sim_track.pt() - seed.Et()) / sim_track.pt() ;
 }
 
-bool reco_energy_consistent(const SimTrack& sim_track, const Seed& seed, float consistency_cut=0.1) {
-    return fabs(reco_energy_rel_err(sim_track, seed)) < consistency_cut;
-}
-
 float reco_energy_rel_err(const SimTrack& sim_track, const Track& track) {
     return (sim_track.pt() - track.pt()) / sim_track.pt() ;
 }
 
-bool reco_energy_consistent(const SimTrack& sim_track, const Track& track, float consistency_cut=0.1) {
+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;
 }
 
+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));
+}
+
 void run(bool silent){
     using namespace std;
     using namespace fv;
@@ -76,7 +91,7 @@ void run(bool silent){
     register_objects(tds);
 
 
-    float hoe_cut = the_config->get("hoe_cut").as<float>(999);
+    auto hoe_cut = the_config->get("hoe_cut").as<float>(999);
 
     bs = {
         tds.track_branch<float>("bsp_x"),
@@ -91,7 +106,7 @@ void run(bool silent){
         NoMatch = 0,
         SeedMatched = 1,
         TrackMatched = 2,
-        SeedAndTrackMatched = 3
+        SeedAndTrackMatched = 3,
     };
 
     std::map<std::tuple<int, int>, ContainerTH2<float>*> residuals_dRz;
@@ -197,32 +212,31 @@ void run(bool silent){
         }
     }
 
-    auto& seed_eff_v_pt  = *tds.register_container<EfficiencyContainer<float>>("seed_eff_v_pt",  THParams::lookup("eff_v_pt"));
-    auto& seed_eff_v_eta = *tds.register_container<EfficiencyContainer<float>>("seed_eff_v_eta", THParams::lookup("eff_v_eta"));
-    auto& seed_eff_v_phi = *tds.register_container<EfficiencyContainer<float>>("seed_eff_v_phi", THParams::lookup("eff_v_phi"));
-    auto& seed_eff_v_eta_pt  = *tds.register_container<EfficiencyContainer2D<float>>("seed_eff_v_eta_pt",  THParams::lookup("eff_v_eta_pt"));
-
-    auto& tracking_eff_v_pt  = *tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_pt",  THParams::lookup("eff_v_pt"));
-    auto& tracking_eff_v_eta = *tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_eta", THParams::lookup("eff_v_eta"));
-    auto& tracking_eff_v_phi = *tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_phi", THParams::lookup("eff_v_phi"));
-    auto& tracking_eff_v_eta_pt  = *tds.register_container<EfficiencyContainer2D<float>>("tracking_eff_v_eta_pt",  THParams::lookup("eff_v_eta_pt"));
-
-    auto& tracking_eff_v_pt2  = *tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_pt2",  THParams::lookup("eff_v_pt"));
-    auto& tracking_eff_v_eta2 = *tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_eta2", THParams::lookup("eff_v_eta"));
-    auto& tracking_eff_v_phi2 = *tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_phi2", THParams::lookup("eff_v_phi"));
-
-    auto& seed_pur_v_pt  = *tds.register_container<EfficiencyContainer<float>>("seed_pur_v_pt",  THParams::lookup("pur_v_pt"));
-    auto& seed_pur_v_eta = *tds.register_container<EfficiencyContainer<float>>("seed_pur_v_eta", THParams::lookup("pur_v_eta"));
-    auto& seed_pur_v_phi = *tds.register_container<EfficiencyContainer<float>>("seed_pur_v_phi", THParams::lookup("pur_v_phi"));
-
-    auto& tracking_pur_v_pt  = *tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_pt",  THParams::lookup("pur_v_pt"));
-    auto& tracking_pur_v_eta = *tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_eta", THParams::lookup("pur_v_eta"));
-    auto& tracking_pur_v_phi = *tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_phi", THParams::lookup("pur_v_phi"));
-
-    auto& tracking_pur_v_pt2  = *tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_pt2",  THParams::lookup("pur_v_pt"));
-    auto& tracking_pur_v_eta2 = *tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_eta2", THParams::lookup("pur_v_eta"));
-    auto& tracking_pur_v_phi2 = *tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_phi2", THParams::lookup("pur_v_phi"));
-
+    KinColl<EfficiencyContainer<float>> seed_eff = {
+            tds.register_container<EfficiencyContainer<float>>("seed_eff_v_pt",  THParams::lookup("eff_v_pt")),
+            tds.register_container<EfficiencyContainer<float>>("seed_eff_v_eta", THParams::lookup("eff_v_eta")),
+            tds.register_container<EfficiencyContainer<float>>("seed_eff_v_phi", THParams::lookup("eff_v_phi"))};
+    KinColl<EfficiencyContainer<float>> seed_pur = {
+            tds.register_container<EfficiencyContainer<float>>("seed_pur_v_pt",  THParams::lookup("pur_v_pt")),
+            tds.register_container<EfficiencyContainer<float>>("seed_pur_v_eta", THParams::lookup("pur_v_eta")),
+            tds.register_container<EfficiencyContainer<float>>("seed_pur_v_phi", THParams::lookup("pur_v_phi"))};
+    KinColl<EfficiencyContainer<float>> tracking_eff = {
+            tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_pt",  THParams::lookup("eff_v_pt")),
+            tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_eta", THParams::lookup("eff_v_eta")),
+            tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_phi", THParams::lookup("eff_v_phi"))};
+    KinColl<EfficiencyContainer<float>> tracking_pur = {
+            tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_pt",  THParams::lookup("pur_v_pt")),
+            tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_eta", THParams::lookup("pur_v_eta")),
+            tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_phi", THParams::lookup("pur_v_phi"))};
+
+    KinColl<EfficiencyContainer<float>> tracking_eff_dR = {
+            tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_pt_dR",  THParams::lookup("eff_v_pt")),
+            tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_eta_dR", THParams::lookup("eff_v_eta")),
+            tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_phi_dR", THParams::lookup("eff_v_phi"))};
+    KinColl<EfficiencyContainer<float>> tracking_pur_dR = {
+            tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_pt_dR",  THParams::lookup("pur_v_pt")),
+            tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_eta_dR", THParams::lookup("pur_v_eta")),
+            tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_phi_dR", THParams::lookup("pur_v_phi"))};
 
 
     auto& hit_vs_layer_barrel = *tds.register_container<ContainerTH2<int>>("hit_vs_layer_barrel", THParams::lookup("hit_vs_layer"));
@@ -251,6 +265,7 @@ void run(bool silent){
         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];
@@ -259,13 +274,12 @@ void run(bool silent){
                     break;
                 }
             }
-            seed_eff_v_eta_pt.fill(sim_track.eta(), sim_track.pt(), matched);
             if (abs(sim_track.eta()) < 2.4)
-                seed_eff_v_pt.fill(sim_track.pt(), matched);
+                seed_eff.pt->fill(sim_track.pt(), matched);
             if (sim_track.pt() > 20.0)
-                seed_eff_v_eta.fill(sim_track.eta(), matched);
+                seed_eff.eta->fill(sim_track.eta(), matched);
             if (abs(sim_track.eta()) < 2.4 and sim_track.pt() > 20)
-                seed_eff_v_phi.fill(sim_track.phi(), matched);
+                seed_eff.phi->fill(sim_track.phi(), matched);
         }
 
         // Tracking Efficiency
@@ -280,30 +294,29 @@ void run(bool silent){
                     break;
                 }
             }
-            tracking_eff_v_eta_pt.fill(sim_track.eta(), sim_track.pt(), matched);
             if (abs(sim_track.eta()) < 2.4)
-                tracking_eff_v_pt.fill(sim_track.pt(), matched);
+                tracking_eff.pt->fill(sim_track.pt(), matched);
             if (sim_track.pt() > 20.0)
-                tracking_eff_v_eta.fill(sim_track.eta(), matched);
+                tracking_eff.eta->fill(sim_track.eta(), matched);
             if (abs(sim_track.eta()) < 2.4 and sim_track.pt() > 20)
-                tracking_eff_v_phi.fill(sim_track.phi(), matched);
-
-
+                tracking_eff.phi->fill(sim_track.phi(), matched);
 
+            // dR-matching
             matched = false;
-            for (const auto& seed_idx : sim_track.seedIdx()) {
-                const Seed& seed = seeds[seed_idx];
-                if (is_good_seed(seed, hoe_cut) and seed.trkIdx()>=0 and reco_energy_consistent(sim_track, seed)) {
-                    matched=true;
+            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_v_pt2.fill(sim_track.pt(), matched);
+                tracking_eff_dR.pt->fill(sim_track.pt(), matched);
             if (sim_track.pt() > 20.0)
-                tracking_eff_v_eta2.fill(sim_track.eta(), matched);
+                tracking_eff_dR.eta->fill(sim_track.eta(), matched);
             if (abs(sim_track.eta()) < 2.4 and sim_track.pt() > 20)
-                tracking_eff_v_phi2.fill(sim_track.phi(), matched);
+                tracking_eff_dR.phi->fill(sim_track.phi(), matched);
         }
 
         // Seeding Purity
@@ -317,11 +330,11 @@ void run(bool silent){
                 }
             }
             if (abs(seed.eta()) < 2.4)
-                seed_pur_v_pt.fill(seed.pt(), match);
+                seed_pur.pt->fill(seed.pt(), match);
             if (seed.pt() > 20)
-                seed_pur_v_eta.fill(seed.eta(), match);
+                seed_pur.eta->fill(seed.eta(), match);
             if (abs(seed.eta()) < 2.4 and seed.pt() > 20)
-                seed_pur_v_phi.fill(seed.phi(), match);
+                seed_pur.phi->fill(seed.phi(), match);
         }
 
         // Tracking Purity
@@ -329,33 +342,36 @@ void run(bool silent){
             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 match = false;
+            bool matched = false;
             for (const auto& sim_track_idx : gsf_track.simTrkIdx()) {
                 if (is_good_sim(sim_tracks[sim_track_idx])) {
-                    match = true;
+                    matched = true;
                     break;
                 }
             }
             if (abs(gsf_track.eta()) < 2.4)
-                tracking_pur_v_pt.fill(gsf_track.pt(), match);
+                tracking_pur.pt->fill(gsf_track.pt(), matched);
             if (gsf_track.pt() > 20)
-                tracking_pur_v_eta.fill(gsf_track.eta(), match);
+                tracking_pur.eta->fill(gsf_track.eta(), matched);
             if (abs(gsf_track.eta()) < 2.4 and gsf_track.pt() > 20)
-                tracking_pur_v_phi.fill(gsf_track.phi(), match);
+                tracking_pur.phi->fill(gsf_track.phi(), matched);
 
-            match = false;
-            for (const auto& sim_track_idx : seed.simTrkIdx()) {
-                if (is_good_sim(sim_tracks[sim_track_idx])) {
-                    match = true;
+
+            // 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_v_pt2.fill(gsf_track.pt(), match);
-            if (gsf_track.pt() > 20)
-                tracking_pur_v_eta2.fill(gsf_track.eta(), match);
+                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_v_phi2.fill(gsf_track.phi(), match);
+                tracking_pur_dR.phi->fill(gsf_track.phi(), matched);
         }
 
         // Hit Residuals

looper/filval

@@ -1 +1 @@
-Subproject commit 9f1f43a93c2cfd7df9ec914723266ed626fe2d3b
+Subproject commit a227be04fcce25f0beb72116839bbf59af36a85e