Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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