Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 #include "HeavyFlavorAnalysis/RecoDecay/test/stubs/TestBPHRecoDecay.h"
0002 
0003 #include "FWCore/Framework/interface/MakerMacros.h"
0004 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0006 
0007 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHRecoBuilder.h"
0008 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHRecoSelect.h"
0009 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHRecoCandidate.h"
0010 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHPlusMinusCandidate.h"
0011 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHMomentumSelect.h"
0012 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHVertexSelect.h"
0013 #include "HeavyFlavorAnalysis/RecoDecay/interface/BPHTrackReference.h"
0014 
0015 #include "DataFormats/PatCandidates/interface/Muon.h"
0016 #include "DataFormats/TrackReco/interface/Track.h"
0017 
0018 #include "DataFormats/PatCandidates/interface/GenericParticle.h"
0019 #include "DataFormats/PatCandidates/interface/CompositeCandidate.h"
0020 
0021 #include "RecoVertex/KinematicFit/interface/TwoTrackMassKinematicConstraint.h"
0022 #include <TH1.h>
0023 #include <TFile.h>
0024 #include <TMath.h>
0025 
0026 #include <set>
0027 #include <string>
0028 #include <iostream>
0029 #include <fstream>
0030 
0031 using namespace std;
0032 
0033 #define SET_LABEL(NAME, PSET) (NAME = PSET.getParameter<string>(#NAME))
0034 // SET_LABEL(xyz,ps);
0035 // is equivalent to
0036 // xyz = ps.getParameter<string>( "xyx" )
0037 
0038 TestBPHRecoDecay::TestBPHRecoDecay(const edm::ParameterSet& ps) {
0039   usePM = (!SET_LABEL(patMuonLabel, ps).empty());
0040   useCC = (!SET_LABEL(ccCandsLabel, ps).empty());
0041   usePF = (!SET_LABEL(pfCandsLabel, ps).empty());
0042   usePC = (!SET_LABEL(pcCandsLabel, ps).empty());
0043   useGP = (!SET_LABEL(gpCandsLabel, ps).empty());
0044 
0045   if (usePM)
0046     consume<pat::MuonCollection>(patMuonToken, patMuonLabel);
0047   if (useCC)
0048     consume<vector<pat::CompositeCandidate>>(ccCandsToken, ccCandsLabel);
0049   if (usePF)
0050     consume<vector<reco::PFCandidate>>(pfCandsToken, pfCandsLabel);
0051   if (usePC)
0052     consume<vector<BPHTrackReference::candidate>>(pcCandsToken, pcCandsLabel);
0053   if (useGP)
0054     consume<vector<pat::GenericParticle>>(gpCandsToken, gpCandsLabel);
0055   SET_LABEL(outDump, ps);
0056   SET_LABEL(outHist, ps);
0057   if (outDump.empty())
0058     fPtr = &cout;
0059   else
0060     fPtr = new ofstream(outDump.c_str());
0061 }
0062 
0063 TestBPHRecoDecay::~TestBPHRecoDecay() {}
0064 
0065 void TestBPHRecoDecay::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0066   edm::ParameterSetDescription desc;
0067   desc.add<string>("patMuonLabel", "");
0068   desc.add<string>("ccCandsLabel", "");
0069   desc.add<string>("pfCandsLabel", "");
0070   desc.add<string>("pcCandsLabel", "");
0071   desc.add<string>("gpCandsLabel", "");
0072   desc.add<string>("outDump", "dump.txt");
0073   desc.add<string>("outHist", "hist.root");
0074   descriptions.add("testBPHRecoDecay", desc);
0075   return;
0076 }
0077 
0078 void TestBPHRecoDecay::beginJob() {
0079   *fPtr << "TestBPHRecoDecay::beginJob" << endl;
0080   createHisto("massJPsi", 50, 2.85, 3.35);     // JPsi mass
0081   createHisto("mcstJPsi", 50, 2.85, 3.35);     // JPsi mass, with constraint
0082   createHisto("massPhi", 50, 0.995, 1.045);    // Phi  mass
0083   createHisto("massBu", 20, 5.0, 5.5);         // Bu   mass
0084   createHisto("mcstBu", 20, 5.0, 5.5);         // Bu   mass, with constraint
0085   createHisto("massBs", 20, 5.1, 5.6);         // Bs   mass
0086   createHisto("mcstBs", 20, 5.1, 5.6);         // Bs   mass, with constraint
0087   createHisto("massBsPhi", 50, 0.995, 1.045);  // Phi  mass in Bs decay
0088   return;
0089 }
0090 
0091 void TestBPHRecoDecay::analyze(const edm::Event& ev, const edm::EventSetup& es) {
0092   ostream& outF = *fPtr;
0093   outF << "--------- event " << ev.id().run() << " / " << ev.id().event() << " ---------" << endl;
0094 
0095   // get object collections
0096   // collections are got through "BPHTokenWrapper" interface to allow
0097   // uniform access in different CMSSW versions
0098 
0099   int nrc = 0;
0100 
0101   // get reco::PFCandidate collection (in full AOD )
0102   edm::Handle<vector<reco::PFCandidate>> pfCands;
0103   if (usePF) {
0104     pfCandsToken.get(ev, pfCands);
0105     nrc = pfCands->size();
0106     if (pfCands.isValid())
0107       outF << nrc << " pfCands found" << endl;
0108     else
0109       outF << "no pfCands" << endl;
0110   }
0111 
0112   // get pat::PackedCandidate collection (in MiniAOD)
0113   // pat::PackedCandidate is not defined in CMSSW_5XY, so a
0114   // typedef (BPHTrackReference::candidate) is used, actually referring
0115   // to pat::PackedCandidate only for CMSSW versions where it's defined
0116   edm::Handle<vector<BPHTrackReference::candidate>> pcCands;
0117   if (usePC) {
0118     pcCandsToken.get(ev, pcCands);
0119     nrc = pcCands->size();
0120     if (pcCands.isValid())
0121       outF << nrc << " pcCands found" << endl;
0122     else
0123       outF << "no pcCands" << endl;
0124   }
0125 
0126   // get pat::GenericParticle collection (in skimmed data)
0127   edm::Handle<vector<pat::GenericParticle>> gpCands;
0128   if (useGP) {
0129     gpCandsToken.get(ev, gpCands);
0130     nrc = gpCands->size();
0131     if (gpCands.isValid())
0132       outF << nrc << " gpCands found" << endl;
0133     else
0134       outF << "no gpCands" << endl;
0135   }
0136 
0137   // get pat::Muon collection (in full AOD and MiniAOD)
0138   edm::Handle<pat::MuonCollection> patMuon;
0139   if (usePM) {
0140     patMuonToken.get(ev, patMuon);
0141     int n = patMuon->size();
0142     if (patMuon.isValid())
0143       outF << n << " muons found" << endl;
0144     else
0145       outF << "no muons" << endl;
0146   }
0147 
0148   // get muons from pat::CompositeCandidate objects describing onia;
0149   // muons from all composite objects are copied to an unique std::vector
0150   vector<const reco::Candidate*> muDaugs;
0151   set<const pat::Muon*> muonSet;
0152   if (useCC) {
0153     edm::Handle<vector<pat::CompositeCandidate>> ccCands;
0154     ccCandsToken.get(ev, ccCands);
0155     int n = ccCands->size();
0156     if (ccCands.isValid())
0157       outF << n << " ccCands found" << endl;
0158     else
0159       outF << "no ccCands" << endl;
0160     muDaugs.clear();
0161     muDaugs.reserve(n);
0162     muonSet.clear();
0163     set<const pat::Muon*>::const_iterator iter;
0164     set<const pat::Muon*>::const_iterator iend;
0165     int i;
0166     for (i = 0; i < n; ++i) {
0167       const pat::CompositeCandidate& cc = ccCands->at(i);
0168       int j;
0169       int m = cc.numberOfDaughters();
0170       for (j = 0; j < m; ++j) {
0171         const reco::Candidate* dp = cc.daughter(j);
0172         const pat::Muon* mp = dynamic_cast<const pat::Muon*>(dp);
0173         iter = muonSet.begin();
0174         iend = muonSet.end();
0175         bool add = (mp != nullptr) && (muonSet.find(mp) == iend);
0176         while (add && (iter != iend)) {
0177           if (BPHRecoBuilder::sameTrack(mp, *iter++, 1.0e-5))
0178             add = false;
0179         }
0180         if (add)
0181           muonSet.insert(mp);
0182       }
0183     }
0184     iter = muonSet.begin();
0185     iend = muonSet.end();
0186     while (iter != iend)
0187       muDaugs.push_back(*iter++);
0188     if ((n = muDaugs.size()))
0189       outF << n << " muons found" << endl;
0190   }
0191 
0192   //
0193   // starting objects selection
0194   //
0195 
0196   // muon selection by charge
0197   class MuonChargeSelect : public BPHRecoSelect {
0198   public:
0199     MuonChargeSelect(int c) : charge(c) {}
0200     ~MuonChargeSelect() override {}
0201     bool accept(const reco::Candidate& cand) const override {
0202       const pat::Muon* p = dynamic_cast<const pat::Muon*>(&cand);
0203       if (p == nullptr)
0204         return false;
0205       return ((charge * cand.charge()) > 0);
0206     }
0207 
0208   private:
0209     int charge;
0210   };
0211 
0212   // muon selection by Pt
0213   class MuonPtSelect : public BPHRecoSelect {
0214   public:
0215     MuonPtSelect(float pt) : ptCut(pt) {}
0216     ~MuonPtSelect() override {}
0217     bool accept(const reco::Candidate& cand) const override {
0218       const pat::Muon* p = dynamic_cast<const pat::Muon*>(&cand);
0219       if (p == nullptr)
0220         return false;
0221       return (p->p4().pt() > ptCut);
0222     }
0223 
0224   private:
0225     float ptCut;
0226   };
0227 
0228   // muon selection by eta
0229   class MuonEtaSelect : public BPHRecoSelect {
0230   public:
0231     MuonEtaSelect(float eta) : etaCut(eta) {}
0232     ~MuonEtaSelect() override {}
0233     bool accept(const reco::Candidate& cand) const override {
0234       const pat::Muon* p = dynamic_cast<const pat::Muon*>(&cand);
0235       if (p == nullptr)
0236         return false;
0237       return (fabs(p->p4().eta()) < etaCut);
0238     }
0239 
0240   private:
0241     float etaCut;
0242   };
0243 
0244   // kaon selection by charge
0245   class KaonChargeSelect : public BPHRecoSelect {
0246   public:
0247     KaonChargeSelect(int c) : charge(c) {}
0248     ~KaonChargeSelect() override {}
0249     bool accept(const reco::Candidate& cand) const override { return ((charge * cand.charge()) > 0); }
0250 
0251   private:
0252     int charge;
0253   };
0254 
0255   class KaonNeutralVeto : public BPHRecoSelect {
0256   public:
0257     KaonNeutralVeto() {}
0258     ~KaonNeutralVeto() override {}
0259     bool accept(const reco::Candidate& cand) const override { return lround(fabs(cand.charge())); }
0260   };
0261 
0262   // kaon selection by Pt
0263   class KaonPtSelect : public BPHRecoSelect {
0264   public:
0265     KaonPtSelect(float pt) : ptCut(pt) {}
0266     ~KaonPtSelect() override {}
0267     bool accept(const reco::Candidate& cand) const override { return (cand.p4().pt() > ptCut); }
0268 
0269   private:
0270     float ptCut;
0271   };
0272 
0273   // kaon selection by eta
0274   class KaonEtaSelect : public BPHRecoSelect {
0275   public:
0276     KaonEtaSelect(float eta) : etaCut(eta) {}
0277     ~KaonEtaSelect() override {}
0278     bool accept(const reco::Candidate& cand) const override { return (fabs(cand.p4().eta()) < etaCut); }
0279 
0280   private:
0281     float etaCut;
0282   };
0283 
0284   //
0285   // reconstructed object selection
0286   //
0287 
0288   // selection by mass
0289   class MassSelect : public BPHMomentumSelect {
0290   public:
0291     MassSelect(double minMass, double maxMass) : mMin(minMass), mMax(maxMass) {}
0292     ~MassSelect() override {}
0293     bool accept(const BPHDecayMomentum& cand) const override {
0294       double mass = cand.composite().mass();
0295       return ((mass > mMin) && (mass < mMax));
0296     }
0297 
0298   private:
0299     double mMin;
0300     double mMax;
0301   };
0302 
0303   // selection by chi^2
0304   class Chi2Select : public BPHVertexSelect {
0305   public:
0306     Chi2Select(double minProb) : mProb(minProb) {}
0307     ~Chi2Select() override {}
0308     bool accept(const BPHDecayVertex& cand) const override {
0309       const reco::Vertex& v = cand.vertex();
0310       if (v.isFake())
0311         return false;
0312       if (!v.isValid())
0313         return false;
0314       return (TMath::Prob(v.chi2(), lround(v.ndof())) > mProb);
0315     }
0316 
0317   private:
0318     double mProb;
0319   };
0320 
0321   // build and dump JPsi
0322 
0323   outF << "build and dump JPsi" << endl;
0324   MuonPtSelect muPt(4.0);
0325   MuonEtaSelect muEta(2.1);
0326   string muPos = "MuPos";
0327   string muNeg = "MuNeg";
0328   BPHRecoBuilder bJPsi(es);
0329   if (usePM) {
0330     bJPsi.add(muPos, BPHRecoBuilder::createCollection(patMuon, "cfmig"), 0.105658);
0331     bJPsi.add(muNeg, BPHRecoBuilder::createCollection(patMuon, "cfmig"), 0.105658);
0332   } else if (useCC) {
0333     bJPsi.add(muPos, BPHRecoBuilder::createCollection(muDaugs, "cfmig"), 0.105658);
0334     bJPsi.add(muNeg, BPHRecoBuilder::createCollection(muDaugs, "cfmig"), 0.105658);
0335   }
0336   bJPsi.filter(muPos, muPt);
0337   bJPsi.filter(muNeg, muPt);
0338   bJPsi.filter(muPos, muEta);
0339   bJPsi.filter(muNeg, muEta);
0340 
0341   MassSelect massJPsi(2.5, 3.7);
0342   Chi2Select chi2Valid(0.0);
0343   bJPsi.filter(massJPsi);
0344   bJPsi.filter(chi2Valid);
0345   vector<BPHPlusMinusConstCandPtr> lJPsi = BPHPlusMinusCandidate::build(bJPsi, muPos, muNeg, 3.096916, 0.00004);
0346   //  //  BPHPlusMinusCandidate::build function has embedded charge selection
0347   //  //  alternatively simple BPHRecoCandidate::build function can be used
0348   //  //  as in the following
0349   //  MuonChargeSelect mqPos( +1 );
0350   //  MuonChargeSelect mqNeg( +1 );
0351   //  bJPsi.filter( nPos, mqPos );
0352   //  bJPsi.filter( nNeg, mqNeg );
0353   //  vector<BPHRecoConstCandPtr> lJPsi = BPHRecoCandidate::build( bJPsi,
0354   //                                      3.096916, 0.00004 );
0355   int iJPsi;
0356   int nJPsi = lJPsi.size();
0357   outF << nJPsi << " JPsi cand found" << endl;
0358   for (iJPsi = 0; iJPsi < nJPsi; ++iJPsi)
0359     dumpRecoCand("JPsi", lJPsi[iJPsi].get());
0360 
0361   // build and dump Phi
0362 
0363   outF << "build and dump Phi" << endl;
0364   BPHRecoBuilder bPhi(es);
0365   KaonChargeSelect tkPos(+1);
0366   KaonChargeSelect tkNeg(-1);
0367   KaonPtSelect tkPt(0.7);
0368   string kPos = "KPos";
0369   string kNeg = "KNeg";
0370   if (usePF) {
0371     bPhi.add(kPos, BPHRecoBuilder::createCollection(pfCands), 0.493677);
0372     bPhi.add(kNeg, BPHRecoBuilder::createCollection(pfCands), 0.493677);
0373   } else if (usePC) {
0374     bPhi.add(kPos, BPHRecoBuilder::createCollection(pcCands), 0.493677);
0375     bPhi.add(kNeg, BPHRecoBuilder::createCollection(pcCands), 0.493677);
0376   } else if (useGP) {
0377     bPhi.add(kPos, BPHRecoBuilder::createCollection(gpCands), 0.493677);
0378     bPhi.add(kNeg, BPHRecoBuilder::createCollection(gpCands), 0.493677);
0379   }
0380   bPhi.filter(kPos, tkPos);
0381   bPhi.filter(kNeg, tkNeg);
0382   bPhi.filter(kPos, tkPt);
0383   bPhi.filter(kNeg, tkPt);
0384 
0385   MassSelect massPhi(1.00, 1.04);
0386   bPhi.filter(massPhi);
0387   vector<BPHRecoConstCandPtr> lPhi = BPHRecoCandidate::build(bPhi);
0388   //  //  BPHRecoCandidate::build function requires explicit charge selection
0389   //  //  alternatively BPHPlusMinusCandidate::build function can be used
0390   //  //  as in the following
0391   //  //  (filter functions with tkPos and tkNeg can be dropped)
0392   //  vector<const BPHPlusMinusConstCandPtr> lPhi =
0393   //               BPHPlusMinusCandidate::build( bPhi, kPos, kNeg );
0394   int iPhi;
0395   int nPhi = lPhi.size();
0396   outF << nPhi << " Phi cand found" << endl;
0397   for (iPhi = 0; iPhi < nPhi; ++iPhi)
0398     dumpRecoCand("Phi", lPhi[iPhi].get());
0399 
0400   // build and dump Bs
0401 
0402   if (nJPsi && nPhi) {
0403     outF << "build and dump Bs" << endl;
0404     BPHRecoBuilder bBs(es);
0405     bBs.setMinPDiffererence(1.0e-5);
0406     bBs.add("JPsi", lJPsi);
0407     bBs.add("Phi", lPhi);
0408     MassSelect mJPsi(2.946916, 3.246916);
0409     MassSelect mPhi(1.009461, 1.029461);
0410     bBs.filter("JPsi", mJPsi);
0411     bBs.filter("Phi", mPhi);
0412     Chi2Select chi2Bs(0.02);
0413     bBs.filter(chi2Bs);
0414     vector<BPHRecoConstCandPtr> lBs = BPHRecoCandidate::build(bBs);
0415     int iBs;
0416     int nBs = lBs.size();
0417     outF << nBs << " Bs cand found" << endl;
0418     // apply kinematic fit
0419     for (iBs = 0; iBs < nBs; ++iBs) {
0420       // get candidate and cast constness away
0421       const BPHRecoCandidate* cptr = lBs[iBs].get();
0422       cptr->kinematicTree("JPsi", 3.096916, 0.000040);
0423     }
0424     for (iBs = 0; iBs < nBs; ++iBs)
0425       dumpRecoCand("Bs", lBs[iBs].get());
0426   }
0427 
0428   // build and dump Bu
0429 
0430   if (nJPsi && nrc) {
0431     outF << "build and dump Bu" << endl;
0432     BPHRecoBuilder bBu(es);
0433     bBu.setMinPDiffererence(1.0e-5);
0434     bBu.add("JPsi", lJPsi);
0435     if (usePF) {
0436       bBu.add("Kaon", BPHRecoBuilder::createCollection(pfCands), 0.493677);
0437     } else if (usePC) {
0438       bBu.add("Kaon", BPHRecoBuilder::createCollection(pcCands), 0.493677);
0439     } else if (useGP) {
0440       bBu.add("Kaon", BPHRecoBuilder::createCollection(gpCands), 0.493677);
0441     }
0442     MassSelect mJPsi(2.946916, 3.246916);
0443     KaonNeutralVeto knv;
0444     bBu.filter("JPsi", mJPsi);
0445     bBu.filter("Kaon", tkPt);
0446     bBu.filter("Kaon", knv);
0447     Chi2Select chi2Bu(0.02);
0448     bBu.filter(chi2Bu);
0449     vector<BPHRecoConstCandPtr> lBu = BPHRecoCandidate::build(bBu);
0450     int iBu;
0451     int nBu = lBu.size();
0452     outF << nBu << " Bu cand found" << endl;
0453     // apply kinematic fit
0454     for (iBu = 0; iBu < nBu; ++iBu) {
0455       // get candidate and cast constness away
0456       const BPHRecoCandidate* cptr = lBu[iBu].get();
0457       cptr->kinematicTree("JPsi", 3.096916, 0.000040);
0458     }
0459     for (iBu = 0; iBu < nBu; ++iBu)
0460       dumpRecoCand("Bu", lBu[iBu].get());
0461   }
0462 
0463   return;
0464 }
0465 
0466 void TestBPHRecoDecay::endJob() {
0467   *fPtr << "TestBPHRecoDecay::endJob" << endl;
0468   TDirectory* currentDir = gDirectory;
0469   TFile file(outHist.c_str(), "RECREATE");
0470   map<string, TH1F*>::iterator iter = histoMap.begin();
0471   map<string, TH1F*>::iterator iend = histoMap.end();
0472   while (iter != iend)
0473     iter++->second->Write();
0474   currentDir->cd();
0475   return;
0476 }
0477 
0478 void TestBPHRecoDecay::dumpRecoCand(const string& name, const BPHRecoCandidate* cand) {
0479   fillHisto(name, cand);
0480 
0481   ostream& outF = *fPtr;
0482 
0483   static string cType = " cowboy";
0484   static string sType = " sailor";
0485   static string dType = "";
0486   string* type;
0487   const BPHPlusMinusCandidate* pmCand = dynamic_cast<const BPHPlusMinusCandidate*>(cand);
0488   if (pmCand != nullptr) {
0489     if (pmCand->isCowboy())
0490       type = &cType;
0491     else
0492       type = &sType;
0493   } else
0494     type = &dType;
0495 
0496   outF << "****** " << name << "   cand mass: " << cand->composite().mass() << " momentum " << cand->composite().px()
0497        << " " << cand->composite().py() << " " << cand->composite().pz() << *type << endl;
0498 
0499   bool validFit = cand->isValidFit();
0500   const RefCountedKinematicTree kt = cand->kinematicTree();
0501   const RefCountedKinematicParticle kp = cand->currentParticle();
0502   if (validFit) {
0503     outF << "****** " << name << " constr mass: " << cand->p4().mass() << " momentum " << cand->p4().px() << " "
0504          << cand->p4().py() << " " << cand->p4().pz() << endl;
0505   }
0506 
0507   const reco::Vertex& vx = cand->vertex();
0508   const reco::Vertex::Point& vp = vx.position();
0509   double chi2 = vx.chi2();
0510   int ndof = lround(vx.ndof());
0511   double prob = TMath::Prob(chi2, ndof);
0512   string tdca = "";
0513   if (pmCand != nullptr) {
0514     stringstream sstr;
0515     sstr << " - " << pmCand->cAppInRPhi().distance();
0516     tdca = sstr.str();
0517   }
0518   outF << "****** " << name << " vertex: " << vx.isFake() << " " << vx.isValid() << " - " << chi2 << " " << ndof << " "
0519        << prob << " - " << vp.X() << " " << vp.Y() << " " << vp.Z() << tdca << endl;
0520 
0521   const vector<string>& dl = cand->daugNames();
0522   int i;
0523   int n = dl.size();
0524   for (i = 0; i < n; ++i) {
0525     const string& name = dl[i];
0526     const reco::Candidate* dp = cand->getDaug(name);
0527     GlobalPoint gp(vp.X(), vp.Y(), vp.Z());
0528     GlobalVector dm(0.0, 0.0, 0.0);
0529     const reco::TransientTrack* tt = cand->getTransientTrack(dp);
0530     if (tt != nullptr) {
0531       TrajectoryStateClosestToPoint tscp = tt->trajectoryStateClosestToPoint(gp);
0532       dm = tscp.momentum();
0533       //      TrajectoryStateOnSurface tsos = tt->stateOnSurface( gp );
0534       //      GlobalVector gv = tsos.globalMomentum();
0535     }
0536     outF << "daughter " << i << " " << name << " " << (dp->charge() > 0 ? '+' : '-') << " momentum: " << dp->px() << " "
0537          << dp->py() << " " << dp->pz() << " - at vertex: " << dm.x() << " " << dm.y() << " " << dm.z() << endl;
0538   }
0539   const vector<string>& dc = cand->compNames();
0540   int j;
0541   int m = dc.size();
0542   for (j = 0; j < m; ++j) {
0543     const string& name = dc[j];
0544     const BPHRecoCandidate* dp = cand->getComp(name).get();
0545     outF << "composite daughter " << j << " " << name << " momentum: " << dp->composite().px() << " "
0546          << dp->composite().py() << " " << dp->composite().pz() << endl;
0547   }
0548 
0549   if (validFit) {
0550     const RefCountedKinematicVertex kv = cand->currentDecayVertex();
0551     GlobalPoint gp = kv->position();
0552     outF << "   kin fit vertex: " << gp.x() << " " << gp.y() << " " << gp.z() << endl;
0553     vector<RefCountedKinematicParticle> dk = kt->finalStateParticles();
0554     int k;
0555     int l = dk.size();
0556     for (k = 0; k < l; ++k) {
0557       const reco::TransientTrack& tt = dk[k]->refittedTransientTrack();
0558       TrajectoryStateClosestToPoint tscp = tt.trajectoryStateClosestToPoint(gp);
0559       GlobalVector dm = tscp.momentum();
0560       //    TrajectoryStateOnSurface tsos = tt.stateOnSurface( gp );
0561       //    GlobalVector dm = tsos.globalMomentum();
0562       outF << "daughter " << k << " refitted: " << dm.x() << " " << dm.y() << " " << dm.z() << endl;
0563     }
0564   }
0565 
0566   return;
0567 }
0568 
0569 void TestBPHRecoDecay::fillHisto(const string& name, const BPHRecoCandidate* cand) {
0570   string mass = "mass";
0571   string mcst = "mcst";
0572   fillHisto(mass + name, cand->composite().mass());
0573   if (cand->isValidFit())
0574     fillHisto(mcst + name, cand->p4().mass());
0575 
0576   const vector<string>& dc = cand->compNames();
0577   int i;
0578   int n = dc.size();
0579   for (i = 0; i < n; ++i) {
0580     const string& daug = dc[i];
0581     const BPHRecoCandidate* dptr = cand->getComp(daug).get();
0582     fillHisto(mass + name + daug, dptr->composite().mass());
0583   }
0584   return;
0585 }
0586 
0587 void TestBPHRecoDecay::fillHisto(const string& name, float x) {
0588   map<string, TH1F*>::iterator iter = histoMap.find(name);
0589   map<string, TH1F*>::iterator iend = histoMap.end();
0590   if (iter == iend)
0591     return;
0592   iter->second->Fill(x);
0593   return;
0594 }
0595 
0596 void TestBPHRecoDecay::createHisto(const string& name, int nbin, float hmin, float hmax) {
0597   histoMap[name] = new TH1F(name.c_str(), name.c_str(), nbin, hmin, hmax);
0598   return;
0599 }
0600 
0601 #include "FWCore/Framework/interface/MakerMacros.h"
0602 
0603 DEFINE_FWK_MODULE(TestBPHRecoDecay);
[ 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!