Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:49:20

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()
0040     : BPHDecayVertex(nullptr),
0041       massConst(-1.0),
0042       massSigma(-1.0),
0043       oldKPs(true),
0044       oldFit(true),
0045       oldMom(true),
0046       kinTree(nullptr) {}
0047 
0048 BPHKinematicFit::BPHKinematicFit(const BPHKinematicFit* ptr)
0049     : BPHDecayVertex(ptr, nullptr),
0050       massConst(-1.0),
0051       massSigma(-1.0),
0052       oldKPs(true),
0053       oldFit(true),
0054       oldMom(true),
0055       kinTree(nullptr) {
0056   map<const reco::Candidate*, const reco::Candidate*> iMap;
0057   const vector<const reco::Candidate*>& daug = daughters();
0058   const vector<Component>& list = ptr->componentList();
0059   int i;
0060   int n = daug.size();
0061   for (i = 0; i < n; ++i) {
0062     const reco::Candidate* cand = daug[i];
0063     iMap[originalReco(cand)] = cand;
0064   }
0065   for (i = 0; i < n; ++i) {
0066     const Component& c = list[i];
0067     dMSig[iMap[c.cand]] = c.msig;
0068   }
0069   const vector<BPHRecoConstCandPtr>& dComp = daughComp();
0070   int j;
0071   int m = dComp.size();
0072   for (j = 0; j < m; ++j) {
0073     const BPHRecoCandidate* rc = dComp[j].get();
0074     const map<const reco::Candidate*, double>& dMap = rc->dMSig;
0075     const map<const BPHRecoCandidate*, FlyingParticle>& cMap = rc->cKinP;
0076     dMSig.insert(dMap.begin(), dMap.end());
0077     cKinP.insert(cMap.begin(), cMap.end());
0078     cKinP[rc];
0079   }
0080 }
0081 
0082 //--------------
0083 // Destructor --
0084 //--------------
0085 BPHKinematicFit::~BPHKinematicFit() {}
0086 
0087 //--------------
0088 // Operations --
0089 //--------------
0090 /// apply a mass constraint
0091 void BPHKinematicFit::setConstraint(double mass, double sigma) {
0092   oldKPs = oldFit = oldMom = true;
0093   massConst = mass;
0094   massSigma = sigma;
0095   return;
0096 }
0097 
0098 /// retrieve the constraint
0099 double BPHKinematicFit::constrMass() const { return massConst; }
0100 
0101 double BPHKinematicFit::constrSigma() const { return massSigma; }
0102 
0103 /// set a decaying daughter as an unique particle fitted independently
0104 void BPHKinematicFit::setIndependentFit(const string& name, bool flag, double mass, double sigma) {
0105   string::size_type pos = name.find('/');
0106   if (pos != string::npos) {
0107     edm::LogPrint("WrongRequest") << "BPHKinematicFit::setIndependentFit: "
0108                                   << "cascade decay specification not admitted " << name;
0109     return;
0110   }
0111   const BPHRecoCandidate* comp = getComp(name).get();
0112   if (comp == nullptr) {
0113     edm::LogPrint("ParticleNotFound") << "BPHKinematicFit::setIndependentFit: " << name << " daughter not found";
0114     return;
0115   }
0116   oldKPs = oldFit = oldMom = true;
0117   FlyingParticle& fp = cKinP[comp];
0118   fp.flag = flag;
0119   fp.mass = mass;
0120   fp.sigma = sigma;
0121   return;
0122 }
0123 
0124 /// get kinematic particles
0125 const vector<RefCountedKinematicParticle>& BPHKinematicFit::kinParticles() const {
0126   if (oldKPs)
0127     buildParticles();
0128   return allParticles;
0129 }
0130 
0131 vector<RefCountedKinematicParticle> BPHKinematicFit::kinParticles(const vector<string>& names) const {
0132   if (oldKPs)
0133     buildParticles();
0134   vector<RefCountedKinematicParticle> plist;
0135   unsigned int m = allParticles.size();
0136   if (m != numParticles())
0137     return plist;
0138   set<RefCountedKinematicParticle> pset;
0139   int i;
0140   int n = names.size();
0141   plist.reserve(m);
0142   for (i = 0; i < n; ++i) {
0143     const string& pname = names[i];
0144     if (pname == "*") {
0145       int j = m;
0146       while (j)
0147         insertParticle(allParticles[--j], plist, pset);
0148       break;
0149     }
0150     string::size_type pos = pname.find('/');
0151     if (pos != string::npos)
0152       getParticles(pname.substr(0, pos), pname.substr(pos + 1), plist, pset);
0153     else
0154       getParticles(pname, "", plist, pset);
0155   }
0156   return plist;
0157 }
0158 
0159 /// perform the kinematic fit and get the result
0160 const RefCountedKinematicTree& BPHKinematicFit::kinematicTree() const {
0161   if (oldFit)
0162     return kinematicTree("", massConst, massSigma);
0163   return kinTree;
0164 }
0165 
0166 const RefCountedKinematicTree& BPHKinematicFit::kinematicTree(const string& name, double mass, double sigma) const {
0167   if (mass < 0)
0168     return kinematicTree(name);
0169   if (sigma < 0)
0170     return kinematicTree(name, mass);
0171   ParticleMass mc = mass;
0172   MassKinematicConstraint kinConst(mc, sigma);
0173   return kinematicTree(name, &kinConst);
0174 }
0175 
0176 const RefCountedKinematicTree& BPHKinematicFit::kinematicTree(const string& name, double mass) const {
0177   if (mass < 0)
0178     return kinematicTree(name);
0179   vector<RefCountedKinematicParticle> kPart;
0180   const BPHKinematicFit* ptr = splitKP(name, &kPart);
0181   if (ptr == nullptr) {
0182     edm::LogPrint("ParticleNotFound") << "BPHKinematicFit::kinematicTree: " << name << " daughter not found";
0183     return kinTree;
0184   }
0185 
0186   int nn = ptr->daughFull().size();
0187   ParticleMass mc = mass;
0188   if (nn == 2) {
0189     TwoTrackMassKinematicConstraint kinConst(mc);
0190     return kinematicTree(kPart, &kinConst);
0191   } else {
0192     MultiTrackMassKinematicConstraint kinConst(mc, nn);
0193     return kinematicTree(kPart, &kinConst);
0194   }
0195 }
0196 
0197 const RefCountedKinematicTree& BPHKinematicFit::kinematicTree(const string& name) const {
0198   KinematicConstraint* kc = nullptr;
0199   return kinematicTree(name, kc);
0200 }
0201 
0202 const RefCountedKinematicTree& BPHKinematicFit::kinematicTree(const string& name, KinematicConstraint* kc) const {
0203   vector<RefCountedKinematicParticle> kComp;
0204   vector<RefCountedKinematicParticle> kTail;
0205   const BPHKinematicFit* ptr = splitKP(name, &kComp, &kTail);
0206   if (ptr == nullptr) {
0207     edm::LogPrint("ParticleNotFound") << "BPHKinematicFit::kinematicTree: " << name << " daughter not found";
0208     return kinTree;
0209   }
0210   try {
0211     KinematicParticleVertexFitter vtxFitter;
0212     RefCountedKinematicTree compTree = vtxFitter.fit(kComp);
0213     if (compTree->isEmpty())
0214       return kinTree;
0215     if (kc != nullptr) {
0216       KinematicParticleFitter kinFitter;
0217       compTree = kinFitter.fit(kc, compTree);
0218       if (compTree->isEmpty())
0219         return kinTree;
0220     }
0221     compTree->movePointerToTheTop();
0222     if (!kTail.empty()) {
0223       RefCountedKinematicParticle compPart = compTree->currentParticle();
0224       if (!compPart->currentState().isValid())
0225         return kinTree;
0226       kTail.push_back(compPart);
0227       kinTree = vtxFitter.fit(kTail);
0228     } else {
0229       kinTree = compTree;
0230     }
0231   } catch (std::exception const&) {
0232     edm::LogPrint("FitFailed") << "BPHKinematicFit::kinematicTree: "
0233                                << "kin fit reset";
0234     kinTree = RefCountedKinematicTree(nullptr);
0235   }
0236   return kinTree;
0237 }
0238 
0239 const RefCountedKinematicTree& BPHKinematicFit::kinematicTree(const string& name,
0240                                                               MultiTrackKinematicConstraint* kc) const {
0241   vector<RefCountedKinematicParticle> kPart;
0242   if (splitKP(name, &kPart) == nullptr) {
0243     edm::LogPrint("ParticleNotFound") << "BPHKinematicFit::kinematicTree: " << name << " daughter not found";
0244     return kinTree;
0245   }
0246   return kinematicTree(kPart, kc);
0247 }
0248 
0249 /// reset the kinematic fit
0250 void BPHKinematicFit::resetKinematicFit() const {
0251   oldKPs = oldFit = oldMom = true;
0252   return;
0253 }
0254 
0255 /// get fit status
0256 bool BPHKinematicFit::isEmpty() const {
0257   kinematicTree();
0258   if (kinTree.get() == nullptr)
0259     return true;
0260   return kinTree->isEmpty();
0261 }
0262 
0263 bool BPHKinematicFit::isValidFit() const {
0264   const RefCountedKinematicParticle kPart = topParticle();
0265   if (kPart.get() == nullptr)
0266     return false;
0267   return kPart->currentState().isValid();
0268 }
0269 
0270 /// get current particle
0271 const RefCountedKinematicParticle BPHKinematicFit::currentParticle() const {
0272   if (isEmpty())
0273     return RefCountedKinematicParticle(nullptr);
0274   return kinTree->currentParticle();
0275 }
0276 
0277 const RefCountedKinematicVertex BPHKinematicFit::currentDecayVertex() const {
0278   if (isEmpty())
0279     return RefCountedKinematicVertex(nullptr);
0280   return kinTree->currentDecayVertex();
0281 }
0282 
0283 /// get top particle
0284 const RefCountedKinematicParticle BPHKinematicFit::topParticle() const {
0285   if (isEmpty())
0286     return RefCountedKinematicParticle(nullptr);
0287   return kinTree->topParticle();
0288 }
0289 
0290 const RefCountedKinematicVertex BPHKinematicFit::topDecayVertex() const {
0291   if (isEmpty())
0292     return RefCountedKinematicVertex(nullptr);
0293   kinTree->movePointerToTheTop();
0294   return kinTree->currentDecayVertex();
0295 }
0296 
0297 ParticleMass BPHKinematicFit::mass() const {
0298   const RefCountedKinematicParticle kPart = topParticle();
0299   if (kPart.get() == nullptr)
0300     return -1.0;
0301   const KinematicState kStat = kPart->currentState();
0302   if (kStat.isValid())
0303     return kStat.mass();
0304   return -1.0;
0305 }
0306 
0307 /// compute total momentum after the fit
0308 const math::XYZTLorentzVector& BPHKinematicFit::p4() const {
0309   if (oldMom)
0310     fitMomentum();
0311   return totalMomentum;
0312 }
0313 
0314 /// retrieve particle mass sigma
0315 double BPHKinematicFit::getMassSigma(const reco::Candidate* cand) const {
0316   map<const reco::Candidate*, double>::const_iterator iter = dMSig.find(cand);
0317   return (iter != dMSig.end() ? iter->second : -1);
0318 }
0319 
0320 /// retrieve independent fit flag
0321 bool BPHKinematicFit::getIndependentFit(const std::string& name) const {
0322   const BPHRecoCandidate* comp = getComp(name).get();
0323   map<const BPHRecoCandidate*, FlyingParticle>::const_iterator iter = cKinP.find(comp);
0324   return (iter != cKinP.end() ? iter->second.flag : 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) && (sigma > 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) const {
0417   const BPHRecoCandidate* cptr = getComp(moth).get();
0418   if (cptr != nullptr) {
0419     if (cKinP.at(cptr).flag) {
0420       insertParticle(kCDMap[cptr], kl, ks);
0421     } else {
0422       vector<string> list;
0423       if (!daug.empty()) {
0424         list.push_back(daug);
0425       } else {
0426         const vector<string>& dNames = cptr->daugNames();
0427         const vector<string>& cNames = cptr->compNames();
0428         list.insert(list.end(), dNames.begin(), dNames.end());
0429         list.insert(list.end(), cNames.begin(), cNames.end());
0430       }
0431       getParticles(moth, list, kl, ks);
0432     }
0433     return;
0434   }
0435   const reco::Candidate* dptr = getDaug(moth);
0436   if (dptr != nullptr) {
0437     insertParticle(kinMap[dptr], kl, ks);
0438     return;
0439   }
0440   edm::LogPrint("ParticleNotFound") << "BPHKinematicFit::getParticles: " << moth << " not found";
0441   return;
0442 }
0443 
0444 void BPHKinematicFit::getParticles(const string& moth,
0445                                    const vector<string>& daug,
0446                                    vector<RefCountedKinematicParticle>& kl,
0447                                    set<RefCountedKinematicParticle>& ks) const {
0448   int i;
0449   int n = daug.size();
0450   for (i = 0; i < n; ++i) {
0451     const string& name = daug[i];
0452     string::size_type pos = name.find('/');
0453     if (pos != string::npos)
0454       getParticles(moth + "/" + name.substr(0, pos), name.substr(pos + 1), kl, ks);
0455     else
0456       getParticles(moth + "/" + name, "", kl, ks);
0457   }
0458   return;
0459 }
0460 
0461 unsigned int BPHKinematicFit::numParticles(const BPHKinematicFit* cand) const {
0462   if (cand == nullptr)
0463     cand = this;
0464   unsigned int n = cand->daughters().size();
0465   const vector<string>& cnames = cand->compNames();
0466   int i = cnames.size();
0467   while (i) {
0468     const BPHRecoCandidate* comp = cand->getComp(cnames[--i]).get();
0469     if (cKinP.at(comp).flag)
0470       ++n;
0471     else
0472       n += numParticles(comp);
0473   }
0474   return n;
0475 }
0476 
0477 void BPHKinematicFit::insertParticle(RefCountedKinematicParticle& kp,
0478                                      vector<RefCountedKinematicParticle>& kl,
0479                                      set<RefCountedKinematicParticle>& ks) {
0480   if (ks.find(kp) != ks.end())
0481     return;
0482   kl.push_back(kp);
0483   ks.insert(kp);
0484   return;
0485 }
0486 
0487 const BPHKinematicFit* BPHKinematicFit::splitKP(const string& name,
0488                                                 vector<RefCountedKinematicParticle>* kComp,
0489                                                 vector<RefCountedKinematicParticle>* kTail) const {
0490   kinTree = RefCountedKinematicTree(nullptr);
0491   oldFit = false;
0492   kinParticles();
0493   if (allParticles.size() != numParticles())
0494     return nullptr;
0495   kComp->clear();
0496   if (kTail == nullptr)
0497     kTail = kComp;
0498   else
0499     kTail->clear();
0500   if (name.empty()) {
0501     *kComp = allParticles;
0502     return this;
0503   }
0504   const BPHRecoCandidate* comp = getComp(name).get();
0505   int ns;
0506   if (comp != nullptr) {
0507     ns = (cKinP.at(comp).flag ? 1 : numParticles(comp));
0508   } else {
0509     edm::LogPrint("ParticleNotFound") << "BPHKinematicFit::splitKP: " << name << " daughter not found";
0510     *kTail = allParticles;
0511     return nullptr;
0512   }
0513   vector<string> nfull(2);
0514   nfull[0] = name;
0515   nfull[1] = "*";
0516   vector<RefCountedKinematicParticle> kPart = kinParticles(nfull);
0517   vector<RefCountedKinematicParticle>::const_iterator iter = kPart.begin();
0518   vector<RefCountedKinematicParticle>::const_iterator imid = iter + ns;
0519   vector<RefCountedKinematicParticle>::const_iterator iend = kPart.end();
0520   kComp->insert(kComp->end(), iter, imid);
0521   kTail->insert(kTail->end(), imid, iend);
0522   return comp;
0523 }
0524 
0525 const RefCountedKinematicTree& BPHKinematicFit::kinematicTree(const vector<RefCountedKinematicParticle>& kPart,
0526                                                               MultiTrackKinematicConstraint* kc) const {
0527   try {
0528     KinematicConstrainedVertexFitter cvf;
0529     kinTree = cvf.fit(kPart, kc);
0530   } catch (std::exception const&) {
0531     edm::LogPrint("FitFailed") << "BPHKinematicFit::kinematicTree: "
0532                                << "kin fit reset";
0533     kinTree = RefCountedKinematicTree(nullptr);
0534   }
0535   return kinTree;
0536 }
0537 
0538 void BPHKinematicFit::fitMomentum() const {
0539   if (isValidFit()) {
0540     const KinematicState& ks = topParticle()->currentState();
0541     GlobalVector tm = ks.globalMomentum();
0542     double x = tm.x();
0543     double y = tm.y();
0544     double z = tm.z();
0545     double m = ks.mass();
0546     double e = sqrt((x * x) + (y * y) + (z * z) + (m * m));
0547     totalMomentum.SetPxPyPzE(x, y, z, e);
0548   } else {
0549     edm::LogPrint("FitNotFound") << "BPHKinematicFit::fitMomentum: "
0550                                  << "simple momentum sum computed";
0551     math::XYZTLorentzVector tm;
0552     const vector<const reco::Candidate*>& daug = daughters();
0553     int n = daug.size();
0554     while (n--)
0555       tm += daug[n]->p4();
0556     const vector<BPHRecoConstCandPtr>& comp = daughComp();
0557     int m = comp.size();
0558     while (m--)
0559       tm += comp[m]->p4();
0560     totalMomentum = tm;
0561   }
0562   oldMom = false;
0563   return;
0564 }