KalmanVertexFit/plugins/KVFTest.cc

0001 #include "RecoVertex/KalmanVertexFit/plugins/KVFTest.h"
0002
0003 #include "DataFormats/VertexReco/interface/Vertex.h"
0004 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0005 #include "DataFormats/TrackReco/interface/Track.h"
0006 #include "DataFormats/Common/interface/Handle.h"
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0010 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0011 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
0012 #include <iostream>
0013 #include <memory>
0014
0015 using namespace reco;
0016 using namespace edm;
0017 using namespace std;
0018
0019 KVFTest::KVFTest(const edm::ParameterSet& iConfig)
0020     : estoken_MF(esConsumes()),
0021       estoken_TTB(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0022       theConfig(iConfig) {
0023   token_tracks = consumes<TrackCollection>(iConfig.getParameter<string>("TrackLabel"));
0024   outputFile_ = iConfig.getUntrackedParameter<std::string>("outputFile");
0025   kvfPSet = iConfig.getParameter<edm::ParameterSet>("KVFParameters");
0026   rootFile_ = TFile::Open(outputFile_.c_str(), "RECREATE");
0027   edm::LogInfo("RecoVertex/KVFTest") << "Initializing KVF TEST analyser  - Output file: " << outputFile_ << "\n";
0028
0029   token_TrackTruth = consumes<TrackingParticleCollection>(edm::InputTag("trackingtruth", "TrackTruth"));
0030   token_VertexTruth = consumes<TrackingVertexCollection>(edm::InputTag("trackingtruth", "VertexTruth"));
0031   token_associatorForParamAtPca =
0032       consumes<reco::TrackToTrackingParticleAssociator>(edm::InputTag("trackAssociatorByChi2"));
0033 }
0034
0035 KVFTest::~KVFTest() { delete rootFile_; }
0036
0037 void KVFTest::beginJob() {}
0038
0039 void KVFTest::endJob() {}
0040
0041 //
0042 // member functions
0043 //
0044
0045 void KVFTest::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0046   edm::Handle<reco::TrackToTrackingParticleAssociator> associatorForParamAtPca;
0047   iEvent.getByToken(token_associatorForParamAtPca, associatorForParamAtPca);
0048
0049   if (not tree) {
0050     tree = std::make_unique<SimpleVertexTree>("VertexFitter", &iSetup.getData(estoken_MF));
0051   }
0052
0053   edm::LogInfo("RecoVertex/KVFTest") << "Reconstructing event number: " << iEvent.id() << "\n";
0054
0055   // get RECO tracks from the event
0056   // `tks` can be used as a ptr to a reco::TrackCollection
0057   edm::Handle<edm::View<reco::Track> > tks;
0058   iEvent.getByToken(token_tracks, tks);
0059   if (!tks.isValid()) {
0060     edm::LogInfo("RecoVertex/KVFTest") << "Exception during event number: " << iEvent.id() << "\n";
0061   } else {
0062     edm::LogInfo("RecoVertex/KVFTest") << "Found: " << (*tks).size() << " reconstructed tracks"
0063                                        << "\n";
0064     edm::LogPrint("RecoVertex/KVFTest") << "got " << (*tks).size() << " tracks " << std::endl;
0065
0066     // Transform Track to TransientTrack
0067
0068     //get the builder:
0069     const auto& theB = &iSetup.getData(estoken_TTB);
0070     //do the conversion:
0071     std::vector<TransientTrack> t_tks = theB->build(tks);
0072
0073     edm::LogInfo("RecoVertex/KVFTest") << "Found: " << t_tks.size() << " reconstructed tracks"
0074                                        << "\n";
0075
0076     // Call the KalmanVertexFitter if more than 1 track
0077     if (t_tks.size() > 1) {
0078       //      KalmanVertexFitter kvf(kvfPSet);
0079       KalmanVertexFitter kvf(true);
0080       TransientVertex tv = kvf.vertex(t_tks);
0081
0082       edm::LogPrint("RecoVertex/KVFTest") << "Position: " << Vertex::Point(tv.position()) << "\n";
0083
0084       // For the analysis: compare to your SimVertex
0085       TrackingVertex sv = getSimVertex(iEvent);
0086       edm::Handle<TrackingParticleCollection> TPCollectionH;
0087       iEvent.getByToken(token_TrackTruth, TPCollectionH);
0088       if (!TPCollectionH.isValid()) {
0089         edm::LogInfo("RecoVertex/KVFTest") << "Exception during event number: " << iEvent.id() << "\n";
0090       } else {
0091         const TrackingParticleCollection tPC = *(TPCollectionH.product());
0092         reco::RecoToSimCollection recSimColl = associatorForParamAtPca->associateRecoToSim(tks, TPCollectionH);
0093         tree->fill(tv, &sv, &recSimColl);
0094       }
0095     }
0096   }
0097 }
0098
0099 //Returns the first vertex in the list.
0100
0101 TrackingVertex KVFTest::getSimVertex(const edm::Event& iEvent) const {
0102   // get the simulated vertices
0103   edm::Handle<TrackingVertexCollection> TVCollectionH;
0104   iEvent.getByToken(token_VertexTruth, TVCollectionH);
0105   const TrackingVertexCollection tPC = *(TVCollectionH.product());
0106
0107   //    Handle<edm::SimVertexContainer> simVtcs;
0108   //    iEvent.getByLabel("g4SimHits", simVtcs);
0109   //    std::cout << "SimVertex " << simVtcs->size() << std::endl;
0110   //    for(edm::SimVertexContainer::const_iterator v=simVtcs->begin();
0111   //        v!=simVtcs->end(); ++v){
0112   //      std::cout << "simvtx "
0113   //           << v->position().x() << " "
0114   //           << v->position().y() << " "
0115   //           << v->position().z() << " "
0116   //           << v->parentIndex() << " "
0117   //           << v->noParent() << " "
0118   //               << std::endl;
0119   //    }
0120   return *(tPC.begin());
0121 }
0122 DEFINE_FWK_MODULE(KVFTest);