KinematicFit/plugins/KineExample.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    KineExample
0004 // Class:      KineExample
0005 //
0006 /**\class KineExample KineExample.cc RecoVertex/KineExample/src/KineExample.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/TrackReco/interface/Track.h"
0021 #include "DataFormats/TrackReco/interface/TrackFwd.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/Framework/interface/MakerMacros.h"
0030 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
0033 #include "RecoVertex/KinematicFit/interface/KinematicConstrainedVertexFitter.h"
0034 #include "RecoVertex/KinematicFit/interface/KinematicParticleFitter.h"
0035 #include "RecoVertex/KinematicFit/interface/KinematicParticleVertexFitter.h"
0036 #include "RecoVertex/KinematicFit/interface/MassKinematicConstraint.h"
0037 #include "RecoVertex/KinematicFit/interface/TwoTrackMassKinematicConstraint.h"
0038 #include "RecoVertex/KinematicFitPrimitives/interface/MultiTrackKinematicConstraint.h"
0039 #include "RecoVertex/KinematicFitPrimitives/interface/ParticleMass.h"
0040 #include "RecoVertex/KinematicFitPrimitives/interface/RefCountedKinematicParticle.h"
0041 #include "RecoVertex/KinematicFitPrimitives/interface/RefCountedKinematicTree.h"
0042 #include "RecoVertex/KinematicFitPrimitives/interface/RefCountedKinematicVertex.h"
0043 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0044 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
0045 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0046 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0047 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0048 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0049 #include <RecoVertex/KinematicFitPrimitives/interface/KinematicParticleFactoryFromTransientTrack.h>
0050 //#include "MagneticField/Engine/interface/MagneticField.h"
0051 //#include "MagneticField/Records/interface/IdealMagneticField.h"
0052
0053 #include <TFile.h>
0054
0055 /**
0056    * This is a very simple test analyzer mean to test the KalmanVertexFitter
0057    */
0058
0059 class KineExample : public edm::one::EDAnalyzer<edm::one::WatchRuns> {
0060 public:
0061   explicit KineExample(const edm::ParameterSet&);
0062   ~KineExample() override;
0063
0064   void analyze(const edm::Event&, const edm::EventSetup&) override;
0065   void beginRun(edm::Run const&, edm::EventSetup const&) override;
0066   void endRun(edm::Run const&, edm::EventSetup const&) override{};
0067
0068 private:
0069   void printout(const RefCountedKinematicVertex& myVertex) const;
0070   void printout(const RefCountedKinematicParticle& myParticle) const;
0071   void printout(const RefCountedKinematicTree& myTree) const;
0072
0073   TrackingVertex getSimVertex(const edm::Event& iEvent) const;
0074
0075   const edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> estoken_ttk;
0076   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> estoken_mf;
0077
0078   edm::ParameterSet theConfig;
0079   edm::ParameterSet kvfPSet;
0080   //   TFile*  rootFile_;
0081
0082   std::string outputFile_;  // output file
0083   edm::EDGetTokenT<reco::TrackCollection> token_tracks;
0084   //   edm::EDGetTokenT<TrackingParticleCollection> token_TrackTruth;
0085   edm::EDGetTokenT<TrackingVertexCollection> token_VertexTruth;
0086 };
0087
0088 using namespace reco;
0089 using namespace edm;
0090 using namespace std;
0091
0092 KineExample::KineExample(const edm::ParameterSet& iConfig)
0093     : estoken_ttk(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0094       estoken_mf(esConsumes<edm::Transition::BeginRun>()) {
0095   token_tracks = consumes<TrackCollection>(iConfig.getParameter<string>("TrackLabel"));
0096   outputFile_ = iConfig.getUntrackedParameter<std::string>("outputFile");
0097   kvfPSet = iConfig.getParameter<edm::ParameterSet>("KVFParameters");
0098   //   rootFile_ = TFile::Open(outputFile_.c_str(),"RECREATE");
0099   edm::LogInfo("RecoVertex/KineExample") << "Initializing KVF TEST analyser  - Output file: " << outputFile_ << "\n";
0100   //   token_TrackTruth = consumes<TrackingParticleCollection>(edm::InputTag("trackingtruth", "TrackTruth"));
0101   token_VertexTruth = consumes<TrackingVertexCollection>(edm::InputTag("trackingtruth", "VertexTruth"));
0102 }
0103
0104 KineExample::~KineExample() = default;
0105
0106 void KineExample::beginRun(Run const& run, EventSetup const& setup) {
0107   //const MagneticField* magField = &setup.getData(estoken_mf);
0108   //tree.reset(new SimpleVertexTree("VertexFitter", magField.product()));
0109 }
0110
0111 //
0112 // member functions
0113 //
0114
0115 void KineExample::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0116   try {
0117     edm::LogPrint("KineExample") << "Reconstructing event number: " << iEvent.id() << "\n";
0118
0119     // get RECO tracks from the event
0120     // `tks` can be used as a ptr to a reco::TrackCollection
0121     edm::Handle<reco::TrackCollection> tks;
0122     iEvent.getByToken(token_tracks, tks);
0123     if (!tks.isValid()) {
0124       edm::LogPrint("KineExample") << "Couln't find track collection: " << iEvent.id() << "\n";
0125     } else {
0126       edm::LogInfo("RecoVertex/KineExample") << "Found: " << (*tks).size() << " reconstructed tracks"
0127                                              << "\n";
0128       edm::LogPrint("KineExample") << "got " << (*tks).size() << " tracks " << endl;
0129
0130       // Transform Track to TransientTrack
0131       //get the builder:
0132       edm::ESHandle<TransientTrackBuilder> theB = iSetup.getHandle(estoken_ttk);
0133       //do the conversion:
0134       vector<TransientTrack> t_tks = (*theB).build(tks);
0135
0136       edm::LogPrint("KineExample") << "Found: " << t_tks.size() << " reconstructed tracks"
0137                                    << "\n";
0138
0139       // Do a KindFit, if >= 4 tracks.
0140       if (t_tks.size() > 3) {
0141         // For a first test, suppose that the first four tracks are the 2 muons,
0142         // then the 2 kaons. Since this will not be true, the result of the fit
0143         // will not be meaningfull, but at least you will get the idea of how to
0144         // do such a fit.
0145
0146         //First, to get started, a simple vertex fit:
0147
0148         vector<TransientTrack> ttv;
0149         ttv.push_back(t_tks[0]);
0150         ttv.push_back(t_tks[1]);
0151         ttv.push_back(t_tks[2]);
0152         ttv.push_back(t_tks[3]);
0153         KalmanVertexFitter kvf(false);
0154         TransientVertex tv = kvf.vertex(ttv);
0155         if (!tv.isValid())
0156           edm::LogPrint("KineExample") << "KVF failed\n";
0157         else
0158           edm::LogPrint("KineExample") << "KVF fit Position: " << Vertex::Point(tv.position()) << "\n";
0159
0160         TransientTrack ttMuPlus = t_tks[0];
0161         TransientTrack ttMuMinus = t_tks[1];
0162         TransientTrack ttKPlus = t_tks[2];
0163         TransientTrack ttKMinus = t_tks[3];
0164
0165         //the final state muons and kaons from the Bs->J/PsiPhi->mumuKK decay
0166         //Creating a KinematicParticleFactory
0167         KinematicParticleFactoryFromTransientTrack pFactory;
0168
0169         //The mass of a muon and the insignificant mass sigma to avoid singularities in the covariance matrix.
0170         ParticleMass muon_mass = 0.1056583;
0171         ParticleMass kaon_mass = 0.493677;
0172         float muon_sigma = 0.0000001;
0173         float kaon_sigma = 0.000016;
0174
0175         //initial chi2 and ndf before kinematic fits. The chi2 of the reconstruction is not considered
0176         float chi = 0.;
0177         float ndf = 0.;
0178
0179         //making particles
0180         vector<RefCountedKinematicParticle> muonParticles;
0181         vector<RefCountedKinematicParticle> phiParticles;
0182         vector<RefCountedKinematicParticle> allParticles;
0183         muonParticles.push_back(pFactory.particle(ttMuPlus, muon_mass, chi, ndf, muon_sigma));
0184         muonParticles.push_back(pFactory.particle(ttMuMinus, muon_mass, chi, ndf, muon_sigma));
0185         allParticles.push_back(pFactory.particle(ttMuPlus, muon_mass, chi, ndf, muon_sigma));
0186         allParticles.push_back(pFactory.particle(ttMuMinus, muon_mass, chi, ndf, muon_sigma));
0187
0188         phiParticles.push_back(pFactory.particle(ttKPlus, kaon_mass, chi, ndf, kaon_sigma));
0189         phiParticles.push_back(pFactory.particle(ttKMinus, kaon_mass, chi, ndf, kaon_sigma));
0190         allParticles.push_back(pFactory.particle(ttKPlus, kaon_mass, chi, ndf, kaon_sigma));
0191         allParticles.push_back(pFactory.particle(ttKMinus, kaon_mass, chi, ndf, kaon_sigma));
0192
0193         /* Example of a simple vertex fit, without other constraints
0194        * The reconstructed decay tree is a result of the kinematic fit
0195        * The KinematicParticleVertexFitter fits the final state particles to their vertex and
0196        * reconstructs the decayed state
0197        */
0198         KinematicParticleVertexFitter fitter;
0199         edm::LogPrint("KineExample") << "Simple vertex fit with KinematicParticleVertexFitter:\n";
0200         RefCountedKinematicTree vertexFitTree = fitter.fit(allParticles);
0201
0202         printout(vertexFitTree);
0203
0204         /////Example of global fit:
0205
0206         //creating the constraint for the J/Psi mass
0207         ParticleMass jpsi = 3.09687;
0208
0209         //creating the two track mass constraint
0210         MultiTrackKinematicConstraint* j_psi_c = new TwoTrackMassKinematicConstraint(jpsi);
0211
0212         //creating the fitter
0213         KinematicConstrainedVertexFitter kcvFitter;
0214
0215         //obtaining the resulting tree
0216         RefCountedKinematicTree myTree = kcvFitter.fit(allParticles, j_psi_c);
0217
0218         edm::LogPrint("KineExample") << "\nGlobal fit done:\n";
0219         printout(myTree);
0220
0221         //creating the vertex fitter
0222         KinematicParticleVertexFitter kpvFitter;
0223
0224         //reconstructing a J/Psi decay
0225         RefCountedKinematicTree jpTree = kpvFitter.fit(muonParticles);
0226
0227         //creating the particle fitter
0228         KinematicParticleFitter csFitter;
0229
0230         // creating the constraint
0231         float jp_m_sigma = 0.00004;
0232         KinematicConstraint* jpsi_c2 = new MassKinematicConstraint(jpsi, jp_m_sigma);
0233
0234         //the constrained fit:
0235         jpTree = csFitter.fit(jpsi_c2, jpTree);
0236
0237         //getting the J/Psi KinematicParticle and putting it together with the kaons.
0238         //The J/Psi KinematicParticle has a pointer to the tree it belongs to
0239         jpTree->movePointerToTheTop();
0240         RefCountedKinematicParticle jpsi_part = jpTree->currentParticle();
0241         phiParticles.push_back(jpsi_part);
0242
0243         //making a vertex fit and thus reconstructing the Bs parameters
0244         // the resulting tree includes all the final state tracks, the J/Psi meson,
0245         // its decay vertex, the Bs meson and its decay vertex.
0246         RefCountedKinematicTree bsTree = kpvFitter.fit(phiParticles);
0247         edm::LogPrint("KineExample") << "Sequential fit done:\n";
0248         printout(bsTree);
0249
0250         //       // For the analysis: compare to your SimVertex
0251         //       TrackingVertex sv = getSimVertex(iEvent);
0252         //   edm::Handle<TrackingParticleCollection>  TPCollectionH ;
0253         //   iEvent.getByToken(token_TrackTruth, TPCollectionH);
0254         //   const TrackingParticleCollection tPC = *(TPCollectionH.product());
0255         //       reco::RecoToSimCollection recSimColl=associatorForParamAtPca->associateRecoToSim(tks,
0256         //                                        TPCollectionH,
0257         //                                        &iEvent,
0258         //                                                                            &iSetup);
0259         //
0260         //       tree->fill(tv, &sv, &recSimColl);
0261         //     }
0262       }
0263     }
0264   } catch (cms::Exception& err) {
0265     edm::LogError("KineExample") << "Exception during event number: " << iEvent.id() << "\n" << err.what() << "\n";
0266   }
0267 }
0268
0269 void KineExample::printout(const RefCountedKinematicVertex& myVertex) const {
0270   if (myVertex->vertexIsValid()) {
0271     edm::LogPrint("KineExample") << "Decay vertex: " << myVertex->position() << myVertex->chiSquared() << " "
0272                                  << myVertex->degreesOfFreedom() << endl;
0273   } else
0274     edm::LogPrint("KineExample") << "Decay vertex Not valid\n";
0275 }
0276
0277 void KineExample::printout(const RefCountedKinematicParticle& myParticle) const {
0278   edm::LogPrint("KineExample") << "Particle: \n";
0279   //accessing the reconstructed Bs meson parameters:
0280   //SK: uncomment if needed  AlgebraicVector7 bs_par = myParticle->currentState().kinematicParameters().vector();
0281
0282   //and their joint covariance matrix:
0283   //SK:uncomment if needed  AlgebraicSymMatrix77 bs_er = myParticle->currentState().kinematicParametersError().matrix();
0284   edm::LogPrint("KineExample") << "Momentum at vertex: " << myParticle->currentState().globalMomentum() << endl;
0285   edm::LogPrint("KineExample") << "Parameters at vertex: " << myParticle->currentState().kinematicParameters().vector()
0286                                << endl;
0287 }
0288
0289 void KineExample::printout(const RefCountedKinematicTree& myTree) const {
0290   if (!myTree->isValid()) {
0291     edm::LogPrint("KineExample") << "Tree is invalid. Fit failed.\n";
0292     return;
0293   }
0294
0295   //accessing the tree components, move pointer to top
0296   myTree->movePointerToTheTop();
0297
0298   //We are now at the top of the decay tree getting the B_s reconstructed KinematicPartlcle
0299   RefCountedKinematicParticle b_s = myTree->currentParticle();
0300   printout(b_s);
0301
0302   // The B_s decay vertex
0303   RefCountedKinematicVertex b_dec_vertex = myTree->currentDecayVertex();
0304   printout(b_dec_vertex);
0305
0306   // Get all the children of Bs:
0307   //In this way, the pointer is not moved
0308   vector<RefCountedKinematicParticle> bs_children = myTree->finalStateParticles();
0309
0310   for (unsigned int i = 0; i < bs_children.size(); ++i) {
0311     printout(bs_children[i]);
0312   }
0313
0314   //Now navigating down the tree , pointer is moved:
0315   bool child = myTree->movePointerToTheFirstChild();
0316
0317   if (child)
0318     while (myTree->movePointerToTheNextChild()) {
0319       RefCountedKinematicParticle aChild = myTree->currentParticle();
0320       printout(aChild);
0321     }
0322 }
0323
0324 //Returns the first vertex in the list.
0325
0326 TrackingVertex KineExample::getSimVertex(const edm::Event& iEvent) const {
0327   // get the simulated vertices
0328   edm::Handle<TrackingVertexCollection> TVCollectionH;
0329   iEvent.getByToken(token_VertexTruth, TVCollectionH);
0330   const TrackingVertexCollection tPC = *(TVCollectionH.product());
0331
0332   //    Handle<edm::SimVertexContainer> simVtcs;
0333   //    iEvent.getByLabel("g4SimHits", simVtcs);
0334   //    edm::LogPrint("KineExample") << "SimVertex " << simVtcs->size() << std::endl;
0335   //    for(edm::SimVertexContainer::const_iterator v=simVtcs->begin();
0336   //        v!=simVtcs->end(); ++v){
0337   //      edm::LogPrint("KineExample") << "simvtx "
0338   //           << v->position().x() << " "
0339   //           << v->position().y() << " "
0340   //           << v->position().z() << " "
0341   //           << v->parentIndex() << " "
0342   //           << v->noParent() << " "
0343   //               << std::endl;
0344   //    }
0345   return *(tPC.begin());
0346 }
0347 DEFINE_FWK_MODULE(KineExample);