File indexing completed on 2023-03-17 10:49:44
0001 #include "DataFormats/EgammaCandidates/interface/Conversion.h"
0002 #include "DataFormats/TrackReco/interface/Track.h"
0003 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0004 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0005
0006 using namespace reco;
0007
0008 Conversion::Conversion(const reco::CaloClusterPtrVector& sc,
0009 const std::vector<reco::TrackRef>& tr,
0010 const std::vector<math::XYZPointF>& trackPositionAtEcal,
0011 const reco::Vertex& convVtx,
0012 const std::vector<reco::CaloClusterPtr>& matchingBC,
0013 const float DCA,
0014 const std::vector<math::XYZPointF>& innPoint,
0015 const std::vector<math::XYZVectorF>& trackPin,
0016 const std::vector<math::XYZVectorF>& trackPout,
0017 const float mva,
0018 ConversionAlgorithm algo)
0019 :
0020
0021 caloCluster_(sc),
0022 trackToBaseRefs_(),
0023 thePositionAtEcal_(trackPositionAtEcal),
0024 theConversionVertex_(convVtx),
0025 theMatchingBCs_(matchingBC),
0026 theTrackInnerPosition_(innPoint),
0027 theTrackPin_(trackPin),
0028 theTrackPout_(trackPout),
0029 theMinDistOfApproach_(DCA),
0030 theMVAout_(mva),
0031 qualityMask_(0),
0032 nSharedHits_(0),
0033 algorithm_(algo) {
0034 trackToBaseRefs_.reserve(tr.size());
0035 for (auto const& track : tr) {
0036 trackToBaseRefs_.emplace_back(track);
0037 }
0038 }
0039
0040 Conversion::Conversion(const reco::CaloClusterPtrVector& sc,
0041 const std::vector<edm::RefToBase<reco::Track> >& tr,
0042 const std::vector<math::XYZPointF>& trackPositionAtEcal,
0043 const reco::Vertex& convVtx,
0044 const std::vector<reco::CaloClusterPtr>& matchingBC,
0045 const float DCA,
0046 const std::vector<math::XYZPointF>& innPoint,
0047 const std::vector<math::XYZVectorF>& trackPin,
0048 const std::vector<math::XYZVectorF>& trackPout,
0049 const std::vector<uint8_t>& nHitsBeforeVtx,
0050 const std::vector<Measurement1DFloat>& dlClosestHitToVtx,
0051 uint8_t nSharedHits,
0052 const float mva,
0053 ConversionAlgorithm algo)
0054 :
0055
0056 caloCluster_(sc),
0057 trackToBaseRefs_(tr),
0058 thePositionAtEcal_(trackPositionAtEcal),
0059 theConversionVertex_(convVtx),
0060 theMatchingBCs_(matchingBC),
0061 theTrackInnerPosition_(innPoint),
0062 theTrackPin_(trackPin),
0063 theTrackPout_(trackPout),
0064 nHitsBeforeVtx_(nHitsBeforeVtx),
0065 dlClosestHitToVtx_(dlClosestHitToVtx),
0066 theMinDistOfApproach_(DCA),
0067 theMVAout_(mva),
0068 qualityMask_(0),
0069 nSharedHits_(nSharedHits),
0070 algorithm_(algo) {}
0071
0072 Conversion::Conversion(const reco::CaloClusterPtrVector& sc,
0073 const std::vector<reco::TrackRef>& tr,
0074 const reco::Vertex& convVtx,
0075 ConversionAlgorithm algo)
0076 : caloCluster_(sc),
0077 trackToBaseRefs_(),
0078 theConversionVertex_(convVtx),
0079 qualityMask_(0),
0080 nSharedHits_(0),
0081 algorithm_(algo) {
0082 trackToBaseRefs_.reserve(tr.size());
0083 for (auto const& track : tr) {
0084 trackToBaseRefs_.emplace_back(track);
0085 }
0086
0087 theMinDistOfApproach_ = 9999.;
0088 theMVAout_ = 9999.;
0089 thePositionAtEcal_.push_back(math::XYZPointF(0., 0., 0.));
0090 thePositionAtEcal_.push_back(math::XYZPointF(0., 0., 0.));
0091 theTrackInnerPosition_.push_back(math::XYZPointF(0., 0., 0.));
0092 theTrackInnerPosition_.push_back(math::XYZPointF(0., 0., 0.));
0093 theTrackPin_.push_back(math::XYZVectorF(0., 0., 0.));
0094 theTrackPin_.push_back(math::XYZVectorF(0., 0., 0.));
0095 theTrackPout_.push_back(math::XYZVectorF(0., 0., 0.));
0096 theTrackPout_.push_back(math::XYZVectorF(0., 0., 0.));
0097 }
0098
0099 Conversion::Conversion(const reco::CaloClusterPtrVector& sc,
0100 const std::vector<edm::RefToBase<reco::Track> >& tr,
0101 const reco::Vertex& convVtx,
0102 ConversionAlgorithm algo)
0103 : caloCluster_(sc),
0104 trackToBaseRefs_(tr),
0105 theConversionVertex_(convVtx),
0106 qualityMask_(0),
0107 nSharedHits_(0),
0108 algorithm_(algo) {
0109 theMinDistOfApproach_ = 9999.;
0110 theMVAout_ = 9999.;
0111 thePositionAtEcal_.push_back(math::XYZPointF(0., 0., 0.));
0112 thePositionAtEcal_.push_back(math::XYZPointF(0., 0., 0.));
0113 theTrackInnerPosition_.push_back(math::XYZPointF(0., 0., 0.));
0114 theTrackInnerPosition_.push_back(math::XYZPointF(0., 0., 0.));
0115 theTrackPin_.push_back(math::XYZVectorF(0., 0., 0.));
0116 theTrackPin_.push_back(math::XYZVectorF(0., 0., 0.));
0117 theTrackPout_.push_back(math::XYZVectorF(0., 0., 0.));
0118 theTrackPout_.push_back(math::XYZVectorF(0., 0., 0.));
0119 }
0120
0121 Conversion::Conversion() {
0122 algorithm_ = 0;
0123 qualityMask_ = 0;
0124 theMinDistOfApproach_ = 9999.;
0125 nSharedHits_ = 0;
0126 theMVAout_ = 9999.;
0127 thePositionAtEcal_.push_back(math::XYZPointF(0., 0., 0.));
0128 thePositionAtEcal_.push_back(math::XYZPointF(0., 0., 0.));
0129 theTrackInnerPosition_.push_back(math::XYZPointF(0., 0., 0.));
0130 theTrackInnerPosition_.push_back(math::XYZPointF(0., 0., 0.));
0131 theTrackPin_.push_back(math::XYZVectorF(0., 0., 0.));
0132 theTrackPin_.push_back(math::XYZVectorF(0., 0., 0.));
0133 theTrackPout_.push_back(math::XYZVectorF(0., 0., 0.));
0134 theTrackPout_.push_back(math::XYZVectorF(0., 0., 0.));
0135 }
0136
0137 std::string const Conversion::algoNames[] = {"undefined", "ecalSeeded", "trackerOnly", "mixed", "pflow"};
0138
0139 Conversion::ConversionAlgorithm Conversion::algoByName(const std::string& name) {
0140 ConversionAlgorithm size = algoSize;
0141 int index = std::find(algoNames, algoNames + size, name) - algoNames;
0142 if (index == size)
0143 return undefined;
0144
0145 return ConversionAlgorithm(index);
0146 }
0147
0148 Conversion* Conversion::clone() const { return new Conversion(*this); }
0149
0150 std::vector<edm::RefToBase<reco::Track> > const& Conversion::tracks() const { return trackToBaseRefs_; }
0151
0152 bool Conversion::isConverted() const {
0153 if (this->nTracks() == 2)
0154 return true;
0155 else
0156 return false;
0157 }
0158
0159 double Conversion::pairInvariantMass() const {
0160 double invMass = -99.;
0161 const float mElec = 0.000511;
0162 if (nTracks() == 2) {
0163 double px = tracksPin()[0].x() + tracksPin()[1].x();
0164 double py = tracksPin()[0].y() + tracksPin()[1].y();
0165 double pz = tracksPin()[0].z() + tracksPin()[1].z();
0166 double mom1 = tracksPin()[0].Mag2();
0167 double mom2 = tracksPin()[1].Mag2();
0168 double e = sqrt(mom1 + mElec * mElec) + sqrt(mom2 + mElec * mElec);
0169 invMass = (e * e - px * px - py * py - pz * pz);
0170 if (invMass > 0)
0171 invMass = sqrt(invMass);
0172 else
0173 invMass = -1;
0174 }
0175
0176 return invMass;
0177 }
0178
0179 double Conversion::pairCotThetaSeparation() const {
0180 double dCotTheta = -99.;
0181
0182 if (nTracks() == 2) {
0183 double theta1 = this->tracksPin()[0].Theta();
0184 double theta2 = this->tracksPin()[1].Theta();
0185 dCotTheta = 1. / tan(theta1) - 1. / tan(theta2);
0186 }
0187
0188 return dCotTheta;
0189 }
0190
0191 math::XYZVectorF Conversion::pairMomentum() const {
0192 if (nTracks() == 2) {
0193 return this->tracksPin()[0] + this->tracksPin()[1];
0194 }
0195 return math::XYZVectorF(0., 0., 0.);
0196 }
0197
0198 math::XYZTLorentzVectorF Conversion::refittedPair4Momentum() const {
0199 math::XYZTLorentzVectorF p4;
0200 if (this->conversionVertex().isValid())
0201 p4 = this->conversionVertex().p4(0.000511, 0.5);
0202
0203 return p4;
0204 }
0205
0206 math::XYZVectorF Conversion::refittedPairMomentum() const {
0207 if (this->conversionVertex().isValid()) {
0208 return this->refittedPair4Momentum().Vect();
0209 }
0210 return math::XYZVectorF(0., 0., 0.);
0211 }
0212
0213 double Conversion::EoverP() const {
0214 double ep = -99.;
0215
0216 if (nTracks() > 0) {
0217 unsigned int size = this->caloCluster().size();
0218 float etot = 0.;
0219 for (unsigned int i = 0; i < size; i++) {
0220 etot += caloCluster()[i]->energy();
0221 }
0222 if (this->pairMomentum().Mag2() != 0)
0223 ep = etot / sqrt(this->pairMomentum().Mag2());
0224 }
0225
0226 return ep;
0227 }
0228
0229 double Conversion::EoverPrefittedTracks() const {
0230 double ep = -99.;
0231
0232 if (nTracks() > 0) {
0233 unsigned int size = this->caloCluster().size();
0234 float etot = 0.;
0235 for (unsigned int i = 0; i < size; i++) {
0236 etot += caloCluster()[i]->energy();
0237 }
0238 if (this->refittedPairMomentum().Mag2() != 0)
0239 ep = etot / sqrt(this->refittedPairMomentum().Mag2());
0240 }
0241
0242 return ep;
0243 }
0244
0245 std::vector<double> Conversion::tracksSigned_d0() const {
0246 std::vector<double> result;
0247 result.reserve(tracks().size());
0248
0249 for (auto const& track : tracks()) {
0250 result.emplace_back(track->d0() * track->charge());
0251 }
0252 return result;
0253 }
0254
0255 double Conversion::dPhiTracksAtVtx() const {
0256 double result = -99;
0257 if (nTracks() == 2) {
0258 result = tracksPin()[0].phi() - tracksPin()[1].phi();
0259 if (result > pi) {
0260 result = result - twopi;
0261 }
0262 if (result < -pi) {
0263 result = result + twopi;
0264 }
0265 }
0266
0267 return result;
0268 }
0269
0270 double Conversion::dPhiTracksAtEcal() const {
0271 double result = -99.;
0272
0273 if (nTracks() == 2 && bcMatchingWithTracks()[0].isNonnull() && bcMatchingWithTracks()[1].isNonnull()) {
0274 float recoPhi1 = ecalImpactPosition()[0].phi();
0275 if (recoPhi1 > pi) {
0276 recoPhi1 = recoPhi1 - twopi;
0277 }
0278 if (recoPhi1 < -pi) {
0279 recoPhi1 = recoPhi1 + twopi;
0280 }
0281
0282 float recoPhi2 = ecalImpactPosition()[1].phi();
0283 if (recoPhi2 > pi) {
0284 recoPhi2 = recoPhi2 - twopi;
0285 }
0286 if (recoPhi2 < -pi) {
0287 recoPhi2 = recoPhi2 + twopi;
0288 }
0289
0290 result = recoPhi1 - recoPhi2;
0291
0292 if (result > pi) {
0293 result = result - twopi;
0294 }
0295 if (result < -pi) {
0296 result = result + twopi;
0297 }
0298 }
0299
0300 return result;
0301 }
0302
0303 double Conversion::dEtaTracksAtEcal() const {
0304 double result = -99.;
0305
0306 if (nTracks() == 2 && bcMatchingWithTracks()[0].isNonnull() && bcMatchingWithTracks()[1].isNonnull()) {
0307 result = ecalImpactPosition()[0].eta() - ecalImpactPosition()[1].eta();
0308 }
0309
0310 return result;
0311 }
0312
0313 double Conversion::dxy(const math::XYZPoint& myBeamSpot) const {
0314 const reco::Vertex& vtx = conversionVertex();
0315 if (!vtx.isValid())
0316 return -9999.;
0317
0318 math::XYZVectorF mom = refittedPairMomentum();
0319
0320 double dxy = (-(vtx.x() - myBeamSpot.x()) * mom.y() + (vtx.y() - myBeamSpot.y()) * mom.x()) / mom.rho();
0321 return dxy;
0322 }
0323
0324 double Conversion::dz(const math::XYZPoint& myBeamSpot) const {
0325 const reco::Vertex& vtx = conversionVertex();
0326 if (!vtx.isValid())
0327 return -9999.;
0328
0329 math::XYZVectorF mom = refittedPairMomentum();
0330
0331 double dz =
0332 (vtx.z() - myBeamSpot.z()) -
0333 ((vtx.x() - myBeamSpot.x()) * mom.x() + (vtx.y() - myBeamSpot.y()) * mom.y()) / mom.rho() * mom.z() / mom.rho();
0334 return dz;
0335 }
0336
0337 double Conversion::lxy(const math::XYZPoint& myBeamSpot) const {
0338 const reco::Vertex& vtx = conversionVertex();
0339 if (!vtx.isValid())
0340 return -9999.;
0341
0342 math::XYZVectorF mom = refittedPairMomentum();
0343
0344 double dbsx = vtx.x() - myBeamSpot.x();
0345 double dbsy = vtx.y() - myBeamSpot.y();
0346 double lxy = (mom.x() * dbsx + mom.y() * dbsy) / mom.rho();
0347 return lxy;
0348 }
0349
0350 double Conversion::lz(const math::XYZPoint& myBeamSpot) const {
0351 const reco::Vertex& vtx = conversionVertex();
0352 if (!vtx.isValid())
0353 return -9999.;
0354
0355 math::XYZVectorF mom = refittedPairMomentum();
0356
0357 double lz = (vtx.z() - myBeamSpot.z()) * mom.z() / std::abs(mom.z());
0358 return lz;
0359 }
