DataScouting/src/AlphaTVarProducer.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/AlphaTVarProducer.h"
0017
0018 #include "TVector3.h"
0019
0020 #include <memory>
0021 #include <vector>
0022
0023 //
0024 // constructors and destructor
0025 //
0026 AlphaTVarProducer::AlphaTVarProducer(const edm::ParameterSet &iConfig)
0027     : inputJetTag_(iConfig.getParameter<edm::InputTag>("inputJetTag")) {
0028   produces<std::vector<double>>();
0029
0030   // set Token(-s)
0031   inputJetTagToken_ = consumes<reco::CaloJetCollection>(iConfig.getParameter<edm::InputTag>("inputJetTag"));
0032
0033   LogDebug("") << "Inputs: " << inputJetTag_.encode() << " ";
0034 }
0035
0036 // ------------ method called to produce the data  ------------
0037 void AlphaTVarProducer::produce(edm::StreamID, edm::Event &iEvent, const edm::EventSetup &iSetup) const {
0038   using namespace std;
0039   using namespace edm;
0040   using namespace reco;
0041
0042   // get hold of collection of objects
0043   edm::Handle<reco::CaloJetCollection> calojet_handle;
0044   iEvent.getByToken(inputJetTagToken_, calojet_handle);
0045
0046   std::unique_ptr<std::vector<double>> result(new std::vector<double>);
0047   // check the the input collections are available
0048   if (calojet_handle.isValid()) {
0049     std::vector<TLorentzVector> myJets;
0050     reco::CaloJetCollection::const_iterator jetIt;
0051     for (jetIt = calojet_handle->begin(); jetIt != calojet_handle->end(); jetIt++) {
0052       TLorentzVector j;
0053       j.SetPtEtaPhiE(jetIt->pt(), jetIt->eta(), jetIt->phi(), jetIt->energy());
0054       myJets.push_back(j);
0055     }
0056
0057     double alphaT = CalcAlphaT(myJets);
0058     double HT = CalcHT(myJets);
0059
0060     result->push_back(alphaT);
0061     result->push_back(HT);
0062   }
0063   iEvent.put(std::move(result));
0064 }
0065
0066 double AlphaTVarProducer::CalcAlphaT(const std::vector<TLorentzVector> &jets) const {
0067   std::vector<double> ETs;
0068   TVector3 MHT{CalcMHT(jets), 0.0, 0.0};
0069   float HT = CalcHT(jets);
0070   // float HT = 0;
0071   for (unsigned int i = 0; i < jets.size(); i++) {
0072     if (jets[i].Et() > 50. && fabs(jets[i].Eta()) < 2.5)
0073       ETs.push_back(jets[i].Et());
0074     // HT += jets[i].Et();
0075   }
0076   if (ETs.size() < 2.)
0077     return 0.0;
0078   if (ETs.size() > 16.)
0079     return 0.0;
0080   float DHT = deltaHt(ETs);
0081
0082   float AlphaT = alphaT(HT, DHT, MHT.Mag());
0083
0084   return AlphaT;
0085 }
0086
0087 double AlphaTVarProducer::deltaHt(const std::vector<double> &ETs) {
0088   if (ETs.size() > 16.)
0089     return 9999999;
0090   std::vector<double> diff(1 << (ETs.size() - 1), 0.);
0091   for (unsigned i = 0; i < diff.size(); i++)
0092     for (unsigned j = 0; j < ETs.size(); j++)
0093       diff[i] += ETs[j] * (1 - 2 * (int(i >> j) & 1));
0094   std::vector<double>::const_iterator it;
0095   double min = 9999999;
0096   for (it = diff.begin(); it != diff.end(); it++)
0097     if (*it < min)
0098       min = *it;
0099   return min;
0100 }
0101
0102 double AlphaTVarProducer::alphaT(const double HT, const double DHT, const double MHT) {
0103   return 0.5 * (HT - DHT) / sqrt(HT * HT - MHT * MHT);
0104 }
0105
0106 double AlphaTVarProducer::CalcHT(const std::vector<TLorentzVector> &jets) {
0107   double HT = 0;
0108   for (unsigned int i = 0; i < jets.size(); i++) {
0109     if (jets[i].Et() > 50. && fabs(jets[i].Eta()) < 2.5)
0110       HT += jets[i].Et();
0111   }
0112
0113   return HT;
0114 }
0115
0116 double AlphaTVarProducer::CalcMHT(const std::vector<TLorentzVector> &jets) {
0117   TVector3 MHT;
0118   for (unsigned int i = 0; i < jets.size(); i++) {
0119     if (jets[i].Et() > 50. && fabs(jets[i].Eta()) < 2.5)
0120       MHT -= jets[i].Vect();
0121   }
0122   MHT.SetZ(0.0);
0123   return MHT.Mag();
0124 }
0125
0126 DEFINE_FWK_MODULE(AlphaTVarProducer);