HiJetAlgos/plugins/HiL1Subtractor.cc

0001 #include "RecoHI/HiJetAlgos/plugins/HiL1Subtractor.h"
0002
0003 #include <vector>
0004
0005 using namespace std;
0006
0007 //
0008 // constants, enums and typedefs
0009 //
0010
0011 //
0012 // static data member definitions
0013 //
0014
0015 //
0016 // constructors and destructor
0017 //
0018 HiL1Subtractor::HiL1Subtractor(const edm::ParameterSet& iConfig)
0019     : jetType_(iConfig.getParameter<std::string>("jetType")),
0020       rhoTagString_(iConfig.getParameter<std::string>("rhoTag"))
0021
0022 {
0023   rhoTag_ = (consumes<std::vector<double> >(rhoTagString_));
0024
0025   if (jetType_ == "CaloJet") {
0026     produces<reco::CaloJetCollection>();
0027     caloJetSrc_ = (consumes<edm::View<reco::CaloJet> >(iConfig.getParameter<edm::InputTag>("src")));
0028   } else if (jetType_ == "PFJet") {
0029     produces<reco::PFJetCollection>();
0030     pfJetSrc_ = (consumes<edm::View<reco::PFJet> >(iConfig.getParameter<edm::InputTag>("src")));
0031   } else if (jetType_ == "GenJet") {
0032     produces<reco::GenJetCollection>();
0033     genJetSrc_ = (consumes<edm::View<reco::GenJet> >(iConfig.getParameter<edm::InputTag>("src")));
0034
0035   } else {
0036     throw cms::Exception("InvalidInput") << "invalid jet type in HiL1Subtractor\n";
0037   }
0038 }
0039
0040 //
0041 // member functions
0042 //
0043
0044 // ------------ method called to produce the data  ------------
0045 void HiL1Subtractor::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0046   // get the input jet collection and create output jet collection
0047
0048   double medianPtkt[12] = {0};
0049   // right now, identical loop for calo and PF jets, should template
0050   if (jetType_ == "GenJet") {
0051     auto jets = std::make_unique<reco::GenJetCollection>();
0052     edm::Handle<edm::View<reco::GenJet> > h_jets;
0053     iEvent.getByToken(genJetSrc_, h_jets);
0054
0055     // Grab appropriate rho, hard coded for the moment
0056     edm::Handle<std::vector<double> > rs;
0057     iEvent.getByToken(rhoTag_, rs);
0058
0059     int rsize = rs->size();
0060
0061     for (int j = 0; j < rsize; j++) {
0062       double medianpt = rs->at(j);
0063       medianPtkt[j] = medianpt;
0064     }
0065
0066     // loop over the jets
0067     int jetsize = h_jets->size();
0068     for (int ijet = 0; ijet < jetsize; ++ijet) {
0069       reco::GenJet jet = ((*h_jets)[ijet]);
0070
0071       double jet_eta = jet.eta();
0072       double jet_et = jet.et();
0073
0074       //std::cout<<" pre-subtracted jet_et "<<jet_et<<std::endl;
0075
0076       if (fabs(jet_eta) <= 3) {
0077         double rho = -999;
0078
0079         if (jet_eta < -2.5 && jet_eta > -3.5)
0080           rho = medianPtkt[2];
0081         if (jet_eta < -1.5 && jet_eta > -2.5)
0082           rho = medianPtkt[3];
0083         if (jet_eta < -0.5 && jet_eta > -1.5)
0084           rho = medianPtkt[4];
0085         if (jet_eta < 0.5 && jet_eta > -0.5)
0086           rho = medianPtkt[5];
0087         if (jet_eta < 1.5 && jet_eta > 0.5)
0088           rho = medianPtkt[6];
0089         if (jet_eta < 2.5 && jet_eta > 1.5)
0090           rho = medianPtkt[7];
0091         if (jet_eta < 3.5 && jet_eta > 2.5)
0092           rho = medianPtkt[8];
0093
0094         double jet_area = jet.jetArea();
0095
0096         double CorrFactor = 0.;
0097         if (rho * jet_area < jet_et)
0098           CorrFactor = 1.0 - rho * jet_area / jet_et;
0099         jet.scaleEnergy(CorrFactor);
0100         jet.setPileup(rho * jet_area);
0101
0102         //std::cout<<"  correction factor "<<1.0 - rho*jet_area/jet_et<<std::endl;
0103       }
0104
0105       //std::cout<<" subtracted jet_et "<<jet.et()<<std::endl;
0106       jets->push_back(jet);
0107     }
0108     iEvent.put(std::move(jets));
0109
0110   } else if (jetType_ == "CaloJet") {
0111     auto jets = std::make_unique<reco::CaloJetCollection>();
0112     edm::Handle<edm::View<reco::CaloJet> > h_jets;
0113     iEvent.getByToken(caloJetSrc_, h_jets);
0114
0115     // Grab appropriate rho, hard coded for the moment
0116     edm::Handle<std::vector<double> > rs;
0117     iEvent.getByToken(rhoTag_, rs);
0118
0119     int rsize = rs->size();
0120
0121     for (int j = 0; j < rsize; j++) {
0122       double medianpt = rs->at(j);
0123       medianPtkt[j] = medianpt;
0124     }
0125
0126     // loop over the jets
0127
0128     int jetsize = h_jets->size();
0129
0130     for (int ijet = 0; ijet < jetsize; ++ijet) {
0131       reco::CaloJet jet = ((*h_jets)[ijet]);
0132
0133       double jet_eta = jet.eta();
0134       double jet_et = jet.et();
0135
0136       //std::cout<<" pre-subtracted jet_et "<<jet_et<<std::endl;
0137
0138       if (fabs(jet_eta) <= 3) {
0139         double rho = -999;
0140
0141         if (jet_eta < -2.5 && jet_eta > -3.5)
0142           rho = medianPtkt[2];
0143         if (jet_eta < -1.5 && jet_eta > -2.5)
0144           rho = medianPtkt[3];
0145         if (jet_eta < -0.5 && jet_eta > -1.5)
0146           rho = medianPtkt[4];
0147         if (jet_eta < 0.5 && jet_eta > -0.5)
0148           rho = medianPtkt[5];
0149         if (jet_eta < 1.5 && jet_eta > 0.5)
0150           rho = medianPtkt[6];
0151         if (jet_eta < 2.5 && jet_eta > 1.5)
0152           rho = medianPtkt[7];
0153         if (jet_eta < 3.5 && jet_eta > 2.5)
0154           rho = medianPtkt[8];
0155
0156         double jet_area = jet.jetArea();
0157
0158         double CorrFactor = 0.;
0159         if (rho * jet_area < jet_et)
0160           CorrFactor = 1.0 - rho * jet_area / jet_et;
0161         jet.scaleEnergy(CorrFactor);
0162         jet.setPileup(rho * jet_area);
0163
0164         //std::cout<<"  correction factor "<<1.0 - rho*jet_area/jet_et<<std::endl;
0165       }
0166
0167       //std::cout<<" subtracted jet_et "<<jet.et()<<std::endl;
0168
0169       jets->push_back(jet);
0170     }
0171     iEvent.put(std::move(jets));
0172
0173   } else if (jetType_ == "PFJet") {
0174     auto jets = std::make_unique<reco::PFJetCollection>();
0175     edm::Handle<edm::View<reco::PFJet> > h_jets;
0176     iEvent.getByToken(pfJetSrc_, h_jets);
0177
0178     // Grab appropriate rho, hard coded for the moment
0179     edm::Handle<std::vector<double> > rs;
0180     iEvent.getByToken(rhoTag_, rs);
0181
0182     int rsize = rs->size();
0183
0184     for (int j = 0; j < rsize; j++) {
0185       double medianpt = rs->at(j);
0186       medianPtkt[j] = medianpt;
0187     }
0188
0189     // loop over the jets
0190
0191     int jetsize = h_jets->size();
0192
0193     for (int ijet = 0; ijet < jetsize; ++ijet) {
0194       reco::PFJet jet = ((*h_jets)[ijet]);
0195
0196       double jet_eta = jet.eta();
0197       double jet_et = jet.et();
0198
0199       //std::cout<<" pre-subtracted jet_et "<<jet_et<<std::endl;
0200
0201       if (fabs(jet_eta) <= 3) {
0202         double rho = -999;
0203
0204         if (jet_eta < -2.5 && jet_eta > -3.5)
0205           rho = medianPtkt[2];
0206         if (jet_eta < -1.5 && jet_eta > -2.5)
0207           rho = medianPtkt[3];
0208         if (jet_eta < -0.5 && jet_eta > -1.5)
0209           rho = medianPtkt[4];
0210         if (jet_eta < 0.5 && jet_eta > -0.5)
0211           rho = medianPtkt[5];
0212         if (jet_eta < 1.5 && jet_eta > 0.5)
0213           rho = medianPtkt[6];
0214         if (jet_eta < 2.5 && jet_eta > 1.5)
0215           rho = medianPtkt[7];
0216         if (jet_eta < 3.5 && jet_eta > 2.5)
0217           rho = medianPtkt[8];
0218
0219         double jet_area = jet.jetArea();
0220
0221         double CorrFactor = 0.;
0222         if (rho * jet_area < jet_et)
0223           CorrFactor = 1.0 - rho * jet_area / jet_et;
0224         jet.scaleEnergy(CorrFactor);
0225         jet.setPileup(rho * jet_area);
0226
0227         //std::cout<<"  correction factor "<<1.0 - rho*jet_area/jet_et<<std::endl;
0228       }
0229
0230       //std::cout<<" subtracted jet_et "<<jet.et()<<std::endl;
0231
0232       jets->push_back(jet);
0233     }
0234     iEvent.put(std::move(jets));
0235   }
0236 }
0237
0238 //define this as a plug-in
0239 DEFINE_FWK_MODULE(HiL1Subtractor);