File indexing completed on 2021-06-16 03:20:14
0001 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0002 #include "DataFormats/Common/interface/RefToBase.h"
0003 #include "DataFormats/ParticleFlowReco/interface/PFConversion.h"
0004 #include "DataFormats/VertexReco/interface/Vertex.h"
0005 #include "FWCore/Framework/interface/ESHandle.h"
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/Framework/interface/stream/EDProducer.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "MagneticField/Engine/interface/MagneticField.h"
0010 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0011 #include "RecoParticleFlow/PFTracking/interface/PFTrackTransformer.h"
0012 #include "TrackingTools/IPTools/interface/IPTools.h"
0013 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0014 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0015 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0016
0017 class PFConversionProducer : public edm::stream::EDProducer<> {
0018 public:
0019
0020 explicit PFConversionProducer(const edm::ParameterSet&);
0021
0022
0023 ~PFConversionProducer() override;
0024
0025 private:
0026 void beginRun(const edm::Run&, const edm::EventSetup&) override;
0027 void endRun(const edm::Run&, const edm::EventSetup&) override;
0028
0029
0030 void produce(edm::Event&, const edm::EventSetup&) override;
0031
0032
0033 PFTrackTransformer* pfTransformer_;
0034 edm::EDGetTokenT<reco::ConversionCollection> pfConversionContainer_;
0035 edm::EDGetTokenT<reco::VertexCollection> vtx_h;
0036
0037 const edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> transientTrackToken_;
0038 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldToken_;
0039 };
0040
0041 #include "FWCore/Framework/interface/MakerMacros.h"
0042 DEFINE_FWK_MODULE(PFConversionProducer);
0043
0044 typedef std::multimap<unsigned, std::vector<unsigned> > BlockMap;
0045 using namespace std;
0046 using namespace edm;
0047
0048 PFConversionProducer::PFConversionProducer(const ParameterSet& iConfig)
0049 : pfTransformer_(nullptr),
0050 transientTrackToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0051 magneticFieldToken_(esConsumes<edm::Transition::BeginRun>()) {
0052 produces<reco::PFRecTrackCollection>();
0053 produces<reco::PFConversionCollection>();
0054
0055 pfConversionContainer_ = consumes<reco::ConversionCollection>(iConfig.getParameter<InputTag>("conversionCollection"));
0056
0057 vtx_h = consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("PrimaryVertexLabel"));
0058 }
0059
0060 PFConversionProducer::~PFConversionProducer() { delete pfTransformer_; }
0061
0062 void PFConversionProducer::produce(Event& iEvent, const EventSetup& iSetup) {
0063
0064 auto pfConversionColl = std::make_unique<reco::PFConversionCollection>();
0065 auto pfRecTrackColl = std::make_unique<reco::PFRecTrackCollection>();
0066
0067 TransientTrackBuilder const& thebuilder = iSetup.getData(transientTrackToken_);
0068
0069 reco::PFRecTrackRefProd pfTrackRefProd = iEvent.getRefBeforePut<reco::PFRecTrackCollection>();
0070 Handle<reco::ConversionCollection> convCollH;
0071 iEvent.getByToken(pfConversionContainer_, convCollH);
0072
0073 const reco::ConversionCollection& convColl = *(convCollH.product());
0074
0075 Handle<reco::VertexCollection> vertex;
0076 iEvent.getByToken(vtx_h, vertex);
0077
0078 reco::Vertex dummy;
0079 const reco::Vertex* pv = &dummy;
0080 if (vertex.isValid()) {
0081 pv = &*vertex->begin();
0082 } else {
0083 reco::Vertex::Error e;
0084 e(0, 0) = 0.0015 * 0.0015;
0085 e(1, 1) = 0.0015 * 0.0015;
0086 e(2, 2) = 15. * 15.;
0087 reco::Vertex::Point p(0, 0, 0);
0088 dummy = reco::Vertex(p, e, 0, 0, 0);
0089 }
0090
0091 int idx = 0;
0092 multimap<unsigned int, unsigned int> trackmap;
0093 std::vector<unsigned int> conv_coll(0);
0094
0095
0096 for (unsigned int icoll1 = 0; icoll1 < convColl.size(); icoll1++) {
0097 if ((!convColl[icoll1].quality(reco::Conversion::arbitratedMergedEcalGeneral)) ||
0098 (!convColl[icoll1].quality(reco::Conversion::highPurity)))
0099 continue;
0100
0101 bool greater_prob = false;
0102 std::vector<edm::RefToBase<reco::Track> > tracksRefColl1 = convColl[icoll1].tracks();
0103 for (unsigned it1 = 0; it1 < tracksRefColl1.size(); it1++) {
0104 reco::TrackRef trackRef1 = (tracksRefColl1[it1]).castTo<reco::TrackRef>();
0105
0106 for (unsigned int icoll2 = 0; icoll2 < convColl.size(); icoll2++) {
0107 if (icoll1 == icoll2)
0108 continue;
0109 if ((!convColl[icoll2].quality(reco::Conversion::arbitratedMergedEcalGeneral)) ||
0110 (!convColl[icoll2].quality(reco::Conversion::highPurity)))
0111 continue;
0112 std::vector<edm::RefToBase<reco::Track> > tracksRefColl2 = convColl[icoll2].tracks();
0113 for (unsigned it2 = 0; it2 < tracksRefColl2.size(); it2++) {
0114 reco::TrackRef trackRef2 = (tracksRefColl2[it2]).castTo<reco::TrackRef>();
0115 double like1 = -999;
0116 double like2 = -999;
0117
0118 int shared = 0;
0119 for (auto const& hit1 : trackRef1->recHits())
0120 if (hit1->isValid()) {
0121
0122 for (auto const& hit2 : trackRef2->recHits()) {
0123 if (hit2->isValid() && hit1->sharesInput(hit2, TrackingRecHit::some))
0124 shared++;
0125 }
0126 }
0127 float frac = 0;
0128
0129 float size1 = trackRef1->found();
0130 float size2 = trackRef2->found();
0131
0132 if (size1 > size2)
0133 frac = (double)shared / size2;
0134 else
0135 frac = (double)shared / size1;
0136 if (frac > 0.9) {
0137 like1 = ChiSquaredProbability(convColl[icoll1].conversionVertex().chi2(),
0138 convColl[icoll1].conversionVertex().ndof());
0139 like2 = ChiSquaredProbability(convColl[icoll2].conversionVertex().chi2(),
0140 convColl[icoll2].conversionVertex().ndof());
0141 }
0142 if (like2 > like1) {
0143 greater_prob = true;
0144 break;
0145 }
0146 }
0147
0148 if (greater_prob)
0149 break;
0150 }
0151 if (greater_prob)
0152 break;
0153 }
0154 if (!greater_prob)
0155 conv_coll.push_back(icoll1);
0156 }
0157
0158
0159 for (unsigned iColl = 0; iColl < conv_coll.size(); iColl++) {
0160 unsigned int collindex = conv_coll[iColl];
0161
0162 std::vector<reco::PFRecTrackRef> pfRecTkcoll;
0163
0164 std::vector<edm::RefToBase<reco::Track> > tracksRefColl = convColl[collindex].tracks();
0165
0166 for (unsigned it = 0; it < tracksRefColl.size(); it++) {
0167 reco::TrackRef trackRef = (tracksRefColl[it]).castTo<reco::TrackRef>();
0168 reco::PFRecTrack pfRecTrack(trackRef->charge(), reco::PFRecTrack::KF, trackRef.key(), trackRef);
0169
0170 Trajectory FakeTraj;
0171 bool valid = pfTransformer_->addPoints(pfRecTrack, *trackRef, FakeTraj);
0172 if (valid) {
0173 double stip = -999;
0174 const reco::PFTrajectoryPoint& atECAL = pfRecTrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance);
0175
0176 if (atECAL.isValid()) {
0177 GlobalVector direction(pfRecTrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance).position().x(),
0178 pfRecTrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance).position().y(),
0179 pfRecTrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance).position().z());
0180 stip = IPTools::signedTransverseImpactParameter(thebuilder.build(*trackRef), direction, *pv)
0181 .second.significance();
0182 }
0183 pfRecTrack.setSTIP(stip);
0184 pfRecTkcoll.push_back(reco::PFRecTrackRef(pfTrackRefProd, idx++));
0185 pfRecTrackColl->push_back(pfRecTrack);
0186 }
0187 }
0188
0189 reco::ConversionRef niRef(convCollH, collindex);
0190 pfConversionColl->push_back(reco::PFConversion(niRef, pfRecTkcoll));
0191 }
0192 iEvent.put(std::move(pfRecTrackColl));
0193 iEvent.put(std::move(pfConversionColl));
0194 }
0195
0196
0197 void PFConversionProducer::beginRun(const edm::Run& run, const EventSetup& iSetup) {
0198 auto const& magneticField = &iSetup.getData(magneticFieldToken_);
0199 pfTransformer_ = new PFTrackTransformer(math::XYZVector(magneticField->inTesla(GlobalPoint(0, 0, 0))));
0200 pfTransformer_->OnlyProp();
0201 }
0202
0203
0204 void PFConversionProducer::endRun(const edm::Run& run, const EventSetup& iSetup) {
0205 delete pfTransformer_;
0206 pfTransformer_ = nullptr;
0207 }