DataScouting/src/DiJetVarProducer.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 "FWCore/Framework/interface/MakerMacros.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 #include "FWCore/Utilities/interface/InputTag.h"
0012
0013 #include "DQM/DataScouting/interface/DiJetVarProducer.h"
0014
0015 #include "TLorentzVector.h"
0016 #include "TVector3.h"
0017
0018 #include <memory>
0019 #include <vector>
0020
0021 //
0022 // constructors and destructor
0023 //
0024 DiJetVarProducer::DiJetVarProducer(const edm::ParameterSet &iConfig)
0025     : inputJetTag_(iConfig.getParameter<edm::InputTag>("inputJetTag")),
0026       wideJetDeltaR_(iConfig.getParameter<double>("wideJetDeltaR")) {
0027   // register your products
0028   // produces<std::vector<double> >("dijetvariables");
0029   produces<std::vector<math::PtEtaPhiMLorentzVector>>("widejets");
0030
0031   // set Token(-s)
0032   inputJetTagToken_ = consumes<reco::CaloJetCollection>(iConfig.getParameter<edm::InputTag>("inputJetTag"));
0033
0034   LogDebug("") << "Input Jet Tag: " << inputJetTag_.encode() << " ";
0035   LogDebug("") << "Radius Parameter Wide Jet: " << wideJetDeltaR_ << ".";
0036 }
0037
0038 DiJetVarProducer::~DiJetVarProducer() {}
0039
0040 // ------------ method called to produce the data  ------------
0041 void DiJetVarProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0042   using namespace std;
0043   using namespace edm;
0044   using namespace reco;
0045
0046   // ## The output collections
0047   // std::unique_ptr<std::vector<double> > dijetvariables(new
0048   // std::vector<double>);
0049   std::unique_ptr<std::vector<math::PtEtaPhiMLorentzVector>> widejets(new std::vector<math::PtEtaPhiMLorentzVector>);
0050
0051   // ## Get jet collection
0052   edm::Handle<reco::CaloJetCollection> calojets_handle;
0053   iEvent.getByToken(inputJetTagToken_, calojets_handle);
0054   // cout << "size: " << calojets_handle->size() << endl;
0055
0056   // ## Wide Jet Algorithm
0057   // At least two jets
0058   if (calojets_handle->size() >= 2) {
0059     TLorentzVector wj1_tmp;
0060     TLorentzVector wj2_tmp;
0061     TLorentzVector wj1;
0062     TLorentzVector wj2;
0063     TLorentzVector wdijet;
0064
0065     //        // Loop over all the input jets
0066     //        for(reco::CaloJetCollection::const_iterator it =
0067     //        calojets_handle->begin(); it != calojets_handle->end(); ++it)
0068     //           {
0069     //     cout << "jet: " << it->pt() << " " << it->eta() << " " << it->phi()
0070     // << endl;
0071     //           }
0072
0073     // Find two leading jets
0074     TLorentzVector jet1, jet2;
0075
0076     reco::CaloJetCollection::const_iterator j1 = calojets_handle->begin();
0077     reco::CaloJetCollection::const_iterator j2 = j1;
0078     ++j2;
0079
0080     jet1.SetPtEtaPhiM(j1->pt(), j1->eta(), j1->phi(), j1->mass());
0081     jet2.SetPtEtaPhiM(j2->pt(), j2->eta(), j2->phi(), j2->mass());
0082
0083     // cout << "j1: " << jet1.Pt() << " " << jet1.Eta() << " " << jet1.Phi() <<
0084     // endl; cout << "j2: " << jet2.Pt() << " " << jet2.Eta() << " " <<
0085     // jet2.Phi() << endl;
0086
0087     // Create wide jets (radiation recovery algorithm)
0088     for (reco::CaloJetCollection::const_iterator it = calojets_handle->begin(); it != calojets_handle->end(); ++it) {
0089       TLorentzVector currentJet;
0090       currentJet.SetPtEtaPhiM(it->pt(), it->eta(), it->phi(), it->mass());
0091
0092       double DeltaR1 = currentJet.DeltaR(jet1);
0093       double DeltaR2 = currentJet.DeltaR(jet2);
0094
0095       if (DeltaR1 < DeltaR2 && DeltaR1 < wideJetDeltaR_) {
0096         wj1_tmp += currentJet;
0097       } else if (DeltaR2 < wideJetDeltaR_) {
0098         wj2_tmp += currentJet;
0099       }
0100     }
0101
0102     // Re-order the wide jets in pT
0103     if (wj1_tmp.Pt() > wj2_tmp.Pt()) {
0104       wj1 = wj1_tmp;
0105       wj2 = wj2_tmp;
0106     } else {
0107       wj1 = wj2_tmp;
0108       wj2 = wj1_tmp;
0109     }
0110
0111     // Create dijet system
0112     wdijet = wj1 + wj2;
0113
0114     //        cout << "j1 wide: " << wj1.Pt() << " " << wj1.Eta() << " " <<
0115     //        wj1.Phi() << " " << wj1.M() << endl; cout << "j2 wide: " <<
0116     //        wj2.Pt() << " " << wj2.Eta() << " " << wj2.Phi() << " " << wj2.M()
0117     //        << endl; cout << "MJJWide: " << wdijet.M() << endl; cout <<
0118     //        "DeltaEtaJJWide: " << fabs(wj1.Eta()-wj2.Eta()) << endl; cout <<
0119     //        "DeltaPhiJJWide: " << fabs(wj1.DeltaPhi(wj2)) << endl;
0120
0121     //        // Put variables in the container
0122     //        dijetvariables->push_back( wdijet.M() );                 //0 =
0123     //        MJJWide dijetvariables->push_back( fabs(wj1.Eta()-wj2.Eta()) );
0124     //        //1 = DeltaEtaJJWide dijetvariables->push_back(
0125     //        fabs(wj1.DeltaPhi(wj2)) );    //2 = DeltaPhiJJWide
0126
0127     // Put widejets in the container
0128     math::PtEtaPhiMLorentzVector wj1math(wj1.Pt(), wj1.Eta(), wj1.Phi(), wj1.M());
0129     math::PtEtaPhiMLorentzVector wj2math(wj2.Pt(), wj2.Eta(), wj2.Phi(), wj2.M());
0130     widejets->push_back(wj1math);
0131     widejets->push_back(wj2math);
0132   }
0133   //    else
0134   //      {
0135   //        // Put variables in the container
0136   //        dijetvariables->push_back( -1 );                //0 = MJJWide
0137   //        dijetvariables->push_back( -1 );                //1 = DeltaEtaJJWide
0138   //        dijetvariables->push_back( -1 );                //2 = DeltaPhiJJWide
0139   //      }
0140
0141   // ## Put objects in the Event
0142   // iEvent.put(std::move(dijetvariables), "dijetvariables");
0143   iEvent.put(std::move(widejets), "widejets");
0144 }
0145
0146 DEFINE_FWK_MODULE(DiJetVarProducer);