Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 #include "HeavyFlavorAnalysis/SpecificDecay/test/stubs/TestBPHSpecificDecay.h"
0002 
0003 #include "FWCore/Framework/interface/MakerMacros.h"
0004 
0005 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHRecoBuilder.h"
0006 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHRecoSelect.h"
0007 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHRecoCandidate.h"
0008 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHPlusMinusCandidate.h"
0009 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHMomentumSelect.h"
0010 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHVertexSelect.h"
0011 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHTrackReference.h"
0012 
0013 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHMuonPtSelect.h"
0014 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHMuonEtaSelect.h"
0015 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHParticlePtSelect.h"
0016 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHParticleNeutralVeto.h"
0017 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHMultiSelect.h"
0018 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHMassSelect.h"
0019 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHChi2Select.h"
0020 
0021 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHOniaToMuMuBuilder.h"
0022 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHKx0ToKPiBuilder.h"
0023 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHPhiToKKBuilder.h"
0024 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHBuToJPsiKBuilder.h"
0025 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHBsToJPsiPhiBuilder.h"
0026 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHBdToJPsiKxBuilder.h"
0027 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHParticleMasses.h"
0028 
0029 #include "DataFormats/PatCandidates/interface/Muon.h"
0030 #include "DataFormats/TrackReco/interface/Track.h"
0031 
0032 #include "DataFormats/PatCandidates/interface/GenericParticle.h"
0033 #include "DataFormats/PatCandidates/interface/CompositeCandidate.h"
0034 
0035 #include <TH1.h>
0036 #include <TFile.h>
0037 
0038 #include <set>
0039 #include <string>
0040 #include <iostream>
0041 #include <fstream>
0042 
0043 using namespace std;
0044 
0045 #define SET_LABEL(NAME, PSET) (NAME = PSET.getParameter<string>(#NAME))
0046 // SET_LABEL(xyz,ps);
0047 // is equivalent to
0048 // xyz = ps.getParameter<string>( "xyx" )
0049 
0050 TestBPHSpecificDecay::TestBPHSpecificDecay(const edm::ParameterSet& ps) {
0051   usePM = (!SET_LABEL(patMuonLabel, ps).empty());
0052   useCC = (!SET_LABEL(ccCandsLabel, ps).empty());
0053   usePF = (!SET_LABEL(pfCandsLabel, ps).empty());
0054   usePC = (!SET_LABEL(pcCandsLabel, ps).empty());
0055   useGP = (!SET_LABEL(gpCandsLabel, ps).empty());
0056 
0057   if (usePM)
0058     consume<pat::MuonCollection>(patMuonToken, patMuonLabel);
0059   if (useCC)
0060     consume<vector<pat::CompositeCandidate>>(ccCandsToken, ccCandsLabel);
0061   if (usePF)
0062     consume<vector<reco::PFCandidate>>(pfCandsToken, pfCandsLabel);
0063   if (usePC)
0064     consume<vector<BPHTrackReference::candidate>>(pcCandsToken, pcCandsLabel);
0065   if (useGP)
0066     consume<vector<pat::GenericParticle>>(gpCandsToken, gpCandsLabel);
0067 
0068   SET_LABEL(outDump, ps);
0069   SET_LABEL(outHist, ps);
0070   if (outDump.empty())
0071     fPtr = &cout;
0072   else
0073     fPtr = new ofstream(outDump.c_str());
0074 }
0075 
0076 TestBPHSpecificDecay::~TestBPHSpecificDecay() {}
0077 
0078 void TestBPHSpecificDecay::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0079   edm::ParameterSetDescription desc;
0080   desc.add<string>("patMuonLabel", "");
0081   desc.add<string>("ccCandsLabel", "");
0082   desc.add<string>("pfCandsLabel", "");
0083   desc.add<string>("pcCandsLabel", "");
0084   desc.add<string>("gpCandsLabel", "");
0085   desc.add<string>("outDump", "dump.txt");
0086   desc.add<string>("outHist", "hist.root");
0087   descriptions.add("testBPHSpecificDecay", desc);
0088   return;
0089 }
0090 
0091 void TestBPHSpecificDecay::beginJob() {
0092   *fPtr << "TestBPHSpecificDecay::beginJob" << endl;
0093   createHisto("massJPsi", 60, 2.95, 3.25);   // JPsi mass
0094   createHisto("mcstJPsi", 60, 2.95, 3.25);   // JPsi mass, with constraint
0095   createHisto("massKx0", 50, 0.80, 1.05);    // Kx0  mass
0096   createHisto("massPhi", 40, 1.01, 1.03);    // Phi  mass
0097   createHisto("massBu", 20, 5.00, 5.50);     // Bu   mass
0098   createHisto("mcstBu", 20, 5.00, 5.50);     // Bu   mass, with constraint
0099   createHisto("massBd", 20, 5.00, 5.50);     // Bd   mass
0100   createHisto("mcstBd", 20, 5.00, 5.50);     // Bd   mass, with constraint
0101   createHisto("massBs", 20, 5.10, 5.60);     // Bs   mass
0102   createHisto("mcstBs", 20, 5.10, 5.60);     // Bs   mass, with constraint
0103   createHisto("massBsPhi", 50, 1.01, 1.03);  // Phi  mass in Bs decay
0104   createHisto("massBdKx0", 50, 0.80, 1.05);  // Kx0  mass in Bd decay
0105 
0106   createHisto("massFull", 200, 2.00, 12.00);  // Full mass
0107   createHisto("massFsel", 200, 2.00, 12.00);  // Full mass
0108   createHisto("massPhi", 70, 0.85, 1.20);     // Psi1 mass
0109   createHisto("massPsi1", 60, 2.95, 3.25);    // Psi1 mass
0110   createHisto("massPsi2", 50, 3.55, 3.80);    // Psi2 mass
0111   createHisto("massUps1", 130, 9.10, 9.75);   // Ups1 mass
0112   createHisto("massUps2", 90, 9.75, 10.20);   // Ups2 mass
0113   createHisto("massUps3", 80, 10.20, 10.60);  // Ups3 mass
0114 
0115   return;
0116 }
0117 
0118 void TestBPHSpecificDecay::analyze(const edm::Event& ev, const edm::EventSetup& es) {
0119   ostream& outF = *fPtr;
0120   outF << "--------- event " << ev.id().run() << " / " << ev.id().event() << " ---------" << endl;
0121 
0122   // get object collections
0123   // collections are got through "BPHTokenWrapper" interface to allow
0124   // uniform access in different CMSSW versions
0125 
0126   int nrc = 0;
0127 
0128   // get reco::PFCandidate collection (in full AOD )
0129   edm::Handle<vector<reco::PFCandidate>> pfCands;
0130   if (usePF) {
0131     pfCandsToken.get(ev, pfCands);
0132     nrc = pfCands->size();
0133   }
0134 
0135   // get pat::PackedCandidate collection (in MiniAOD)
0136   // pat::PackedCandidate is not defined in CMSSW_5XY, so a
0137   // typedef (BPHTrackReference::candidate) is used, actually referring
0138   // to pat::PackedCandidate only for CMSSW versions where it's defined
0139   edm::Handle<vector<BPHTrackReference::candidate>> pcCands;
0140   if (usePC) {
0141     pcCandsToken.get(ev, pcCands);
0142     nrc = pcCands->size();
0143   }
0144 
0145   // get pat::GenericParticle collection (in skimmed data)
0146   edm::Handle<vector<pat::GenericParticle>> gpCands;
0147   if (useGP) {
0148     gpCandsToken.get(ev, gpCands);
0149     nrc = gpCands->size();
0150   }
0151 
0152   // get pat::Muon collection (in full AOD and MiniAOD)
0153   edm::Handle<pat::MuonCollection> patMuon;
0154   if (usePM) {
0155     patMuonToken.get(ev, patMuon);
0156   }
0157 
0158   // get muons from pat::CompositeCandidate objects describing onia;
0159   // muons from all composite objects are copied to an unique std::vector
0160   vector<const reco::Candidate*> muDaugs;
0161   set<const pat::Muon*> muonSet;
0162   if (useCC) {
0163     edm::Handle<vector<pat::CompositeCandidate>> ccCands;
0164     ccCandsToken.get(ev, ccCands);
0165     int n = ccCands->size();
0166     muDaugs.clear();
0167     muDaugs.reserve(n);
0168     muonSet.clear();
0169     set<const pat::Muon*>::const_iterator iter;
0170     set<const pat::Muon*>::const_iterator iend;
0171     int i;
0172     for (i = 0; i < n; ++i) {
0173       const pat::CompositeCandidate& cc = ccCands->at(i);
0174       int j;
0175       int m = cc.numberOfDaughters();
0176       for (j = 0; j < m; ++j) {
0177         const reco::Candidate* dp = cc.daughter(j);
0178         const pat::Muon* mp = dynamic_cast<const pat::Muon*>(dp);
0179         iter = muonSet.begin();
0180         iend = muonSet.end();
0181         bool add = (mp != nullptr) && (muonSet.find(mp) == iend);
0182         while (add && (iter != iend)) {
0183           if (BPHRecoBuilder::sameTrack(mp, *iter++, 1.0e-5))
0184             add = false;
0185         }
0186         if (add)
0187           muonSet.insert(mp);
0188       }
0189     }
0190     iter = muonSet.begin();
0191     iend = muonSet.end();
0192     while (iter != iend)
0193       muDaugs.push_back(*iter++);
0194   }
0195 
0196   // reconstruct resonances
0197 
0198   outF << "build and dump full onia" << endl;
0199   BPHOniaToMuMuBuilder* onia = nullptr;
0200   if (usePM)
0201     onia = new BPHOniaToMuMuBuilder(
0202         es, BPHRecoBuilder::createCollection(patMuon, "cfmig"), BPHRecoBuilder::createCollection(patMuon, "cfmig"));
0203   else if (useCC)
0204     onia = new BPHOniaToMuMuBuilder(
0205         es, BPHRecoBuilder::createCollection(muDaugs, "cfmig"), BPHRecoBuilder::createCollection(muDaugs, "cfmig"));
0206 
0207   vector<BPHPlusMinusConstCandPtr> lFull = onia->build();
0208   int iFull;
0209   int nFull = lFull.size();
0210   outF << nFull << " onia cand found" << endl;
0211   for (iFull = 0; iFull < nFull; ++iFull)
0212     dumpRecoCand("Full", lFull[iFull].get());
0213 
0214   BPHMuonPtSelect ptSel1(3.0);
0215   BPHMuonPtSelect ptSel2(3.5);
0216   BPHMuonPtSelect ptSel3(4.5);
0217   BPHMuonEtaSelect etaSel1(1.6);
0218   BPHMuonEtaSelect etaSel2(1.4);
0219   BPHMuonEtaSelect etaSel3(1.2);
0220   BPHMultiSelect<BPHSlimSelect<BPHRecoSelect>> select1(BPHSelectOperation::and_mode);
0221   select1.include(ptSel1);
0222   select1.include(etaSel1);
0223   select1.include(etaSel2, false);
0224   BPHMultiSelect<BPHSlimSelect<BPHRecoSelect>> select2(BPHSelectOperation::and_mode);
0225   select2.include(ptSel2);
0226   select2.include(etaSel2);
0227   select2.include(etaSel3, false);
0228   BPHMultiSelect<BPHSlimSelect<BPHRecoSelect>> select3(BPHSelectOperation::and_mode);
0229   select3.include(ptSel3);
0230   select3.include(etaSel3);
0231   BPHMultiSelect<BPHSlimSelect<BPHRecoSelect>> muoSel(BPHSelectOperation::or_mode);
0232   muoSel.include(select1);
0233   muoSel.include(select2);
0234   muoSel.include(select3);
0235 
0236   BPHMassSelect massPh(0.85, 1.20);
0237   BPHMassSelect massP1(2.95, 3.25);
0238   BPHMassSelect massP2(3.55, 3.80);
0239   BPHMassSelect massU1(9.10, 9.75);
0240   BPHMassSelect massU2(9.75, 10.20);
0241   BPHMassSelect massU3(10.20, 10.60);
0242 
0243   outF << "extract and dump JPsi" << endl;
0244   vector<BPHPlusMinusConstCandPtr> lJPsi = onia->getList(BPHOniaToMuMuBuilder::Psi1);
0245   int iJPsi;
0246   int nJPsi = lJPsi.size();
0247   outF << nJPsi << " JPsi cand found" << endl;
0248   for (iJPsi = 0; iJPsi < nJPsi; ++iJPsi)
0249     dumpRecoCand("JPsi", lJPsi[iJPsi].get());
0250 
0251   outF << "extract and dump specific onia" << endl;
0252   vector<BPHPlusMinusConstCandPtr> lPmm = onia->getList(BPHOniaToMuMuBuilder::Phi, &muoSel, &massPh);
0253   vector<BPHPlusMinusConstCandPtr> lPsi1 = onia->getList(BPHOniaToMuMuBuilder::Psi1, &muoSel, &massP1);
0254   vector<BPHPlusMinusConstCandPtr> lPsi2 = onia->getList(BPHOniaToMuMuBuilder::Psi2, &muoSel, &massP2);
0255   vector<BPHPlusMinusConstCandPtr> lUps1 = onia->getList(BPHOniaToMuMuBuilder::Ups1, &muoSel, &massU1);
0256   vector<BPHPlusMinusConstCandPtr> lUps2 = onia->getList(BPHOniaToMuMuBuilder::Ups2, &muoSel, &massU2);
0257   vector<BPHPlusMinusConstCandPtr> lUps3 = onia->getList(BPHOniaToMuMuBuilder::Ups3, &muoSel, &massU3);
0258   int iPhi;
0259   int nPhi = lPmm.size();
0260   outF << nPhi << " PhiMuMu cand found" << endl;
0261   for (iPhi = 0; iPhi < nPhi; ++iPhi)
0262     dumpRecoCand("PhiMuMu", lPmm[iPhi].get());
0263   int iPsi1;
0264   int nPsi1 = lPsi1.size();
0265   outF << nPsi1 << " Psi1 cand found" << endl;
0266   for (iPsi1 = 0; iPsi1 < nPsi1; ++iPsi1)
0267     dumpRecoCand("Psi1", lPsi1[iPsi1].get());
0268   int iPsi2;
0269   int nPsi2 = lPsi2.size();
0270   outF << nPsi2 << " Psi2 cand found" << endl;
0271   for (iPsi2 = 0; iPsi2 < nPsi2; ++iPsi2)
0272     dumpRecoCand("Psi2", lPsi2[iPsi2].get());
0273   int iUps1;
0274   int nUps1 = lUps1.size();
0275   outF << nUps1 << " Ups1 cand found" << endl;
0276   for (iUps1 = 0; iUps1 < nUps1; ++iUps1)
0277     dumpRecoCand("Ups1", lUps1[iUps1].get());
0278   int iUps2;
0279   int nUps2 = lUps2.size();
0280   outF << nUps2 << " Ups2 cand found" << endl;
0281   for (iUps2 = 0; iUps2 < nUps2; ++iUps2)
0282     dumpRecoCand("Ups2", lUps2[iUps2].get());
0283   int iUps3;
0284   int nUps3 = lUps3.size();
0285   outF << nUps3 << " Ups3 cand found" << endl;
0286   for (iUps3 = 0; iUps3 < nUps3; ++iUps3)
0287     dumpRecoCand("Ups3", lUps3[iUps3].get());
0288   delete onia;
0289 
0290   if (!nPsi1)
0291     return;
0292   if (!nrc)
0293     return;
0294 
0295   BPHParticlePtSelect tkPt(0.7);
0296   BPHMassSelect mJPsi(3.00, 3.17);
0297   BPHChi2Select chi2Bs(0.02);
0298 
0299   // build and dump Bu
0300 
0301   outF << "build and dump Bu" << endl;
0302   BPHBuToJPsiKBuilder* bu = nullptr;
0303   if (usePF)
0304     bu = new BPHBuToJPsiKBuilder(es,
0305                                  lJPsi,  //lFull,//lPsi1,
0306                                  BPHRecoBuilder::createCollection(pfCands));
0307   else if (usePC)
0308     bu = new BPHBuToJPsiKBuilder(es,
0309                                  lJPsi,  //lFull,//lPsi1,
0310                                  BPHRecoBuilder::createCollection(pcCands));
0311   else if (useGP)
0312     bu = new BPHBuToJPsiKBuilder(es,
0313                                  lJPsi,  //lFull,//lPsi1,
0314                                  BPHRecoBuilder::createCollection(gpCands));
0315 
0316   vector<BPHRecoConstCandPtr> lBu = bu->build();
0317 
0318   int iBu;
0319   int nBu = lBu.size();
0320   outF << nBu << " Bu cand found" << endl;
0321   for (iBu = 0; iBu < nBu; ++iBu)
0322     dumpRecoCand("Bu", lBu[iBu].get());
0323   // the following is an example of decay reconstruction starting from
0324   // specific reco::Candidates
0325   // here the final decay products are taken from already reconstructed B+,
0326   // so there's no physical sense in the operation
0327   for (iBu = 0; iBu < nBu; ++iBu) {
0328     const BPHRecoCandidate* bu = lBu[iBu].get();
0329     const reco::Candidate* mPos = bu->originalReco(bu->getDaug("JPsi/MuPos"));
0330     const reco::Candidate* mNeg = bu->originalReco(bu->getDaug("JPsi/MuNeg"));
0331     const reco::Candidate* kaon = bu->originalReco(bu->getDaug("Kaon"));
0332     BPHRecoCandidatePtr njp = BPHPlusMinusCandidateWrap::create(&es);
0333     njp->add("MuPos", mPos, BPHParticleMasses::muonMass, BPHParticleMasses::muonMSigma);
0334     njp->add("MuNeg", mNeg, BPHParticleMasses::muonMass, BPHParticleMasses::muonMSigma);
0335     BPHRecoCandidate nbu(&es);
0336     nbu.add("JPsi", njp);
0337     nbu.add("Kaon", kaon, BPHParticleMasses::kaonMass, BPHParticleMasses::kaonMSigma);
0338     nbu.kinematicTree("JPsi", BPHParticleMasses::jPsiMass, BPHParticleMasses::jPsiMWidth);
0339     dumpRecoCand("nBu", &nbu);
0340   }
0341 
0342   // build and dump Kx0
0343 
0344   BPHKx0ToKPiBuilder* kx0 = nullptr;
0345   if (usePF)
0346     kx0 = new BPHKx0ToKPiBuilder(
0347         es, BPHRecoBuilder::createCollection(pfCands), BPHRecoBuilder::createCollection(pfCands));
0348   else if (usePC)
0349     kx0 = new BPHKx0ToKPiBuilder(
0350         es, BPHRecoBuilder::createCollection(pcCands), BPHRecoBuilder::createCollection(pcCands));
0351   else if (useGP)
0352     kx0 = new BPHKx0ToKPiBuilder(
0353         es, BPHRecoBuilder::createCollection(gpCands), BPHRecoBuilder::createCollection(gpCands));
0354 
0355   vector<BPHPlusMinusConstCandPtr> lKx0 = kx0->build();
0356 
0357   int iKx0;
0358   int nKx0 = lKx0.size();
0359   outF << nKx0 << " Kx0 cand found" << endl;
0360   for (iKx0 = 0; iKx0 < nKx0; ++iKx0)
0361     dumpRecoCand("Kx0", lKx0[iKx0].get());
0362 
0363   delete kx0;
0364 
0365   // build and dump Bd
0366 
0367   outF << "build and dump Bd" << endl;
0368   if (nKx0) {
0369     BPHBdToJPsiKxBuilder* bd = new BPHBdToJPsiKxBuilder(es, lJPsi, lKx0);
0370     vector<BPHRecoConstCandPtr> lBd = bd->build();
0371     int iBd;
0372     int nBd = lBd.size();
0373     outF << nBd << " Bd cand found" << endl;
0374     for (iBd = 0; iBd < nBd; ++iBd)
0375       dumpRecoCand("Bd", lBd[iBd].get());
0376   }
0377 
0378   // build and dump Phi
0379 
0380   BPHPhiToKKBuilder* phi = nullptr;
0381   if (usePF)
0382     phi =
0383         new BPHPhiToKKBuilder(es, BPHRecoBuilder::createCollection(pfCands), BPHRecoBuilder::createCollection(pfCands));
0384   else if (usePC)
0385     phi =
0386         new BPHPhiToKKBuilder(es, BPHRecoBuilder::createCollection(pcCands), BPHRecoBuilder::createCollection(pcCands));
0387   else if (useGP)
0388     phi =
0389         new BPHPhiToKKBuilder(es, BPHRecoBuilder::createCollection(gpCands), BPHRecoBuilder::createCollection(gpCands));
0390 
0391   vector<BPHPlusMinusConstCandPtr> lPkk = phi->build();
0392 
0393   int iPkk;
0394   int nPkk = lPkk.size();
0395   outF << nPkk << " PhiKK cand found" << endl;
0396   for (iPkk = 0; iPkk < nPkk; ++iPkk)
0397     dumpRecoCand("PhiKK", lPkk[iPkk].get());
0398 
0399   delete phi;
0400 
0401   // build and dump Bs
0402 
0403   outF << "build and dump Bs" << endl;
0404   if (nPkk) {
0405     BPHBsToJPsiPhiBuilder* bs = new BPHBsToJPsiPhiBuilder(es, lJPsi, lPkk);
0406     vector<BPHRecoConstCandPtr> lBs = bs->build();
0407     int iBs;
0408     int nBs = lBs.size();
0409     outF << nBs << " Bs cand found" << endl;
0410     for (iBs = 0; iBs < nBs; ++iBs)
0411       dumpRecoCand("Bs", lBs[iBs].get());
0412   }
0413 
0414   return;
0415 }
0416 
0417 void TestBPHSpecificDecay::endJob() {
0418   *fPtr << "TestBPHSpecificDecay::endJob" << endl;
0419   TDirectory* currentDir = gDirectory;
0420   TFile file(outHist.c_str(), "RECREATE");
0421   map<string, TH1F*>::iterator iter = histoMap.begin();
0422   map<string, TH1F*>::iterator iend = histoMap.end();
0423   while (iter != iend)
0424     iter++->second->Write();
0425   currentDir->cd();
0426   return;
0427 }
0428 
0429 void TestBPHSpecificDecay::dumpRecoCand(const string& name, const BPHRecoCandidate* cand) {
0430   fillHisto(name, cand);
0431   if ((name == "PhiMuMu") || (name == "Psi1") || (name == "Psi2") || (name == "Ups1") || (name == "Ups2") ||
0432       (name == "Ups3"))
0433     fillHisto("Fsel", cand);
0434 
0435   ostream& outF = *fPtr;
0436 
0437   static string cType = " cowboy";
0438   static string sType = " sailor";
0439   static string dType = "";
0440   string* type;
0441   const BPHPlusMinusCandidate* pmCand = dynamic_cast<const BPHPlusMinusCandidate*>(cand);
0442   if (pmCand != nullptr) {
0443     if (pmCand->isCowboy())
0444       type = &cType;
0445     else
0446       type = &sType;
0447   } else
0448     type = &dType;
0449 
0450   outF << "****** " << name << "   cand mass: " << cand->composite().mass() << " momentum " << cand->composite().px()
0451        << " " << cand->composite().py() << " " << cand->composite().pz() << *type << endl;
0452 
0453   bool validFit = cand->isValidFit();
0454   const RefCountedKinematicTree kt = cand->kinematicTree();
0455   const RefCountedKinematicParticle kp = cand->currentParticle();
0456   if (validFit) {
0457     outF << "****** " << name << " constr mass: " << cand->p4().mass() << " momentum " << cand->p4().px() << " "
0458          << cand->p4().py() << " " << cand->p4().pz() << endl;
0459   }
0460 
0461   const reco::Vertex& vx = cand->vertex();
0462   const reco::Vertex::Point& vp = vx.position();
0463   double chi2 = vx.chi2();
0464   int ndof = lround(vx.ndof());
0465   double prob = TMath::Prob(chi2, ndof);
0466   string tdca = "";
0467   if (pmCand != nullptr) {
0468     stringstream sstr;
0469     sstr << " - " << pmCand->cAppInRPhi().distance();
0470     tdca = sstr.str();
0471   }
0472   outF << "****** " << name << " vertex: " << vx.isFake() << " " << vx.isValid() << " - " << chi2 << " " << ndof << " "
0473        << prob << " - " << vp.X() << " " << vp.Y() << " " << vp.Z() << tdca << endl;
0474 
0475   const vector<string>& dl = cand->daugNames();
0476   int i;
0477   int n = dl.size();
0478   for (i = 0; i < n; ++i) {
0479     const string& name = dl[i];
0480     const reco::Candidate* dp = cand->getDaug(name);
0481     GlobalPoint gp(vp.X(), vp.Y(), vp.Z());
0482     GlobalVector dm(0.0, 0.0, 0.0);
0483     const reco::TransientTrack* tt = cand->getTransientTrack(dp);
0484     if (tt != nullptr) {
0485       TrajectoryStateClosestToPoint tscp = tt->trajectoryStateClosestToPoint(gp);
0486       dm = tscp.momentum();
0487       //      TrajectoryStateOnSurface tsos = tt->stateOnSurface( gp );
0488       //      GlobalVector gv = tsos.globalMomentum();
0489     }
0490     outF << "daughter " << i << " " << name << " " << (dp->charge() > 0 ? '+' : '-') << " momentum: " << dp->px() << " "
0491          << dp->py() << " " << dp->pz() << " - at vertex: " << dm.x() << " " << dm.y() << " " << dm.z() << endl;
0492   }
0493   const vector<string>& dc = cand->compNames();
0494   int j;
0495   int m = dc.size();
0496   for (j = 0; j < m; ++j) {
0497     const string& name = dc[j];
0498     const BPHRecoCandidate* dp = cand->getComp(name).get();
0499     outF << "composite daughter " << j << " " << name << " momentum: " << dp->composite().px() << " "
0500          << dp->composite().py() << " " << dp->composite().pz() << endl;
0501   }
0502 
0503   if (validFit) {
0504     const RefCountedKinematicVertex kv = cand->currentDecayVertex();
0505     GlobalPoint gp = kv->position();
0506     outF << "   kin fit vertex: " << gp.x() << " " << gp.y() << " " << gp.z() << endl;
0507     vector<RefCountedKinematicParticle> dk = kt->finalStateParticles();
0508     int k;
0509     int l = dk.size();
0510     for (k = 0; k < l; ++k) {
0511       const reco::TransientTrack& tt = dk[k]->refittedTransientTrack();
0512       TrajectoryStateClosestToPoint tscp = tt.trajectoryStateClosestToPoint(gp);
0513       GlobalVector dm = tscp.momentum();
0514       //    TrajectoryStateOnSurface tsos = tt.stateOnSurface( gp );
0515       //    GlobalVector dm = tsos.globalMomentum();
0516       outF << "daughter " << k << " refitted: " << dm.x() << " " << dm.y() << " " << dm.z() << endl;
0517     }
0518   }
0519 
0520   return;
0521 }
0522 
0523 void TestBPHSpecificDecay::fillHisto(const string& name, const BPHRecoCandidate* cand) {
0524   string mass = "mass";
0525   string mcst = "mcst";
0526   fillHisto(mass + name, cand->composite().mass());
0527   if (cand->isValidFit())
0528     fillHisto(mcst + name, cand->p4().mass());
0529 
0530   const vector<string>& dc = cand->compNames();
0531   int i;
0532   int n = dc.size();
0533   for (i = 0; i < n; ++i) {
0534     const string& daug = dc[i];
0535     const BPHRecoCandidate* dptr = cand->getComp(daug).get();
0536     fillHisto(mass + name + daug, dptr->composite().mass());
0537   }
0538   return;
0539 }
0540 
0541 void TestBPHSpecificDecay::fillHisto(const string& name, float x) {
0542   map<string, TH1F*>::iterator iter = histoMap.find(name);
0543   map<string, TH1F*>::iterator iend = histoMap.end();
0544   if (iter == iend)
0545     return;
0546   iter->second->Fill(x);
0547   return;
0548 }
0549 
0550 void TestBPHSpecificDecay::createHisto(const string& name, int nbin, float hmin, float hmax) {
0551   histoMap[name] = new TH1F(name.c_str(), name.c_str(), nbin, hmin, hmax);
0552   return;
0553 }
0554 
0555 #include "FWCore/Framework/interface/MakerMacros.h"
0556 
0557 DEFINE_FWK_MODULE(TestBPHSpecificDecay);