GaussianSumVertexFit/test/GsfTest.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    GsfTest
0004 // Class:      GsfTest
0005 //
0006 /**\class GsfTest GsfTest.cc RecoVertex/GsfTest/src/GsfTest.cc
0007
0008  Description: steers tracker primary vertex reconstruction and storage
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013
0014 // system include files
0015 #include <memory>
0016 #include <iostream>
0017
0018 // user include files
0019 #include "DataFormats/Common/interface/Handle.h"
0020 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0021 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
0022 #include "DataFormats/VertexReco/interface/Vertex.h"
0023 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0024 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0025 #include "FWCore/Framework/interface/ESHandle.h"
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/EventSetup.h"
0028 #include "FWCore/Framework/interface/Frameworkfwd.h"
0029 #include "FWCore/PluginManager/interface/ModuleDef.h"
0030 #include "FWCore/Framework/interface/MakerMacros.h"
0031 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0032 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0033 #include "MagneticField/Engine/interface/MagneticField.h"
0034 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0035 #include "RecoVertex/GaussianSumVertexFit/interface/GsfVertexFitter.h"
0036 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
0037 #include "RecoVertex/KalmanVertexFit/interface/SimpleVertexTree.h"
0038 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0039 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
0040 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0041 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0042 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0043 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0044
0045 #include <TFile.h>
0046
0047 /**
0048    * This is a very simple test analyzer mean to test the KalmanVertexFitter
0049    */
0050
0051 class GsfTest : public edm::one::EDAnalyzer<> {
0052 public:
0053   explicit GsfTest(const edm::ParameterSet&);
0054   ~GsfTest();
0055
0056   virtual void analyze(const edm::Event&, const edm::EventSetup&);
0057
0058 private:
0059   TrackingVertex getSimVertex(const edm::Event& iEvent) const;
0060
0061   const edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> estoken_ttk;
0062   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> estoken_mf;
0063
0064   edm::ParameterSet gsfPSet;
0065
0066   std::unique_ptr<SimpleVertexTree> tree;
0067   TFile* rootFile_;
0068
0069   std::string outputFile_;  // output file
0070   edm::EDGetTokenT<reco::TrackCollection> token_tracks;
0071   //   edm::EDGetTokenT<TrackingParticleCollection> token_TrackTruth;
0072   edm::EDGetTokenT<TrackingVertexCollection> token_VertexTruth;
0073 };
0074
0075 using namespace reco;
0076 using namespace edm;
0077 using namespace std;
0078
0079 GsfTest::GsfTest(const edm::ParameterSet& iConfig)
0080     : estoken_ttk(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))), estoken_mf(esConsumes()) {
0081   token_tracks = consumes<TrackCollection>(iConfig.getParameter<string>("TrackLabel"));
0082   outputFile_ = iConfig.getUntrackedParameter<std::string>("outputFile");
0083   gsfPSet = iConfig.getParameter<edm::ParameterSet>("GsfParameters");
0084   rootFile_ = TFile::Open(outputFile_.c_str(), "RECREATE");
0085   edm::LogInfo("RecoVertex/GsfTest") << "Initializing KVF TEST analyser  - Output file: " << outputFile_ << "\n";
0086   //   token_TrackTruth = consumes<TrackingParticleCollection>(edm::InputTag("trackingtruth", "TrackTruth"));
0087   token_VertexTruth = consumes<TrackingVertexCollection>(edm::InputTag("trackingtruth", "VertexTruth"));
0088 }
0089
0090 GsfTest::~GsfTest() { delete rootFile_; }
0091
0092 //
0093 // member functions
0094 //
0095
0096 void GsfTest::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0097   if (not tree) {
0098     const MagneticField* magField = &iSetup.getData(estoken_mf);
0099     tree.reset(new SimpleVertexTree("VertexFitter", magField));
0100   }
0101
0102   try {
0103     edm::LogInfo("RecoVertex/GsfTest") << "Reconstructing event number: " << iEvent.id() << "\n";
0104
0105     // get RECO tracks from the event
0106     // `tks` can be used as a ptr to a reco::TrackCollection
0107     edm::Handle<edm::View<reco::Track> > tks;
0108     iEvent.getByToken(token_tracks, tks);
0109
0110     edm::LogInfo("RecoVertex/GsfTest") << "Found: " << (*tks).size() << " reconstructed tracks"
0111                                        << "\n";
0112
0113     // Transform Track to TransientTrack
0114
0115     //get the builder:
0116     edm::ESHandle<TransientTrackBuilder> theB = iSetup.getHandle(estoken_ttk);
0117     //do the conversion:
0118     vector<TransientTrack> t_tks = (*theB).build(tks);
0119
0120     edm::LogInfo("RecoVertex/GsfTest") << "Found: " << t_tks.size() << " reconstructed tracks"
0121                                        << "\n";
0122
0123     // Call the KalmanVertexFitter if more than 1 track
0124     if (t_tks.size() > 1) {
0125       //      KalmanVertexFitter kvf(kvfPSet);
0126       // For the analysis: compare to your SimVertex
0127       //       TrackingVertex sv = getSimVertex(iEvent);
0128       //       edm::LogPrint("GsfTest") << "SimV Position: " << Vertex::Point(sv.position()) << "\n";
0129
0130       KalmanVertexFitter kvf;
0131       TransientVertex tv2 = kvf.vertex(t_tks);
0132       edm::LogPrint("GsfTest") << "KVF Position:  " << Vertex::Point(tv2.position()) << tv2.normalisedChiSquared()
0133                                << " " << tv2.degreesOfFreedom() << "\n";
0134
0135       GsfVertexFitter gsf(gsfPSet);
0136       TransientVertex tv = gsf.vertex(t_tks);
0137       edm::LogPrint("GsfTest") << "Position: " << Vertex::Point(tv.position()) << "\n";
0138
0139       //   edm::Handle<TrackingParticleCollection>  TPCollectionH ;
0140       //   iEvent.getByToken(token_TrackTruth, TPCollectionH);
0141       //   const TrackingParticleCollection tPC = *(TPCollectionH.product());
0142       //       reco::RecoToSimCollection recSimColl=associatorForParamAtPca->associateRecoToSim(tks,
0143       //                                          TPCollectionH,
0144       //                                          &iEvent);
0145
0146       //       tree->fill(tv, &sv, 0, 0.0);
0147     }
0148
0149   }
0150
0151   catch (cms::Exception& err) {
0152     edm::LogError("GsfTest") << "Exception during event number: " << iEvent.id() << "\n" << err.what() << "\n";
0153   }
0154 }
0155
0156 //Returns the first vertex in the list.
0157
0158 TrackingVertex GsfTest::getSimVertex(const edm::Event& iEvent) const {
0159   // get the simulated vertices
0160   edm::Handle<TrackingVertexCollection> TVCollectionH;
0161   iEvent.getByToken(token_VertexTruth, TVCollectionH);
0162   const TrackingVertexCollection tPC = *(TVCollectionH.product());
0163
0164   //    Handle<edm::SimVertexContainer> simVtcs;
0165   //    iEvent.getByLabel("g4SimHits", simVtcs);
0166   //    edm::LogPrint("GsfTest") << "SimVertex " << simVtcs->size() << std::endl;
0167   //    for(edm::SimVertexContainer::const_iterator v=simVtcs->begin();
0168   //        v!=simVtcs->end(); ++v){
0169   //      edm::LogPrint("GsfTest") << "simvtx "
0170   //           << v->position().x() << " "
0171   //           << v->position().y() << " "
0172   //           << v->position().z() << " "
0173   //           << v->parentIndex() << " "
0174   //           << v->noParent() << " "
0175   //               << std::endl;
0176   //    }
0177   return *(tPC.begin());
0178 }
0179
0180 DEFINE_FWK_MODULE(GsfTest);