File indexing completed on 2023-05-08 23:19:13
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0012 #include "DataFormats/Common/interface/Handle.h"
0013 #include "DataFormats/EgammaTrackReco/interface/ConversionTrack.h"
0014 #include "DataFormats/EgammaTrackReco/interface/ConversionTrackFwd.h"
0015 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
0016 #include "DataFormats/TrackReco/interface/Track.h"
0017 #include "DataFormats/TrackReco/interface/TrackBase.h"
0018 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0019 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0020 #include "FWCore/Framework/interface/ESHandle.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/Framework/interface/EventSetup.h"
0023 #include "FWCore/Framework/interface/stream/EDProducer.h"
0024 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0025 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0026 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0027 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0028 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0029
0030 #include <vector>
0031
0032 class ConversionTrackMerger : public edm::stream::EDProducer<> {
0033 public:
0034 explicit ConversionTrackMerger(const edm::ParameterSet& conf);
0035
0036 ~ConversionTrackMerger() override;
0037
0038 void produce(edm::Event& e, const edm::EventSetup& c) override;
0039
0040 private:
0041 edm::ParameterSet conf_;
0042
0043 edm::EDGetTokenT<reco::ConversionTrackCollection> trackProducer1;
0044 edm::EDGetTokenT<reco::ConversionTrackCollection> trackProducer2;
0045
0046 std::unique_ptr<reco::ConversionTrackCollection> outputTrks;
0047 };
0048
0049 #include "FWCore/Framework/interface/MakerMacros.h"
0050 DEFINE_FWK_MODULE(ConversionTrackMerger);
0051
0052 ConversionTrackMerger::ConversionTrackMerger(edm::ParameterSet const& conf) : conf_(conf) {
0053
0054 trackProducer1 = consumes<reco::ConversionTrackCollection>(conf_.getParameter<edm::InputTag>("TrackProducer1"));
0055 trackProducer2 = consumes<reco::ConversionTrackCollection>(conf_.getParameter<edm::InputTag>("TrackProducer2"));
0056
0057 produces<reco::ConversionTrackCollection>();
0058 }
0059
0060
0061 ConversionTrackMerger::~ConversionTrackMerger() {}
0062
0063
0064 void ConversionTrackMerger::produce(edm::Event& e, const edm::EventSetup& es) {
0065 double shareFrac = conf_.getParameter<double>("ShareFrac");
0066 bool allowFirstHitShare = conf_.getParameter<bool>("allowFirstHitShare");
0067 bool checkCharge = conf_.getParameter<bool>("checkCharge");
0068 double minProb = conf_.getParameter<double>("minProb");
0069
0070 int outputPreferCollection = conf_.getParameter<int>("outputPreferCollection");
0071 int trackerOnlyPreferCollection = conf_.getParameter<int>("trackerOnlyPreferCollection");
0072 int arbitratedEcalSeededPreferCollection = conf_.getParameter<int>("arbitratedEcalSeededPreferCollection");
0073 int arbitratedMergedPreferCollection = conf_.getParameter<int>("arbitratedMergedPreferCollection");
0074 int arbitratedMergedEcalGeneralPreferCollection =
0075 conf_.getParameter<int>("arbitratedMergedEcalGeneralPreferCollection");
0076
0077
0078
0079
0080
0081
0082 const reco::ConversionTrackCollection* TC1 = nullptr;
0083 static const reco::ConversionTrackCollection s_empty1, s_empty2;
0084 edm::Handle<reco::ConversionTrackCollection> trackCollection1;
0085 e.getByToken(trackProducer1, trackCollection1);
0086 if (trackCollection1.isValid()) {
0087 TC1 = trackCollection1.product();
0088
0089 } else {
0090 TC1 = &s_empty1;
0091 edm::LogWarning("ConversionTrackMerger") << "1st TrackCollection not found;"
0092 << " will only clean 2nd TrackCollection ";
0093 }
0094 reco::ConversionTrackCollection tC1 = *TC1;
0095
0096 const reco::ConversionTrackCollection* TC2 = nullptr;
0097 edm::Handle<reco::ConversionTrackCollection> trackCollection2;
0098 e.getByToken(trackProducer2, trackCollection2);
0099 if (trackCollection2.isValid()) {
0100 TC2 = trackCollection2.product();
0101
0102 } else {
0103 TC2 = &s_empty2;
0104 edm::LogWarning("ConversionTrackMerger") << "2nd TrackCollection not found;"
0105 << " will only clean 1st TrackCollection ";
0106 }
0107 reco::ConversionTrackCollection tC2 = *TC2;
0108
0109
0110 outputTrks = std::make_unique<reco::ConversionTrackCollection>();
0111 int i;
0112
0113 std::vector<int> selected1;
0114 for (unsigned int i = 0; i < tC1.size(); ++i) {
0115 selected1.push_back(1);
0116 }
0117 std::vector<int> selected2;
0118 for (unsigned int i = 0; i < tC2.size(); ++i) {
0119 selected2.push_back(1);
0120 }
0121
0122 std::map<reco::ConversionTrackCollection::const_iterator, std::vector<const TrackingRecHit*> > rh1;
0123 std::map<reco::ConversionTrackCollection::const_iterator, std::vector<const TrackingRecHit*> > rh2;
0124 for (reco::ConversionTrackCollection::const_iterator track = tC1.begin(); track != tC1.end(); ++track) {
0125 trackingRecHit_iterator itB = track->track()->recHitsBegin();
0126 trackingRecHit_iterator itE = track->track()->recHitsEnd();
0127 for (trackingRecHit_iterator it = itB; it != itE; ++it) {
0128 const TrackingRecHit* hit = &(**it);
0129 rh1[track].push_back(hit);
0130 }
0131 }
0132 for (reco::ConversionTrackCollection::const_iterator track = tC2.begin(); track != tC2.end(); ++track) {
0133 trackingRecHit_iterator jtB = track->track()->recHitsBegin();
0134 trackingRecHit_iterator jtE = track->track()->recHitsEnd();
0135 for (trackingRecHit_iterator jt = jtB; jt != jtE; ++jt) {
0136 const TrackingRecHit* hit = &(**jt);
0137 rh2[track].push_back(hit);
0138 }
0139 }
0140
0141 if ((!tC1.empty()) && (!tC2.empty())) {
0142 i = -1;
0143 for (reco::ConversionTrackCollection::iterator track = tC1.begin(); track != tC1.end(); ++track) {
0144 i++;
0145
0146
0147 selected1[i] = selected1[i] && outputPreferCollection != 0;
0148 track->setIsTrackerOnly(track->isTrackerOnly() && trackerOnlyPreferCollection != 0);
0149 track->setIsArbitratedEcalSeeded(track->isArbitratedEcalSeeded() && arbitratedEcalSeededPreferCollection != 0);
0150 track->setIsArbitratedMerged(track->isArbitratedMerged() && arbitratedMergedPreferCollection != 0);
0151 track->setIsArbitratedMergedEcalGeneral(track->isArbitratedMergedEcalGeneral() &&
0152 arbitratedMergedEcalGeneralPreferCollection != 0);
0153
0154 std::vector<const TrackingRecHit*>& iHits = rh1[track];
0155 unsigned nh1 = iHits.size();
0156 int j = -1;
0157 for (reco::ConversionTrackCollection::iterator track2 = tC2.begin(); track2 != tC2.end(); ++track2) {
0158 j++;
0159
0160
0161 selected2[j] = selected2[j] && outputPreferCollection != 0;
0162 track2->setIsTrackerOnly(track2->isTrackerOnly() && trackerOnlyPreferCollection != 0);
0163 track2->setIsArbitratedEcalSeeded(track2->isArbitratedEcalSeeded() &&
0164 arbitratedEcalSeededPreferCollection != 0);
0165 track2->setIsArbitratedMerged(track2->isArbitratedMerged() && arbitratedMergedPreferCollection != 0);
0166 track2->setIsArbitratedMergedEcalGeneral(track2->isArbitratedMergedEcalGeneral() &&
0167 arbitratedMergedEcalGeneralPreferCollection != 0);
0168
0169 std::vector<const TrackingRecHit*>& jHits = rh2[track2];
0170 unsigned nh2 = jHits.size();
0171 int noverlap = 0;
0172 int firstoverlap = 0;
0173 for (unsigned ih = 0; ih < nh1; ++ih) {
0174 const TrackingRecHit* it = iHits[ih];
0175 if (it->isValid()) {
0176 for (unsigned jh = 0; jh < nh2; ++jh) {
0177 const TrackingRecHit* jt = jHits[jh];
0178 if (jt->isValid()) {
0179 if (it->sharesInput(jt, TrackingRecHit::some)) {
0180 noverlap++;
0181 if (allowFirstHitShare && (ih == 0) && (jh == 0))
0182 firstoverlap = 1;
0183 }
0184 }
0185 }
0186 }
0187 }
0188 int nhit1 = track->track()->numberOfValidHits();
0189 int nhit2 = track2->track()->numberOfValidHits();
0190
0191
0192 if ((noverlap - firstoverlap) > (std::min(nhit1, nhit2) - firstoverlap) * shareFrac &&
0193 (!checkCharge || track->track()->charge() * track2->track()->charge() > 0)) {
0194
0195
0196
0197 double probFirst = ChiSquaredProbability(track->track()->chi2(), track->track()->ndof());
0198 double probSecond = ChiSquaredProbability(track2->track()->chi2(), track2->track()->ndof());
0199
0200
0201 bool keepFirst = (nhit1 > nhit2 ||
0202 (nhit1 == nhit2 && track->track()->normalizedChi2() < track2->track()->normalizedChi2()));
0203
0204
0205 keepFirst |= (probFirst > minProb && probSecond <= minProb);
0206 keepFirst &= !(probFirst <= minProb && probSecond > minProb);
0207
0208 bool keepSecond = !keepFirst;
0209
0210
0211
0212 selected1[i] = selected1[i] && ((keepFirst && outputPreferCollection == 3) || outputPreferCollection == -1 ||
0213 outputPreferCollection == 1);
0214 track->setIsTrackerOnly(track->isTrackerOnly() &&
0215 ((keepFirst && trackerOnlyPreferCollection == 3) ||
0216 trackerOnlyPreferCollection == -1 || trackerOnlyPreferCollection == 1));
0217 track->setIsArbitratedEcalSeeded(track->isArbitratedEcalSeeded() &&
0218 ((keepFirst && arbitratedEcalSeededPreferCollection == 3) ||
0219 arbitratedEcalSeededPreferCollection == -1 ||
0220 arbitratedEcalSeededPreferCollection == 1));
0221 track->setIsArbitratedMerged(track->isArbitratedMerged() &&
0222 ((keepFirst && arbitratedMergedPreferCollection == 3) ||
0223 arbitratedMergedPreferCollection == -1 ||
0224 arbitratedMergedPreferCollection == 1));
0225 track->setIsArbitratedMergedEcalGeneral(track->isArbitratedMergedEcalGeneral() &&
0226 ((keepFirst && arbitratedMergedEcalGeneralPreferCollection == 3) ||
0227 arbitratedMergedEcalGeneralPreferCollection == -1 ||
0228 arbitratedMergedEcalGeneralPreferCollection == 1));
0229
0230 selected2[j] = selected2[j] && ((keepSecond && outputPreferCollection == 3) || outputPreferCollection == -1 ||
0231 outputPreferCollection == 2);
0232 track2->setIsTrackerOnly(track2->isTrackerOnly() &&
0233 ((keepSecond && trackerOnlyPreferCollection == 3) ||
0234 trackerOnlyPreferCollection == -1 || trackerOnlyPreferCollection == 2));
0235 track2->setIsArbitratedEcalSeeded(track2->isArbitratedEcalSeeded() &&
0236 ((keepSecond && arbitratedEcalSeededPreferCollection == 3) ||
0237 arbitratedEcalSeededPreferCollection == -1 ||
0238 arbitratedEcalSeededPreferCollection == 2));
0239 track2->setIsArbitratedMerged(track2->isArbitratedMerged() &&
0240 ((keepSecond && arbitratedMergedPreferCollection == 3) ||
0241 arbitratedMergedPreferCollection == -1 ||
0242 arbitratedMergedPreferCollection == 2));
0243 track2->setIsArbitratedMergedEcalGeneral(track2->isArbitratedMergedEcalGeneral() &&
0244 ((keepSecond && arbitratedMergedEcalGeneralPreferCollection == 3) ||
0245 arbitratedMergedEcalGeneralPreferCollection == -1 ||
0246 arbitratedMergedEcalGeneralPreferCollection == 2));
0247
0248 }
0249 }
0250 }
0251 }
0252
0253
0254
0255
0256
0257 if (!tC1.empty()) {
0258 i = 0;
0259 for (reco::ConversionTrackCollection::const_iterator track = tC1.begin(); track != tC1.end(); ++track, ++i) {
0260
0261 if (!selected1[i]) {
0262 continue;
0263 }
0264
0265 outputTrks->push_back(*track);
0266 }
0267 }
0268
0269
0270
0271 if (!tC2.empty()) {
0272 i = 0;
0273 for (reco::ConversionTrackCollection::const_iterator track = tC2.begin(); track != tC2.end(); ++track, ++i) {
0274
0275 if (!selected2[i]) {
0276 continue;
0277 }
0278
0279 outputTrks->push_back(*track);
0280 }
0281 }
0282
0283 e.put(std::move(outputTrks));
0284 return;
0285
0286 }