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 // ------------ method called to produce the data  ------------
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   // get hold of collection of objects
0045   Handle<vector<math::XYZTLorentzVector>> hemispheres;
0046   iEvent.getByToken(inputTagToken_, hemispheres);
0047
0048   // get hold of the MET Collection
0049   Handle<CaloMETCollection> inputMet;
0050   iEvent.getByToken(inputMetTagToken_, inputMet);
0051
0052   std::unique_ptr<std::vector<double>> result(new std::vector<double>);
0053   // check the the input collections are available
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   // use gamma times MRstar
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   // now we can calculate MTR
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   // filter events
0119   return float(MTR) / float(MR);  // R
0120 }
0121
0122 DEFINE_FWK_MODULE(RazorVarProducer);