Adds hoe info to ntuple

looper/analysis/TrackingNtupleMod.h

@@ -146,6 +146,7 @@ public :
    std::vector<std::vector<int> > *see_simTrkIdx;
    std::vector<unsigned int> *see_offset;
    std::vector<float> *see_Et;
+   std::vector<float> *see_hoe;
    std::vector<std::vector<int> > *see_isBarrel;
    std::vector<std::vector<int> > *see_layerOrDiskNr;
    std::vector<std::vector<int> > *see_isValidPos;
@@ -314,6 +315,7 @@ public :
    TBranch        *b_see_simTrkIdx;   //!
    TBranch        *b_see_offset;   //!
    TBranch        *b_see_Et;   //!
+   TBranch        *b_see_hoe;   //!
    TBranch        *b_see_isBarrel;   //!
    TBranch        *b_see_layerOrDiskNr;   //!
    TBranch        *b_see_isValidPos;   //!
@@ -534,6 +536,7 @@ void TrackingNtuple::Init(TTree *tree)
    see_simTrkIdx = 0;
    see_offset = 0;
    see_Et = 0;
+   see_hoe = 0;
    see_isBarrel = 0;
    see_layerOrDiskNr = 0;
    see_isValidPos = 0;
@@ -706,6 +709,7 @@ void TrackingNtuple::Init(TTree *tree)
    fChain->SetBranchAddress("see_simTrkIdx", &see_simTrkIdx, &b_see_simTrkIdx);
    fChain->SetBranchAddress("see_offset", &see_offset, &b_see_offset);
    fChain->SetBranchAddress("see_Et", &see_Et, &b_see_Et);
+   fChain->SetBranchAddress("see_hoe", &see_hoe, &b_see_hoe);
    fChain->SetBranchAddress("see_isBarrel", &see_isBarrel, &b_see_isBarrel);
    fChain->SetBranchAddress("see_layerOrDiskNr", &see_layerOrDiskNr, &b_see_layerOrDiskNr);
    fChain->SetBranchAddress("see_isValidPos", &see_isValidPos, &b_see_isValidPos);

looper/analysis/TrackingNtupleObjs.hpp

@@ -1,4 +1,4 @@
-/** TrackingNtupleObjs.hpp created on 2018-01-29 17:54:20.680281 by generate_class.py
+/** TrackingNtupleObjs.hpp created on 2018-03-08 10:59:46.474049 by generate_class.py
  * AVOID EDITING THIS FILE BY HAND!! Instead edit TrackingNtupleObjs.yaml and re-run
  * generate_class.py
  */
@@ -56,6 +56,7 @@ class SeedCollection {
     Value<vector<vector<int>>>* val_simTrkIdx;
     Value<vector<unsigned int>>* val_offset;
     Value<vector<float>>* val_Et;
+    Value<vector<float>>* val_hoe;
     Value<vector<vector<int>>>* val_isBarrel;
     Value<vector<vector<int>>>* val_layerOrDiskNr;
     Value<vector<vector<int>>>* val_isValidPos;
@@ -95,6 +96,7 @@ class SeedCollection {
         val_simTrkIdx = tds.track_branch_obj<vector<vector<int>>>("see_simTrkIdx");
         val_offset = tds.track_branch_obj<vector<unsigned int>>("see_offset");
         val_Et = tds.track_branch_obj<vector<float>>("see_Et");
+        val_hoe = tds.track_branch_obj<vector<float>>("see_hoe");
         val_isBarrel = tds.track_branch_obj<vector<vector<int>>>("see_isBarrel");
         val_layerOrDiskNr = tds.track_branch_obj<vector<vector<int>>>("see_layerOrDiskNr");
         val_isValidPos = tds.track_branch_obj<vector<vector<int>>>("see_isValidPos");
@@ -145,6 +147,7 @@ struct Seed {
     const vector<int>& simTrkIdx() const {return (*collection->val_simTrkIdx)().at(idx);}
     const unsigned int& offset() const {return (*collection->val_offset)().at(idx);}
     const float& Et() const {return (*collection->val_Et)().at(idx);}
+    const float& hoe() const {return (*collection->val_hoe)().at(idx);}
     const vector<int>& isBarrel() const {return (*collection->val_isBarrel)().at(idx);}
     const vector<int>& layerOrDiskNr() const {return (*collection->val_layerOrDiskNr)().at(idx);}
     const vector<int>& isValidPos() const {return (*collection->val_isValidPos)().at(idx);}

looper/analysis/TrackingNtupleObjs.yaml

@@ -53,6 +53,8 @@ Seed:
         type: unsigned int
       - name: Et
         type: float
+      - name: hoe
+        type: float
       - name: isBarrel
         type: vector<int>
       - name: layerOrDiskNr

looper/analysis/config.yaml

@@ -1,32 +1,25 @@
 #max-events: 50000
 debug: false
 
-output-file: ../hists/extra-wide-window.root
-source-file-key: extra-wide-window
+output-file: ../hists/old-seeding.root
+source-file-key: old-seeding
+
+hoe_cut: 0.15
 
 extra-narrow-window:
-    files: ../data/new_trees_extra_narrow_windows/trackingNtuple_*.root
+    files: ../data/EGamma_18-03-07_trackingNtuple_extra-narrow/180307_215057/0000/trackingNtuple_*.root
 
 narrow-window:
-    files: ../data/new_trees_narrow_windows/trackingNtuple_*.root
+    files: ../data/EGamma_18-03-07_trackingNtuple_narrow/180307_215909/0000/trackingNtuple_*.root
 
 wide-window:
-    files: ../data/new_trees_wide_windows/trackingNtuple_*.root
+    files: ../data/EGamma_18-03-07_trackingNtuple_wide/180307_220112/0000/trackingNtuple_*.root
 
 extra-wide-window:
-    files: ../data/new_trees_extra_wide_windows/trackingNtuple_*.root
-
-nwp-tight-window:
-    files: ../data/new_trees_nwp_tight_windows/trackingNtuple_*.root
-
-nwp-window:
-    files: ../data/new_trees_nwp_windows/trackingNtuple_*.root
+    files: ../data/EGamma_18-03-07_trackingNtuple_extra-wide/180307_220241/0000/trackingNtuple_*.root
 
 old-seeding:
-    files: ../data/new_trees_old_seeding/trackingNtuple_*.root
-
-nwp-eta-breakdown:
-    files: ../data/new_trees_nwp_eta_breakdown/trackingNtuple_*.root
+    files: ../data/EGamma_18-03-08_trackingNtuple_old-seeding/180308_171441/0000/trackingNtuple_*.root
 
 hist-params:
 
@@ -181,4 +174,10 @@ hist-params:
         label_x: "delta E"
         nbins_x: 200
         low_x: -1
-        low_x: 1
+        high_x: 1
+
+    n_seeds:
+        label_x: "Number of Seeds in Event"
+        nbins_x: 201
+        low_x: -0.5
+        high_x: 200.5

looper/analysis/tracking_eff.cpp

@@ -46,6 +46,10 @@ bool is_good_sim(const SimTrack& sim_track) {
     return abs(sim_track.pdgId()) == 11 and in_lum_region(vertex);
 };
 
+bool is_good_seed(const Seed& seed, float hoe_cut) {
+    return seed.isECALDriven() and seed.hoe() < hoe_cut;
+}
+
 float reco_energy_rel_err(const SimTrack& sim_track, const Seed& seed) {
     return (sim_track.pt() - seed.Et()) / sim_track.pt() ;
 }
@@ -54,6 +58,14 @@ bool reco_energy_consistent(const SimTrack& sim_track, const Seed& seed, float c
     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) {
+    return fabs(reco_energy_rel_err(sim_track, track)) < consistency_cut;
+}
+
 void run(bool silent){
     using namespace std;
     using namespace fv;
@@ -63,12 +75,8 @@ void run(bool silent){
     TrackingDataSet tds(output_filename, file_list, "trackingNtuple/tree");
     register_objects(tds);
 
-    int x = 12;
-    int y = 13;
 
-    for (auto z : vector<pair<int,int>>{{x, x}, {y, y}}) {
-        cout << z.first << endl;
-    }
+    float hoe_cut = the_config->get("hoe_cut").as<float>(999);
 
     bs = {
         tds.track_branch<float>("bsp_x"),
@@ -79,37 +87,70 @@ void run(bool silent){
         tds.track_branch<float>("bsp_sigmaz")
     };
 
-    std::map<std::tuple<int, int, bool>, ContainerTH2<float>*> BPIX_residuals_dRz;
-    std::map<std::tuple<int, int, bool>, ContainerTH2<float>*> BPIX_residuals_dPhi;
-    std::map<std::tuple<int, int, bool>, ContainerTH2<float>*> FPIX_residuals_dRz;
-    std::map<std::tuple<int, int, bool>, ContainerTH2<float>*> FPIX_residuals_dPhi;
+    enum TMType {
+        NoMatch = 0,
+        SeedMatched = 1,
+        TrackMatched = 2,
+        SeedAndTrackMatched = 3
+    };
+
+    std::map<std::tuple<int, int>, ContainerTH2<float>*> residuals_dRz;
+    std::map<std::tuple<int, int>, ContainerTH2<float>*> residuals_dPhi;
+
+    std::map<std::tuple<int, int, int>, ContainerTH2<float>*> BPIX_residuals_dRz;
+    std::map<std::tuple<int, int, int>, ContainerTH2<float>*> BPIX_residuals_dPhi;
+    std::map<std::tuple<int, int, int>, ContainerTH2<float>*> FPIX_residuals_dRz;
+    std::map<std::tuple<int, int, int>, ContainerTH2<float>*> FPIX_residuals_dPhi;
 
-    int layer, hit, isFake;
     stringstream name;
-    auto set_name = [&name, &layer, &hit, &isFake](const std::string& var, const std::string& region) {
+    auto set_name = [&name](const std::string& var, const std::string& region,
+                            const int& layer, const int& hit, const int& tm_type) {
         name.str("");
-        name << var << "_" << region << "_L" << layer << "_H" << hit << "_v_Et";
-        if (isFake)
-            name << "_fake";
-
+        name << var << "_" << region;
+        if (layer)
+            name << "_L" << layer;
+        name << "_H" << hit << "_v_Et";
+        switch (tm_type) {
+            case NoMatch:
+                name << "_NoMatch";
+                break;
+            case SeedMatched:
+                name << "_SeedMatched";
+                break;
+            case TrackMatched:
+                name << "_TrackMatched";
+                break;
+            case SeedAndTrackMatched:
+                name << "_SeedAndTrackMatched";
+                break;
+            default:
+                break;
+        }
     };
 
     THParams hist_params;
-    for (layer=1; layer<=4; layer++) {
-        for (hit=1; hit<=3; hit++) {
-            for (isFake=0; isFake<=1; isFake++) {
+    for (int hit=1; hit<=3; hit++) {
+        for (int tm_type = NoMatch; tm_type <= SeedAndTrackMatched; tm_type++) {
+            for (int layer=1; layer<=4; layer++) {
                 hist_params = (hit == 1) ? THParams::lookup("dRz_v_Et") : THParams::lookup("dRz_v_Et_outer_hits");
-                set_name("dRz", "BPIX");
-                BPIX_residuals_dRz[{layer, hit, isFake}] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
-                set_name("dRz", "FPIX");
-                FPIX_residuals_dRz[{layer, hit, isFake}] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
+                set_name("dRz", "BPIX", layer, hit, tm_type);
+                BPIX_residuals_dRz[{layer, hit, tm_type}] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
+                set_name("dRz", "FPIX", layer, hit, tm_type);
+                FPIX_residuals_dRz[{layer, hit, tm_type}] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
 
                 hist_params = (hit == 1) ? THParams::lookup("dPhi_v_Et") : THParams::lookup("dPhi_v_Et_outer_hits");
-                set_name("dPhi", "BPIX");
-                BPIX_residuals_dPhi[{layer, hit, isFake}] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
-                set_name("dPhi", "FPIX");
-                FPIX_residuals_dPhi[{layer, hit, isFake}] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
+                set_name("dPhi", "BPIX", layer, hit, tm_type);
+                BPIX_residuals_dPhi[{layer, hit, tm_type}] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
+                set_name("dPhi", "FPIX", layer, hit, tm_type);
+                FPIX_residuals_dPhi[{layer, hit, tm_type}] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
             }
+            hist_params = (hit == 1) ? THParams::lookup("dRz_v_Et") : THParams::lookup("dRz_v_Et_outer_hits");
+            set_name("dRz", "ALL", 0, hit, tm_type);
+            residuals_dRz[{hit, tm_type}] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
+
+            hist_params = (hit == 1) ? THParams::lookup("dPhi_v_Et") : THParams::lookup("dPhi_v_Et_outer_hits");
+            set_name("dPhi", "ALL", 0, hit, tm_type);
+            residuals_dPhi[{hit, tm_type}] = tds.register_container<ContainerTH2<float>>(name.str(), hist_params);
         }
     }
 
@@ -192,6 +233,8 @@ void run(bool silent){
 
     auto& ecal_energy_resolution = *tds.register_container<ContainerTH1<float>>("ecal_energy_resolution ", THParams::lookup("ecal_energy_resolution"));
 
+    auto& n_seeds = *tds.register_container<ContainerTH1<size_t>>("n_seeds", THParams::lookup("n_seeds"));
+
     while (tds.next(!silent)) {
 
         for (const auto& sim_track : sim_tracks) {
@@ -199,7 +242,7 @@ void run(bool silent){
             if (seeds.size() == 0) continue;
             for (const auto &seed_idx : sim_track.seedIdx()) {
                 const Seed& seed = seeds[seed_idx];
-                if (not seed.isECALDriven()) continue;
+                if (!is_good_seed(seed, hoe_cut)) continue;
                 ecal_energy_resolution.fill(reco_energy_rel_err(sim_track, seed));
             }
         }
@@ -211,7 +254,7 @@ void run(bool silent){
             bool matched = false;
             for (const auto& seed_idx : sim_track.seedIdx()) {
                 const Seed& seed = seeds[seed_idx];
-                if (seed.isECALDriven()) {
+                if (is_good_seed(seed, hoe_cut)) {
                     matched=true;
                     break;
                 }
@@ -232,7 +275,7 @@ void run(bool silent){
             bool matched = false;
             for (const auto& track_idx : sim_track.trkIdx()) {
                 const Seed& seed = seeds[gsf_tracks[track_idx].seedIdx()];
-                if (seed.isECALDriven()) {
+                if (is_good_seed(seed, hoe_cut)) {
                     matched=true;
                     break;
                 }
@@ -250,7 +293,7 @@ void run(bool silent){
             matched = false;
             for (const auto& seed_idx : sim_track.seedIdx()) {
                 const Seed& seed = seeds[seed_idx];
-                if (seed.isECALDriven() and seed.trkIdx()>=0 and reco_energy_consistent(sim_track, seed)) {
+                if (is_good_seed(seed, hoe_cut) and seed.trkIdx()>=0 and reco_energy_consistent(sim_track, seed)) {
                     matched=true;
                     break;
                 }
@@ -265,7 +308,7 @@ void run(bool silent){
 
         // Seeding Purity
         for (const auto& seed : seeds) {
-            if (!seed.isECALDriven()) continue;
+            if (!is_good_seed(seed, hoe_cut)) continue;
             bool match = false;
             for (const auto& sim_track_idx : seed.simTrkIdx()) {
                 if(is_good_sim(sim_tracks[sim_track_idx])) {
@@ -284,6 +327,8 @@ void run(bool silent){
         // Tracking Purity
         for (const auto& gsf_track : gsf_tracks) {
             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;
             for (const auto& sim_track_idx : gsf_track.simTrkIdx()) {
                 if (is_good_sim(sim_tracks[sim_track_idx])) {
@@ -299,7 +344,7 @@ void run(bool silent){
                 tracking_pur_v_phi.fill(gsf_track.phi(), match);
 
             match = false;
-            for (const auto& sim_track_idx : seeds[gsf_track.seedIdx()].simTrkIdx()) {
+            for (const auto& sim_track_idx : seed.simTrkIdx()) {
                 if (is_good_sim(sim_tracks[sim_track_idx])) {
                     match = true;
                     break;
@@ -316,18 +361,37 @@ void run(bool silent){
         // Hit Residuals
         for (const auto& seed : seeds) {
             if (seed.trkIdx() < 0) continue;  // require that seed produced gsf-track
-            if (!seed.isECALDriven()) continue;
-            bool is_sim_matched = false;
+            if (!is_good_seed(seed, hoe_cut)) continue;
+            bool is_seed_sim_matched = false;
             for (const auto& sim_track_idx : seed.simTrkIdx()) {
                 const SimTrack& sim_track = sim_tracks[sim_track_idx];
                 if (is_good_sim(sim_track) and reco_energy_consistent(sim_track, seed)) {
-                    is_sim_matched = true;
+                    is_seed_sim_matched = true;
                     break;
-//                    std::cout << "Matched to non-electron: " << sim_tracks[sim_track_idx].pdgId()
-//                              << " (" << seed.simTrkIdx().size() << ")" << std::endl;
                 }
             }
+
+            bool is_track_sim_matched = false;
             const auto the_track = gsf_tracks[seed.trkIdx()];
+            for (const auto& sim_track_idx : the_track.simTrkIdx()) {
+                const SimTrack& sim_track = sim_tracks[sim_track_idx];
+                if (is_good_sim(sim_track) and reco_energy_consistent(sim_track, the_track)) {
+                    is_track_sim_matched = true;
+                    break;
+                }
+            }
+
+            std::vector<TMType> tm_types;
+            if (is_seed_sim_matched and is_track_sim_matched) {
+                tm_types = {SeedAndTrackMatched, SeedMatched, TrackMatched};
+            } else if (is_seed_sim_matched) {
+                tm_types = {SeedMatched};
+            } else if (is_track_sim_matched) {
+                tm_types = {TrackMatched};
+            } else {
+                tm_types = {NoMatch};
+            }
+
             vector<int>   hits_valid;
             vector<float> hits_dRz;
             vector<float> hits_dPhi;
@@ -343,41 +407,45 @@ void run(bool silent){
             const vector<int>& hits_isBarrel = seed.isBarrel();
             const vector<int>& hits_layerOrDiskNr = seed.layerOrDiskNr();
 
-            size_t nHits = hits_valid.size();
-            for (size_t hit_idx=0; hit_idx<nHits; hit_idx++) {
-                if (!hits_valid[hit_idx]) continue;
-                bool isBarrel = hits_isBarrel[hit_idx] == 1;
-                int layerOrDiskNr = hits_layerOrDiskNr[hit_idx];
-                float dRz = abs(hits_dRz[hit_idx]);
-                float dPhi = abs(hits_dPhi[hit_idx]);
-                int eta_region = get_region(static_cast<int>(hit_idx + 1), seed.eta());
-
-
-                if (is_sim_matched) {
-                    dRz_residuals_v_eta[hit_idx + 1]->fill(dRz, abs(seed.eta()));
-                    dPhi_residuals_v_eta[hit_idx + 1]->fill(dPhi, abs(seed.eta()));
-
-                    dRz_residuals_v_region[{hit_idx + 1, eta_region}]->fill(dRz, seed.Et());
-                    dPhi_residuals_v_region[{hit_idx + 1, eta_region}]->fill(dPhi, seed.Et());
-
-                }
-
-                if (layerOrDiskNr == -1) continue;  // subDet is not set w/ old seeding
-
-                if (isBarrel)
-                    hit_vs_layer_barrel.fill(layerOrDiskNr, static_cast<const int &>(hit_idx + 1));
-                else
-                    hit_vs_layer_forward.fill(layerOrDiskNr, static_cast<const int &>(hit_idx + 1));
-
-                if (isBarrel) {
-                    BPIX_residuals_dRz[{layerOrDiskNr, hit_idx+1, !is_sim_matched}]->fill(dRz, seed.Et());
-                    BPIX_residuals_dPhi[{layerOrDiskNr, hit_idx+1, !is_sim_matched}]->fill(dPhi, seed.Et());
-                } else {
-                    FPIX_residuals_dRz[{layerOrDiskNr, hit_idx+1, !is_sim_matched}]->fill(dRz, seed.Et());
-                    FPIX_residuals_dPhi[{layerOrDiskNr, hit_idx+1, !is_sim_matched}]->fill(dPhi, seed.Et());
+            for (const auto& tm_type : tm_types) {
+                size_t nHits = hits_valid.size();
+                for (size_t hit_idx = 0; hit_idx < nHits; hit_idx++) {
+                    if (!hits_valid[hit_idx]) continue;
+                    bool isBarrel = hits_isBarrel[hit_idx] == 1;
+                    int layerOrDiskNr = hits_layerOrDiskNr[hit_idx];
+                    float dRz = abs(hits_dRz[hit_idx]);
+                    float dPhi = abs(hits_dPhi[hit_idx]);
+                    int eta_region = get_region(static_cast<int>(hit_idx + 1), seed.eta());
+
+
+                    if (is_seed_sim_matched) {
+                        dRz_residuals_v_eta[hit_idx + 1]->fill(dRz, abs(seed.eta()));
+                        dPhi_residuals_v_eta[hit_idx + 1]->fill(dPhi, abs(seed.eta()));
+
+                        dRz_residuals_v_region[{hit_idx + 1, eta_region}]->fill(dRz, seed.Et());
+                        dPhi_residuals_v_region[{hit_idx + 1, eta_region}]->fill(dPhi, seed.Et());
+                    }
+
+                    residuals_dRz[{hit_idx + 1, tm_type}]->fill(dRz, seed.Et());
+                    residuals_dPhi[{hit_idx + 1, tm_type}]->fill(dPhi, seed.Et());
+                    if (layerOrDiskNr == -1) continue;  // subDet is not set w/ old seeding
+
+                    if (isBarrel)
+                        hit_vs_layer_barrel.fill(layerOrDiskNr, static_cast<const int &>(hit_idx + 1));
+                    else
+                        hit_vs_layer_forward.fill(layerOrDiskNr, static_cast<const int &>(hit_idx + 1));
+
+                    if (isBarrel) {
+                        BPIX_residuals_dRz[{layerOrDiskNr, hit_idx + 1, tm_type}]->fill(dRz, seed.Et());
+                        BPIX_residuals_dPhi[{layerOrDiskNr, hit_idx + 1, tm_type}]->fill(dPhi, seed.Et());
+                    } else {
+                        FPIX_residuals_dRz[{layerOrDiskNr, hit_idx + 1, tm_type}]->fill(dRz, seed.Et());
+                        FPIX_residuals_dPhi[{layerOrDiskNr, hit_idx + 1, tm_type}]->fill(dPhi, seed.Et());
+                    }
                 }
             }
         }
+        n_seeds.fill(seeds.size());
     }
 
     tds.save_all();