DataScouting/src/RazorVarProducer.cc

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 // constructors and destructor
0025 //
0026 RazorVarProducer::RazorVarProducer(const edm::ParameterSet &iConfig)
0027     : inputTag_(iConfig.getParameter<edm::InputTag>("inputTag")),
0028       inputMetTag_(iConfig.getParameter<edm::InputTag>("inputMetTag")) {
0029   // set Token(-s)
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 // ------------ method called to produce the data  ------------
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   // get hold of collection of objects
0047   Handle<vector<math::XYZTLorentzVector>> hemispheres;
0048   iEvent.getByToken(inputTagToken_, hemispheres);
0049
0050   // get hold of the MET Collection
0051   Handle<CaloMETCollection> inputMet;
0052   iEvent.getByToken(inputMetTagToken_, inputMet);
0053
0054   std::unique_ptr<std::vector<double>> result(new std::vector<double>);
0055   // check the the input collections are available
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   // use gamma times MRstar
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   // now we can calculate MTR
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   // filter events
0121   return float(MTR) / float(MR);  // R
0122 }
0123
0124 DEFINE_FWK_MODULE(RazorVarProducer);