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 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 = '' diff --git a/Recon/src/Recon/ParticleFlow.cxx b/Recon/src/Recon/ParticleFlow.cxx index 2b4b30b1cd..3cc247aef2 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, @@ -72,7 +75,7 @@ void ParticleFlow::fillCandEMCalo(ldmx::PFCandidate& cand, corr = eCorr_->GetY()[0]; } else if (e > eCorr_->GetX()[eCorr_->GetN() - 1]) { corr = eCorr_->GetY()[eCorr_->GetN() - 1]; - } else { // else look up calibration factor + } else { corr = eCorr_->Eval(e); } cand.setEcalEnergy(e * corr); @@ -135,7 +138,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 // @@ -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,10 +270,12 @@ 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{}; + std::map HadEMPairs{}; + for (int i = 0; i < ecalClusters.size(); i++) { if (emHadCaloMap.count(i)) { // pick first (highest-energy) unused matching cluster @@ -257,6 +284,27 @@ void ParticleFlow::produce(framework::Event& event) { HadIsEMLinked[had_idx] = true; EMIsHadLinked[i] = true; EMHadPairs[i] = had_idx; + HadEMPairs[had_idx] = i; + break; + } + } + } + } + + // + // 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; } } @@ -271,21 +319,31 @@ void ParticleFlow::produce(framework::Event& event) { // Begin building pf candidates from tracks // - // std::vector chargedMatch; - // std::vector chargedUnmatch; + // loop through 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); } @@ -293,12 +351,12 @@ void ParticleFlow::produce(framework::Event& event) { // std::vector emMatch; // std::vector emUnmatch; 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,14 +364,17 @@ void ParticleFlow::produce(framework::Event& event) { } pfCands.push_back(cand); } - std::vector hadOnly; + + // std::vector hadOnly; 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); pfCands.push_back(cand); } + // // track / ecal cluster arbitration // std::vector caloMatchedTks; @@ -386,6 +447,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 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 +