00001 #include <Math/VectorUtil.h>
00002 #include "AnalysisDataFormats/TopObjects/interface/TtGenEvent.h"
00003 #include "TopQuarkAnalysis/TopJetCombination/plugins/TtSemiLepHypGeom.h"
00004
00005
00006 TtSemiLepHypGeom::TtSemiLepHypGeom(const edm::ParameterSet& cfg):
00007 TtSemiLepHypothesis( cfg ),
00008 maxNJets_ (cfg.getParameter<int> ("maxNJets" )),
00009 useDeltaR_(cfg.getParameter<bool>("useDeltaR"))
00010 {
00011 if(maxNJets_<4 && maxNJets_!=-1)
00012 throw cms::Exception("WrongConfig")
00013 << "Parameter maxNJets can not be set to " << maxNJets_ << ". \n"
00014 << "It has to be larger than 4 or can be set to -1 to take all jets.";
00015 }
00016
00017 TtSemiLepHypGeom::~TtSemiLepHypGeom() { }
00018
00019 void
00020 TtSemiLepHypGeom::buildHypo(edm::Event& evt,
00021 const edm::Handle<edm::View<reco::RecoCandidate> >& leps,
00022 const edm::Handle<std::vector<pat::MET> >& mets,
00023 const edm::Handle<std::vector<pat::Jet> >& jets,
00024 std::vector<int>& match, const unsigned int iComb)
00025 {
00026 if(leps->empty() || mets->empty() || jets->size()<4){
00027
00028 return;
00029 }
00030
00031 unsigned maxNJets = maxNJets_;
00032 if(maxNJets_ == -1 || (int)jets->size() < maxNJets_) maxNJets = jets->size();
00033
00034 match.clear();
00035 for(unsigned int i=0; i<5; ++i)
00036 match.push_back(-1);
00037
00038
00039
00040
00041
00042 double minDist=-1.;
00043 int lightQ =-1;
00044 int lightQBar=-1;
00045 for(unsigned int idx=0; idx<maxNJets; ++idx){
00046 for(unsigned int jdx=(idx+1); jdx<maxNJets; ++jdx){
00047 double dist = distance((*jets)[idx].p4(), (*jets)[jdx].p4());
00048 if( minDist<0. || dist<minDist ){
00049 minDist=dist;
00050 lightQ =idx;
00051 lightQBar=jdx;
00052 }
00053 }
00054 }
00055
00056 reco::Particle::LorentzVector wHad = (*jets)[lightQ].p4() + (*jets)[lightQBar].p4();
00057
00058
00059
00060
00061
00062 minDist=-1.;
00063 int hadB=-1;
00064 for(unsigned int idx=0; idx<maxNJets; ++idx){
00065
00066 if( (int)idx!=lightQ && (int)idx!=lightQBar ) {
00067 double dist = distance((*jets)[idx].p4(), wHad);
00068 if( minDist<0. || dist<minDist ){
00069 minDist=dist;
00070 hadB=idx;
00071 }
00072 }
00073 }
00074
00075
00076
00077
00078
00079 minDist=-1.;
00080 int lepB=-1;
00081 for(unsigned int idx=0; idx<maxNJets; ++idx){
00082
00083 if( (int)idx!=lightQ && (int)idx!=lightQBar && (int)idx!=hadB ){
00084 double dist = distance((*jets)[idx].p4(), (*leps)[0].p4());
00085 if( minDist<0. || dist<minDist ){
00086 minDist=dist;
00087 lepB=idx;
00088 }
00089 }
00090 }
00091
00092
00093
00094
00095 if( isValid(lightQ, jets) ){
00096 setCandidate(jets, lightQ, lightQ_);
00097 match[TtSemiLepEvtPartons::LightQ] = lightQ;
00098 }
00099
00100 if( isValid(lightQBar, jets) ){
00101 setCandidate(jets, lightQBar, lightQBar_);
00102 match[TtSemiLepEvtPartons::LightQBar] = lightQBar;
00103 }
00104
00105 if( isValid(hadB, jets) ){
00106 setCandidate(jets, hadB, hadronicB_);
00107 match[TtSemiLepEvtPartons::HadB] = hadB;
00108 }
00109
00110 if( isValid(lepB, jets) ){
00111 setCandidate(jets, lepB, leptonicB_);
00112 match[TtSemiLepEvtPartons::LepB] = lepB;
00113 }
00114
00115
00116
00117
00118 setCandidate(leps, 0, lepton_);
00119 match[TtSemiLepEvtPartons::Lepton] = 0;
00120
00121
00122
00123
00124 setCandidate(mets, 0, neutrino_);
00125 }
00126
00127 double
00128 TtSemiLepHypGeom::distance(const math::XYZTLorentzVector& v1, const math::XYZTLorentzVector& v2)
00129 {
00130
00131
00132 if(useDeltaR_) return ROOT::Math::VectorUtil::DeltaR(v1, v2);
00133 return fabs(v1.theta() - v2.theta());
00134 }