RecoUtils/src/CandKinematicVertexFitter.cc

0001 #include "PhysicsTools/RecoUtils/interface/CandKinematicVertexFitter.h"
0002 #include "RecoVertex/KinematicFit/interface/KinematicParticleFitter.h"
0003 #include "RecoVertex/KinematicFit/interface/MassKinematicConstraint.h"
0004 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0005 #include "DataFormats/Candidate/interface/VertexCompositeCandidate.h"
0006 #include "FWCore/Utilities/interface/EDMException.h"
0007 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0008 #include "FWCore/Framework/interface/ESHandle.h"
0009 #include <sstream>
0010 #include <iostream>
0011 using namespace reco;
0012 using namespace std;
0013
0014 // perform the kinematic fit

0015 bool CandKinematicVertexFitter::fit(const vector<RefCountedKinematicParticle> &particles) const {
0016   try {
0017     tree_ = fitter_.fit(particles);
0018   } catch (std::exception &err) {
0019     std::cerr << ">>> exception thrown by KinematicParticleVertexFitter:\n"
0020               << err.what() << "\n"
0021               << ">>> candidate not fitted to common vertex" << std::endl;
0022     return false;
0023   }
0024   //check tree_ is valid here!

0025   if (tree_->isValid())
0026     return true;
0027   else
0028     return false;
0029 }
0030
0031 // main method called by CandProducer sets the VertexCompositeCandidate

0032 void CandKinematicVertexFitter::set(VertexCompositeCandidate &c) const {
0033   if (bField_ == nullptr)
0034     throw edm::Exception(edm::errors::InvalidReference)
0035         << "B-Field was not set up CandKinematicVertexFitter.\n"
0036         << "the following method must be called before fitting a candidate:\n"
0037         << " CandKinematicVertexFitter:.set( const MagneticField * )" << endl;
0038   vector<RefCountedKinematicParticle> particles;
0039   vector<Candidate *> daughters;
0040   vector<RecoCandidate::TrackType> trackTypes;
0041   // fill particles with KinematicParticles and daughters with Candidates of the daughters of c

0042   fill(particles, daughters, trackTypes, c);
0043   assert(particles.size() == daughters.size());
0044
0045   // attempt to fit the KinematicParticles, particles

0046   if (fit(particles)) {
0047     // after the fit, tree_ contains the KinematicTree from the fit

0048     tree_->movePointerToTheTop();
0049     // set the kinematic properties of the daughters from the fit

0050     RefCountedKinematicVertex vertex = tree_->currentDecayVertex();
0051     if (vertex->vertexIsValid()) {
0052       Candidate::Point vtx(vertex->position());
0053       c.setVertex(vtx);
0054       vector<RefCountedKinematicParticle> treeParticles = tree_->daughterParticles();
0055       vector<RefCountedKinematicParticle>::const_iterator particleIt = treeParticles.begin();
0056       vector<Candidate *>::const_iterator daughterIt = daughters.begin(), daughtersEnd = daughters.end();
0057       vector<RecoCandidate::TrackType>::const_iterator trackTypeIt = trackTypes.begin();
0058       Candidate::LorentzVector mp4(0, 0, 0, 0);
0059       for (; daughterIt != daughtersEnd; ++particleIt, ++daughterIt, ++trackTypeIt) {
0060         Candidate &daughter = **daughterIt;
0061         GlobalVector p3 = (*particleIt)->currentState().globalMomentum();
0062         double px = p3.x(), py = p3.y(), pz = p3.z(), p = p3.mag();
0063         double energy;
0064
0065         if (!daughter.longLived())
0066           daughter.setVertex(vtx);
0067         double scale;
0068         switch (*trackTypeIt) {
0069           case RecoCandidate::gsfTrackType:
0070             //gsf used for electron tracks

0071             energy = daughter.energy();
0072             scale = energy / p;
0073             px *= scale;
0074             py *= scale;
0075             pz *= scale;
0076             [[fallthrough]];
0077           default:
0078             double mass = (*particleIt)->currentState().mass();
0079             energy = sqrt(p * p + mass * mass);
0080         };
0081         Candidate::LorentzVector dp4(px, py, pz, energy);
0082         daughter.setP4(dp4);
0083         mp4 += dp4;
0084       }
0085       c.setP4(mp4);
0086       c.setChi2AndNdof(chi2_ = vertex->chiSquared(), ndof_ = vertex->degreesOfFreedom());
0087       GlobalError err = vertex->error();
0088       cov_(0, 0) = err.cxx();
0089       cov_(0, 1) = err.cyx();
0090       cov_(0, 2) = err.czx();
0091       cov_(1, 2) = err.czy();
0092       cov_(1, 1) = err.cyy();
0093       cov_(2, 2) = err.czz();
0094       c.setCovariance(cov_);
0095     }
0096   } else {
0097     c.setChi2AndNdof(chi2_ = -1, ndof_ = 0);
0098     c.setCovariance(cov_ = CovarianceMatrix(ROOT::Math::SMatrixIdentity()));
0099   }
0100 }
0101
0102 // methond to fill the properties of a CompositeCandidate's daughters

0103 void CandKinematicVertexFitter::fill(vector<RefCountedKinematicParticle> &particles,
0104                                      vector<Candidate *> &daughters,
0105                                      vector<RecoCandidate::TrackType> &trackTypes,
0106                                      Candidate &c) const {
0107   size_t nDau = c.numberOfDaughters();
0108   // loop through CompositeCandidate daughters

0109   for (unsigned int j = 0; j < nDau; ++j) {
0110     Candidate *d = c.daughter(j);
0111     if (d == nullptr) {
0112       ostringstream message;
0113       message << "Can't access in write mode candidate daughters. "
0114               << "pdgId = " << c.pdgId() << ".\n";
0115       const Candidate *d1 = c.daughter(j);
0116       if (d1 == nullptr)
0117         message << "Null daughter also found in read-only mode\n";
0118       else
0119         message << "Daughter found in read-only mode with id: " << d1->pdgId() << "\n";
0120       throw edm::Exception(edm::errors::InvalidReference) << message.str();
0121     }
0122     //check for a daughter which itself is a composite

0123     if (d->numberOfDaughters() > 0) {
0124       //try to cast to VertexCompositeCandiate

0125       VertexCompositeCandidate *vtxDau = dynamic_cast<VertexCompositeCandidate *>(d);
0126       if (vtxDau != nullptr && vtxDau->vertexChi2() > 0) {
0127         // if VertexCompositeCandidate refit vtxDau via the set method

0128         (*this).set(*vtxDau);
0129         // if mass constraint is desired, do it here

0130         if (vtxDau->massConstraint()) {
0131           KinematicParticleFitter csFitter;
0132           //get particle mass from pdg table via pdgid number

0133           const ParticleData *data = pdt_->particle(vtxDau->pdgId());
0134           ParticleMass mass = data->mass();
0135           float mass_sigma = mass * 0.000001;  //needs a sigma for the fit

0136           // create a KinematicConstraint and refit the tree with it

0137           //KinematicConstraint * mass_c = new MassKinematicConstraint(mass,mass_sigma);

0138           MassKinematicConstraint mkc(mass, mass_sigma);
0139           KinematicConstraint *mass_c(&mkc);
0140           tree_ = csFitter.fit(mass_c, tree_);
0141           //CHECK THIS! the following works, but might not be safe

0142           //tree_ = csFitter.fit(&(MassKinematicConstraint(mass,mass_sigma)),tree_);

0143         }
0144         // add the kinematic particle from the fit to particles

0145         RefCountedKinematicParticle current = (*this).currentParticle();
0146         particles.push_back(current);
0147         daughters.push_back(d);
0148         trackTypes.push_back(RecoCandidate::noTrackType);
0149       } else {
0150         fill(particles, daughters, trackTypes, *d);
0151       }
0152     } else {
0153       //get track, make KinematicParticle and add to particles so it can be fit

0154       TrackRef trk = d->get<TrackRef>();
0155       RecoCandidate::TrackType type = d->get<RecoCandidate::TrackType>();
0156       if (!trk.isNull()) {
0157         TransientTrack trTrk(trk, bField_);
0158         float chi2 = 0, ndof = 0;
0159         ParticleMass mass = d->mass();
0160         float sigma = mass * 1.e-6;
0161         particles.push_back(factory_.particle(trTrk, mass, chi2, ndof, sigma));
0162         daughters.push_back(d);
0163         trackTypes.push_back(type);
0164       } else {
0165         cerr << ">>> warning: candidate of type " << d->pdgId() << " has no track reference." << endl;
0166       }
0167     }
0168   }
0169 }