File indexing completed on 2021-02-14 13:09:36
0001 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0002
0003 #include "DataFormats/Common/interface/Handle.h"
0004
0005 #include "DataFormats/JetReco/interface/CaloJet.h"
0006 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0007
0008 #include "DataFormats/METReco/interface/CaloMET.h"
0009 #include "DataFormats/METReco/interface/CaloMETCollection.h"
0010
0011 #include "FWCore/Framework/interface/MakerMacros.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013
0014 #include "FWCore/Utilities/interface/InputTag.h"
0015
0016 #include "DQM/DataScouting/interface/RazorVarProducer.h"
0017
0018 #include "TVector3.h"
0019
0020 #include <memory>
0021 #include <vector>
0022
0023
0024
0025
0026 RazorVarProducer::RazorVarProducer(const edm::ParameterSet &iConfig)
0027 : inputTag_(iConfig.getParameter<edm::InputTag>("inputTag")),
0028 inputMetTag_(iConfig.getParameter<edm::InputTag>("inputMetTag")) {
0029
0030 inputTagToken_ = consumes<std::vector<math::XYZTLorentzVector>>(iConfig.getParameter<edm::InputTag>("inputTag"));
0031 inputMetTagToken_ = consumes<reco::CaloMETCollection>(iConfig.getParameter<edm::InputTag>("inputMetTag"));
0032
0033 produces<std::vector<double>>();
0034
0035 LogDebug("") << "Inputs: " << inputTag_.encode() << " " << inputMetTag_.encode() << ".";
0036 }
0037
0038 RazorVarProducer::~RazorVarProducer() {}
0039
0040
0041 void RazorVarProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0042 using namespace std;
0043 using namespace edm;
0044 using namespace reco;
0045
0046
0047 Handle<vector<math::XYZTLorentzVector>> hemispheres;
0048 iEvent.getByToken(inputTagToken_, hemispheres);
0049
0050
0051 Handle<CaloMETCollection> inputMet;
0052 iEvent.getByToken(inputMetTagToken_, inputMet);
0053
0054 std::unique_ptr<std::vector<double>> result(new std::vector<double>);
0055
0056 if (hemispheres.isValid() && inputMet.isValid() && hemispheres->size() > 1) {
0057 TLorentzVector ja(hemispheres->at(0).x(), hemispheres->at(0).y(), hemispheres->at(0).z(), hemispheres->at(0).t());
0058 TLorentzVector jb(hemispheres->at(1).x(), hemispheres->at(1).y(), hemispheres->at(1).z(), hemispheres->at(1).t());
0059
0060 std::vector<math::XYZTLorentzVector> muonVec;
0061 const double MR = CalcMR(ja, jb);
0062 const double R = CalcR(MR, ja, jb, inputMet, muonVec);
0063 result->push_back(MR);
0064 result->push_back(R);
0065 }
0066 iEvent.put(std::move(result));
0067 }
0068
0069 double RazorVarProducer::CalcMR(TLorentzVector ja, TLorentzVector jb) {
0070 if (ja.Pt() <= 0.1)
0071 return -1;
0072
0073 ja.SetPtEtaPhiM(ja.Pt(), ja.Eta(), ja.Phi(), 0.0);
0074 jb.SetPtEtaPhiM(jb.Pt(), jb.Eta(), jb.Phi(), 0.0);
0075
0076 if (ja.Pt() > jb.Pt()) {
0077 TLorentzVector temp = ja;
0078 ja = jb;
0079 jb = temp;
0080 }
0081
0082 double A = ja.P();
0083 double B = jb.P();
0084 double az = ja.Pz();
0085 double bz = jb.Pz();
0086 TVector3 jaT, jbT;
0087 jaT.SetXYZ(ja.Px(), ja.Py(), 0.0);
0088 jbT.SetXYZ(jb.Px(), jb.Py(), 0.0);
0089 double ATBT = (jaT + jbT).Mag2();
0090
0091 double MR = sqrt((A + B) * (A + B) - (az + bz) * (az + bz) -
0092 (jbT.Dot(jbT) - jaT.Dot(jaT)) * (jbT.Dot(jbT) - jaT.Dot(jaT)) / (jaT + jbT).Mag2());
0093
0094 double mybeta = (jbT.Dot(jbT) - jaT.Dot(jaT)) / sqrt(ATBT * ((A + B) * (A + B) - (az + bz) * (az + bz)));
0095
0096 double mygamma = 1. / sqrt(1. - mybeta * mybeta);
0097
0098
0099 return MR * mygamma;
0100 }
0101
0102 double RazorVarProducer::CalcR(double MR,
0103 const TLorentzVector &ja,
0104 const TLorentzVector &jb,
0105 edm::Handle<reco::CaloMETCollection> inputMet,
0106 const std::vector<math::XYZTLorentzVector> &muons) {
0107
0108 TVector3 met;
0109 met.SetPtEtaPhi((inputMet->front()).pt(), 0.0, (inputMet->front()).phi());
0110
0111 std::vector<math::XYZTLorentzVector>::const_iterator muonIt;
0112 for (muonIt = muons.begin(); muonIt != muons.end(); muonIt++) {
0113 TVector3 tmp;
0114 tmp.SetPtEtaPhi(muonIt->pt(), 0, muonIt->phi());
0115 met -= tmp;
0116 }
0117
0118 double MTR = sqrt(0.5 * (met.Mag() * (ja.Pt() + jb.Pt()) - met.Dot(ja.Vect() + jb.Vect())));
0119
0120
0121 return float(MTR) / float(MR);
0122 }
0123
0124 DEFINE_FWK_MODULE(RazorVarProducer);