ConversionSeedGenerators/plugins/ConversionSeedFilterCharge.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    ConversionSeedFilterCharge
0004 // Class:      ConversionSeedFilterCharge
0005 //
0006 /**\class ConversionSeedFilterCharge ConversionSeedFilterCharge.cc RecoTracker/TkSeedGenerator/src/ConversionSeedFilterCharge.cc
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Giuseppe Cerati
0015 //         Created:  Thu Mar 11 10:48:48 CET 2010
0016 //
0017 //
0018
0019 #include <memory>
0020 #include "DataFormats/Math/interface/deltaPhi.h"
0021 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0022 #include "FWCore/Framework/interface/Event.h"
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/MakerMacros.h"
0025 #include "FWCore/Framework/interface/global/EDProducer.h"
0026 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0027 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0028 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0029 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0030 #include "MagneticField/Engine/interface/MagneticField.h"
0031 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0032 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0033 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0034 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0035
0036 class ConversionSeedFilterCharge : public edm::global::EDProducer<> {
0037 public:
0038   explicit ConversionSeedFilterCharge(const edm::ParameterSet&);
0039   ~ConversionSeedFilterCharge() override = default;
0040
0041 private:
0042   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0043   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken;
0044   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> mfToken;
0045
0046   edm::EDGetTokenT<TrajectorySeedCollection> inputCollPos;
0047   edm::EDGetTokenT<TrajectorySeedCollection> inputCollNeg;
0048   const double deltaPhiCut, deltaCotThetaCut, deltaRCut, deltaZCut;
0049
0050   const uint32_t maxInputSeeds;
0051 };
0052
0053 ConversionSeedFilterCharge::ConversionSeedFilterCharge(const edm::ParameterSet& cfg)
0054     : geomToken(esConsumes()),
0055       mfToken(esConsumes()),
0056       inputCollPos(consumes<TrajectorySeedCollection>(cfg.getParameter<edm::InputTag>("seedCollectionPos"))),
0057       inputCollNeg(consumes<TrajectorySeedCollection>(cfg.getParameter<edm::InputTag>("seedCollectionNeg"))),
0058       deltaPhiCut(cfg.getParameter<double>("deltaPhiCut")),
0059       deltaCotThetaCut(cfg.getParameter<double>("deltaCotThetaCut")),
0060       deltaRCut(cfg.getParameter<double>("deltaRCut")),
0061       deltaZCut(cfg.getParameter<double>("deltaZCut")),
0062       maxInputSeeds(cfg.getParameter<uint32_t>("maxInputSeeds")) {
0063   produces<TrajectorySeedCollection>();
0064 }
0065
0066 void ConversionSeedFilterCharge::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0067   using namespace edm;
0068   using namespace std;
0069   Handle<TrajectorySeedCollection> pInPos;
0070   iEvent.getByToken(inputCollPos, pInPos);
0071   Handle<TrajectorySeedCollection> pInNeg;
0072   iEvent.getByToken(inputCollNeg, pInNeg);
0073
0074   const TrackerGeometry* theG = &iSetup.getData(geomToken);
0075   const MagneticField* theMF = &iSetup.getData(mfToken);
0076
0077   auto result = std::make_unique<TrajectorySeedCollection>();
0078   result->reserve(pInPos->size());
0079
0080   if (pInPos->size() < maxInputSeeds && pInNeg->size() < maxInputSeeds) {
0081     edm::LogInfo("ConversionSeedFilterCharge")
0082         << "New Event \t Pos " << pInPos->size() << " \t Neg " << pInNeg->size() << std::endl;
0083
0084     std::vector<int> inResult;
0085     for (TrajectorySeedCollection::const_iterator iS1 = pInPos->begin(); iS1 != pInPos->end(); ++iS1) {
0086       PTrajectoryStateOnDet state1 = iS1->startingState();
0087       DetId detId1(state1.detId());
0088       TrajectoryStateOnSurface tsos1 =
0089           trajectoryStateTransform::transientState(state1, &(theG->idToDet(detId1)->surface()), theMF);
0090       double phi1 = tsos1.globalMomentum().phi();
0091       double cotTheta1 = 1 / tan(tsos1.globalMomentum().theta());
0092       double r1 = tsos1.globalPosition().perp();
0093       double z1 = tsos1.globalPosition().z();
0094       //cout << "detId1=" << detId1 << " phi1=" << phi1 << " cotTheta1=" << cotTheta1 << " r1=" << r1 << " z1=" << z1 << endl;
0095
0096       bool pushed = false;
0097       for (TrajectorySeedCollection::const_iterator iS2 = pInNeg->begin(); iS2 != pInNeg->end(); ++iS2) {
0098         PTrajectoryStateOnDet state2 = iS2->startingState();
0099         DetId detId2(state2.detId());
0100         TrajectoryStateOnSurface tsos2 =
0101             trajectoryStateTransform::transientState(state2, &(theG->idToDet(detId2)->surface()), theMF);
0102
0103         double deltaPhi = fabs(reco::deltaPhi(phi1, tsos2.globalMomentum().phi()));
0104         double deltaCotTheta = fabs(cotTheta1 - 1 / tan(tsos2.globalMomentum().theta()));
0105         double deltaR = fabs(r1 - tsos2.globalPosition().perp());
0106         double deltaZ = fabs(z1 - tsos2.globalPosition().z());
0107
0108         //   double phi2 = tsos2.globalMomentum().phi();
0109         //   double cotTheta2 = 2/tan(tsos2.globalMomentum().theta());
0110         //   double r2 = tsos2.globalPosition().perp();
0111         //   double z2 = tsos2.globalPosition().z();
0112         //cout << "j=" << j << " detId2=" << detId2 << " phi2=" << phi2 << " cotTheta2=" << cotTheta2 << " r2=" << r2 << " z2=" << z2 << endl;
0113
0114         if (deltaPhi < deltaPhiCut && deltaCotTheta < deltaCotThetaCut && deltaR < deltaRCut && deltaZ < deltaZCut) {
0115           edm::LogInfo("ConversionSeedFilterCharge") << "[SearchAmongSeeds] match in pos " << iS1 - pInPos->begin()
0116                                                      << " " << iS2 - pInNeg->begin() << std::endl;
0117           //cout << "match" << endl;
0118           if (!pushed) {
0119             result->push_back(*iS1);
0120             pushed = true;
0121           }
0122           if (std::find(inResult.begin(), inResult.end(), iS2 - pInNeg->begin()) == inResult.end()) {
0123             result->push_back(*iS2);
0124             inResult.push_back(iS2 - pInNeg->begin());
0125           }
0126         }
0127       }
0128     }
0129   }
0130
0131   edm::LogInfo("ConversionSeedFilterCharge") << "\nNew Event : result size " << result->size() << std::endl;
0132
0133   //cout << "sizes: pInPos=" << pInPos->size() << " pInNeg=" << pInNeg->size() << " result=" << result->size() << endl;
0134   iEvent.put(std::move(result));
0135 }
0136
0137 DEFINE_FWK_MODULE(ConversionSeedFilterCharge);