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 // ------------ method called to produce the data  ------------
0039 void DiJetVarProducer::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   // ## The output collections
0045   // std::unique_ptr<std::vector<double> > dijetvariables(new
0046   // std::vector<double>);
0047   std::unique_ptr<std::vector<math::PtEtaPhiMLorentzVector>> widejets(new std::vector<math::PtEtaPhiMLorentzVector>);
0048
0049   // ## Get jet collection
0050   edm::Handle<reco::CaloJetCollection> calojets_handle;
0051   iEvent.getByToken(inputJetTagToken_, calojets_handle);
0052   // cout << "size: " << calojets_handle->size() << endl;
0053
0054   // ## Wide Jet Algorithm
0055   // At least two jets
0056   if (calojets_handle->size() >= 2) {
0057     TLorentzVector wj1_tmp;
0058     TLorentzVector wj2_tmp;
0059     TLorentzVector wj1;
0060     TLorentzVector wj2;
0061     TLorentzVector wdijet;
0062
0063     //        // Loop over all the input jets
0064     //        for(reco::CaloJetCollection::const_iterator it =
0065     //        calojets_handle->begin(); it != calojets_handle->end(); ++it)
0066     //           {
0067     //     cout << "jet: " << it->pt() << " " << it->eta() << " " << it->phi()
0068     // << endl;
0069     //           }
0070
0071     // Find two leading jets
0072     TLorentzVector jet1, jet2;
0073
0074     reco::CaloJetCollection::const_iterator j1 = calojets_handle->begin();
0075     reco::CaloJetCollection::const_iterator j2 = j1;
0076     ++j2;
0077
0078     jet1.SetPtEtaPhiM(j1->pt(), j1->eta(), j1->phi(), j1->mass());
0079     jet2.SetPtEtaPhiM(j2->pt(), j2->eta(), j2->phi(), j2->mass());
0080
0081     // cout << "j1: " << jet1.Pt() << " " << jet1.Eta() << " " << jet1.Phi() <<
0082     // endl; cout << "j2: " << jet2.Pt() << " " << jet2.Eta() << " " <<
0083     // jet2.Phi() << endl;
0084
0085     // Create wide jets (radiation recovery algorithm)
0086     for (reco::CaloJetCollection::const_iterator it = calojets_handle->begin(); it != calojets_handle->end(); ++it) {
0087       TLorentzVector currentJet;
0088       currentJet.SetPtEtaPhiM(it->pt(), it->eta(), it->phi(), it->mass());
0089
0090       double DeltaR1 = currentJet.DeltaR(jet1);
0091       double DeltaR2 = currentJet.DeltaR(jet2);
0092
0093       if (DeltaR1 < DeltaR2 && DeltaR1 < wideJetDeltaR_) {
0094         wj1_tmp += currentJet;
0095       } else if (DeltaR2 < wideJetDeltaR_) {
0096         wj2_tmp += currentJet;
0097       }
0098     }
0099
0100     // Re-order the wide jets in pT
0101     if (wj1_tmp.Pt() > wj2_tmp.Pt()) {
0102       wj1 = wj1_tmp;
0103       wj2 = wj2_tmp;
0104     } else {
0105       wj1 = wj2_tmp;
0106       wj2 = wj1_tmp;
0107     }
0108
0109     // Create dijet system
0110     wdijet = wj1 + wj2;
0111
0112     //        cout << "j1 wide: " << wj1.Pt() << " " << wj1.Eta() << " " <<
0113     //        wj1.Phi() << " " << wj1.M() << endl; cout << "j2 wide: " <<
0114     //        wj2.Pt() << " " << wj2.Eta() << " " << wj2.Phi() << " " << wj2.M()
0115     //        << endl; cout << "MJJWide: " << wdijet.M() << endl; cout <<
0116     //        "DeltaEtaJJWide: " << fabs(wj1.Eta()-wj2.Eta()) << endl; cout <<
0117     //        "DeltaPhiJJWide: " << fabs(wj1.DeltaPhi(wj2)) << endl;
0118
0119     //        // Put variables in the container
0120     //        dijetvariables->push_back( wdijet.M() );                 //0 =
0121     //        MJJWide dijetvariables->push_back( fabs(wj1.Eta()-wj2.Eta()) );
0122     //        //1 = DeltaEtaJJWide dijetvariables->push_back(
0123     //        fabs(wj1.DeltaPhi(wj2)) );    //2 = DeltaPhiJJWide
0124
0125     // Put widejets in the container
0126     math::PtEtaPhiMLorentzVector wj1math(wj1.Pt(), wj1.Eta(), wj1.Phi(), wj1.M());
0127     math::PtEtaPhiMLorentzVector wj2math(wj2.Pt(), wj2.Eta(), wj2.Phi(), wj2.M());
0128     widejets->push_back(wj1math);
0129     widejets->push_back(wj2math);
0130   }
0131   //    else
0132   //      {
0133   //        // Put variables in the container
0134   //        dijetvariables->push_back( -1 );                //0 = MJJWide
0135   //        dijetvariables->push_back( -1 );                //1 = DeltaEtaJJWide
0136   //        dijetvariables->push_back( -1 );                //2 = DeltaPhiJJWide
0137   //      }
0138
0139   // ## Put objects in the Event
0140   // iEvent.put(std::move(dijetvariables), "dijetvariables");
0141   iEvent.put(std::move(widejets), "widejets");
0142 }
0143
0144 DEFINE_FWK_MODULE(DiJetVarProducer);