RecoDecay/src/BPHKinematicFit.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author Paolo Ronchese INFN Padova
0005  *
0006  */
0007
0008 //-----------------------
0009 // This Class' Header --
0010 //-----------------------
0011 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHKinematicFit.h"
0012
0013 //-------------------------------
0014 // Collaborating Class Headers --
0015 //-------------------------------
0016 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHRecoCandidate.h"
0017 #include "RecoVertex/KinematicFitPrimitives/interface/KinematicParticleFactoryFromTransientTrack.h"
0018 #include "RecoVertex/KinematicFitPrimitives/interface/RefCountedKinematicParticle.h"
0019 #include "RecoVertex/KinematicFit/interface/KinematicParticleVertexFitter.h"
0020 #include "RecoVertex/KinematicFit/interface/KinematicConstrainedVertexFitter.h"
0021 #include "RecoVertex/KinematicFit/interface/KinematicParticleFitter.h"
0022 #include "RecoVertex/KinematicFit/interface/MassKinematicConstraint.h"
0023 #include "RecoVertex/KinematicFit/interface/TwoTrackMassKinematicConstraint.h"
0024 #include "RecoVertex/KinematicFit/interface/MultiTrackMassKinematicConstraint.h"
0025 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0026
0027 //---------------
0028 // C++ Headers --
0029 //---------------
0030 using namespace std;
0031
0032 //-------------------
0033 // Initializations --
0034 //-------------------
0035
0036 //----------------
0037 // Constructors --
0038 //----------------
0039 BPHKinematicFit::BPHKinematicFit(int daugNum, int compNum)
0040     : BPHDecayMomentum(daugNum, compNum),
0041       BPHDecayVertex(nullptr),
0042       massConst(-1.0),
0043       massSigma(-1.0),
0044       oldKPs(true),
0045       oldFit(true),
0046       oldMom(true),
0047       kinTree(nullptr) {}
0048
0049 BPHKinematicFit::BPHKinematicFit(const BPHKinematicFit* ptr)
0050     : BPHDecayVertex(ptr, nullptr),
0051       massConst(-1.0),
0052       massSigma(-1.0),
0053       oldKPs(true),
0054       oldFit(true),
0055       oldMom(true),
0056       kinTree(nullptr) {
0057   map<const reco::Candidate*, const reco::Candidate*> iMap;
0058   const vector<const reco::Candidate*>& daug = daughters();
0059   const vector<Component>& list = ptr->componentList();
0060   int i;
0061   int n = daug.size();
0062   for (i = 0; i < n; ++i) {
0063     const reco::Candidate* cand = daug[i];
0064     iMap[originalReco(cand)] = cand;
0065   }
0066   for (i = 0; i < n; ++i) {
0067     const Component& c = list[i];
0068     dMSig[iMap[c.cand]] = c.msig;
0069   }
0070   const vector<BPHRecoConstCandPtr>& dComp = daughComp();
0071   int j;
0072   int m = dComp.size();
0073   for (j = 0; j < m; ++j) {
0074     const BPHRecoCandidate* rc = dComp[j].get();
0075     const map<const reco::Candidate*, double>& dMap = rc->dMSig;
0076     const map<const BPHRecoCandidate*, FlyingParticle>& cMap = rc->cKinP;
0077     dMSig.insert(dMap.begin(), dMap.end());
0078     cKinP.insert(cMap.begin(), cMap.end());
0079     cKinP[rc];
0080   }
0081 }
0082
0083 //--------------
0084 // Operations --
0085 //--------------
0086 /// apply a mass constraint
0087 void BPHKinematicFit::setConstraint(double mass, double sigma) {
0088   oldKPs = oldFit = oldMom = true;
0089   massConst = mass;
0090   massSigma = sigma;
0091   return;
0092 }
0093
0094 /// retrieve the constraint
0095 double BPHKinematicFit::constrMass() const { return massConst; }
0096
0097 double BPHKinematicFit::constrSigma() const { return massSigma; }
0098
0099 /// set a decaying daughter as an unique particle fitted independently
0100 void BPHKinematicFit::setIndependentFit(const string& name, bool flag, double mass, double sigma) {
0101   string::size_type pos = name.find('/');
0102   if (pos != string::npos) {
0103     edm::LogPrint("WrongRequest") << "BPHKinematicFit::setIndependentFit: "
0104                                   << "cascade decay specification not admitted " << name;
0105     return;
0106   }
0107   const BPHRecoCandidate* comp = getComp(name).get();
0108   if (comp == nullptr) {
0109     edm::LogPrint("ParticleNotFound") << "BPHKinematicFit::setIndependentFit: " << name << " daughter not found";
0110     return;
0111   }
0112   oldKPs = oldFit = oldMom = true;
0113   FlyingParticle& fp = cKinP[comp];
0114   fp.flag = flag;
0115   fp.mass = mass;
0116   fp.sigma = sigma;
0117   return;
0118 }
0119
0120 /// get kinematic particles
0121 const vector<RefCountedKinematicParticle>& BPHKinematicFit::kinParticles() const {
0122   if (oldKPs)
0123     buildParticles();
0124   return allParticles;
0125 }
0126
0127 vector<RefCountedKinematicParticle> BPHKinematicFit::kinParticles(const vector<string>& names) const {
0128   if (oldKPs)
0129     buildParticles();
0130   vector<RefCountedKinematicParticle> plist;
0131   unsigned int m = allParticles.size();
0132   if (m != numParticles())
0133     return plist;
0134   set<RefCountedKinematicParticle> pset;
0135   int i;
0136   int n = names.size();
0137   plist.reserve(m);
0138   for (i = 0; i < n; ++i) {
0139     const string& pname = names[i];
0140     if (pname == "*") {
0141       int j = m;
0142       while (j)
0143         insertParticle(allParticles[--j], plist, pset);
0144       break;
0145     }
0146     string::size_type pos = pname.find('/');
0147     if (pos != string::npos)
0148       getParticles(pname.substr(0, pos), pname.substr(pos + 1), plist, pset, this);
0149     else
0150       getParticles(pname, "", plist, pset, this);
0151   }
0152   return plist;
0153 }
0154
0155 /// perform the kinematic fit and get the result
0156 const RefCountedKinematicTree& BPHKinematicFit::kinematicTree() const {
0157   if (oldFit)
0158     return kinematicTree("", massConst, massSigma);
0159   return kinTree;
0160 }
0161
0162 const RefCountedKinematicTree& BPHKinematicFit::kinematicTree(const string& name, double mass, double sigma) const {
0163   if (mass < 0)
0164     return kinematicTree(name);
0165   if (sigma < 0)
0166     return kinematicTree(name, mass);
0167   ParticleMass mc = mass;
0168   MassKinematicConstraint kinConst(mc, sigma);
0169   return kinematicTree(name, &kinConst);
0170 }
0171
0172 const RefCountedKinematicTree& BPHKinematicFit::kinematicTree(const string& name, double mass) const {
0173   if (mass < 0)
0174     return kinematicTree(name);
0175   vector<RefCountedKinematicParticle> kPart;
0176   const BPHKinematicFit* ptr = splitKP(name, &kPart);
0177   if (ptr == nullptr) {
0178     edm::LogPrint("ParticleNotFound") << "BPHKinematicFit::kinematicTree: " << name << " daughter not found";
0179     return kinTree;
0180   }
0181   int nn = ptr->daughFull().size();
0182   ParticleMass mc = mass;
0183   if (nn == 2) {
0184     TwoTrackMassKinematicConstraint kinConst(mc);
0185     return kinematicTree(kPart, &kinConst);
0186   } else {
0187     MultiTrackMassKinematicConstraint kinConst(mc, nn);
0188     return kinematicTree(kPart, &kinConst);
0189   }
0190 }
0191
0192 const RefCountedKinematicTree& BPHKinematicFit::kinematicTree(const string& name) const {
0193   KinematicConstraint* kc = nullptr;
0194   return kinematicTree(name, kc);
0195 }
0196
0197 const RefCountedKinematicTree& BPHKinematicFit::kinematicTree(const string& name, KinematicConstraint* kc) const {
0198   vector<RefCountedKinematicParticle> kComp;
0199   vector<RefCountedKinematicParticle> kTail;
0200   const BPHKinematicFit* ptr = splitKP(name, &kComp, &kTail);
0201   if (ptr == nullptr) {
0202     edm::LogPrint("ParticleNotFound") << "BPHKinematicFit::kinematicTree: " << name << " daughter not found";
0203     return kinTree;
0204   }
0205   try {
0206     KinematicParticleVertexFitter vtxFitter;
0207     RefCountedKinematicTree compTree = vtxFitter.fit(kComp);
0208     if (compTree->isEmpty())
0209       return kinTree;
0210     if (kc != nullptr) {
0211       KinematicParticleFitter kinFitter;
0212       compTree = kinFitter.fit(kc, compTree);
0213       if (compTree->isEmpty())
0214         return kinTree;
0215     }
0216     compTree->movePointerToTheTop();
0217     if (!kTail.empty()) {
0218       RefCountedKinematicParticle compPart = compTree->currentParticle();
0219       if (!compPart->currentState().isValid())
0220         return kinTree;
0221       kTail.push_back(compPart);
0222       kinTree = vtxFitter.fit(kTail);
0223     } else {
0224       kinTree = compTree;
0225     }
0226   } catch (std::exception const&) {
0227     edm::LogPrint("FitFailed") << "BPHKinematicFit::kinematicTree: "
0228                                << "kin fit reset";
0229     kinTree = RefCountedKinematicTree(nullptr);
0230   }
0231   return kinTree;
0232 }
0233
0234 const RefCountedKinematicTree& BPHKinematicFit::kinematicTree(const string& name,
0235                                                               MultiTrackKinematicConstraint* kc) const {
0236   vector<RefCountedKinematicParticle> kPart;
0237   if (splitKP(name, &kPart) == nullptr) {
0238     edm::LogPrint("ParticleNotFound") << "BPHKinematicFit::kinematicTree: " << name << " daughter not found";
0239     return kinTree;
0240   }
0241   return kinematicTree(kPart, kc);
0242 }
0243
0244 /// reset the kinematic fit
0245 void BPHKinematicFit::resetKinematicFit() const {
0246   oldKPs = oldFit = oldMom = true;
0247   return;
0248 }
0249
0250 /// get fit status
0251 bool BPHKinematicFit::isEmpty() const {
0252   kinematicTree();
0253   if (kinTree.get() == nullptr)
0254     return true;
0255   return kinTree->isEmpty();
0256 }
0257
0258 bool BPHKinematicFit::isValidFit() const {
0259   const RefCountedKinematicParticle kPart = topParticle();
0260   if (kPart.get() == nullptr)
0261     return false;
0262   return kPart->currentState().isValid();
0263 }
0264
0265 /// get current particle
0266 const RefCountedKinematicParticle BPHKinematicFit::currentParticle() const {
0267   if (isEmpty())
0268     return RefCountedKinematicParticle(nullptr);
0269   return kinTree->currentParticle();
0270 }
0271
0272 const RefCountedKinematicVertex BPHKinematicFit::currentDecayVertex() const {
0273   if (isEmpty())
0274     return RefCountedKinematicVertex(nullptr);
0275   return kinTree->currentDecayVertex();
0276 }
0277
0278 /// get top particle
0279 const RefCountedKinematicParticle BPHKinematicFit::topParticle() const {
0280   if (isEmpty())
0281     return RefCountedKinematicParticle(nullptr);
0282   return kinTree->topParticle();
0283 }
0284
0285 const RefCountedKinematicVertex BPHKinematicFit::topDecayVertex() const {
0286   if (isEmpty())
0287     return RefCountedKinematicVertex(nullptr);
0288   kinTree->movePointerToTheTop();
0289   return kinTree->currentDecayVertex();
0290 }
0291
0292 ParticleMass BPHKinematicFit::mass() const {
0293   const RefCountedKinematicParticle kPart = topParticle();
0294   if (kPart.get() == nullptr)
0295     return -1.0;
0296   const KinematicState kStat = kPart->currentState();
0297   if (kStat.isValid())
0298     return kStat.mass();
0299   return -1.0;
0300 }
0301
0302 /// compute total momentum after the fit
0303 const math::XYZTLorentzVector& BPHKinematicFit::p4() const {
0304   if (oldMom)
0305     fitMomentum();
0306   return totalMomentum;
0307 }
0308
0309 /// retrieve particle mass sigma
0310 double BPHKinematicFit::getMassSigma(const reco::Candidate* cand) const {
0311   map<const reco::Candidate*, double>::const_iterator iter = dMSig.find(cand);
0312   return (iter != dMSig.end() ? iter->second : -1);
0313 }
0314
0315 /// retrieve independent fit flag
0316 bool BPHKinematicFit::getIndependentFit(const string& name, double& mass, double& sigma) const {
0317   const BPHRecoCandidate* comp = getComp(name).get();
0318   map<const BPHRecoCandidate*, FlyingParticle>::const_iterator iter = cKinP.find(comp);
0319   if ((iter != cKinP.end()) && iter->second.flag) {
0320     mass = iter->second.mass;
0321     sigma = iter->second.sigma;
0322     return true;
0323   }
0324   return false;
0325 }
0326
0327 /// add a simple particle and specify a criterion to search for
0328 /// the associated track
0329 void BPHKinematicFit::addK(const string& name, const reco::Candidate* daug, double mass, double sigma) {
0330   addK(name, daug, "cfhpmig", mass, sigma);
0331   return;
0332 }
0333
0334 /// add a previously reconstructed particle giving it a name
0335 void BPHKinematicFit::addK(
0336     const string& name, const reco::Candidate* daug, const string& searchList, double mass, double sigma) {
0337   addV(name, daug, searchList, mass);
0338   dMSig[daughters().back()] = sigma;
0339   return;
0340 }
0341
0342 /// add a previously reconstructed particle giving it a name
0343 void BPHKinematicFit::addK(const string& name, const BPHRecoConstCandPtr& comp) {
0344   addV(name, comp);
0345   const map<const reco::Candidate*, double>& dMap = comp->dMSig;
0346   const map<const BPHRecoCandidate*, FlyingParticle>& cMap = comp->cKinP;
0347   dMSig.insert(dMap.begin(), dMap.end());
0348   cKinP.insert(cMap.begin(), cMap.end());
0349   cKinP[comp.get()];
0350   return;
0351 }
0352
0353 // utility function used to cash reconstruction results
0354 void BPHKinematicFit::setNotUpdated() const {
0355   BPHDecayVertex::setNotUpdated();
0356   resetKinematicFit();
0357   return;
0358 }
0359
0360 // build kin particles, perform the fit and compute the total momentum
0361 void BPHKinematicFit::buildParticles() const {
0362   kinMap.clear();
0363   kCDMap.clear();
0364   allParticles.clear();
0365   allParticles.reserve(daughFull().size());
0366   addParticles(allParticles, kinMap, kCDMap);
0367   oldKPs = false;
0368   return;
0369 }
0370
0371 void BPHKinematicFit::addParticles(vector<RefCountedKinematicParticle>& kl,
0372                                    map<const reco::Candidate*, RefCountedKinematicParticle>& km,
0373                                    map<const BPHRecoCandidate*, RefCountedKinematicParticle>& cm) const {
0374   const vector<const reco::Candidate*>& daug = daughters();
0375   KinematicParticleFactoryFromTransientTrack pFactory;
0376   int n = daug.size();
0377   float chi = 0.0;
0378   float ndf = 0.0;
0379   while (n--) {
0380     const reco::Candidate* cand = daug[n];
0381     ParticleMass mass = cand->mass();
0382     float sigma = dMSig.find(cand)->second;
0383     if (sigma < 0)
0384       sigma = 1.0e-7;
0385     reco::TransientTrack* tt = getTransientTrack(cand);
0386     if (tt != nullptr)
0387       kl.push_back(km[cand] = pFactory.particle(*tt, mass, chi, ndf, sigma));
0388   }
0389   const vector<BPHRecoConstCandPtr>& comp = daughComp();
0390   int m = comp.size();
0391   while (m--) {
0392     const BPHRecoCandidate* cptr = comp[m].get();
0393     const FlyingParticle& fp = cKinP.at(cptr);
0394     if (fp.flag) {
0395       BPHRecoCandidate* tptr = cptr->clone();
0396       double mass = fp.mass;
0397       double sigma = fp.sigma;
0398       if (mass > 0.0)
0399         tptr->setConstraint(mass, sigma);
0400       tmpList.push_back(BPHRecoConstCandPtr(tptr));
0401       if (tptr->isEmpty())
0402         return;
0403       if (!tptr->isValidFit())
0404         return;
0405       kl.push_back(cm[cptr] = tptr->topParticle());
0406     } else {
0407       cptr->addParticles(kl, km, cm);
0408     }
0409   }
0410   return;
0411 }
0412
0413 void BPHKinematicFit::getParticles(const string& moth,
0414                                    const string& daug,
0415                                    vector<RefCountedKinematicParticle>& kl,
0416                                    set<RefCountedKinematicParticle>& ks,
0417                                    const BPHKinematicFit* curr) const {
0418   const BPHRecoCandidate* cptr = getComp(moth).get();
0419   if (cptr != nullptr) {
0420     if (curr->cKinP.at(cptr).flag) {
0421       insertParticle(kCDMap[cptr], kl, ks);
0422     } else {
0423       vector<string> list;
0424       if (!daug.empty()) {
0425         list.push_back(daug);
0426       } else {
0427         const vector<string>& dNames = cptr->daugNames();
0428         const vector<string>& cNames = cptr->compNames();
0429         list.insert(list.end(), dNames.begin(), dNames.end());
0430         list.insert(list.end(), cNames.begin(), cNames.end());
0431       }
0432       getParticles(moth, list, kl, ks, cptr);
0433     }
0434     return;
0435   }
0436   const reco::Candidate* dptr = getDaug(moth);
0437   if (dptr != nullptr) {
0438     insertParticle(kinMap[dptr], kl, ks);
0439     return;
0440   }
0441   edm::LogPrint("ParticleNotFound") << "BPHKinematicFit::getParticles: " << moth << " not found";
0442   return;
0443 }
0444
0445 void BPHKinematicFit::getParticles(const string& moth,
0446                                    const vector<string>& daug,
0447                                    vector<RefCountedKinematicParticle>& kl,
0448                                    set<RefCountedKinematicParticle>& ks,
0449                                    const BPHKinematicFit* curr) const {
0450   int i;
0451   int n = daug.size();
0452   for (i = 0; i < n; ++i) {
0453     const string& name = daug[i];
0454     string::size_type pos = name.find('/');
0455     if (pos != string::npos)
0456       getParticles(moth + "/" + name.substr(0, pos), name.substr(pos + 1), kl, ks, curr);
0457     else
0458       getParticles(moth + "/" + name, "", kl, ks, curr);
0459   }
0460   return;
0461 }
0462
0463 unsigned int BPHKinematicFit::numParticles(const BPHKinematicFit* cand) const {
0464   if (cand == nullptr)
0465     cand = this;
0466   unsigned int n = cand->daughters().size();
0467   const vector<string>& cnames = cand->compNames();
0468   int i = cnames.size();
0469   while (i) {
0470     const BPHRecoCandidate* comp = cand->getComp(cnames[--i]).get();
0471     if (cand->cKinP.at(comp).flag)
0472       ++n;
0473     else
0474       n += numParticles(comp);
0475   }
0476   return n;
0477 }
0478
0479 void BPHKinematicFit::insertParticle(RefCountedKinematicParticle& kp,
0480                                      vector<RefCountedKinematicParticle>& kl,
0481                                      set<RefCountedKinematicParticle>& ks) {
0482   if (ks.find(kp) != ks.end())
0483     return;
0484   kl.push_back(kp);
0485   ks.insert(kp);
0486   return;
0487 }
0488
0489 const BPHKinematicFit* BPHKinematicFit::splitKP(const string& name,
0490                                                 vector<RefCountedKinematicParticle>* kComp,
0491                                                 vector<RefCountedKinematicParticle>* kTail) const {
0492   kinTree = RefCountedKinematicTree(nullptr);
0493   oldFit = false;
0494   kinParticles();
0495   if (allParticles.size() != numParticles())
0496     return nullptr;
0497   kComp->clear();
0498   if (kTail == nullptr)
0499     kTail = kComp;
0500   else
0501     kTail->clear();
0502   if (name.empty()) {
0503     *kComp = allParticles;
0504     return this;
0505   }
0506   const BPHRecoCandidate* comp = getComp(name).get();
0507   int ns;
0508   if (comp != nullptr) {
0509     ns = (cKinP.at(comp).flag ? 1 : numParticles(comp));
0510   } else {
0511     edm::LogPrint("ParticleNotFound") << "BPHKinematicFit::splitKP: " << name << " daughter not found";
0512     *kTail = allParticles;
0513     return nullptr;
0514   }
0515   vector<string> nfull(2);
0516   nfull[0] = name;
0517   nfull[1] = "*";
0518   vector<RefCountedKinematicParticle> kPart = kinParticles(nfull);
0519   vector<RefCountedKinematicParticle>::const_iterator iter = kPart.begin();
0520   vector<RefCountedKinematicParticle>::const_iterator imid = iter + ns;
0521   vector<RefCountedKinematicParticle>::const_iterator iend = kPart.end();
0522   kComp->insert(kComp->end(), iter, imid);
0523   kTail->insert(kTail->end(), imid, iend);
0524   return comp;
0525 }
0526
0527 const RefCountedKinematicTree& BPHKinematicFit::kinematicTree(const vector<RefCountedKinematicParticle>& kPart,
0528                                                               MultiTrackKinematicConstraint* kc) const {
0529   try {
0530     KinematicConstrainedVertexFitter cvf;
0531     kinTree = cvf.fit(kPart, kc);
0532   } catch (std::exception const&) {
0533     edm::LogPrint("FitFailed") << "BPHKinematicFit::kinematicTree: "
0534                                << "kin fit reset";
0535     kinTree = RefCountedKinematicTree(nullptr);
0536   }
0537   return kinTree;
0538 }
0539
0540 void BPHKinematicFit::fitMomentum() const {
0541   if (isValidFit()) {
0542     const KinematicState& ks = topParticle()->currentState();
0543     GlobalVector tm = ks.globalMomentum();
0544     double x = tm.x();
0545     double y = tm.y();
0546     double z = tm.z();
0547     double m = ks.mass();
0548     double e = sqrt((x * x) + (y * y) + (z * z) + (m * m));
0549     totalMomentum.SetPxPyPzE(x, y, z, e);
0550   } else {
0551     edm::LogPrint("FitNotFound") << "BPHKinematicFit::fitMomentum: "
0552                                  << "simple momentum sum computed";
0553     math::XYZTLorentzVector tm;
0554     const vector<const reco::Candidate*>& daug = daughters();
0555     int n = daug.size();
0556     while (n--)
0557       tm += daug[n]->p4();
0558     const vector<BPHRecoConstCandPtr>& comp = daughComp();
0559     int m = comp.size();
0560     while (m--)
0561       tm += comp[m]->p4();
0562     totalMomentum = tm;
0563   }
0564   oldMom = false;
0565   return;
0566 }