Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-07-08 01:45:56 

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