diff --git a/ElectronID/Boost/Beam.h b/ElectronID/Boost/Beam.h new file mode 100644 index 0000000..8843381 --- /dev/null +++ b/ElectronID/Boost/Beam.h @@ -0,0 +1,74 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Wouter Deconinck + +#pragma once + +#include +// #include +// #include +#include + +using ROOT::Math::PxPyPzEVector; + +namespace eicrecon { + +// template auto find_first_with_pdg(const T* parts, const std::set& pdg) { +// T c; +// c.setSubsetCollection(); +// for (const auto& p : *parts) { +// if (pdg.count(p.getPDG()) > 0) { +// c.push_back(p); +// break; +// } +// } +// return c; +// } + +// template +// auto find_first_with_status_pdg(const T* parts, const std::set& status, +// const std::set& pdg) { +// T c; +// c.setSubsetCollection(); +// for (const auto& p : *parts) { +// if (status.count(p.getGeneratorStatus()) > 0 && pdg.count(p.getPDG()) > 0) { +// c.push_back(p); +// break; +// } +// } +// return c; +// } + +// inline auto find_first_beam_electron(const edm4hep::MCParticleCollection* mcparts) { +// return find_first_with_status_pdg(mcparts, {4}, {11}); +// } + +// inline auto find_first_beam_hadron(const edm4hep::MCParticleCollection* mcparts) { +// return find_first_with_status_pdg(mcparts, {4}, {2212, 2112}); +// } + +// inline auto find_first_scattered_electron(const edm4hep::MCParticleCollection* mcparts) { +// return find_first_with_status_pdg(mcparts, {1}, {11}); +// } + +// inline auto find_first_scattered_electron(const edm4eic::ReconstructedParticleCollection* rcparts) { +// return find_first_with_pdg(rcparts, {11}); +// } + +template +PxPyPzEVector round_beam_four_momentum(const Vector3& p_in, const float mass, + const std::vector& pz_set, + const float crossing_angle = 0.0) { + PxPyPzEVector p_out; + for (const auto& pz : pz_set) { + if (fabs(p_in.Z() / pz - 1) < 0.1) { + p_out.SetPz(pz); + break; + } + } + p_out.SetPx(p_out.Pz() * sin(crossing_angle)); + p_out.SetPz(p_out.Pz() * cos(crossing_angle)); + p_out.SetE(std::hypot(p_out.Px(), p_out.Pz(), mass)); + return p_out; +} + +} // namespace eicrecon diff --git a/ElectronID/Boost/Boost.h b/ElectronID/Boost/Boost.h new file mode 100644 index 0000000..1dc9eca --- /dev/null +++ b/ElectronID/Boost/Boost.h @@ -0,0 +1,101 @@ +// Modified from EICrecon Boost.h + +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Wouter Deconinck, Barak Schmookler + +#pragma once + +#include +#include +#include +#include +#include +#include + +using ROOT::Math::PxPyPzEVector; + +namespace eicrecon { + +using ROOT::Math::LorentzRotation; + +inline LorentzRotation determine_boost(PxPyPzEVector ei, PxPyPzEVector pi) { + + using ROOT::Math::Boost; + using ROOT::Math::RotationX; + using ROOT::Math::RotationY; + + // Step 1: Find the needed boosts and rotations from the incoming lepton and hadron beams + // (note, this will give you a perfect boost, in principle you will not know the beam momenta exactly and should use an average) + + PxPyPzEVector eo = ei; + PxPyPzEVector po = pi; + + // cout << "lab e and p: " << Form("(%f, %f, %f, %f) (%f, %f, %f, %f)", ei.Px(), ei.Py(), ei.Pz(), ei.E(), pi.Px(), pi.Py(), pi.Pz(), pi.E()) << endl; + + // Define the Boost to make beams back-to-back + const auto cmBoost = (ei + pi).BoostToCM(); + const Boost boost_to_cm(cmBoost); + + // Boost to COM frame + pi = boost_to_cm(pi); + ei = boost_to_cm(ei); + // cout << "CM boosted e and p: " << Form("(%f, %f, %f, %f) (%f, %f, %f, %f)", ei.Px(), ei.Py(), ei.Pz(), ei.E(), pi.Px(), pi.Py(), pi.Pz(), pi.E()) << endl; + + // This will boost beams from a center of momentum frame back to (nearly) their original energies + PxPyPzEVector eh(0, 0, -1*sqrt(pow(eo.E(),2)-pow(eo.M(),2)), eo.E()); + PxPyPzEVector ph(0, 0, sqrt(pow(po.E(),2)-pow(po.M(),2)), po.E()); + // cout << "headon e and p: " << Form("(%f, %f, %f, %f) (%f, %f, %f, %f)", eh.Px(), eh.Py(), eh.Pz(), eh.E(), ph.Px(), ph.Py(), ph.Pz(), ph.E()) << endl; + + const auto hoBoost = (eh + ph).BoostToCM(); + // const Boost headon_to_cm(hoBoost); + const Boost boost_to_headon(-hoBoost); + + // PxPyPzEVector et = headon_to_cm(eh); + // PxPyPzEVector pt = headon_to_cm(ph); + // cout << "headon boosted to CM e and p: " << Form("(%f, %f, %f, %f) (%f, %f, %f, %f)", et.Px(), et.Py(), et.Pz(), et.E(), pt.Px(), pt.Py(), pt.Pz(), pt.E()) << endl; + + // PxPyPzEVector erb = boost_to_headon(et); + // PxPyPzEVector prb = boost_to_headon(pt); + // cout << "reversed boost e and p: " << Form("(%f, %f, %f, %f) (%f, %f, %f, %f)", erb.Px(), erb.Py(), erb.Pz(), erb.E(), prb.Px(), prb.Py(), prb.Pz(), prb.E()) << endl; + + // Boost and rotate the incoming beams to find the proper rotations TLorentzVector + + // Rotate to head-on + RotationY rotAboutY(-1.0 * atan2(pi.Px(), pi.Pz())); // Rotate to remove x component of beams + RotationX rotAboutX(+1.0 * atan2(pi.Py(), pi.Pz())); // Rotate to remove y component of beams + + // PxPyPzEVector er = rotAboutX(rotAboutY(ei)); + // PxPyPzEVector pr = rotAboutX(rotAboutY(pi)); + // cout << "rotation e and p: " << Form("(%f, %f, %f, %f) (%f, %f, %f, %f)", er.Px(), er.Py(), er.Pz(), er.E(), pr.Px(), pr.Py(), pr.Pz(), pr.E()) << endl; + + // PxPyPzEVector ef = boost_to_headon(er); + // PxPyPzEVector pf = boost_to_headon(pr); + // cout << "final boost e and p: " << Form("(%f, %f, %f, %f) (%f, %f, %f, %f)", ef.Px(), ef.Py(), ef.Pz(), ef.E(), pf.Px(), pf.Py(), pf.Pz(), pf.E()) << endl; + + // cout << "**" << endl; + + // final matrix: P' = [BtoH][RX][RY][BtoCM]P <-- Matrix multi. goes from R to L + LorentzRotation tf; + tf *= boost_to_headon; + tf *= rotAboutX; + tf *= rotAboutY; + tf *= boost_to_cm; + + // PxPyPzEVector em = tf(eo); + // PxPyPzEVector pm = tf(po); + // cout << "boost matrix e and p: " << Form("(%f, %f, %f, %f) (%f, %f, %f, %f)", em.Px(), em.Py(), em.Pz(), em.E(), pm.Px(), pm.Py(), pm.Pz(), pm.E()) << endl; + + return tf; +} + +inline PxPyPzEVector apply_boost(const LorentzRotation& tf, PxPyPzEVector part) { + + // Step 2: Apply boosts and rotations to any particle 4-vector + // (here too, choices will have to be made as to what the 4-vector is for reconstructed particles) + + // Boost and rotate particle 4-momenta into the headon frame + tf(part); + return part; +} + +} // namespace eicrecon diff --git a/ElectronID/Boost/getBoost.h b/ElectronID/Boost/getBoost.h new file mode 100644 index 0000000..f08665a --- /dev/null +++ b/ElectronID/Boost/getBoost.h @@ -0,0 +1,47 @@ +#pragma once + +#include "Beam.h" +#include "Boost.h" + +#include +using ROOT::Math::PxPyPzEVector; + +#include +using ROOT::Math::LorentzRotation; + +LorentzRotation getBoost(double eE, double eN, double mE, double mN) { + + TVector3 ve(0,0,-sqrt(eE*eE-mE*mE)); + TVector3 vn(0,0,sqrt(eN*eN-mN*mN)); + + const PxPyPzEVector ei( + eicrecon::round_beam_four_momentum( + ve, + mE, + {-1*eE}, + 0.0) + ); + + const PxPyPzEVector ni( + eicrecon::round_beam_four_momentum( + vn, + mN, + {eN}, + -0.025) + ); + + // std::cout << "initial electron beam: " << Form("(%f, %f, %f, %f)", ei.Px(), ei.Py(), ei.Pz(), ei.E()) << std::endl; + // std::cout << "initial nucleon beam: " << Form("(%f, %f, %f, %f)", ni.Px(), ni.Py(), ni.Pz(), ni.E()) << std::endl; + // std::cout << "Created vector, e, p (mrad) " << ei.Theta() * 1000 << " " << ni.Theta() * 1000 << std::endl; + + LorentzRotation boost = eicrecon::determine_boost(ei, ni); + + PxPyPzEVector boosted_e = boost(ei); + PxPyPzEVector boosted_n = boost(ni); + + // std::cout << "Boosted electron beam: " << Form("(%f, %f, %f, %f)", boosted_e.Px(), boosted_e.Py(), boosted_e.Pz(), boosted_e.E()) << std::endl; + // std::cout << "Boosted nucleon beam: " << Form("(%f, %f, %f, %f)", boosted_n.Px(), boosted_n.Py(), boosted_n.Pz(), boosted_n.E()) << std::endl; + // std::cout << "Incoming headon, e, p (mrad) " << boosted_e.Theta() * 1000 << " " << boosted_n.Theta() * 1000 << std::endl; + + return boost; +} diff --git a/ElectronID/ElectronID.cc b/ElectronID/ElectronID.cc index 2b50165..508a21c 100644 --- a/ElectronID/ElectronID.cc +++ b/ElectronID/ElectronID.cc @@ -1,27 +1,31 @@ #include "ElectronID.hh" #include "edm4hep/utils/vector_utils.h" -#include "edm4hep/utils/kinematics.h" #include "edm4eic/ClusterCollection.h" #include +#include +using ROOT::Math::PxPyPzEVector; + ElectronID::ElectronID() { mEe = 10.; mEh = 100.; std::cout << "!!! ElectronID: You have not specified beam energies...defaulting to 10x100 GeV !!!" << std::endl; - mEoP_min = 0.9; + mEoP_min = 0.8; mEoP_max = 1.2; mDeltaH_min = 0.*mEe; mDeltaH_min = 5.*mEe; - mIsoR = 0.4; + mIsoR = 0.7; mIsoE = 0.9; + minTrackPoints = 3; + boost = LorentzRotation(); // Initialize to identity } ElectronID::ElectronID(double Ee, double Eh) { @@ -29,120 +33,287 @@ ElectronID::ElectronID(double Ee, double Eh) { mEe = Ee; mEh = Eh; - mEoP_min = 0.9; + mEoP_min = 0.8; mEoP_max = 1.2; mDeltaH_min = 0.*mEe; mDeltaH_min = 5.*mEe; - mIsoR = 0.4; + // mIsoR = 2; // tested for theta > 150 + mIsoR = 0.7; mIsoE = 0.9; + minTrackPoints = 3; + + boost = LorentzRotation(); // Initialize to identity } ElectronID::~ElectronID() { } - void ElectronID::SetEvent(const podio::Frame* event) { + // std::cout << "** Setting event in ElectronID... " << std::endl; mEvent = event; + eScatIndex = -1; + hfs_dpt.clear(); + hfs_dpz.clear(); + hfs_de.clear(); + hfs_theta.clear(); + e_det.clear(); + jet_e_det.clear(); + pi_det.clear(); + else_det.clear(); + det_val.clear(); + // std::cout << "** Event set in ElectronID. " << std::endl; + return; +} + +edm4hep::MCParticle ElectronID::GetMC(edm4eic::ReconstructedParticle e_rec) { + + // std::cout << "Getting MC particle associated with reconstructed particle..." << std::endl; + + const auto& RecoMC = static_cast(*(mEvent->get("ReconstructedParticleAssociations"))); + for(const auto& assoc : RecoMC) { + if(assoc.getRec() == e_rec) + return assoc.getSim(); + } + return edm4hep::MCParticle(); } +int ElectronID::Check_eID(edm4eic::ReconstructedParticle e_rec) { + + edm4hep::MCParticleCollection meMC = GetMCElectron(); + if ( meMC.size() == 0 ) + return 86; // No MC electron found + + const auto& RecoMC = static_cast(*(mEvent->get("ReconstructedParticleAssociations"))); + for(const auto& assoc : RecoMC) { + if(assoc.getRec() == e_rec) + { + if (assoc.getSim() == meMC[0]) + return 0; + else + return assoc.getSim().getPDG(); + } + } + + return 86; +} + +void ElectronID::CheckClusters() { + + const auto& EcalEndcapNClusters = static_cast(*(mEvent->get("EcalEndcapNClusters"))); + const auto& EcalBarrelScFiClusters = static_cast(*(mEvent->get("EcalBarrelScFiClusters"))); + const auto& EcalEndcapPClusters = static_cast(*(mEvent->get("EcalEndcapPClusters"))); + + std::cout << " Number of clusters in EcalEndcapN: " << EcalEndcapNClusters.size() << std::endl; + std::cout << " Number of clusters in EcalBarrelScFi: " << EcalBarrelScFiClusters.size() << std::endl; + std::cout << " Number of clusters in EcalEndcapP: " << EcalEndcapPClusters.size() << std::endl; + + return; +} + +edm4eic::ReconstructedParticleCollection ElectronID::FindHadronicFinalState(int object_id) { + + // edm4eic::HadronicFinalStateCollection meRecon; + edm4eic::ReconstructedParticleCollection meRecon; + meRecon.setSubsetCollection(); + + // auto& rcparts = mEvent->get("HadronicFinalState"); + const auto& rcparts = static_cast(*(mEvent->get("ReconstructedParticles"))); + + for(const auto& mcp : rcparts) { + if ( mcp.getObjectID().index != object_id ) + meRecon.push_back(mcp); + } + + return meRecon; +} edm4eic::ReconstructedParticleCollection ElectronID::FindScatteredElectron() { + // std::cout << "\nFinding scattered electron candidates..." << std::endl; + // Get all the edm4eic objects needed for electron ID - auto& rcparts = mEvent->get("ReconstructedParticles"); + const auto& rcparts = static_cast(*(mEvent->get("ReconstructedParticles"))); // Create collection for storing scattered electron candidates // (subset collection of ReconstructedParticleCollection) edm4eic::ReconstructedParticleCollection scatteredElectronCandidates; - scatteredElectronCandidates->setSubsetCollection(); + scatteredElectronCandidates.setSubsetCollection(); + + edm4eic::ReconstructedParticleCollection scatteredElectronCandidates_noEoP; + scatteredElectronCandidates_noEoP.setSubsetCollection(); // Loop over all ReconstructedParticles for this event for (const auto& reconPart : rcparts) { - // Require negative particle - if(reconPart.getCharge() >= 0) continue; - // Require at least one track and one cluster - if(reconPart.getClusters().size() == 0 || reconPart.getTracks().size() == 0) continue; + if(reconPart.getClusters().size() > 0 && reconPart.getTracks().size() > 0) + { + // Require negative particle + if(reconPart.getCharge() >= 0) + continue; - // Calculate rcpart_ member variables for this event - CalculateParticleValues(reconPart, rcparts); + int n_track_points = reconPart.getTracks()[0].measurements_size(); - // Calculate E/p and isolation fraction for this event - // Note that the rcpart_ variables are set in CalculateParticleValues - double recon_EoP = rcpart_sum_cluster_E / edm4hep::utils::magnitude(reconPart.getMomentum()); - double recon_isoE = rcpart_sum_cluster_E / rcpart_isolation_E; - - // Apply scattered electron ID cuts - if(recon_EoP < mEoP_min || recon_EoP > mEoP_max) continue; - if(recon_isoE < mIsoE) continue; + // Calculate rcpart_ member variables for this event + CalculateParticleValues(reconPart, rcparts); - // If particle passes cuts, add to output collection - scatteredElectronCandidates.push_back(reconPart); + // Calculate isolation fraction for this event + double recon_isoE = rcpart_sum_cluster_E / rcpart_isolation_E; - } + // Calculate E/p for this event + double recon_EoP = rcpart_sum_cluster_E / edm4hep::utils::magnitude(reconPart.getMomentum()); + + det_val.push_back({Check_eID(reconPart), n_track_points, recon_EoP, recon_isoE}); + + if ( n_track_points < minTrackPoints ) + continue; - return scatteredElectronCandidates; + if(recon_isoE < mIsoE) + continue; + double trackTheta = edm4hep::utils::anglePolar(reconPart.getMomentum())*(180./M_PI); + double clusterTheta = GetMomentumVectorFromCluster(reconPart, 0).Theta()*(180./M_PI); + + if ( recon_EoP > mEoP_min && recon_EoP < mEoP_max ) + { + scatteredElectronCandidates.push_back(reconPart); + continue; + } + else if ( (trackTheta > 158 && trackTheta < 162) || (clusterTheta > 22 && clusterTheta < 33) ) + { + scatteredElectronCandidates_noEoP.push_back(reconPart); + continue; + } + } + } + + // If EoP is found use that, otherwise loosen cuts + if (scatteredElectronCandidates.size() > 0) + return scatteredElectronCandidates; + else + return scatteredElectronCandidates_noEoP; } -edm4hep::MCParticleCollection ElectronID::GetMCElectron() { +edm4hep::MCParticleCollection ElectronID::GetMCHadronicFinalState() { - edm4hep::MCParticleCollection meMC; - meMC->setSubsetCollection(); + edm4hep::MCParticleCollection mhMC; + mhMC.setSubsetCollection(); - auto& mcparts = mEvent->get("MCParticles"); + const auto& mcparts = static_cast(*(mEvent->get("MCParticles"))); - std::vector mc_electrons; - + std::vector mc_hadronic; + edm4hep::MCParticleCollection meMC = GetMCElectron(); + + bool found_scattered_e = false; for(const auto& mcp : mcparts) { - if(mcp.getPDG() == 11 && mcp.getGeneratorStatus() == 1) { - mc_electrons.push_back(mcp); + if (mcp.getGeneratorStatus() == 1) + { + if ( meMC.size() == 0 ) + mhMC.push_back(mcp); + else if (mcp.getObjectID().index != meMC[0].getObjectID().index ) + mhMC.push_back(mcp); } } - // For now, just take first electron - if(mc_electrons.size() > 0) { - meMC.push_back(mc_electrons[0]); + return mhMC; +} + +edm4hep::MCParticleCollection ElectronID::GetMCElectron() { + + edm4hep::MCParticleCollection meMC; + meMC.setSubsetCollection(); + + const auto& mcparts = static_cast(*(mEvent->get("MCParticles"))); + if ( eScatIndex != -1 ) + meMC.push_back(mcparts[eScatIndex]); + + //// + + for (const auto& mcp : mcparts) + { + if (mcp.getPDG() != 11 || mcp.getGeneratorStatus() != 4) + continue; + + std::vector stack; + stack.insert(stack.end(), mcp.getDaughters().begin(), mcp.getDaughters().end()); + + int shortest_gen = 999; + int generations = 0; + edm4hep::MCParticle meMC_candidates; + bool found_candidate = false; + + while (!stack.empty() ) { + generations++; + auto cur = stack.back(); + stack.pop_back(); + + if (cur.getPDG() == 11 && cur.getGeneratorStatus() == 1) { + + if ( generations < shortest_gen ) + { + shortest_gen = generations; + meMC_candidates = cur; + found_candidate = true; + } + break; + } + + const auto& kids = cur.getDaughters(); + if (!kids.empty()) { + stack.insert(stack.end(), kids.begin(), kids.end()); + } + } + + if ( found_candidate ) + meMC.push_back(meMC_candidates); } return meMC; - } edm4eic::ReconstructedParticleCollection ElectronID::GetTruthReconElectron() { - edm4hep::MCParticleCollection meMC = GetMCElectron(); + // cout << "New process " << endl; + + const edm4hep::MCParticleCollection meMC = GetMCElectron(); edm4eic::ReconstructedParticleCollection meRecon; - meRecon->setSubsetCollection(); + meRecon.setSubsetCollection(); - auto& RecoMC = mEvent->get("ReconstructedParticleAssociations"); + if ( meMC.size() == 0 ) + return meRecon; // No MC electron found + + const auto& RecoMC = static_cast(*(mEvent->get("ReconstructedParticleAssociations"))); + + for(const auto& assoc : RecoMC) + { + auto e_candidat = assoc.getSim(); - for(const auto& assoc : RecoMC) { if(assoc.getSim() == meMC[0]) { meRecon.push_back(assoc.getRec()); - break; + // break; } } - return meRecon; + const auto& rcparts = static_cast(*(mEvent->get("ReconstructedParticles"))); + edm4hep::Vector3f lead_pos(meMC[0].getEndpoint().x, meMC[0].getEndpoint().y, meMC[0].getEndpoint().z); + CheckSurroundingClusters(lead_pos, rcparts); + // std::cout << meRecon.size() << " reconstructed particles associated with the MC electron." << std::endl; + + return meRecon; } - - void ElectronID::CalculateParticleValues(const edm4eic::ReconstructedParticle& rcp, const edm4eic::ReconstructedParticleCollection& rcparts) { rcpart_sum_cluster_E = 0.; rcpart_lead_cluster_E = 0.; rcpart_isolation_E = 0.; - rcpart_deltaH = 0.; const edm4eic::Cluster* lead_cluster = nullptr; @@ -156,13 +327,21 @@ void ElectronID::CalculateParticleValues(const edm4eic::ReconstructedParticle& r if(!lead_cluster) return; - const auto& lead_pos = lead_cluster->getPosition(); + CheckSurroundingClusters(lead_cluster->getPosition(), rcparts); + + return; +} + +void ElectronID::CheckSurroundingClusters(const edm4hep::Vector3f& lead_pos, + const edm4eic::ReconstructedParticleCollection& rcparts) { + + rcpart_isolation_E = 0.; + rcpart_n_clusters = 0; + double lead_eta = edm4hep::utils::eta(lead_pos); double lead_phi = edm4hep::utils::angleAzimuthal(lead_pos); for (const auto& other_rcp : rcparts) { - if (&other_rcp == &rcp) continue; // Skip the same particle - for (const auto& other_cluster : other_rcp.getClusters()) { const auto& other_pos = other_cluster.getPosition(); @@ -181,6 +360,7 @@ void ElectronID::CalculateParticleValues(const edm4eic::ReconstructedParticle& r // Check if the cluster is within the isolation cone if (dR < mIsoR) { rcpart_isolation_E += other_cluster.getEnergy(); + rcpart_n_clusters ++; } } } @@ -188,13 +368,54 @@ void ElectronID::CalculateParticleValues(const edm4eic::ReconstructedParticle& r return; } +void ElectronID::GetEminusPzSum(double &TrackEminusPzSum, double &CalEminusPzSum) { + + const auto& rcparts = static_cast(*(mEvent->get("ReconstructedParticles"))); + + for (const auto& reconPart : rcparts) { + + // Require at least one track and one cluster + // if(reconPart.getClusters().size() == 0 || reconPart.getTracks().size() == 0) continue; + + if(reconPart.getTracks().size() > 0 ) + { + int n_track_points = reconPart.getTracks()[0].measurements_size(); + if ( n_track_points < minTrackPoints ) continue; + + PxPyPzEVector vT(reconPart.getMomentum().x, reconPart.getMomentum().y, reconPart.getMomentum().z, reconPart.getEnergy()); + vT = boost(vT); + TrackEminusPzSum += (vT.E() - vT.Pz()); + + if ( reconPart.getClusters().size() > 0 ) + { + PxPyPzEVector vC(reconPart.getMomentum().x, reconPart.getMomentum().y, reconPart.getMomentum().z, GetCalorimeterEnergy(reconPart)); + // PxPyPzEVector vC = GetMomentumVectorFromCluster(reconPart, reconPart.getMass()); + // if ( reconPart.getTracks().size() > 0 && n_track_points >= minTrackPoints ) + // vC.SetPxPyPzE(reconPart.getMomentum().x, reconPart.getMomentum().y, reconPart.getMomentum().z, GetCalorimeterEnergy(reconPart)); + vC = boost(vC); + CalEminusPzSum += (vC.E() - vC.Pz()); + } + } + else if (reconPart.getClusters().size() > 0 ) + { + PxPyPzEVector vC = GetMomentumVectorFromCluster(reconPart, 0); + vC = boost(vC); + CalEminusPzSum += (vC.E() - vC.E()*std::cos(vC.Theta())); + } + } + + // std::cout << " recon E - Pz sum: " << reconEminusPzSum << std::endl; + + return; +} + edm4eic::ReconstructedParticle ElectronID::SelectHighestPT(const edm4eic::ReconstructedParticleCollection& ecandidates) { edm4eic::ReconstructedParticle erec; double max_pT = 0.; - for(auto ecand : ecandidates) { - double e_pT = edm4hep::utils::magnitudeTransverse(ecand.getMomentum()); + for(const auto& ecand : ecandidates) { + double e_pT = ecand.getTracks().size() > 0 ? edm4hep::utils::magnitudeTransverse(ecand.getMomentum()) : GetMomentumVectorFromCluster(ecand, 0.000511).Pt(); if(e_pT > max_pT) { erec = ecand; max_pT = e_pT; @@ -202,7 +423,45 @@ edm4eic::ReconstructedParticle ElectronID::SelectHighestPT(const edm4eic::Recons } return erec; +} +double ElectronID::GetClusterCone(const edm4eic::ReconstructedParticle& rcp, double frac=1) +{ + const edm4eic::Cluster* lead_cluster = nullptr; + double sum_cluster_E = 0.; + double lead_cluster_E = 0.; + + for (const auto& cluster : rcp.getClusters()) { + sum_cluster_E += cluster.getEnergy(); + if(cluster.getEnergy() > lead_cluster_E) { + lead_cluster = &cluster; + lead_cluster_E = cluster.getEnergy(); + } + } + + if (!lead_cluster) + return -1.; + + std::vector> dr_de; + + for (const auto& cluster : rcp.getClusters()) + { + if (lead_cluster == &cluster) + continue; + + double d_eta = edm4hep::utils::eta(cluster.getPosition()) - edm4hep::utils::eta(lead_cluster->getPosition()); + double d_phi = edm4hep::utils::angleAzimuthal(cluster.getPosition()) - edm4hep::utils::angleAzimuthal(lead_cluster->getPosition()); + + if (d_phi > M_PI) d_phi -= 2 * M_PI; + if (d_phi < -M_PI) d_phi += 2 * M_PI; + + double dR = std::sqrt(std::pow(d_eta, 2) + std::pow(d_phi, 2)); + dr_de.emplace_back(dR, cluster.getEnergy()); + } + + std::sort(dr_de.begin(), dr_de.end(), [](const auto& left, const auto& right) {return left.first < right.first;}); + + return dr_de.back().first; // Return the dR of the farthest cluster } double ElectronID::GetCalorimeterEnergy(const edm4eic::ReconstructedParticle& rcp) { @@ -215,4 +474,46 @@ double ElectronID::GetCalorimeterEnergy(const edm4eic::ReconstructedParticle& rc } +double ElectronID::GetClusterTheta(const edm4eic::ReconstructedParticle& rcp) { + + const edm4eic::Cluster* lead_cluster = nullptr; + double lead_E = 0.; + + for (const auto& cluster : rcp.getClusters()) { + if(cluster.getEnergy() > lead_E) { + lead_cluster = &cluster; + lead_E = cluster.getEnergy(); + } + } + + if(!lead_cluster) return -999.; + + return lead_cluster->getIntrinsicTheta(); +} + +PxPyPzEVector ElectronID::GetMomentumVectorFromCluster(const edm4eic::ReconstructedParticle& rcp, double mass) { + const edm4eic::Cluster* lead_cluster = nullptr; + double sum_cluster_E = 0.; + double lead_E = 0.; + + for (const auto& cluster : rcp.getClusters()) { + sum_cluster_E += cluster.getEnergy(); + if(cluster.getEnergy() > lead_E) { + lead_cluster = &cluster; + lead_E = cluster.getEnergy(); + } + } + + if(!lead_cluster) + cout << "Warning: No clusters found for this reconstructed particle!" << endl; + + if(!lead_cluster) return PxPyPzEVector(0, 0, 0, 0); + + double p = std::sqrt(std::pow(sum_cluster_E, 2) - std::pow(mass, 2)); + double px = p * (lead_cluster->getPosition().x / edm4hep::utils::magnitude(lead_cluster->getPosition())); + double py = p * (lead_cluster->getPosition().y / edm4hep::utils::magnitude(lead_cluster->getPosition())); + double pz = p * (lead_cluster->getPosition().z / edm4hep::utils::magnitude(lead_cluster->getPosition())); + + return PxPyPzEVector(px, py, pz, sum_cluster_E); +} diff --git a/ElectronID/ElectronID.hh b/ElectronID/ElectronID.hh index ad5e047..0633c11 100644 --- a/ElectronID/ElectronID.hh +++ b/ElectronID/ElectronID.hh @@ -4,7 +4,18 @@ #include "podio/Frame.h" #include "edm4eic/ReconstructedParticleCollection.h" +#include "edm4eic/MCRecoParticleAssociationCollection.h" #include "edm4hep/MCParticleCollection.h" +#include "edm4hep/utils/vector_utils.h" +#include "edm4eic/ClusterCollection.h" + +#include +using ROOT::Math::LorentzRotation; + +#include +using ROOT::Math::PxPyPzEVector; + +#include class ElectronID{ @@ -18,14 +29,53 @@ public: inline void SetEoPMin(double eopmin) {mEoP_min = eopmin;} inline void SetDeltaHMin(double deltahmin) {mDeltaH_min = deltahmin;} inline void SetIsolation(double isor, double isoe) {mIsoR = isor; mIsoE = isoe;} + inline void SetMinTrackPoints(int minPoints) { minTrackPoints = minPoints; } void SetEvent(const podio::Frame* event); + void SetBoost(LorentzRotation fboost) { boost = fboost; } + int Check_eID(edm4eic::ReconstructedParticle e_rec); + edm4hep::MCParticle GetMC(edm4eic::ReconstructedParticle e_rec); + edm4eic::ReconstructedParticleCollection FindHadronicFinalState(int object_id); edm4eic::ReconstructedParticleCollection FindScatteredElectron(); edm4eic::ReconstructedParticleCollection GetTruthReconElectron(); edm4hep::MCParticleCollection GetMCElectron(); + edm4hep::MCParticleCollection GetMCHadronicFinalState(); edm4eic::ReconstructedParticle SelectHighestPT(const edm4eic::ReconstructedParticleCollection& rcparts); double GetCalorimeterEnergy(const edm4eic::ReconstructedParticle& rcp); + double GetClusterCone(const edm4eic::ReconstructedParticle& rcp, double frac=1); + PxPyPzEVector GetMomentumVectorFromCluster(const edm4eic::ReconstructedParticle& rcp, double mass); + double GetClusterTheta(const edm4eic::ReconstructedParticle& rcp); + void GetEminusPzSum(double &TrackEminusPzSum, double &CalEminusPzSum); + void CheckClusters(); + + double get_mEoP_min() const { return mEoP_min; } + double get_mEoP_max() const { return mEoP_max; } + double get_mDeltaH_min() const { return mDeltaH_min; } + double get_mDeltaH_max() const { return mDeltaH_max; } + double get_mIsoR() const { return mIsoR; } + double get_mIsoE() const { return mIsoE; } + int GetMinTrackPoints() const { return minTrackPoints; } + + // for HFS QA + vector hfs_dpt; + vector hfs_dpz; + vector hfs_de; + vector hfs_theta; + + struct DetValues { + int parType; // 0 for dis electron, rest follow pdg code + int nTrackPoints; + double recon_EoP; + double recon_isoE; + }; + vector det_val; + vector e_det; + vector jet_e_det; + vector pi_det; + vector else_det; + + double rcpart_n_clusters; private: @@ -33,6 +83,7 @@ private: double mEe; double mEh; + LorentzRotation boost; double mEoP_min; double mEoP_max; @@ -40,16 +91,18 @@ private: double mDeltaH_max; double mIsoR; double mIsoE; + int minTrackPoints = 3; void CalculateParticleValues(const edm4eic::ReconstructedParticle& rcp, const edm4eic::ReconstructedParticleCollection& rcparts); + void CheckSurroundingClusters(const edm4hep::Vector3f& lead_pos, + const edm4eic::ReconstructedParticleCollection& rcparts); double rcpart_sum_cluster_E; double rcpart_lead_cluster_E; double rcpart_isolation_E; - double rcpart_deltaH; - + int eScatIndex; }; #endif diff --git a/ElectronID/InclusiveSkim.C b/ElectronID/InclusiveSkim.C index 6788d41..e99bf0f 100644 --- a/ElectronID/InclusiveSkim.C +++ b/ElectronID/InclusiveSkim.C @@ -1,12 +1,8 @@ +#include "preLoadLib.hh" + // Data model headers -#include "edm4eic/ReconstructedParticleCollection.h" -#include "edm4hep/MCParticleCollection.h" -#include "edm4hep/utils/vector_utils.h" -#include "edm4hep/utils/kinematics.h" -#include "edm4eic/ClusterCollection.h" -#include "edm4eic/MCRecoParticleAssociationCollection.h" #include "podio/Frame.h" -#include "podio/ROOTFrameReader.h" +#include "podio/ROOTReader.h" // ROOT headers #include "TTree.h" @@ -18,22 +14,32 @@ // Analysis headers #include "InclusiveSkim.h" #include "ElectronID.cc" +#include "Boost/getBoost.h" void InclusiveSkim() { double Ee = 10.; double Eh = 100.; - vector inFiles = {"pythia8NCDIS_10x100_minQ2=10_beamEffects_xAngle=-0.025_hiDiv_1.0001.eicrecon.tree.edm4eic.root"}; + // access local file + // vector inFiles = {"pythia8NCDIS_10x100_minQ2=10_beamEffects_xAngle=-0.025_hiDiv_1.0001.eicrecon.tree.edm4eic.root"}; - auto reader = podio::ROOTFrameReader(); - reader.openFiles(inFiles); + // access remote file + vector inFiles = {"root://dtn-eic.jlab.org:1094//volatile/eic/EPIC//RECO/26.02.0/epic_craterlake/DIS/NC/18x275/minQ2=10/pythia8NCDIS_18x275_minQ2=10_beamEffects_xAngle=-0.025_hiDiv_1.0000.eicrecon.edm4eic.root"}; - ElectronID* eFinder = new ElectronID(Ee, Eh); + auto reader = podio::ROOTReader(); + reader.openFiles(inFiles); TString outFileName = Form("inclusive_skim_%.0fx%.0fGeV.root", Ee, Eh); CreateOutputTree(outFileName); + ElectronID* eFinder = new ElectronID(Ee, Eh); + eFinder->SetEoPMin(0.8); // set E/p cut + eFinder->SetIsolation(0.7, 0.9); // set isolation parameters (cone size, energy fraction) + eFinder->SetMinTrackPoints(4); // set number of minimum points required for a track + + LorentzRotation boost = getBoost( Ee, Eh, MASS_ELECTRON, MASS_PROTON); // if you have your own boost calculation, you can replace this line. + eFinder->SetBoost(boost); for(size_t ev = 0; ev < reader.getEntries("events"); ev++) { @@ -188,4 +194,4 @@ TLorentzVector GetHadronBeam(double fEh) { hadron_beam.SetE(std::hypot(fEh, Mp)); return hadron_beam; -} +} \ No newline at end of file diff --git a/ElectronID/InclusiveSkim.h b/ElectronID/InclusiveSkim.h index 928b6ff..bdf3177 100644 --- a/ElectronID/InclusiveSkim.h +++ b/ElectronID/InclusiveSkim.h @@ -11,6 +11,9 @@ void ResetVariables(); void CalculateElectronKinematics(double fEe, double fEh, TLorentzVector kf, float& xB, float& Q2, float& W, float& y, float& nu); TLorentzVector GetHadronBeam(double fEh); +#define MASS_ELECTRON 0.000511 +#define MASS_PROTON 0.93827 + using namespace std; int positive_eID; @@ -40,4 +43,3 @@ float e_clust_Q2; float e_clust_W; float e_clust_y; float e_clust_nu; - diff --git a/ElectronID/README b/ElectronID/README new file mode 100644 index 0000000..b0c7f76 --- /dev/null +++ b/ElectronID/README @@ -0,0 +1,10 @@ +DIS Electron Finder: +* Example codes assume the use of eic-shell and PODIO +* ElectronID.cc and ElectronID.hh provides a list of functions for DIS electron analysis. Use the function FindScatteredElectron() for the standard electron finding method. +* InclusiveSkim.C is an example on how to use the ElectronFinder. This can be ran with the most updated eic-shell for simulation campaign after 25.05. + +The current algorithm apply cuts as follows to select the DIS electron: +* Particle has associated track with negative curvature (negative charge). +* Track with momentum $p$ can be matched to ECAL cluster with compatible energy $E_\text{clus}$ such that $0.9 < E_\text{clus}/p < 1.2$. +* Isolation of the matched ECAL cluster, such that the matched cluster constitutes more than $90~\%$ of the energy contained within an $\eta-\phi$ cone centered on the matched cluster: $E_\text{clus} / \sum (E_{\Delta R < 0.4}) > 0.9$, with $\Delta R = \sqrt{\Delta \eta^2 - \Delta \phi^2}$. +* If more than one candidate remains after this selection, the highest pT candidate is chosen as the DIS electron. \ No newline at end of file diff --git a/ElectronID/preLoadLib.hh b/ElectronID/preLoadLib.hh new file mode 100644 index 0000000..186f421 --- /dev/null +++ b/ElectronID/preLoadLib.hh @@ -0,0 +1,8 @@ +R__LOAD_LIBRARY(libpodio.so) +R__LOAD_LIBRARY(libpodioRootIO.so) +R__LOAD_LIBRARY(/opt/local/lib/libedm4hep.so) +R__LOAD_LIBRARY(/opt/local/lib/libedm4hepDict.so) +R__LOAD_LIBRARY(/opt/local/lib/libedm4hepUtils.so) +R__LOAD_LIBRARY(/opt/local/lib/libedm4hepRDF.so) +R__LOAD_LIBRARY(/opt/local/lib/libedm4eic.so) +R__LOAD_LIBRARY(/opt/local/lib/libedm4eicDict.so) \ No newline at end of file