EgammaPhotonAlgos/src/ConversionVertexFinder.cc

0001 #include <iostream>
0002 #include <memory>
0003 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionVertexFinder.h"
0004 // Framework
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 //
0008 //Kinematic constraint vertex fitter
0009 #include <RecoVertex/KinematicFitPrimitives/interface/KinematicParticleFactoryFromTransientTrack.h>
0010
0011 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0012
0013 // new templated one
0014 #include "RecoVertex/KinematicFit/interface/KinematicConstrainedVertexFitterT.h"
0015 #include "RecoVertex/KinematicFit/interface/ColinearityKinematicConstraintT.h"
0016
0017 //
0018
0019 ConversionVertexFinder::ConversionVertexFinder(const edm::ParameterSet& config) : conf_(config) {
0020   LogDebug("ConversionVertexFinder") << "ConversionVertexFinder CTOR  "
0021                                      << "\n";
0022   maxDelta_ = conf_.getParameter<double>("maxDelta");
0023   maxReducedChiSq_ = conf_.getParameter<double>("maxReducedChiSq");
0024   minChiSqImprovement_ = conf_.getParameter<double>("minChiSqImprovement");
0025   maxNbrOfIterations_ = conf_.getParameter<int>("maxNbrOfIterations");
0026   kcvFitter_ = new KinematicConstrainedVertexFitter();
0027   kcvFitter_->setParameters(conf_);
0028 }
0029
0030 ConversionVertexFinder::~ConversionVertexFinder() {
0031   LogDebug("ConversionVertexFinder") << "ConversionVertexFinder DTOR "
0032                                      << "\n";
0033   delete kcvFitter_;
0034 }
0035
0036 bool ConversionVertexFinder::run(const std::vector<reco::TransientTrack>& _pair, reco::Vertex& the_vertex) {
0037   std::vector<reco::TransientTrack> pair = _pair;
0038   bool found = false;
0039
0040   if (pair.size() < 2)
0041     return found;
0042
0043   float sigma = 0.00000001;
0044   float chi = 0.;
0045   float ndf = 0.;
0046   float mass = 0.000511;
0047
0048   KinematicParticleFactoryFromTransientTrack pFactory;
0049
0050   std::vector<RefCountedKinematicParticle> particles;
0051
0052   particles.push_back(
0053       pFactory.particle(pair[0], mass, chi, ndf, sigma, *pair[0].innermostMeasurementState().freeState()));
0054   particles.push_back(
0055       pFactory.particle(pair[1], mass, chi, ndf, sigma, *pair[1].innermostMeasurementState().freeState()));
0056
0057 #ifdef OldKineFit
0058   ColinearityKinematicConstraint constr(ColinearityKinematicConstraint::PhiTheta);
0059
0060   RefCountedKinematicTree myTree = kcvFitter_->fit(particles, &constr);
0061 #else
0062
0063   // bizzare way to the get the field...
0064   const MagneticField* mf = pair[0].field();
0065
0066   ColinearityKinematicConstraintT<colinearityKinematic::PhiTheta> constr;
0067   KinematicConstrainedVertexFitterT<2, 2> kcvFitter(mf);
0068   kcvFitter.setParameters(conf_);
0069   RefCountedKinematicTree myTree = kcvFitter.fit(particles, &constr);
0070
0071 #ifdef KineFitDebug
0072
0073   ColinearityKinematicConstraint oldconstr(ColinearityKinematicConstraint::PhiTheta);
0074
0075   RefCountedKinematicTree oldTree = kcvFitter_->fit(particles, &oldconstr);
0076
0077   if (oldTree->isValid()) {
0078     std::cout << "old " << kcvFitter_->getNit() << std::endl;
0079     RefCountedKinematicVertex gamma_dec_vertex = oldTree->currentDecayVertex();
0080     if (gamma_dec_vertex->vertexIsValid())
0081       std::cout << gamma_dec_vertex->chiSquared() << " " << gamma_dec_vertex->degreesOfFreedom() << std::endl;
0082     std::cout << oldTree->currentParticle()->currentState().globalMomentum()
0083               << oldTree->currentParticle()->currentState().globalPosition() << std::endl;
0084     std::vector<RefCountedKinematicParticle> fStates = oldTree->finalStateParticles();
0085     for (unsigned int kk = 0; kk < fStates.size(); kk++) {
0086       std::cout << fStates[kk]->currentState().globalMomentum() << fStates[kk]->currentState().globalPosition()
0087                 << std::endl;
0088       std::cout << fStates[kk]->currentState().kinematicParametersError().matrix() << std::endl;
0089     }
0090   } else
0091     std::cout << "old invalid " << kcvFitter_->getNit() << std::endl;
0092
0093   if (myTree->isValid()) {
0094     std::cout << "new " << kcvFitter.getNit() << std::endl;
0095     RefCountedKinematicVertex gamma_dec_vertex = myTree->currentDecayVertex();
0096     if (gamma_dec_vertex->vertexIsValid())
0097       std::cout << gamma_dec_vertex->chiSquared() << " " << gamma_dec_vertex->degreesOfFreedom() << std::endl;
0098     std::cout << myTree->currentParticle()->currentState().globalMomentum()
0099               << myTree->currentParticle()->currentState().globalPosition() << std::endl;
0100     std::vector<RefCountedKinematicParticle> fStates = myTree->finalStateParticles();
0101     for (unsigned int kk = 0; kk < fStates.size(); kk++) {
0102       std::cout << fStates[kk]->currentState().globalMomentum() << fStates[kk]->currentState().globalPosition()
0103                 << std::endl;
0104       std::cout << fStates[kk]->currentState().kinematicParametersError().matrix() << std::endl;
0105     }
0106   } else
0107     std::cout << "new invalid " << kcvFitter.getNit() << std::endl;
0108
0109 #endif  // TemplateKineFitDebug
0110
0111 #endif
0112
0113   if (myTree->isValid()) {
0114     myTree->movePointerToTheTop();
0115     RefCountedKinematicParticle the_photon = myTree->currentParticle();
0116     if (the_photon->currentState().isValid()) {
0117       //const ParticleMass photon_mass = the_photon->currentState().mass();
0118       RefCountedKinematicVertex gamma_dec_vertex;
0119       gamma_dec_vertex = myTree->currentDecayVertex();
0120       if (gamma_dec_vertex->vertexIsValid()) {
0121         const float chi2Prob =
0122             ChiSquaredProbability(gamma_dec_vertex->chiSquared(), gamma_dec_vertex->degreesOfFreedom());
0123         if (chi2Prob > 0.) {  // no longer cut here, only ask positive probability here
0124           //const math::XYZPoint vtxPos(gamma_dec_vertex->position());
0125           the_vertex = *gamma_dec_vertex;
0126           found = true;
0127         }
0128       }
0129     }
0130   }
0131
0132   return found;
0133 }
0134
0135 TransientVertex ConversionVertexFinder::run(const std::vector<reco::TransientTrack>& pair) {
0136   LogDebug("ConversionVertexFinder") << "ConversionVertexFinder run pair size " << pair.size() << "\n";
0137
0138   //for ( std::vector<reco::TransientTrack>::const_iterator iTk=pair.begin(); iTk!=pair.end(); ++iTk) {
0139   // LogDebug("ConversionVertexFinder") << "  ConversionVertexFinder  Tracks in the pair  charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner momentum " << iTk->track().innerMomentum() << "\n";
0140   //}
0141
0142   reco::Vertex theVertex;
0143   KalmanVertexFitter fitter(true);
0144   TransientVertex transientVtx;
0145
0146   const std::string metname = "ConversionVertexFinder| ConversionVertexFinder";
0147   try {
0148     transientVtx = fitter.vertex(pair);
0149
0150   } catch (cms::Exception& e) {
0151     edm::LogWarning(metname) << "cms::Exception caught in ConversionVertexFinder::run\n" << e.explainSelf();
0152   }
0153
0154   return transientVtx;
0155 }