EgammaCandidates/src/Conversion.cc

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 }