Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:09:35

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