Back to home page

Project CMSSW displayed by LXR

 
 

    


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

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/SpecificDecay/interface/BPHDecayToTkpTknSymChargeBuilder.h"
0012 
0013 //-------------------------------
0014 // Collaborating Class Headers --
0015 //-------------------------------
0016 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHRecoBuilder.h"
0017 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHPlusMinusCandidate.h"
0018 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHTrackReference.h"
0019 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHParticlePtSelect.h"
0020 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHParticleEtaSelect.h"
0021 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHMassSelect.h"
0022 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHMassSymSelect.h"
0023 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHChi2Select.h"
0024 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHParticleMasses.h"
0025 #include "DataFormats/Candidate/interface/Candidate.h"
0026 #include "DataFormats/Math/interface/LorentzVector.h"
0027 
0028 //---------------
0029 // C++ Headers --
0030 //---------------
0031 using namespace std;
0032 
0033 //-------------------
0034 // Initializations --
0035 //-------------------
0036 
0037 //----------------
0038 // Constructors --
0039 //----------------
0040 BPHDecayToTkpTknSymChargeBuilder::BPHDecayToTkpTknSymChargeBuilder(
0041     const edm::EventSetup& es,
0042     const std::string& daug1Name,
0043     double daug1Mass,
0044     double daug1Sigma,
0045     const std::string& daug2Name,
0046     double daug2Mass,
0047     double daug2Sigma,
0048     const BPHRecoBuilder::BPHGenericCollection* posCollection,
0049     const BPHRecoBuilder::BPHGenericCollection* negCollection,
0050     double expectedMass)
0051     : BPHDecayGenericBuilder(es),
0052       d1Name(daug1Name),
0053       d1Mass(daug1Mass),
0054       d1Sigma(daug1Sigma),
0055       d2Name(daug2Name),
0056       d2Mass(daug2Mass),
0057       d2Sigma(daug2Sigma),
0058       eMass(expectedMass),
0059       pCollection(posCollection),
0060       nCollection(negCollection),
0061       pt1Min(-1.0),
0062       pt2Min(-1.0),
0063       eta1Max(10.0),
0064       eta2Max(10.0),
0065       dzMax(1.0) {}
0066 
0067 //--------------
0068 // Destructor --
0069 //--------------
0070 BPHDecayToTkpTknSymChargeBuilder::~BPHDecayToTkpTknSymChargeBuilder() {}
0071 
0072 //--------------
0073 // Operations --
0074 //--------------
0075 vector<BPHPlusMinusConstCandPtr> BPHDecayToTkpTknSymChargeBuilder::build() {
0076   if (updated)
0077     return recList;
0078 
0079   recList.clear();
0080 
0081   // extract basic informations from input collections
0082 
0083   vector<Particle*> pList;
0084   vector<Particle*> nList;
0085 
0086   addParticle(pCollection, +1, pList);
0087   addParticle(nCollection, -1, nList);
0088   int iPos;
0089   int iNeg;
0090   int nPos = pList.size();
0091   int nNeg = nList.size();
0092   double massMin = getMassMin();
0093   double massMax = getMassMax();
0094   double mSqMin = massMin * massMin * 0.9;
0095   double mSqMax = massMax * massMax * 1.1;
0096   if (mSqMin < 0.0)
0097     mSqMin = 0.0;
0098 
0099   for (iPos = 0; iPos < nPos; ++iPos) {
0100     Particle* pc = pList[iPos];
0101     const reco::Track* pt = pc->track;
0102     double px = pc->px;
0103     double py = pc->py;
0104     double pz = pc->pz;
0105     double p1 = pc->e1;
0106     double p2 = pc->e2;
0107     for (iNeg = 0; iNeg < nNeg; ++iNeg) {
0108       Particle* nc = nList[iNeg];
0109       const reco::Track* nt = nc->track;
0110       if (fabs(nt->dz() - pt->dz()) > dzMax)
0111         continue;
0112       double nx = nc->px;
0113       double ny = nc->py;
0114       double nz = nc->pz;
0115       double n1 = nc->e1;
0116       double n2 = nc->e2;
0117       const float tx = px + nx;
0118       const float ty = py + ny;
0119       const float tz = pz + nz;
0120       const float ta = ((p1 > 0.0) && (n2 > 0.0) ? p1 + n2 : -1.0);
0121       const float tb = ((p2 > 0.0) && (n1 > 0.0) ? p2 + n1 : -1.0);
0122       float ma = (ta > 0 ? (ta * ta) - ((tx * tx) + (ty * ty) + (tz * tz)) : -1.0);
0123       float mb = (tb > 0 ? (tb * tb) - ((tx * tx) + (ty * ty) + (tz * tz)) : -1.0);
0124       //      if ( ma > 0.0 ) ma = sqrt( ma );
0125       //      else            ma = -1.0;
0126       //      if ( mb > 0.0 ) mb = sqrt( mb );
0127       //      else            mb = -1.0;
0128       if (((ma < mSqMin) || (ma > mSqMax)) && ((mb < mSqMin) || (mb > mSqMax)))
0129         continue;
0130       BPHPlusMinusCandidatePtr rc(nullptr);
0131       float rcMass = -1.0;
0132       if (ma > 0.0) {
0133         rc = BPHPlusMinusCandidateWrap::create(evSetup);
0134         rc->add(d1Name, pc->cand, d1Mass, d1Sigma);
0135         rc->add(d2Name, nc->cand, d2Mass, d2Sigma);
0136         rcMass = rc->composite().mass();
0137       }
0138       BPHPlusMinusCandidatePtr rb(nullptr);
0139       float rbMass = -1.0;
0140       if (mb > 0.0) {
0141         rb = BPHPlusMinusCandidateWrap::create(evSetup);
0142         rb->add(d1Name, nc->cand, d1Mass, d1Sigma);
0143         rb->add(d2Name, pc->cand, d2Mass, d2Sigma);
0144         rbMass = rb->composite().mass();
0145       }
0146       BPHPlusMinusCandidatePtr* rp(nullptr);
0147       double mass = -1.0;
0148       if ((rc.get() != nullptr) && ((rb.get() == nullptr) || (fabs(rcMass - eMass) < fabs(rbMass - eMass)))) {
0149         mass = rcMass;
0150         rp = &rc;
0151       } else {
0152         mass = rbMass;
0153         rp = &rb;
0154       }
0155       BPHPlusMinusCandidate* rr = rp->get();
0156       if (mass < massMin)
0157         continue;
0158       if (mass > massMax)
0159         continue;
0160       if (!chi2Sel->accept(*rr))
0161         continue;
0162       recList.push_back(*rp);
0163     }
0164   }
0165 
0166   for (iPos = 0; iPos < nPos; ++iPos)
0167     delete pList[iPos];
0168   for (iNeg = 0; iNeg < nNeg; ++iNeg)
0169     delete nList[iNeg];
0170 
0171   updated = true;
0172   return recList;
0173 }
0174 
0175 /// set cuts
0176 void BPHDecayToTkpTknSymChargeBuilder::setTrk1PtMin(double pt) {
0177   updated = false;
0178   pt1Min = pt;
0179   return;
0180 }
0181 
0182 void BPHDecayToTkpTknSymChargeBuilder::setTrk2PtMin(double pt) {
0183   updated = false;
0184   pt2Min = pt;
0185   return;
0186 }
0187 
0188 void BPHDecayToTkpTknSymChargeBuilder::setTrk1EtaMax(double eta) {
0189   updated = false;
0190   eta1Max = eta;
0191   return;
0192 }
0193 
0194 void BPHDecayToTkpTknSymChargeBuilder::setTrk2EtaMax(double eta) {
0195   updated = false;
0196   eta2Max = eta;
0197   return;
0198 }
0199 
0200 void BPHDecayToTkpTknSymChargeBuilder::setDzMax(double dz) {
0201   updated = false;
0202   dzMax = dz;
0203   return;
0204 }
0205 
0206 void BPHDecayToTkpTknSymChargeBuilder::addParticle(const BPHRecoBuilder::BPHGenericCollection* collection,
0207                                                    int charge,
0208                                                    vector<Particle*>& list) {
0209   int i;
0210   int n = collection->size();
0211   list.reserve(n);
0212   for (i = 0; i < n; ++i) {
0213     const reco::Candidate& cand = collection->get(i);
0214     int q = cand.charge();
0215     if ((charge > 0) && (q <= 0))
0216       continue;
0217     if ((charge < 0) && (q >= 0))
0218       continue;
0219     const reco::Candidate::LorentzVector p4 = cand.p4();
0220     const reco::Track* tk = BPHTrackReference::getTrack(cand, "cfhp");
0221     if (tk == nullptr)
0222       continue;
0223     double px = p4.px();
0224     double py = p4.py();
0225     double pz = p4.pz();
0226     double p2 = (px * px) + (py * py) + (pz * pz);
0227     double e1 = ((p4.pt() >= pt1Min) && (p4.eta() <= eta1Max) ? sqrt(p2 + (d1Mass * d1Mass)) : -1.0);
0228     double e2 = ((p4.pt() >= pt2Min) && (p4.eta() <= eta2Max) ? sqrt(p2 + (d2Mass * d2Mass)) : -1.0);
0229     if ((e1 > 0.0) || (e2 > 0.0))
0230       list.push_back(new Particle(&cand, tk, px, py, pz, e1, e2));
0231   }
0232   return;
0233 }