Onia2MuMu/src/OniaPhotonConversionProducer.cc

0001 #include "HeavyFlavorAnalysis/Onia2MuMu/interface/OniaPhotonConversionProducer.h"
0002 #include "DataFormats/EgammaCandidates/interface/Conversion.h"
0003 #include "DataFormats/PatCandidates/interface/CompositeCandidate.h"
0004 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0005 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0006 #include "DataFormats/TrackReco/interface/Track.h"
0007 #include "DataFormats/TrackReco/interface/TrackBase.h"
0008
0009 #include "CommonTools/Utils/interface/StringToEnumValue.h"
0010 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0011
0012 #include "RecoVertex/VertexTools/interface/VertexDistanceXY.h"
0013 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0014
0015 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0016
0017 #include <TMath.h>
0018 #include <vector>
0019
0020 // to order from high to low ProbChi2
0021 bool ConversionLessByChi2(const reco::Conversion& c1, const reco::Conversion& c2) {
0022   return TMath::Prob(c1.conversionVertex().chi2(), c1.conversionVertex().ndof()) >
0023          TMath::Prob(c2.conversionVertex().chi2(), c2.conversionVertex().ndof());
0024 }
0025
0026 bool ConversionEqualByTrack(const reco::Conversion& c1, const reco::Conversion& c2) {
0027   bool atLeastOneInCommon = false;
0028   for (auto const& tk1 : c1.tracks()) {
0029     for (auto const& tk2 : c2.tracks()) {
0030       if (tk1 == tk2) {
0031         atLeastOneInCommon = true;
0032         break;
0033       }
0034     }
0035   }
0036   return atLeastOneInCommon;
0037 }
0038
0039 bool lt_(std::pair<double, short> a, std::pair<double, short> b) { return a.first < b.first; }
0040
0041 // define operator== for conversions, those with at least one track in common
0042 namespace reco {
0043   bool operator==(const reco::Conversion& c1, const reco::Conversion& c2) {
0044     return c1.tracks()[0] == c2.tracks()[0] || c1.tracks()[1] == c2.tracks()[1] || c1.tracks()[1] == c2.tracks()[0] ||
0045            c1.tracks()[0] == c2.tracks()[1];
0046   }
0047 }  // namespace reco
0048
0049 OniaPhotonConversionProducer::OniaPhotonConversionProducer(const edm::ParameterSet& ps) {
0050   convCollectionToken_ = consumes<reco::ConversionCollection>(ps.getParameter<edm::InputTag>("conversions"));
0051   thePVsToken_ = consumes<reco::VertexCollection>(ps.getParameter<edm::InputTag>("primaryVertexTag"));
0052
0053   wantTkVtxCompatibility_ = ps.getParameter<bool>("wantTkVtxCompatibility");
0054   sigmaTkVtxComp_ = ps.getParameter<uint32_t>("sigmaTkVtxComp");
0055   wantCompatibleInnerHits_ = ps.getParameter<bool>("wantCompatibleInnerHits");
0056   TkMinNumOfDOF_ = ps.getParameter<uint32_t>("TkMinNumOfDOF");
0057
0058   wantHighpurity_ = ps.getParameter<bool>("wantHighpurity");
0059   _vertexChi2ProbCut = ps.getParameter<double>("vertexChi2ProbCut");
0060   _trackchi2Cut = ps.getParameter<double>("trackchi2Cut");
0061   _minDistanceOfApproachMinCut = ps.getParameter<double>("minDistanceOfApproachMinCut");
0062   _minDistanceOfApproachMaxCut = ps.getParameter<double>("minDistanceOfApproachMaxCut");
0063
0064   pfCandidateCollectionToken_ = consumes<reco::PFCandidateCollection>(ps.getParameter<edm::InputTag>("pfcandidates"));
0065
0066   pi0OnlineSwitch_ = ps.getParameter<bool>("pi0OnlineSwitch");
0067   pi0SmallWindow_ = ps.getParameter<std::vector<double>>("pi0SmallWindow");
0068   pi0LargeWindow_ = ps.getParameter<std::vector<double>>("pi0LargeWindow");
0069
0070   std::string algo = ps.getParameter<std::string>("convAlgo");
0071   convAlgo_ = StringToEnumValue<reco::Conversion::ConversionAlgorithm>(algo);
0072
0073   std::vector<std::string> qual = ps.getParameter<std::vector<std::string>>("convQuality");
0074   if (!qual[0].empty())
0075     convQuality_ = StringToEnumValue<reco::Conversion::ConversionQuality>(qual);
0076
0077   convSelectionCuts_ = ps.getParameter<std::string>("convSelection");
0078   convSelection_ = std::make_unique<StringCutObjectSelector<reco::Conversion>>(convSelectionCuts_);
0079   produces<pat::CompositeCandidateCollection>("conversions");
0080 }
0081
0082 void OniaPhotonConversionProducer::produce(edm::Event& event, const edm::EventSetup& esetup) {
0083   std::unique_ptr<reco::ConversionCollection> outCollection(new reco::ConversionCollection);
0084   std::unique_ptr<pat::CompositeCandidateCollection> patoutCollection(new pat::CompositeCandidateCollection);
0085
0086   edm::Handle<reco::VertexCollection> priVtxs;
0087   event.getByToken(thePVsToken_, priVtxs);
0088
0089   edm::Handle<reco::ConversionCollection> pConv;
0090   event.getByToken(convCollectionToken_, pConv);
0091
0092   edm::Handle<reco::PFCandidateCollection> pfcandidates;
0093   event.getByToken(pfCandidateCollectionToken_, pfcandidates);
0094
0095   const reco::PFCandidateCollection pfphotons = selectPFPhotons(*pfcandidates);
0096
0097   for (reco::ConversionCollection::const_iterator conv = pConv->begin(); conv != pConv->end(); ++conv) {
0098     if (!(*convSelection_)(*conv)) {
0099       continue;  // selection string
0100     }
0101     if (convAlgo_ != 0 && conv->algo() != convAlgo_) {
0102       continue;  // select algorithm
0103     }
0104     if (!convQuality_.empty()) {
0105       bool flagsok = true;
0106       for (std::vector<int>::const_iterator flag = convQuality_.begin(); flag != convQuality_.end(); ++flag) {
0107         reco::Conversion::ConversionQuality q = (reco::Conversion::ConversionQuality)(*flag);
0108         if (!conv->quality(q)) {
0109           flagsok = false;
0110           break;
0111         }
0112       }
0113       if (!flagsok) {
0114         continue;
0115       }
0116     }
0117     outCollection->push_back(*conv);
0118   }
0119
0120   removeDuplicates(*outCollection);
0121
0122   for (reco::ConversionCollection::const_iterator conv = outCollection->begin(); conv != outCollection->end(); ++conv) {
0123     bool flag1 = true;
0124     bool flag2 = true;
0125     bool flag3 = true;
0126     bool flag4 = true;
0127
0128     // The logic implies that by default this flags are true and if the check are not wanted conversions are saved.
0129     // If checks are required and failed then don't save the conversion.
0130
0131     bool flagTkVtxCompatibility = true;
0132     if (!checkTkVtxCompatibility(*conv, *priVtxs.product())) {
0133       flagTkVtxCompatibility = false;
0134       if (wantTkVtxCompatibility_) {
0135         flag1 = false;
0136       }
0137     }
0138     bool flagCompatibleInnerHits = false;
0139     if (conv->tracks().size() == 2) {
0140       reco::HitPattern hitPatA = conv->tracks().at(0)->hitPattern();
0141       reco::HitPattern hitPatB = conv->tracks().at(1)->hitPattern();
0142       if (foundCompatibleInnerHits(hitPatA, hitPatB) && foundCompatibleInnerHits(hitPatB, hitPatA))
0143         flagCompatibleInnerHits = true;
0144     }
0145     if (wantCompatibleInnerHits_ && !flagCompatibleInnerHits) {
0146       flag2 = false;
0147     }
0148     bool flagHighpurity = true;
0149     if (!HighpuritySubset(*conv, *priVtxs.product())) {
0150       flagHighpurity = false;
0151       if (wantHighpurity_) {
0152         flag3 = false;
0153       }
0154     }
0155     bool pizero_rejected = false;
0156     bool large_pizero_window = CheckPi0(*conv, pfphotons, pizero_rejected);
0157     if (pi0OnlineSwitch_ && pizero_rejected) {
0158       flag4 = false;
0159     }
0160
0161     int flags = 0;
0162     if (flag1 && flag2 && flag3 && flag4) {
0163       flags = PackFlags(
0164           *conv, flagTkVtxCompatibility, flagCompatibleInnerHits, flagHighpurity, pizero_rejected, large_pizero_window);
0165       std::unique_ptr<pat::CompositeCandidate> pat_conv(makePhotonCandidate(*conv));
0166       pat_conv->addUserInt("flags", flags);
0167       patoutCollection->push_back(*pat_conv);
0168     }
0169   }
0170   event.put(std::move(patoutCollection), "conversions");
0171 }
0172
0173 int OniaPhotonConversionProducer::PackFlags(const reco::Conversion& conv,
0174                                             bool flagTkVtxCompatibility,
0175                                             bool flagCompatibleInnerHits,
0176                                             bool flagHighpurity,
0177                                             bool pizero_rejected,
0178                                             bool large_pizero_window) {
0179   int flags = 0;
0180   if (flagTkVtxCompatibility)
0181     flags += 1;
0182   if (flagCompatibleInnerHits)
0183     flags += 2;
0184   if (flagHighpurity)
0185     flags += 4;
0186   if (pizero_rejected)
0187     flags += 8;
0188   if (large_pizero_window)
0189     flags += 16;
0190
0191   flags += (conv.algo() * 32);
0192   int q_mask = 0;
0193   std::vector<std::string> s_quals;
0194   s_quals.push_back("generalTracksOnly");
0195   s_quals.push_back("arbitratedEcalSeeded");
0196   s_quals.push_back("arbitratedMerged");
0197   s_quals.push_back("arbitratedMergedEcalGeneral");
0198   s_quals.push_back("highPurity");
0199   s_quals.push_back("highEfficiency");
0200   s_quals.push_back("ecalMatched1Track");
0201   s_quals.push_back("ecalMatched2Track");
0202   std::vector<int> i_quals = StringToEnumValue<reco::Conversion::ConversionQuality>(s_quals);
0203   for (std::vector<int>::const_iterator qq = i_quals.begin(); qq != i_quals.end(); ++qq) {
0204     reco::Conversion::ConversionQuality q = (reco::Conversion::ConversionQuality)(*qq);
0205     if (conv.quality(q))
0206       q_mask = *qq;
0207   }
0208   flags += (q_mask * 32 * 8);
0209   return flags;
0210 }
0211
0212 /** Put in out collection only those conversion candidates that are not sharing tracks.
0213     If sharing, keep the one with the best chi2.
0214  */
0215 void OniaPhotonConversionProducer::removeDuplicates(reco::ConversionCollection& c) {
0216   // first sort from high to low chi2 prob
0217   std::sort(c.begin(), c.end(), ConversionLessByChi2);
0218   int iter1 = 0;
0219   // Cycle over all the elements of the collection and compare to all the following,
0220   // if two elements have at least one track in common delete the element with the lower chi2
0221   while (iter1 < (((int)c.size()) - 1)) {
0222     int iter2 = iter1 + 1;
0223     while (iter2 < (int)c.size())
0224       if (ConversionEqualByTrack(c[iter1], c[iter2])) {
0225         c.erase(c.begin() + iter2);
0226       } else {
0227         iter2++;  // Increment index only if this element is no duplicate.
0228             // If it is duplicate check again the same index since the vector rearranged elements index after erasing
0229       }
0230     iter1++;
0231   }
0232 }
0233
0234 bool OniaPhotonConversionProducer::checkTkVtxCompatibility(const reco::Conversion& conv,
0235                                                            const reco::VertexCollection& priVtxs) {
0236   std::vector<std::pair<double, short>> idx[2];
0237   short ik = -1;
0238   for (auto const& tk : conv.tracks()) {
0239     ik++;
0240     short count = -1;
0241     for (auto const& vtx : priVtxs) {
0242       count++;
0243
0244       double dz_ = tk->dz(vtx.position());
0245       double dzError_ = tk->dzError();
0246       dzError_ = sqrt(dzError_ * dzError_ + vtx.covariance(2, 2));
0247       if (fabs(dz_) / dzError_ > sigmaTkVtxComp_)
0248         continue;
0249       idx[ik].push_back(std::pair<double, short>(fabs(dz_), count));
0250     }
0251     if (idx[ik].empty())
0252       return false;
0253     if (idx[ik].size() > 1)
0254       std::stable_sort(idx[ik].begin(), idx[ik].end(), lt_);
0255   }
0256   if (ik != 1)
0257     return false;
0258   if (idx[0][0].second == idx[1][0].second)
0259     return true;
0260   if (idx[0].size() > 1 && idx[0][1].second == idx[1][0].second)
0261     return true;
0262   if (idx[1].size() > 1 && idx[0][0].second == idx[1][1].second)
0263     return true;
0264
0265   return false;
0266 }
0267
0268 bool OniaPhotonConversionProducer::foundCompatibleInnerHits(const reco::HitPattern& hitPatA,
0269                                                             const reco::HitPattern& hitPatB) {
0270   size_t count = 0;
0271   uint32_t oldSubStr = 0;
0272   for (int i = 0; i < hitPatA.numberOfAllHits(reco::HitPattern::HitCategory::TRACK_HITS) && count < 2; i++) {
0273     uint32_t hitA = hitPatA.getHitPattern(reco::HitPattern::HitCategory::TRACK_HITS, i);
0274     if (!hitPatA.validHitFilter(hitA) || !hitPatA.trackerHitFilter(hitA))
0275       continue;
0276
0277     if (hitPatA.getSubStructure(hitA) == oldSubStr && hitPatA.getLayer(hitA) == oldSubStr)
0278       continue;
0279
0280     if (hitPatB.getTrackerMonoStereo(
0281             reco::HitPattern::HitCategory::TRACK_HITS, hitPatA.getSubStructure(hitA), hitPatA.getLayer(hitA)) != 0)
0282       return true;
0283
0284     oldSubStr = hitPatA.getSubStructure(hitA);
0285     count++;
0286   }
0287   return false;
0288 }
0289
0290 bool OniaPhotonConversionProducer::HighpuritySubset(const reco::Conversion& conv,
0291                                                     const reco::VertexCollection& priVtxs) {
0292   // select high purity conversions our way:
0293   // vertex chi2 cut
0294   if (ChiSquaredProbability(conv.conversionVertex().chi2(), conv.conversionVertex().ndof()) < _vertexChi2ProbCut)
0295     return false;
0296
0297   // d0 cut
0298   // Find closest primary vertex
0299   int closest_pv_index = 0;
0300   int i = 0;
0301   for (auto const& vtx : priVtxs) {
0302     if (conv.zOfPrimaryVertexFromTracks(vtx.position()) <
0303         conv.zOfPrimaryVertexFromTracks(priVtxs[closest_pv_index].position()))
0304       closest_pv_index = i;
0305     i++;
0306   }
0307   // Now check impact parameter wtr with the just found closest primary vertex
0308   for (auto const& tk : conv.tracks())
0309     if (-tk->dxy(priVtxs[closest_pv_index].position()) * tk->charge() / tk->dxyError() < 0)
0310       return false;
0311
0312   // chi2 of single tracks
0313   for (auto const& tk : conv.tracks())
0314     if (tk->normalizedChi2() > _trackchi2Cut)
0315       return false;
0316
0317   // dof for each track
0318   for (auto const& tk : conv.tracks())
0319     if (tk->ndof() < TkMinNumOfDOF_)
0320       return false;
0321
0322   // distance of approach cut
0323   if (conv.distOfMinimumApproach() < _minDistanceOfApproachMinCut ||
0324       conv.distOfMinimumApproach() > _minDistanceOfApproachMaxCut)
0325     return false;
0326
0327   return true;
0328 }
0329
0330 pat::CompositeCandidate* OniaPhotonConversionProducer::makePhotonCandidate(const reco::Conversion& conv) {
0331   pat::CompositeCandidate* photonCand = new pat::CompositeCandidate();
0332   photonCand->setP4(convertVector(conv.refittedPair4Momentum()));
0333   photonCand->setVertex(conv.conversionVertex().position());
0334
0335   photonCand->addUserData<reco::Track>("track0", *conv.tracks()[0]);
0336   photonCand->addUserData<reco::Track>("track1", *conv.tracks()[1]);
0337
0338   return photonCand;
0339 }
0340
0341 // create a collection of PF photons
0342 const reco::PFCandidateCollection OniaPhotonConversionProducer::selectPFPhotons(
0343     const reco::PFCandidateCollection& pfcandidates) {
0344   reco::PFCandidateCollection pfphotons;
0345   for (reco::PFCandidateCollection::const_iterator cand = pfcandidates.begin(); cand != pfcandidates.end(); ++cand) {
0346     if (cand->particleId() == reco::PFCandidate::gamma)
0347       pfphotons.push_back(*cand);
0348   }
0349   return pfphotons;
0350 }
0351
0352 bool OniaPhotonConversionProducer::CheckPi0(const reco::Conversion& conv,
0353                                             const reco::PFCandidateCollection& photons,
0354                                             bool& pizero_rejected) {
0355   // 2 windows are defined for Pi0 rejection, Conversions that, paired with others photons from the event, have an
0356   // invariant mass inside the "small" window will be pizero_rejected and those that falls in the large window will
0357   // be CheckPi0.
0358   bool check_small = false;
0359   bool check_large = false;
0360
0361   float small1 = pi0SmallWindow_[0];
0362   float small2 = pi0SmallWindow_[1];
0363   float large1 = pi0LargeWindow_[0];
0364   float large2 = pi0LargeWindow_[1];
0365   for (reco::PFCandidateCollection::const_iterator photon = photons.begin(); photon != photons.end(); ++photon) {
0366     float inv = (conv.refittedPair4Momentum() + photon->p4()).M();
0367     if (inv > large1 && inv < large2) {
0368       check_large = true;
0369       if (inv > small1 && inv < small2) {
0370         check_small = true;
0371         break;
0372       }
0373     }
0374   }
0375   pizero_rejected = check_small;
0376   return check_large;
0377 }
0378
0379 reco::Candidate::LorentzVector OniaPhotonConversionProducer::convertVector(const math::XYZTLorentzVectorF& v) {
0380   return reco::Candidate::LorentzVector(v.x(), v.y(), v.z(), v.t());
0381 }
0382
0383 void OniaPhotonConversionProducer::endStream() {}
0384 //define this as a plug-in
0385 DEFINE_FWK_MODULE(OniaPhotonConversionProducer);