Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 #include "HeavyFlavorAnalysis/SpecificDecay/plugins/BPHWriteSpecificDecay.h"
0002 
0003 #include "FWCore/Framework/interface/MakerMacros.h"
0004 
0005 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHMuonPtSelect.h"
0006 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHMuonEtaSelect.h"
0007 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHParticlePtSelect.h"
0008 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHParticleNeutralVeto.h"
0009 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHMassSelect.h"
0010 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHChi2Select.h"
0011 
0012 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHOniaToMuMuBuilder.h"
0013 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHKx0ToKPiBuilder.h"
0014 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHPhiToKKBuilder.h"
0015 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHBuToJPsiKBuilder.h"
0016 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHBuToPsi2SKBuilder.h"
0017 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHBsToJPsiPhiBuilder.h"
0018 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHBdToJPsiKxBuilder.h"
0019 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHK0sToPiPiBuilder.h"
0020 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHLambda0ToPPiBuilder.h"
0021 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHParticleMasses.h"
0022 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHBdToJPsiKsBuilder.h"
0023 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHLbToJPsiL0Builder.h"
0024 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHBcToJPsiPiBuilder.h"
0025 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHPsi2SToJPsiPiPiBuilder.h"
0026 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHX3872ToJPsiPiPiBuilder.h"
0027 
0028 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHAnalyzerTokenWrapper.h"
0029 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHRecoBuilder.h"
0030 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHRecoSelect.h"
0031 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHRecoCandidate.h"
0032 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHPlusMinusCandidate.h"
0033 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHMomentumSelect.h"
0034 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHVertexSelect.h"
0035 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHTrackReference.h"
0036 
0037 #include "DataFormats/PatCandidates/interface/Muon.h"
0038 #include "DataFormats/TrackReco/interface/Track.h"
0039 
0040 #include "DataFormats/PatCandidates/interface/GenericParticle.h"
0041 #include "DataFormats/PatCandidates/interface/CompositeCandidate.h"
0042 
0043 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0044 #include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistance.h"
0045 
0046 #include "FWCore/Framework/interface/ESHandle.h"
0047 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0048 
0049 #include <algorithm>
0050 #include <set>
0051 #include <string>
0052 #include <iostream>
0053 using namespace std;
0054 
0055 #define SET_PAR(TYPE, NAME, PSET) (NAME = PSET.getParameter<TYPE>(#NAME))
0056 // SET_PAR(string,xyz,ps);
0057 // is equivalent to
0058 // ( xyz = ps.getParameter< string >( "xyx" ) )
0059 
0060 BPHWriteSpecificDecay::BPHWriteSpecificDecay(const edm::ParameterSet& ps) {
0061   usePV = (!SET_PAR(string, pVertexLabel, ps).empty());
0062   usePM = (!SET_PAR(string, patMuonLabel, ps).empty());
0063   useCC = (!SET_PAR(string, ccCandsLabel, ps).empty());
0064   usePF = (!SET_PAR(string, pfCandsLabel, ps).empty());
0065   usePC = (!SET_PAR(string, pcCandsLabel, ps).empty());
0066   useGP = (!SET_PAR(string, gpCandsLabel, ps).empty());
0067   useK0 = (!SET_PAR(string, k0CandsLabel, ps).empty());
0068   useL0 = (!SET_PAR(string, l0CandsLabel, ps).empty());
0069   useKS = (!SET_PAR(string, kSCandsLabel, ps).empty());
0070   useLS = (!SET_PAR(string, lSCandsLabel, ps).empty());
0071   SET_PAR(string, oniaName, ps);
0072   SET_PAR(string, sdName, ps);
0073   SET_PAR(string, ssName, ps);
0074   SET_PAR(string, buName, ps);
0075   SET_PAR(string, bpName, ps);
0076   SET_PAR(string, bdName, ps);
0077   SET_PAR(string, bsName, ps);
0078   SET_PAR(string, k0Name, ps);
0079   SET_PAR(string, l0Name, ps);
0080   SET_PAR(string, b0Name, ps);
0081   SET_PAR(string, lbName, ps);
0082   SET_PAR(string, bcName, ps);
0083   SET_PAR(string, psi2SName, ps);
0084   SET_PAR(string, x3872Name, ps);
0085 
0086   SET_PAR(bool, writeMomentum, ps);
0087   SET_PAR(bool, writeVertex, ps);
0088 
0089   rMap["Onia"] = Onia;
0090   rMap["PhiMuMu"] = Pmm;
0091   rMap["Psi1"] = Psi1;
0092   rMap["Psi2"] = Psi2;
0093   rMap["Ups"] = Ups;
0094   rMap["Ups1"] = Ups1;
0095   rMap["Ups2"] = Ups2;
0096   rMap["Ups3"] = Ups3;
0097   rMap["Kx0"] = Kx0;
0098   rMap["PhiKK"] = Pkk;
0099   rMap["Bu"] = Bu;
0100   rMap["Bp"] = Bp;
0101   rMap["Bd"] = Bd;
0102   rMap["Bs"] = Bs;
0103   rMap["K0s"] = K0s;
0104   rMap["Lambda0"] = Lambda0;
0105   rMap["B0"] = B0;
0106   rMap["Lambdab"] = Lambdab;
0107   rMap["Bc"] = Bc;
0108   rMap["Psi2S"] = Psi2S;
0109   rMap["X3872"] = X3872;
0110 
0111   pMap["ptMin"] = ptMin;
0112   pMap["etaMax"] = etaMax;
0113   pMap["mJPsiMin"] = mPsiMin;
0114   pMap["mJPsiMax"] = mPsiMax;
0115   pMap["mKx0Min"] = mKx0Min;
0116   pMap["mKx0Max"] = mKx0Max;
0117   pMap["mPhiMin"] = mPhiMin;
0118   pMap["mPhiMax"] = mPhiMax;
0119   pMap["mK0sMin"] = mK0sMin;
0120   pMap["mK0sMax"] = mK0sMax;
0121   pMap["mLambda0Min"] = mLambda0Min;
0122   pMap["mLambda0Max"] = mLambda0Max;
0123   pMap["massMin"] = massMin;
0124   pMap["massMax"] = massMax;
0125   pMap["probMin"] = probMin;
0126   pMap["massFitMin"] = mFitMin;
0127   pMap["massFitMax"] = mFitMax;
0128   pMap["constrMass"] = constrMass;
0129   pMap["constrSigma"] = constrSigma;
0130 
0131   fMap["requireJPsi"] = requireJPsi;
0132   fMap["constrMJPsi"] = constrMJPsi;
0133   fMap["constrMPsi2"] = constrMPsi2;
0134 
0135   fMap["writeCandidate"] = writeCandidate;
0136 
0137   recoOnia = recoKx0 = writeKx0 = recoPkk = writePkk = recoBu = writeBu = recoBp = writeBp = recoBd = writeBd = recoBs =
0138       writeBs = recoK0s = writeK0s = recoLambda0 = writeLambda0 = recoB0 = writeB0 = recoLambdab = writeLambdab =
0139           recoBc = writeBc = recoPsi2S = writePsi2S = recoX3872 = writeX3872 = false;
0140 
0141   writeOnia = true;
0142   const vector<edm::ParameterSet> recoSelect = ps.getParameter<vector<edm::ParameterSet> >("recoSelect");
0143   int iSel;
0144   int nSel = recoSelect.size();
0145   for (iSel = 0; iSel < nSel; ++iSel)
0146     setRecoParameters(recoSelect[iSel]);
0147   if (!recoOnia)
0148     writeOnia = false;
0149 
0150   if (recoBu)
0151     recoOnia = true;
0152   if (recoBd)
0153     recoOnia = recoKx0 = true;
0154   if (recoBs)
0155     recoOnia = recoPkk = true;
0156   if (recoB0)
0157     recoOnia = recoK0s = true;
0158   if (recoLambdab)
0159     recoOnia = recoLambda0 = true;
0160   if (recoBc)
0161     recoOnia = true;
0162   if (recoPsi2S)
0163     recoOnia = true;
0164   if (recoX3872)
0165     recoOnia = true;
0166   if (writeBu)
0167     writeOnia = true;
0168   if (writeBd)
0169     writeOnia = writeKx0 = true;
0170   if (writeBs)
0171     writeOnia = writePkk = true;
0172   if (writeB0)
0173     writeOnia = writeK0s = true;
0174   if (writeLambdab)
0175     writeOnia = writeLambda0 = true;
0176   if (writeBc)
0177     writeOnia = true;
0178   if (writePsi2S)
0179     writeOnia = true;
0180   if (writeX3872)
0181     writeOnia = true;
0182   if (recoBp && !recoPsi2S && !recoX3872)
0183     recoPsi2S = true;
0184   if (writeBp && !writePsi2S && !writeX3872)
0185     writePsi2S = true;
0186   allKx0 = (parMap[Kx0][requireJPsi] < 0);
0187   allPkk = (parMap[Pkk][requireJPsi] < 0);
0188   allK0s = (parMap[K0s][requireJPsi] < 0);
0189   allLambda0 = (parMap[Lambda0][requireJPsi] < 0);
0190 
0191   esConsume<MagneticField, IdealMagneticFieldRecord>(magFieldToken);
0192   esConsume<TransientTrackBuilder, TransientTrackRecord>(ttBToken, "TransientTrackBuilder");
0193   if (usePV)
0194     consume<vector<reco::Vertex> >(pVertexToken, pVertexLabel);
0195   if (usePM)
0196     consume<pat::MuonCollection>(patMuonToken, patMuonLabel);
0197   if (useCC)
0198     consume<vector<pat::CompositeCandidate> >(ccCandsToken, ccCandsLabel);
0199   if (usePF)
0200     consume<vector<reco::PFCandidate> >(pfCandsToken, pfCandsLabel);
0201   if (usePC)
0202     consume<vector<BPHTrackReference::candidate> >(pcCandsToken, pcCandsLabel);
0203   if (useGP)
0204     consume<vector<pat::GenericParticle> >(gpCandsToken, gpCandsLabel);
0205   if (useK0)
0206     consume<vector<reco::VertexCompositeCandidate> >(k0CandsToken, k0CandsLabel);
0207   if (useL0)
0208     consume<vector<reco::VertexCompositeCandidate> >(l0CandsToken, l0CandsLabel);
0209   if (useKS)
0210     consume<vector<reco::VertexCompositePtrCandidate> >(kSCandsToken, kSCandsLabel);
0211   if (useLS)
0212     consume<vector<reco::VertexCompositePtrCandidate> >(lSCandsToken, lSCandsLabel);
0213 
0214   if (writeOnia)
0215     produces<pat::CompositeCandidateCollection>(oniaName);
0216   if (writeKx0)
0217     produces<pat::CompositeCandidateCollection>(sdName);
0218   if (writePkk)
0219     produces<pat::CompositeCandidateCollection>(ssName);
0220   if (writeBu)
0221     produces<pat::CompositeCandidateCollection>(buName);
0222   if (writeBp)
0223     produces<pat::CompositeCandidateCollection>(bpName);
0224   if (writeBd)
0225     produces<pat::CompositeCandidateCollection>(bdName);
0226   if (writeBs)
0227     produces<pat::CompositeCandidateCollection>(bsName);
0228   if (writeK0s)
0229     produces<pat::CompositeCandidateCollection>(k0Name);
0230   if (writeLambda0)
0231     produces<pat::CompositeCandidateCollection>(l0Name);
0232   if (writeB0)
0233     produces<pat::CompositeCandidateCollection>(b0Name);
0234   if (writeLambdab)
0235     produces<pat::CompositeCandidateCollection>(lbName);
0236   if (writeBc)
0237     produces<pat::CompositeCandidateCollection>(bcName);
0238   if (writePsi2S)
0239     produces<pat::CompositeCandidateCollection>(psi2SName);
0240   if (writeX3872)
0241     produces<pat::CompositeCandidateCollection>(x3872Name);
0242 }
0243 
0244 void BPHWriteSpecificDecay::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0245   edm::ParameterSetDescription desc;
0246   desc.add<string>("pVertexLabel", "");
0247   desc.add<string>("patMuonLabel", "");
0248   desc.add<string>("ccCandsLabel", "");
0249   desc.add<string>("pfCandsLabel", "");
0250   desc.add<string>("pcCandsLabel", "");
0251   desc.add<string>("gpCandsLabel", "");
0252   desc.add<string>("k0CandsLabel", "");
0253   desc.add<string>("l0CandsLabel", "");
0254   desc.add<string>("kSCandsLabel", "");
0255   desc.add<string>("lSCandsLabel", "");
0256   desc.add<string>("oniaName", "oniaCand");
0257   desc.add<string>("sdName", "kx0Cand");
0258   desc.add<string>("ssName", "phiCand");
0259   desc.add<string>("buName", "buFitted");
0260   desc.add<string>("bpName", "bpFitted");
0261   desc.add<string>("bdName", "bdFitted");
0262   desc.add<string>("bsName", "bsFitted");
0263   desc.add<string>("k0Name", "k0Fitted");
0264   desc.add<string>("l0Name", "l0Fitted");
0265   desc.add<string>("b0Name", "b0Fitted");
0266   desc.add<string>("lbName", "lbFitted");
0267   desc.add<string>("bcName", "bcFitted");
0268   desc.add<string>("psi2SName", "psi2SFitted");
0269   desc.add<string>("x3872Name", "x3872Fitted");
0270   desc.add<bool>("writeVertex", true);
0271   desc.add<bool>("writeMomentum", true);
0272   edm::ParameterSetDescription dpar;
0273   dpar.add<string>("name");
0274   dpar.add<double>("ptMin", -2.0e35);
0275   dpar.add<double>("etaMax", -2.0e35);
0276   dpar.add<double>("mJPsiMin", -2.0e35);
0277   dpar.add<double>("mJPsiMax", -2.0e35);
0278   dpar.add<double>("mKx0Min", -2.0e35);
0279   dpar.add<double>("mKx0Max", -2.0e35);
0280   dpar.add<double>("mPhiMin", -2.0e35);
0281   dpar.add<double>("mPhiMax", -2.0e35);
0282   dpar.add<double>("mK0sMin", -2.0e35);
0283   dpar.add<double>("mK0sMax", -2.0e35);
0284   dpar.add<double>("mLambda0Min", -2.0e35);
0285   dpar.add<double>("mLambda0Max", -2.0e35);
0286   dpar.add<double>("massMin", -2.0e35);
0287   dpar.add<double>("massMax", -2.0e35);
0288   dpar.add<double>("probMin", -2.0e35);
0289   dpar.add<double>("massFitMin", -2.0e35);
0290   dpar.add<double>("massFitMax", -2.0e35);
0291   dpar.add<double>("constrMass", -2.0e35);
0292   dpar.add<double>("constrSigma", -2.0e35);
0293   dpar.add<bool>("requireJPsi", true);
0294   dpar.add<bool>("constrMJPsi", true);
0295   dpar.add<bool>("constrMPsi2", true);
0296   dpar.add<bool>("writeCandidate", true);
0297   vector<edm::ParameterSet> rpar;
0298   desc.addVPSet("recoSelect", dpar, rpar);
0299   descriptions.add("bphWriteSpecificDecay", desc);
0300   return;
0301 }
0302 
0303 void BPHWriteSpecificDecay::produce(edm::Event& ev, const edm::EventSetup& es) {
0304   BPHEventSetupWrapper ew(es, BPHRecoCandidate::transientTrackBuilder, &ttBToken);
0305   fill(ev, ew);
0306   if (writeOnia)
0307     write(ev, lFull, oniaName);
0308   if (writeKx0)
0309     write(ev, lSd, sdName);
0310   if (writePkk)
0311     write(ev, lSs, ssName);
0312   if (writeBu)
0313     write(ev, lBu, buName);
0314   if (writeBp)
0315     write(ev, lBp, bpName);
0316   if (writeBd)
0317     write(ev, lBd, bdName);
0318   if (writeBs)
0319     write(ev, lBs, bsName);
0320   if (writeK0s)
0321     write(ev, lK0, k0Name);
0322   if (writeLambda0)
0323     write(ev, lL0, l0Name);
0324   if (writeB0)
0325     write(ev, lB0, b0Name);
0326   if (writeLambdab)
0327     write(ev, lLb, lbName);
0328   if (writeBc)
0329     write(ev, lBc, bcName);
0330   if (writePsi2S)
0331     write(ev, lPsi2S, psi2SName);
0332   if (writeX3872)
0333     write(ev, lX3872, x3872Name);
0334   return;
0335 }
0336 
0337 void BPHWriteSpecificDecay::fill(edm::Event& ev, const BPHEventSetupWrapper& es) {
0338   lFull.clear();
0339   lJPsi.clear();
0340   lSd.clear();
0341   lSs.clear();
0342   lBu.clear();
0343   lBp.clear();
0344   lBd.clear();
0345   lBs.clear();
0346   lK0.clear();
0347   lL0.clear();
0348   lB0.clear();
0349   lLb.clear();
0350   lBc.clear();
0351   lPsi2S.clear();
0352   lX3872.clear();
0353   jPsiOMap.clear();
0354   daughMap.clear();
0355   pvRefMap.clear();
0356   ccRefMap.clear();
0357 
0358   // get magnetic field
0359   // data are got through "BPHESTokenWrapper" interface to allow
0360   // uniform access in different CMSSW versions
0361   edm::ESHandle<MagneticField> magneticField;
0362   magFieldToken.get(*es.get(), magneticField);
0363 
0364   // get object collections
0365   // collections are got through "BPHTokenWrapper" interface to allow
0366   // uniform access in different CMSSW versions
0367 
0368   edm::Handle<std::vector<reco::Vertex> > pVertices;
0369   pVertexToken.get(ev, pVertices);
0370   int npv = pVertices->size();
0371 
0372   int nrc = 0;
0373 
0374   // get reco::PFCandidate collection (in full AOD )
0375   edm::Handle<vector<reco::PFCandidate> > pfCands;
0376   if (usePF) {
0377     pfCandsToken.get(ev, pfCands);
0378     nrc = pfCands->size();
0379   }
0380 
0381   // get pat::PackedCandidate collection (in MiniAOD)
0382   // pat::PackedCandidate is not defined in CMSSW_5XY, so a
0383   // typedef (BPHTrackReference::candidate) is used, actually referring
0384   // to pat::PackedCandidate only for CMSSW versions where it's defined
0385   edm::Handle<vector<BPHTrackReference::candidate> > pcCands;
0386   if (usePC) {
0387     pcCandsToken.get(ev, pcCands);
0388     nrc = pcCands->size();
0389   }
0390 
0391   // get pat::GenericParticle collection (in skimmed data)
0392   edm::Handle<vector<pat::GenericParticle> > gpCands;
0393   if (useGP) {
0394     gpCandsToken.get(ev, gpCands);
0395     nrc = gpCands->size();
0396   }
0397 
0398   // get pat::Muon collection (in full AOD and MiniAOD)
0399   edm::Handle<pat::MuonCollection> patMuon;
0400   if (usePM) {
0401     patMuonToken.get(ev, patMuon);
0402   }
0403 
0404   // get K0 reco::VertexCompositeCandidate collection (in full AOD)
0405   edm::Handle<std::vector<reco::VertexCompositeCandidate> > k0Cand;
0406   if (useK0) {
0407     k0CandsToken.get(ev, k0Cand);
0408   }
0409 
0410   // get Lambda0 reco::VertexCompositeCandidate collection (in full AOD)
0411   edm::Handle<std::vector<reco::VertexCompositeCandidate> > l0Cand;
0412   if (useL0) {
0413     l0CandsToken.get(ev, l0Cand);
0414   }
0415 
0416   // get K0 reco::VertexCompositePtrCandidate collection (in MiniAOD)
0417   edm::Handle<std::vector<reco::VertexCompositePtrCandidate> > kSCand;
0418   if (useKS) {
0419     kSCandsToken.get(ev, kSCand);
0420   }
0421 
0422   // get Lambda0 reco::VertexCompositePtrCandidate collection (in MiniAOD)
0423   edm::Handle<std::vector<reco::VertexCompositePtrCandidate> > lSCand;
0424   if (useLS) {
0425     lSCandsToken.get(ev, lSCand);
0426   }
0427 
0428   // get muons from pat::CompositeCandidate objects describing onia;
0429   // muons from all composite objects are copied to an unique std::vector
0430   vector<const reco::Candidate*> muDaugs;
0431   set<const pat::Muon*> muonSet;
0432   typedef multimap<const reco::Candidate*, const pat::CompositeCandidate*> mu_cc_map;
0433   mu_cc_map muCCMap;
0434   if (useCC) {
0435     edm::Handle<vector<pat::CompositeCandidate> > ccCands;
0436     ccCandsToken.get(ev, ccCands);
0437     int n = ccCands->size();
0438     muDaugs.clear();
0439     muDaugs.reserve(n);
0440     muonSet.clear();
0441     set<const pat::Muon*>::const_iterator iter;
0442     set<const pat::Muon*>::const_iterator iend;
0443     int i;
0444     for (i = 0; i < n; ++i) {
0445       const pat::CompositeCandidate& cc = ccCands->at(i);
0446       int j;
0447       int m = cc.numberOfDaughters();
0448       for (j = 0; j < m; ++j) {
0449         const reco::Candidate* dp = cc.daughter(j);
0450         const pat::Muon* mp = dynamic_cast<const pat::Muon*>(dp);
0451         iter = muonSet.begin();
0452         iend = muonSet.end();
0453         bool add = (mp != nullptr) && (muonSet.find(mp) == iend);
0454         while (add && (iter != iend)) {
0455           if (BPHRecoBuilder::sameTrack(mp, *iter++, 1.0e-5))
0456             add = false;
0457         }
0458         if (add)
0459           muonSet.insert(mp);
0460         // associate muon to the CompositeCandidate containing it
0461         muCCMap.insert(pair<const reco::Candidate*, const pat::CompositeCandidate*>(dp, &cc));
0462       }
0463     }
0464     iter = muonSet.begin();
0465     iend = muonSet.end();
0466     while (iter != iend)
0467       muDaugs.push_back(*iter++);
0468   }
0469 
0470   map<recoType, map<parType, double> >::const_iterator rIter = parMap.begin();
0471   map<recoType, map<parType, double> >::const_iterator rIend = parMap.end();
0472 
0473   // reconstruct quarkonia
0474 
0475   BPHOniaToMuMuBuilder* onia = nullptr;
0476   if (recoOnia) {
0477     if (usePM)
0478       onia = new BPHOniaToMuMuBuilder(
0479           es, BPHRecoBuilder::createCollection(patMuon, "ingmcf"), BPHRecoBuilder::createCollection(patMuon, "ingmcf"));
0480     else if (useCC)
0481       onia = new BPHOniaToMuMuBuilder(
0482           es, BPHRecoBuilder::createCollection(muDaugs, "ingmcf"), BPHRecoBuilder::createCollection(muDaugs, "ingmcf"));
0483   }
0484 
0485   if (onia != nullptr) {
0486     while (rIter != rIend) {
0487       const map<recoType, map<parType, double> >::value_type& rEntry = *rIter++;
0488       recoType rType = rEntry.first;
0489       const map<parType, double>& pMap = rEntry.second;
0490       BPHOniaToMuMuBuilder::oniaType type;
0491       switch (rType) {
0492         case Pmm:
0493           type = BPHOniaToMuMuBuilder::Phi;
0494           break;
0495         case Psi1:
0496           type = BPHOniaToMuMuBuilder::Psi1;
0497           break;
0498         case Psi2:
0499           type = BPHOniaToMuMuBuilder::Psi2;
0500           break;
0501         case Ups:
0502           type = BPHOniaToMuMuBuilder::Ups;
0503           break;
0504         case Ups1:
0505           type = BPHOniaToMuMuBuilder::Ups1;
0506           break;
0507         case Ups2:
0508           type = BPHOniaToMuMuBuilder::Ups2;
0509           break;
0510         case Ups3:
0511           type = BPHOniaToMuMuBuilder::Ups3;
0512           break;
0513         default:
0514           continue;
0515       }
0516       map<parType, double>::const_iterator pIter = pMap.begin();
0517       map<parType, double>::const_iterator pIend = pMap.end();
0518       while (pIter != pIend) {
0519         const map<parType, double>::value_type& pEntry = *pIter++;
0520         parType id = pEntry.first;
0521         double pv = pEntry.second;
0522         switch (id) {
0523           case ptMin:
0524             onia->setPtMin(type, pv);
0525             break;
0526           case etaMax:
0527             onia->setEtaMax(type, pv);
0528             break;
0529           case massMin:
0530             onia->setMassMin(type, pv);
0531             break;
0532           case massMax:
0533             onia->setMassMax(type, pv);
0534             break;
0535           case probMin:
0536             onia->setProbMin(type, pv);
0537             break;
0538           case constrMass:
0539             onia->setConstr(type, pv, onia->getConstrSigma(type));
0540             break;
0541           case constrSigma:
0542             onia->setConstr(type, onia->getConstrMass(type), pv);
0543             break;
0544           default:
0545             break;
0546         }
0547       }
0548     }
0549     lFull = onia->build();
0550   }
0551 
0552   // associate onia to primary vertex
0553 
0554   int iFull;
0555   int nFull = lFull.size();
0556   map<const BPHRecoCandidate*, const reco::Vertex*> oniaVtxMap;
0557 
0558   typedef mu_cc_map::const_iterator mu_cc_iter;
0559   for (iFull = 0; iFull < nFull; ++iFull) {
0560     const reco::Vertex* pVtx = nullptr;
0561     int pvId = 0;
0562     const BPHPlusMinusCandidate* ptr = lFull[iFull].get();
0563     const std::vector<const reco::Candidate*>& daugs = ptr->daughters();
0564 
0565     // try to recover primary vertex association in skim data:
0566     // get the CompositeCandidate containing both muons
0567     pair<mu_cc_iter, mu_cc_iter> cc0 = muCCMap.equal_range(ptr->originalReco(daugs[0]));
0568     pair<mu_cc_iter, mu_cc_iter> cc1 = muCCMap.equal_range(ptr->originalReco(daugs[1]));
0569     mu_cc_iter iter0 = cc0.first;
0570     mu_cc_iter iend0 = cc0.second;
0571     mu_cc_iter iter1 = cc1.first;
0572     mu_cc_iter iend1 = cc1.second;
0573     while ((iter0 != iend0) && (pVtx == nullptr)) {
0574       const pat::CompositeCandidate* ccp = iter0++->second;
0575       while (iter1 != iend1) {
0576         if (ccp != iter1++->second)
0577           continue;
0578         pVtx = ccp->userData<reco::Vertex>("PVwithmuons");
0579         const reco::Vertex* sVtx = nullptr;
0580         const reco::Vertex::Point& pPos = pVtx->position();
0581         float dMin = 999999.;
0582         int ipv;
0583         for (ipv = 0; ipv < npv; ++ipv) {
0584           const reco::Vertex* tVtx = &pVertices->at(ipv);
0585           const reco::Vertex::Point& tPos = tVtx->position();
0586           float dist = pow(pPos.x() - tPos.x(), 2) + pow(pPos.y() - tPos.y(), 2) + pow(pPos.z() - tPos.z(), 2);
0587           if (dist < dMin) {
0588             dMin = dist;
0589             sVtx = tVtx;
0590             pvId = ipv;
0591           }
0592         }
0593         pVtx = sVtx;
0594         break;
0595       }
0596     }
0597 
0598     // if not found, as for other type of input data,
0599     // try to get the nearest primary vertex in z direction
0600     if (pVtx == nullptr) {
0601       const reco::Vertex::Point& sVtp = ptr->vertex().position();
0602       GlobalPoint cPos(sVtp.x(), sVtp.y(), sVtp.z());
0603       const pat::CompositeCandidate& sCC = ptr->composite();
0604       GlobalVector cDir(sCC.px(), sCC.py(), sCC.pz());
0605       GlobalPoint bPos(0.0, 0.0, 0.0);
0606       GlobalVector bDir(0.0, 0.0, 1.0);
0607       TwoTrackMinimumDistance ttmd;
0608       bool state = ttmd.calculate(GlobalTrajectoryParameters(cPos, cDir, TrackCharge(0), &(*magneticField)),
0609                                   GlobalTrajectoryParameters(bPos, bDir, TrackCharge(0), &(*magneticField)));
0610       float minDz = 999999.;
0611       float extrapZ = (state ? ttmd.points().first.z() : -9e20);
0612       int ipv;
0613       for (ipv = 0; ipv < npv; ++ipv) {
0614         const reco::Vertex& tVtx = pVertices->at(ipv);
0615         float deltaZ = fabs(extrapZ - tVtx.position().z());
0616         if (deltaZ < minDz) {
0617           minDz = deltaZ;
0618           pVtx = &tVtx;
0619           pvId = ipv;
0620         }
0621       }
0622     }
0623 
0624     oniaVtxMap[ptr] = pVtx;
0625     pvRefMap[ptr] = vertex_ref(pVertices, pvId);
0626   }
0627   pVertexToken.get(ev, pVertices);
0628 
0629   // get JPsi subsample and associate JPsi candidate to original
0630   // generic onia candidate
0631   if (nFull)
0632     lJPsi = onia->getList(BPHOniaToMuMuBuilder::Psi1);
0633 
0634   bool jPsiFound = !lJPsi.empty();
0635   delete onia;
0636 
0637   if (!nrc)
0638     return;
0639 
0640   int ij;
0641   int io;
0642   int nj = lJPsi.size();
0643   int no = lFull.size();
0644   for (ij = 0; ij < nj; ++ij) {
0645     const BPHRecoCandidate* jp = lJPsi[ij].get();
0646     for (io = 0; io < no; ++io) {
0647       const BPHRecoCandidate* oc = lFull[io].get();
0648       if ((jp->originalReco(jp->getDaug("MuPos")) == oc->originalReco(oc->getDaug("MuPos"))) &&
0649           (jp->originalReco(jp->getDaug("MuNeg")) == oc->originalReco(oc->getDaug("MuNeg")))) {
0650         jPsiOMap[jp] = oc;
0651         break;
0652       }
0653     }
0654   }
0655 
0656   // build and dump Bu
0657 
0658   BPHBuToJPsiKBuilder* bu = nullptr;
0659   if (recoBu && jPsiFound) {
0660     if (usePF)
0661       bu = new BPHBuToJPsiKBuilder(es, lJPsi, BPHRecoBuilder::createCollection(pfCands, "f"));
0662     else if (usePC)
0663       bu = new BPHBuToJPsiKBuilder(es, lJPsi, BPHRecoBuilder::createCollection(pcCands, "p"));
0664     else if (useGP)
0665       bu = new BPHBuToJPsiKBuilder(es, lJPsi, BPHRecoBuilder::createCollection(gpCands, "h"));
0666   }
0667 
0668   if (bu != nullptr) {
0669     rIter = parMap.find(Bu);
0670     if (rIter != rIend) {
0671       const map<parType, double>& pMap = rIter->second;
0672       map<parType, double>::const_iterator pIter = pMap.begin();
0673       map<parType, double>::const_iterator pIend = pMap.end();
0674       while (pIter != pIend) {
0675         const map<parType, double>::value_type& pEntry = *pIter++;
0676         parType id = pEntry.first;
0677         double pv = pEntry.second;
0678         switch (id) {
0679           case ptMin:
0680             bu->setKPtMin(pv);
0681             break;
0682           case etaMax:
0683             bu->setKEtaMax(pv);
0684             break;
0685           case mPsiMin:
0686             bu->setJPsiMassMin(pv);
0687             break;
0688           case mPsiMax:
0689             bu->setJPsiMassMax(pv);
0690             break;
0691           case massMin:
0692             bu->setMassMin(pv);
0693             break;
0694           case massMax:
0695             bu->setMassMax(pv);
0696             break;
0697           case probMin:
0698             bu->setProbMin(pv);
0699             break;
0700           case mFitMin:
0701             bu->setMassFitMin(pv);
0702             break;
0703           case mFitMax:
0704             bu->setMassFitMax(pv);
0705             break;
0706           case constrMJPsi:
0707             bu->setConstr(pv > 0);
0708             break;
0709           case writeCandidate:
0710             writeBu = (pv > 0);
0711             break;
0712           default:
0713             break;
0714         }
0715       }
0716     }
0717     lBu = bu->build();
0718     delete bu;
0719   }
0720 
0721   // build and dump Kx0
0722 
0723   vector<BPHPlusMinusConstCandPtr> lKx0;
0724   BPHKx0ToKPiBuilder* kx0 = nullptr;
0725   if (recoKx0 && (jPsiFound || allKx0)) {
0726     if (usePF)
0727       kx0 = new BPHKx0ToKPiBuilder(
0728           es, BPHRecoBuilder::createCollection(pfCands, "f"), BPHRecoBuilder::createCollection(pfCands, "f"));
0729     else if (usePC)
0730       kx0 = new BPHKx0ToKPiBuilder(
0731           es, BPHRecoBuilder::createCollection(pcCands, "p"), BPHRecoBuilder::createCollection(pcCands, "p"));
0732     else if (useGP)
0733       kx0 = new BPHKx0ToKPiBuilder(
0734           es, BPHRecoBuilder::createCollection(gpCands, "h"), BPHRecoBuilder::createCollection(gpCands, "h"));
0735   }
0736 
0737   set<BPHRecoConstCandPtr> sKx0;
0738 
0739   if (kx0 != nullptr) {
0740     rIter = parMap.find(Kx0);
0741     if (rIter != rIend) {
0742       const map<parType, double>& pMap = rIter->second;
0743       map<parType, double>::const_iterator pIter = pMap.begin();
0744       map<parType, double>::const_iterator pIend = pMap.end();
0745       while (pIter != pIend) {
0746         const map<parType, double>::value_type& pEntry = *pIter++;
0747         parType id = pEntry.first;
0748         double pv = pEntry.second;
0749         switch (id) {
0750           case ptMin:
0751             kx0->setPtMin(pv);
0752             break;
0753           case etaMax:
0754             kx0->setEtaMax(pv);
0755             break;
0756           case massMin:
0757             kx0->setMassMin(pv);
0758             break;
0759           case massMax:
0760             kx0->setMassMax(pv);
0761             break;
0762           case probMin:
0763             kx0->setProbMin(pv);
0764             break;
0765           case writeCandidate:
0766             writeKx0 = (pv > 0);
0767             break;
0768           default:
0769             break;
0770         }
0771       }
0772     }
0773     lKx0 = kx0->build();
0774     if (allKx0)
0775       sKx0.insert(lKx0.begin(), lKx0.end());
0776     delete kx0;
0777   }
0778 
0779   bool kx0Found = !lKx0.empty();
0780 
0781   // build and dump Bd -> JPsi Kx0
0782 
0783   if (recoBd && jPsiFound && kx0Found) {
0784     BPHBdToJPsiKxBuilder* bd = new BPHBdToJPsiKxBuilder(es, lJPsi, lKx0);
0785     rIter = parMap.find(Bd);
0786     if (rIter != rIend) {
0787       const map<parType, double>& pMap = rIter->second;
0788       map<parType, double>::const_iterator pIter = pMap.begin();
0789       map<parType, double>::const_iterator pIend = pMap.end();
0790       while (pIter != pIend) {
0791         const map<parType, double>::value_type& pEntry = *pIter++;
0792         parType id = pEntry.first;
0793         double pv = pEntry.second;
0794         switch (id) {
0795           case mPsiMin:
0796             bd->setJPsiMassMin(pv);
0797             break;
0798           case mPsiMax:
0799             bd->setJPsiMassMax(pv);
0800             break;
0801           case mKx0Min:
0802             bd->setKxMassMin(pv);
0803             break;
0804           case mKx0Max:
0805             bd->setKxMassMax(pv);
0806             break;
0807           case massMin:
0808             bd->setMassMin(pv);
0809             break;
0810           case massMax:
0811             bd->setMassMax(pv);
0812             break;
0813           case probMin:
0814             bd->setProbMin(pv);
0815             break;
0816           case mFitMin:
0817             bd->setMassFitMin(pv);
0818             break;
0819           case mFitMax:
0820             bd->setMassFitMax(pv);
0821             break;
0822           case constrMJPsi:
0823             bd->setConstr(pv > 0);
0824             break;
0825           case writeCandidate:
0826             writeBd = (pv > 0);
0827             break;
0828           default:
0829             break;
0830         }
0831       }
0832     }
0833 
0834     lBd = bd->build();
0835     delete bd;
0836 
0837     int iBd;
0838     int nBd = lBd.size();
0839     for (iBd = 0; iBd < nBd; ++iBd)
0840       sKx0.insert(lBd[iBd]->getComp("Kx0"));
0841   }
0842   set<BPHRecoConstCandPtr>::const_iterator kx0_iter = sKx0.begin();
0843   set<BPHRecoConstCandPtr>::const_iterator kx0_iend = sKx0.end();
0844   lSd.reserve(sKx0.size());
0845   while (kx0_iter != kx0_iend)
0846     lSd.push_back(*kx0_iter++);
0847 
0848   // build and dump Phi
0849 
0850   vector<BPHPlusMinusConstCandPtr> lPhi;
0851   BPHPhiToKKBuilder* phi = nullptr;
0852   if (recoPkk && (jPsiFound || allPkk)) {
0853     if (usePF)
0854       phi = new BPHPhiToKKBuilder(
0855           es, BPHRecoBuilder::createCollection(pfCands, "f"), BPHRecoBuilder::createCollection(pfCands, "f"));
0856     else if (usePC)
0857       phi = new BPHPhiToKKBuilder(
0858           es, BPHRecoBuilder::createCollection(pcCands, "p"), BPHRecoBuilder::createCollection(pcCands, "p"));
0859     else if (useGP)
0860       phi = new BPHPhiToKKBuilder(
0861           es, BPHRecoBuilder::createCollection(gpCands, "h"), BPHRecoBuilder::createCollection(gpCands, "h"));
0862   }
0863 
0864   set<BPHRecoConstCandPtr> sPhi;
0865 
0866   if (phi != nullptr) {
0867     rIter = parMap.find(Pkk);
0868     if (rIter != rIend) {
0869       const map<parType, double>& pMap = rIter->second;
0870       map<parType, double>::const_iterator pIter = pMap.begin();
0871       map<parType, double>::const_iterator pIend = pMap.end();
0872       while (pIter != pIend) {
0873         const map<parType, double>::value_type& pEntry = *pIter++;
0874         parType id = pEntry.first;
0875         double pv = pEntry.second;
0876         switch (id) {
0877           case ptMin:
0878             phi->setPtMin(pv);
0879             break;
0880           case etaMax:
0881             phi->setEtaMax(pv);
0882             break;
0883           case massMin:
0884             phi->setMassMin(pv);
0885             break;
0886           case massMax:
0887             phi->setMassMax(pv);
0888             break;
0889           case probMin:
0890             phi->setProbMin(pv);
0891             break;
0892           case writeCandidate:
0893             writePkk = (pv > 0);
0894             break;
0895           default:
0896             break;
0897         }
0898       }
0899     }
0900     lPhi = phi->build();
0901     if (allPkk)
0902       sPhi.insert(lPhi.begin(), lPhi.end());
0903     delete phi;
0904   }
0905 
0906   bool phiFound = !lPhi.empty();
0907 
0908   // build and dump Bs
0909 
0910   if (recoBs && jPsiFound && phiFound) {
0911     BPHBsToJPsiPhiBuilder* bs = new BPHBsToJPsiPhiBuilder(es, lJPsi, lPhi);
0912     rIter = parMap.find(Bs);
0913     if (rIter != rIend) {
0914       const map<parType, double>& pMap = rIter->second;
0915       map<parType, double>::const_iterator pIter = pMap.begin();
0916       map<parType, double>::const_iterator pIend = pMap.end();
0917       while (pIter != pIend) {
0918         const map<parType, double>::value_type& pEntry = *pIter++;
0919         parType id = pEntry.first;
0920         double pv = pEntry.second;
0921         switch (id) {
0922           case mPsiMin:
0923             bs->setJPsiMassMin(pv);
0924             break;
0925           case mPsiMax:
0926             bs->setJPsiMassMax(pv);
0927             break;
0928           case mPhiMin:
0929             bs->setPhiMassMin(pv);
0930             break;
0931           case mPhiMax:
0932             bs->setPhiMassMax(pv);
0933             break;
0934           case massMin:
0935             bs->setMassMin(pv);
0936             break;
0937           case massMax:
0938             bs->setMassMax(pv);
0939             break;
0940           case probMin:
0941             bs->setProbMin(pv);
0942             break;
0943           case mFitMin:
0944             bs->setMassFitMin(pv);
0945             break;
0946           case mFitMax:
0947             bs->setMassFitMax(pv);
0948             break;
0949           case constrMJPsi:
0950             bs->setConstr(pv > 0);
0951             break;
0952           case writeCandidate:
0953             writeBs = (pv > 0);
0954             break;
0955           default:
0956             break;
0957         }
0958       }
0959     }
0960 
0961     lBs = bs->build();
0962     delete bs;
0963 
0964     int iBs;
0965     int nBs = lBs.size();
0966     for (iBs = 0; iBs < nBs; ++iBs)
0967       sPhi.insert(lBs[iBs]->getComp("Phi"));
0968   }
0969   set<BPHRecoConstCandPtr>::const_iterator phi_iter = sPhi.begin();
0970   set<BPHRecoConstCandPtr>::const_iterator phi_iend = sPhi.end();
0971   lSs.reserve(sPhi.size());
0972   while (phi_iter != phi_iend)
0973     lSs.push_back(*phi_iter++);
0974 
0975   // build K0
0976 
0977   BPHK0sToPiPiBuilder* k0s = nullptr;
0978   if (recoK0s && (jPsiFound || allK0s)) {
0979     if (useK0)
0980       k0s = new BPHK0sToPiPiBuilder(es, k0Cand.product(), "cfp");
0981     else if (useKS)
0982       k0s = new BPHK0sToPiPiBuilder(es, kSCand.product(), "cfp");
0983   }
0984   if (k0s != nullptr) {
0985     rIter = parMap.find(K0s);
0986     if (rIter != rIend) {
0987       const map<parType, double>& pMap = rIter->second;
0988       map<parType, double>::const_iterator pIter = pMap.begin();
0989       map<parType, double>::const_iterator pIend = pMap.end();
0990       while (pIter != pIend) {
0991         const map<parType, double>::value_type& pEntry = *pIter++;
0992         parType id = pEntry.first;
0993         double pv = pEntry.second;
0994         switch (id) {
0995           case ptMin:
0996             k0s->setPtMin(pv);
0997             break;
0998           case etaMax:
0999             k0s->setEtaMax(pv);
1000             break;
1001           case massMin:
1002             k0s->setMassMin(pv);
1003             break;
1004           case massMax:
1005             k0s->setMassMax(pv);
1006             break;
1007           case probMin:
1008             k0s->setProbMin(pv);
1009             break;
1010           case writeCandidate:
1011             writeK0s = (pv > 0);
1012             break;
1013           default:
1014             break;
1015         }
1016       }
1017     }
1018     lK0 = k0s->build();
1019     delete k0s;
1020   }
1021 
1022   bool k0Found = !lK0.empty();
1023 
1024   // build Lambda0
1025 
1026   BPHLambda0ToPPiBuilder* l0s = nullptr;
1027   if (recoLambda0 && (jPsiFound || allLambda0)) {
1028     if (useL0)
1029       l0s = new BPHLambda0ToPPiBuilder(es, l0Cand.product(), "cfp");
1030     else if (useLS)
1031       l0s = new BPHLambda0ToPPiBuilder(es, lSCand.product(), "cfp");
1032   }
1033   if (l0s != nullptr) {
1034     rIter = parMap.find(Lambda0);
1035     if (rIter != rIend) {
1036       const map<parType, double>& pMap = rIter->second;
1037       map<parType, double>::const_iterator pIter = pMap.begin();
1038       map<parType, double>::const_iterator pIend = pMap.end();
1039       while (pIter != pIend) {
1040         const map<parType, double>::value_type& pEntry = *pIter++;
1041         parType id = pEntry.first;
1042         double pv = pEntry.second;
1043         switch (id) {
1044           case ptMin:
1045             l0s->setPtMin(pv);
1046             break;
1047           case etaMax:
1048             l0s->setEtaMax(pv);
1049             break;
1050           case massMin:
1051             l0s->setMassMin(pv);
1052             break;
1053           case massMax:
1054             l0s->setMassMax(pv);
1055             break;
1056           case probMin:
1057             l0s->setProbMin(pv);
1058             break;
1059           case writeCandidate:
1060             writeLambda0 = (pv > 0);
1061             break;
1062           default:
1063             break;
1064         }
1065       }
1066     }
1067     lL0 = l0s->build();
1068     delete l0s;
1069   }
1070 
1071   bool l0Found = !lL0.empty();
1072 
1073   // build and dump Bd -> JPsi K0s
1074 
1075   if (recoB0 && jPsiFound && k0Found) {
1076     BPHBdToJPsiKsBuilder* b0 = new BPHBdToJPsiKsBuilder(es, lJPsi, lK0);
1077     rIter = parMap.find(B0);
1078     if (rIter != rIend) {
1079       const map<parType, double>& pMap = rIter->second;
1080       map<parType, double>::const_iterator pIter = pMap.begin();
1081       map<parType, double>::const_iterator pIend = pMap.end();
1082       while (pIter != pIend) {
1083         const map<parType, double>::value_type& pEntry = *pIter++;
1084         parType id = pEntry.first;
1085         double pv = pEntry.second;
1086         switch (id) {
1087           case mPsiMin:
1088             b0->setJPsiMassMin(pv);
1089             break;
1090           case mPsiMax:
1091             b0->setJPsiMassMax(pv);
1092             break;
1093           case mK0sMin:
1094             b0->setK0MassMin(pv);
1095             break;
1096           case mK0sMax:
1097             b0->setK0MassMax(pv);
1098             break;
1099           case massMin:
1100             b0->setMassMin(pv);
1101             break;
1102           case massMax:
1103             b0->setMassMax(pv);
1104             break;
1105           case probMin:
1106             b0->setProbMin(pv);
1107             break;
1108           case mFitMin:
1109             b0->setMassFitMin(pv);
1110             break;
1111           case mFitMax:
1112             b0->setMassFitMax(pv);
1113             break;
1114           case constrMJPsi:
1115             b0->setConstr(pv > 0);
1116             break;
1117           case writeCandidate:
1118             writeB0 = (pv > 0);
1119             break;
1120           default:
1121             break;
1122         }
1123       }
1124     }
1125 
1126     lB0 = b0->build();
1127     const map<const BPHRecoCandidate*, const BPHRecoCandidate*>& b0Map = b0->daughMap();
1128     daughMap.insert(b0Map.begin(), b0Map.end());
1129     delete b0;
1130   }
1131 
1132   // build and dump Lambdab -> JPsi Lambda0
1133 
1134   if (recoLambdab && jPsiFound && l0Found) {
1135     BPHLbToJPsiL0Builder* lb = new BPHLbToJPsiL0Builder(es, lJPsi, lL0);
1136     rIter = parMap.find(Lambdab);
1137     if (rIter != rIend) {
1138       const map<parType, double>& pMap = rIter->second;
1139       map<parType, double>::const_iterator pIter = pMap.begin();
1140       map<parType, double>::const_iterator pIend = pMap.end();
1141       while (pIter != pIend) {
1142         const map<parType, double>::value_type& pEntry = *pIter++;
1143         parType id = pEntry.first;
1144         double pv = pEntry.second;
1145         switch (id) {
1146           case mPsiMin:
1147             lb->setJPsiMassMin(pv);
1148             break;
1149           case mPsiMax:
1150             lb->setJPsiMassMax(pv);
1151             break;
1152           case mLambda0Min:
1153             lb->setLambda0MassMin(pv);
1154             break;
1155           case mLambda0Max:
1156             lb->setLambda0MassMax(pv);
1157             break;
1158           case massMin:
1159             lb->setMassMin(pv);
1160             break;
1161           case massMax:
1162             lb->setMassMax(pv);
1163             break;
1164           case probMin:
1165             lb->setProbMin(pv);
1166             break;
1167           case mFitMin:
1168             lb->setMassFitMin(pv);
1169             break;
1170           case mFitMax:
1171             lb->setMassFitMax(pv);
1172             break;
1173           case constrMJPsi:
1174             lb->setConstr(pv > 0);
1175             break;
1176           case writeCandidate:
1177             writeLambdab = (pv > 0);
1178             break;
1179           default:
1180             break;
1181         }
1182       }
1183     }
1184 
1185     lLb = lb->build();
1186     const map<const BPHRecoCandidate*, const BPHRecoCandidate*>& ldMap = lb->daughMap();
1187     daughMap.insert(ldMap.begin(), ldMap.end());
1188     delete lb;
1189   }
1190 
1191   // build and dump Bc
1192 
1193   BPHBcToJPsiPiBuilder* bc = nullptr;
1194   if (recoBc && jPsiFound) {
1195     if (usePF)
1196       bc = new BPHBcToJPsiPiBuilder(es, lJPsi, BPHRecoBuilder::createCollection(pfCands, "f"));
1197     else if (usePC)
1198       bc = new BPHBcToJPsiPiBuilder(es, lJPsi, BPHRecoBuilder::createCollection(pcCands, "p"));
1199     else if (useGP)
1200       bc = new BPHBcToJPsiPiBuilder(es, lJPsi, BPHRecoBuilder::createCollection(gpCands, "h"));
1201   }
1202 
1203   if (bc != nullptr) {
1204     rIter = parMap.find(Bc);
1205     if (rIter != rIend) {
1206       const map<parType, double>& pMap = rIter->second;
1207       map<parType, double>::const_iterator pIter = pMap.begin();
1208       map<parType, double>::const_iterator pIend = pMap.end();
1209       while (pIter != pIend) {
1210         const map<parType, double>::value_type& pEntry = *pIter++;
1211         parType id = pEntry.first;
1212         double pv = pEntry.second;
1213         switch (id) {
1214           case ptMin:
1215             bc->setPiPtMin(pv);
1216             break;
1217           case etaMax:
1218             bc->setPiEtaMax(pv);
1219             break;
1220           case mPsiMin:
1221             bc->setJPsiMassMin(pv);
1222             break;
1223           case mPsiMax:
1224             bc->setJPsiMassMax(pv);
1225             break;
1226           case massMin:
1227             bc->setMassMin(pv);
1228             break;
1229           case massMax:
1230             bc->setMassMax(pv);
1231             break;
1232           case probMin:
1233             bc->setProbMin(pv);
1234             break;
1235           case mFitMin:
1236             bc->setMassFitMin(pv);
1237             break;
1238           case mFitMax:
1239             bc->setMassFitMax(pv);
1240             break;
1241           case constrMJPsi:
1242             bc->setConstr(pv > 0);
1243             break;
1244           case writeCandidate:
1245             writeBc = (pv > 0);
1246             break;
1247           default:
1248             break;
1249         }
1250       }
1251     }
1252     lBc = bc->build();
1253     delete bc;
1254   }
1255 
1256   // build and dump Psi2S
1257 
1258   BPHPsi2SToJPsiPiPiBuilder* psi2S = nullptr;
1259   if (recoPsi2S && jPsiFound) {
1260     if (usePF)
1261       psi2S = new BPHPsi2SToJPsiPiPiBuilder(
1262           es, lJPsi, BPHRecoBuilder::createCollection(pfCands, "f"), BPHRecoBuilder::createCollection(pfCands, "f"));
1263     else if (usePC)
1264       psi2S = new BPHPsi2SToJPsiPiPiBuilder(
1265           es, lJPsi, BPHRecoBuilder::createCollection(pcCands, "p"), BPHRecoBuilder::createCollection(pcCands, "p"));
1266     else if (useGP)
1267       psi2S = new BPHPsi2SToJPsiPiPiBuilder(
1268           es, lJPsi, BPHRecoBuilder::createCollection(gpCands, "h"), BPHRecoBuilder::createCollection(gpCands, "h"));
1269   }
1270 
1271   if (psi2S != nullptr) {
1272     rIter = parMap.find(Psi2S);
1273     if (rIter != rIend) {
1274       const map<parType, double>& pMap = rIter->second;
1275       map<parType, double>::const_iterator pIter = pMap.begin();
1276       map<parType, double>::const_iterator pIend = pMap.end();
1277       while (pIter != pIend) {
1278         const map<parType, double>::value_type& pEntry = *pIter++;
1279         parType id = pEntry.first;
1280         double pv = pEntry.second;
1281         switch (id) {
1282           case ptMin:
1283             psi2S->setPiPtMin(pv);
1284             break;
1285           case etaMax:
1286             psi2S->setPiEtaMax(pv);
1287             break;
1288           case mPsiMin:
1289             psi2S->setJPsiMassMin(pv);
1290             break;
1291           case mPsiMax:
1292             psi2S->setJPsiMassMax(pv);
1293             break;
1294           case massMin:
1295             psi2S->setMassMin(pv);
1296             break;
1297           case massMax:
1298             psi2S->setMassMax(pv);
1299             break;
1300           case probMin:
1301             psi2S->setProbMin(pv);
1302             break;
1303           case mFitMin:
1304             psi2S->setMassFitMin(pv);
1305             break;
1306           case mFitMax:
1307             psi2S->setMassFitMax(pv);
1308             break;
1309           case constrMJPsi:
1310             psi2S->setConstr(pv > 0);
1311             break;
1312           case writeCandidate:
1313             writePsi2S = (pv > 0);
1314             break;
1315           default:
1316             break;
1317         }
1318       }
1319     }
1320     lPsi2S = psi2S->build();
1321     delete psi2S;
1322   }
1323 
1324   // build and dump X3872
1325 
1326   BPHX3872ToJPsiPiPiBuilder* x3872 = nullptr;
1327   if (recoX3872 && jPsiFound) {
1328     if (usePF)
1329       x3872 = new BPHX3872ToJPsiPiPiBuilder(
1330           es, lJPsi, BPHRecoBuilder::createCollection(pfCands, "f"), BPHRecoBuilder::createCollection(pfCands, "f"));
1331     else if (usePC)
1332       x3872 = new BPHX3872ToJPsiPiPiBuilder(
1333           es, lJPsi, BPHRecoBuilder::createCollection(pcCands, "p"), BPHRecoBuilder::createCollection(pcCands, "p"));
1334     else if (useGP)
1335       x3872 = new BPHX3872ToJPsiPiPiBuilder(
1336           es, lJPsi, BPHRecoBuilder::createCollection(gpCands, "h"), BPHRecoBuilder::createCollection(gpCands, "h"));
1337   }
1338 
1339   if (x3872 != nullptr) {
1340     rIter = parMap.find(X3872);
1341     if (rIter != rIend) {
1342       const map<parType, double>& pMap = rIter->second;
1343       map<parType, double>::const_iterator pIter = pMap.begin();
1344       map<parType, double>::const_iterator pIend = pMap.end();
1345       while (pIter != pIend) {
1346         const map<parType, double>::value_type& pEntry = *pIter++;
1347         parType id = pEntry.first;
1348         double pv = pEntry.second;
1349         switch (id) {
1350           case ptMin:
1351             x3872->setPiPtMin(pv);
1352             break;
1353           case etaMax:
1354             x3872->setPiEtaMax(pv);
1355             break;
1356           case mPsiMin:
1357             x3872->setJPsiMassMin(pv);
1358             break;
1359           case mPsiMax:
1360             x3872->setJPsiMassMax(pv);
1361             break;
1362           case massMin:
1363             x3872->setMassMin(pv);
1364             break;
1365           case massMax:
1366             x3872->setMassMax(pv);
1367             break;
1368           case probMin:
1369             x3872->setProbMin(pv);
1370             break;
1371           case mFitMin:
1372             x3872->setMassFitMin(pv);
1373             break;
1374           case mFitMax:
1375             x3872->setMassFitMax(pv);
1376             break;
1377           case constrMJPsi:
1378             x3872->setConstr(pv > 0);
1379             break;
1380           case writeCandidate:
1381             writeX3872 = (pv > 0);
1382             break;
1383           default:
1384             break;
1385         }
1386       }
1387     }
1388     lX3872 = x3872->build();
1389     delete x3872;
1390   }
1391 
1392   // merge Psi2S and X3872
1393   class ResTrkTrkCompare {
1394   public:
1395     bool operator()(const BPHRecoConstCandPtr& l, const BPHRecoConstCandPtr& r) const {
1396       vector<const reco::Track*> tl = l->tracks();
1397       vector<const reco::Track*> tr = r->tracks();
1398       if (tl.size() < tr.size())
1399         return true;
1400       sort(tl.begin(), tl.end());
1401       sort(tr.begin(), tr.end());
1402       int n = tr.size();
1403       int i;
1404       for (i = 0; i < n; ++i) {
1405         if (tl[i] < tr[i])
1406           return true;
1407         if (tl[i] > tr[i])
1408           return false;
1409       }
1410       return false;
1411     }
1412   } rttc;
1413   set<BPHRecoConstCandPtr, ResTrkTrkCompare> sjpPiPi(rttc);
1414   sjpPiPi.insert(lPsi2S.begin(), lPsi2S.end());
1415   sjpPiPi.insert(lX3872.begin(), lX3872.end());
1416   vector<BPHRecoConstCandPtr> ljpPiPi;
1417   ljpPiPi.insert(ljpPiPi.end(), sjpPiPi.begin(), sjpPiPi.end());
1418   bool jpPiPiFound = !ljpPiPi.empty();
1419 
1420   // build and dump Bp
1421 
1422   BPHBuToPsi2SKBuilder* bp = nullptr;
1423   if (recoBp && jpPiPiFound) {
1424     if (usePF)
1425       bp = new BPHBuToPsi2SKBuilder(es, ljpPiPi, BPHRecoBuilder::createCollection(pfCands, "f"));
1426     else if (usePC)
1427       bp = new BPHBuToPsi2SKBuilder(es, ljpPiPi, BPHRecoBuilder::createCollection(pcCands, "p"));
1428     else if (useGP)
1429       bp = new BPHBuToPsi2SKBuilder(es, ljpPiPi, BPHRecoBuilder::createCollection(gpCands, "h"));
1430   }
1431 
1432   if (bp != nullptr) {
1433     class BPHBuToPsi2SSelect : public BPHMassFitSelect {
1434     public:
1435       BPHBuToPsi2SSelect()
1436           : BPHMassFitSelect("Psi2S", BPHParticleMasses::psi2Mass, BPHParticleMasses::psi2MWidth, 5.0, 6.0) {}
1437       ~BPHBuToPsi2SSelect() override = default;
1438       bool accept(const BPHKinematicFit& cand) const override {
1439         const_cast<BPHRecoCandidate*>(cand.getComp("Psi2S").get())
1440             ->setIndependentFit("JPsi", true, BPHParticleMasses::jPsiMass, BPHParticleMasses::jPsiMWidth);
1441         return BPHMassFitSelect::accept(cand);
1442       }
1443     };
1444     bool mcJPsi = false;
1445     bool mcPsi2 = true;
1446     rIter = parMap.find(Bp);
1447     if (rIter != rIend) {
1448       const map<parType, double>& pMap = rIter->second;
1449       map<parType, double>::const_iterator pIter = pMap.begin();
1450       map<parType, double>::const_iterator pIend = pMap.end();
1451       while (pIter != pIend) {
1452         const map<parType, double>::value_type& pEntry = *pIter++;
1453         parType id = pEntry.first;
1454         double pv = pEntry.second;
1455         switch (id) {
1456           case ptMin:
1457             bp->setKPtMin(pv);
1458             break;
1459           case etaMax:
1460             bp->setKEtaMax(pv);
1461             break;
1462           case mPsiMin:
1463             bp->setPsi2SMassMin(pv);
1464             break;
1465           case mPsiMax:
1466             bp->setPsi2SMassMax(pv);
1467             break;
1468           case massMin:
1469             bp->setMassMin(pv);
1470             break;
1471           case massMax:
1472             bp->setMassMax(pv);
1473             break;
1474           case probMin:
1475             bp->setProbMin(pv);
1476             break;
1477           case mFitMin:
1478             bp->setMassFitMin(pv);
1479             break;
1480           case mFitMax:
1481             bp->setMassFitMax(pv);
1482             break;
1483           case constrMJPsi:
1484             mcJPsi = (pv > 0);
1485             break;
1486           case constrMPsi2:
1487             mcPsi2 = (pv > 0);
1488             break;
1489           case writeCandidate:
1490             writeBp = (pv > 0);
1491             break;
1492           default:
1493             break;
1494         }
1495       }
1496     }
1497     if (mcJPsi)
1498       bp->setMassFitSelect(mcPsi2 ? new BPHBuToPsi2SSelect
1499                                   : new BPHMassFitSelect("Psi2S/JPsi",
1500                                                          BPHParticleMasses::jPsiMass,
1501                                                          BPHParticleMasses::jPsiMWidth,
1502                                                          bp->getMassFitMin(),
1503                                                          bp->getMassFitMax()));
1504     else
1505       bp->setConstr(mcPsi2);
1506     lBp = bp->build();
1507     const map<const BPHRecoCandidate*, const BPHRecoCandidate*>& bpMap = bp->daughMap();
1508     daughMap.insert(bpMap.begin(), bpMap.end());
1509     delete bp;
1510   }
1511 
1512   return;
1513 }
1514 
1515 void BPHWriteSpecificDecay::setRecoParameters(const edm::ParameterSet& ps) {
1516   const string& name = ps.getParameter<string>("name");
1517   bool writeCandidate = ps.getParameter<bool>("writeCandidate");
1518   switch (rMap[name]) {
1519     case Onia:
1520       recoOnia = true;
1521       writeOnia = writeCandidate;
1522       break;
1523     case Pmm:
1524     case Psi1:
1525     case Psi2:
1526     case Ups:
1527     case Ups1:
1528     case Ups2:
1529     case Ups3:
1530       recoOnia = true;
1531       break;
1532     case Kx0:
1533       recoKx0 = true;
1534       allKx0 = false;
1535       writeKx0 = writeCandidate;
1536       break;
1537     case Pkk:
1538       recoPkk = true;
1539       allPkk = false;
1540       writePkk = writeCandidate;
1541       break;
1542     case Bu:
1543       recoBu = true;
1544       writeBu = writeCandidate;
1545       break;
1546     case Bp:
1547       recoBp = true;
1548       writeBp = writeCandidate;
1549       break;
1550     case Bd:
1551       recoBd = true;
1552       writeBd = writeCandidate;
1553       break;
1554     case Bs:
1555       recoBs = true;
1556       writeBs = writeCandidate;
1557       break;
1558     case K0s:
1559       recoK0s = true;
1560       allK0s = false;
1561       writeK0s = writeCandidate;
1562       break;
1563     case Lambda0:
1564       recoLambda0 = true;
1565       allLambda0 = false;
1566       writeLambda0 = writeCandidate;
1567       break;
1568     case B0:
1569       recoB0 = true;
1570       writeB0 = writeCandidate;
1571       break;
1572     case Lambdab:
1573       recoLambdab = true;
1574       writeLambdab = writeCandidate;
1575       break;
1576     case Bc:
1577       recoBc = true;
1578       writeBc = writeCandidate;
1579       break;
1580     case Psi2S:
1581       recoPsi2S = true;
1582       writePsi2S = writeCandidate;
1583       break;
1584     case X3872:
1585       recoX3872 = true;
1586       writeX3872 = writeCandidate;
1587       break;
1588   }
1589 
1590   map<string, parType>::const_iterator pIter = pMap.begin();
1591   map<string, parType>::const_iterator pIend = pMap.end();
1592   while (pIter != pIend) {
1593     const map<string, parType>::value_type& entry = *pIter++;
1594     const string& pn = entry.first;
1595     parType id = entry.second;
1596     double pv = ps.getParameter<double>(pn);
1597     if (pv > -1.0e35)
1598       edm::LogVerbatim("Configuration") << "BPHWriteSpecificDecay::setRecoParameters: set " << pn << " for " << name
1599                                         << " : " << (parMap[rMap[name]][id] = pv);
1600   }
1601 
1602   map<string, parType>::const_iterator fIter = fMap.begin();
1603   map<string, parType>::const_iterator fIend = fMap.end();
1604   while (fIter != fIend) {
1605     const map<string, parType>::value_type& entry = *fIter++;
1606     const string& fn = entry.first;
1607     parType id = entry.second;
1608     double pv = (ps.getParameter<bool>(fn) ? 1 : -1);
1609     if (pv > -1.0e35)
1610       edm::LogVerbatim("Configuration") << "BPHWriteSpecificDecay::setRecoParameters: set " << fn << " for " << name
1611                                         << " : " << (parMap[rMap[name]][id] = pv);
1612   }
1613 
1614   return;
1615 }
1616 
1617 void BPHWriteSpecificDecay::addTrackModes(const std::string& name,
1618                                           const BPHRecoCandidate& cand,
1619                                           std::string& modes,
1620                                           bool& count) {
1621   for (const std::map<std::string, const reco::Candidate*>::value_type& entry : cand.daugMap()) {
1622     if (count)
1623       modes += "#";
1624     modes += (name + entry.first + ":" + cand.getTMode(entry.second));
1625     count = true;
1626   }
1627   for (const std::map<std::string, BPHRecoConstCandPtr>::value_type& entry : cand.compMap()) {
1628     addTrackModes(entry.first + "/", *entry.second, modes, count);
1629   }
1630   return;
1631 }
1632 
1633 void BPHWriteSpecificDecay::addTrackModes(const std::string& name,
1634                                           const BPHRecoCandidate& cand,
1635                                           pat::CompositeCandidate& cc) {
1636   for (const std::map<std::string, const reco::Candidate*>::value_type& entry : cand.daugMap())
1637     cc.addUserData(name + entry.first, string(1, cand.getTMode(entry.second)), true);
1638   for (const std::map<std::string, BPHRecoConstCandPtr>::value_type& entry : cand.compMap())
1639     addTrackModes(name + entry.first + "/", *entry.second, cc);
1640   return;
1641 }
1642 
1643 #include "FWCore/Framework/interface/MakerMacros.h"
1644 
1645 DEFINE_FWK_MODULE(BPHWriteSpecificDecay);
[ 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!