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 private:
0026   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0027   void endRun(const edm::Run&, const edm::EventSetup&) override;
0028
0029   ///Produce the PFRecTrack collection
0030   void produce(edm::Event&, const edm::EventSetup&) override;
0031
0032   ///PFTrackTransformer
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   //create the empty collections
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   //Find PV for IP calculation, if there is no PV in collection than use dummy
0078   reco::Vertex dummy;
0079   const reco::Vertex* pv = &dummy;
0080   if (vertex.isValid()) {
0081     pv = &*vertex->begin();
0082   } else {  // create a dummy PV
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;                                    //index of track in PFRecTrack collection
0092   multimap<unsigned int, unsigned int> trackmap;  //Map of Collections and tracks
0093   std::vector<unsigned int> conv_coll(0);
0094
0095   // CLEAN CONVERSION COLLECTION FOR DUPLICATES
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           //number of shared hits
0118           int shared = 0;
0119           for (auto const& hit1 : trackRef1->recHits())
0120             if (hit1->isValid()) {
0121               //count number of shared hits
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           //number of valid hits in tracks that are duplicates
0129           float size1 = trackRef1->found();
0130           float size2 = trackRef2->found();
0131           //divide number of shared hits by the total number of hits for the track with less hits
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         }  //end loop over tracks in collection 2
0147
0148         if (greater_prob)
0149           break;  //if a duplicate track is found in a collection with greater Chi^2 probability for Vertex fit then break out of comparison loop
0150       }           //end loop over collection 2 checking
0151       if (greater_prob)
0152         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
0153     }           //end loop over tracks in collection 1
0154     if (!greater_prob)
0155       conv_coll.push_back(icoll1);
0156   }  //end loop over collection 1
0157
0158   //Finally fill empty collections
0159   for (unsigned iColl = 0; iColl < conv_coll.size(); iColl++) {
0160     unsigned int collindex = conv_coll[iColl];
0161     //std::cout<<"Filling this collection"<<collindex<<endl;
0162     std::vector<reco::PFRecTrackRef> pfRecTkcoll;
0163
0164     std::vector<edm::RefToBase<reco::Track> > tracksRefColl = convColl[collindex].tracks();
0165     // convert the secondary tracks
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       //std::cout<<"Track Pt "<<trackRef->pt()<<std::endl;
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         //if extrapolation to ECAL is valid then calculate STIP
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     }  //end loop over tracks
0188     //store reference to the Conversion collection
0189     reco::ConversionRef niRef(convCollH, collindex);
0190     pfConversionColl->push_back(reco::PFConversion(niRef, pfRecTkcoll));
0191   }  //end loop over collections
0192   iEvent.put(std::move(pfRecTrackColl));
0193   iEvent.put(std::move(pfConversionColl));
0194 }
0195
0196 // ------------ method called once each job just before starting event loop  ------------
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 // ------------ method called once each job just after ending the event loop  ------------
0204 void PFConversionProducer::endRun(const edm::Run& run, const EventSetup& iSetup) {
0205   delete pfTransformer_;
0206   pfTransformer_ = nullptr;
0207 }