Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:15:38 

0001 #include "HeavyFlavorAnalysis/SpecificDecay/test/stubs/TestBPHSpecificDecay.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/BPHBsToJPsiPhiBuilder.h"
0017 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHBdToJPsiKxBuilder.h"
0018 #include "HeavyFlavorAnalysis/SpecificDecay/interface/BPHParticleMasses.h"
0019 
0020 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHRecoBuilder.h"
0021 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHRecoSelect.h"
0022 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHRecoCandidate.h"
0023 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHPlusMinusCandidate.h"
0024 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHMomentumSelect.h"
0025 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHVertexSelect.h"
0026 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHTrackReference.h"
0027 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHMultiSelect.h"
0028 
0029 #include "DataFormats/TrackReco/interface/Track.h"
0030 
0031 #include "DataFormats/PatCandidates/interface/Muon.h"
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 void TestBPHSpecificDecay::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0077   edm::ParameterSetDescription desc;
0078   desc.add<string>("patMuonLabel", "");
0079   desc.add<string>("ccCandsLabel", "");
0080   desc.add<string>("pfCandsLabel", "");
0081   desc.add<string>("pcCandsLabel", "");
0082   desc.add<string>("gpCandsLabel", "");
0083   desc.add<string>("outDump", "dump.txt");
0084   desc.add<string>("outHist", "hist.root");
0085   descriptions.add("testBPHSpecificDecay", desc);
0086   return;
0087 }
0088 
0089 void TestBPHSpecificDecay::beginJob() {
0090   *fPtr << "TestBPHSpecificDecay::beginJob" << endl;
0091   createHisto("massJPsi", 60, 2.95, 3.25);   // JPsi mass
0092   createHisto("mcstJPsi", 60, 2.95, 3.25);   // JPsi mass, with constraint
0093   createHisto("massKx0", 50, 0.80, 1.05);    // Kx0  mass
0094   createHisto("massPhi", 40, 1.01, 1.03);    // Phi  mass
0095   createHisto("massBu", 20, 5.00, 5.50);     // Bu   mass
0096   createHisto("mcstBu", 20, 5.00, 5.50);     // Bu   mass, with constraint
0097   createHisto("massBd", 20, 5.00, 5.50);     // Bd   mass
0098   createHisto("mcstBd", 20, 5.00, 5.50);     // Bd   mass, with constraint
0099   createHisto("massBs", 20, 5.10, 5.60);     // Bs   mass
0100   createHisto("mcstBs", 20, 5.10, 5.60);     // Bs   mass, with constraint
0101   createHisto("massBsPhi", 50, 1.01, 1.03);  // Phi  mass in Bs decay
0102   createHisto("massBdKx0", 50, 0.80, 1.05);  // Kx0  mass in Bd decay
0103 
0104   createHisto("massFull", 200, 2.00, 12.00);  // Full mass
0105   createHisto("massFsel", 200, 2.00, 12.00);  // Full mass
0106   createHisto("massPhi", 70, 0.85, 1.20);     // Psi1 mass
0107   createHisto("massPsi1", 60, 2.95, 3.25);    // Psi1 mass
0108   createHisto("massPsi2", 50, 3.55, 3.80);    // Psi2 mass
0109   createHisto("massUps1", 130, 9.10, 9.75);   // Ups1 mass
0110   createHisto("massUps2", 90, 9.75, 10.20);   // Ups2 mass
0111   createHisto("massUps3", 80, 10.20, 10.60);  // Ups3 mass
0112 
0113   return;
0114 }
0115 
0116 void TestBPHSpecificDecay::analyze(const edm::Event& ev, const edm::EventSetup& es) {
0117   ostream& outF = *fPtr;
0118   outF << "--------- event " << ev.id().run() << " / " << ev.id().event() << " ---------" << endl;
0119 
0120   BPHEventSetupWrapper ew(es, BPHRecoCandidate::transientTrackBuilder, &ttBToken);
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         ew, BPHRecoBuilder::createCollection(patMuon, "cfmig"), BPHRecoBuilder::createCollection(patMuon, "cfmig"));
0203   else if (useCC)
0204     onia = new BPHOniaToMuMuBuilder(
0205         ew, 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(ew, lJPsi, BPHRecoBuilder::createCollection(pfCands));
0305   else if (usePC)
0306     bu = new BPHBuToJPsiKBuilder(ew, lJPsi, BPHRecoBuilder::createCollection(pcCands));
0307   else if (useGP)
0308     bu = new BPHBuToJPsiKBuilder(ew, lJPsi, BPHRecoBuilder::createCollection(gpCands));
0309 
0310   vector<BPHRecoConstCandPtr> lBu = bu->build();
0311 
0312   int iBu;
0313   int nBu = lBu.size();
0314   outF << nBu << " Bu cand found" << endl;
0315   for (iBu = 0; iBu < nBu; ++iBu)
0316     dumpRecoCand("Bu", lBu[iBu].get());
0317   // the following is an example of decay reconstruction starting from
0318   // specific reco::Candidates
0319   // here the final decay products are taken from already reconstructed B+,
0320   // so there's no physical sense in the operation
0321   for (iBu = 0; iBu < nBu; ++iBu) {
0322     const BPHRecoCandidate* bu = lBu[iBu].get();
0323     const reco::Candidate* mPos = bu->originalReco(bu->getDaug("JPsi/MuPos"));
0324     const reco::Candidate* mNeg = bu->originalReco(bu->getDaug("JPsi/MuNeg"));
0325     const reco::Candidate* kaon = bu->originalReco(bu->getDaug("Kaon"));
0326     BPHRecoCandidatePtr njp = BPHPlusMinusCandidateWrap::create(&ew);
0327     njp->add("MuPos", mPos, BPHParticleMasses::muonMass, BPHParticleMasses::muonMSigma);
0328     njp->add("MuNeg", mNeg, BPHParticleMasses::muonMass, BPHParticleMasses::muonMSigma);
0329     BPHRecoCandidate nbu(&ew);
0330     nbu.add("JPsi", njp);
0331     nbu.add("Kaon", kaon, BPHParticleMasses::kaonMass, BPHParticleMasses::kaonMSigma);
0332     nbu.kinematicTree("JPsi", BPHParticleMasses::jPsiMass, BPHParticleMasses::jPsiMWidth);
0333     dumpRecoCand("nBu", &nbu);
0334   }
0335 
0336   // build and dump Kx0
0337 
0338   BPHKx0ToKPiBuilder* kx0 = nullptr;
0339   if (usePF)
0340     kx0 = new BPHKx0ToKPiBuilder(
0341         ew, BPHRecoBuilder::createCollection(pfCands), BPHRecoBuilder::createCollection(pfCands));
0342   else if (usePC)
0343     kx0 = new BPHKx0ToKPiBuilder(
0344         ew, BPHRecoBuilder::createCollection(pcCands), BPHRecoBuilder::createCollection(pcCands));
0345   else if (useGP)
0346     kx0 = new BPHKx0ToKPiBuilder(
0347         ew, BPHRecoBuilder::createCollection(gpCands), BPHRecoBuilder::createCollection(gpCands));
0348 
0349   vector<BPHPlusMinusConstCandPtr> lKx0 = kx0->build();
0350 
0351   int iKx0;
0352   int nKx0 = lKx0.size();
0353   outF << nKx0 << " Kx0 cand found" << endl;
0354   for (iKx0 = 0; iKx0 < nKx0; ++iKx0)
0355     dumpRecoCand("Kx0", lKx0[iKx0].get());
0356 
0357   delete kx0;
0358 
0359   // build and dump Bd
0360 
0361   outF << "build and dump Bd" << endl;
0362   if (nKx0) {
0363     BPHBdToJPsiKxBuilder* bd = new BPHBdToJPsiKxBuilder(ew, lJPsi, lKx0);
0364     vector<BPHRecoConstCandPtr> lBd = bd->build();
0365     int iBd;
0366     int nBd = lBd.size();
0367     outF << nBd << " Bd cand found" << endl;
0368     for (iBd = 0; iBd < nBd; ++iBd)
0369       dumpRecoCand("Bd", lBd[iBd].get());
0370   }
0371 
0372   // build and dump Phi
0373 
0374   BPHPhiToKKBuilder* phi = nullptr;
0375   if (usePF)
0376     phi =
0377         new BPHPhiToKKBuilder(ew, BPHRecoBuilder::createCollection(pfCands), BPHRecoBuilder::createCollection(pfCands));
0378   else if (usePC)
0379     phi =
0380         new BPHPhiToKKBuilder(ew, BPHRecoBuilder::createCollection(pcCands), BPHRecoBuilder::createCollection(pcCands));
0381   else if (useGP)
0382     phi =
0383         new BPHPhiToKKBuilder(ew, BPHRecoBuilder::createCollection(gpCands), BPHRecoBuilder::createCollection(gpCands));
0384 
0385   vector<BPHPlusMinusConstCandPtr> lPkk = phi->build();
0386 
0387   int iPkk;
0388   int nPkk = lPkk.size();
0389   outF << nPkk << " PhiKK cand found" << endl;
0390   for (iPkk = 0; iPkk < nPkk; ++iPkk)
0391     dumpRecoCand("PhiKK", lPkk[iPkk].get());
0392 
0393   delete phi;
0394 
0395   // build and dump Bs
0396 
0397   outF << "build and dump Bs" << endl;
0398   if (nPkk) {
0399     BPHBsToJPsiPhiBuilder* bs = new BPHBsToJPsiPhiBuilder(ew, lJPsi, lPkk);
0400     vector<BPHRecoConstCandPtr> lBs = bs->build();
0401     int iBs;
0402     int nBs = lBs.size();
0403     outF << nBs << " Bs cand found" << endl;
0404     for (iBs = 0; iBs < nBs; ++iBs)
0405       dumpRecoCand("Bs", lBs[iBs].get());
0406   }
0407 
0408   return;
0409 }
0410 
0411 void TestBPHSpecificDecay::endJob() {
0412   *fPtr << "TestBPHSpecificDecay::endJob" << endl;
0413   TDirectory* currentDir = gDirectory;
0414   TFile file(outHist.c_str(), "RECREATE");
0415   map<string, TH1F*>::iterator iter = histoMap.begin();
0416   map<string, TH1F*>::iterator iend = histoMap.end();
0417   while (iter != iend)
0418     iter++->second->Write();
0419   currentDir->cd();
0420   return;
0421 }
0422 
0423 void TestBPHSpecificDecay::dumpRecoCand(const string& name, const BPHRecoCandidate* cand) {
0424   fillHisto(name, cand);
0425   if ((name == "PhiMuMu") || (name == "Psi1") || (name == "Psi2") || (name == "Ups1") || (name == "Ups2") ||
0426       (name == "Ups3"))
0427     fillHisto("Fsel", cand);
0428 
0429   ostream& outF = *fPtr;
0430 
0431   static string cType = " cowboy";
0432   static string sType = " sailor";
0433   static string dType = "";
0434   string* type;
0435   const BPHPlusMinusCandidate* pmCand = dynamic_cast<const BPHPlusMinusCandidate*>(cand);
0436   if (pmCand != nullptr) {
0437     if (pmCand->isCowboy())
0438       type = &cType;
0439     else
0440       type = &sType;
0441   } else
0442     type = &dType;
0443 
0444   outF << "****** " << name << "   cand mass: " << cand->composite().mass() << " momentum " << cand->composite().px()
0445        << " " << cand->composite().py() << " " << cand->composite().pz() << *type << endl;
0446 
0447   bool validFit = cand->isValidFit();
0448   const RefCountedKinematicTree kt = cand->kinematicTree();
0449   const RefCountedKinematicParticle kp = cand->currentParticle();
0450   if (validFit) {
0451     outF << "****** " << name << " constr mass: " << cand->p4().mass() << " momentum " << cand->p4().px() << " "
0452          << cand->p4().py() << " " << cand->p4().pz() << endl;
0453   }
0454 
0455   const reco::Vertex& vx = cand->vertex();
0456   const reco::Vertex::Point& vp = vx.position();
0457   double chi2 = vx.chi2();
0458   int ndof = lround(vx.ndof());
0459   double prob = TMath::Prob(chi2, ndof);
0460   string tdca = "";
0461   if (pmCand != nullptr) {
0462     stringstream sstr;
0463     sstr << " - " << pmCand->cAppInRPhi().distance();
0464     tdca = sstr.str();
0465   }
0466   outF << "****** " << name << " vertex: " << vx.isFake() << " " << vx.isValid() << " - " << chi2 << " " << ndof << " "
0467        << prob << " - " << vp.X() << " " << vp.Y() << " " << vp.Z() << tdca << endl;
0468 
0469   const vector<string>& dl = cand->daugNames();
0470   int i;
0471   int n = dl.size();
0472   for (i = 0; i < n; ++i) {
0473     const string& name = dl[i];
0474     const reco::Candidate* dp = cand->getDaug(name);
0475     GlobalPoint gp(vp.X(), vp.Y(), vp.Z());
0476     GlobalVector dm(0.0, 0.0, 0.0);
0477     const reco::TransientTrack* tt = cand->getTransientTrack(dp);
0478     if (tt != nullptr) {
0479       TrajectoryStateClosestToPoint tscp = tt->trajectoryStateClosestToPoint(gp);
0480       dm = tscp.momentum();
0481       //      TrajectoryStateOnSurface tsos = tt->stateOnSurface( gp );
0482       //      GlobalVector gv = tsos.globalMomentum();
0483     }
0484     outF << "daughter " << i << " " << name << " " << (dp->charge() > 0 ? '+' : '-') << " momentum: " << dp->px() << " "
0485          << dp->py() << " " << dp->pz() << " - at vertex: " << dm.x() << " " << dm.y() << " " << dm.z() << endl;
0486   }
0487   const vector<string>& dc = cand->compNames();
0488   int j;
0489   int m = dc.size();
0490   for (j = 0; j < m; ++j) {
0491     const string& name = dc[j];
0492     const BPHRecoCandidate* dp = cand->getComp(name).get();
0493     outF << "composite daughter " << j << " " << name << " momentum: " << dp->composite().px() << " "
0494          << dp->composite().py() << " " << dp->composite().pz() << endl;
0495   }
0496 
0497   if (validFit) {
0498     const RefCountedKinematicVertex kv = cand->currentDecayVertex();
0499     GlobalPoint gp = kv->position();
0500     outF << "   kin fit vertex: " << gp.x() << " " << gp.y() << " " << gp.z() << endl;
0501     vector<RefCountedKinematicParticle> dk = kt->finalStateParticles();
0502     int k;
0503     int l = dk.size();
0504     for (k = 0; k < l; ++k) {
0505       const reco::TransientTrack& tt = dk[k]->refittedTransientTrack();
0506       TrajectoryStateClosestToPoint tscp = tt.trajectoryStateClosestToPoint(gp);
0507       GlobalVector dm = tscp.momentum();
0508       //    TrajectoryStateOnSurface tsos = tt.stateOnSurface( gp );
0509       //    GlobalVector dm = tsos.globalMomentum();
0510       outF << "daughter " << k << " refitted: " << dm.x() << " " << dm.y() << " " << dm.z() << endl;
0511     }
0512   }
0513 
0514   return;
0515 }
0516 
0517 void TestBPHSpecificDecay::fillHisto(const string& name, const BPHRecoCandidate* cand) {
0518   string mass = "mass";
0519   string mcst = "mcst";
0520   fillHisto(mass + name, cand->composite().mass());
0521   if (cand->isValidFit())
0522     fillHisto(mcst + name, cand->p4().mass());
0523 
0524   const vector<string>& dc = cand->compNames();
0525   int i;
0526   int n = dc.size();
0527   for (i = 0; i < n; ++i) {
0528     const string& daug = dc[i];
0529     const BPHRecoCandidate* dptr = cand->getComp(daug).get();
0530     fillHisto(mass + name + daug, dptr->composite().mass());
0531   }
0532   return;
0533 }
0534 
0535 void TestBPHSpecificDecay::fillHisto(const string& name, float x) {
0536   map<string, TH1F*>::iterator iter = histoMap.find(name);
0537   map<string, TH1F*>::iterator iend = histoMap.end();
0538   if (iter == iend)
0539     return;
0540   iter->second->Fill(x);
0541   return;
0542 }
0543 
0544 void TestBPHSpecificDecay::createHisto(const string& name, int nbin, float hmin, float hmax) {
0545   histoMap[name] = new TH1F(name.c_str(), name.c_str(), nbin, hmin, hmax);
0546   return;
0547 }
0548 
0549 #include "FWCore/Framework/interface/MakerMacros.h"
0550 
0551 DEFINE_FWK_MODULE(TestBPHSpecificDecay);
[ 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!