Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-04-07 05:50:32 

0001 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0002 #include "DataFormats/PatCandidates/interface/Muon.h"
0003 #include "FWCore/Common/interface/TriggerNames.h"
0004 #include "FWCore/Framework/interface/ESHandle.h"
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHParticleMasses.h"
0008 #include "HeavyFlavorAnalysis/SpecificDecay/plugins/BPHHistoSpecificDecay.h"
0009 #include "RecoVertex/VertexTools/interface/VertexDistanceXY.h"
0010 
0011 #include "TFile.h"
0012 #include "TH1.h"
0013 #include "TMath.h"
0014 #include "TTree.h"
0015 #include "TVector3.h"
0016 
0017 #include "Math/VectorUtil.h"
0018 #include <set>
0019 #include <string>
0020 #include <iostream>
0021 #include <fstream>
0022 #include <cmath>
0023 
0024 using namespace std;
0025 
0026 #define SET_LABEL(NAME, PSET) (NAME = PSET.getParameter<string>(#NAME))
0027 // SET_LABEL(xyz,ps);
0028 // is equivalent to
0029 // xyz = ps.getParameter<string>( "xyx" );
0030 
0031 #define CAT3(A, B, C) A##B##C
0032 #define STRING_NX(A) #A
0033 #define STRING(A) STRING_NX(A)
0034 #define CHK_TRIG(RESULTS, NAMES, INDEX, PATH)                 \
0035   if (NAMES[INDEX].find(STRING(CAT3(HLT_, PATH, _v))) == 0) { \
0036     flag_##PATH = RESULTS->accept(INDEX);                     \
0037     continue;                                                 \
0038   }
0039 // CHK_TRIG( trigResults, names, i, xyz );
0040 // is equivalent to
0041 // if ( names[i].find( "HLT_xyz_v" ) == 0 ) { flag_xyz = trigResults->accept( i ); break; }
0042 
0043 class VertexAnalysis {
0044 public:
0045   static double cAlpha(const reco::Vertex* pvtx, const reco::Vertex* svtx, float px, float py) {
0046     TVector3 disp(svtx->x() - pvtx->x(), svtx->y() - pvtx->y(), 0);
0047     TVector3 cmom(px, py, 0);
0048     return disp.Dot(cmom) / (disp.Perp() * cmom.Perp());
0049   }
0050   static double cAlpha(const reco::Vertex* pvtx, const reco::Vertex* svtx, const TVector3& cmom) {
0051     TVector3 disp(svtx->x() - pvtx->x(), svtx->y() - pvtx->y(), 0);
0052     return disp.Dot(cmom) / (disp.Perp() * cmom.Perp());
0053   }
0054   static void dist2D(const reco::Vertex* pvtx,
0055                      const reco::Vertex* svtx,
0056                      float px,
0057                      float py,
0058                      float mass,
0059                      double& ctauPV,
0060                      double& ctauErrPV) {
0061     dist2D(pvtx, svtx, px, py, cAlpha(pvtx, svtx, px, py), mass, ctauPV, ctauErrPV);
0062     return;
0063   }
0064   static void dist2D(const reco::Vertex* pvtx,
0065                      const reco::Vertex* svtx,
0066                      float px,
0067                      float py,
0068                      double cosAlpha,
0069                      float mass,
0070                      double& ctauPV,
0071                      double& ctauErrPV) {
0072     TVector3 cmom(px, py, 0);
0073     AlgebraicVector3 vmom(px, py, 0);
0074     VertexDistanceXY vdistXY;
0075     Measurement1D distXY = vdistXY.distance(*svtx, *pvtx);
0076     ctauPV = distXY.value() * cosAlpha * mass / cmom.Perp();
0077     GlobalError sve = svtx->error();
0078     GlobalError pve = pvtx->error();
0079     AlgebraicSymMatrix33 vXYe = sve.matrix() + pve.matrix();
0080     ctauErrPV = sqrt(ROOT::Math::Similarity(vmom, vXYe)) * mass / cmom.Perp2();
0081     return;
0082   }
0083 };
0084 
0085 class BPHUserData {
0086 public:
0087   template <class T>
0088   static const T* get(const pat::CompositeCandidate& cand, const string& name) {
0089     if (cand.hasUserData(name))
0090       return cand.userData<T>(name);
0091     return nullptr;
0092   }
0093   static float get(const pat::CompositeCandidate& cand, const string& name, float d = 0.0) {
0094     if (cand.hasUserFloat(name))
0095       return cand.userFloat(name);
0096     return d;
0097   }
0098   template <class T>
0099   static const T* getByRef(const pat::CompositeCandidate& cand, const string& name) {
0100     if (cand.hasUserData(name)) {
0101       typedef edm::Ref<std::vector<T>> objRef;
0102       const objRef* ref = cand.userData<objRef>(name);
0103       if (ref == nullptr)
0104         return nullptr;
0105       if (ref->isNull())
0106         return nullptr;
0107       return ref->get();
0108     }
0109     return nullptr;
0110   }
0111 };
0112 
0113 class BPHDaughters {
0114 public:
0115   static vector<const reco::Candidate*> get(const pat::CompositeCandidate& cand, float massMin, float massMax) {
0116     int i;
0117     int n = cand.numberOfDaughters();
0118     vector<const reco::Candidate*> v;
0119     v.reserve(n);
0120     for (i = 0; i < n; ++i) {
0121       const reco::Candidate* dptr = cand.daughter(i);
0122       float mass = dptr->mass();
0123       if ((mass > massMin) && (mass < massMax))
0124         v.push_back(dptr);
0125     }
0126     return v;
0127   }
0128 };
0129 
0130 class BPHSoftMuonSelect {
0131 public:
0132   BPHSoftMuonSelect(int cutTrackerLayers = 5,
0133                     int cutPixelLayers = 0,
0134                     float maxDxy = 0.3,
0135                     float maxDz = 20.0,
0136                     bool goodMuon = true,
0137                     bool highPurity = true)
0138       : cutTL(cutTrackerLayers), cutPL(cutPixelLayers), maxXY(maxDxy), maxZ(maxDz), gM(goodMuon), hP(highPurity) {}
0139   ~BPHSoftMuonSelect() {}
0140   bool accept(const reco::Candidate& cand, const reco::Vertex* pv) const {
0141     const pat::Muon* p = dynamic_cast<const pat::Muon*>(&cand);
0142     if (p == nullptr)
0143       return false;
0144     if (gM && !muon::isGoodMuon(*p, muon::TMOneStationTight))
0145       return false;
0146     if (p->innerTrack()->hitPattern().trackerLayersWithMeasurement() <= cutTL)
0147       return false;
0148     if (p->innerTrack()->hitPattern().pixelLayersWithMeasurement() <= cutPL)
0149       return false;
0150     if (hP && !p->innerTrack()->quality(reco::TrackBase::highPurity))
0151       return false;
0152     if (pv == nullptr)
0153       return true;
0154     const reco::Vertex::Point& pos = pv->position();
0155     if (fabs(p->innerTrack()->dxy(pos)) >= maxXY)
0156       return false;
0157     if (fabs(p->innerTrack()->dz(pos)) >= maxZ)
0158       return false;
0159     return true;
0160   }
0161 
0162 private:
0163   const reco::Vertex* pv;
0164   int cutTL;
0165   int cutPL;
0166   float maxXY;
0167   float maxZ;
0168   bool gM;
0169   bool hP;
0170 };
0171 
0172 class BPHDaughterSelect : public BPHHistoSpecificDecay::CandidateSelect {
0173 public:
0174   BPHDaughterSelect(float ptMinLoose,
0175                     float ptMinTight,
0176                     float etaMaxLoose,
0177                     float etaMaxTight,
0178                     const BPHSoftMuonSelect* softMuonselector = nullptr)
0179       : pLMin(ptMinLoose), pTMin(ptMinTight), eLMax(etaMaxLoose), eTMax(etaMaxTight), sms(softMuonselector) {}
0180   ~BPHDaughterSelect() override {}
0181   bool accept(const pat::CompositeCandidate& cand, const reco::Vertex* pv = nullptr) const override {
0182     return accept(cand, pLMin, pTMin, eLMax, eTMax, pv, sms);
0183   }
0184   static bool accept(const pat::CompositeCandidate& cand,
0185                      float ptMinLoose,
0186                      float ptMinTight,
0187                      float etaMaxLoose,
0188                      float etaMaxTight,
0189                      const reco::Vertex* pv = nullptr,
0190                      const BPHSoftMuonSelect* softMuonselector = nullptr) {
0191     const reco::Candidate* dptr0 = cand.daughter(0);
0192     const reco::Candidate* dptr1 = cand.daughter(1);
0193     if (dptr0 == nullptr)
0194       return false;
0195     if (dptr1 == nullptr)
0196       return false;
0197     float pt0 = dptr0->pt();
0198     float pt1 = dptr1->pt();
0199     if ((pt0 < ptMinLoose) || (pt1 < ptMinLoose))
0200       return false;
0201     if ((pt0 < ptMinTight) && (pt1 < ptMinTight))
0202       return false;
0203     float eta0 = fabs(dptr0->eta());
0204     float eta1 = fabs(dptr1->eta());
0205     if ((etaMaxLoose > 0) && ((eta0 > etaMaxLoose) || (eta1 > etaMaxLoose)))
0206       return false;
0207     if ((etaMaxTight > 0) && ((eta0 > etaMaxTight) && (eta1 > etaMaxTight)))
0208       return false;
0209     if (softMuonselector != nullptr) {
0210       const reco::Vertex* pvtx = BPHUserData::getByRef<reco::Vertex>(cand, "primaryVertex");
0211       if (pvtx == nullptr)
0212         return false;
0213       if (!softMuonselector->accept(*dptr0, pvtx))
0214         return false;
0215       if (!softMuonselector->accept(*dptr1, pvtx))
0216         return false;
0217     }
0218     return true;
0219   }
0220 
0221 private:
0222   float pLMin;
0223   float pTMin;
0224   float eLMax;
0225   float eTMax;
0226   const BPHSoftMuonSelect* sms;
0227 };
0228 
0229 class BPHCompositeBasicSelect : public BPHHistoSpecificDecay::CandidateSelect {
0230 public:
0231   BPHCompositeBasicSelect(
0232       float massMin, float massMax, float ptMin = -1.0, float etaMax = -1.0, float rapidityMax = -1.0)
0233       : mMin(massMin), mMax(massMax), pMin(ptMin), eMax(etaMax), yMax(rapidityMax) {}
0234   ~BPHCompositeBasicSelect() override {}
0235   bool accept(const pat::CompositeCandidate& cand, const reco::Vertex* pv = nullptr) const override {
0236     if (((mMin > 0) && (mMax < 0)) || ((mMin < 0) && (mMax > 0)) || ((mMin > 0) && (mMax > 0) && (mMin < mMax))) {
0237       float mass = cand.mass();
0238       if (mass < mMin)
0239         return false;
0240       if ((mMax > 0) && (mass > mMax))
0241         return false;
0242     }
0243     if (cand.pt() < pMin)
0244       return false;
0245     if ((eMax > 0) && (fabs(cand.eta()) > eMax))
0246       return false;
0247     if ((yMax > 0) && (fabs(cand.rapidity()) > yMax))
0248       return false;
0249     return true;
0250   }
0251 
0252 private:
0253   float mMin;
0254   float mMax;
0255   float pMin;
0256   float eMax;
0257   float yMax;
0258 };
0259 
0260 class BPHFittedBasicSelect : public BPHHistoSpecificDecay::CandidateSelect {
0261 public:
0262   BPHFittedBasicSelect(float massMin, float massMax, float ptMin = -1.0, float etaMax = -1.0, float rapidityMax = -1.0)
0263       : mMin(massMin), mMax(massMax), pMin(ptMin), eMax(etaMax), yMax(rapidityMax) {}
0264   ~BPHFittedBasicSelect() override {}
0265   bool accept(const pat::CompositeCandidate& cand, const reco::Vertex* pv = nullptr) const override {
0266     if (!cand.hasUserFloat("fitMass"))
0267       return false;
0268     float mass = cand.userFloat("fitMass");
0269     if (((mMin > 0) && (mMax < 0)) || ((mMin < 0) && (mMax > 0)) || ((mMin > 0) && (mMax > 0) && (mMin < mMax))) {
0270       if (mass < mMin)
0271         return false;
0272       if ((mMax > 0) && (mass > mMax))
0273         return false;
0274     }
0275     const Vector3DBase<float, GlobalTag>* fmom = BPHUserData::get<Vector3DBase<float, GlobalTag>>(cand, "fitMomentum");
0276     if (fmom == nullptr)
0277       return false;
0278     if (pMin > 0) {
0279       if (fmom->transverse() < pMin)
0280         return false;
0281     }
0282     if (eMax > 0) {
0283       if (fabs(fmom->eta()) > eMax)
0284         return false;
0285     }
0286     if (yMax > 0) {
0287       float x = fmom->x();
0288       float y = fmom->y();
0289       float z = fmom->z();
0290       float e = sqrt((x * x) + (y * y) + (z * z) + (mass * mass));
0291       float r = log((e + z) / (e - z)) / 2;
0292       if (fabs(r) > yMax)
0293         return false;
0294     }
0295     return true;
0296   }
0297 
0298 private:
0299   float mMin;
0300   float mMax;
0301   float pMin;
0302   float eMax;
0303   float yMax;
0304 };
0305 
0306 class BPHGenericVertexSelect : public BPHHistoSpecificDecay::CandidateSelect {
0307 public:
0308   BPHGenericVertexSelect(char vType, float probMin, float cosMin = -2.0, float sigMin = -1.0, char dMode = 'r')
0309       : type(vType), pMin(probMin), cMin(cosMin), sMin(sigMin), mode(dMode) {}
0310   ~BPHGenericVertexSelect() override {}
0311   bool accept(const pat::CompositeCandidate& cand, const reco::Vertex* pvtx) const override {
0312     if (pvtx == nullptr)
0313       return false;
0314     const reco::Vertex* svtx = nullptr;
0315     float px;
0316     float py;
0317     float mass;
0318     switch (type) {
0319       case 'c':
0320         svtx = BPHUserData::get<reco::Vertex>(cand, "vertex");
0321         px = cand.px();
0322         py = cand.py();
0323         mass = cand.mass();
0324         break;
0325       case 'f':
0326         svtx = BPHUserData::get<reco::Vertex>(cand, "fitVertex");
0327         {
0328           const Vector3DBase<float, GlobalTag>* fmom =
0329               BPHUserData::get<Vector3DBase<float, GlobalTag>>(cand, "fitMomentum");
0330           if (fmom == nullptr)
0331             return false;
0332           px = fmom->x();
0333           py = fmom->y();
0334         }
0335         if (!cand.hasUserFloat("fitMass"))
0336           return false;
0337         mass = cand.userFloat("fitMass");
0338         break;
0339       default:
0340         return false;
0341     }
0342     if (svtx == nullptr)
0343       return false;
0344     if (pMin > 0) {
0345       if (ChiSquaredProbability(svtx->chi2(), svtx->ndof()) < pMin)
0346         return false;
0347     }
0348     if ((cMin > -1.0) || (sMin > 0)) {
0349       float cosAlpha = VertexAnalysis::cAlpha(pvtx, svtx, px, py);
0350       if (cosAlpha < cMin)
0351         return false;
0352       if (sMin < 0)
0353         return true;
0354       double ctauPV;
0355       double ctauErrPV;
0356       VertexAnalysis::dist2D(pvtx, svtx, px, py, cosAlpha, mass, ctauPV, ctauErrPV);
0357       float dTest;
0358       switch (mode) {
0359         case 'a':
0360         case 'd':
0361           dTest = ctauPV;
0362           break;
0363         case 'r':
0364         case 's':
0365         default:
0366           dTest = ctauPV / ctauErrPV;
0367           break;
0368       }
0369       if (dTest < sMin)
0370         return false;
0371     }
0372     return true;
0373   }
0374 
0375 private:
0376   char type;
0377   float pMin;
0378   float cMin;
0379   float sMin;
0380   char mode;
0381 };
0382 
0383 BPHHistoSpecificDecay::BPHHistoSpecificDecay(const edm::ParameterSet& ps) {
0384   useTrig = (!SET_LABEL(trigResultsLabel, ps).empty());
0385   useOnia = (!SET_LABEL(oniaCandsLabel, ps).empty());
0386   useSd = (!SET_LABEL(sdCandsLabel, ps).empty());
0387   useSs = (!SET_LABEL(ssCandsLabel, ps).empty());
0388   useBu = (!SET_LABEL(buCandsLabel, ps).empty());
0389   useBd = (!SET_LABEL(bdCandsLabel, ps).empty());
0390   useBs = (!SET_LABEL(bsCandsLabel, ps).empty());
0391   useK0 = (!SET_LABEL(k0CandsLabel, ps).empty());
0392   useL0 = (!SET_LABEL(l0CandsLabel, ps).empty());
0393   useB0 = (!SET_LABEL(b0CandsLabel, ps).empty());
0394   useLb = (!SET_LABEL(lbCandsLabel, ps).empty());
0395   useBc = (!SET_LABEL(bcCandsLabel, ps).empty());
0396   useX3872 = (!SET_LABEL(x3872CandsLabel, ps).empty());
0397   if (useTrig)
0398     consume<edm::TriggerResults>(trigResultsToken, trigResultsLabel);
0399   if (useOnia)
0400     consume<vector<pat::CompositeCandidate>>(oniaCandsToken, oniaCandsLabel);
0401   if (useSd)
0402     consume<vector<pat::CompositeCandidate>>(sdCandsToken, sdCandsLabel);
0403   if (useSs)
0404     consume<vector<pat::CompositeCandidate>>(ssCandsToken, ssCandsLabel);
0405   if (useBu)
0406     consume<vector<pat::CompositeCandidate>>(buCandsToken, buCandsLabel);
0407   if (useBd)
0408     consume<vector<pat::CompositeCandidate>>(bdCandsToken, bdCandsLabel);
0409   if (useBs)
0410     consume<vector<pat::CompositeCandidate>>(bsCandsToken, bsCandsLabel);
0411   if (useK0)
0412     consume<vector<pat::CompositeCandidate>>(k0CandsToken, k0CandsLabel);
0413   if (useL0)
0414     consume<vector<pat::CompositeCandidate>>(l0CandsToken, l0CandsLabel);
0415   if (useB0)
0416     consume<vector<pat::CompositeCandidate>>(b0CandsToken, b0CandsLabel);
0417   if (useLb)
0418     consume<vector<pat::CompositeCandidate>>(lbCandsToken, lbCandsLabel);
0419   if (useBc)
0420     consume<vector<pat::CompositeCandidate>>(bcCandsToken, bcCandsLabel);
0421   if (useX3872)
0422     consume<vector<pat::CompositeCandidate>>(x3872CandsToken, x3872CandsLabel);
0423 
0424   static const BPHSoftMuonSelect* sms = new BPHSoftMuonSelect;
0425 
0426   //////////// quarkonia selection ////////////
0427 
0428   double phiIMassMin = 0.85;
0429   double phiIMassMax = 1.20;
0430   double phiIPtMin = 18.0;
0431   double phiIEtaMax = -1.0;
0432   double phiIYMax = -1.0;
0433   double phiBMassMin = 0.85;
0434   double phiBMassMax = 1.20;
0435   double phiBPtMin = 14.0;
0436   double phiBEtaMax = -1.0;
0437   double phiBYMax = 1.25;
0438   double jPsiIMassMin = 2.80;
0439   double jPsiIMassMax = 3.40;
0440   double jPsiIPtMin = 25.0;
0441   double jPsiIEtaMax = -1.0;
0442   double jPsiIYMax = -1.0;
0443   double jPsiBMassMin = 2.80;
0444   double jPsiBMassMax = 3.40;
0445   double jPsiBPtMin = 20.0;
0446   double jPsiBEtaMax = -1.0;
0447   double jPsiBYMax = 1.25;
0448   double psi2IMassMin = 3.40;
0449   double psi2IMassMax = 4.00;
0450   double psi2IPtMin = 18.0;
0451   double psi2IEtaMax = -1.0;
0452   double psi2IYMax = -1.0;
0453   double psi2BMassMin = 3.40;
0454   double psi2BMassMax = 4.00;
0455   double psi2BPtMin = 10.0;
0456   double psi2BEtaMax = -1.0;
0457   double psi2BYMax = 1.25;
0458   double upsIMassMin = 8.50;
0459   double upsIMassMax = 11.0;
0460   double upsIPtMin = 15.0;
0461   double upsIEtaMax = -1.0;
0462   double upsIYMax = -1.0;
0463   double upsBMassMin = 8.50;
0464   double upsBMassMax = 11.0;
0465   double upsBPtMin = 12.0;
0466   //  double  upsBEtaMax  =  1.5; // 2017
0467   //  double  upsBYMax    = -1.0; // 2017
0468   double upsBEtaMax = -1.0;  // 2018
0469   double upsBYMax = 1.4;     // 2018
0470 
0471   double oniaProbMin = 0.005;
0472   double oniaCosMin = -2.0;
0473   double oniaSigMin = -1.0;
0474 
0475   double oniaMuPtMinLoose = 2.0;
0476   double oniaMuPtMinTight = -1.0;
0477   double oniaMuEtaMaxLoose = -1.0;
0478   double oniaMuEtaMaxTight = -1.0;
0479 
0480   phiIBasicSelect = new BPHCompositeBasicSelect(phiIMassMin, phiIMassMax, phiIPtMin, phiIEtaMax, phiIYMax);
0481   phiBBasicSelect = new BPHCompositeBasicSelect(phiBMassMin, phiBMassMax, phiBPtMin, phiBEtaMax, phiBYMax);
0482   jPsiIBasicSelect = new BPHCompositeBasicSelect(jPsiIMassMin, jPsiIMassMax, jPsiIPtMin, jPsiIEtaMax, jPsiIYMax);
0483   jPsiBBasicSelect = new BPHCompositeBasicSelect(jPsiBMassMin, jPsiBMassMax, jPsiBPtMin, jPsiBEtaMax, jPsiBYMax);
0484   psi2IBasicSelect = new BPHCompositeBasicSelect(psi2IMassMin, psi2IMassMax, psi2IPtMin, psi2IEtaMax, psi2IYMax);
0485   psi2BBasicSelect = new BPHCompositeBasicSelect(psi2BMassMin, psi2BMassMax, psi2BPtMin, psi2BEtaMax, psi2BYMax);
0486   upsIBasicSelect = new BPHCompositeBasicSelect(upsIMassMin, upsIMassMax, upsIPtMin, upsIEtaMax, upsIYMax);
0487   upsBBasicSelect = new BPHCompositeBasicSelect(upsBMassMin, upsBMassMax, upsBPtMin, upsBEtaMax, upsBYMax);
0488   //  oniaVertexSelect   = new BPHCompositeVertexSelect(
0489   oniaVertexSelect = new BPHGenericVertexSelect('c', oniaProbMin, oniaCosMin, oniaSigMin);
0490   oniaDaughterSelect =
0491       new BPHDaughterSelect(oniaMuPtMinLoose, oniaMuPtMinTight, oniaMuEtaMaxLoose, oniaMuEtaMaxTight, sms);
0492 
0493   double npJPsiMassMin = BPHParticleMasses::jPsiMass - 0.150;
0494   double npJPsiMassMax = BPHParticleMasses::jPsiMass + 0.150;
0495   double npJPsiPtMin = 8.0;
0496   double npJPsiEtaMax = -1.0;
0497   double npJPsiYMax = -1.0;
0498   double npMuPtMinLoose = 4.0;
0499   double npMuPtMinTight = -1.0;
0500   double npMuEtaMaxLoose = 2.2;
0501   double npMuEtaMaxTight = -1.0;
0502 
0503   npJPsiBasicSelect = new BPHCompositeBasicSelect(npJPsiMassMin, npJPsiMassMax, npJPsiPtMin, npJPsiEtaMax, npJPsiYMax);
0504   npJPsiDaughterSelect = new BPHDaughterSelect(npMuPtMinLoose, npMuPtMinTight, npMuEtaMaxLoose, npMuEtaMaxTight, sms);
0505 
0506   //////////// Bu selection ////////////
0507 
0508   double buIMassMin = 0.0;
0509   double buIMassMax = 999999.0;
0510   double buIPtMin = 27.0;
0511   double buIEtaMax = -1.0;
0512   double buIYMax = -1.0;
0513   double buIJPsiMassMin = BPHParticleMasses::jPsiMass - 0.150;
0514   double buIJPsiMassMax = BPHParticleMasses::jPsiMass + 0.150;
0515   double buIJPsiPtMin = 25.0;
0516   double buIJPsiEtaMax = -1.0;
0517   double buIJPsiYMax = -1.0;
0518   double buIProbMin = 0.15;
0519   double buICosMin = -2.0;
0520   double buISigMin = -1.0;
0521   //  double buIMuPtMinLoose  = -1.0;
0522   //  double buIMuPtMinTight  = -1.0;
0523   //  double buIMuEtaMaxLoose = -1.0;
0524   //  double buIMuEtaMaxTight = -1.0;
0525 
0526   buIKPtMin = 2.0;
0527 
0528   buIBasicSelect = new BPHFittedBasicSelect(buIMassMin, buIMassMax, buIPtMin, buIEtaMax, buIYMax);
0529   buIJPsiBasicSelect =
0530       new BPHCompositeBasicSelect(buIJPsiMassMin, buIJPsiMassMax, buIJPsiPtMin, buIJPsiEtaMax, buIJPsiYMax);
0531   buIVertexSelect = new BPHGenericVertexSelect('f', buIProbMin, buICosMin, buISigMin);
0532   buIJPsiDaughterSelect = nullptr;
0533   //  buIJPsiDaughterSelect = new BPHDaughterSelect(
0534   //                              buIMuPtMinLoose , buIMuPtMinTight ,
0535   //                              buIMuEtaMaxLoose, buMuEtaMaxTight, sms );
0536 
0537   double buDMassMin = 0.0;
0538   double buDMassMax = 999999.0;
0539   double buDPtMin = 10.0;
0540   double buDEtaMax = -1.0;
0541   double buDYMax = -1.0;
0542   double buDJPsiMassMin = BPHParticleMasses::jPsiMass - 0.150;
0543   double buDJPsiMassMax = BPHParticleMasses::jPsiMass + 0.150;
0544   double buDJPsiPtMin = 8.0;
0545   double buDJPsiEtaMax = -1.0;
0546   double buDJPsiYMax = -1.0;
0547   double buDProbMin = 0.10;
0548   double buDCosMin = 0.99;
0549   double buDSigMin = 3.0;
0550   //  double buDMuPtMinLoose  = -1.0;
0551   //  double buDMuPtMinTight  = -1.0;
0552   //  double buDMuEtaMaxLoose = -1.0;
0553   //  double buDMuEtaMaxTight = -1.0;
0554 
0555   buDKPtMin = 1.6;
0556 
0557   buDBasicSelect = new BPHFittedBasicSelect(buDMassMin, buDMassMax, buDPtMin, buDEtaMax, buDYMax);
0558   buDJPsiBasicSelect =
0559       new BPHCompositeBasicSelect(buDJPsiMassMin, buDJPsiMassMax, buDJPsiPtMin, buDJPsiEtaMax, buDJPsiYMax);
0560   buDVertexSelect = new BPHGenericVertexSelect('f', buDProbMin, buDCosMin, buDSigMin);
0561   buDJPsiDaughterSelect = nullptr;
0562   //  buDJPsiDaughterSelect = new BPHDaughterSelect(
0563   //                              buDMuPtMinLoose , buDMuPtMinTight ,
0564   //                              buDMuEtaMaxLoose, buDMuEtaMaxTight, sms );
0565 
0566   //////////// Bd -> JPsi Kx0 selection ////////////
0567 
0568   double bdIMassMin = 0.0;
0569   double bdIMassMax = 999999.0;
0570   double bdIPtMin = 27.0;
0571   double bdIEtaMax = -1.0;
0572   double bdIYMax = -1.0;
0573   double bdIJPsiMassMin = BPHParticleMasses::jPsiMass - 0.150;
0574   double bdIJPsiMassMax = BPHParticleMasses::jPsiMass + 0.150;
0575   double bdIJPsiPtMin = 25.0;
0576   double bdIJPsiEtaMax = -1.0;
0577   double bdIJPsiYMax = -1.0;
0578   double bdIKx0MassMin = BPHParticleMasses::kx0Mass - 0.100;
0579   double bdIKx0MassMax = BPHParticleMasses::kx0Mass + 0.100;
0580   double bdIKx0PtMin = -1.0;
0581   double bdIKx0EtaMax = -1.0;
0582   double bdIKx0YMax = -1.0;
0583   double bdIProbMin = 0.15;
0584   double bdICosMin = -2.0;
0585   double bdISigMin = -1.0;
0586   //  double bdIMuPtMinLoose  =  -1.0;
0587   //  double bdIMuPtMinTight  =  -1.0;
0588   //  double bdIMuEtaMaxLoose =  -1.0;
0589   //  double bdIMuEtaMaxTight =  -1.0;
0590 
0591   bdIBasicSelect = new BPHFittedBasicSelect(bdIMassMin, bdIMassMax, bdIPtMin, bdIEtaMax, bdIYMax);
0592   bdIJPsiBasicSelect =
0593       new BPHCompositeBasicSelect(bdIJPsiMassMin, bdIJPsiMassMax, bdIJPsiPtMin, bdIJPsiEtaMax, bdIJPsiYMax);
0594   bdIKx0BasicSelect = new BPHCompositeBasicSelect(bdIKx0MassMin, bdIKx0MassMax, bdIKx0PtMin, bdIKx0EtaMax, bdIKx0YMax);
0595   bdIVertexSelect = new BPHGenericVertexSelect('f', bdIProbMin, bdICosMin, bdISigMin);
0596   bdIJPsiDaughterSelect = nullptr;
0597   //  bdIJPsiDaughterSelect = new BPHDaughterSelect(
0598   //                              bdIMuPtMinLoose , bdIMuPtMinTight ,
0599   //                              bdIMuEtaMaxLoose, bdIMuEtaMaxTight, sms );
0600 
0601   double bdDMassMin = 0.0;
0602   double bdDMassMax = 999999.0;
0603   double bdDPtMin = 10.0;
0604   double bdDEtaMax = -1.0;
0605   double bdDYMax = -1.0;
0606   double bdDJPsiMassMin = BPHParticleMasses::jPsiMass - 0.150;
0607   double bdDJPsiMassMax = BPHParticleMasses::jPsiMass + 0.150;
0608   double bdDJPsiPtMin = 8.0;
0609   double bdDJPsiEtaMax = -1.0;
0610   double bdDJPsiYMax = -1.0;
0611   double bdDKx0MassMin = BPHParticleMasses::kx0Mass - 0.100;
0612   double bdDKx0MassMax = BPHParticleMasses::kx0Mass + 0.100;
0613   double bdDKx0PtMin = -1.0;
0614   double bdDKx0EtaMax = -1.0;
0615   double bdDKx0YMax = -1.0;
0616   double bdDProbMin = 0.10;
0617   double bdDCosMin = 0.99;
0618   double bdDSigMin = 3.0;
0619   //  double bdDMuPtMinLoose  = -1.0;
0620   //  double bdDMuPtMinTight  = -1.0;
0621   //  double bdDMuEtaMaxLoose = -1.0;
0622   //  double bdDMuEtaMaxTight = -1.0;
0623 
0624   bdDBasicSelect = new BPHFittedBasicSelect(bdDMassMin, bdDMassMax, bdDPtMin, bdDEtaMax, bdDYMax);
0625   bdDJPsiBasicSelect =
0626       new BPHCompositeBasicSelect(bdDJPsiMassMin, bdDJPsiMassMax, bdDJPsiPtMin, bdDJPsiEtaMax, bdDJPsiYMax);
0627   bdDKx0BasicSelect = new BPHCompositeBasicSelect(bdDKx0MassMin, bdDKx0MassMax, bdDKx0PtMin, bdDKx0EtaMax, bdDKx0YMax);
0628   bdDVertexSelect = new BPHGenericVertexSelect('f', bdDProbMin, bdDCosMin, bdDSigMin);
0629   bdDJPsiDaughterSelect = nullptr;
0630   //  bdDJPsiDaughterSelect = new BPHDaughterSelect(
0631   //                              bdDMuPtMinLoose , bdDMuPtMinTight ,
0632   //                              bdDMuEtaMaxLoose, bdDMuEtaMaxTight, sms );
0633 
0634   //////////// Bs selection ////////////
0635 
0636   double bsIMassMin = 0.0;
0637   double bsIMassMax = 999999.0;
0638   double bsIPtMin = 27.0;
0639   double bsIEtaMax = -1.0;
0640   double bsIYMax = -1.0;
0641   double bsIJPsiMassMin = BPHParticleMasses::jPsiMass - 0.150;
0642   double bsIJPsiMassMax = BPHParticleMasses::jPsiMass + 0.150;
0643   double bsIJPsiPtMin = 25.0;
0644   double bsIJPsiEtaMax = -1.0;
0645   double bsIJPsiYMax = -1.0;
0646   double bsIPhiMassMin = BPHParticleMasses::phiMass - 0.010;
0647   double bsIPhiMassMax = BPHParticleMasses::phiMass + 0.010;
0648   double bsIPhiPtMin = -1.0;
0649   double bsIPhiEtaMax = -1.0;
0650   double bsIPhiYMax = -1.0;
0651   double bsIProbMin = 0.15;
0652   double bsICosMin = -2.0;
0653   double bsISigMin = -1.0;
0654   //  double bsIMuPtMinLoose  = -1.0;
0655   //  double bsIMuPtMinTight  = -1.0;
0656   //  double bsIMuEtaMaxLoose = -1.0;
0657   //  double bsIMuEtaMaxTight = -1.0;
0658 
0659   bsIBasicSelect = new BPHFittedBasicSelect(bsIMassMin, bsIMassMax, bsIPtMin, bsIEtaMax, bsIYMax);
0660   bsIJPsiBasicSelect =
0661       new BPHCompositeBasicSelect(bsIJPsiMassMin, bsIJPsiMassMax, bsIJPsiPtMin, bsIJPsiEtaMax, bsIJPsiYMax);
0662   bsIPhiBasicSelect = new BPHCompositeBasicSelect(bsIPhiMassMin, bsIPhiMassMax, bsIPhiPtMin, bsIPhiEtaMax, bsIPhiYMax);
0663   bsIVertexSelect = new BPHGenericVertexSelect('f', bsIProbMin, bsICosMin, bsISigMin);
0664   bsIJPsiDaughterSelect = nullptr;
0665   //  bsIJPsiDaughterSelect = new BPHDaughterSelect(
0666   //                              bsIMuPtMinLoose , bsIMuPtMinTight ,
0667   //                              bsIMuEtaMaxLoose, bsIMuEtaMaxTight, sms );
0668 
0669   double bsDMassMin = 0.0;
0670   double bsDMassMax = 999999.0;
0671   double bsDPtMin = 10.0;
0672   double bsDEtaMax = -1.0;
0673   double bsDYMax = -1.0;
0674   double bsDJPsiMassMin = BPHParticleMasses::jPsiMass - 0.150;
0675   double bsDJPsiMassMax = BPHParticleMasses::jPsiMass + 0.150;
0676   double bsDJPsiPtMin = 8.0;
0677   double bsDJPsiEtaMax = -1.0;
0678   double bsDJPsiYMax = -1.0;
0679   double bsDPhiMassMin = BPHParticleMasses::phiMass - 0.010;
0680   double bsDPhiMassMax = BPHParticleMasses::phiMass + 0.010;
0681   double bsDPhiPtMin = -1.0;
0682   double bsDPhiEtaMax = -1.0;
0683   double bsDPhiYMax = -1.0;
0684   double bsDProbMin = 0.10;
0685   double bsDCosMin = 0.99;
0686   double bsDSigMin = 3.0;
0687   //  double bsDMuPtMinLoose  = -1.0;
0688   //  double bsDMuPtMinTight  = -1.0;
0689   //  double bsDMuEtaMaxLoose = -1.0;
0690   //  double bsDMuEtaMaxTight = -1.0;
0691 
0692   bsDBasicSelect = new BPHFittedBasicSelect(bsDMassMin, bsDMassMax, bsDPtMin, bsDEtaMax, bsDYMax);
0693   bsDJPsiBasicSelect =
0694       new BPHCompositeBasicSelect(bsDJPsiMassMin, bsDJPsiMassMax, bsDJPsiPtMin, bsDJPsiEtaMax, bsDJPsiYMax);
0695   bsDPhiBasicSelect = new BPHCompositeBasicSelect(bsDPhiMassMin, bsDPhiMassMax, bsDPhiPtMin, bsDPhiEtaMax, bsDPhiYMax);
0696   bsDVertexSelect = new BPHGenericVertexSelect('f', bsDProbMin, bsDCosMin, bsDSigMin);
0697   bsDJPsiDaughterSelect = nullptr;
0698   //  bsDJPsiDaughterSelect = new BPHDaughterSelect(
0699   //                              bsDMuPtMinLoose , bsDMuPtMinTight ,
0700   //                              bsDMuEtaMaxLoose, bsDMuEtaMaxTight, sms );
0701 
0702   //////////// Bd -> JPsi K0s selection ////////////
0703 
0704   double b0IMassMin = 0.0;
0705   double b0IMassMax = 999999.0;
0706   double b0IPtMin = 27.0;
0707   double b0IEtaMax = -1.0;
0708   double b0IYMax = -1.0;
0709   double b0IJPsiMassMin = BPHParticleMasses::jPsiMass - 0.150;
0710   double b0IJPsiMassMax = BPHParticleMasses::jPsiMass + 0.150;
0711   double b0IJPsiPtMin = 25.0;
0712   double b0IJPsiEtaMax = -1.0;
0713   double b0IJPsiYMax = -1.0;
0714   double b0IK0sMassMin = BPHParticleMasses::k0sMass - 0.010;
0715   double b0IK0sMassMax = BPHParticleMasses::k0sMass + 0.010;
0716   double b0IK0sPtMin = -1.0;
0717   double b0IK0sEtaMax = -1.0;
0718   double b0IK0sYMax = -1.0;
0719   double b0IProbMin = 0.15;
0720   double b0ICosMin = -2.0;
0721   double b0ISigMin = -1.0;
0722   //  double b0IMuPtMinLoose  =  -1.0;
0723   //  double b0IMuPtMinTight  =  -1.0;
0724   //  double b0IMuEtaMaxLoose =  -1.0;
0725   //  double b0IMuEtaMaxTight =  -1.0;
0726 
0727   b0IBasicSelect = new BPHFittedBasicSelect(b0IMassMin, b0IMassMax, b0IPtMin, b0IEtaMax, b0IYMax);
0728   b0IJPsiBasicSelect =
0729       new BPHCompositeBasicSelect(b0IJPsiMassMin, b0IJPsiMassMax, b0IJPsiPtMin, b0IJPsiEtaMax, b0IJPsiYMax);
0730   b0IK0sBasicSelect = new BPHFittedBasicSelect(b0IK0sMassMin, b0IK0sMassMax, b0IK0sPtMin, b0IK0sEtaMax, b0IK0sYMax);
0731   b0IVertexSelect = new BPHGenericVertexSelect('f', b0IProbMin, b0ICosMin, b0ISigMin);
0732   b0IJPsiDaughterSelect = nullptr;
0733   //  b0IJPsiDaughterSelect = new BPHDaughterSelect(
0734   //                              b0IMuPtMinLoose , b0IMuPtMinTight ,
0735   //                              b0IMuEtaMaxLoose, b0IMuEtaMaxTight, sms );
0736 
0737   double b0DMassMin = 0.0;
0738   double b0DMassMax = 999999.0;
0739   double b0DPtMin = 10.0;
0740   double b0DEtaMax = -1.0;
0741   double b0DYMax = -1.0;
0742   double b0DJPsiMassMin = BPHParticleMasses::jPsiMass - 0.150;
0743   double b0DJPsiMassMax = BPHParticleMasses::jPsiMass + 0.150;
0744   double b0DJPsiPtMin = 8.0;
0745   double b0DJPsiEtaMax = -1.0;
0746   double b0DJPsiYMax = -1.0;
0747   double b0DK0sMassMin = BPHParticleMasses::k0sMass - 0.010;
0748   double b0DK0sMassMax = BPHParticleMasses::k0sMass + 0.010;
0749   double b0DK0sPtMin = -1.0;
0750   double b0DK0sEtaMax = -1.0;
0751   double b0DK0sYMax = -1.0;
0752   double b0DProbMin = 0.10;
0753   double b0DCosMin = 0.99;
0754   double b0DSigMin = 3.0;
0755   //  double b0DMuPtMinLoose  = -1.0;
0756   //  double b0DMuPtMinTight  = -1.0;
0757   //  double b0DMuEtaMaxLoose = -1.0;
0758   //  double b0DMuEtaMaxTight = -1.0;
0759 
0760   b0DBasicSelect = new BPHFittedBasicSelect(b0DMassMin, b0DMassMax, b0DPtMin, b0DEtaMax, b0DYMax);
0761   b0DJPsiBasicSelect =
0762       new BPHCompositeBasicSelect(b0DJPsiMassMin, b0DJPsiMassMax, b0DJPsiPtMin, b0DJPsiEtaMax, b0DJPsiYMax);
0763   b0DK0sBasicSelect = new BPHFittedBasicSelect(b0DK0sMassMin, b0DK0sMassMax, b0DK0sPtMin, b0DK0sEtaMax, b0DK0sYMax);
0764   b0DVertexSelect = new BPHGenericVertexSelect('f', b0DProbMin, b0DCosMin, b0DSigMin);
0765   b0DJPsiDaughterSelect = nullptr;
0766   //  b0DJPsiDaughterSelect = new BPHDaughterSelect(
0767   //                              b0DMuPtMinLoose , b0DMuPtMinTight ,
0768   //                              b0DMuEtaMaxLoose, b0DMuEtaMaxTight, sms );
0769 
0770   //////////// Lambdab -> JPsi Lambda0 selection ////////////
0771 
0772   double lbIMassMin = 0.0;
0773   double lbIMassMax = 999999.0;
0774   double lbIPtMin = 27.0;
0775   double lbIEtaMax = -1.0;
0776   double lbIYMax = -1.0;
0777   double lbIJPsiMassMin = BPHParticleMasses::jPsiMass - 0.150;
0778   double lbIJPsiMassMax = BPHParticleMasses::jPsiMass + 0.150;
0779   double lbIJPsiPtMin = 25.0;
0780   double lbIJPsiEtaMax = -1.0;
0781   double lbIJPsiYMax = -1.0;
0782   double lbILambda0MassMin = BPHParticleMasses::lambda0Mass - 0.006;
0783   double lbILambda0MassMax = BPHParticleMasses::lambda0Mass + 0.006;
0784   double lbILambda0PtMin = -1.0;
0785   double lbILambda0EtaMax = -1.0;
0786   double lbILambda0YMax = -1.0;
0787   double lbIProbMin = 0.10;
0788   double lbICosMin = -2.0;
0789   double lbISigMin = -1.0;
0790   //  double lbIMuPtMinLoose   =  -1.0;
0791   //  double lbIMuPtMinTight   =  -1.0;
0792   //  double lbIMuEtaMaxLoose  =  -1.0;
0793   //  double lbIMuEtaMaxTight  =  -1.0;
0794 
0795   lbIBasicSelect = new BPHFittedBasicSelect(lbIMassMin, lbIMassMax, lbIPtMin, lbIEtaMax, lbIYMax);
0796   lbIJPsiBasicSelect =
0797       new BPHCompositeBasicSelect(lbIJPsiMassMin, lbIJPsiMassMax, lbIJPsiPtMin, lbIJPsiEtaMax, lbIJPsiYMax);
0798   lbILambda0BasicSelect =
0799       new BPHFittedBasicSelect(lbILambda0MassMin, lbILambda0MassMax, lbILambda0PtMin, lbILambda0EtaMax, lbILambda0YMax);
0800   lbIVertexSelect = new BPHGenericVertexSelect('f', lbIProbMin, lbICosMin, lbISigMin);
0801   lbIJPsiDaughterSelect = nullptr;
0802   //  lbIJPsiDaughterSelect = new BPHDaughterSelect(
0803   //                              lbIMuPtMinLoose , lbIMuPtMinTight ,
0804   //                              lbIMuEtaMaxLoose, lbIMuEtaMaxTight, sms );
0805 
0806   double lbDMassMin = 0.0;
0807   double lbDMassMax = 999999.0;
0808   double lbDPtMin = 10.0;
0809   double lbDEtaMax = -1.0;
0810   double lbDYMax = -1.0;
0811   double lbDJPsiMassMin = BPHParticleMasses::jPsiMass - 0.150;
0812   double lbDJPsiMassMax = BPHParticleMasses::jPsiMass + 0.150;
0813   double lbDJPsiPtMin = 8.0;
0814   double lbDJPsiEtaMax = -1.0;
0815   double lbDJPsiYMax = -1.0;
0816   double lbDLambda0MassMin = BPHParticleMasses::lambda0Mass - 0.006;
0817   double lbDLambda0MassMax = BPHParticleMasses::lambda0Mass + 0.006;
0818   double lbDLambda0PtMin = -1.0;
0819   double lbDLambda0EtaMax = -1.0;
0820   double lbDLambda0YMax = -1.0;
0821   double lbDProbMin = 0.10;
0822   double lbDCosMin = 0.99;
0823   double lbDSigMin = 3.0;
0824   //  double lbDMuPtMinLoose   = -1.0;
0825   //  double lbDMuPtMinTight   = -1.0;
0826   //  double lbDMuEtaMaxLoose  = -1.0;
0827   //  double lbDMuEtaMaxTight  = -1.0;
0828 
0829   lbDBasicSelect = new BPHFittedBasicSelect(lbDMassMin, lbDMassMax, lbDPtMin, lbDEtaMax, lbDYMax);
0830   lbDJPsiBasicSelect =
0831       new BPHCompositeBasicSelect(lbDJPsiMassMin, lbDJPsiMassMax, lbDJPsiPtMin, lbDJPsiEtaMax, lbDJPsiYMax);
0832   lbDLambda0BasicSelect =
0833       new BPHFittedBasicSelect(lbDLambda0MassMin, lbDLambda0MassMax, lbDLambda0PtMin, lbDLambda0EtaMax, lbDLambda0YMax);
0834   lbDVertexSelect = new BPHGenericVertexSelect('f', lbDProbMin, lbDCosMin, lbDSigMin);
0835   lbDJPsiDaughterSelect = nullptr;
0836   //  lbDJPsiDaughterSelect = new BPHDaughterSelect(
0837   //                              lbDMuPtMinLoose , lbDMuPtMinTight ,
0838   //                              lbDMuEtaMaxLoose, lbDMuEtaMaxTight, sms );
0839 
0840   //////////// Bc selection ////////////
0841 
0842   double bcIMassMin = 0.0;
0843   double bcIMassMax = 999999.0;
0844   double bcIPtMin = 27.0;
0845   double bcIEtaMax = -1.0;
0846   double bcIYMax = -1.0;
0847   double bcIJPsiMassMin = BPHParticleMasses::jPsiMass - 0.150;
0848   double bcIJPsiMassMax = BPHParticleMasses::jPsiMass + 0.150;
0849   double bcIJPsiPtMin = 25.0;
0850   double bcIJPsiEtaMax = -1.0;
0851   double bcIJPsiYMax = -1.0;
0852   double bcIJPsiProbMin = 0.005;
0853   double bcIProbMin = 0.10;
0854   double bcICosMin = -2.0;
0855   double bcISigMin = -1.0;
0856   double bcIDistMin = 0.01;
0857   //  double bcIMuPtMinLoose  = -1.0;
0858   //  double bcIMuPtMinTight  = -1.0;
0859   //  double bcIMuEtaMaxLoose = -1.0;
0860   //  double bcIMuEtaMaxTight = -1.0;
0861 
0862   bcIPiPtMin = 3.5;
0863 
0864   bcIBasicSelect = new BPHFittedBasicSelect(bcIMassMin, bcIMassMax, bcIPtMin, bcIEtaMax, bcIYMax);
0865   bcIJPsiBasicSelect =
0866       new BPHCompositeBasicSelect(bcIJPsiMassMin, bcIJPsiMassMax, bcIJPsiPtMin, bcIJPsiEtaMax, bcIJPsiYMax);
0867   bcIJPsiVertexSelect = new BPHGenericVertexSelect('c', bcIJPsiProbMin);
0868   bcIVertexSelect = new BPHGenericVertexSelect('f', bcIProbMin, bcICosMin, bcISigMin, bcIDistMin);
0869   bcIJPsiDaughterSelect = nullptr;
0870   //  bcIJPsiDaughterSelect = new BPHDaughterSelect(
0871   //                              bcIMuPtMinLoose , bcIMuPtMinTight ,
0872   //                              bcIMuEtaMaxLoose, bcMuEtaMaxTight, sms );
0873 
0874   double bcDMassMin = 0.0;
0875   double bcDMassMax = 999999.0;
0876   double bcDPtMin = 8.0;
0877   double bcDEtaMax = -1.0;
0878   double bcDYMax = -1.0;
0879   double bcDJPsiMassMin = BPHParticleMasses::jPsiMass - 0.150;
0880   double bcDJPsiMassMax = BPHParticleMasses::jPsiMass + 0.150;
0881   double bcDJPsiPtMin = 7.0;
0882   double bcDJPsiEtaMax = -1.0;
0883   double bcDJPsiYMax = -1.0;
0884   double bcDJPsiProbMin = 0.005;
0885   double bcDProbMin = 0.10;
0886   double bcDCosMin = 0.99;
0887   double bcDSigMin = 3.0;
0888   //  double bcDMuPtMinLoose  = -1.0;
0889   //  double bcDMuPtMinTight  = -1.0;
0890   //  double bcDMuEtaMaxLoose = -1.0;
0891   //  double bcDMuEtaMaxTight = -1.0;
0892 
0893   bcJPsiDcaMax = 0.5;
0894   bcDPiPtMin = 3.5;
0895 
0896   bcDBasicSelect = new BPHFittedBasicSelect(bcDMassMin, bcDMassMax, bcDPtMin, bcDEtaMax, bcDYMax);
0897   bcDJPsiBasicSelect =
0898       new BPHCompositeBasicSelect(bcDJPsiMassMin, bcDJPsiMassMax, bcDJPsiPtMin, bcDJPsiEtaMax, bcDJPsiYMax);
0899   bcDJPsiVertexSelect = new BPHGenericVertexSelect('c', bcDJPsiProbMin);
0900   bcDVertexSelect = new BPHGenericVertexSelect('f', bcDProbMin, bcDCosMin, bcDSigMin);
0901   bcDJPsiDaughterSelect = nullptr;
0902   //  bcDJPsiDaughterSelect = new BPHDaughterSelect(
0903   //                              bcDMuPtMinLoose , bcDMuPtMinTight ,
0904   //                              bcDMuEtaMaxLoose, bcDMuEtaMaxTight, sms );
0905 
0906   //////////// X3872 selection ////////////
0907 
0908   double x3872IMassMin = 0.0;
0909   double x3872IMassMax = 999999.0;
0910   double x3872IPtMin = 27.0;
0911   double x3872IEtaMax = -1.0;
0912   double x3872IYMax = -1.0;
0913   double x3872IJPsiMassMin = BPHParticleMasses::jPsiMass - 0.150;
0914   double x3872IJPsiMassMax = BPHParticleMasses::jPsiMass + 0.150;
0915   double x3872IJPsiPtMin = 25.0;
0916   double x3872IJPsiEtaMax = -1.0;
0917   double x3872IJPsiYMax = -1.0;
0918   double x3872IJPsiProbMin = 0.10;
0919   double x3872IProbMin = 0.10;
0920   double x3872ICosMin = -2.0;
0921   double x3872ISigMin = -1.0;
0922   double x3872IDistMin = 0.01;
0923   //  double x3872IMuPtMinLoose  = -1.0;
0924   //  double x3872IMuPtMinTight  = -1.0;
0925   //  double x3872IMuEtaMaxLoose = -1.0;
0926   //  double x3872IMuEtaMaxTight = -1.0;
0927 
0928   x3872JPsiDcaMax = 0.5;
0929   x3872IPiPtMin = 1.2;
0930 
0931   x3872IBasicSelect = new BPHFittedBasicSelect(x3872IMassMin, x3872IMassMax, x3872IPtMin, x3872IEtaMax, x3872IYMax);
0932   x3872IJPsiBasicSelect = new BPHCompositeBasicSelect(
0933       x3872IJPsiMassMin, x3872IJPsiMassMax, x3872IJPsiPtMin, x3872IJPsiEtaMax, x3872IJPsiYMax);
0934   x3872IJPsiVertexSelect = new BPHGenericVertexSelect('c', x3872IJPsiProbMin);
0935   x3872IVertexSelect = new BPHGenericVertexSelect('f', x3872IProbMin, x3872ICosMin, x3872ISigMin, x3872IDistMin);
0936   x3872IJPsiDaughterSelect = nullptr;
0937   //  x3872IJPsiDaughterSelect = new BPHDaughterSelect(
0938   //                                 x3872IMuPtMinLoose , x3872IMuPtMinTight,
0939   //                                 x3872IMuEtaMaxLoose, x3872MuEtaMaxTight,
0940   //                                 sms );
0941 
0942   double x3872DMassMin = 0.0;
0943   double x3872DMassMax = 999999.0;
0944   double x3872DPtMin = 8.0;
0945   double x3872DEtaMax = -1.0;
0946   double x3872DYMax = -1.0;
0947   double x3872DJPsiMassMin = BPHParticleMasses::jPsiMass - 0.150;
0948   double x3872DJPsiMassMax = BPHParticleMasses::jPsiMass + 0.150;
0949   double x3872DJPsiPtMin = 7.0;
0950   double x3872DJPsiEtaMax = -1.0;
0951   double x3872DJPsiYMax = -1.0;
0952   double x3872DJPsiProbMin = 0.10;
0953   double x3872DProbMin = 0.10;
0954   double x3872DCosMin = 0.99;
0955   double x3872DSigMin = 3.0;
0956   //  double x3872DMuPtMinLoose  = -1.0;
0957   //  double x3872DMuPtMinTight  = -1.0;
0958   //  double x3872DMuEtaMaxLoose = -1.0;
0959   //  double x3872DMuEtaMaxTight = -1.0;
0960 
0961   x3872DPiPtMin = 1.2;
0962 
0963   x3872DBasicSelect = new BPHFittedBasicSelect(x3872DMassMin, x3872DMassMax, x3872DPtMin, x3872DEtaMax, x3872DYMax);
0964   x3872DJPsiBasicSelect = new BPHCompositeBasicSelect(
0965       x3872DJPsiMassMin, x3872DJPsiMassMax, x3872DJPsiPtMin, x3872DJPsiEtaMax, x3872DJPsiYMax);
0966   x3872DJPsiVertexSelect = new BPHGenericVertexSelect('c', x3872DJPsiProbMin);
0967   x3872DVertexSelect = new BPHGenericVertexSelect('f', x3872DProbMin, x3872DCosMin, x3872DSigMin);
0968   x3872DJPsiDaughterSelect = nullptr;
0969   //  x3872DJPsiDaughterSelect = new BPHDaughterSelect(
0970   //                                 x3872DMuPtMinLoose , x3872DMuPtMinTight ,
0971   //                                 x3872DMuEtaMaxLoose, x3872DMuEtaMaxTight,,
0972   //                                 sms );
0973 }
0974 
0975 BPHHistoSpecificDecay::~BPHHistoSpecificDecay() {
0976   delete phiIBasicSelect;
0977   delete phiBBasicSelect;
0978   delete jPsiIBasicSelect;
0979   delete jPsiBBasicSelect;
0980   delete psi2IBasicSelect;
0981   delete psi2BBasicSelect;
0982   delete upsIBasicSelect;
0983   delete upsBBasicSelect;
0984   delete oniaVertexSelect;
0985   delete oniaDaughterSelect;
0986 
0987   delete npJPsiBasicSelect;
0988   delete npJPsiDaughterSelect;
0989 
0990   delete buIBasicSelect;
0991   delete buIJPsiBasicSelect;
0992   delete buIVertexSelect;
0993   delete buIJPsiDaughterSelect;
0994   delete buDBasicSelect;
0995   delete buDJPsiBasicSelect;
0996   delete buDVertexSelect;
0997   delete buDJPsiDaughterSelect;
0998 
0999   delete bdIBasicSelect;
1000   delete bdIJPsiBasicSelect;
1001   delete bdIKx0BasicSelect;
1002   delete bdIVertexSelect;
1003   delete bdIJPsiDaughterSelect;
1004   delete bdDBasicSelect;
1005   delete bdDJPsiBasicSelect;
1006   delete bdDKx0BasicSelect;
1007   delete bdDVertexSelect;
1008   delete bdDJPsiDaughterSelect;
1009 
1010   delete bsIBasicSelect;
1011   delete bsIJPsiBasicSelect;
1012   delete bsIPhiBasicSelect;
1013   delete bsIVertexSelect;
1014   delete bsIJPsiDaughterSelect;
1015   delete bsDBasicSelect;
1016   delete bsDJPsiBasicSelect;
1017   delete bsDPhiBasicSelect;
1018   delete bsDVertexSelect;
1019   delete bsDJPsiDaughterSelect;
1020 
1021   delete b0IBasicSelect;
1022   delete b0IJPsiBasicSelect;
1023   delete b0IK0sBasicSelect;
1024   delete b0IVertexSelect;
1025   delete b0IJPsiDaughterSelect;
1026   delete b0DBasicSelect;
1027   delete b0DJPsiBasicSelect;
1028   delete b0DK0sBasicSelect;
1029   delete b0DVertexSelect;
1030   delete b0DJPsiDaughterSelect;
1031 
1032   delete lbIBasicSelect;
1033   delete lbIJPsiBasicSelect;
1034   delete lbILambda0BasicSelect;
1035   delete lbIVertexSelect;
1036   delete lbIJPsiDaughterSelect;
1037   delete lbDBasicSelect;
1038   delete lbDJPsiBasicSelect;
1039   delete lbDLambda0BasicSelect;
1040   delete lbDVertexSelect;
1041   delete lbDJPsiDaughterSelect;
1042 
1043   delete bcIBasicSelect;
1044   delete bcIJPsiBasicSelect;
1045   delete bcIJPsiVertexSelect;
1046   delete bcIVertexSelect;
1047   delete bcIJPsiDaughterSelect;
1048   delete bcDBasicSelect;
1049   delete bcDJPsiBasicSelect;
1050   delete bcDJPsiVertexSelect;
1051   delete bcDVertexSelect;
1052   delete bcDJPsiDaughterSelect;
1053 
1054   delete x3872IBasicSelect;
1055   delete x3872IJPsiBasicSelect;
1056   delete x3872IJPsiVertexSelect;
1057   delete x3872IVertexSelect;
1058   delete x3872IJPsiDaughterSelect;
1059   delete x3872DBasicSelect;
1060   delete x3872DJPsiBasicSelect;
1061   delete x3872DJPsiVertexSelect;
1062   delete x3872DVertexSelect;
1063   delete x3872DJPsiDaughterSelect;
1064 }
1065 
1066 void BPHHistoSpecificDecay::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1067   edm::ParameterSetDescription desc;
1068   desc.add<string>("trigResultsLabel", "");
1069   desc.add<string>("oniaCandsLabel", "");
1070   desc.add<string>("sdCandsLabel", "");
1071   desc.add<string>("ssCandsLabel", "");
1072   desc.add<string>("buCandsLabel", "");
1073   desc.add<string>("bdCandsLabel", "");
1074   desc.add<string>("bsCandsLabel", "");
1075   desc.add<string>("k0CandsLabel", "");
1076   desc.add<string>("l0CandsLabel", "");
1077   desc.add<string>("b0CandsLabel", "");
1078   desc.add<string>("lbCandsLabel", "");
1079   desc.add<string>("bcCandsLabel", "");
1080   desc.add<string>("x3872CandsLabel", "");
1081   descriptions.add("process.bphHistoSpecificDecay", desc);
1082   return;
1083 }
1084 
1085 void BPHHistoSpecificDecay::beginJob() {
1086   createHisto("massDIPhi", 50, 0.90, 1.15);      // Phi  mass inclusive
1087   createHisto("massTIPhi", 50, 0.90, 1.15);      // Phi  mass inclusive
1088   createHisto("massDBPhi", 50, 0.90, 1.15);      // Phi  mass barrel
1089   createHisto("massTBPhi", 50, 0.90, 1.15);      // Phi  mass barrel
1090   createHisto("massDIJPsi", 35, 2.95, 3.30);     // JPsi mass inclusive
1091   createHisto("massTIJPsi", 35, 2.95, 3.30);     // JPsi mass inclusive
1092   createHisto("massDBJPsi", 35, 2.95, 3.30);     // JPsi mass barrel
1093   createHisto("massTBJPsi", 35, 2.95, 3.30);     // JPsi mass barrel
1094   createHisto("massDIPsi2", 60, 3.40, 4.00);     // Psi2 mass inclusive
1095   createHisto("massTIPsi2", 60, 3.40, 4.00);     // Psi2 mass inclusive
1096   createHisto("massDBPsi2", 60, 3.40, 4.00);     // Psi2 mass barrel
1097   createHisto("massTBPsi2", 60, 3.40, 4.00);     // Psi2 mass barrel
1098   createHisto("massDIUps123", 115, 8.70, 11.0);  // Ups  mass inclusive
1099   createHisto("massTIUps123", 115, 8.70, 11.0);  // Ups  mass inclusive
1100   createHisto("massDBUps123", 115, 8.70, 11.0);  // Ups  mass barrel
1101   createHisto("massTBUps123", 115, 8.70, 11.0);  // Ups  mass barrel
1102   createHisto("massDIBu", 100, 5.00, 6.00);      // Bu   mass inclusive
1103   createHisto("massTIBu", 100, 5.00, 6.00);      // Bu   mass inclusive
1104   createHisto("massDDBu", 100, 5.00, 6.00);      // Bu   mass displaced
1105   createHisto("massTDBu", 100, 5.00, 6.00);      // Bu   mass displaced
1106   createHisto("massDIBd", 100, 5.00, 6.00);      // Bd   mass inclusive
1107   createHisto("massTIBd", 100, 5.00, 6.00);      // Bd   mass inclusive
1108   createHisto("massDDBd", 100, 5.00, 6.00);      // Bd   mass displaced
1109   createHisto("massTDBd", 100, 5.00, 6.00);      // Bd   mass displaced
1110   createHisto("massDIBs", 100, 5.00, 6.00);      // Bs   mass inclusive
1111   createHisto("massTIBs", 100, 5.00, 6.00);      // Bs   mass inclusive
1112   createHisto("massDDBs", 100, 5.00, 6.00);      // Bs   mass displaced
1113   createHisto("massTDBs", 100, 5.00, 6.00);      // Bs   mass displaced
1114   createHisto("massDIBc", 100, 6.00, 7.00);      // Bc   mass inclusive
1115   createHisto("massTIBc", 100, 6.00, 7.00);      // Bc   mass inclusive
1116   createHisto("massDDBc", 100, 6.00, 7.00);      // Bc   mass displaced
1117   createHisto("massTDBc", 100, 6.00, 7.00);      // Bc   mass displaced
1118   createHisto("massDIX3872", 40, 3.60, 4.00);    // X3872 mass inclusive
1119   createHisto("massTIX3872", 40, 3.60, 4.00);    // X3872 mass inclusive
1120   createHisto("massDDX3872", 40, 3.60, 4.00);    // X3872 mass displaced
1121   createHisto("massTDX3872", 40, 3.60, 4.00);    // X3872 mass displaced
1122   createHisto("mfitDIBu", 100, 5.00, 6.00);      // Bu   mass, with constr.
1123   createHisto("mfitTIBu", 100, 5.00, 6.00);      // Bu   mass, with constr.
1124   createHisto("mfitDDBu", 100, 5.00, 6.00);      // Bu   mass, with constr.
1125   createHisto("mfitTDBu", 100, 5.00, 6.00);      // Bu   mass, with constr.
1126   createHisto("mfitDIBd", 100, 5.00, 6.00);      // Bd   mass, with constr.
1127   createHisto("mfitTIBd", 100, 5.00, 6.00);      // Bd   mass, with constr.
1128   createHisto("mfitDDBd", 100, 5.00, 6.00);      // Bd   mass, with constr.
1129   createHisto("mfitTDBd", 100, 5.00, 6.00);      // Bd   mass, with constr.
1130   createHisto("mfitDIBs", 100, 5.00, 6.00);      // Bs   mass, with constr.
1131   createHisto("mfitTIBs", 100, 5.00, 6.00);      // Bs   mass, with constr.
1132   createHisto("mfitDDBs", 100, 5.00, 6.00);      // Bs   mass, with constr.
1133   createHisto("mfitTDBs", 100, 5.00, 6.00);      // Bs   mass, with constr.
1134   createHisto("mfitDIBc", 100, 6.00, 7.00);      // Bc   mass, with constr.
1135   createHisto("mfitTIBc", 100, 6.00, 7.00);      // Bc   mass, with constr.
1136   createHisto("mfitDDBc", 100, 6.00, 7.00);      // Bc   mass, with constr.
1137   createHisto("mfitTDBc", 100, 6.00, 7.00);      // Bc   mass, with constr.
1138   createHisto("mfitDIX3872", 40, 3.60, 4.00);    // X3872 mass, with constr.
1139   createHisto("mfitTIX3872", 40, 3.60, 4.00);    // X3872 mass, with constr.
1140   createHisto("mfitDDX3872", 40, 3.60, 4.00);    // X3872 mass, with constr.
1141   createHisto("mfitTDX3872", 40, 3.60, 4.00);    // X3872 mass, with constr.
1142   createHisto("massDIBuJPsi", 35, 2.95, 3.30);   // JPsi mass in Bu decay
1143   createHisto("massTIBuJPsi", 35, 2.95, 3.30);   // JPsi mass in Bu decay
1144   createHisto("massDDBuJPsi", 35, 2.95, 3.30);   // JPsi mass in Bu decay
1145   createHisto("massTDBuJPsi", 35, 2.95, 3.30);   // JPsi mass in Bu decay
1146   createHisto("massDIBdJPsi", 35, 2.95, 3.30);   // JPsi mass in Bd decay
1147   createHisto("massTIBdJPsi", 35, 2.95, 3.30);   // JPsi mass in Bd decay
1148   createHisto("massDDBdJPsi", 35, 2.95, 3.30);   // JPsi mass in Bd decay
1149   createHisto("massTDBdJPsi", 35, 2.95, 3.30);   // JPsi mass in Bd decay
1150   createHisto("massDIBdKx0", 50, 0.80, 1.05);    // Kx0  mass in Bd decay
1151   createHisto("massTIBdKx0", 50, 0.80, 1.05);    // Kx0  mass in Bd decay
1152   createHisto("massDDBdKx0", 50, 0.80, 1.05);    // Kx0  mass in Bd decay
1153   createHisto("massTDBdKx0", 50, 0.80, 1.05);    // Kx0  mass in Bd decay
1154   createHisto("massDIBsJPsi", 35, 2.95, 3.30);   // JPsi mass in Bs decay
1155   createHisto("massTIBsJPsi", 35, 2.95, 3.30);   // JPsi mass in Bs decay
1156   createHisto("massDDBsJPsi", 35, 2.95, 3.30);   // JPsi mass in Bs decay
1157   createHisto("massTDBsJPsi", 35, 2.95, 3.30);   // JPsi mass in Bs decay
1158   createHisto("massDIBsPhi", 50, 1.01, 1.03);    // Phi  mass in Bs decay
1159   createHisto("massTIBsPhi", 50, 1.01, 1.03);    // Phi  mass in Bs decay
1160   createHisto("massDDBsPhi", 50, 1.01, 1.03);    // Phi  mass in Bs decay
1161   createHisto("massTDBsPhi", 50, 1.01, 1.03);    // Phi  mass in Bs decay
1162   createHisto("massDIBcJPsi", 35, 2.95, 3.30);   // JPsi mass in Bc decay
1163   createHisto("massTIBcJPsi", 35, 2.95, 3.30);   // JPsi mass in Bc decay
1164   createHisto("massDDBcJPsi", 35, 2.95, 3.30);   // JPsi mass in Bc decay
1165   createHisto("massTDBcJPsi", 35, 2.95, 3.30);   // JPsi mass in Bc decay
1166   createHisto("massDIX3JPsi", 35, 2.95, 3.30);   // JPsi mass in X3872 decay
1167   createHisto("massTIX3JPsi", 35, 2.95, 3.30);   // JPsi mass in X3872 decay
1168   createHisto("massDDX3JPsi", 35, 2.95, 3.30);   // JPsi mass in X3872 decay
1169   createHisto("massTDX3JPsi", 35, 2.95, 3.30);   // JPsi mass in X3872 decay
1170   createHisto("massDK0s", 50, 0.40, 0.60);       // K0s  mass
1171   createHisto("mfitDK0s", 50, 0.40, 0.60);       // K0s  mass
1172   createHisto("massDLambda0", 60, 1.00, 1.30);   // Lambda0 mass
1173   createHisto("mfitDLambda0", 60, 1.00, 1.30);   // Lambda0 mass
1174   createHisto("massDIB0", 50, 5.00, 6.00);       // B0   mass inclusive
1175   createHisto("massTIB0", 50, 5.00, 6.00);       // B0   mass inclusive
1176   createHisto("massDDB0", 50, 5.00, 6.00);       // B0   mass displaced
1177   createHisto("massTDB0", 50, 5.00, 6.00);       // B0   mass displaced
1178   createHisto("mfitDIB0", 50, 5.00, 6.00);       // B0   mass, with constr.
1179   createHisto("mfitTIB0", 50, 5.00, 6.00);       // B0   mass, with constr.
1180   createHisto("mfitDDB0", 50, 5.00, 6.00);       // B0   mass, with constr.
1181   createHisto("mfitTDB0", 50, 5.00, 6.00);       // B0   mass, with constr.
1182   createHisto("massDIB0JPsi", 35, 2.95, 3.30);   // JPsi mass in B0 decay
1183   createHisto("massTIB0JPsi", 35, 2.95, 3.30);   // JPsi mass in B0 decay
1184   createHisto("massDDB0JPsi", 35, 2.95, 3.30);   // JPsi mass in B0 decay
1185   createHisto("massTDB0JPsi", 35, 2.95, 3.30);   // JPsi mass in B0 decay
1186   createHisto("massDIB0K0s", 50, 0.40, 0.60);    // K0s  mass in B0 decay
1187   createHisto("massTIB0K0s", 50, 0.40, 0.60);    // K0s  mass in B0 decay
1188   createHisto("massDDB0K0s", 50, 0.40, 0.60);    // K0s  mass in B0 decay
1189   createHisto("massTDB0K0s", 50, 0.40, 0.60);    // K0s  mass in B0 decay
1190   createHisto("mfitDIB0K0s", 50, 0.40, 0.60);    // K0s  mass in B0 decay
1191   createHisto("mfitTIB0K0s", 50, 0.40, 0.60);    // K0s  mass in B0 decay
1192   createHisto("mfitDDB0K0s", 50, 0.40, 0.60);    // K0s  mass in B0 decay
1193   createHisto("mfitTDB0K0s", 50, 0.40, 0.60);    // K0s  mass in B0 decay
1194   createHisto("massDILambdab", 25, 5.00, 6.00);  // Lambdab mass inclusive
1195   createHisto("massTILambdab", 25, 5.00, 6.00);  // Lambdab mass inclusive
1196   createHisto("massDDLambdab", 25, 5.00, 6.00);  // Lambdab mass displaced
1197   createHisto("massTDLambdab", 25, 5.00, 6.00);  // Lambdab mass displaced
1198   createHisto("mfitDILambdab", 25, 5.00, 6.00);  // Lambdab mass, with constr.
1199   createHisto("mfitTILambdab", 25, 5.00, 6.00);  // Lambdab mass, with constr.
1200   createHisto("mfitDDLambdab", 25, 5.00, 6.00);  // Lambdab mass, with constr.
1201   createHisto("mfitTDLambdab", 25, 5.00, 6.00);  // Lambdab mass, with constr.
1202   createHisto("massDILbJPsi", 35, 2.95, 3.30);   // JPsi mass in Lambdab decay
1203   createHisto("massTILbJPsi", 35, 2.95, 3.30);   // JPsi mass in Lambdab decay
1204   createHisto("massDDLbJPsi", 35, 2.95, 3.30);   // JPsi mass in Lambdab decay
1205   createHisto("massTDLbJPsi", 35, 2.95, 3.30);   // JPsi mass in Lambdab decay
1206   createHisto("massDILbL0", 60, 1.00, 1.30);     // L0   mass in Lambdab decay
1207   createHisto("massTILbL0", 60, 1.00, 1.30);     // L0   mass in Lambdab decay
1208   createHisto("massDDLbL0", 60, 1.00, 1.30);     // L0   mass in Lambdab decay
1209   createHisto("massTDLbL0", 60, 1.00, 1.30);     // L0   mass in Lambdab decay
1210   createHisto("mfitDILbL0", 60, 1.00, 1.30);     // L0   mass in Lambdab decay
1211   createHisto("mfitTILbL0", 60, 1.00, 1.30);     // L0   mass in Lambdab decay
1212   createHisto("mfitDDLbL0", 60, 1.00, 1.30);     // L0   mass in Lambdab decay
1213   createHisto("mfitTDLbL0", 60, 1.00, 1.30);     // L0   mass in Lambdab decay
1214 
1215   createHisto("massFull", 200, 2.00, 12.0);  // Full onia mass
1216 
1217   createHisto("ctauDIJPsi", 60, -0.05, 0.25);     // JPsi ctau inclusive
1218   createHisto("ctauTIJPsi", 60, -0.05, 0.25);     // JPsi ctau inclusive
1219   createHisto("ctauDBJPsi", 60, -0.05, 0.25);     // JPsi ctau barrel
1220   createHisto("ctauTBJPsi", 60, -0.05, 0.25);     // JPsi ctau barrel
1221   createHisto("ctauDIBu", 60, -0.05, 0.25);       // Bu   ctau inclusive
1222   createHisto("ctauTIBu", 60, -0.05, 0.25);       // Bu   ctau inclusive
1223   createHisto("ctauDDBu", 60, -0.05, 0.25);       // Bu   ctau displaced
1224   createHisto("ctauTDBu", 60, -0.05, 0.25);       // Bu   ctau displaced
1225   createHisto("ctauDIBd", 60, -0.05, 0.25);       // Bd   ctau inclusive
1226   createHisto("ctauTIBd", 60, -0.05, 0.25);       // Bd   ctau inclusive
1227   createHisto("ctauDDBd", 60, -0.05, 0.25);       // Bd   ctau displaced
1228   createHisto("ctauTDBd", 60, -0.05, 0.25);       // Bd   ctau displaced
1229   createHisto("ctauDIBs", 60, -0.05, 0.25);       // Bs   ctau inclusive
1230   createHisto("ctauTIBs", 60, -0.05, 0.25);       // Bs   ctau inclusive
1231   createHisto("ctauDDBs", 60, -0.05, 0.25);       // Bs   ctau displaced
1232   createHisto("ctauTDBs", 60, -0.05, 0.25);       // Bs   ctau displaced
1233   createHisto("ctauDIB0", 60, -0.05, 0.25);       // B0   ctau inclusive
1234   createHisto("ctauTIB0", 60, -0.05, 0.25);       // B0   ctau inclusive
1235   createHisto("ctauDDB0", 60, -0.05, 0.25);       // B0   ctau displaced
1236   createHisto("ctauTDB0", 60, -0.05, 0.25);       // B0   ctau displaced
1237   createHisto("ctauDILambdab", 60, -0.05, 0.25);  // Lambdab ctau inclusive
1238   createHisto("ctauTILambdab", 60, -0.05, 0.25);  // Lambdab ctau inclusive
1239   createHisto("ctauDDLambdab", 60, -0.05, 0.25);  // Lambdab ctau displaced
1240   createHisto("ctauTDLambdab", 60, -0.05, 0.25);  // Lambdab ctau displaced
1241 
1242   recoName = new string;
1243   tree = fs->make<TTree>("BPHReco", "BPHReco");
1244   b_runNumber = tree->Branch("runNumber", &runNumber, "runNumber/i");
1245   b_lumiSection = tree->Branch("lumiSection", &lumiSection, "lumiSection/i");
1246   b_eventNumber = tree->Branch("eventNumber", &eventNumber, "eventNumber/i");
1247   b_recoName = tree->Branch("recoName", &recoName, 8192, 99);
1248   b_recoMass = tree->Branch("recoMass", &recoMass, "recoMass/F");
1249   b_recoTime = tree->Branch("recoTime", &recoTime, "recoTime/F");
1250   b_recoErrT = tree->Branch("recoErrT", &recoErrT, "recoErrT/F");
1251 
1252   return;
1253 }
1254 
1255 void BPHHistoSpecificDecay::analyze(const edm::Event& ev, const edm::EventSetup& es) {
1256   //  if ( ev.id().run  () !=    316239 ) return;
1257   //  if ( ev.id().event() != 170736782 ) return;
1258   static map<string, ofstream*> ofMap;
1259   if (ofMap.empty()) {
1260     ofMap["BarPhi"] = nullptr;
1261     ofMap["IncJPsi"] = nullptr;
1262     ofMap["BarJPsi"] = nullptr;
1263     ofMap["IncPsi2"] = nullptr;
1264     ofMap["BarPsi2"] = nullptr;
1265     ofMap["BarUpsilon123"] = nullptr;
1266     ofMap["InclusiveBu"] = nullptr;
1267     ofMap["DisplacedBu"] = nullptr;
1268     ofMap["InclusiveBd"] = nullptr;
1269     ofMap["DisplacedBd"] = nullptr;
1270     ofMap["InclusiveBs"] = nullptr;
1271     ofMap["DisplacedBs"] = nullptr;
1272     ofMap["K0s"] = nullptr;
1273     ofMap["Lambda0"] = nullptr;
1274     ofMap["InclusiveB0"] = nullptr;
1275     ofMap["DisplacedB0"] = nullptr;
1276     ofMap["InclusiveLambdab"] = nullptr;
1277     ofMap["DisplacedLambdab"] = nullptr;
1278     ofMap["InclusiveBc"] = nullptr;
1279     ofMap["DisplacedBc"] = nullptr;
1280     ofMap["InclusiveX3872"] = nullptr;
1281     ofMap["DisplacedX3872"] = nullptr;
1282     map<string, ofstream*>::iterator iter = ofMap.begin();
1283     map<string, ofstream*>::iterator iend = ofMap.end();
1284     string name = "list";
1285     while (iter != iend) {
1286       iter->second = new ofstream(name + iter->first);
1287       ++iter;
1288     }
1289   }
1290 
1291   // event number
1292   runNumber = ev.id().run();
1293   lumiSection = ev.id().luminosityBlock();
1294   eventNumber = ev.id().event();
1295 
1296   // get object collections
1297   // collections are got through "BPHTokenWrapper" interface to allow
1298   // uniform access in different CMSSW versions
1299 
1300   //////////// trigger results ////////////
1301 
1302   edm::Handle<edm::TriggerResults> trigResults;
1303   const edm::TriggerNames* trigNames = nullptr;
1304   if (useTrig) {
1305     trigResultsToken.get(ev, trigResults);
1306     if (trigResults.isValid())
1307       trigNames = &ev.triggerNames(*trigResults);
1308   }
1309 
1310   bool flag_Dimuon25_Jpsi = false;
1311   bool flag_Dimuon20_Jpsi_Barrel_Seagulls = false;
1312   bool flag_Dimuon14_Phi_Barrel_Seagulls = false;
1313   bool flag_Dimuon18_PsiPrime = false;
1314   bool flag_Dimuon10_PsiPrime_Barrel_Seagulls = false;
1315   bool flag_Dimuon12_Upsilon_eta1p5 = false;
1316   bool flag_Dimuon12_Upsilon_y1p4 = false;
1317   bool flag_DoubleMu4_JpsiTrk_Displaced = false;
1318   if (trigNames != nullptr) {
1319     const edm::TriggerNames::Strings& names = trigNames->triggerNames();
1320     int iObj;
1321     int nObj = names.size();
1322     for (iObj = 0; iObj < nObj; ++iObj) {
1323       //      cout << names[iObj] << endl;
1324       CHK_TRIG(trigResults, names, iObj, Dimuon25_Jpsi)
1325       CHK_TRIG(trigResults, names, iObj, Dimuon20_Jpsi_Barrel_Seagulls)
1326       CHK_TRIG(trigResults, names, iObj, Dimuon14_Phi_Barrel_Seagulls)
1327       CHK_TRIG(trigResults, names, iObj, Dimuon18_PsiPrime)
1328       CHK_TRIG(trigResults, names, iObj, Dimuon10_PsiPrime_Barrel_Seagulls)
1329       CHK_TRIG(trigResults, names, iObj, Dimuon12_Upsilon_eta1p5)
1330       CHK_TRIG(trigResults, names, iObj, Dimuon12_Upsilon_y1p4)
1331       CHK_TRIG(trigResults, names, iObj, DoubleMu4_JpsiTrk_Displaced)
1332     }
1333   }
1334 
1335   //  cout <<       "Dimuon25_Jpsi "
1336   //       << ( flag_Dimuon25_Jpsi ? 'A' : 'R' ) << endl;
1337   //  cout <<       "Dimuon20_Jpsi_Barrel_Seagulls "
1338   //       << ( flag_Dimuon20_Jpsi_Barrel_Seagulls ? 'A' : 'R' ) << endl;
1339   //  cout <<       "Dimuon14_Phi_Barrel_Seagulls "
1340   //       << ( flag_Dimuon14_Phi_Barrel_Seagulls ? 'A' : 'R' ) << endl;
1341   //  cout <<       "Dimuon18_PsiPrime "
1342   //       << ( flag_Dimuon18_PsiPrime ? 'A' : 'R' ) << endl;
1343   //  cout <<       "Dimuon10_PsiPrime_Barrel_Seagulls "
1344   //       << ( flag_Dimuon10_PsiPrime_Barrel_Seagulls ? 'A' : 'R' ) << endl;
1345   //  cout <<       "Dimuon12_Upsilon_eta1p5 "
1346   //       << ( flag_Dimuon12_Upsilon_eta1p5 ? 'A' : 'R' ) << endl;
1347   //  cout <<       "DoubleMu4_JpsiTrk_Displaced "
1348   //       << ( flag_DoubleMu4_JpsiTrk_Displaced ? 'A' : 'R' ) << endl;
1349 
1350   //////////// quarkonia ////////////
1351 
1352   edm::Handle<vector<pat::CompositeCandidate>> oniaCands;
1353   int iqo;
1354   int nqo = 0;
1355   if (useOnia) {
1356     oniaCandsToken.get(ev, oniaCands);
1357     nqo = oniaCands->size();
1358   }
1359 
1360   for (iqo = 0; iqo < nqo; ++iqo) {
1361     LogTrace("DataDump") << "*********** quarkonium " << iqo << "/" << nqo << " ***********";
1362     const pat::CompositeCandidate& cand = oniaCands->at(iqo);
1363     if (!oniaVertexSelect->accept(cand, BPHUserData::getByRef<reco::Vertex>(cand, "primaryVertex")))
1364       continue;
1365     if (!oniaDaughterSelect->accept(cand))
1366       continue;
1367     fillHisto("Full", cand, 'c');
1368     if (phiBBasicSelect->accept(cand)) {
1369       fillHisto("DBPhi", cand, 'c');
1370       if (flag_Dimuon14_Phi_Barrel_Seagulls)
1371         fillHisto("TBPhi", cand, 'c');
1372       if (flag_Dimuon14_Phi_Barrel_Seagulls)
1373         *ofMap["BarPhi"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1374                          << cand.mass() << endl;
1375     }
1376     if (jPsiIBasicSelect->accept(cand)) {
1377       fillHisto("DIJPsi", cand, 'c');
1378       if (flag_Dimuon25_Jpsi)
1379         fillHisto("TIJPsi", cand, 'c');
1380       if (flag_Dimuon25_Jpsi)
1381         *ofMap["IncJPsi"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1382                           << cand.mass() << endl;
1383     }
1384     if (jPsiBBasicSelect->accept(cand)) {
1385       fillHisto("DBJPsi", cand, 'c');
1386       if (flag_Dimuon20_Jpsi_Barrel_Seagulls)
1387         fillHisto("TBJPsi", cand, 'c');
1388       if (flag_Dimuon20_Jpsi_Barrel_Seagulls)
1389         *ofMap["BarJPsi"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1390                           << cand.mass() << endl;
1391     }
1392     if (psi2IBasicSelect->accept(cand)) {
1393       fillHisto("DIPsi2", cand, 'c');
1394       if (flag_Dimuon18_PsiPrime)
1395         fillHisto("TIPsi2", cand, 'c');
1396       if (flag_Dimuon18_PsiPrime)
1397         *ofMap["IncPsi2"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1398                           << cand.mass() << endl;
1399     }
1400     if (psi2BBasicSelect->accept(cand)) {
1401       fillHisto("DBPsi2", cand, 'c');
1402       if (flag_Dimuon10_PsiPrime_Barrel_Seagulls)
1403         fillHisto("TBPsi2", cand, 'c');
1404       if (flag_Dimuon10_PsiPrime_Barrel_Seagulls)
1405         *ofMap["BarPsi2"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1406                           << cand.mass() << endl;
1407     }
1408     if (upsBBasicSelect->accept(cand)) {
1409       fillHisto("DBUps123", cand, 'c');
1410       if (flag_Dimuon12_Upsilon_eta1p5 || flag_Dimuon12_Upsilon_y1p4)
1411         fillHisto("TBUps123", cand, 'c');
1412       if (flag_Dimuon12_Upsilon_eta1p5 || flag_Dimuon12_Upsilon_y1p4)
1413         *ofMap["BarUpsilon123"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1414                                 << cand.mass() << endl;
1415     }
1416   }
1417 
1418   //////////// Bu ////////////
1419 
1420   edm::Handle<vector<pat::CompositeCandidate>> buCands;
1421   int ibu;
1422   int nbu = 0;
1423   if (useBu) {
1424     buCandsToken.get(ev, buCands);
1425     nbu = buCands->size();
1426   }
1427 
1428   for (ibu = 0; ibu < nbu; ++ibu) {
1429     LogTrace("DataDump") << "*********** Bu " << ibu << "/" << nbu << " ***********";
1430     const pat::CompositeCandidate& cand = buCands->at(ibu);
1431     const pat::CompositeCandidate* jPsi = BPHUserData::getByRef<pat::CompositeCandidate>(cand, "refToJPsi");
1432     LogTrace("DataDump") << "JPsi: " << jPsi;
1433     if (jPsi == nullptr)
1434       continue;
1435     if (!npJPsiBasicSelect->accept(*jPsi))
1436       continue;
1437     if (!npJPsiDaughterSelect->accept(*jPsi))
1438       continue;
1439     const reco::Candidate* kptr = BPHDaughters::get(cand, 0.49, 0.50).front();
1440     if (kptr == nullptr)
1441       continue;
1442     if (buIBasicSelect->accept(cand) && buIJPsiBasicSelect->accept(*jPsi) &&
1443         //         buIJPsiDaughterSelect->accept( *jPsi ) &&
1444         buIVertexSelect->accept(cand, BPHUserData::getByRef<reco::Vertex>(*jPsi, "primaryVertex")) &&
1445         (kptr->pt() > buIKPtMin)) {
1446       fillHisto("DIBu", cand, 'f');
1447       fillHisto("DIBuJPsi", *jPsi, 'c');
1448       if (flag_Dimuon25_Jpsi) {
1449         fillHisto("TIBu", cand, 'f');
1450         fillHisto("TIBuJPsi", *jPsi, 'c');
1451         *ofMap["InclusiveBu"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1452                               << (cand.hasUserFloat("fitMass") ? cand.userFloat("fitMass") : -1) << endl;
1453       }
1454     }
1455     if (buDBasicSelect->accept(cand) && buDJPsiBasicSelect->accept(*jPsi) &&
1456         //         buDJPsiDaughterSelect->accept( *jPsi ) &&
1457         buDVertexSelect->accept(cand, BPHUserData::getByRef<reco::Vertex>(*jPsi, "primaryVertex")) &&
1458         (kptr->pt() > buDKPtMin)) {
1459       fillHisto("DDBu", cand, 'f');
1460       fillHisto("DDBuJPsi", *jPsi, 'c');
1461       if (flag_DoubleMu4_JpsiTrk_Displaced) {
1462         fillHisto("TDBu", cand, 'f');
1463         fillHisto("TDBuJPsi", *jPsi, 'c');
1464         *ofMap["DisplacedBu"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1465                               << (cand.hasUserFloat("fitMass") ? cand.userFloat("fitMass") : -1) << endl;
1466       }
1467     }
1468   }
1469 
1470   //////////// Bd -> JPsi Kx0 ////////////
1471 
1472   edm::Handle<vector<pat::CompositeCandidate>> bdCands;
1473   int ibd;
1474   int nbd = 0;
1475   if (useBd) {
1476     bdCandsToken.get(ev, bdCands);
1477     nbd = bdCands->size();
1478   }
1479 
1480   for (ibd = 0; ibd < nbd; ++ibd) {
1481     LogTrace("DataDump") << "*********** Bd " << ibd << "/" << nbd << " ***********";
1482     const pat::CompositeCandidate& cand = bdCands->at(ibd);
1483     const pat::CompositeCandidate* jPsi = BPHUserData::getByRef<pat::CompositeCandidate>(cand, "refToJPsi");
1484     LogTrace("DataDump") << "JPsi: " << jPsi;
1485     if (jPsi == nullptr)
1486       continue;
1487     if (!npJPsiBasicSelect->accept(*jPsi))
1488       continue;
1489     if (!npJPsiDaughterSelect->accept(*jPsi))
1490       continue;
1491     const pat::CompositeCandidate* kx0 = BPHUserData::getByRef<pat::CompositeCandidate>(cand, "refToKx0");
1492     LogTrace("DataDump") << "Kx0: " << kx0;
1493     if (kx0 == nullptr)
1494       continue;
1495     if (bdIBasicSelect->accept(cand) && bdIJPsiBasicSelect->accept(*jPsi) && bdIKx0BasicSelect->accept(*kx0) &&
1496         //         bdIJPsiDaughterSelect->accept( *jPsi ) &&
1497         bdIVertexSelect->accept(cand, BPHUserData::getByRef<reco::Vertex>(*jPsi, "primaryVertex"))) {
1498       fillHisto("DIBd", cand, 'f');
1499       fillHisto("DIBdJPsi", *jPsi, 'c');
1500       fillHisto("DIBdKx0", *kx0, 'c');
1501       if (flag_Dimuon25_Jpsi) {
1502         fillHisto("TIBd", cand, 'f');
1503         fillHisto("TIBdJPsi", *jPsi, 'c');
1504         fillHisto("TIBdKx0", *kx0, 'c');
1505         *ofMap["InclusiveBd"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1506                               << (cand.hasUserFloat("fitMass") ? cand.userFloat("fitMass") : -1) << endl;
1507       }
1508     }
1509     if (bdDBasicSelect->accept(cand) && bdDJPsiBasicSelect->accept(*jPsi) && bdDKx0BasicSelect->accept(*kx0) &&
1510         //         bdDJPsiDaughterSelect->accept( *jPsi ) &&
1511         bdDVertexSelect->accept(cand, BPHUserData::getByRef<reco::Vertex>(*jPsi, "primaryVertex"))) {
1512       fillHisto("DDBd", cand, 'f');
1513       fillHisto("DDBdJPsi", *jPsi, 'c');
1514       fillHisto("DDBdKx0", *kx0, 'c');
1515       if (flag_DoubleMu4_JpsiTrk_Displaced) {
1516         fillHisto("TDBd", cand, 'f');
1517         fillHisto("TDBdJPsi", *jPsi, 'c');
1518         fillHisto("TDBdKx0", *kx0, 'c');
1519         *ofMap["DisplacedBd"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1520                               << (cand.hasUserFloat("fitMass") ? cand.userFloat("fitMass") : -1) << endl;
1521       }
1522     }
1523   }
1524 
1525   //////////// Bs ////////////
1526 
1527   edm::Handle<vector<pat::CompositeCandidate>> bsCands;
1528   int ibs;
1529   int nbs = 0;
1530   if (useBs) {
1531     bsCandsToken.get(ev, bsCands);
1532     nbs = bsCands->size();
1533   }
1534 
1535   for (ibs = 0; ibs < nbs; ++ibs) {
1536     LogTrace("DataDump") << "*********** Bs " << ibs << "/" << nbs << " ***********";
1537     const pat::CompositeCandidate& cand = bsCands->at(ibs);
1538     const pat::CompositeCandidate* jPsi = BPHUserData::getByRef<pat::CompositeCandidate>(cand, "refToJPsi");
1539     LogTrace("DataDump") << "JPsi: " << jPsi;
1540     if (jPsi == nullptr)
1541       continue;
1542     if (!npJPsiBasicSelect->accept(*jPsi))
1543       continue;
1544     if (!npJPsiDaughterSelect->accept(*jPsi))
1545       continue;
1546     const pat::CompositeCandidate* phi = BPHUserData::getByRef<pat::CompositeCandidate>(cand, "refToPhi");
1547     LogTrace("DataDump") << "Phi: " << phi;
1548     if (phi == nullptr)
1549       continue;
1550     if (bsIBasicSelect->accept(cand) && bsIJPsiBasicSelect->accept(*jPsi) && bsIPhiBasicSelect->accept(*phi) &&
1551         //         bsIJPsiDaughterSelect->accept( *jPsi ) &&
1552         bsIVertexSelect->accept(cand, BPHUserData::getByRef<reco::Vertex>(*jPsi, "primaryVertex"))) {
1553       fillHisto("DIBs", cand, 'f');
1554       fillHisto("DIBsJPsi", *jPsi, 'c');
1555       fillHisto("DIBsPhi", *phi, 'c');
1556       if (flag_Dimuon25_Jpsi) {
1557         fillHisto("TIBs", cand, 'f');
1558         fillHisto("TIBsJPsi", *jPsi, 'c');
1559         fillHisto("TIBsPhi", *phi, 'c');
1560         *ofMap["InclusiveBs"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1561                               << (cand.hasUserFloat("fitMass") ? cand.userFloat("fitMass") : -1) << endl;
1562       }
1563     }
1564     if (bsDBasicSelect->accept(cand) && bsDJPsiBasicSelect->accept(*jPsi) && bsDPhiBasicSelect->accept(*phi) &&
1565         //         bsDJPsiDaughterSelect->accept( *jPsi ) &&
1566         bsDVertexSelect->accept(cand, BPHUserData::getByRef<reco::Vertex>(*jPsi, "primaryVertex"))) {
1567       fillHisto("DDBs", cand, 'f');
1568       fillHisto("DDBsJPsi", *jPsi, 'c');
1569       fillHisto("DDBsPhi", *phi, 'c');
1570       if (flag_DoubleMu4_JpsiTrk_Displaced) {
1571         fillHisto("TDBs", cand, 'f');
1572         fillHisto("TDBsJPsi", *jPsi, 'c');
1573         fillHisto("TDBsPhi", *phi, 'c');
1574         *ofMap["DisplacedBs"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1575                               << (cand.hasUserFloat("fitMass") ? cand.userFloat("fitMass") : -1) << endl;
1576       }
1577     }
1578   }
1579 
1580   //////////// K0s ////////////
1581 
1582   edm::Handle<vector<pat::CompositeCandidate>> k0Cands;
1583   int ik0;
1584   int nk0 = 0;
1585   if (useK0) {
1586     k0CandsToken.get(ev, k0Cands);
1587     nk0 = k0Cands->size();
1588   }
1589 
1590   for (ik0 = 0; ik0 < nk0; ++ik0) {
1591     LogTrace("DataDump") << "*********** K0 " << ik0 << "/" << nk0 << " ***********";
1592     const pat::CompositeCandidate& cand = k0Cands->at(ik0);
1593     fillHisto("DK0s", cand, 'f');
1594     *ofMap["K0s"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1595                   << (cand.hasUserFloat("fitMass") ? cand.userFloat("fitMass") : -1) << endl;
1596   }
1597 
1598   //////////// Lambda0 ////////////
1599 
1600   edm::Handle<vector<pat::CompositeCandidate>> l0Cands;
1601   int il0;
1602   int nl0 = 0;
1603   if (useL0) {
1604     l0CandsToken.get(ev, l0Cands);
1605     nl0 = l0Cands->size();
1606   }
1607 
1608   for (il0 = 0; il0 < nl0; ++il0) {
1609     LogTrace("DataDump") << "*********** Lambda0 " << il0 << "/" << nl0 << " ***********";
1610     const pat::CompositeCandidate& cand = l0Cands->at(il0);
1611     fillHisto("DLambda0", cand, 'f');
1612     *ofMap["Lambda0"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1613                       << (cand.hasUserFloat("fitMass") ? cand.userFloat("fitMass") : -1) << endl;
1614   }
1615 
1616   //////////// Bd -> JPsi K0s ////////////
1617 
1618   edm::Handle<vector<pat::CompositeCandidate>> b0Cands;
1619   int ib0;
1620   int nb0 = 0;
1621   if (useB0) {
1622     b0CandsToken.get(ev, b0Cands);
1623     nb0 = b0Cands->size();
1624   }
1625 
1626   //  cout << nb0 << ' ' << ev.id().run() << ' ' << ev.id().event();
1627   //  if ( nb0 ) cout << " *************************";
1628   //  cout << endl;
1629 
1630   for (ib0 = 0; ib0 < nb0; ++ib0) {
1631     LogTrace("DataDump") << "*********** B0 " << ib0 << "/" << nb0 << " ***********";
1632     const pat::CompositeCandidate& cand = b0Cands->at(ib0);
1633     const pat::CompositeCandidate* jPsi = BPHUserData::getByRef<pat::CompositeCandidate>(cand, "refToJPsi");
1634     LogTrace("DataDump") << "JPsi: " << jPsi;
1635     if (jPsi == nullptr)
1636       continue;
1637     if (!npJPsiBasicSelect->accept(*jPsi))
1638       continue;
1639     if (!npJPsiDaughterSelect->accept(*jPsi))
1640       continue;
1641     const pat::CompositeCandidate* k0s = BPHUserData::getByRef<pat::CompositeCandidate>(cand, "refToK0s");
1642     LogTrace("DataDump") << "K0s: " << k0s;
1643     if (k0s == nullptr)
1644       continue;
1645     if (b0IBasicSelect->accept(cand) && b0IJPsiBasicSelect->accept(*jPsi) && b0IK0sBasicSelect->accept(*k0s) &&
1646         //         b0IJPsiDaughterSelect->accept( *jPsi ) &&
1647         b0IVertexSelect->accept(cand, BPHUserData::getByRef<reco::Vertex>(*jPsi, "primaryVertex"))) {
1648       fillHisto("DIB0", cand, 'f');
1649       fillHisto("DIB0JPsi", *jPsi, 'c');
1650       fillHisto("DIB0K0s", *k0s, 'c');
1651       if (flag_Dimuon25_Jpsi) {
1652         fillHisto("TIB0", cand, 'f');
1653         fillHisto("TIB0JPsi", *jPsi, 'c');
1654         fillHisto("TIB0K0s", *k0s, 'c');
1655         *ofMap["InclusiveB0"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1656                               << (cand.hasUserFloat("fitMass") ? cand.userFloat("fitMass") : -1) << endl;
1657       }
1658     }
1659     if (b0DBasicSelect->accept(cand) && b0DJPsiBasicSelect->accept(*jPsi) && b0DK0sBasicSelect->accept(*k0s) &&
1660         //         b0DJPsiDaughterSelect->accept( *jPsi ) &&
1661         b0DVertexSelect->accept(cand, BPHUserData::getByRef<reco::Vertex>(*jPsi, "primaryVertex"))) {
1662       fillHisto("DDB0", cand, 'f');
1663       fillHisto("DDB0JPsi", *jPsi, 'c');
1664       fillHisto("DDB0K0s", *k0s, 'c');
1665       if (flag_DoubleMu4_JpsiTrk_Displaced) {
1666         fillHisto("TDB0", cand, 'f');
1667         fillHisto("TDB0JPsi", *jPsi, 'c');
1668         fillHisto("TDB0K0s", *k0s, 'c');
1669         *ofMap["DisplacedB0"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1670                               << (cand.hasUserFloat("fitMass") ? cand.userFloat("fitMass") : -1) << endl;
1671       }
1672     }
1673   }
1674 
1675   //////////// Lambdab -> JPsi Lambda0///////////
1676 
1677   edm::Handle<vector<pat::CompositeCandidate>> lbCands;
1678   int ilb;
1679   int nlb = 0;
1680   if (useLb) {
1681     lbCandsToken.get(ev, lbCands);
1682     nlb = lbCands->size();
1683   }
1684 
1685   for (ilb = 0; ilb < nlb; ++ilb) {
1686     LogTrace("DataDump") << "*********** Lambdab " << ilb << "/" << nlb << " ***********";
1687     const pat::CompositeCandidate& cand = lbCands->at(ilb);
1688     const pat::CompositeCandidate* jPsi = BPHUserData::getByRef<pat::CompositeCandidate>(cand, "refToJPsi");
1689     LogTrace("DataDump") << "JPsi: " << jPsi;
1690     if (jPsi == nullptr)
1691       continue;
1692     if (!npJPsiBasicSelect->accept(*jPsi))
1693       continue;
1694     if (!npJPsiDaughterSelect->accept(*jPsi))
1695       continue;
1696     const pat::CompositeCandidate* l0 = BPHUserData::getByRef<pat::CompositeCandidate>(cand, "refToLambda0");
1697     LogTrace("DataDump") << "Lambda0: " << l0;
1698     if (l0 == nullptr)
1699       continue;
1700     if (lbIBasicSelect->accept(cand) && lbIJPsiBasicSelect->accept(*jPsi) && lbILambda0BasicSelect->accept(*l0) &&
1701         //         lbIJPsiDaughterSelect->accept( *jPsi ) &&
1702         lbIVertexSelect->accept(cand, BPHUserData::getByRef<reco::Vertex>(*jPsi, "primaryVertex"))) {
1703       fillHisto("DILambdab", cand, 'f');
1704       fillHisto("DILbJPsi", *jPsi, 'c');
1705       fillHisto("DILbL0", *l0, 'c');
1706       if (flag_Dimuon25_Jpsi) {
1707         fillHisto("TILambdab", cand, 'f');
1708         fillHisto("TILbJPsi", *jPsi, 'c');
1709         fillHisto("TILbL0", *l0, 'c');
1710         *ofMap["InclusiveLambdab"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1711                                    << (cand.hasUserFloat("fitMass") ? cand.userFloat("fitMass") : -1) << endl;
1712       }
1713     }
1714     if (lbDBasicSelect->accept(cand) && lbDJPsiBasicSelect->accept(*jPsi) && lbDLambda0BasicSelect->accept(*l0) &&
1715         //         lbDJPsiDaughterSelect->accept( *jPsi ) &&
1716         lbDVertexSelect->accept(cand, BPHUserData::getByRef<reco::Vertex>(*jPsi, "primaryVertex"))) {
1717       fillHisto("DDLambdab", cand, 'f');
1718       fillHisto("DDLbJPsi", *jPsi, 'c');
1719       fillHisto("DDLbL0", *l0, 'c');
1720       if (flag_DoubleMu4_JpsiTrk_Displaced) {
1721         fillHisto("TDLambdab", cand, 'f');
1722         fillHisto("TDLbJPsi", *jPsi, 'c');
1723         fillHisto("TDLbL0", *l0, 'c');
1724         *ofMap["DisplacedLambdab"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1725                                    << (cand.hasUserFloat("fitMass") ? cand.userFloat("fitMass") : -1) << endl;
1726       }
1727     }
1728   }
1729 
1730   //////////// Bc ////////////
1731 
1732   edm::Handle<vector<pat::CompositeCandidate>> bcCands;
1733   int ibc;
1734   int nbc = 0;
1735   if (useBc) {
1736     bcCandsToken.get(ev, bcCands);
1737     nbc = bcCands->size();
1738   }
1739 
1740   for (ibc = 0; ibc < nbc; ++ibc) {
1741     LogTrace("DataDump") << "*********** Bc " << ibc << "/" << nbc << " ***********";
1742     const pat::CompositeCandidate& cand = bcCands->at(ibc);
1743     const pat::CompositeCandidate* jPsi = BPHUserData::getByRef<pat::CompositeCandidate>(cand, "refToJPsi");
1744     LogTrace("DataDump") << "JPsi: " << jPsi;
1745     if (jPsi == nullptr)
1746       continue;
1747     //    if ( BPHUserData::get( *jPsi, "dca", -1.0 ) < bcJPsiDcaMax ) continue;
1748     if (!npJPsiBasicSelect->accept(*jPsi))
1749       continue;
1750     if (!npJPsiDaughterSelect->accept(*jPsi))
1751       continue;
1752     const reco::Candidate* pptr = BPHDaughters::get(cand, 0.13, 0.14).front();
1753     if (pptr == nullptr)
1754       continue;
1755 
1756     if (bcIBasicSelect->accept(cand) && bcIJPsiBasicSelect->accept(*jPsi) &&
1757         //         bcIJPsiDaughterSelect->accept( *jPsi ) &&
1758         bcIJPsiVertexSelect->accept(*jPsi, BPHUserData::getByRef<reco::Vertex>(*jPsi, "primaryVertex")) &&
1759         bcIVertexSelect->accept(cand, BPHUserData::getByRef<reco::Vertex>(*jPsi, "primaryVertex")) &&
1760         (pptr->pt() > bcIPiPtMin)) {
1761       fillHisto("DIBc", cand, 'f');
1762       fillHisto("DIBcJPsi", *jPsi, 'c');
1763       if (flag_Dimuon25_Jpsi) {
1764         fillHisto("TIBc", cand, 'f');
1765         fillHisto("TIBcJPsi", *jPsi, 'c');
1766         *ofMap["InclusiveBc"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1767                               << (cand.hasUserFloat("fitMass") ? cand.userFloat("fitMass") : -1) << endl;
1768       }
1769     }
1770     if (bcDBasicSelect->accept(cand) && bcDJPsiBasicSelect->accept(*jPsi) &&
1771         //         bcDJPsiDaughterSelect->accept( *jPsi ) &&
1772         bcDVertexSelect->accept(cand, BPHUserData::getByRef<reco::Vertex>(*jPsi, "primaryVertex")) &&
1773         bcDVertexSelect->accept(cand, BPHUserData::getByRef<reco::Vertex>(*jPsi, "primaryVertex")) &&
1774         (pptr->pt() > bcDPiPtMin)) {
1775       fillHisto("DDBc", cand, 'f');
1776       fillHisto("DDBcJPsi", *jPsi, 'c');
1777       if (flag_DoubleMu4_JpsiTrk_Displaced) {
1778         fillHisto("TDBc", cand, 'f');
1779         fillHisto("TDBcJPsi", *jPsi, 'c');
1780         *ofMap["DisplacedBc"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1781                               << (cand.hasUserFloat("fitMass") ? cand.userFloat("fitMass") : -1) << endl;
1782       }
1783     }
1784   }
1785 
1786   //////////// X3872 ////////////
1787 
1788   edm::Handle<vector<pat::CompositeCandidate>> x3872Cands;
1789   int ix3872;
1790   int nx3872 = 0;
1791   if (useX3872) {
1792     x3872CandsToken.get(ev, x3872Cands);
1793     nx3872 = x3872Cands->size();
1794   }
1795 
1796   for (ix3872 = 0; ix3872 < nx3872; ++ix3872) {
1797     LogTrace("DataDump") << "*********** X3872 " << ix3872 << "/" << nx3872 << " ***********";
1798     const pat::CompositeCandidate& cand = x3872Cands->at(ix3872);
1799     const pat::CompositeCandidate* jPsi = BPHUserData::getByRef<pat::CompositeCandidate>(cand, "refToJPsi");
1800     LogTrace("DataDump") << "JPsi: " << jPsi;
1801     if (jPsi == nullptr)
1802       continue;
1803     //    if ( BPHUserData::get( *jPsi, "dca", -1.0 ) < x3872JPsiDcaMax ) continue;
1804     if (!npJPsiBasicSelect->accept(*jPsi))
1805       continue;
1806     if (!npJPsiDaughterSelect->accept(*jPsi))
1807       continue;
1808     const reco::Candidate* ppt1 = BPHDaughters::get(cand, 0.13, 0.14)[0];
1809     const reco::Candidate* ppt2 = BPHDaughters::get(cand, 0.13, 0.14)[1];
1810     if (ppt1 == nullptr)
1811       continue;
1812     if (ppt2 == nullptr)
1813       continue;
1814     if (x3872IBasicSelect->accept(cand) && x3872IJPsiBasicSelect->accept(*jPsi) &&
1815         //         x3872IJPsiDaughterSelect->accept( *jPsi ) &&
1816         x3872IJPsiVertexSelect->accept(*jPsi, BPHUserData::getByRef<reco::Vertex>(*jPsi, "primaryVertex")) &&
1817         x3872IVertexSelect->accept(cand, BPHUserData::getByRef<reco::Vertex>(*jPsi, "primaryVertex")) &&
1818         (ppt1->pt() > x3872IPiPtMin) && (ppt2->pt() > x3872IPiPtMin)) {
1819       fillHisto("DIX3872", cand, 'f');
1820       fillHisto("DIX3872JPsi", *jPsi, 'c');
1821       if (flag_Dimuon25_Jpsi) {
1822         fillHisto("TIX3872", cand, 'f');
1823         fillHisto("TIX3872JPsi", *jPsi, 'c');
1824         *ofMap["InclusiveX3872"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1825                                  << (cand.hasUserFloat("fitMass") ? cand.userFloat("fitMass") : -1) << endl;
1826       }
1827     }
1828     if (x3872DBasicSelect->accept(cand) && x3872DJPsiBasicSelect->accept(*jPsi) &&
1829         //         x3872DJPsiDaughterSelect->accept( *jPsi ) &&
1830         x3872DVertexSelect->accept(cand, BPHUserData::getByRef<reco::Vertex>(*jPsi, "primaryVertex")) &&
1831         x3872DVertexSelect->accept(cand, BPHUserData::getByRef<reco::Vertex>(*jPsi, "primaryVertex")) &&
1832         (ppt1->pt() > x3872DPiPtMin) && (ppt2->pt() > x3872DPiPtMin)) {
1833       fillHisto("DDX3872", cand, 'f');
1834       fillHisto("DDX3872JPsi", *jPsi, 'c');
1835       if (flag_DoubleMu4_JpsiTrk_Displaced) {
1836         fillHisto("TDX3872", cand, 'f');
1837         fillHisto("TDX3872JPsi", *jPsi, 'c');
1838         *ofMap["DisplacedX3872"] << ev.id().run() << ' ' << ev.id().luminosityBlock() << ' ' << ev.id().event() << ' '
1839                                  << (cand.hasUserFloat("fitMass") ? cand.userFloat("fitMass") : -1) << endl;
1840       }
1841     }
1842   }
1843 
1844   return;
1845 }
1846 
1847 void BPHHistoSpecificDecay::endJob() {
1848   //  tree->Write();
1849   return;
1850 }
1851 
1852 void BPHHistoSpecificDecay::fillHisto(const string& name, const pat::CompositeCandidate& cand, char svType) {
1853   *recoName = name;
1854   float mass = (cand.hasUserFloat("fitMass") ? cand.userFloat("fitMass") : -1);
1855   fillHisto("mass" + name, cand.mass());
1856   fillHisto("mfit" + name, mass);
1857   recoMass = mass;
1858   recoTime = -999999.0;
1859   recoErrT = -999999.0;
1860 
1861   const reco::Vertex* pvtx = BPHUserData::getByRef<reco::Vertex>(cand, "primaryVertex");
1862   if (pvtx == nullptr) {
1863     const pat::CompositeCandidate* jPsi = BPHUserData::getByRef<pat::CompositeCandidate>(cand, "refToJPsi");
1864     if (jPsi == nullptr)
1865       return;
1866     pvtx = BPHUserData::getByRef<reco::Vertex>(*jPsi, "primaryVertex");
1867   }
1868   //  if ( pvtx == nullptr ) return;
1869 
1870   if (pvtx != nullptr) {
1871     const reco::Vertex* svtx = nullptr;
1872     if (svType == 'f')
1873       svtx = BPHUserData::get<reco::Vertex>(cand, "fitVertex");
1874     if (svtx == nullptr)
1875       svtx = BPHUserData::get<reco::Vertex>(cand, "vertex");
1876     //  if ( svtx == nullptr ) return;
1877     if (svtx != nullptr) {
1878       float px;
1879       float py;
1880       const Vector3DBase<float, GlobalTag>* fmom =
1881           BPHUserData::get<Vector3DBase<float, GlobalTag>>(cand, "fitMomentum");
1882       if (fmom != nullptr) {
1883         px = fmom->x();
1884         py = fmom->y();
1885       } else {
1886         px = cand.px();
1887         py = cand.py();
1888       }
1889       double ctauPV2;
1890       double ctauErrPV2;
1891       VertexAnalysis::dist2D(pvtx, svtx, px, py, mass, ctauPV2, ctauErrPV2);
1892       recoTime = ctauPV2;
1893       recoErrT = ctauErrPV2;
1894       fillHisto("ctau" + name, ctauPV2);
1895     }
1896   }
1897   tree->Fill();
1898   return;
1899 }
1900 
1901 void BPHHistoSpecificDecay::fillHisto(const string& name, float x) {
1902   map<string, TH1F*>::iterator iter = histoMap.find(name);
1903   map<string, TH1F*>::iterator iend = histoMap.end();
1904   if (iter == iend)
1905     return;
1906   iter->second->Fill(x);
1907   return;
1908 }
1909 
1910 void BPHHistoSpecificDecay::createHisto(const string& name, int nbin, float hmin, float hmax) {
1911   histoMap[name] = fs->make<TH1F>(name.c_str(), name.c_str(), nbin, hmin, hmax);
1912   return;
1913 }
1914 
1915 #include "FWCore/Framework/interface/MakerMacros.h"
1916 
1917 DEFINE_FWK_MODULE(BPHHistoSpecificDecay);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!