Przeglądaj źródła

Implements additional efficiency/purity measurements for electron seeds

Caleb Fangmeier 7 lat temu
rodzic
commit
a5fa101f1c

+ 1 - 0
.gitignore

@@ -1,6 +1,7 @@
 
 legacy/
 data/
+output/
 build/
 env/
 analysis/output/

+ 11 - 0
.idea/EGamma.iml

@@ -0,0 +1,11 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<module type="PYTHON_MODULE" version="4">
+  <component name="NewModuleRootManager">
+    <content url="file://$MODULE_DIR$" />
+    <orderEntry type="inheritedJdk" />
+    <orderEntry type="sourceFolder" forTests="false" />
+  </component>
+  <component name="TestRunnerService">
+    <option name="PROJECT_TEST_RUNNER" value="Unittests" />
+  </component>
+</module>

+ 4 - 0
.idea/misc.xml

@@ -0,0 +1,4 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="ProjectRootManager" version="2" project-jdk-name="Python 3.6 (TTTT2)" project-jdk-type="Python SDK" />
+</project>

+ 8 - 0
.idea/modules.xml

@@ -0,0 +1,8 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="ProjectModuleManager">
+    <modules>
+      <module fileurl="file://$PROJECT_DIR$/.idea/EGamma.iml" filepath="$PROJECT_DIR$/.idea/EGamma.iml" />
+    </modules>
+  </component>
+</project>

+ 6 - 0
.idea/vcs.xml

@@ -0,0 +1,6 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="VcsDirectoryMappings">
+    <mapping directory="$PROJECT_DIR$" vcs="Git" />
+  </component>
+</project>

+ 0 - 8
analysis/TrackingNtuple928.h

@@ -338,8 +338,6 @@ public :
    vector<float>   *scl_y;
    vector<float>   *scl_z;
    vector<vector<int> > *scl_charge;
-   vector<vector<int> > *scl_lay1;
-   vector<vector<int> > *scl_lay2;
    vector<vector<int> > *scl_subDet1;
    vector<vector<int> > *scl_subDet2;
    vector<vector<float> > *scl_dRz1;
@@ -665,8 +663,6 @@ public :
    TBranch        *b_scl_y;   //!
    TBranch        *b_scl_z;   //!
    TBranch        *b_scl_charge;   //!
-   TBranch        *b_scl_lay1;   //!
-   TBranch        *b_scl_lay2;   //!
    TBranch        *b_scl_subDet1;   //!
    TBranch        *b_scl_subDet2;   //!
    TBranch        *b_scl_dRz1;   //!
@@ -1044,8 +1040,6 @@ void TrackingNtuple::Init(TTree *tree)
    scl_y = 0;
    scl_z = 0;
    scl_charge = 0;
-   scl_lay1 = 0;
-   scl_lay2 = 0;
    scl_subDet1 = 0;
    scl_subDet2 = 0;
    scl_dRz1 = 0;
@@ -1375,8 +1369,6 @@ void TrackingNtuple::Init(TTree *tree)
    fChain->SetBranchAddress("scl_y", &scl_y, &b_scl_y);
    fChain->SetBranchAddress("scl_z", &scl_z, &b_scl_z);
    fChain->SetBranchAddress("scl_charge", &scl_charge, &b_scl_charge);
-   fChain->SetBranchAddress("scl_lay1", &scl_lay1, &b_scl_lay1);
-   fChain->SetBranchAddress("scl_lay2", &scl_lay2, &b_scl_lay2);
    fChain->SetBranchAddress("scl_subDet1", &scl_subDet1, &b_scl_subDet1);
    fChain->SetBranchAddress("scl_subDet2", &scl_subDet2, &b_scl_subDet2);
    fChain->SetBranchAddress("scl_dRz1", &scl_dRz1, &b_scl_dRz1);

+ 69 - 5
analysis/config.yaml

@@ -1,4 +1,49 @@
-max-events: 100
+max-events: 0
+
+#output-file: old_seeding.root
+#source-files:
+#    - filename: ../data/old_seeding/trackingNtuple_1.root
+#    - filename: ../data/old_seeding/trackingNtuple_2.root
+#    - filename: ../data/old_seeding/trackingNtuple_3.root
+#    - filename: ../data/old_seeding/trackingNtuple_4.root
+#    - filename: ../data/old_seeding/trackingNtuple_5.root
+#    - filename: ../data/old_seeding/trackingNtuple_6.root
+#    - filename: ../data/old_seeding/trackingNtuple_7.root
+#    - filename: ../data/old_seeding/trackingNtuple_8.root
+#    - filename: ../data/old_seeding/trackingNtuple_9.root
+#    - filename: ../data/old_seeding/trackingNtuple_10.root
+#    - filename: ../data/old_seeding/trackingNtuple_11.root
+#    - filename: ../data/old_seeding/trackingNtuple_12.root
+#    - filename: ../data/old_seeding/trackingNtuple_13.root
+#    - filename: ../data/old_seeding/trackingNtuple_14.root
+#    - filename: ../data/old_seeding/trackingNtuple_15.root
+#    - filename: ../data/old_seeding/trackingNtuple_16.root
+#    - filename: ../data/old_seeding/trackingNtuple_17.root
+#    - filename: ../data/old_seeding/trackingNtuple_18.root
+#    - filename: ../data/old_seeding/trackingNtuple_19.root
+
+output-file: new_seeding.root
+source-files:
+    - filename: ../data/new_seeding/trackingNtuple_1.root
+    - filename: ../data/new_seeding/trackingNtuple_2.root
+    - filename: ../data/new_seeding/trackingNtuple_3.root
+    - filename: ../data/new_seeding/trackingNtuple_4.root
+    - filename: ../data/new_seeding/trackingNtuple_5.root
+    - filename: ../data/new_seeding/trackingNtuple_6.root
+    - filename: ../data/new_seeding/trackingNtuple_7.root
+    - filename: ../data/new_seeding/trackingNtuple_8.root
+    - filename: ../data/new_seeding/trackingNtuple_9.root
+    - filename: ../data/new_seeding/trackingNtuple_10.root
+    - filename: ../data/new_seeding/trackingNtuple_11.root
+    - filename: ../data/new_seeding/trackingNtuple_12.root
+    - filename: ../data/new_seeding/trackingNtuple_13.root
+    - filename: ../data/new_seeding/trackingNtuple_14.root
+    - filename: ../data/new_seeding/trackingNtuple_15.root
+    - filename: ../data/new_seeding/trackingNtuple_16.root
+    - filename: ../data/new_seeding/trackingNtuple_17.root
+    - filename: ../data/new_seeding/trackingNtuple_18.root
+    - filename: ../data/new_seeding/trackingNtuple_19.root
+
 hist-params:
     sc_dphi2_v_dz1_L1_L2:
         label_x: "$\\Delta \\phi_2$(rad)"
@@ -76,19 +121,38 @@ hist-params:
         label_y: ""
     eff_v_pt:
         label_x: "eff_v_pt"
-        nbins: 100
+        nbins: 20
         low: 0
-        high: 500
+        high: 300
         label_y: ""
     eff_v_eta:
         label_x: "eff_v_eta"
-        nbins: 100
+        nbins: 20
         low: -4.0
         high: 4.0
         label_y: ""
     eff_v_phi:
         label_x: "eff_v_phi"
-        nbins: 100
+        nbins: 20
+        low: -3.2
+        high: 3.2
+        label_y: ""
+
+    pur_v_pt:
+        label_x: "pur_v_pt"
+        nbins: 20
+        low: 0
+        high: 300
+        label_y: ""
+    pur_v_eta:
+        label_x: "pur_v_eta"
+        nbins: 20
+        low: -4.0
+        high: 4.0
+        label_y: ""
+    pur_v_phi:
+        label_x: "pur_v_phi"
+        nbins: 20
         low: -3.2
         high: 3.2
         label_y: ""

+ 0 - 69
analysis/templates/dashboard.htm

@@ -1,69 +0,0 @@
-<!DOCTYPE html>
-<html lang="en">
-<head>
-  <title>EGamma Dashboard</title>
-  <meta charset="utf-8">
-  <meta name="viewport" content="width=device-width, initial-scale=2">
-  <link rel="stylesheet" href="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.7/css/bootstrap.min.css">
-  <script src="https://ajax.googleapis.com/ajax/libs/jquery/3.2.1/jquery.min.js"></script>
-  <script src="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.7/js/bootstrap.min.js"></script>
-
-<script type="text/x-mathjax-config">
-MathJax.Hub.Config({
-  config: ["MMLorHTML.js"],
-  jax: ["input/TeX", "output/HTML-CSS", "output/NativeMML"],
-  extensions: ["MathMenu.js", "MathZoom.js"]
-});
-</script>
-<script type="text/javascript" src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS_HTML-full"> </script>
-</head>
-<body>
-<div class="container-fluid">
-{% for r, plot_row in enumerate(plots) %}
-  <div class="row">
-  {% for c, plot in enumerate(plot_row) %}
-    <div class="col-md-4">
-      <div class="well">
-    {% if plot.title %}
-        <h3 class="text-center">{{ plot.title }}</h3>
-    {% endif %}
-        <a href="{{ outdir+plot.name+".png" }}"><img src="{{ outdir+plot.name+".png" }}" style="width:100%"></a>
-        <div class="caption">
-          <p class="text-center"> {{ plot.name }} </p>
-          <div class="panel-group" id="accordion">
-    {% for id, (i,j) in enumerate(plot.docs.keys()) %}
-            <div class="panel-heading">
-              <h4 class="panel-title">
-                <button data-toggle="collapse" data-parent="#accordion" class="btn btn-info" href="#collapse{{r}}-{{c}}-{{id}}">
-                  Plot at ({{i+1}}, {{j+1}})</button>
-              </h4>
-            </div>
-            <div id="collapse{{r}}-{{c}}-{{id}}" class="panel-collapse collapse">
-              <div class="panel-body">
-      {% for doc, argdict in zip(plot.docs[(i,j)], plot.argdicts[(i,j)]) %}
-                <p class="text-left">{{ doc|safe }}</p>
-                <hr>
-                <p class="text-left"><strong>Plot Arguments</strong></p>
-                <table class="table table-hover">
-                  <tbody>
-        {% for key, val in argdict.items() %}
-                    <tr>
-                      <td>{{ key }}</td> <td>{{ val }}</td>
-                    </tr>
-        {% endfor %}
-                  </tbody>
-                </table>
-      {% endfor %}
-              </div>
-            </div>
-    {% endfor %}
-          </div>
-        </div>
-      </div>
-    </div>
-  {% endfor %}
-  </div>
-{% endfor %}
-</div>
-</body>
-</html>

+ 55 - 29
analysis/tracking_eff.cpp

@@ -83,16 +83,7 @@ float rad(const float& x, const float& y) {
 }
 
 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(11); */
     register_objects(tds);
 
     // Indices of prompt gen electrons
@@ -103,7 +94,8 @@ void run_analysis(const vector<fv::util::DataFileDescriptor>& dfds, const string
                 if (abs(sim_track.pdgId()) != 11) continue;
                 const SimVertex& sim_vertex = sim_vertices[sim_track.parentVtxIdx()];
                 float rho = rad(sim_vertex.x(), sim_vertex.y());
-                if (rho > 1) continue;
+                float z = abs(sim_vertex.z());
+                if (rho > 1 or z > 10) continue;
                 idxs.push_back(sim_track.idx);
             }
             return idxs;
@@ -151,12 +143,50 @@ void run_analysis(const vector<fv::util::DataFileDescriptor>& dfds, const string
     tds.register_container<EfficiencyContainer<size_t>>("track_eff_v_phi",
         prompt_electrons, TH1Params::lookup("eff_v_phi"), nullptr, &with_track, &e_phi);
 
-    /* std::ofstream printout("track_eff.txt"); */
-    /* tds.register_container<Printer<vector<size_t>>>("prompt_electron_idx", prompt_electrons, &printout); */
-    /* tds.register_container<Printer<vector<size_t>>>("with track:        ", apply(with_seed, prompt_electrons), &printout); */
+    // Indices of ecal-driven seeds
+    auto ecal_seeds = func_value<vector<size_t>>("ecal_seeds",
+        FUNC(([](){
+            vector<size_t> idxs;
+            for(const auto& seed : seeds) {
+                if (seed.ecalDriven())
+                    idxs.push_back(seed.idx);
+            }
+            return idxs;
+        })));
 
-    /* tds.register_container<ContainerTH1Many<float>>("dxys", sim_tracks.val_pca_dxy, "dxys", */
-    /*                                                 TH1Params::lookup("dxys")); */
+    function<bool(size_t)> with_sim_track = [](const size_t& seed_idx) {
+        for (const int& simTrkIdx : seeds[seed_idx].simTrkIdx()) {
+            const SimTrack& sim_track = sim_tracks[simTrkIdx];
+            if (abs(sim_track.pdgId()) != 11) continue;
+            const SimVertex& sim_vertex = sim_vertices[sim_track.parentVtxIdx()];
+            float rho = rad(sim_vertex.x(), sim_vertex.y());
+            float z = abs(sim_vertex.z());
+            if (rho > 1 or z > 10) continue;
+            return true;
+        }
+        return false;
+    };
+
+    function<float(size_t)> seed_pt = [](const size_t& seed_idx) {
+        return seeds[seed_idx].pt();
+    };
+
+    function<float(size_t)> seed_eta = [](const size_t& seed_idx) {
+        return seeds[seed_idx].eta();
+    };
+
+    function<float(size_t)> seed_phi = [](const size_t& seed_idx) {
+        return seeds[seed_idx].phi();
+    };
+
+    tds.register_container<EfficiencyContainer<size_t>>("seed_pur_v_pt",
+        ecal_seeds, TH1Params::lookup("pur_v_pt"), nullptr, &with_sim_track, &seed_pt);
+
+    tds.register_container<EfficiencyContainer<size_t>>("seed_pur_v_eta",
+        ecal_seeds, TH1Params::lookup("pur_v_eta"), nullptr, &with_sim_track, &seed_eta);
+
+    tds.register_container<EfficiencyContainer<size_t>>("seed_pur_v_phi",
+        ecal_seeds, TH1Params::lookup("pur_v_phi"), nullptr, &with_sim_track, &seed_phi);
 
     tds.process(silent);
     tds.save_all();
@@ -165,22 +195,18 @@ void run_analysis(const vector<fv::util::DataFileDescriptor>& dfds, const string
 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("-c")){
+    if(args.cmdOptionExists("-c")) {
         fv::util::init_config(args.getCmdOption("-c"));
-    }
-    if(args.cmdOptionExists("-F")){
-        auto file_list = fv::util::read_input_list(args.getCmdOption("-F"));
+        auto output_filename = fv::util::the_config->get_output_filename();
+        if (output_filename == "") return -1;
+
+        fv::util::init_log(fv::util::LogPriority::kLogInfo);
+        gSystem->Load("libfilval.so");
+
+        auto file_list = fv::util::the_config->get_source_files();
         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") : "";
-        string data_category = args.cmdOptionExists("-c") ? args.getCmdOption("-c") : "";
-        fv::util::DataFileDescriptor dfd(input_filename, data_label, data_category);
-        run_analysis(std::vector<fv::util::DataFileDescriptor>({dfd}), output_filename, silent);
-    }else {
-        cout << "Usage: ./tracking_eff (-s) {-o output_filename} -F datafiles.txt" << endl;
-        cout << "       ./tracking_eff (-s) {-l DATA_LABEL} {-c DATA_CATEGORY} {-o output_filename} -f treefile.root" << endl;
+    } else {
+        cout << "Usage: ./tracking_eff (-s) -c config_file.yaml" << endl;
     }
     return 0;
 }

+ 0 - 7
analysis/tracking_validation.cpp

@@ -861,10 +861,3 @@ int main(int argc, char * argv[]){
     }
     return 0;
 }
-
-/* int main(int argc, char * argv[]){ */
-/*     using namespace fv::util; */
-/*     init_config(argv[1]); */
-/*     int max_events = the_config.get("max-events").as<int>(); */
-/*     cout << max_events << endl; */
-/* } */

+ 120 - 0
eff_plots.py

@@ -0,0 +1,120 @@
+#!/usr/bin/env python
+import numpy as np
+import matplotlib.pyplot as plt
+
+from filval.result_set import ResultSet
+from filval.histogram_utils import hist, hist_add, hist_normalize, hist_scale
+from filval.plotter import (decl_plot, render_plots, hist_plot, hist_plot_stack, Plot, generate_dashboard)
+
+
+@decl_plot
+def plot_seed_eff(old_seeds, new_seeds):
+    r"""## Seeding Efficiency
+
+    The proportion of gen-level electrons origination in $\rho<1$cm and $|z|<10$cm from the beam spot that have
+    an associated Seed, matched via rechit-simhit associations in the pixel detector.
+    """
+    _, ((ax_pt, ax_eta), (ax_phi, _)) = plt.subplots(2, 2)
+
+    errors = True
+    plt.sca(ax_pt)
+    hist_plot(hist(new_seeds.seed_eff_v_pt), include_errors=errors, label='New')
+    hist_plot(hist(old_seeds.seed_eff_v_pt), include_errors=errors, title='Efficiency vs. Pt', label='Old')
+    plt.xlabel("Pt")
+    plt.ylim((0, 1.1))
+
+    plt.sca(ax_eta)
+    hist_plot(hist(new_seeds.seed_eff_v_eta), include_errors=errors, label='New')
+    hist_plot(hist(old_seeds.seed_eff_v_eta), include_errors=errors, title='Efficiency vs. Eta', label='Old')
+    plt.xlabel("Eta")
+    plt.ylim((0, 1.1))
+    plt.legend()
+
+    plt.sca(ax_phi)
+    hist_plot(hist(new_seeds.seed_eff_v_phi), include_errors=errors, label='New')
+    hist_plot(hist(old_seeds.seed_eff_v_phi), include_errors=errors, title='Efficiency vs. Phi', label='Old')
+    plt.xlabel("Phi")
+    plt.ylim((0, 1.1))
+
+@decl_plot
+def plot_track_eff(old_seeds, new_seeds):
+    r"""## Tracking Efficiency
+
+    The proportion of gen-level electrons origination in $\rho<1$cm and $|z|<10$cm from the beam spot that have
+    an associated GSF reconstructed electron.
+    """
+    _, ((ax_pt, ax_eta), (ax_phi, _)) = plt.subplots(2, 2)
+
+    errors = True
+    plt.sca(ax_pt)
+    hist_plot(hist(new_seeds.track_eff_v_pt), include_errors=errors, label='New')
+    hist_plot(hist(old_seeds.track_eff_v_pt), include_errors=errors, title='Efficiency vs. Pt', label='Old')
+    plt.xlabel("Pt")
+    plt.ylim((0, 1.1))
+
+    plt.sca(ax_eta)
+    hist_plot(hist(new_seeds.track_eff_v_eta), include_errors=errors, label='New')
+    hist_plot(hist(old_seeds.track_eff_v_eta), include_errors=errors, title='Efficiency vs. Eta', label='Old')
+    plt.xlabel("Eta")
+    plt.ylim((0, 1.1))
+    plt.legend()
+
+    plt.sca(ax_phi)
+    hist_plot(hist(new_seeds.track_eff_v_phi), include_errors=errors, label='New')
+    hist_plot(hist(old_seeds.track_eff_v_phi), include_errors=errors, title='Efficiency vs. Phi', label='Old')
+    plt.xlabel("Phi")
+    plt.ylim((0, 1.1))
+
+@decl_plot
+def plot_seed_purity(old_seeds, new_seeds):
+    r"""## Seed Purity
+
+    The proportion of ECAL-driven seeds that have a matched gen-level electron originating in
+    $\rho<1$cm and $|z|<10$cm from the beam spot.
+    """
+    _, ((ax_pt, ax_eta), (ax_phi, _)) = plt.subplots(2, 2)
+
+    errors = True
+    plt.sca(ax_pt)
+    hist_plot(hist(new_seeds.seed_pur_v_pt), include_errors=errors, label='New')
+    hist_plot(hist(old_seeds.seed_pur_v_pt), include_errors=errors, title='Purity vs. Pt', label='Old')
+    plt.xlabel("Pt")
+    plt.ylim((0, 1.1))
+
+    plt.sca(ax_eta)
+    hist_plot(hist(new_seeds.seed_pur_v_eta), include_errors=errors, label='New')
+    hist_plot(hist(old_seeds.seed_pur_v_eta), include_errors=errors, title='Purity vs. Eta', label='Old')
+    plt.xlabel("Eta")
+    plt.ylim((0, 1.1))
+    plt.legend()
+
+    plt.sca(ax_phi)
+    hist_plot(hist(new_seeds.seed_pur_v_phi), include_errors=errors, label='New')
+    hist_plot(hist(old_seeds.seed_pur_v_phi), include_errors=errors, title='Purity vs. Phi', label='Old')
+    plt.xlabel("Phi")
+    plt.ylim((0, 1.1))
+
+if __name__ == '__main__':
+    old_seeds = ResultSet("old_seeds", 'build/old_seeding.root')
+    new_seeds = ResultSet("new_seeds", 'build/new_seeding.root')
+
+    # Next, declare all of the (sub)plots that will be assembled into full
+    # figures later
+    seed_eff = (plot_seed_eff, (old_seeds, new_seeds), {})
+    track_eff = (plot_track_eff, (old_seeds, new_seeds), {})
+    seed_pur = (plot_seed_purity, (old_seeds, new_seeds), {})
+
+    # Now assemble the plots into figures.
+    plots = [
+        Plot([[seed_eff]],
+             'Seeding Efficiency'),
+        Plot([[track_eff]],
+             'Tracking Efficiency'),
+        Plot([[seed_pur]],
+             'ECAL-Driven Seed Purity'),
+    ]
+
+    # Finally, render and save the plots and generate the html+bootstrap
+    # dashboard to view them
+    render_plots(plots, to_disk=False)
+    generate_dashboard(plots, 'Seeding Efficiency', output='eff_plots.html', source_file=__file__)

+ 1 - 1
filval

@@ -1 +1 @@
-Subproject commit 7824f845e477f6b37413bed3bfcc2df78eb84d51
+Subproject commit 74e26df8bf18ef7f23a7110248136b8161f5acd0

analysis/plots.py → plots.py