Adds info on prompt reconstruction

looper/analysis/TrackingNtupleObjs.hpp

@@ -1,4 +1,4 @@
-/** TrackingNtupleObjs.hpp created on 2018-05-22 17:07:59.784846 by generate_class.py
+/** TrackingNtupleObjs.hpp created on 2018-06-06 16:43:59.886647 by generate_class.py
  * AVOID EDITING THIS FILE BY HAND!! Instead edit TrackingNtupleObjs.yaml and re-run
  * generate_class.py
  */
@@ -430,6 +430,10 @@ struct Seed {
 
 };
 
+bool operator==(const Seed& obj1, const Seed& obj2) {
+    return obj1.idx == obj2.idx;
+}
+
 Seed SeedCollection::iter::operator*() {
     return {collection, idx};
 }
@@ -992,12 +996,202 @@ struct Track {
 
 };
 
+bool operator==(const Track& obj1, const Track& obj2) {
+    return obj1.idx == obj2.idx;
+}
+
 Track TrackCollection::iter::operator*() {
     return {collection, idx};
 }
 Track TrackCollection::operator[](size_t idx) {
     return {this, idx};
 }
+struct Gen;
+
+class GenCollection {
+  public:
+    class iter {
+      public:
+        iter(GenCollection* collection, size_t idx)
+          :collection(collection), idx(idx) { }
+        iter operator++() { ++idx; return *this; }
+        bool operator!=(const iter & other) { return idx != other.idx; }
+        Gen operator*();
+      private:
+        GenCollection* collection;
+        size_t idx;
+    };
+    TrackingDataSet* tds;
+
+    Value<vector<float>>* val_px;
+    bool px_loaded;
+    Value<vector<float>>* val_py;
+    bool py_loaded;
+    Value<vector<float>>* val_pz;
+    bool pz_loaded;
+    Value<vector<float>>* val_charge;
+    bool charge_loaded;
+    Value<vector<int>>* val_pdgId;
+    bool pdgId_loaded;
+    Value<vector<float>>* val_vx;
+    bool vx_loaded;
+    Value<vector<float>>* val_vy;
+    bool vy_loaded;
+    Value<vector<float>>* val_vz;
+    bool vz_loaded;
+    Value<vector<int>>* val_status;
+    bool status_loaded;
+    Value<vector<int>>* val_mother;
+    bool mother_loaded;
+    Value<vector<bool>>* val_isTauDecayProduct;
+    bool isTauDecayProduct_loaded;
+    Value<vector<bool>>* val_isDirectHadronDecayProduct;
+    bool isDirectHadronDecayProduct_loaded;
+    Value<vector<bool>>* val_isPrompt;
+    bool isPrompt_loaded;
+
+    GenCollection() { }
+
+    void init(TrackingDataSet* tds){
+        this->tds = tds;
+    }
+
+    size_t size() {
+        if (!this->px_loaded) {
+            this->val_px = this->tds->track_branch_obj<vector<float>>("gen_px");
+            this->px_loaded = true;
+        }
+        return (*this->val_px)().size();
+    }
+
+
+    Gen operator[](size_t);
+    iter begin() { return iter(this, 0); }
+    iter end() { return iter(this, size()); }
+};
+
+struct Gen {
+    GenCollection* collection;
+    size_t idx;
+    Gen(GenCollection* collection, const size_t idx)
+      :collection(collection), idx(idx) { }
+
+    const float& px() const {
+        if (!collection->px_loaded) {
+            collection->val_px = collection->tds->track_branch_obj<vector<float>>("gen_px");
+            collection->px_loaded = true;
+        }
+        return (*collection->val_px)().at(idx);
+    }
+
+    const float& py() const {
+        if (!collection->py_loaded) {
+            collection->val_py = collection->tds->track_branch_obj<vector<float>>("gen_py");
+            collection->py_loaded = true;
+        }
+        return (*collection->val_py)().at(idx);
+    }
+
+    const float& pz() const {
+        if (!collection->pz_loaded) {
+            collection->val_pz = collection->tds->track_branch_obj<vector<float>>("gen_pz");
+            collection->pz_loaded = true;
+        }
+        return (*collection->val_pz)().at(idx);
+    }
+
+    const float& charge() const {
+        if (!collection->charge_loaded) {
+            collection->val_charge = collection->tds->track_branch_obj<vector<float>>("gen_charge");
+            collection->charge_loaded = true;
+        }
+        return (*collection->val_charge)().at(idx);
+    }
+
+    const int& pdgId() const {
+        if (!collection->pdgId_loaded) {
+            collection->val_pdgId = collection->tds->track_branch_obj<vector<int>>("gen_pdgId");
+            collection->pdgId_loaded = true;
+        }
+        return (*collection->val_pdgId)().at(idx);
+    }
+
+    const float& vx() const {
+        if (!collection->vx_loaded) {
+            collection->val_vx = collection->tds->track_branch_obj<vector<float>>("gen_vx");
+            collection->vx_loaded = true;
+        }
+        return (*collection->val_vx)().at(idx);
+    }
+
+    const float& vy() const {
+        if (!collection->vy_loaded) {
+            collection->val_vy = collection->tds->track_branch_obj<vector<float>>("gen_vy");
+            collection->vy_loaded = true;
+        }
+        return (*collection->val_vy)().at(idx);
+    }
+
+    const float& vz() const {
+        if (!collection->vz_loaded) {
+            collection->val_vz = collection->tds->track_branch_obj<vector<float>>("gen_vz");
+            collection->vz_loaded = true;
+        }
+        return (*collection->val_vz)().at(idx);
+    }
+
+    const int& status() const {
+        if (!collection->status_loaded) {
+            collection->val_status = collection->tds->track_branch_obj<vector<int>>("gen_status");
+            collection->status_loaded = true;
+        }
+        return (*collection->val_status)().at(idx);
+    }
+
+    const int& mother() const {
+        if (!collection->mother_loaded) {
+            collection->val_mother = collection->tds->track_branch_obj<vector<int>>("gen_mother");
+            collection->mother_loaded = true;
+        }
+        return (*collection->val_mother)().at(idx);
+    }
+
+    const bool& isTauDecayProduct() const {
+        if (!collection->isTauDecayProduct_loaded) {
+            collection->val_isTauDecayProduct = collection->tds->track_branch_obj<vector<bool>>("gen_isTauDecayProduct");
+            collection->isTauDecayProduct_loaded = true;
+        }
+        return (*collection->val_isTauDecayProduct)().at(idx);
+    }
+
+    const bool& isDirectHadronDecayProduct() const {
+        if (!collection->isDirectHadronDecayProduct_loaded) {
+            collection->val_isDirectHadronDecayProduct = collection->tds->track_branch_obj<vector<bool>>("gen_isDirectHadronDecayProduct");
+            collection->isDirectHadronDecayProduct_loaded = true;
+        }
+        return (*collection->val_isDirectHadronDecayProduct)().at(idx);
+    }
+
+    const bool& isPrompt() const {
+        if (!collection->isPrompt_loaded) {
+            collection->val_isPrompt = collection->tds->track_branch_obj<vector<bool>>("gen_isPrompt");
+            collection->isPrompt_loaded = true;
+        }
+        return (*collection->val_isPrompt)().at(idx);
+    }
+
+};
+
+bool operator==(const Gen& obj1, const Gen& obj2) {
+    return obj1.idx == obj2.idx;
+}
+
+Gen GenCollection::iter::operator*() {
+    return {collection, idx};
+}
+Gen GenCollection::operator[](size_t idx) {
+    return {this, idx};
+}
 struct SimTrack;
 
 class SimTrackCollection {
@@ -1324,6 +1518,10 @@ struct SimTrack {
 
 };
 
+bool operator==(const SimTrack& obj1, const SimTrack& obj2) {
+    return obj1.idx == obj2.idx;
+}
+
 SimTrack SimTrackCollection::iter::operator*() {
     return {collection, idx};
 }
@@ -1456,6 +1654,10 @@ struct SimVertex {
 
 };
 
+bool operator==(const SimVertex& obj1, const SimVertex& obj2) {
+    return obj1.idx == obj2.idx;
+}
+
 SimVertex SimVertexCollection::iter::operator*() {
     return {collection, idx};
 }
@@ -1588,6 +1790,10 @@ struct SuperCluster {
 
 };
 
+bool operator==(const SuperCluster& obj1, const SuperCluster& obj2) {
+    return obj1.idx == obj2.idx;
+}
+
 SuperCluster SuperClusterCollection::iter::operator*() {
     return {collection, idx};
 }

looper/analysis/TrackingNtupleObjs.yaml

@@ -182,6 +182,38 @@ Track:
       - name: genDR
         type: float
 
+
+Gen:
+    treename_prefix: gen
+    fields:
+      - name: px
+        type: float
+      - name: py
+        type: float
+      - name: pz
+        type: float
+      - name: charge
+        type: float
+      - name: pdgId
+        type: int
+      - name: vx
+        type: float
+      - name: vy
+        type: float
+      - name: vz
+        type: float
+      - name: status
+        type: int
+      - name: mother
+        type: int
+      - name: isTauDecayProduct
+        type: bool
+      - name: isDirectHadronDecayProduct
+        type: bool
+      - name: isPrompt
+        type: bool
+
+
 SimTrack:
     treename_prefix: sim
     fields:

looper/analysis/config.yaml

@@ -1,4 +1,4 @@
-#max-events: 50000
+#max-events: 20
 debug: false
 
 #output-file: ../hists/ttbar-old-seeding.root
@@ -26,7 +26,7 @@ tt-old-default:
 
 
 local:
-    files: /home/caleb/Sources/EGamma/data/*.root
+    files: /home/caleb/Sources/EGamma/data/zee.root
 
 hist-params:
 

looper/analysis/tracking_eff.cpp

@@ -14,6 +14,7 @@
 typedef ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<float>> Vec4;
 
 SeedCollection seeds;
+GenCollection gens;
 SimTrackCollection sim_tracks;
 SimVertexCollection sim_vertices;
 TrackCollection gsf_tracks;
@@ -23,6 +24,7 @@ std::vector<SimTrack> sim_els;
 
 void register_objects(TrackingDataSet& tds){
     seeds.init(&tds);
+    gens.init(&tds);
     sim_tracks.init(&tds);
     sim_vertices.init(&tds);
     gsf_tracks.init(&tds);
@@ -85,23 +87,50 @@ Vec4 scl_p4(const SuperCluster& scl) {
     return Vec4(scl.px(), scl.py(), scl.pz(), scl.e());
 }
 
-double deltaR(float scl_eta, float scl_phi, const SimTrack& sim_track) {
-    return sqrt(pow(sim_track.eta() - scl_eta, 2.0) +
-                pow(sim_track.phi() - scl_phi, 2.0));
-}
-
-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));
+double deltaR(double eta1, double phi1, double eta2, double phi2) {
+    return sqrt(pow(eta1 - eta2, 2.0) +
+                pow(phi1 - phi2, 2.0));
 }
 
 void update_sim_els() {
     sim_els.clear();
     for (const SimTrack sim_track: sim_tracks) {
-        if (is_good_sim(sim_track)) sim_els.push_back(sim_track);
+        if (is_good_sim(sim_track)) {
+            sim_els.push_back(sim_track);
+        }
     }
 }
 
+std::vector<SimTrack> get_prompt_sims() {
+    /*
+     * Returns sim tracks that pass is_good_sim as well as match well with a
+     * prompt gen electron
+     */
+    std::vector<Gen> prompt_gen;
+    for (const auto& gen : gens) {
+        if (abs(gen.pdgId()) == 11 and gen.isPrompt()) {
+            prompt_gen.push_back(gen);
+        }
+    }
+
+    std::vector<SimTrack> prompt_sim;
+    for (const auto& gen : prompt_gen) {
+        float gen_phi = atan2(gen.py(), gen.px());
+        float gen_eta = pseudorapidityP(gen.px(), gen.py(), gen.pz());
+        int gen_id = gen.pdgId();
+        for (const auto& sim_el : sim_els) {
+            float sim_phi = sim_el.phi();
+            float sim_eta = sim_el.eta();
+            int sim_id = sim_el.pdgId();
+            if (gen_id == sim_id and deltaR(gen_eta, gen_phi, sim_eta, sim_phi) < 0.05) {
+                prompt_sim.push_back(sim_el);
+                break;
+            }
+        }
+    }
+    return prompt_sim;
+}
+
 
 auto get_match_stats(bool dRMatch) {
     int nSimTrack = 0;
@@ -121,30 +150,34 @@ auto get_match_stats(bool dRMatch) {
     for(const auto& gsf_track : gsf_tracks)
         gsf_matches[gsf_track.idx] = {};
 
-    nSimTrack = sim_matches.size();
-    nGSFTrack = gsf_matches.size();
+    nSimTrack = static_cast<int>(sim_matches.size());
+    nGSFTrack = static_cast<int>(gsf_matches.size());
 
     for (const auto& sim_track : sim_els) {
         if (dRMatch) {
             for (const auto& gsf_track : gsf_tracks) {
-                float dR = deltaR(sim_track, gsf_track);
+                double dR = deltaR(sim_track.eta(), sim_track.phi(), gsf_track.eta(), gsf_track.phi());
                 if (dR < 0.2) {
                     sim_matches[sim_track.idx].push_back(gsf_track.idx);
                     gsf_matches[gsf_track.idx].push_back(sim_track.idx);
                     matched_dR.push_back(dR);
                     matched_dpT.push_back((sim_track.pt() - gsf_track.pt())/sim_track.pt());
-                    if (gsf_track.q() != sim_track.q()) nFlipped++;
+                    if (gsf_track.q() != sim_track.q()) {
+                        nFlipped++;
+                    }
                 }
             }
         } else {
             for (const auto& gsfIdx : sim_track.trkIdx()) {
                 const Track& gsf_track = gsf_tracks[gsfIdx];
-                float dR = deltaR(sim_track, gsf_track);
+                double dR = deltaR(sim_track.eta(), sim_track.phi(), gsf_track.eta(), gsf_track.phi());
                 sim_matches[sim_track.idx].push_back(gsf_track.idx);
                 gsf_matches[gsf_track.idx].push_back(sim_track.idx);
-                matched_dR.push_back(dR);
+                matched_dR.push_back(static_cast<float &&>(dR));
                 matched_dpT.push_back((sim_track.pt() - gsf_track.pt())/sim_track.pt());
-                if (gsf_track.q() != sim_track.q()) nFlipped++;
+                if (gsf_track.q() != sim_track.q()) {
+                    nFlipped++;
+                }
             }
         }
     }
@@ -204,7 +237,7 @@ void run(){
                             const int& layer, const int& hit, const int& tm_type) {
         name.str("");
         name << var << "_" << region;
-        if (layer)
+        if (layer>0)
             name << "_L" << layer;
         name << "_H" << hit << "_v_Et";
         switch (tm_type) {
@@ -314,6 +347,11 @@ void run(){
             tds.register_container<ContainerTH1<float>>("scl_v_eta", THParams::lookup("eta")),
             tds.register_container<ContainerTH1<float>>("scl_v_phi", THParams::lookup("phi"))};
 
+    KinColl<ContainerTH1<float>> prompt_kinems = {
+            tds.register_container<ContainerTH1<float>>("prompt_v_pt",  THParams::lookup("pt")),
+            tds.register_container<ContainerTH1<float>>("prompt_v_eta", THParams::lookup("eta")),
+            tds.register_container<ContainerTH1<float>>("prompt_v_phi", THParams::lookup("phi"))};
+
 
     KinColl<EfficiencyContainer<float>> seed_eff = {
             tds.register_container<EfficiencyContainer<float>>("seed_eff_v_pt",  THParams::lookup("pt")),
@@ -323,6 +361,7 @@ void run(){
             tds.register_container<EfficiencyContainer<float>>("seed_pur_v_pt",  THParams::lookup("pt")),
             tds.register_container<EfficiencyContainer<float>>("seed_pur_v_eta", THParams::lookup("eta")),
             tds.register_container<EfficiencyContainer<float>>("seed_pur_v_phi", THParams::lookup("phi"))};
+
     KinColl<EfficiencyContainer<float>> tracking_eff = {
             tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_pt",  THParams::lookup("pt")),
             tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_eta", THParams::lookup("eta")),
@@ -331,6 +370,14 @@ void run(){
             tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_pt",  THParams::lookup("pt")),
             tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_eta", THParams::lookup("eta")),
             tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_phi", THParams::lookup("phi"))};
+    KinColl<EfficiencyContainer<float>> prompt_eff = {
+            tds.register_container<EfficiencyContainer<float>>("prompt_eff_v_pt",  THParams::lookup("pt")),
+            tds.register_container<EfficiencyContainer<float>>("prompt_eff_v_eta", THParams::lookup("eta")),
+            tds.register_container<EfficiencyContainer<float>>("prompt_eff_v_phi", THParams::lookup("phi"))};
+    KinColl<EfficiencyContainer<float>> prompt_pur = {
+            tds.register_container<EfficiencyContainer<float>>("prompt_pur_v_pt",  THParams::lookup("pt")),
+            tds.register_container<EfficiencyContainer<float>>("prompt_pur_v_eta", THParams::lookup("eta")),
+            tds.register_container<EfficiencyContainer<float>>("prompt_pur_v_phi", THParams::lookup("phi"))};
 
     KinColl<EfficiencyContainer<float>> tracking_eff_dR = {
             tds.register_container<EfficiencyContainer<float>>("tracking_eff_v_pt_dR",  THParams::lookup("pt")),
@@ -340,6 +387,14 @@ void run(){
             tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_pt_dR",  THParams::lookup("pt")),
             tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_eta_dR", THParams::lookup("eta")),
             tds.register_container<EfficiencyContainer<float>>("tracking_pur_v_phi_dR", THParams::lookup("phi"))};
+    KinColl<EfficiencyContainer<float>> prompt_eff_dR = {
+            tds.register_container<EfficiencyContainer<float>>("prompt_eff_v_pt_dR",  THParams::lookup("pt")),
+            tds.register_container<EfficiencyContainer<float>>("prompt_eff_v_eta_dR", THParams::lookup("eta")),
+            tds.register_container<EfficiencyContainer<float>>("prompt_eff_v_phi_dR", THParams::lookup("phi"))};
+    KinColl<EfficiencyContainer<float>> prompt_pur_dR = {
+            tds.register_container<EfficiencyContainer<float>>("prompt_pur_v_pt_dR",  THParams::lookup("pt")),
+            tds.register_container<EfficiencyContainer<float>>("prompt_pur_v_eta_dR", THParams::lookup("eta")),
+            tds.register_container<EfficiencyContainer<float>>("prompt_pur_v_phi_dR", THParams::lookup("phi"))};
 
     KinColl<EfficiencyContainer<float>> fake_rate = {
             tds.register_container<EfficiencyContainer<float>>("fake_rate_v_pt",  THParams::lookup("pt")),
@@ -392,6 +447,7 @@ void run(){
     auto& n_seeds = *tds.register_container<ContainerTH1<size_t>>("n_seeds", THParams::lookup("n_seeds"));
     auto& n_good_seeds = *tds.register_container<ContainerTH1<size_t>>("n_good_seeds", THParams::lookup("n_seeds"));
     auto& n_good_sim = *tds.register_container<ContainerTH1<int>>("n_good_sim", THParams::lookup("n_seeds"));
+    auto& n_prompt = *tds.register_container<ContainerTH1<int>>("n_prompt", THParams::lookup("n_seeds"));
     auto& n_gsf_tracks = *tds.register_container<ContainerTH1<int>>("n_gsf_track", THParams::lookup("n_seeds"));
     auto& n_scl = *tds.register_container<ContainerTH1<size_t>>("n_scl", THParams::lookup("n_seeds"));
     auto& n_good_scl = *tds.register_container<ContainerTH1<size_t>>("n_good_scl", THParams::lookup("n_seeds"));
@@ -430,6 +486,8 @@ void run(){
     while (tds.next()) {
 
         update_sim_els();
+        vector<SimTrack> prompt_sims = get_prompt_sims();
+
 
         size_t _n_good_seeds = 0;
         for (const auto& seed : seeds) {
@@ -449,6 +507,13 @@ void run(){
             good_sim_kinems.phi->fill(sim_track.phi());
         }
 
+        n_prompt.fill(prompt_sims.size());
+        for (const auto& prompt : prompt_sims) {
+            prompt_kinems.pt->fill(prompt.pt());
+            prompt_kinems.eta->fill(prompt.eta());
+            prompt_kinems.phi->fill(prompt.phi());
+        }
+
         for (const auto& gsf_track : gsf_tracks) {
             gsf_kinems.pt->fill(gsf_track.pt());
             gsf_kinems.eta->fill(gsf_track.eta());
@@ -534,6 +599,7 @@ void run(){
         for (const auto& sim_track : sim_els) {
             if (gsf_tracks.size() == 0) continue;
             bool matched = false;
+            bool prompt_matched = false;
             for (const auto& track_idx : sim_track.trkIdx()) {
                 const Seed& seed = seeds[gsf_tracks[track_idx].seedIdx()];
                 if (is_good_seed(seed, hoe_cut)) {
@@ -548,11 +614,20 @@ void run(){
             if (abs(sim_track.eta()) < 2.4 and sim_track.pt() > 20)
                 tracking_eff.phi->fill(sim_track.phi(), matched);
 
+            if (std::find(prompt_sims.begin(), prompt_sims.end(), sim_track) != prompt_sims.end()) {
+                if (abs(sim_track.eta()) < 2.4)
+                    prompt_eff.pt->fill(sim_track.pt(), matched);
+                if (sim_track.pt() > 20.0)
+                    prompt_eff.eta->fill(sim_track.eta(), matched);
+                if (abs(sim_track.eta()) < 2.4 and sim_track.pt() > 20)
+                    prompt_eff.phi->fill(sim_track.phi(), matched);
+            }
+
             // dR-matching
             matched = false;
             for (const auto& gsf_track : gsf_tracks) {
                 const Seed& seed = seeds[gsf_track.seedIdx()];
-                double dR = deltaR(sim_track, gsf_track);
+                double dR = deltaR(sim_track.eta(), sim_track.phi(), gsf_track.eta(), gsf_track.phi());
                 if (is_good_seed(seed, hoe_cut) and dR < 0.2) {
                     matched = true;
                     break;
@@ -564,6 +639,15 @@ void run(){
                 tracking_eff_dR.eta->fill(sim_track.eta(), matched);
             if (abs(sim_track.eta()) < 2.4 and sim_track.pt() > 20)
                 tracking_eff_dR.phi->fill(sim_track.phi(), matched);
+
+            if (find(prompt_sims.begin(), prompt_sims.end(), sim_track) != prompt_sims.end()) {
+                if (abs(sim_track.eta()) < 2.4)
+                    prompt_eff_dR.pt->fill(sim_track.pt(), matched);
+                if (sim_track.pt() > 20.0)
+                    prompt_eff_dR.eta->fill(sim_track.eta(), matched);
+                if (abs(sim_track.eta()) < 2.4 and sim_track.pt() > 20)
+                    prompt_eff_dR.phi->fill(sim_track.phi(), matched);
+            }
         }
 
         // Seeding Purity
@@ -590,10 +674,13 @@ void run(){
             const Seed& seed = seeds[gsf_track.seedIdx()];
             if (!is_good_seed(seed, hoe_cut)) continue;
             bool matched = false;
+            bool prompt_matched = false;
             for (const auto& sim_track_idx : gsf_track.simTrkIdx()) {
-                if (is_good_sim(sim_tracks[sim_track_idx])) {
-                    matched = true;
-                    break;
+                const auto& sim_track = sim_tracks[sim_track_idx];
+                if (!is_good_sim(sim_track)) continue;
+                matched = true;
+                if (find(prompt_sims.begin(), prompt_sims.end(), sim_track) != prompt_sims.end()) {
+                    prompt_matched = true;
                 }
             }
             if (abs(gsf_track.eta()) < 2.4)
@@ -603,13 +690,23 @@ void run(){
             if (abs(gsf_track.eta()) < 2.4 and gsf_track.pt() > 20)
                 tracking_pur.phi->fill(gsf_track.phi(), matched);
 
+            if (abs(gsf_track.eta()) < 2.4)
+                prompt_pur.pt->fill(gsf_track.pt(), prompt_matched);
+            if (gsf_track.pt() > 20)
+                prompt_pur.eta->fill(gsf_track.eta(), prompt_matched);
+            if (abs(gsf_track.eta()) < 2.4 and gsf_track.pt() > 20)
+                prompt_pur.phi->fill(gsf_track.phi(), prompt_matched);
+
             // dR-matching
             matched = false;
+            prompt_matched = false;
             for (const auto& sim_track : sim_els) {
-                double dR = deltaR(sim_track, gsf_track);
+                double dR = deltaR(sim_track.eta(), sim_track.phi(), gsf_track.eta(), gsf_track.phi());
                 if (dR < 0.2) {
                     matched = true;
-                    break;
+                    if (find(prompt_sims.begin(), prompt_sims.end(), sim_track) != prompt_sims.end()) {
+                        prompt_matched = true;
+                    }
                 }
             }
             if (abs(gsf_track.eta()) < 2.4)
@@ -618,6 +715,13 @@ void run(){
                 tracking_pur_dR.eta->fill(gsf_track.eta(), matched);
             if (abs(gsf_track.eta()) < 2.4 and gsf_track.pt() > 20)
                 tracking_pur_dR.phi->fill(gsf_track.phi(), matched);
+
+            if (abs(gsf_track.eta()) < 2.4)
+                prompt_pur_dR.pt->fill(gsf_track.pt(), prompt_matched);
+            if (gsf_track.pt() > 20)
+                prompt_pur_dR.eta->fill(gsf_track.eta(), prompt_matched);
+            if (abs(gsf_track.eta()) < 2.4 and gsf_track.pt() > 20)
+                prompt_pur_dR.phi->fill(gsf_track.phi(), prompt_matched);
         }
 
         // Fake-Rate
@@ -637,7 +741,7 @@ void run(){
         for (const auto &scl : scls) {
             Vec4 p4 = scl_p4(scl);
             for(const auto& sim_track : sim_els) {
-                if (deltaR(p4.eta(), p4.phi(), sim_track) > 0.2) continue;
+                if (deltaR(p4.eta(), p4.phi(), sim_track.eta(), sim_track.phi()) > 0.2) continue;
                 if (((p4.Et() - sim_track.pt()) / p4.Et()) > 0.1) continue;
                 tm_scls.insert(scl.idx);
                 break;

looper/filval

@@ -1 +1 @@
-Subproject commit 1316fbe1dba82784d3ba92f2005f2abf5c1dc054
+Subproject commit de09db9ee970db8e263c285e0e2d4128b545e9a5