EgammaTools/src/LowPtConversion.cc

0001 #include "RecoEgamma/EgammaTools/interface/LowPtConversion.h"
0002
0003 ////////////////////////////////////////////////////////////////////////////////
0004 // Matched to any conversion (without selections)
0005 //
0006 bool LowPtConversion::wpOpen() const { return matched_; }
0007
0008 ////////////////////////////////////////////////////////////////////////////////
0009 // Nancy's baseline selections for conversions
0010 // Based on: https://github.com/CMSBParking/BParkingNANO/blob/b2664ed/BParkingNano/plugins/ConversionSelector.cc#L253-L300
0011 bool LowPtConversion::wpLoose() const {
0012   return (wpOpen() && ntracks_ == 2 && valid_ && quality_high_purity_ && chi2prob_ > 0.0005);
0013 }
0014
0015 ////////////////////////////////////////////////////////////////////////////////
0016 // Nancy's selection for analysis of conversions
0017 // Based on: slide 20 of https://indico.cern.ch/event/814402/contributions/3401312/
0018 bool LowPtConversion::wpTight() const {
0019   return (wpLoose() && sum_nhits_before_vtx_ <= 1 && l_xy_ > 0. && mass_from_conv_ > 0. &&  // sanity check
0020           mass_from_conv_ < 0.05);
0021 }
0022
0023 ////////////////////////////////////////////////////////////////////////////////
0024 // adds minimal set of flags to electron userData
0025 void LowPtConversion::addUserVars(pat::Electron& ele) const {
0026   ele.addUserInt("convOpen", matched_ ? 1 : 0);
0027   ele.addUserInt("convLoose", wpLoose() ? 1 : 0);
0028   ele.addUserInt("convTight", wpTight() ? 1 : 0);
0029   ele.addUserInt("convLead", matched_lead_.isNonnull() ? 1 : 0);
0030   ele.addUserInt("convTrail", matched_trail_.isNonnull() ? 1 : 0);
0031   if (ele.hasUserInt("convExtra") == false) {
0032     ele.addUserInt("convExtra", 0);
0033   }
0034 }
0035
0036 ////////////////////////////////////////////////////////////////////////////////
0037 // adds all variables to electron userData
0038 void LowPtConversion::addExtraUserVars(pat::Electron& ele) const {
0039   // Flag that indicates if extra variables are added to electron userData
0040   ele.addUserInt("convExtra", 1, true);  // overwrite
0041
0042   // quality
0043   ele.addUserInt("convValid", valid_ ? 1 : 0);
0044   ele.addUserFloat("convChi2Prob", chi2prob_);
0045   ele.addUserInt("convQualityHighPurity", quality_high_purity_ ? 1 : 0);
0046   ele.addUserInt("convQualityHighEff", quality_high_efficiency_ ? 1 : 0);
0047
0048   // tracks
0049   ele.addUserInt("convTracksN", ntracks_);
0050   ele.addUserFloat("convMinTrkPt", min_trk_pt_);
0051   ele.addUserInt("convLeadIdx", ilead_);
0052   ele.addUserInt("convTrailIdx", itrail_);
0053
0054   // displacement
0055   ele.addUserFloat("convLxy", l_xy_);
0056   ele.addUserFloat("convVtxRadius", vtx_radius_);
0057
0058   // invariant mass
0059   ele.addUserFloat("convMass", mass_from_conv_);
0060   ele.addUserFloat("convMassFromPin", mass_from_Pin_);
0061   ele.addUserFloat("convMassBeforeFit", mass_before_fit_);
0062   ele.addUserFloat("convMassAfterFit", mass_after_fit_);
0063
0064   // hits before vertex
0065   ele.addUserInt("convLeadNHitsBeforeVtx", lead_nhits_before_vtx_);
0066   ele.addUserInt("convTrailNHitsBeforeVtx", trail_nhits_before_vtx_);
0067   ele.addUserInt("convMaxNHitsBeforeVtx", max_nhits_before_vtx_);
0068   ele.addUserInt("convSumNHitsBeforeVtx", sum_nhits_before_vtx_);
0069   ele.addUserInt("convDeltaExpectedNHitsInner", delta_expected_nhits_inner_);
0070
0071   // opening angle
0072   ele.addUserFloat("convDeltaCotFromPin", delta_cot_from_Pin_);
0073 }
0074
0075 ////////////////////////////////////////////////////////////////////////////////
0076 //
0077 bool LowPtConversion::match(const reco::BeamSpot& beamSpot,
0078                             const reco::ConversionCollection& conversions,
0079                             const pat::Electron& ele) {
0080   // Iterate through conversions and calculate quantities (requirement from Nancy)
0081   for (const auto& conv : conversions) {
0082     // Filter
0083     if (conv.tracks().size() != 2) {
0084       continue;
0085     }
0086
0087     // Quality
0088     valid_ = conv.conversionVertex().isValid();                                                         // (=true)
0089     chi2prob_ = ChiSquaredProbability(conv.conversionVertex().chi2(), conv.conversionVertex().ndof());  // (<0.005)
0090     quality_high_purity_ = conv.quality(reco::Conversion::highPurity);                                  // (=true)
0091     quality_high_efficiency_ = conv.quality(reco::Conversion::highEfficiency);                          // (none)
0092
0093     // Tracks
0094     ntracks_ = conv.tracks().size();  // (=2)
0095     min_trk_pt_ = -1.;                // (>0.5)
0096     for (const auto& trk : conv.tracks()) {
0097       if (trk.isNonnull() && trk.isAvailable() && (min_trk_pt_ < 0. || trk->pt() < min_trk_pt_)) {
0098         min_trk_pt_ = trk->pt();
0099       }
0100     }
0101     ilead_ = -1;
0102     itrail_ = -1;
0103     if (conv.tracks().size() == 2) {
0104       const edm::RefToBase<reco::Track>& trk1 = conv.tracks().front();
0105       const edm::RefToBase<reco::Track>& trk2 = conv.tracks().back();
0106       if (trk1.isNonnull() && trk1.isAvailable() && trk2.isNonnull() && trk2.isAvailable()) {
0107         if (trk1->pt() > trk2->pt()) {
0108           ilead_ = 0;
0109           itrail_ = 1;
0110         } else {
0111           ilead_ = 1;
0112           itrail_ = 0;
0113         }
0114       }
0115     }
0116
0117     // Transverse displacement (with respect to beamspot) and vertex radius
0118     math::XYZVectorF p_refitted = conv.refittedPairMomentum();
0119     float dx = conv.conversionVertex().x() - beamSpot.x0();
0120     float dy = conv.conversionVertex().y() - beamSpot.y0();
0121     l_xy_ = (p_refitted.x() * dx + p_refitted.y() * dy) / p_refitted.rho();
0122     vtx_radius_ = sqrt(conv.conversionVertex().position().perp2());  // (1.5<r<4.)
0123
0124     // invariant mass from track pair from conversion
0125     mass_from_conv_ = conv.pairInvariantMass();
0126
0127     // Invariant mass from Pin before fit to common vertex
0128     if (conv.tracksPin().size() >= 2 && ilead_ > -1 && itrail_ > -1) {
0129       math::XYZVectorF lead_Pin = conv.tracksPin().at(ilead_);
0130       math::XYZVectorF trail_Pin = conv.tracksPin().at(itrail_);
0131       mass_from_Pin_ = mee(lead_Pin.x(), lead_Pin.y(), lead_Pin.z(), trail_Pin.x(), trail_Pin.y(), trail_Pin.z());
0132       // Opening angle
0133       delta_cot_from_Pin_ = 1. / tan(trail_Pin.theta()) - 1. / tan(lead_Pin.theta());
0134     }
0135
0136     // Invariant mass before fit to common vertex
0137     if (conv.tracks().size() >= 2 && ilead_ > -1 && itrail_ > -1) {
0138       auto lead_before_vtx_fit = conv.tracks().at(ilead_)->momentum();
0139       auto trail_before_vtx_fit = conv.tracks().at(itrail_)->momentum();
0140       mass_before_fit_ = mee(lead_before_vtx_fit.x(),
0141                              lead_before_vtx_fit.y(),
0142                              lead_before_vtx_fit.z(),
0143                              trail_before_vtx_fit.x(),
0144                              trail_before_vtx_fit.y(),
0145                              trail_before_vtx_fit.z());
0146     }
0147
0148     // Invariant mass after the fit to common vertex
0149     if (conv.conversionVertex().refittedTracks().size() >= 2 && ilead_ > -1 && itrail_ > -1) {
0150       auto const& lead_after_vtx_fit = conv.conversionVertex().refittedTracks().at(ilead_);
0151       auto const& trail_after_vtx_fit = conv.conversionVertex().refittedTracks().at(itrail_);
0152       mass_after_fit_ = mee(lead_after_vtx_fit.px(),
0153                             lead_after_vtx_fit.py(),
0154                             lead_after_vtx_fit.pz(),
0155                             trail_after_vtx_fit.px(),
0156                             trail_after_vtx_fit.py(),
0157                             trail_after_vtx_fit.pz());
0158       // Difference in expeted hits
0159       delta_expected_nhits_inner_ =
0160           lead_after_vtx_fit.hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS) -
0161           trail_after_vtx_fit.hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS);
0162     }
0163
0164     // Hits prior to vertex
0165     if (ilead_ > -1 && itrail_ > -1) {
0166       auto const& nHits = conv.nHitsBeforeVtx();
0167       bool enoughTracks = nHits.size() > 1;
0168       lead_nhits_before_vtx_ = enoughTracks ? nHits.at(ilead_) : 0;
0169       trail_nhits_before_vtx_ = enoughTracks ? nHits.at(itrail_) : 0;
0170       max_nhits_before_vtx_ = enoughTracks ? std::max(nHits[0], nHits[1]) : 0;
0171       sum_nhits_before_vtx_ = enoughTracks ? nHits[0] + nHits[1] : 0;
0172     }
0173
0174     // Attempt to match conversion track to electron
0175     for (uint itrk = 0; itrk < conv.tracks().size(); ++itrk) {
0176       const edm::RefToBase<reco::Track> trk = conv.tracks()[itrk];
0177       if (trk.isNull()) {
0178         continue;
0179       }
0180       reco::GsfTrackRef ref = ele.core()->gsfTrack();
0181       reco::GsfTrackRef gsf = ele.gsfTrack();
0182       if (gsf.isNull()) {
0183         continue;
0184       }
0185       if (ref.id() == trk.id() && ref.key() == trk.key()) {
0186         matched_ = true;
0187         if (static_cast<int>(itrk) == ilead_) {
0188           matched_lead_ = trk;
0189         }
0190         if (static_cast<int>(itrk) == itrail_) {
0191           matched_trail_ = trk;
0192         }
0193       }
0194     }  // track loop
0195   }  // conversions loop
0196
0197   return matched_;
0198 }
0199
0200 ////////////////////////////////////////////////////////////////////////////////
0201 //
0202 float LowPtConversion::mee(float px1, float py1, float pz1, float px2, float py2, float pz2) {
0203   const float m = 0.000511;
0204   const float px = px1 + px2;
0205   const float py = py1 + py2;
0206   const float pz = pz1 + pz2;
0207   const float p1 = px1 * px1 + py1 * py1 + pz1 * pz1;
0208   const float p2 = px2 * px2 + py2 * py2 + pz2 * pz2;
0209   const float e = sqrt(p1 + m * m) + sqrt(p2 + m * m);
0210   const float mass = (e * e - px * px - py * py - pz * pz);
0211   return mass > 0. ? sqrt(mass) : -1.;
0212 }