File indexing completed on 2024-04-06 12:15:33
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
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
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 }
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;
0100 }
0101 if (convAlgo_ != 0 && conv->algo() != convAlgo_) {
0102 continue;
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
0129
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
0213
0214
0215 void OniaPhotonConversionProducer::removeDuplicates(reco::ConversionCollection& c) {
0216
0217 std::sort(c.begin(), c.end(), ConversionLessByChi2);
0218 int iter1 = 0;
0219
0220
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++;
0228
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
0293
0294 if (ChiSquaredProbability(conv.conversionVertex().chi2(), conv.conversionVertex().ndof()) < _vertexChi2ProbCut)
0295 return false;
0296
0297
0298
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
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
0313 for (auto const& tk : conv.tracks())
0314 if (tk->normalizedChi2() > _trackchi2Cut)
0315 return false;
0316
0317
0318 for (auto const& tk : conv.tracks())
0319 if (tk->ndof() < TkMinNumOfDOF_)
0320 return false;
0321
0322
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
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
0356
0357
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
0385 DEFINE_FWK_MODULE(OniaPhotonConversionProducer);