BPHNano/plugins/helper.h

0001 // original author: RK18 team
0002
0003 #ifndef PhysicsTools_BPHNano_helpers
0004 #define PhysicsTools_BPHNano_helpers
0005
0006 #include <algorithm>
0007 #include <limits>
0008 #include <memory>
0009 #include <vector>
0010
0011 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0012 #include "DataFormats/Candidate/interface/Candidate.h"
0013 #include "DataFormats/GeometryCommonDetAlgo/interface/GlobalError.h"
0014 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
0015 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0016 #include "DataFormats/GeometryVector/interface/PV3DBase.h"
0017 #include "DataFormats/Math/interface/deltaR.h"
0018 #include "DataFormats/PatCandidates/interface/CompositeCandidate.h"
0019 #include "Math/LorentzVector.h"
0020 #include "RecoVertex/KinematicFitPrimitives/interface/RefCountedKinematicTree.h"
0021 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
0022 #include "TVector3.h"
0023 #include "TrackingTools/GeomPropagators/interface/AnalyticalImpactPointExtrapolator.h"
0024 #include "TrackingTools/IPTools/interface/IPTools.h"
0025 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0026
0027 typedef std::vector<reco::TransientTrack> TransientTrackCollection;
0028
0029 namespace bph {
0030
0031   constexpr float PROT_MASS = 0.938272;
0032   constexpr float K_MASS = 0.493677;
0033   constexpr float PI_MASS = 0.139571;
0034   constexpr float LEP_SIGMA = 0.0000001;
0035   constexpr float K_SIGMA = 0.000016;
0036   constexpr float PI_SIGMA = 0.000016;
0037   constexpr float PROT_SIGMA = 0.000016;
0038   constexpr float MUON_MASS = 0.10565837;
0039   constexpr float ELECTRON_MASS = 0.000511;
0040
0041   inline std::pair<float, float> min_max_dr(const std::vector<edm::Ptr<reco::Candidate>>& cands) {
0042     float min_dr = std::numeric_limits<float>::max();
0043     float max_dr = 0.;
0044     for (size_t i = 0; i < cands.size(); ++i) {
0045       for (size_t j = i + 1; j < cands.size(); ++j) {
0046         float dr = reco::deltaR(*cands.at(i), *cands.at(j));
0047         min_dr = std::min(min_dr, dr);
0048         max_dr = std::max(max_dr, dr);
0049       }
0050     }
0051     return std::make_pair(min_dr, max_dr);
0052   }
0053
0054   template <typename FITTER, typename LORENTZ_VEC>
0055   inline double cos_theta_2D(const FITTER& fitter, const reco::BeamSpot& bs, const LORENTZ_VEC& p4) {
0056     if (!fitter.success())
0057       return -2;
0058     GlobalPoint point = fitter.fitted_vtx();
0059     auto bs_pos = bs.position(point.z());
0060     math::XYZVector delta(point.x() - bs_pos.x(), point.y() - bs_pos.y(), 0.);
0061     math::XYZVector pt(p4.px(), p4.py(), 0.);
0062     double den = (delta.R() * pt.R());
0063     return (den != 0.) ? delta.Dot(pt) / den : -2;
0064   }
0065
0066   template <typename FITTER>
0067   inline Measurement1D l_xy(const FITTER& fitter, const reco::BeamSpot& bs) {
0068     if (!fitter.success())
0069       return {-2, -2};
0070     GlobalPoint point = fitter.fitted_vtx();
0071     GlobalError err = fitter.fitted_vtx_uncertainty();
0072     auto bs_pos = bs.position(point.z());
0073     GlobalPoint delta(point.x() - bs_pos.x(), point.y() - bs_pos.y(), 0.);
0074     return {delta.perp(), sqrt(err.rerr(delta))};
0075   }
0076
0077   /*
0078 inline GlobalPoint FlightDistVector (const reco::BeamSpot & bm, GlobalPoint
0079 Bvtx)
0080 {
0081    GlobalPoint Dispbeamspot(-1*( (bm.x0()-Bvtx.x()) + (Bvtx.z()-bm.z0()) *
0082 bm.dxdz()), -1*( (bm.y0()-Bvtx.y()) + (Bvtx.z()-bm.z0()) * bm.dydz()), 0);
0083    return std::move(Dispbeamspot);
0084 }
0085 */
0086
0087   inline float CosA(GlobalPoint& dist, ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double>>& Bp4) {
0088     math::XYZVector vperp(dist.x(), dist.y(), 0);
0089     math::XYZVector pperp(Bp4.Px(), Bp4.Py(), 0);
0090     return vperp.Dot(pperp) / (vperp.R() * pperp.R());
0091   }
0092
0093   inline std::pair<double, double> computeDCA(const reco::TransientTrack& trackTT, const reco::BeamSpot& beamSpot) {
0094     double DCABS = -1.;
0095     double DCABSErr = -1.;
0096
0097     TrajectoryStateClosestToPoint theDCAXBS = trackTT.trajectoryStateClosestToPoint(
0098         GlobalPoint(beamSpot.position().x(), beamSpot.position().y(), beamSpot.position().z()));
0099     if (theDCAXBS.isValid()) {
0100       DCABS = theDCAXBS.perigeeParameters().transverseImpactParameter();
0101       DCABSErr = theDCAXBS.perigeeError().transverseImpactParameterError();
0102     }
0103
0104     return std::make_pair(DCABS, DCABSErr);
0105   }
0106
0107   inline bool track_to_lepton_match(edm::Ptr<reco::Candidate> l_ptr, auto iso_tracks_id, unsigned int iTrk) {
0108     for (unsigned int i = 0; i < l_ptr->numberOfSourceCandidatePtrs(); ++i) {
0109       if (!((l_ptr->sourceCandidatePtr(i)).isNonnull() && (l_ptr->sourceCandidatePtr(i)).isAvailable()))
0110         continue;
0111       const edm::Ptr<reco::Candidate>& source = l_ptr->sourceCandidatePtr(i);
0112       if (source.id() == iso_tracks_id && source.key() == iTrk) {
0113         return true;
0114       }
0115     }
0116     return false;
0117   }
0118
0119   inline std::pair<bool, Measurement1D> absoluteImpactParameter(const TrajectoryStateOnSurface& tsos,
0120                                                                 RefCountedKinematicVertex vertex,
0121                                                                 VertexDistance& distanceComputer) {
0122     if (!tsos.isValid()) {
0123       return std::pair<bool, Measurement1D>(false, Measurement1D(0., 0.));
0124     }
0125     GlobalPoint refPoint = tsos.globalPosition();
0126     GlobalError refPointErr = tsos.cartesianError().position();
0127     GlobalPoint vertexPosition = vertex->vertexState().position();
0128     GlobalError vertexPositionErr = RecoVertex::convertError(vertex->vertexState().error());
0129     return std::pair<bool, Measurement1D>(
0130         true,
0131         distanceComputer.distance(VertexState(vertexPosition, vertexPositionErr), VertexState(refPoint, refPointErr)));
0132   }
0133
0134   inline std::pair<bool, Measurement1D> absoluteImpactParameter3D(const TrajectoryStateOnSurface& tsos,
0135                                                                   RefCountedKinematicVertex vertex) {
0136     VertexDistance3D dist;
0137     return absoluteImpactParameter(tsos, vertex, dist);
0138   }
0139
0140   inline std::pair<bool, Measurement1D> absoluteTransverseImpactParameter(const TrajectoryStateOnSurface& tsos,
0141                                                                           RefCountedKinematicVertex vertex) {
0142     VertexDistanceXY dist;
0143     return absoluteImpactParameter(tsos, vertex, dist);
0144   }
0145
0146   inline std::pair<bool, Measurement1D> signedImpactParameter3D(const TrajectoryStateOnSurface& tsos,
0147                                                                 RefCountedKinematicVertex vertex,
0148                                                                 const reco::BeamSpot& bs,
0149                                                                 double pv_z) {
0150     VertexDistance3D dist;
0151
0152     std::pair<bool, Measurement1D> result = absoluteImpactParameter(tsos, vertex, dist);
0153     if (!result.first)
0154       return result;
0155
0156     // Compute Sign
0157     auto bs_pos = bs.position(vertex->vertexState().position().z());
0158     GlobalPoint impactPoint = tsos.globalPosition();
0159     GlobalVector IPVec(impactPoint.x() - vertex->vertexState().position().x(),
0160                        impactPoint.y() - vertex->vertexState().position().y(),
0161                        impactPoint.z() - vertex->vertexState().position().z());
0162
0163     GlobalVector direction(vertex->vertexState().position().x() - bs_pos.x(),
0164                            vertex->vertexState().position().y() - bs_pos.y(),
0165                            vertex->vertexState().position().z() - pv_z);
0166
0167     double prod = IPVec.dot(direction);
0168     double sign = (prod >= 0) ? 1. : -1.;
0169
0170     // Apply sign to the result
0171     return std::pair<bool, Measurement1D>(result.first,
0172                                           Measurement1D(sign * result.second.value(), result.second.error()));
0173   }
0174
0175   inline std::pair<bool, Measurement1D> signedTransverseImpactParameter(const TrajectoryStateOnSurface& tsos,
0176                                                                         RefCountedKinematicVertex vertex,
0177                                                                         const reco::BeamSpot& bs) {
0178     VertexDistanceXY dist;
0179
0180     std::pair<bool, Measurement1D> result = absoluteImpactParameter(tsos, vertex, dist);
0181     if (!result.first)
0182       return result;
0183
0184     // Compute Sign
0185     auto bs_pos = bs.position(vertex->vertexState().position().z());
0186     GlobalPoint impactPoint = tsos.globalPosition();
0187     GlobalVector IPVec(impactPoint.x() - vertex->vertexState().position().x(),
0188                        impactPoint.y() - vertex->vertexState().position().y(),
0189                        0.);
0190     GlobalVector direction(
0191         vertex->vertexState().position().x() - bs_pos.x(), vertex->vertexState().position().y() - bs_pos.y(), 0);
0192
0193     double prod = IPVec.dot(direction);
0194     double sign = (prod >= 0) ? 1. : -1.;
0195
0196     // Apply sign to the result
0197     return std::pair<bool, Measurement1D>(result.first,
0198                                           Measurement1D(sign * result.second.value(), result.second.error()));
0199   }
0200
0201   inline std::vector<float> TrackerIsolation(edm::Handle<pat::CompositeCandidateCollection>& tracks,
0202                                              pat::CompositeCandidate& B,
0203                                              std::vector<std::string>& dnames) {
0204     std::vector<float> iso(dnames.size(), 0);
0205     for (size_t k_idx = 0; k_idx < tracks->size(); ++k_idx) {
0206       edm::Ptr<pat::CompositeCandidate> trk_ptr(tracks, k_idx);
0207       for (size_t iname = 0; iname < dnames.size(); ++iname) {
0208         float dr = deltaR(B.userFloat("fitted_" + dnames[iname] + "_eta"),
0209                           B.userFloat("fitted_" + dnames[iname] + "_phi"),
0210                           trk_ptr->eta(),
0211                           trk_ptr->phi());
0212         if (dr > 0 && dr < 0.4)
0213           iso[iname] += trk_ptr->pt();
0214       }
0215     }
0216     return iso;
0217   }
0218
0219 }  // namespace bph
0220 #endif