RecoUtils/src/CandMassKinFitter.cc

0001 #include "PhysicsTools/RecoUtils/interface/CandMassKinFitter.h"
0002 #include "PhysicsTools/KinFitter/interface/TKinFitter.h"
0003 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0004 #include "PhysicsTools/KinFitter/interface/TFitConstraintM.h"
0005 #include "PhysicsTools/KinFitter/interface/TFitParticleMCCart.h"
0006 #include "TMatrixD.h"
0007 #include <iostream>
0008 using namespace reco;
0009 using namespace std;
0010
0011 FitQuality CandMassKinFitter::set(Candidate& c) const {
0012   TKinFitter fitter("CandMassFit", "CandMassFit");
0013   TString name("dau0");
0014   size_t daus = c.numberOfDaughters();
0015   vector<TMatrixD> errors(daus, TMatrix(3, 3));
0016   vector<TVector3> momenta(daus);
0017   vector<TFitParticleMCCart*> particles(daus, nullptr);
0018   TFitConstraintM constraint("MassConstraint", "MassConstraint", nullptr, nullptr, mass_);
0019   for (size_t i = 0; i < daus; ++i) {
0020     const Candidate& dau = *c.daughter(i);
0021     const Particle::LorentzVector& p4 = dau.p4();
0022     TMatrixD& err = errors[i];
0023     TVector3& mom = momenta[i];
0024     mom = TVector3(p4.px(), p4.py(), p4.pz());
0025     TrackRef trk = dau.get<TrackRef>();
0026     // dummy errors for now...

0027     // should play with track parametrization...

0028     err.Zero();
0029     err(0, 0) = 0.1;
0030     err(1, 1) = 0.1;
0031     err(2, 2) = 0.1;
0032     fitter.addMeasParticle(particles[i] = new TFitParticleMCCart(name, name, &mom, dau.mass(), &err));
0033     name[3]++;
0034     constraint.addParticle1(particles[i]);
0035   }
0036   fitter.addConstraint(&constraint);
0037   fitter.setMaxNbIter(30);
0038   fitter.setMaxDeltaS(1e-2);
0039   fitter.setMaxF(1e-1);
0040   fitter.setVerbosity(0);
0041   fitter.fit();
0042   // if (  fitter->getStatus() != 0 ) throw ...

0043   TLorentzVector sum(0, 0, 0, 0);
0044   for (size_t i = 0; i < daus; ++i) {
0045     Candidate& dau = *c.daughter(i);
0046     TFitParticleMCCart* part = particles[i];
0047     const TLorentzVector* p4 = part->getCurr4Vec();
0048     dau.setP4(Particle::LorentzVector(p4->X(), p4->Y(), p4->Z(), p4->T()));
0049     sum += *p4;
0050     delete particles[i];
0051   }
0052   c.setP4(Particle::LorentzVector(sum.X(), sum.Y(), sum.Z(), sum.T()));
0053   return FitQuality(fitter.getS(), fitter.getNDF());
0054 }