CMS 3D CMS Logo
HiL1Subtractor.cc
Go to the documentation of this file.
1 #include "RecoHI/HiJetAlgos/plugins/HiL1Subtractor.h"
2 
3 #include <vector>
4 
5 using namespace std;
6 
7 //
8 // constants, enums and typedefs
9 //
10 double puCent[11] = {-5,-4,-3,-2,-1,0,1,2,3,4,5};
11 double medianPtkt[12];
12 
13 //
14 // static data member definitions
15 //
16 
17 //
18 // constructors and destructor
19 //
20 HiL1Subtractor::HiL1Subtractor(const edm::ParameterSet& iConfig) :
21  jetType_ (iConfig.getParameter<std::string>("jetType") ),
22  rhoTagString_ (iConfig.getParameter<std::string>("rhoTag") )
23 
24 {
25  rhoTag_ = (consumes<std::vector<double> > (rhoTagString_));
26 
27  if(jetType_ == "CaloJet"){
28  produces<reco::CaloJetCollection >();
29  caloJetSrc_ = (consumes< edm::View<reco::CaloJet> >(iConfig.getParameter<edm::InputTag>("src")));
30  }
31  else if(jetType_ == "PFJet"){
32  produces<reco::PFJetCollection >();
33  pfJetSrc_ = (consumes< edm::View<reco::PFJet> >(iConfig.getParameter<edm::InputTag>("src")));
34  }
35  else if(jetType_ == "GenJet"){
36  produces<reco::GenJetCollection >();
37  genJetSrc_ = (consumes< edm::View<reco::GenJet> >(iConfig.getParameter<edm::InputTag>("src")));
38 
39  }
40  else{
41  throw cms::Exception("InvalidInput") << "invalid jet type in HiL1Subtractor\n";
42  }
43 
44 }
45 
46 
47 HiL1Subtractor::~HiL1Subtractor()
48 {
49 }
50 
51 
52 //
53 // member functions
54 //
55 
56 // ------------ method called to produce the data ------------
57 void
58 HiL1Subtractor::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
59 {
60 
61  // get the input jet collection and create output jet collection
62 
63  // right now, identical loop for calo and PF jets, should template
64  if(jetType_ == "GenJet"){
65  auto jets = std::make_unique<reco::GenJetCollection>();
66  edm::Handle< edm::View<reco::GenJet> > h_jets;
67  iEvent.getByToken( genJetSrc_, h_jets );
68 
69  // Grab appropriate rho, hard coded for the moment
70  edm::Handle<std::vector<double> > rs;
71  iEvent.getByToken(rhoTag_, rs);
72 
73  int rsize = rs->size();
74 
75  for(int j = 0; j < rsize; j++){
76  double medianpt=rs->at(j);
77  medianPtkt[j]=medianpt;
78  }
79 
80  // loop over the jets
81  int jetsize = h_jets->size();
82  for(int ijet = 0; ijet < jetsize; ++ijet){
83 
84  reco::GenJet jet = ((*h_jets)[ijet]);
85 
86  double jet_eta = jet.eta();
87  double jet_et = jet.et();
88 
89  //std::cout<<" pre-subtracted jet_et "<<jet_et<<std::endl;
90 
91  if(fabs(jet_eta)<=3){
92 
93  double rho=-999;
94 
95  if (jet_eta<-2.5 && jet_eta>-3.5)rho=medianPtkt[2];
96  if (jet_eta<-1.5 && jet_eta>-2.5)rho=medianPtkt[3];
97  if (jet_eta<-0.5 && jet_eta>-1.5)rho=medianPtkt[4];
98  if (jet_eta<0.5 && jet_eta>-0.5)rho=medianPtkt[5];
99  if (jet_eta<1.5 && jet_eta>0.5)rho=medianPtkt[6];
100  if (jet_eta<2.5 && jet_eta>1.5)rho=medianPtkt[7];
101  if (jet_eta<3.5 && jet_eta>2.5)rho=medianPtkt[8];
102 
103  double jet_area = jet.jetArea();
104 
105  double CorrFactor =0.;
106  if(rho*jet_area<jet_et) CorrFactor = 1.0 - rho*jet_area/jet_et;
107  jet.scaleEnergy( CorrFactor );
108  jet.setPileup(rho*jet_area);
109 
110  //std::cout<<" correction factor "<<1.0 - rho*jet_area/jet_et<<std::endl;
111  }
112 
113  //std::cout<<" subtracted jet_et "<<jet.et()<<std::endl;
114  jets->push_back(jet);
115 
116 
117  }
118  iEvent.put(std::move(jets));
119 
120  }else if(jetType_ == "CaloJet"){
121  auto jets = std::make_unique<reco::CaloJetCollection>();
122  edm::Handle< edm::View<reco::CaloJet> > h_jets;
123  iEvent.getByToken( caloJetSrc_, h_jets );
124 
125  // Grab appropriate rho, hard coded for the moment
126  edm::Handle<std::vector<double> > rs;
127  iEvent.getByToken(rhoTag_, rs);
128 
129 
130  int rsize = rs->size();
131 
132  for(int j = 0; j < rsize; j++){
133  double medianpt=rs->at(j);
134  medianPtkt[j]=medianpt;
135  }
136 
137  // loop over the jets
138 
139  int jetsize = h_jets->size();
140 
141  for(int ijet = 0; ijet < jetsize; ++ijet){
142 
143  reco::CaloJet jet = ((*h_jets)[ijet]);
144 
145  double jet_eta = jet.eta();
146  double jet_et = jet.et();
147 
148  //std::cout<<" pre-subtracted jet_et "<<jet_et<<std::endl;
149 
150  if(fabs(jet_eta)<=3){
151 
152  double rho=-999;
153 
154  if (jet_eta<-2.5 && jet_eta>-3.5)rho=medianPtkt[2];
155  if (jet_eta<-1.5 && jet_eta>-2.5)rho=medianPtkt[3];
156  if (jet_eta<-0.5 && jet_eta>-1.5)rho=medianPtkt[4];
157  if (jet_eta<0.5 && jet_eta>-0.5)rho=medianPtkt[5];
158  if (jet_eta<1.5 && jet_eta>0.5)rho=medianPtkt[6];
159  if (jet_eta<2.5 && jet_eta>1.5)rho=medianPtkt[7];
160  if (jet_eta<3.5 && jet_eta>2.5)rho=medianPtkt[8];
161 
162  double jet_area = jet.jetArea();
163 
164  double CorrFactor =0.;
165  if(rho*jet_area<jet_et) CorrFactor = 1.0 - rho*jet_area/jet_et;
166  jet.scaleEnergy( CorrFactor );
167  jet.setPileup(rho*jet_area);
168 
169  //std::cout<<" correction factor "<<1.0 - rho*jet_area/jet_et<<std::endl;
170  }
171 
172  //std::cout<<" subtracted jet_et "<<jet.et()<<std::endl;
173 
174  jets->push_back(jet);
175 
176 
177  }
178  iEvent.put(std::move(jets));
179 
180  }
181  else if(jetType_ == "PFJet"){
182  auto jets = std::make_unique<reco::PFJetCollection>();
183  edm::Handle< edm::View<reco::PFJet> > h_jets;
184  iEvent.getByToken( pfJetSrc_, h_jets );
185 
186  // Grab appropriate rho, hard coded for the moment
187  edm::Handle<std::vector<double> > rs;
188  iEvent.getByToken(rhoTag_,rs);
189 
190 
191  int rsize = rs->size();
192 
193  for(int j = 0; j < rsize; j++){
194  double medianpt=rs->at(j);
195  medianPtkt[j]=medianpt;
196  }
197 
198  // loop over the jets
199 
200  int jetsize = h_jets->size();
201 
202  for(int ijet = 0; ijet < jetsize; ++ijet){
203 
204  reco::PFJet jet = ((*h_jets)[ijet]);
205 
206  double jet_eta = jet.eta();
207  double jet_et = jet.et();
208 
209  //std::cout<<" pre-subtracted jet_et "<<jet_et<<std::endl;
210 
211  if(fabs(jet_eta)<=3){
212 
213  double rho=-999;
214 
215  if (jet_eta<-2.5 && jet_eta>-3.5)rho=medianPtkt[2];
216  if (jet_eta<-1.5 && jet_eta>-2.5)rho=medianPtkt[3];
217  if (jet_eta<-0.5 && jet_eta>-1.5)rho=medianPtkt[4];
218  if (jet_eta<0.5 && jet_eta>-0.5)rho=medianPtkt[5];
219  if (jet_eta<1.5 && jet_eta>0.5)rho=medianPtkt[6];
220  if (jet_eta<2.5 && jet_eta>1.5)rho=medianPtkt[7];
221  if (jet_eta<3.5 && jet_eta>2.5)rho=medianPtkt[8];
222 
223  double jet_area = jet.jetArea();
224 
225  double CorrFactor =0.;
226  if(rho*jet_area<jet_et) CorrFactor = 1.0 - rho*jet_area/jet_et;
227  jet.scaleEnergy( CorrFactor );
228  jet.setPileup(rho*jet_area);
229 
230  //std::cout<<" correction factor "<<1.0 - rho*jet_area/jet_et<<std::endl;
231  }
232 
233  //std::cout<<" subtracted jet_et "<<jet.et()<<std::endl;
234 
235  jets->push_back(jet);
236 
237 
238  }
239  iEvent.put(std::move(jets));
240 
241  }
242 
243 
244 
245 
246 
247 
248 
249 }
250 
251 // ------------ method called once each job just before starting event loop ------------
252 void
253 HiL1Subtractor::beginJob()
254 {
255 }
256 
257 // ------------ method called once each job just after ending the event loop ------------
258 void
259 HiL1Subtractor::endJob() {
260 }
261 
262 //define this as a plug-in
263 DEFINE_FWK_MODULE(HiL1Subtractor);
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
HiL1Subtractor::endJob
virtual void endJob()
Definition: HiL1Subtractor.cc:259
edm::Event::put
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:122
reco::CaloJet
Jets made from CaloTowers.
Definition: CaloJet.h:29
reco::Jet::scaleEnergy
virtual void scaleEnergy(double fScale)
scale energy of the jet
reco::LeafCandidate::eta
virtual double eta() const final
momentum pseudorapidity
Definition: LeafCandidate.h:137
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:457
reco::Jet::setPileup
virtual void setPileup(float fEnergy)
Set pileup energy contribution as calculated by algorithm.
Definition: Jet.h:108
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
puCent
double puCent[11]
Definition: HiL1Subtractor.cc:10
HiL1Subtractor::HiL1Subtractor
HiL1Subtractor(const edm::ParameterSet &)
Definition: HiL1Subtractor.cc:20
reco::PFJet
Jets made from PFObjects.
Definition: PFJet.h:21
HiL1Subtractor::rhoTagString_
std::string rhoTagString_
Definition: HiL1Subtractor.h:70
reco::LeafCandidate::et
virtual double et() const final
transverse energy
Definition: LeafCandidate.h:112
HiL1Subtractor::jetType_
std::string jetType_
Definition: HiL1Subtractor.h:69
iEvent
int iEvent
Definition: GenABIO.cc:230
HiL1Subtractor::genJetSrc_
edm::EDGetTokenT< edm::View< reco::GenJet > > genJetSrc_
Definition: HiL1Subtractor.h:64
fwrapper::jets
vector< PseudoJet > jets
Definition: fastjetfortran_madfks.cc:42
medianPtkt
double medianPtkt[12]
Definition: HiL1Subtractor.cc:11
reco::GenJet
Jets made from MC generator particles.
Definition: GenJet.h:24
HiL1Subtractor::caloJetSrc_
edm::EDGetTokenT< edm::View< reco::CaloJet > > caloJetSrc_
Definition: HiL1Subtractor.h:65
HiL1Subtractor::produce
virtual void produce(edm::Event &, const edm::EventSetup &)
Definition: HiL1Subtractor.cc:58
HiL1Subtractor::pfJetSrc_
edm::EDGetTokenT< edm::View< reco::PFJet > > pfJetSrc_
Definition: HiL1Subtractor.h:66
reco::Jet::jetArea
virtual float jetArea() const
get jet area
Definition: Jet.h:105
HiL1Subtractor::rhoTag_
edm::EDGetTokenT< std::vector< double > > rhoTag_
Definition: HiL1Subtractor.h:72
eostools.move
def move(src, dest)
Definition: eostools.py:510
HiL1Subtractor::beginJob
virtual void beginJob()
Definition: HiL1Subtractor.cc:253
Generated for CMSSW Reference Manual by doxygen 1.8.11