File indexing completed on 2024-04-06 12:07:02
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
0039 void RazorVarProducer::produce(edm::StreamID, edm::Event &iEvent, const edm::EventSetup &iSetup) const {
0040 using namespace std;
0041 using namespace edm;
0042 using namespace reco;
0043
0044
0045 Handle<vector<math::XYZTLorentzVector>> hemispheres;
0046 iEvent.getByToken(inputTagToken_, hemispheres);
0047
0048
0049 Handle<CaloMETCollection> inputMet;
0050 iEvent.getByToken(inputMetTagToken_, inputMet);
0051
0052 std::unique_ptr<std::vector<double>> result(new std::vector<double>);
0053
0054 if (hemispheres.isValid() && inputMet.isValid() && hemispheres->size() > 1) {
0055 TLorentzVector ja(hemispheres->at(0).x(), hemispheres->at(0).y(), hemispheres->at(0).z(), hemispheres->at(0).t());
0056 TLorentzVector jb(hemispheres->at(1).x(), hemispheres->at(1).y(), hemispheres->at(1).z(), hemispheres->at(1).t());
0057
0058 std::vector<math::XYZTLorentzVector> muonVec;
0059 const double MR = CalcMR(ja, jb);
0060 const double R = CalcR(MR, ja, jb, inputMet, muonVec);
0061 result->push_back(MR);
0062 result->push_back(R);
0063 }
0064 iEvent.put(std::move(result));
0065 }
0066
0067 double RazorVarProducer::CalcMR(TLorentzVector ja, TLorentzVector jb) const {
0068 if (ja.Pt() <= 0.1)
0069 return -1;
0070
0071 ja.SetPtEtaPhiM(ja.Pt(), ja.Eta(), ja.Phi(), 0.0);
0072 jb.SetPtEtaPhiM(jb.Pt(), jb.Eta(), jb.Phi(), 0.0);
0073
0074 if (ja.Pt() > jb.Pt()) {
0075 TLorentzVector temp = ja;
0076 ja = jb;
0077 jb = temp;
0078 }
0079
0080 double A = ja.P();
0081 double B = jb.P();
0082 double az = ja.Pz();
0083 double bz = jb.Pz();
0084 TVector3 jaT, jbT;
0085 jaT.SetXYZ(ja.Px(), ja.Py(), 0.0);
0086 jbT.SetXYZ(jb.Px(), jb.Py(), 0.0);
0087 double ATBT = (jaT + jbT).Mag2();
0088
0089 double MR = sqrt((A + B) * (A + B) - (az + bz) * (az + bz) -
0090 (jbT.Dot(jbT) - jaT.Dot(jaT)) * (jbT.Dot(jbT) - jaT.Dot(jaT)) / (jaT + jbT).Mag2());
0091
0092 double mybeta = (jbT.Dot(jbT) - jaT.Dot(jaT)) / sqrt(ATBT * ((A + B) * (A + B) - (az + bz) * (az + bz)));
0093
0094 double mygamma = 1. / sqrt(1. - mybeta * mybeta);
0095
0096
0097 return MR * mygamma;
0098 }
0099
0100 double RazorVarProducer::CalcR(double MR,
0101 const TLorentzVector &ja,
0102 const TLorentzVector &jb,
0103 edm::Handle<reco::CaloMETCollection> inputMet,
0104 const std::vector<math::XYZTLorentzVector> &muons) const {
0105
0106 TVector3 met;
0107 met.SetPtEtaPhi((inputMet->front()).pt(), 0.0, (inputMet->front()).phi());
0108
0109 std::vector<math::XYZTLorentzVector>::const_iterator muonIt;
0110 for (muonIt = muons.begin(); muonIt != muons.end(); muonIt++) {
0111 TVector3 tmp;
0112 tmp.SetPtEtaPhi(muonIt->pt(), 0, muonIt->phi());
0113 met -= tmp;
0114 }
0115
0116 double MTR = sqrt(0.5 * (met.Mag() * (ja.Pt() + jb.Pt()) - met.Dot(ja.Vect() + jb.Vect())));
0117
0118
0119 return float(MTR) / float(MR);
0120 }
0121
0122 DEFINE_FWK_MODULE(RazorVarProducer);