Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:34:44

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