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