From 44ac16a1b43c49840fab9c5d479f2f8c72853a89 Mon Sep 17 00:00:00 2001 From: tillrue Date: Sun, 8 Jun 2025 11:23:13 +0200 Subject: [PATCH 1/6] Add wide angle events in generators.py Added wide angle event generator. Shots single (mono-energetic) electrons downstream at a configurable angle with respect to the z axis. --- SimCore/python/generators.py | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/SimCore/python/generators.py b/SimCore/python/generators.py index 3b7abde24f..03194bc4e7 100644 --- a/SimCore/python/generators.py +++ b/SimCore/python/generators.py @@ -359,3 +359,32 @@ def single_backwards_positron(energy: float): beam.energy = energy return beam +def single_e_wide_angle_downstream_target(minTheta = 30, maxTheta = 70, minPhi = 0, maxPhi = 360): + """A general particle source configured to shoot electrons downstream from after the target at wide angles such that they hit the side hcal. + + This generator is helpful to study the matching criteria between tracks and signals in the (side) hcal. Theta is the angle with respect to the downstream z axis. Note that the angular distribution and energy of the beam electrons can easily be modified. The particle source starts behind the target. + + Returns + ------- + gun: + configured general particle source to shoot electrons downstream the target at wide angles + """ + myGPS = gps( 'myGPS' , [ + "/gps/particle e-", + "/gps/number 1", + "/gps/pos/type Plane", + "/gps/pos/shape Rectangle", + "/gps/pos/centre 0 0 0 mm", + "/gps/pos/halfx 10 mm", + "/gps/pos/halfy 40 mm", + "/gps/ang/type iso", + f"/gps/ang/mintheta {minTheta} deg", + f"/gps/ang/maxtheta {maxTheta} deg", + f"/gps/ang/minphi {minPhi} deg", + f"/gps/ang/maxphi {maxPhi} deg", + f"/gps/ang/rot1 0 1 0", # These have been determined by trial and error so that theta is the angle from the ldmx z axis and + f"/gps/ang/rot2 1 0 0", # phi is rotates clockwise in the XY plane starting from the negative Y axis + "/gps/ene/mono 8 GeV", + ] ) + return myGPS + From 14a1445594be11a3a89664068625321e6c6402bb Mon Sep 17 00:00:00 2001 From: tillrue Date: Sun, 8 Jun 2025 11:32:22 +0200 Subject: [PATCH 2/6] Add configurable matching criteria for track + (side) hcal cluster matching pfReco.py --- Recon/python/pfReco.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/Recon/python/pfReco.py b/Recon/python/pfReco.py index 2376d00244..39505bd220 100644 --- a/Recon/python/pfReco.py +++ b/Recon/python/pfReco.py @@ -56,6 +56,11 @@ def __init__(self, name='PFlow') : self.inputTrackCollName = 'PFTracks' self.outputCollName = 'PFCandidates' self.singleParticle = False + + # Matching criteria for Track + (side) HCal cluster matching (need to be configured) + self.tkHadCaloMatchDist = 2.0 + self.tkHadCaloMinEnergyRatio = 0.3 + self.tkHadCaloMaxEnergyRatio = 2.0 self.input_ecal_passname = '' self.input_hcal_passname = '' @@ -79,4 +84,4 @@ def __init__(self, name='PFTruth') : self.sim_particles_event_passname = '' self.ecal_sp_hits_event_passname = '' self.target_sp_hits_event_passname = '' - self.target_sp_passname = '' \ No newline at end of file + self.target_sp_passname = '' From a35e0fc95e5e2499e5733cde89966fddcaf7ff88 Mon Sep 17 00:00:00 2001 From: tillrue Date: Sun, 8 Jun 2025 11:34:30 +0200 Subject: [PATCH 3/6] Add configurable matching criteria variables for track + (side) hcal cluster matching ParticleFlow.h --- Recon/include/Recon/ParticleFlow.h | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/Recon/include/Recon/ParticleFlow.h b/Recon/include/Recon/ParticleFlow.h index 225b7991fa..11bf3b3d6c 100644 --- a/Recon/include/Recon/ParticleFlow.h +++ b/Recon/include/Recon/ParticleFlow.h @@ -61,6 +61,11 @@ class ParticleFlow : public framework::Producer { std::string outputCollName_; // configuration bool singleParticle_; + + // configurable matching criteria for Track + (side) HCal cluster matching + double tkHadCaloMatchDist_; + double tkHadCaloMinEnergyRatio_; + double tkHadCaloMaxEnergyRatio_; }; } // namespace recon From 828a75d117b998a213a01134d73f827890c612d7 Mon Sep 17 00:00:00 2001 From: tillrue Date: Tue, 10 Jun 2025 11:59:18 +0200 Subject: [PATCH 4/6] Update ParticleFlow.cxx Added linking between Track and HCal cluster Added and refined arbitration for Track and HCal cluster matching I've implemented these algorithms as part of my bachelor degree project and I have written them very similar to the ones for the track + ecal cluster matching. I tried to stick to the existing variable "convention" but think that they might need some refinement at a later stage as "Calo" and "HadCalo" sound confusingly similar to describe the ECal and HCal respectively. --- Recon/src/Recon/ParticleFlow.cxx | 119 +++++++++++++++++++++++-------- 1 file changed, 89 insertions(+), 30 deletions(-) diff --git a/Recon/src/Recon/ParticleFlow.cxx b/Recon/src/Recon/ParticleFlow.cxx index 2b4b30b1cd..8c029c9c99 100644 --- a/Recon/src/Recon/ParticleFlow.cxx +++ b/Recon/src/Recon/ParticleFlow.cxx @@ -25,6 +25,9 @@ void ParticleFlow::configure(framework::config::Parameters& ps) { // Algorithm configuration singleParticle_ = ps.getParameter("singleParticle"); + tkHadCaloMatchDist_ = ps.getParameter("tkHadCaloMatchDist"); + tkHadCaloMinEnergyRatio_ = ps.getParameter("tkHadCaloMinEnergyRatio"); + tkHadCaloMaxEnergyRatio_ = ps.getParameter("tkHadCaloMaxEnergyRatio"); // Calibration factors, from jason, temperary std::vector em1{250.0, 750.0, 1250.0, 1750.0, 2250.0, 2750.0, @@ -209,21 +212,43 @@ void ParticleFlow::produce(framework::Event& event) { } } - // NOT YET IMPLEMENTED... - // tk-hadcalo linking (Side HCal) + // tk-hadcalo linking (Side HCal) BY TILLRUE std::map > tkHadCaloMap; - // for(int i=0; i > HadCaloTkMap; + std::map, float> tkHadCaloDist; + for(int i=0; i xyz = tk.getPosition(); + const std::vector pxyz = tk.getMomentum(); + const float p = sqrt(pow(pxyz[0], 2) + pow(pxyz[1], 2) + pow(pxyz[2], 2)); + + for(int j=0; j tkHadCaloMinEnergyRatio_ * p && + hcal.getEnergy() < tkHadCaloMaxEnergyRatio_ * p); // matching criteria, need to be configured + if (isMatch) { + if (tkHadCaloMap.count(i)) + tkHadCaloMap[i].push_back(j); + else + tkHadCaloMap[i] = {j}; + if (HadCaloTkMap.count(j)) + HadCaloTkMap[j].push_back(i); + else + HadCaloTkMap[j] = {i}; + } + } + } // // track / ecal cluster arbitration @@ -245,7 +270,7 @@ void ParticleFlow::produce(framework::Event& event) { } } - // track / hcal cluster arbitration + // ecal / hcal cluster arbitration std::vector EMIsHadLinked(ecalClusters.size(), false); std::vector HadIsEMLinked(hcalClusters.size(), false); std::map EMHadPairs{}; @@ -263,6 +288,26 @@ void ParticleFlow::produce(framework::Event& event) { } } + // + // track / hcal cluster arbitration BY TILLRUE + // + std::vector tkIsHadLinked(tracks.size(), false); + std::vector HadIsTkLinked(hcalClusters.size(), false); + std::map tkHadPairs{}; + for (int i = 0; i < tracks.size(); i++) { + if (tkHadCaloMap.count(i)) { + // pick first (highest-energy) unused matching cluster + for (int had_idx : tkHadCaloMap[i]) { + if (!HadIsTkLinked[had_idx]) { + HadIsTkLinked[had_idx] = true; + tkIsHadLinked[i] = true; + tkHadPairs[i] = had_idx; + break; + } + } + } + } + // can consider combining satellite clusters here... // define some "primary cluster" ID criterion // and can add fails to the primaries @@ -271,34 +316,45 @@ void ParticleFlow::produce(framework::Event& event) { // Begin building pf candidates from tracks // - // std::vector chargedMatch; - // std::vector chargedUnmatch; + // starting from tracks for (int i = 0; i < tracks.size(); i++) { ldmx::PFCandidate cand; fillCandTrack(cand, tracks[i]); // append track info to candidate - - if (!tkIsEMLinked[i]) { - // chargedUnmatch.push_back(cand); - } else { // if track is linked with ECal cluster + if (tkIsEMLinked[i]){ // if track is linked with ECal cluster fillCandEMCalo(cand, ecalClusters[tkEMPairs[i]]); - if (EMIsHadLinked[tkEMPairs[i]]) { // if ECal is linked with HCal - // cluster - fillCandHadCalo(cand, hcalClusters[EMHadPairs[tkEMPairs[i]]]); + if (EMIsHadLinked[tkEMPairs[i]]) { // if ECal is linked with HCal cluster + if (!tkIsHadLinked[i]){ // if track is NOT also linked with Hcal + fillCandHadCalo(cand, hcalClusters[EMHadPairs[tkEMPairs[i]]]);} + } + } + + if(tkIsHadLinked[i]){ // if track is linked with Hcal + fillCandHadCalo(cand, hcalClusters[tkHadPairs[i]]); + if(HadIsEMLinked[tkHadPairs[i]]){ // if hcal is linked with ecal + if(!tkIsEMLinked[i]){ // if ecal is NOT also linked with track + fillCandEMCalo(cand,ecalClusters[HadEMPairs[tkHadPairs[i]]]); + } + } else if(!tkIsEMLinked[i]){ + cand.setDistTkHcalMatch(tkHadCaloDist[{i,tkHadPairs[i]}]); + if(hcalClusters.size() > 1){ + cand.setHcalSecondEnergy(hcalClusters[tkHadPairs[i]+1].getEnergy()); + cand.setDistTkHcalSecondCluster(tkHadCaloDist[{i,tkHadPairs[i]+1}]); + } } - // chargedMatch.push_back(cand); } pfCands.push_back(cand); } // std::vector emMatch; // std::vector emUnmatch; + // ECal clusters for (int i = 0; i < ecalClusters.size(); i++) { - // already linked with ECal in the previous step - if (EMIsTkLinked[i]) continue; + if (EMIsTkLinked[i]) continue; // already linked with ECal with Track in the previous step + if (EMIsHadLinked[i] && HadIsTkLinked[EMHadPairs[i]]) continue; // already linked with track through Hcal ldmx::PFCandidate cand; fillCandEMCalo(cand, ecalClusters[i]); - if (EMIsHadLinked[tkEMPairs[i]]) { + if (EMIsHadLinked[i]) { // is Ecal is linked with Hcal fillCandHadCalo(cand, hcalClusters[EMHadPairs[i]]); // emMatch.push_back(cand); } else { @@ -306,9 +362,12 @@ void ParticleFlow::produce(framework::Event& event) { } pfCands.push_back(cand); } - std::vector hadOnly; + + // std::vector hadOnly; + // HCal clusters for (int i = 0; i < hcalClusters.size(); i++) { - if (HadIsEMLinked[i]) continue; + if (HadIsEMLinked[i]) continue; // already linked with ecal + if (HadIsTkLinked[i]) continue; // already linked with track ldmx::PFCandidate cand; fillCandHadCalo(cand, hcalClusters[i]); // hadOnly.push_back(cand); From d0d2227487edddbc98fced53f41bdb7d2575f5bb Mon Sep 17 00:00:00 2001 From: tillrue Date: Tue, 10 Jun 2025 12:18:37 +0200 Subject: [PATCH 5/6] Update ParticleFlow.cxx Added linking between Track and HCal cluster Added and refined arbitration for Track and HCal cluster matching I've implemented these algorithms as part of my bachelor degree project and I have written them very similar to the ones for the track + ecal cluster matching. I tried to stick to the existing variable "convention" but think that they might need some refinement at a later stage as "Calo" and "HadCalo" sound confusingly similar to describe the ECal and HCal respectively. --- Recon/src/Recon/ParticleFlow.cxx | 55 +++++++++++++------------------- 1 file changed, 22 insertions(+), 33 deletions(-) diff --git a/Recon/src/Recon/ParticleFlow.cxx b/Recon/src/Recon/ParticleFlow.cxx index 8c029c9c99..db9af63f19 100644 --- a/Recon/src/Recon/ParticleFlow.cxx +++ b/Recon/src/Recon/ParticleFlow.cxx @@ -10,19 +10,6 @@ void ParticleFlow::configure(framework::config::Parameters& ps) { inputHcalCollName_ = ps.getParameter("inputHcalCollName"); inputTrackCollName_ = ps.getParameter("inputTrackCollName"); outputCollName_ = ps.getParameter("outputCollName"); - - input_ecal_passname_ = ps.getParameter("input_ecal_passname"); - input_hcal_passname_ = ps.getParameter("input_hcal_passname"); - input_tracks_passname_ = - ps.getParameter("input_tracks_passname"); - - input_track_event_passname_ = - ps.getParameter("input_track_event_passname"); - input_ecal_event_passname_ = - ps.getParameter("input_ecal_event_passname"); - input_hcal_event_passname_ = - ps.getParameter("input_hcal_event_passname"); - // Algorithm configuration singleParticle_ = ps.getParameter("singleParticle"); tkHadCaloMatchDist_ = ps.getParameter("tkHadCaloMatchDist"); @@ -70,12 +57,11 @@ void ParticleFlow::fillCandEMCalo(ldmx::PFCandidate& cand, const ldmx::CaloCluster& em) { float corr = 1.; float e = em.getEnergy(); - // update energy: use min or max factor if outside calibration range if (e < eCorr_->GetX()[0]) { - corr = eCorr_->GetY()[0]; + corr = eCorr_->GetX()[0]; } else if (e > eCorr_->GetX()[eCorr_->GetN() - 1]) { - corr = eCorr_->GetY()[eCorr_->GetN() - 1]; - } else { // else look up calibration factor + corr = eCorr_->GetX()[eCorr_->GetN() - 1]; + } else { corr = eCorr_->Eval(e); } cand.setEcalEnergy(e * corr); @@ -95,9 +81,9 @@ void ParticleFlow::fillCandHadCalo(ldmx::PFCandidate& cand, float corr = 1.; float e = had.getEnergy(); if (e < hCorr_->GetX()[0]) { - corr = hCorr_->GetY()[0]; + corr = hCorr_->GetX()[0]; } else if (e > hCorr_->GetX()[hCorr_->GetN() - 1]) { - corr = hCorr_->GetY()[hCorr_->GetN() - 1]; + corr = hCorr_->GetX()[hCorr_->GetN() - 1]; } else { corr = hCorr_->Eval(e); } @@ -115,16 +101,16 @@ void ParticleFlow::fillCandHadCalo(ldmx::PFCandidate& cand, // produce track, ecal, and hcal linking void ParticleFlow::produce(framework::Event& event) { - if (!event.exists(inputTrackCollName_, input_track_event_passname_)) return; - if (!event.exists(inputEcalCollName_, input_ecal_event_passname_)) return; - if (!event.exists(inputHcalCollName_, input_hcal_event_passname_)) return; + if (!event.exists(inputTrackCollName_)) return; + if (!event.exists(inputEcalCollName_)) return; + if (!event.exists(inputHcalCollName_)) return; // get the track and clustering info - const auto ecalClusters = event.getCollection( - inputEcalCollName_, input_ecal_passname_); - const auto hcalClusters = event.getCollection( - inputHcalCollName_, input_hcal_passname_); - const auto tracks = event.getCollection( - inputTrackCollName_, input_tracks_passname_); + const auto ecalClusters = + event.getCollection(inputEcalCollName_); + const auto hcalClusters = + event.getCollection(inputHcalCollName_); + const auto tracks = + event.getCollection(inputTrackCollName_); std::vector pfCands; // multi-particle case @@ -138,7 +124,7 @@ void ParticleFlow::produce(framework::Event& event) { 4c. (Upstream?) Categorize hcal clusters as: EM/Had-like 5. Build candidates by category, moving from Tk-Ecal-Hcal */ - + std::cout << "PF Multiple Particles for Event: "<< event.getEventNumber() << std::endl; // // track-calo linking // @@ -274,6 +260,8 @@ void ParticleFlow::produce(framework::Event& event) { std::vector EMIsHadLinked(ecalClusters.size(), false); std::vector HadIsEMLinked(hcalClusters.size(), false); std::map EMHadPairs{}; + std::map HadEMPairs{}; + for (int i = 0; i < ecalClusters.size(); i++) { if (emHadCaloMap.count(i)) { // pick first (highest-energy) unused matching cluster @@ -282,6 +270,7 @@ void ParticleFlow::produce(framework::Event& event) { HadIsEMLinked[had_idx] = true; EMIsHadLinked[i] = true; EMHadPairs[i] = had_idx; + HadEMPairs[had_idx] = i; break; } } @@ -316,7 +305,7 @@ void ParticleFlow::produce(framework::Event& event) { // Begin building pf candidates from tracks // - // starting from tracks + // loop through tracks for (int i = 0; i < tracks.size(); i++) { ldmx::PFCandidate cand; fillCandTrack(cand, tracks[i]); // append track info to candidate @@ -347,7 +336,6 @@ void ParticleFlow::produce(framework::Event& event) { // std::vector emMatch; // std::vector emUnmatch; - // ECal clusters for (int i = 0; i < ecalClusters.size(); i++) { if (EMIsTkLinked[i]) continue; // already linked with ECal with Track in the previous step if (EMIsHadLinked[i] && HadIsTkLinked[EMHadPairs[i]]) continue; // already linked with track through Hcal @@ -364,7 +352,6 @@ void ParticleFlow::produce(framework::Event& event) { } // std::vector hadOnly; - // HCal clusters for (int i = 0; i < hcalClusters.size(); i++) { if (HadIsEMLinked[i]) continue; // already linked with ecal if (HadIsTkLinked[i]) continue; // already linked with track @@ -373,6 +360,7 @@ void ParticleFlow::produce(framework::Event& event) { // hadOnly.push_back(cand); pfCands.push_back(cand); } + // // track / ecal cluster arbitration // std::vector caloMatchedTks; @@ -445,6 +433,7 @@ void ParticleFlow::produce(framework::Event& event) { // } } else { + std::cout << "Matching Single Particle Signals" << std::endl; // Single-particle builder ldmx::PFCandidate pf; int pid = 0; // initialize pid to add @@ -478,4 +467,4 @@ void ParticleFlow::onProcessEnd() { } // namespace recon -DECLARE_PRODUCER(recon::ParticleFlow); +DECLARE_PRODUCER_NS(recon, ParticleFlow); From ee2cd72dc5f7a33d54d37157548f54c3cd5269e6 Mon Sep 17 00:00:00 2001 From: tillrue Date: Wed, 17 Sep 2025 11:52:58 +0200 Subject: [PATCH 6/6] Update ParticleFlow.cxx Added updates according to current version of trunk --- Recon/src/Recon/ParticleFlow.cxx | 42 +++++++++++++++++++++----------- 1 file changed, 28 insertions(+), 14 deletions(-) diff --git a/Recon/src/Recon/ParticleFlow.cxx b/Recon/src/Recon/ParticleFlow.cxx index db9af63f19..3cc247aef2 100644 --- a/Recon/src/Recon/ParticleFlow.cxx +++ b/Recon/src/Recon/ParticleFlow.cxx @@ -10,6 +10,19 @@ void ParticleFlow::configure(framework::config::Parameters& ps) { inputHcalCollName_ = ps.getParameter("inputHcalCollName"); inputTrackCollName_ = ps.getParameter("inputTrackCollName"); outputCollName_ = ps.getParameter("outputCollName"); + + input_ecal_passname_ = ps.getParameter("input_ecal_passname"); + input_hcal_passname_ = ps.getParameter("input_hcal_passname"); + input_tracks_passname_ = + ps.getParameter("input_tracks_passname"); + + input_track_event_passname_ = + ps.getParameter("input_track_event_passname"); + input_ecal_event_passname_ = + ps.getParameter("input_ecal_event_passname"); + input_hcal_event_passname_ = + ps.getParameter("input_hcal_event_passname"); + // Algorithm configuration singleParticle_ = ps.getParameter("singleParticle"); tkHadCaloMatchDist_ = ps.getParameter("tkHadCaloMatchDist"); @@ -57,10 +70,11 @@ void ParticleFlow::fillCandEMCalo(ldmx::PFCandidate& cand, const ldmx::CaloCluster& em) { float corr = 1.; float e = em.getEnergy(); + // update energy: use min or max factor if outside calibration range if (e < eCorr_->GetX()[0]) { - corr = eCorr_->GetX()[0]; + corr = eCorr_->GetY()[0]; } else if (e > eCorr_->GetX()[eCorr_->GetN() - 1]) { - corr = eCorr_->GetX()[eCorr_->GetN() - 1]; + corr = eCorr_->GetY()[eCorr_->GetN() - 1]; } else { corr = eCorr_->Eval(e); } @@ -81,9 +95,9 @@ void ParticleFlow::fillCandHadCalo(ldmx::PFCandidate& cand, float corr = 1.; float e = had.getEnergy(); if (e < hCorr_->GetX()[0]) { - corr = hCorr_->GetX()[0]; + corr = hCorr_->GetY()[0]; } else if (e > hCorr_->GetX()[hCorr_->GetN() - 1]) { - corr = hCorr_->GetX()[hCorr_->GetN() - 1]; + corr = hCorr_->GetY()[hCorr_->GetN() - 1]; } else { corr = hCorr_->Eval(e); } @@ -101,16 +115,16 @@ void ParticleFlow::fillCandHadCalo(ldmx::PFCandidate& cand, // produce track, ecal, and hcal linking void ParticleFlow::produce(framework::Event& event) { - if (!event.exists(inputTrackCollName_)) return; - if (!event.exists(inputEcalCollName_)) return; - if (!event.exists(inputHcalCollName_)) return; + if (!event.exists(inputTrackCollName_, input_track_event_passname_)) return; + if (!event.exists(inputEcalCollName_, input_ecal_event_passname_)) return; + if (!event.exists(inputHcalCollName_, input_hcal_event_passname_)) return; // get the track and clustering info - const auto ecalClusters = - event.getCollection(inputEcalCollName_); - const auto hcalClusters = - event.getCollection(inputHcalCollName_); - const auto tracks = - event.getCollection(inputTrackCollName_); + const auto ecalClusters = event.getCollection( + inputEcalCollName_, input_ecal_passname_); + const auto hcalClusters = event.getCollection( + inputHcalCollName_, input_hcal_passname_); + const auto tracks = event.getCollection( + inputTrackCollName_, input_tracks_passname_); std::vector pfCands; // multi-particle case @@ -467,4 +481,4 @@ void ParticleFlow::onProcessEnd() { } // namespace recon -DECLARE_PRODUCER_NS(recon, ParticleFlow); +DECLARE_PRODUCER(recon::ParticleFlow);