PFTracking/plugins/PFConversionProducer.cc

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   ///Constructor
0020   explicit PFConversionProducer(const edm::ParameterSet&);
0021
0022   ///Destructor
0023   ~PFConversionProducer() override;
0024
0025   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0026
0027 private:
0028   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0029   void endRun(const edm::Run&, const edm::EventSetup&) override;
0030
0031   ///Produce the PFRecTrack collection
0032   void produce(edm::Event&, const edm::EventSetup&) override;
0033
0034   ///PFTrackTransformer
0035   PFTrackTransformer* pfTransformer_;
0036   edm::EDGetTokenT<reco::ConversionCollection> pfConversionContainer_;
0037   edm::EDGetTokenT<reco::VertexCollection> vtx_h;
0038
0039   const edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> transientTrackToken_;
0040   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldToken_;
0041 };
0042
0043 #include "FWCore/Framework/interface/MakerMacros.h"
0044 DEFINE_FWK_MODULE(PFConversionProducer);
0045
0046 void PFConversionProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0047   edm::ParameterSetDescription desc;
0048   desc.add<edm::InputTag>("conversionCollection", {"allConversions", ""});
0049   desc.add<edm::InputTag>("PrimaryVertexLabel", {"offlinePrimaryVertices"});
0050   descriptions.add("pfConversions", desc);
0051 }
0052
0053 typedef std::multimap<unsigned, std::vector<unsigned> > BlockMap;
0054 using namespace std;
0055 using namespace edm;
0056
0057 PFConversionProducer::PFConversionProducer(const ParameterSet& iConfig)
0058     : pfTransformer_(nullptr),
0059       transientTrackToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0060       magneticFieldToken_(esConsumes<edm::Transition::BeginRun>()) {
0061   produces<reco::PFRecTrackCollection>();
0062   produces<reco::PFConversionCollection>();
0063
0064   pfConversionContainer_ = consumes<reco::ConversionCollection>(iConfig.getParameter<InputTag>("conversionCollection"));
0065
0066   vtx_h = consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("PrimaryVertexLabel"));
0067 }
0068
0069 PFConversionProducer::~PFConversionProducer() { delete pfTransformer_; }
0070
0071 void PFConversionProducer::produce(Event& iEvent, const EventSetup& iSetup) {
0072   //create the empty collections
0073   auto pfConversionColl = std::make_unique<reco::PFConversionCollection>();
0074   auto pfRecTrackColl = std::make_unique<reco::PFRecTrackCollection>();
0075
0076   TransientTrackBuilder const& thebuilder = iSetup.getData(transientTrackToken_);
0077
0078   reco::PFRecTrackRefProd pfTrackRefProd = iEvent.getRefBeforePut<reco::PFRecTrackCollection>();
0079   Handle<reco::ConversionCollection> convCollH;
0080   iEvent.getByToken(pfConversionContainer_, convCollH);
0081
0082   const reco::ConversionCollection& convColl = *(convCollH.product());
0083
0084   Handle<reco::VertexCollection> vertex;
0085   iEvent.getByToken(vtx_h, vertex);
0086   //Find PV for IP calculation, if there is no PV in collection than use dummy
0087   reco::Vertex dummy;
0088   const reco::Vertex* pv = &dummy;
0089   if (vertex.isValid()) {
0090     pv = &*vertex->begin();
0091   } else {  // create a dummy PV
0092     reco::Vertex::Error e;
0093     e(0, 0) = 0.0015 * 0.0015;
0094     e(1, 1) = 0.0015 * 0.0015;
0095     e(2, 2) = 15. * 15.;
0096     reco::Vertex::Point p(0, 0, 0);
0097     dummy = reco::Vertex(p, e, 0, 0, 0);
0098   }
0099
0100   int idx = 0;                                    //index of track in PFRecTrack collection
0101   multimap<unsigned int, unsigned int> trackmap;  //Map of Collections and tracks
0102   std::vector<unsigned int> conv_coll(0);
0103
0104   // CLEAN CONVERSION COLLECTION FOR DUPLICATES
0105   for (unsigned int icoll1 = 0; icoll1 < convColl.size(); icoll1++) {
0106     if ((!convColl[icoll1].quality(reco::Conversion::arbitratedMergedEcalGeneral)) ||
0107         (!convColl[icoll1].quality(reco::Conversion::highPurity)))
0108       continue;
0109
0110     bool greater_prob = false;
0111     std::vector<edm::RefToBase<reco::Track> > tracksRefColl1 = convColl[icoll1].tracks();
0112     for (unsigned it1 = 0; it1 < tracksRefColl1.size(); it1++) {
0113       reco::TrackRef trackRef1 = (tracksRefColl1[it1]).castTo<reco::TrackRef>();
0114
0115       for (unsigned int icoll2 = 0; icoll2 < convColl.size(); icoll2++) {
0116         if (icoll1 == icoll2)
0117           continue;
0118         if ((!convColl[icoll2].quality(reco::Conversion::arbitratedMergedEcalGeneral)) ||
0119             (!convColl[icoll2].quality(reco::Conversion::highPurity)))
0120           continue;
0121         std::vector<edm::RefToBase<reco::Track> > tracksRefColl2 = convColl[icoll2].tracks();
0122         for (unsigned it2 = 0; it2 < tracksRefColl2.size(); it2++) {
0123           reco::TrackRef trackRef2 = (tracksRefColl2[it2]).castTo<reco::TrackRef>();
0124           double like1 = -999;
0125           double like2 = -999;
0126           //number of shared hits
0127           int shared = 0;
0128           for (auto const& hit1 : trackRef1->recHits())
0129             if (hit1->isValid()) {
0130               //count number of shared hits
0131               for (auto const& hit2 : trackRef2->recHits()) {
0132                 if (hit2->isValid() && hit1->sharesInput(hit2, TrackingRecHit::some))
0133                   shared++;
0134               }
0135             }
0136           float frac = 0;
0137           //number of valid hits in tracks that are duplicates
0138           float size1 = trackRef1->found();
0139           float size2 = trackRef2->found();
0140           //divide number of shared hits by the total number of hits for the track with less hits
0141           if (size1 > size2)
0142             frac = (double)shared / size2;
0143           else
0144             frac = (double)shared / size1;
0145           if (frac > 0.9) {
0146             like1 = ChiSquaredProbability(convColl[icoll1].conversionVertex().chi2(),
0147                                           convColl[icoll1].conversionVertex().ndof());
0148             like2 = ChiSquaredProbability(convColl[icoll2].conversionVertex().chi2(),
0149                                           convColl[icoll2].conversionVertex().ndof());
0150           }
0151           if (like2 > like1) {
0152             greater_prob = true;
0153             break;
0154           }
0155         }  //end loop over tracks in collection 2
0156
0157         if (greater_prob)
0158           break;  //if a duplicate track is found in a collection with greater Chi^2 probability for Vertex fit then break out of comparison loop
0159       }  //end loop over collection 2 checking
0160       if (greater_prob)
0161         break;  //if a duplicate track is found in a collection with greater Chi^2 probability for Vertex fit then one does not need to check the other track the collection will not be stored
0162     }  //end loop over tracks in collection 1
0163     if (!greater_prob)
0164       conv_coll.push_back(icoll1);
0165   }  //end loop over collection 1
0166
0167   //Finally fill empty collections
0168   for (unsigned iColl = 0; iColl < conv_coll.size(); iColl++) {
0169     unsigned int collindex = conv_coll[iColl];
0170     //std::cout<<"Filling this collection"<<collindex<<endl;
0171     std::vector<reco::PFRecTrackRef> pfRecTkcoll;
0172
0173     std::vector<edm::RefToBase<reco::Track> > tracksRefColl = convColl[collindex].tracks();
0174     // convert the secondary tracks
0175     for (unsigned it = 0; it < tracksRefColl.size(); it++) {
0176       reco::TrackRef trackRef = (tracksRefColl[it]).castTo<reco::TrackRef>();
0177       reco::PFRecTrack pfRecTrack(trackRef->charge(), reco::PFRecTrack::KF, trackRef.key(), trackRef);
0178       //std::cout<<"Track Pt "<<trackRef->pt()<<std::endl;
0179       Trajectory FakeTraj;
0180       bool valid = pfTransformer_->addPoints(pfRecTrack, *trackRef, FakeTraj);
0181       if (valid) {
0182         double stip = -999;
0183         const reco::PFTrajectoryPoint& atECAL = pfRecTrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance);
0184         //if extrapolation to ECAL is valid then calculate STIP
0185         if (atECAL.isValid()) {
0186           GlobalVector direction(pfRecTrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance).position().x(),
0187                                  pfRecTrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance).position().y(),
0188                                  pfRecTrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance).position().z());
0189           stip = IPTools::signedTransverseImpactParameter(thebuilder.build(*trackRef), direction, *pv)
0190                      .second.significance();
0191         }
0192         pfRecTrack.setSTIP(stip);
0193         pfRecTkcoll.push_back(reco::PFRecTrackRef(pfTrackRefProd, idx++));
0194         pfRecTrackColl->push_back(pfRecTrack);
0195       }
0196     }  //end loop over tracks
0197     //store reference to the Conversion collection
0198     reco::ConversionRef niRef(convCollH, collindex);
0199     pfConversionColl->push_back(reco::PFConversion(niRef, pfRecTkcoll));
0200   }  //end loop over collections
0201   iEvent.put(std::move(pfRecTrackColl));
0202   iEvent.put(std::move(pfConversionColl));
0203 }
0204
0205 // ------------ method called once each job just before starting event loop  ------------
0206 void PFConversionProducer::beginRun(const edm::Run& run, const EventSetup& iSetup) {
0207   auto const& magneticField = &iSetup.getData(magneticFieldToken_);
0208   pfTransformer_ = new PFTrackTransformer(math::XYZVector(magneticField->inTesla(GlobalPoint(0, 0, 0))));
0209   pfTransformer_->OnlyProp();
0210 }
0211
0212 // ------------ method called once each job just after ending the event loop  ------------
0213 void PFConversionProducer::endRun(const edm::Run& run, const EventSetup& iSetup) {
0214   delete pfTransformer_;
0215   pfTransformer_ = nullptr;
0216 }