diff --git a/.github/validation_samples/it_pileup/config.py b/.github/validation_samples/it_pileup/config.py index 07c840e0e0..d85645d242 100644 --- a/.github/validation_samples/it_pileup/config.py +++ b/.github/validation_samples/it_pileup/config.py @@ -104,6 +104,43 @@ hcal_veto = hcal.HcalVetoProcessor() hcal_veto.input_hit_pass_name = thisPassName +# Load and configure particle flow sequence. +# Here we use PF "tracking" and CLUE Ecal clustering +from LDMX.Recon import pfReco +trackPF = pfReco.pfTrackProducer() +trackPF.inputTrackCollName=trackPF.inputTrackCollName+overlayStr #"EcalScoringPlaneHitsOverlay" # +trackPF.input_pass_name=inputPassName +trackPF.doElectronTracking=True +# reference info +truthPF = pfReco.pfTruthProducer() + +# CLUE +import LDMX.Ecal.ecalClusters as cl +cluster = cl.EcalClusterProducer() +cluster.seedThreshold = 350. cluster.dc = 0.3 +cluster.nbrOfLayers = 1 +cluster.nbr_of_layers = cluster.nbrOfLayers #just here to highlight confusion +cluster.reclustering = True +cluster.rec_hit_pass_name=thisPassName #run on process+pileup + +# particle flow: +pfComb=pfReco.pfProducer() +pfComb.inputEcalCollName = cluster.cluster_coll_name # use CLUE +pfComb.input_ecal_passname = thisPassName +# trigger recasting existing CLUE to caloclusters +pfComb.use_existing_ecal_clusters = True + +from LDMX.Recon import pileupFinder +puFinder = pileupFinder.pileupFinder() +puFinder.rec_hit_pass_name=thisPassName +#needs recast caloclusters, not (CLUE) ecalclusters +puFinder.cluster_coll_name=pfComb.inputEcalCollName+"Cast" +puFinder.pf_cand_coll_name=pfComb.outputCollName +puFinder.min_momentum=3000. + +# Load pileup finder +from LDMX.Recon.pfReco import pileupFinder + # Load the DQM modules from LDMX.DQM import dqm @@ -191,7 +228,6 @@ ecalVeto, ecalMip, ecal_veto_pnet, - ecal_cluster.EcalClusterProducer(), hcal_digi_reco, hcal_veto, *ts_digis, @@ -204,6 +240,15 @@ p.sequence.extend(dqm_with_overlay) +# Add PFlow + pileup finding sequence +p.sequence.extend([ + cluster, + trackPF, + truthPF, + pfComb, + puFinder +)] + p.inputFiles = ['ecal_pn.root'] p.outputFiles= ['events.root'] p.histogramFile = 'hist.root' diff --git a/DQM/include/DQM/EcalClusterAnalyzer.h b/DQM/include/DQM/EcalClusterAnalyzer.h index 4cfb6cd917..68270e5585 100644 --- a/DQM/include/DQM/EcalClusterAnalyzer.h +++ b/DQM/include/DQM/EcalClusterAnalyzer.h @@ -61,6 +61,13 @@ class EcalClusterAnalyzer : public framework::Analyzer { // Pass Name for clusters std::string cluster_pass_name_; + // Collection Name for Scoring Plane hits + std::string ecal_sp_hits_coll_name_; + // Pass Name for Scoring Plane hits + std::string ecal_sp_hits_pass_name_; + + // min energy fraction from smaller contributor to consider hit "mixed" + double mixed_hit_cutoff_; std::string ecal_sp_hits_passname_; // Inverse skim flag diff --git a/DQM/python/dqm.py b/DQM/python/dqm.py index 99ff66ce28..3f3651f878 100644 --- a/DQM/python/dqm.py +++ b/DQM/python/dqm.py @@ -1130,7 +1130,9 @@ def __init__(self,name='EcalClusterAnalyzer') : self.cluster_coll_name = 'EcalClusters' self.cluster_pass_name = '' - self.ecal_sp_hits_passname = '' + self.ecal_sp_hits_coll_name = 'EcalScoringPlaneHits' + self.ecal_sp_hits_pass_name = '' + self.mixed_hit_cutoff = 0.05 self.inverse_skim = False @@ -1149,20 +1151,27 @@ def __init__(self,name='EcalClusterAnalyzer') : self.build1DHistogram("same_ancestor", "Percentage of hits in cluster coming from the electron that produced most hits", 21, 0., 105.) self.build1DHistogram("energy_percentage", "Percentage of energy in cluster coming from the electron that produced most of energy", 21, 0., 105.) self.build1DHistogram("mixed_hit_energy", "Percentage of total energy coming from hits with energy contributions from more than one electron", 21, 0., 105.) - self.build1DHistogram("clusterless_hits", "Number of hits not in a cluster", 10, 0., 200.) - self.build1DHistogram("clusterless_hits_percentage", "Percentage of hits not in a cluster", 21, 0., 105.) - self.build1DHistogram("total_rechits_in_event", "Number of RecHits", 20, 0., 500.) - + self.build1DHistogram("unclustered_hits", "Number of hits not in a cluster", 10, 0., 200.) + self.build1DHistogram("unclustered_hits_percentage", "Percentage of hits not in a cluster", 21, 0., 105.) + self.build1DHistogram("total_rechits_in_event", "RecHits per event", 20, 0., 500.) - self.build2DHistogram("total_energy_vs_hits", "Total energy (edep) [MeV]", 30, 0, 150, "Hits in cluster", 20, 0, 200) - self.build2DHistogram("total_energy_vs_purity", "Total energy (edep) [MeV]", 30, 0, 150, "Energy purity %", 21, 0, 105) - self.build2DHistogram("sp_ele_distance_vs_purity", "SP ele distance in xy-plane [mm]", 50, 0, 250, "Energy purity %", 21, 0, 105) + self.build2DHistogram("total_energy_vs_hits", "Total energy (edep) [MeV]", 30, 0., 150., "Hits in cluster", 30, 0., 300.) + self.build2DHistogram("total_energy_vs_purity", "Total energy (edep) [MeV]", 30, 0., 150., "Energy purity %", 21, 0, 105.) + self.build2DHistogram("sp_ele_distance_vs_purity", "SP ele distance in xy-plane [mm]", 50, 0, 250, "Energy purity %", 21, 0., 105.) self.build2DHistogram("sp_clue_distance_vs_layer", "CLUE centroid to SP ele distance in xy-plane [mm]", 125, 0., 250., "Layer", 33, -0.5, 32.5) self.build1DHistogram("sp_clue_distance", "CLUE centroid to SP ele distance in xy-plane [mm]", 125, 0., 250.) self.build1DHistogram("sp_clue_x_residual", "CLUE centroid X - SP ele X [mm]", 250, -250., 250.) self.build1DHistogram("sp_clue_y_residual", "CLUE centroid Y - SP ele Y [mm]", 250, -250., 250.) + self.build1DHistogram("sp_distance", "dR(SPhit_1, SPhit_2)", 100, -1, 202) + self.build1DHistogram("cluster_distance", "dR(cl_1, cl_2)", 100, -1, 202) + self.build1DHistogram("cluster_RMSX", "RMS(hits in cluster) X", 100, -1, 202) + self.build1DHistogram("cluster_RMSY", "RMS(hits in cluster) Y", 100, -1, 202) + self.build2DHistogram("dE_cl2_vs_cl1", "E_{cl}-E_{true}^{SP}, cluster 1 [MeV]", 100, -10000, 10000, "E_{cl}-E_{true}^{SP}, cluster 2 [MeV]", 100, -10000, 10000) + + self.build2DHistogram("tag0frac_vs_SPdist", "dR(SPhit_1, SPhit_2)", 251, -1, 250, "Fraction of mixed (purity less than {int(100*(1.-self.mixed_hit_cutoff))}%) ancestors", 200, 0, 1) + ecal_dqm = [ EcalDigiVerify(), EcalShowerFeatures(), diff --git a/DQM/src/DQM/EcalClusterAnalyzer.cxx b/DQM/src/DQM/EcalClusterAnalyzer.cxx index 4d3527e0b7..3a64e54dd2 100644 --- a/DQM/src/DQM/EcalClusterAnalyzer.cxx +++ b/DQM/src/DQM/EcalClusterAnalyzer.cxx @@ -16,6 +16,12 @@ void EcalClusterAnalyzer::configure(framework::config::Parameters& ps) { cluster_coll_name_ = ps.get("cluster_coll_name"); cluster_pass_name_ = ps.get("cluster_pass_name"); + ecal_sp_hits_coll_name_ = + ps.getParameter("ecal_sp_hits_coll_name"); + ecal_sp_hits_pass_name_ = + ps.getParameter("ecal_sp_hits_pass_name"); + mixed_hit_cutoff_ = ps.getParameter("mixed_hit_cutoff"); + ecal_sp_hits_passname_ = ps.get("ecal_sp_hits_passname"); inverse_skim_ = ps.get("inverse_skim"); n_ecal_clusters_min_ = ps.get("n_ecal_clusters_min"); @@ -79,7 +85,7 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { // Determine the truth information for the recoil electron std::vector> sp_electron_positions; const auto& ecal_sp_hits{event.getCollection( - "EcalScoringPlaneHits", ecal_sp_hits_passname_)}; + "EcalScoringPlaneHits", ecal_sp_hits_pass_name_)}; std::vector sorted_sp_hits = ecal_sp_hits; std::sort(sorted_sp_hits.begin(), sorted_sp_hits.end(), @@ -114,6 +120,8 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { << ", TS counted electrons: " << nbr_of_electrons << ", SP electrons: " << sp_electron_positions.size(); + std::map true_energy; + std::map delta_energy; double sp_ele_dist{9999.}; if (nbr_of_electrons == 2 && sp_electron_positions.size() > 1) { // Measures sp_ele_distance between two electrons in the ECal scoring plane @@ -125,11 +133,18 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { sp_ele_dist = std::sqrt((pos1[0] - pos2[0]) * (pos1[0] - pos2[0]) + (pos1[1] - pos2[1]) * (pos1[1] - pos2[1])); + histograms_.fill("sp_distance", sp_ele_dist); + } // end block about the scoring plane hits ldmx_log(trace) << "Distance between the two e- in the ECal scoring plane: " << sp_ele_dist << " mm"; + double tot_event_energy = 0; + std::vector tot_origin_edep; + tot_origin_edep.resize(nbr_of_electrons_ + 1); + int n_mixed = 0; + // Loop over the rechits and find the matching simhits ldmx_log(trace) << "Loop over the rechits and find the matching simhits"; for (const auto& hit : ecal_rec_hits) { @@ -138,31 +153,65 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { [&hit](const auto& sim_hit) { return sim_hit.getID() == hit.getID(); }); if (it != ecal_sim_hits.end()) { // if found a simhit matching this rechit + ldmx_log(trace) << "\tFound simhit matching rechit with ID" + << hit.getID(); int ancestor = 0; int prev_ancestor = 0; bool tagged = false; int tag = 0; std::vector edep; - edep.resize(nbr_of_electrons + 1); + edep.resize(nbr_of_electrons_ + 1); + double e_tot = 0; // keep track of total from all counted ancestors + ldmx_log(trace) << "\t\tIt has " << it->getNumberOfContribs() + << " contribs. "; for (int i = 0; i < it->getNumberOfContribs(); i++) { // for each contrib in this simhit const auto& contrib = it->getContrib(i); // get origin electron ID ancestor = contrib.origin_id_; + ldmx_log(trace) << "\t\t\tAncestor ID " << ancestor << " with edep " + << contrib.edep_; + tot_event_energy += contrib.edep_; // store energy from this contrib at index = origin electron ID - if (ancestor <= nbr_of_electrons) edep[ancestor] += contrib.edep_; + if (ancestor <= nbr_of_electrons) { + edep[ancestor] += contrib.edep_; + tot_origin_edep[ancestor] += contrib.edep_; + e_tot += contrib.edep_; + } if (!tagged && i != 0 && prev_ancestor != ancestor) { // if origin electron ID does not match previous origin electron ID // this hit has contributions from several electrons, ie mixed case tag = 0; tagged = true; + ldmx_log(trace) << "\t\t\t\tMixed hit! Ancestor ID changed to " + << ancestor; } prev_ancestor = ancestor; + } // over contribs + // now check if mixed really means mixed, i.e. more than small fraction + // from a second electron. + if (tagged) { + for (int i = 1; i < nbr_of_electrons_ + 1; i++) { + if (edep[i] / e_tot > + 1 - mixed_hit_cutoff_) { // one ancestor contributes at least the + // complement to the allowed mixing + // fraction + tagged = false; + ancestor = + distance(edep.begin(), max_element(edep.begin(), edep.end())); + ldmx_log(trace) + << "\t\t\t\tUndid mixed hit tagging, now ancestor = " + << ancestor; + break; + } + } } if (!tagged) { - // if not tagged, hit was from a single electron - tag = prev_ancestor; - } + // if not tagged, hit was from a single electron (within acceptable + // purity) + tag = ancestor; // prev_ancestor; + } else + n_mixed++; histograms_.fill("ancestors", tag); hit_info.insert({hit.getID(), std::make_pair(tag, edep)}); } // end if simhit found @@ -171,6 +220,21 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { // Loop over the clusters int clustered_hits = 0; ldmx_log(trace) << "Loop over the clusters, N = " << n_ecal_clusters; + histograms_.fill("tag0frac_vs_SPdist", sp_ele_dist, + (float)n_mixed / ecal_rec_hits.size()); + ldmx_log(debug) << "Got " << n_mixed << " mixed hits, a fraction of " + << (float)n_mixed / ecal_rec_hits.size(); + + if (ecal_clusters.size() >= 2) { + float d_x = + ecal_clusters[0].getCentroidX() - ecal_clusters[1].getCentroidX(); + float d_y = + ecal_clusters[0].getCentroidY() - ecal_clusters[1].getCentroidY(); + float d_r = std::sqrt(d_x * d_x + d_y * d_y); + histograms_.fill("cluster_distance", d_r); + ldmx_log(trace) << "Gt cluster distance (0,1) = " << d_r; + } + for (const auto& cl : ecal_clusters) { auto layer = cl.getLayer(); ldmx_log(trace) << "Cluster in layer " << layer @@ -217,7 +281,7 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { energy_from_electron.resize(nbr_of_electrons + 1); double energy_sum = 0.; double n_sum = 0.; - + ldmx_log(trace) << "Looping over hits in the cluster"; const auto& hit_ids = cl.getHitIDs(); for (const auto& id : hit_ids) { // for each hit in cluster, find previously stored info @@ -255,6 +319,11 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { // get largest energy contribution double max_energy_contribution = *max_element( energy_from_electron.begin(), energy_from_electron.end()); + std::string to_log; + for (auto nb : energy_from_electron) + to_log.append(std::to_string(nb) + " "); + ldmx_log(debug) << "Energies vector is " << to_log; + // energy purity = largest contribution / all energy histograms_.fill("energy_percentage", 100. * (max_energy_contribution / energy_sum)); @@ -277,12 +346,25 @@ void EcalClusterAnalyzer::analyze(const framework::Event& event) { n_hits_from_electron.end()); histograms_.fill("same_ancestor", 100. * (n_max / n_sum)); } + // find the main contributor + auto elt = distance( + energy_from_electron.begin(), + max_element(energy_from_electron.begin(), energy_from_electron.end())); + ldmx_log(debug) << "Found that the maximum contributing trackID is " << elt; + delta_energy[elt] = cl.getEnergy() - true_energy[elt]; + // delta_energy[2]=e[2]-true_energy[2]; + histograms_.fill("cluster_RMSX", cl.getRMSX()); } // end loop on the clusters + std::string more_log; + for (auto nb : tot_origin_edep) more_log.append(std::to_string(nb) + " "); + ldmx_log(debug) << "Edep per ancestor in event is " << more_log; + ldmx_log(debug) << "Total energy deposited in event: " << tot_event_energy; - histograms_.fill("clusterless_hits", (ecal_rec_hits.size() - clustered_hits)); + histograms_.fill("dE_cl2_vs_cl1", delta_energy[1], delta_energy[2]); + histograms_.fill("unclustered_hits", (ecal_rec_hits.size() - clustered_hits)); histograms_.fill("total_rechits_in_event", ecal_rec_hits.size()); histograms_.fill( - "clusterless_hits_percentage", + "unclustered_hits_percentage", 100. * (ecal_rec_hits.size() - clustered_hits) / ecal_rec_hits.size()); if (inverse_skim_) { diff --git a/Ecal/include/Ecal/CLUE.h b/Ecal/include/Ecal/CLUE.h index 3dc7b14ae8..4c2ab6ac81 100644 --- a/Ecal/include/Ecal/CLUE.h +++ b/Ecal/include/Ecal/CLUE.h @@ -79,6 +79,9 @@ class CLUE { std::vector> setup( const std::vector& hits); + // get distance between clusters in the first layers, proxy for electron sep. + void electronSeparation(std::vector hits); + // connectingLayers marks if we're currently doing 3D clustering (i.e. // connecting seeds between layers) otherwise, layerTag tells us which layer // number we're working on diff --git a/Ecal/src/Ecal/CLUE.cxx b/Ecal/src/Ecal/CLUE.cxx index 169a9001d6..60c4849a2c 100644 --- a/Ecal/src/Ecal/CLUE.cxx +++ b/Ecal/src/Ecal/CLUE.cxx @@ -34,64 +34,66 @@ T CLUE::dist(T x1, T y1, T z1, T x2, T y2, T z2) { likely merged => recluster Did not quite work and I don't remember the idea anymore but leaving the code here for inspo */ -// void electronSeparation(std::vector hits_) { -// std::vector layerThickness = -// { 2., 3.5, 5.3, 5.3, 5.3, 5.3, 5.3, 5.3, 5.3, 5.3, 5.3, 10.5, 10.5, 10.5, 10.5, -// 10.5 }; double air = 10.; std::sort(hits_.begin(), hits_.end(), [](const -// ldmx::EcalHit& a, const ldmx::EcalHit& b) { -// return a.getZPos() < b.getZPos(); -// }); -// std::vector firstLayers; -// std::vector firstLayerClusters; -// int layer_index = 0; -// double layerZ = hits_[0].getZPos(); -// for (const auto& hit : hits_) { -// if (hit.getZPos() > layerZ + layerThickness[layer_index] + air) { -// layer_index++; -// // if (layer_index > limit) break; -// break; -// } -// firstLayers.push_back(hit); -// firstLayerClusters.push_back(IntermediateEcalCluster(hit, layer_index)); - -// } -// bool merge = false; -// do { -// merge = false; -// for (int i = 0; i < firstLayerClusters.size(); i++) { -// if (firstLayerClusters[i].empty()) continue; -// // if (firstLayerClusters[i].centroid().E() >= seed_threshold_) { -// for (int j = i + 1; j < firstLayerClusters.size(); j++) { -// if (firstLayerClusters[j].empty()) continue; -// if (dist(firstLayerClusters[i].centroid().Px(), -// firstLayerClusters[i].centroid().Py(), -// firstLayerClusters[j].centroid().Px(), -// firstLayerClusters[j].centroid().Py()) < 8.) { -// firstLayerClusters[i].add(firstLayerClusters[j]); -// firstLayerClusters[j].clear(); -// merge = true; -// } -// } -// // } else break; -// } -// } while (merge); -// ldmx_log(trace) << "--- ELECTRON SEPARATION ---"; -// for (int i = 0; i < firstLayerClusters.size(); i++) { -// if (firstLayerClusters[i].empty()) continue; -// ldmx_log(trace) << " Cluster " << i << " x_: " -// << firstLayerClusters[i].centroid().Px() << " y_: " -// << firstLayerClusters[i].centroid().Py(); -// for (int j = i + 1; -// j < firstLayerClusters.size(); j++) { -// if (firstLayerClusters[j].empty()) continue; -// auto d = dist(firstLayerClusters[i].centroid().Px(), -// firstLayerClusters[i].centroid().Py(), -// firstLayerClusters[j].centroid().Px(), -// firstLayerClusters[j].centroid().Py()); -// ldmx_log(trace) << "Dist to cluster " << j << ": " << d; -// } -// } -// } +void CLUE::electronSeparation(std::vector hits) { + std::vector layer_thickness = {2., 3.5, 5.3, 5.3, 5.3, 5.3, + 5.3, 5.3, 5.3, 5.3, 5.3, 10.5, + 10.5, 10.5, 10.5, 10.5}; + double air = 10.; + // sort hits in z + std::sort(hits.begin(), hits.end(), + [](const ldmx::EcalHit& a, const ldmx::EcalHit& b) { + return a.getZPos() < b.getZPos(); + }); + + std::vector first_layers; + std::vector first_layer_clusters; + int layer_tag = 0; + double layer_z = hits[0].getZPos(); + for (const auto& hit : hits) { + if (hit.getZPos() > layer_z + layer_thickness[layer_tag] + air) { + layer_tag++; + // if (layerTag > limit) break; + break; + } + first_layers.push_back(hit); + first_layer_clusters.push_back(IntermediateCluster(hit, layer_tag)); + } + bool merge = false; + do { + merge = false; + for (int i = 0; i < first_layer_clusters.size(); i++) { + if (first_layer_clusters[i].empty()) continue; + // if (firstLayerClusters[i].centroid().E() >= seedThreshold_) { + for (int j = i + 1; j < first_layer_clusters.size(); j++) { + if (first_layer_clusters[j].empty()) continue; + if (dist(first_layer_clusters[i].centroid().Px(), + first_layer_clusters[i].centroid().Py(), + first_layer_clusters[j].centroid().Px(), + first_layer_clusters[j].centroid().Py()) < 8.) { + first_layer_clusters[i].add(first_layer_clusters[j]); + first_layer_clusters[j].clear(); + merge = true; + } + } + // } else break; + } + } while (merge); + ldmx_log(trace) << "--- ELECTRON SEPARATION ---"; + for (int i = 0; i < first_layer_clusters.size(); i++) { + if (first_layer_clusters[i].empty()) continue; + ldmx_log(trace) << " Cluster " << i + << " x: " << first_layer_clusters[i].centroid().Px() + << " y: " << first_layer_clusters[i].centroid().Py(); + for (int j = i + 1; j < first_layer_clusters.size(); j++) { + if (first_layer_clusters[j].empty()) continue; + auto d = dist(first_layer_clusters[i].centroid().Px(), + first_layer_clusters[i].centroid().Py(), + first_layer_clusters[j].centroid().Px(), + first_layer_clusters[j].centroid().Py()); + ldmx_log(trace) << "Dist to cluster " << j << ": " << d; + } + } +} // Function to create layers from hits based on their layer number std::vector> CLUE::createLayers( @@ -162,7 +164,9 @@ std::vector> CLUE::setup( // collapse z dimension float x = roundToDecimal(hit->getXPos(), 4); float y = roundToDecimal(hit->getYPos(), 4); - ldmx_log(trace) << " New hit { x: " << x << " y: " << y << "}"; + float z = roundToDecimal(hit->getZPos(), 4); + ldmx_log(trace) << " New hit { x: " << x << " y: " << y << "}" + << " (and z: " << z << ")"; std::pair coords; if (dc_ != 0 && nbr_of_layers_ > 1) { // if more than one layer, divide hit into densities with side dc @@ -188,11 +192,13 @@ std::vector> CLUE::setup( // density coords = {x, y}; } + if (density_map.find(coords) == density_map.end()) { density_map.emplace(coords, std::make_shared(x, y)); - ldmx_log(trace) << " New density created"; + ldmx_log(trace) << " * New density created"; } else { - ldmx_log(trace) << " Found density with x: " << density_map[coords]->x_ + ldmx_log(trace) << " --> Found density with x: " + << density_map[coords]->x_ << " y: " << density_map[coords]->y_; } density_map[coords]->hits_.push_back(hit); @@ -334,8 +340,13 @@ std::vector> CLUE::clustering( density->total_energy_ > rhoc_ && density->delta_ > delta_c_mod; } else { is_seed = density->total_energy_ > rhoc_ && density->delta_ > deltac_; + if (is_seed) { + ldmx_log(trace) << " Distance to event centroid: " + << dist(density->x_, density->y_, + event_centroid_.centroid().x(), + event_centroid_.centroid().y()); + } } - bool is_outlier = (density->total_energy_ < rhoc_) && (density->delta_ > deltao_); density->cluster_id_ = -1; @@ -359,6 +370,10 @@ std::vector> CLUE::clustering( int& parent_index = density->follower_of_; if (parent_index != -1) followers[parent_index].push_back(density->index_); + else + ldmx_log(error) + << " Somehow found a follower with parent index -1: id = " + << density->index_; } else { ldmx_log(trace) << " This is an Outlier"; } @@ -446,7 +461,9 @@ std::vector> CLUE::setupForClue3D() { if (current_seed->total_energy_ > highest_energy) highest_energy = current_seed->total_energy_; ldmx_log(trace) << " Density with index " << current_seed->index_ - << ", energy: " << current_seed->total_energy_; + << ", energy: " << current_seed->total_energy_ + << " position (x,y)= {" << current_seed->x_ << "," + << current_seed->y_ << ")"; int depth = 1; // decide delta_ and followerof from seeds in previous and next layer_ // do { @@ -454,7 +471,9 @@ std::vector> CLUE::setupForClue3D() { if ((layer - depth >= 0) && (layer - depth < seeds_.size())) { ldmx_log(trace) << " Looking at pre-layer: " << layer - depth; // look at previous layer + ldmx_log(trace) << " In previous layer... "; auto& previous_layer = seeds_[layer - depth]; + ldmx_log(trace) << " Got " << previous_layer.size() << " seeds"; for (const auto& prev_seed : previous_layer) { // for each seed in previous layer auto distance_2d_prev = dist(current_seed->x_, current_seed->y_, @@ -484,9 +503,11 @@ std::vector> CLUE::setupForClue3D() { } } } + if (layer + depth < nbr_of_layers_ && layer + depth < seeds_.size()) { ldmx_log(trace) << " Looking at post-layer: " << layer + depth; auto& next_layer = seeds_[layer + depth]; + ldmx_log(trace) << " Got " << next_layer.size() << " seeds"; for (const auto& next_seed : next_layer) { auto distance_2d_next = dist(current_seed->x_, current_seed->y_, next_seed->x_, next_seed->y_); @@ -517,6 +538,7 @@ std::vector> CLUE::setupForClue3D() { } // end of looking at next layer // } while (currentLayer[i]->layerFollowerOf == -1 && (layer - depth >= // 0 || layer + depth < nbr_of_layers_)); + ldmx_log(trace) << " Done setting parents."; densities.push_back(current_seed); } // TODO: This 2 needs to be configurable @@ -601,7 +623,7 @@ void CLUE::cluster(const std::vector& unsorted_hits, double dc, return a->getZPos() < b->getZPos(); }); - seeds_.resize(nbr_of_layers); + seeds_.resize(nbr_of_layers_); if (nbr_of_layers_ > 1) { ldmx_log(debug) << "Creating layers"; diff --git a/Recon/include/Recon/Event/CaloCluster.h b/Recon/include/Recon/Event/CaloCluster.h index 3056c5becf..648b7ae30d 100644 --- a/Recon/include/Recon/Event/CaloCluster.h +++ b/Recon/include/Recon/Event/CaloCluster.h @@ -135,7 +135,7 @@ class CaloCluster { /// @brief Delta unc on unc in y-z plane double getEDYDZ() const { return err_dydz_; } - // get hit rawIDs (unused) + // get hit rawIDs const std::vector& getHitIDs() const { return hit_ids_; } // ability to store limited hit info diff --git a/Recon/include/Recon/Event/PFCandidate.h b/Recon/include/Recon/Event/PFCandidate.h index ce01f1871a..e6bd1af43d 100644 --- a/Recon/include/Recon/Event/PFCandidate.h +++ b/Recon/include/Recon/Event/PFCandidate.h @@ -9,6 +9,9 @@ // ROOT #include "TObject.h" //For ClassDef +// ldmx-sw objects +// #include "Ecal/Event/EcalHit.h" +// #include "Hcal/Event/HcalHit.h" namespace ldmx { @@ -49,12 +52,15 @@ class PFCandidate { pos_hcal_y_ = y_; pos_hcal_z_ = z_; } - void setTrackPxPyPz(float x_, float y_, float z_) { track_px_ = x_; track_py_ = y_; track_pz_ = z_; } + // associate component indices to the pf candidate + void setTrackIndex(int x) { track_idx_ = x; } + void setEcalIndex(int x) { ecal_idx_ = x; } + void setHcalIndex(int x) { hcal_idx_ = x; } void setEcalEnergy(float x_) { ecal_energy_ = x_; } void setEcalRawEnergy(float x_) { ecal_raw_energy_ = x_; } @@ -120,6 +126,10 @@ class PFCandidate { std::vector getHcalPositionXYZ() const { return {pos_hcal_x_, pos_hcal_y_, pos_hcal_z_}; } + // associate component indices to the pf candidate + int getTrackIndex() const { return track_idx_; } + int getEcalIndex() const { return ecal_idx_; } + int getHcalIndex() const { return hcal_idx_; } std::vector getTrackPxPyPz() const { return {track_px_, track_py_, track_pz_}; @@ -224,6 +234,11 @@ class PFCandidate { double truth_energy_{0}; int truth_pdg_id_{0}; + /* Indices of the components making up the PFlow object */ + int track_idx_{-1}; + int ecal_idx_{-1}; + int hcal_idx_{-1}; + /* The ROOT class definition. */ ClassDef(PFCandidate, 2); }; diff --git a/Recon/include/Recon/PileupFinder.h b/Recon/include/Recon/PileupFinder.h new file mode 100644 index 0000000000..1fae41151a --- /dev/null +++ b/Recon/include/Recon/PileupFinder.h @@ -0,0 +1,55 @@ +/** + * @file PileupFinder.h + * @brief Simple PFlow algorithm + * @author Christian Herwig, Fermilab + */ + +#ifndef PARTICLEFLOW_H +#define PARTICLEFLOW_H + +// LDMX Framework +#include "Ecal/Event/EcalCluster.h" +#include "Framework/Configure/Parameters.h" // Needed to import parameters from configuration file +#include "Framework/Event.h" +#include "Framework/EventProcessor.h" //Needed to declare processor +#include "Hcal/Event/HcalCluster.h" +#include "Recon/Event/CaloCluster.h" +#include "Recon/Event/PFCandidate.h" +#include "SimCore/Event/SimParticle.h" +#include "SimCore/Event/SimTrackerHit.h" +#include "TGraph.h" + +namespace recon { + +/** + * @class PileupFinder + * @brief + */ +class PileupFinder : public framework::Producer { + public: + PileupFinder(const std::string& name, framework::Process& process) + : framework::Producer(name, process) {} + + virtual void configure(framework::config::Parameters& ps); + + virtual void produce(framework::Event& event); + + virtual void onProcessEnd(); + + private: + // name of collections for PF input object to be passed + std::string rec_hit_coll_name_; + std::string rec_hit_pass_name_; + std::string pf_cand_coll_name_; + std::string pf_cand_pass_name_; + std::string cluster_coll_name_; + std::string cluster_pass_name_; + // name of collection for pileup-free output hit coll + std::string output_rec_hit_coll_name_; + // configuration + + double min_mom_{0.}; // MeV +}; +} // namespace recon + +#endif /* PILEUPFINDER_H */ diff --git a/Recon/python/pileupFinder.py b/Recon/python/pileupFinder.py new file mode 100644 index 0000000000..e9d28dc873 --- /dev/null +++ b/Recon/python/pileupFinder.py @@ -0,0 +1,24 @@ +"""Configuration for pfReco producers + +Sets all parameters to reasonable defaults. + +Examples +-------- + from LDMX.Recon.pfReco import pileupFinder + p.sequence.append( pileupFinder ) +""" + +from LDMX.Framework import ldmxcfg + +class pileupFinder(ldmxcfg.Producer) : + """Configuration for pileup finding from particle flow objects""" + def __init__(self, name='PileupFinder') : + super().__init__(name, 'recon::PileupFinder','Recon') + self.rec_hit_coll_name = 'EcalRecHits' + self.rec_hit_pass_name = '' + self.cluster_coll_name = 'PFEcalClusters' + self.cluster_pass_name = '' + self.pf_cand_coll_name = 'PFCandidates' + self.pf_cand_pass_name = '' + self.output_rec_hit_coll_name = 'EcalRecHitsNoPileup' + self.min_momentum = 4000. diff --git a/Recon/src/Recon/ParticleFlow.cxx b/Recon/src/Recon/ParticleFlow.cxx index aec3a2ea72..b5cb7c0495 100644 --- a/Recon/src/Recon/ParticleFlow.cxx +++ b/Recon/src/Recon/ParticleFlow.cxx @@ -313,10 +313,12 @@ void ParticleFlow::produce(framework::Event& event) { ldmx::PFCandidate cand; fillCandTrack(cand, tracks[i]); // append track info to candidate + cand.setTrackIndex(i); if (!tk_is_em_linked[i]) { // chargedUnmatch.push_back(cand); } else { // if track is linked with ECal cluster fillCandEMCalo(cand, ecal_clusters[tk_em_pairs[i]]); + cand.setEcalIndex(tk_em_pairs[i]); if (em_is_had_linked[tk_em_pairs[i]]) { // if ECal is linked with HCal // cluster fillCandHadCalo(cand, hcal_clusters[em_had_pairs[tk_em_pairs[i]]]); @@ -331,9 +333,9 @@ void ParticleFlow::produce(framework::Event& event) { for (int i = 0; i < ecal_clusters.size(); i++) { // already linked with ECal in the previous step if (em_is_tk_linked[i]) continue; - ldmx::PFCandidate cand; fillCandEMCalo(cand, ecal_clusters[i]); + cand.setEcalIndex(i); if (em_is_had_linked[tk_em_pairs[i]]) { fillCandHadCalo(cand, hcal_clusters[em_had_pairs[i]]); // emMatch.push_back(cand); @@ -347,6 +349,7 @@ void ParticleFlow::produce(framework::Event& event) { if (had_is_em_linked[i]) continue; ldmx::PFCandidate cand; fillCandHadCalo(cand, hcal_clusters[i]); + cand.setHcalIndex(i); // hadOnly.push_back(cand); pf_cands.push_back(cand); } diff --git a/Recon/src/Recon/PileupFinder.cxx b/Recon/src/Recon/PileupFinder.cxx new file mode 100644 index 0000000000..66f480b651 --- /dev/null +++ b/Recon/src/Recon/PileupFinder.cxx @@ -0,0 +1,104 @@ +#include "Recon/PileupFinder.h" + +#include + +namespace recon { + +void PileupFinder::configure(framework::config::Parameters& ps) { + // I/O + rec_hit_coll_name_ = ps.getParameter("rec_hit_coll_name"); + rec_hit_pass_name_ = ps.getParameter("rec_hit_pass_name"); + pf_cand_coll_name_ = ps.getParameter("pf_cand_coll_name"); + pf_cand_pass_name_ = ps.getParameter("pf_cand_pass_name"); + cluster_coll_name_ = ps.getParameter("cluster_coll_name"); + cluster_pass_name_ = ps.getParameter("cluster_pass_name"); + output_rec_hit_coll_name_ = + ps.getParameter("output_rec_hit_coll_name"); + // Algorithm configuration + min_mom_ = ps.getParameter("min_momentum"); +} + +// get pileup candidates from PFlow and make a cleaned-up hit collection +void PileupFinder::produce(framework::Event& event) { + if (!event.exists(rec_hit_coll_name_, rec_hit_pass_name_)) { // ecal rechits + ldmx_log(error) << "Unable to find (one) collection named " + << rec_hit_coll_name_ << "_" << rec_hit_pass_name_; + return; + } + if (!event.exists(pf_cand_coll_name_, pf_cand_pass_name_)) { + ldmx_log(error) << "Unable to find (one) collection named " + << pf_cand_coll_name_ << "_" << pf_cand_pass_name_; + return; + } + if (!event.exists(cluster_coll_name_, cluster_pass_name_)) { + ldmx_log(error) << "Unable to find (one) collection named " + << cluster_coll_name_ << "_" << cluster_pass_name_; + return; + } + + const auto& ecal_hits{event.getCollection(rec_hit_coll_name_, + rec_hit_pass_name_)}; + + const auto& pf_cands{event.getCollection( + pf_cand_coll_name_, pf_cand_pass_name_)}; + + const auto& clusters{event.getCollection( + cluster_coll_name_, cluster_pass_name_)}; + // get PID 3 and 7 -- the ones with track and ecal matching + // get the high-momentum track ones from there -- pileup candidate! + // get the list of hits associated with pileup candidates + // if a rechit is not on that list, add to output collection. + std::vector output_hits; + std::vector pileup_hit_i_ds; + + // this needs to be a two-step procedure: loop over all clusters deemed to be + // pileup to find all their associated hits + // then loop over that list to make a collection of output hits that doesn't + // contain them + for (const auto& pf_cand : pf_cands) { + if (pf_cand.getPID() == 3 || pf_cand.getPID() == 7) { + // we have both ecal cluster and track + std::vector mom_vec = pf_cand.getTrackPxPyPz(); + float mom = mom_vec[0] * mom_vec[0] + mom_vec[1] * mom_vec[1] + + mom_vec[2] * mom_vec[2]; + mom = sqrt(mom); + + if (mom < min_mom_) continue; + ldmx_log(trace) << "Got pileup candidate with PID = " << pf_cand.getPID() + << " and momentum = " << mom << " MeV."; + + // now! use the hit-candidate association to get the associated ecal hits. + int pf_cl_idx = pf_cand.getEcalIndex(); + ldmx_log(trace) << "Got Ecal cluster with index " << pf_cl_idx + << " while cluster array length is " << clusters.size(); + if (pf_cl_idx < 0) // was never set + continue; + auto cl = clusters[pf_cl_idx]; + auto hit_i_ds = cl.getHitIDs(); + // add to collection of pileup hits + pileup_hit_i_ds.insert(pileup_hit_i_ds.end(), hit_i_ds.begin(), + hit_i_ds.end()); + } // if trk/ecal matched + } // over PF objects + + for (auto hit : ecal_hits) { + auto found_index = std::find(std::begin(pileup_hit_i_ds), + std::end(pileup_hit_i_ds), hit.getID()); + // When the element is not found, std::find returns the end of the range + if (found_index == + std::end(pileup_hit_i_ds)) { // hit not part of any pileup cluster + output_hits.emplace_back(hit); // keep it + ldmx_log(trace) << "Got no-pileup hit! "; + } + } + event.add(output_rec_hit_coll_name_, output_hits); +} +void PileupFinder::onProcessEnd() { + ldmx_log(debug) << "Process ends!"; + + return; +} + +} // namespace recon + +DECLARE_PRODUCER(recon::PileupFinder); diff --git a/TrigScint/include/TrigScint/TruthHitProducer.h b/TrigScint/include/TrigScint/TruthHitProducer.h index e900d9de32..e7cb6b1090 100644 --- a/TrigScint/include/TrigScint/TruthHitProducer.h +++ b/TrigScint/include/TrigScint/TruthHitProducer.h @@ -52,25 +52,20 @@ class TruthHitProducer : public framework::Producer { */ void produce(framework::Event &event) override; - /// Class to set the verbosity level. - // TODO: Make use of the global verbose parameter. - bool verbose_{false}; - + private: /// Name of the input collection containing the sim hits std::string input_collection_; /// Name of the pass that the input collection is on (empty string means take /// any pass) std::string input_pass_name_; + /// Name of the pass that the input simparticles collection is on + std::string sim_particles_pass_name_; /// Name of the output collection that will be used to store the /// selected sim hits std::string output_collection_; - private: - std::string sim_particles_passname_; - std::string input_collection_events_passname_; - }; // TruthHitProducer } // namespace trigscint diff --git a/TrigScint/python/truthHits.py b/TrigScint/python/truthHits.py index e9d9023493..933f9dd355 100644 --- a/TrigScint/python/truthHits.py +++ b/TrigScint/python/truthHits.py @@ -20,10 +20,7 @@ def __init__(self,name) : self.input_collection="TriggerPadUpSimHits" self.input_pass_name="" #take any pass self.output_collection="truthBeamElectronsUp" - self.verbose = False - - self.sim_particles_passname = "" - self.input_collection_events_passname = "" + self.sim_particles_pass_name = "" def up() : """Get the beam electron truth hits for the trigger pad upstream of target""" diff --git a/TrigScint/src/TrigScint/TruthHitProducer.cxx b/TrigScint/src/TrigScint/TruthHitProducer.cxx index dbff995750..80ead80902 100644 --- a/TrigScint/src/TrigScint/TruthHitProducer.cxx +++ b/TrigScint/src/TrigScint/TruthHitProducer.cxx @@ -8,38 +8,38 @@ TruthHitProducer::TruthHitProducer(const std::string &name, : Producer(name, process) {} void TruthHitProducer::configure(framework::config::Parameters ¶meters) { - input_collection_ = parameters.get("input_collection"); - input_pass_name_ = parameters.get("input_pass_name"); - output_collection_ = parameters.get("output_collection"); - sim_particles_passname_ = - parameters.get("sim_particles_passname"); - input_collection_events_passname_ = - parameters.get("input_collection_events_passname"); + input_collection_ = parameters.getParameter("input_collection"); + input_pass_name_ = parameters.getParameter("input_pass_name"); + output_collection_ = + parameters.getParameter("output_collection"); + sim_particles_pass_name_ = + parameters.getParameter("sim_particles_pass_name"); - verbose_ = parameters.get("verbose"); - - if (verbose_) { - ldmx_log(info) << "In TruthHitProducer: configure done!"; - ldmx_log(info) << "Got parameters: " << "\nInput collection: " - << input_collection_ - << "\nInput pass name: " << input_pass_name_ - << "\nOutput collection: " << output_collection_ - << "\nVerbose: " << verbose_; - } + ldmx_log(info) << "In TruthHitProducer: configure done!"; + ldmx_log(info) << "Got parameters: " << "\nInput collection: " + << input_collection_ + << "\nInput pass name: " << input_pass_name_ + << "\nOutput collection: " << output_collection_; } void TruthHitProducer::produce(framework::Event &event) { // Check if the collection exists. If not, don't bother processing the event. - if (!event.exists(input_collection_, input_collection_events_passname_)) { - ldmx_log(error) << "No input collection called " << input_collection_ - << " found; skipping!"; + if (!event.exists(input_collection_, input_pass_name_)) { + ldmx_log(error) << "No input collection called " << input_collection_ << "_" + << input_pass_name_ << " found; skipping!"; + return; + } + if (!event.exists("SimParticles", sim_particles_pass_name_)) { + ldmx_log(error) << "No input SimParticle collection with pass " + << sim_particles_pass_name_ << " found; skipping!"; return; } + // looper over sim hits and aggregate energy depositions for each detID const auto sim_hits{event.getCollection( input_collection_, input_pass_name_)}; auto particle_map{event.getMap( - "SimParticles", sim_particles_passname_)}; + "SimParticles", sim_particles_pass_name_)}; std::vector truth_beam_electrons; @@ -49,15 +49,14 @@ void TruthHitProducer::produce(framework::Event &event) { // check if hit is from beam electron and, if so, add to output collection for (int i = 0; i < sim_hit.getNumberOfContribs(); i++) { auto contrib = sim_hit.getContrib(i); - if (verbose_) { - ldmx_log(debug) << "contrib " << i << " trackID: " << contrib.track_id_ - << " pdgID: " << contrib.pdg_code_ - << " edep: " << contrib.edep_; - ldmx_log(debug) << "\t particle id: " - << particle_map[contrib.track_id_].getPdgID() - << " particle status: " - << particle_map[contrib.track_id_].getGenStatus(); - } + ldmx_log(trace) << "contrib " << i << " track_id: " << contrib.track_id_ + << " pdgID: " << contrib.pdg_code_ + << " edep: " << contrib.edep_; + ldmx_log(trace) << "\t particle id: " + << particle_map[contrib.track_id_].getPdgID() + << " particle status: " + << particle_map[contrib.track_id_].getGenStatus(); + // if the trackID is in the map if (particle_map.find(contrib.track_id_) != particle_map.end()) { // beam electron (PDGID = 11, genStatus == 1) @@ -67,8 +66,8 @@ void TruthHitProducer::produce(framework::Event &event) { } } if (keep) truth_beam_electrons.push_back(sim_hit); - } - } + } // over simhit contribs + } // over simhits event.add(output_collection_, truth_beam_electrons); } } // namespace trigscint