File indexing completed on 2023-03-17 11:22:07
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 #include <memory>
0020 #include "DataFormats/Math/interface/deltaPhi.h"
0021 #include "DataFormats/TrackReco/interface/Track.h"
0022 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0023 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/Frameworkfwd.h"
0026 #include "FWCore/Framework/interface/MakerMacros.h"
0027 #include "FWCore/Framework/interface/stream/EDProducer.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0030 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0031 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0032 #include "MagneticField/Engine/interface/MagneticField.h"
0033 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0034 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0035 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0036 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0037 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0038 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0039
0040 class ConversionSeedFilter : public edm::stream::EDProducer<> {
0041 public:
0042 explicit ConversionSeedFilter(const edm::ParameterSet&);
0043 ~ConversionSeedFilter() override = default;
0044
0045 private:
0046 void produce(edm::Event&, const edm::EventSetup&) override;
0047 bool isCompatible(double* vars1, double* vars2);
0048 void getKine(const TrajectoryStateOnSurface& tsos, double* vars);
0049 void SearchAmongSeeds(const TrajectorySeedCollection* pInPos,
0050 const TrajectorySeedCollection* pInNeg,
0051 TrajectorySeedCollection& selectedColl,
0052 std::vector<bool>& idxPosColl1,
0053 std::vector<bool>& idxPosColl2);
0054 void SearchAmongTracks(const TrajectorySeedCollection* pInSeed,
0055 const reco::TrackCollection* pInTk,
0056 TrajectorySeedCollection& selectedColl,
0057 std::vector<bool>& idxPosColl);
0058 void SearchAmongTrajectories(const TrajectorySeedCollection* pInSeed,
0059 const Trajectory* InTj,
0060 TrajectorySeedCollection& selectedColl,
0061 std::vector<bool>& idxPosColl);
0062
0063 TrajectoryStateOnSurface getTSOS(const TrajectorySeed& ts);
0064 TrajectoryStateOnSurface getTSOS(const reco::Track& tk);
0065 TrajectoryStateOnSurface getTSOS(const Trajectory& tj, const TrajectorySeed& ts);
0066
0067 const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken;
0068 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> mfToken;
0069 const edm::EDGetTokenT<TrajectorySeedCollection> inputCollSeedPos;
0070 const edm::EDGetTokenT<TrajectorySeedCollection> inputCollSeedNeg;
0071 const edm::EDGetTokenT<TrajTrackAssociationCollection> inputTrajectory;
0072 const double deltaPhiCut, deltaCotThetaCut, deltaRCut, deltaZCut;
0073
0074 const TrackerGeometry* theG;
0075 const MagneticField* theMF;
0076 uint32_t maxInputSeeds;
0077 bool takeAll;
0078 };
0079
0080 ConversionSeedFilter::ConversionSeedFilter(const edm::ParameterSet& cfg)
0081 : geomToken(esConsumes()),
0082 mfToken(esConsumes()),
0083
0084 inputCollSeedPos(consumes<TrajectorySeedCollection>(cfg.getParameter<edm::InputTag>("seedCollectionPos"))),
0085
0086 inputCollSeedNeg(consumes<TrajectorySeedCollection>(cfg.getParameter<edm::InputTag>("seedCollectionNeg"))),
0087 inputTrajectory(consumes<TrajTrackAssociationCollection>(cfg.getParameter<edm::InputTag>("inputTrajectory"))),
0088 deltaPhiCut(cfg.getParameter<double>("deltaPhiCut")),
0089 deltaCotThetaCut(cfg.getParameter<double>("deltaCotThetaCut")),
0090 deltaRCut(cfg.getParameter<double>("deltaRCut")),
0091 deltaZCut(cfg.getParameter<double>("deltaZCut")),
0092 maxInputSeeds(cfg.getParameter<uint32_t>("maxInputSeeds")),
0093 takeAll(cfg.getParameter<bool>("takeAll")) {
0094 produces<TrajectorySeedCollection>();
0095 }
0096
0097 void ConversionSeedFilter::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0098 using namespace edm;
0099 using namespace std;
0100 Handle<TrajectorySeedCollection> pInPos;
0101 iEvent.getByToken(inputCollSeedPos, pInPos);
0102 Handle<TrajectorySeedCollection> pInNeg;
0103 iEvent.getByToken(inputCollSeedNeg, pInNeg);
0104
0105
0106
0107
0108 edm::Handle<TrajTrackAssociationCollection> trajTrackAssociations;
0109 iEvent.getByToken(inputTrajectory, trajTrackAssociations);
0110
0111 theG = &iSetup.getData(geomToken);
0112 theMF = &iSetup.getData(mfToken);
0113
0114 auto result = std::make_unique<TrajectorySeedCollection>();
0115
0116 TrajectorySeedCollection selectedColl;
0117
0118 if (takeAll) {
0119 result->insert(result->end(), pInPos->begin(), pInPos->end());
0120 result->insert(result->end(), pInNeg->begin(), pInNeg->end());
0121 } else {
0122 edm::LogInfo("ConversionSeedFilter") << "takeAll " << takeAll;
0123 if (pInPos->size() < maxInputSeeds && pInNeg->size() < maxInputSeeds) {
0124 std::vector<bool> idxPosColl1(pInPos->size(), false);
0125 std::vector<bool> idxPosColl2(pInNeg->size(), false);
0126 selectedColl.reserve(pInPos->size());
0127
0128 edm::LogInfo("ConversionSeedFilter")
0129 << "New Event \t Pos " << pInPos->size() << " \t Neg " << pInNeg->size() << std::endl;
0130
0131 SearchAmongSeeds(pInPos.product(), pInNeg.product(), selectedColl, idxPosColl1, idxPosColl2);
0132
0133
0134
0135 if (trajTrackAssociations.isValid()) {
0136 edm::LogInfo("ConversionSeedFilter") << "Reconstructed tracks " << trajTrackAssociations->size() << std::endl;
0137 for (TrajTrackAssociationCollection::const_iterator association = trajTrackAssociations->begin();
0138 association != trajTrackAssociations->end();
0139 association++) {
0140 const Trajectory* traj = association->key.get();
0141 const reco::Track* track = association->val.get();
0142
0143
0144
0145 if (track->charge() < 0) {
0146 SearchAmongTrajectories(pInPos.product(), traj, selectedColl, idxPosColl1);
0147 } else {
0148 SearchAmongTrajectories(pInNeg.product(), traj, selectedColl, idxPosColl2);
0149 }
0150 }
0151 }
0152 }
0153 result->insert(result->end(), selectedColl.begin(), selectedColl.end());
0154 }
0155
0156 edm::LogInfo("ConversionSeedFilter") << "\nNew Event : result size " << result->size() << std::endl;
0157
0158 iEvent.put(std::move(result));
0159 }
0160
0161 void ConversionSeedFilter::SearchAmongSeeds(const TrajectorySeedCollection* pInPos,
0162 const TrajectorySeedCollection* pInNeg,
0163 TrajectorySeedCollection& selectedColl,
0164 std::vector<bool>& idxPosColl1,
0165 std::vector<bool>& idxPosColl2) {
0166 for (TrajectorySeedCollection::const_iterator iS1 = pInPos->begin(); iS1 != pInPos->end(); ++iS1) {
0167 bool pushed1 = false;
0168
0169 double vars1[4];
0170 getKine(getTSOS(*iS1), vars1);
0171
0172 for (TrajectorySeedCollection::const_iterator iS2 = pInNeg->begin(); iS2 != pInNeg->end(); ++iS2) {
0173 double vars2[4];
0174 getKine(getTSOS(*iS2), vars2);
0175
0176 if (isCompatible(vars1, vars2)) {
0177 edm::LogInfo("ConversionSeedFilter")
0178 << "[SearchAmongSeeds] match in pos " << iS1 - pInPos->begin() << " " << iS2 - pInNeg->begin() << std::endl;
0179 if (!pushed1) {
0180 idxPosColl1[iS1 - pInPos->begin()] = true;
0181 selectedColl.push_back(*iS1);
0182 pushed1 = true;
0183 }
0184 if (!idxPosColl2[iS2 - pInNeg->begin()]) {
0185 selectedColl.push_back(*iS2);
0186 idxPosColl2[iS2 - pInNeg->begin()] = true;
0187 }
0188 }
0189 }
0190 }
0191 }
0192
0193 void ConversionSeedFilter::SearchAmongTracks(const TrajectorySeedCollection* pInSeed,
0194 const reco::TrackCollection* pInTk,
0195 TrajectorySeedCollection& selectedColl,
0196 std::vector<bool>& idxPosColl) {
0197 for (TrajectorySeedCollection::const_iterator iS1 = pInSeed->begin(); iS1 != pInSeed->end(); ++iS1) {
0198 if (idxPosColl[iS1 - pInSeed->begin()])
0199 continue;
0200
0201 double vars1[4];
0202 getKine(getTSOS(*iS1), vars1);
0203
0204 for (reco::TrackCollection::const_iterator iS2 = pInTk->begin(); iS2 != pInTk->end(); ++iS2) {
0205 double vars2[4];
0206 getKine(getTSOS(*iS2), vars2);
0207
0208 if (isCompatible(vars1, vars2)) {
0209 edm::LogInfo("ConversionSeedFilter")
0210 << "[SearchAmongTracks] match in pos " << iS1 - pInSeed->begin() << std::endl;
0211 idxPosColl[iS1 - pInSeed->begin()] = true;
0212 selectedColl.push_back(*iS1);
0213 continue;
0214 }
0215 }
0216 }
0217 }
0218
0219 void ConversionSeedFilter::SearchAmongTrajectories(const TrajectorySeedCollection* pInSeed,
0220 const Trajectory* InTj,
0221 TrajectorySeedCollection& selectedColl,
0222 std::vector<bool>& idxPosColl) {
0223 for (TrajectorySeedCollection::const_iterator iS1 = pInSeed->begin(); iS1 != pInSeed->end(); ++iS1) {
0224 if (idxPosColl[iS1 - pInSeed->begin()])
0225 continue;
0226
0227 double vars1[4];
0228 getKine(getTSOS(*iS1), vars1);
0229
0230 double vars2[4];
0231 getKine(getTSOS(*InTj, *iS1), vars2);
0232
0233 if (isCompatible(vars1, vars2)) {
0234 edm::LogInfo("ConversionSeedFilter")
0235 << "[SearchAmongTrajectories] match seed in pos " << iS1 - pInSeed->begin() << " of " << pInSeed->size()
0236 << " seed charge " << iS1->startingState().parameters().charge() << std::endl;
0237 idxPosColl[iS1 - pInSeed->begin()] = true;
0238 selectedColl.push_back(*iS1);
0239 }
0240 }
0241 }
0242
0243 TrajectoryStateOnSurface ConversionSeedFilter::getTSOS(const TrajectorySeed& ts) {
0244 PTrajectoryStateOnDet state = ts.startingState();
0245 DetId detId(state.detId());
0246 return trajectoryStateTransform::transientState(state, &(theG->idToDet(detId)->surface()), theMF);
0247 }
0248
0249 TrajectoryStateOnSurface ConversionSeedFilter::getTSOS(const reco::Track& tk) {
0250 return trajectoryStateTransform::innerStateOnSurface(tk, *theG, theMF);
0251 }
0252
0253 TrajectoryStateOnSurface ConversionSeedFilter::getTSOS(const Trajectory& tj, const TrajectorySeed& ts) {
0254 PTrajectoryStateOnDet state = ts.startingState();
0255 DetId detId(state.detId());
0256 GlobalPoint p = theG->idToDet(detId)->surface().toGlobal(state.parameters().position());
0257 return tj.closestMeasurement(p).updatedState();
0258 }
0259
0260 bool ConversionSeedFilter::isCompatible(double* vars1, double* vars2) {
0261 double deltaPhi = fabs(reco::deltaPhi(vars1[0], vars2[0]));
0262 double deltaCotTheta = fabs(vars1[1] - vars2[1]);
0263 double deltaR = fabs(vars1[2] - vars2[2]);
0264 double deltaZ = fabs(vars1[3] - vars2[3]);
0265
0266
0267
0268
0269
0270
0271
0272
0273 return deltaPhi < deltaPhiCut && deltaCotTheta < deltaCotThetaCut && deltaR < deltaRCut && deltaZ < deltaZCut;
0274 }
0275
0276 void ConversionSeedFilter::getKine(const TrajectoryStateOnSurface& tsos, double* vars) {
0277 vars[0] = tsos.globalMomentum().phi();
0278 vars[1] = 2 / tan(tsos.globalMomentum().theta());
0279 vars[2] = tsos.globalPosition().perp();
0280 vars[3] = tsos.globalPosition().z();
0281 }
0282
0283 DEFINE_FWK_MODULE(ConversionSeedFilter);