PatUtils/src/RazorComputer.cc

0001 #include "PhysicsTools/PatUtils/interface/RazorComputer.h"
0002 #include "TLorentzVector.h"
0003
0004 RazorBox::RazorBox(const CachingVariable::CachingVariableFactoryArg& arg, edm::ConsumesCollector& iC)
0005     : CachingVariable("RazorBox", arg.n, arg.iConfig, iC) {}
0006
0007 void RazorBox::compute(const edm::Event& iEvent) const {}
0008
0009 RazorComputer::RazorComputer(const CachingVariable::CachingVariableFactoryArg& arg, edm::ConsumesCollector& iC)
0010     : VariableComputer(arg, iC) {
0011   jet_ = edm::Service<InputTagDistributorService>()->retrieve("jet", arg.iConfig);
0012   met_ = edm::Service<InputTagDistributorService>()->retrieve("met", arg.iConfig);
0013   jetToken_ = iC.consumes<std::vector<pat::Jet>>(jet_);
0014   metToken_ = iC.consumes<std::vector<pat::MET>>(met_);
0015   pt_ = 40.;
0016   eta_ = 2.4;
0017
0018   declare("MRT", iC);
0019   declare("MR", iC);
0020   declare("R2", iC);
0021   declare("R", iC);
0022 }
0023
0024 namespace razor {
0025   typedef reco::Candidate::LorentzVector LorentzVector;
0026
0027   double CalcMR(const LorentzVector& ja, const LorentzVector& jb) {
0028     double A = ja.P();
0029     double B = jb.P();
0030     double az = ja.Pz();
0031     double bz = jb.Pz();
0032
0033     double temp = sqrt((A + B) * (A + B) - (az + bz) * (az + bz));
0034
0035     return temp;
0036   }
0037   double CalcMTR(const LorentzVector& j1, const LorentzVector& j2, const pat::MET& met) {
0038     double temp = met.et() * (j1.Pt() + j2.Pt()) - met.px() * (j1.X() + j2.X()) - met.py() * (j1.Y() + j2.Y());
0039     temp /= 2.;
0040
0041     temp = sqrt(temp);
0042
0043     return temp;
0044   }
0045 };  // namespace razor
0046
0047 void RazorComputer::compute(const edm::Event& iEvent) const {
0048   //std::cout<<"getting into computation"<<std::endl;
0049
0050   edm::Handle<std::vector<pat::Jet>> jetH;
0051   edm::Handle<std::vector<pat::MET>> metH;
0052   iEvent.getByToken(jetToken_, jetH);
0053   iEvent.getByToken(metToken_, metH);
0054
0055   typedef reco::Candidate::LorentzVector LorentzVector;
0056
0057   std::vector<LorentzVector> jets;
0058   jets.reserve(jetH.product()->size());
0059   for (std::vector<pat::Jet>::const_iterator jetit = jetH.product()->begin(); jetit != jetH.product()->end(); ++jetit) {
0060     if (jetit->et() > pt_ && fabs(jetit->eta()) < eta_) {
0061       jets.push_back(jetit->p4());
0062     }
0063   }
0064
0065   reco::Candidate::LorentzVector HEM_1, HEM_2;
0066
0067   HEM_1.SetPxPyPzE(0.0, 0.0, 0.0, 0.0);
0068   HEM_2.SetPxPyPzE(0.0, 0.0, 0.0, 0.0);
0069
0070   if (jets.size() < 2) {
0071   } else {
0072     unsigned int N_comb = 1;
0073     for (unsigned int i = 0; i < jets.size(); i++) {
0074       N_comb *= 2;
0075     }
0076
0077     double M_temp;
0078
0079     double M_min = 9999999999.0;
0080     int j_count;
0081
0082     for (unsigned int i = 1; i < N_comb - 1; i++) {
0083       LorentzVector j_temp1, j_temp2;
0084       int itemp = i;
0085       j_count = N_comb / 2;
0086       int count = 0;
0087       while (j_count > 0) {
0088         if (itemp / j_count == 1) {
0089           j_temp1 += jets[count];
0090         } else {
0091           j_temp2 += jets[count];
0092         }
0093         itemp -= j_count * (itemp / j_count);
0094         j_count /= 2;
0095         count++;
0096       }
0097
0098       M_temp = j_temp1.M2() + j_temp2.M2();
0099
0100       if (M_temp < M_min) {
0101         M_min = M_temp;
0102         HEM_1 = j_temp1;
0103         HEM_2 = j_temp2;
0104       }
0105     }
0106
0107     if (HEM_2.Pt() > HEM_1.Pt()) {
0108       LorentzVector temp = HEM_1;
0109       HEM_1 = HEM_2;
0110       HEM_2 = temp;
0111     }
0112   }
0113
0114   double mrt = razor::CalcMTR(HEM_1, HEM_2, metH.product()->at(0));
0115   double mr = razor::CalcMR(HEM_1, HEM_2);
0116
0117   assign("MRT", mrt);
0118   assign("MR", mr);
0119   double r = -1, r2 = -1;
0120   if (mr != 0) {
0121     r = mrt / mr;
0122     r2 = r * r;
0123   }
0124   assign("R", r);
0125   assign("R2", r2);
0126
0127   //std::cout<<"MR,R2 "<<mr<<" , "<<r2<<std::endl;
0128 }