CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
JetAnalyzer_HeavyIons_matching.cc
Go to the documentation of this file.
1 //
2 // Jet Analyzer class for heavy ion jets. for DQM jet analysis monitoring
3 // For CMSSW_7_5_X, especially reading background subtracted jets - jet matching 1 - 1
4 // author: Raghav Kunnawalkam Elayavalli,
5 // April 6th 2015
6 // Rutgers University, email: raghav.k.e at CERN dot CH
7 //
8 
9 
11 
12 using namespace edm;
13 using namespace reco;
14 using namespace std;
15 
16 // declare the constructors:
17 
18 
20  mInputJet1Collection (iConfig.getParameter<edm::InputTag> ("src_Jet1")),
21  mInputJet2Collection (iConfig.getParameter<edm::InputTag> ("src_Jet2")),
22  JetType1 (iConfig.getUntrackedParameter<std::string>("Jet1")),
23  JetType2 (iConfig.getUntrackedParameter<std::string>("Jet2")),
24  mRecoJetPtThreshold (iConfig.getParameter<double> ("recoJetPtThreshold")),
25  mRecoDelRMatch (iConfig.getParameter<double> ("recoDelRMatch")),
26  mRecoJetEtaCut (iConfig.getParameter<double> ("recoJetEtaCut"))
27 {
28  std::string inputCollectionLabelJet1(mInputJet1Collection.label());
29  std::string inputCollectionLabelJet2(mInputJet2Collection.label());
30 
31  //consumes
32 
33  if(std::string("VsCalo") == JetType1) caloJet1Token_ = consumes<reco::CaloJetCollection>(mInputJet1Collection);
34  if(std::string("VsPF") == JetType1) pfJetsToken_ = consumes<reco::PFJetCollection>(mInputJet1Collection);
35  if(std::string("PuCalo") == JetType1) caloJet2Token_ = consumes<reco::CaloJetCollection>(mInputJet1Collection);
36  if(std::string("PuPF") == JetType1) basicJetsToken_ = consumes<reco::BasicJetCollection>(mInputJet1Collection);
37 
38  if(std::string("VsCalo") == JetType2) caloJet1Token_ = consumes<reco::CaloJetCollection>(mInputJet2Collection);
39  if(std::string("VsPF") == JetType2) pfJetsToken_ = consumes<reco::PFJetCollection>(mInputJet2Collection);
40  if(std::string("PuCalo") == JetType2) caloJet2Token_ = consumes<reco::CaloJetCollection>(mInputJet2Collection);
41  if(std::string("PuPF") == JetType2) basicJetsToken_ = consumes<reco::BasicJetCollection>(mInputJet2Collection);
42 
43  // initialize the Jet matching histograms
44 
46  mpT_Jet1_matched = 0;
47  mpT_Jet2_matched = 0;
50 
51  // we need to add histograms which will hold the hadronic and electromagnetic energy content for the unmatched Jets.
52  if(std::string("VsCalo") == JetType1 || std::string("PuCalo") == JetType1){
55  }
56  if(std::string("VsCalo") == JetType2 || std::string("PuCalo") == JetType2) {
59  }
60 
61  if(std::string("VsPF") == JetType1){
67  }
68 
69  if(std::string("VsPF") == JetType2){
75  }
76 
77 
78 }
79 
81  {
82 
83  ibooker.setCurrentFolder("JetMET/HIJetValidation/"+mInputJet1Collection.label()+mInputJet2Collection.label());
84 
85  mpT_ratio_Jet1Jet2 = ibooker.book1D("Ratio_Jet1pT_vs_Jet2pT","",100, 0, 10);
86  mpT_Jet1_matched = ibooker.book1D("Jet1_matched_jet_Spectra","",100, 0, 1000);
87  mpT_Jet2_matched = ibooker.book1D("Jet2_matched_jet_Spectra","",100, 0, 1000);
88  mpT_Jet1_unmatched = ibooker.book1D("Jet1_unmatched_jet_Spectra","",100, 0, 1000);
89  mpT_Jet2_unmatched = ibooker.book1D("Jet2_unmatched_jet_Spectra","",100, 0, 1000);
90 
91  if(std::string("VsCalo") == JetType1 || std::string("PuCalo") == JetType1){
92  mHadEnergy_Jet1_unmatched = ibooker.book1D("HadEnergy_Jet1_unmatched", "HadEnergy_Jet1_unmatched", 50, 0, 200);
93  mEmEnergy_Jet1_unmatched = ibooker.book1D("EmEnergy_Jet1_unmatched", "EmEnergy_Jet1_unmatched", 50, 0, 200);
94  }
95 
96  if(std::string("VsCalo") == JetType2 || std::string("PuCalo") == JetType2){
97  mHadEnergy_Jet2_unmatched = ibooker.book1D("HadEnergy_Jet2_unmatched", "HadEnergy_Jet2_unmatched", 50, 0, 200);
98  mEmEnergy_Jet2_unmatched = ibooker.book1D("EmEnergy_Jet2_unmatched", "EmEnergy_Jet2_unmatched", 50, 0, 200);
99  }
100 
101  if(std::string("VsPF") == JetType1){
102  mChargedHadronEnergy_Jet1_unmatched = ibooker.book1D("ChargedHadronEnergy_Jet1_unmatched", "charged HAD energy unmatched Jet1", 50, 0, 100);
103  mNeutralHadronEnergy_Jet1_unmatched = ibooker.book1D("ChargedHadronEnergy_Jet1_unmatched", "charged HAD energy unmatched Jet1", 50, 0, 100);
104  mChargedEmEnergy_Jet1_unmatched = ibooker.book1D("ChargedEmEnergy_Jet1_unmatched", "charged EM energy unmatched Jet1", 50, 0, 100);
105  mNeutralEmEnergy_Jet1_unmatched = ibooker.book1D("ChargedEmEnergy_Jet1_unmatched", "charged EM energy unmatched Jet1", 50, 0, 100);
106  mChargedMuEnergy_Jet1_unmatched = ibooker.book1D("ChargedMuEnergy_Jet1_unmatched", "charged Mu energy unmatched Jet1", 50, 0, 100);
107  }
108 
109  if(std::string("VsPF") == JetType2){
110  mChargedHadronEnergy_Jet2_unmatched = ibooker.book1D("ChargedHadronEnergy_Jet2_unmatched", "charged HAD energy unmatched Jet2", 50, 0, 100);
111  mNeutralHadronEnergy_Jet2_unmatched = ibooker.book1D("ChargedHadronEnergy_Jet2_unmatched", "charged HAD energy unmatched Jet2", 50, 0, 100);
112  mChargedEmEnergy_Jet2_unmatched = ibooker.book1D("ChargedEmEnergy_Jet2_unmatched", "charged EM energy unmatched Jet2", 50, 0, 100);
113  mNeutralEmEnergy_Jet2_unmatched = ibooker.book1D("ChargedEmEnergy_Jet2_unmatched", "charged EM energy unmatched Jet2", 50, 0, 100);
114  mChargedMuEnergy_Jet2_unmatched = ibooker.book1D("ChargedMuEnergy_Jet2_unmatched", "charged Mu energy unmatched Jet2", 50, 0, 100);
115  }
116 
117  if (mOutputFile.empty ())
118  LogInfo("OutputInfo") << " Histograms will NOT be saved";
119  else
120  LogInfo("OutputInfo") << " Histograms will be saved to file:" << mOutputFile;
121 
122 
123  }
124 
125 
126 
127 //------------------------------------------------------------------------------
128 // ~JetAnalyzer_HeavyIons
129 //------------------------------------------------------------------------------
131 
132 
133 //------------------------------------------------------------------------------
134 // beginJob
135 //------------------------------------------------------------------------------
136 //void JetAnalyzer_HeavyIons_matching::beginJob() {
137 // std::cout<<"inside the begin job function"<<endl;
138 //}
139 
140 
141 //------------------------------------------------------------------------------
142 // endJob
143 //------------------------------------------------------------------------------
144 //void JetAnalyzer_HeavyIons_matching::endJob()
145 //{
146 // if (!mOutputFile.empty() && &*edm::Service<DQMStore>())
147 // {
148 // edm::Service<DQMStore>()->save(mOutputFile);
149 // }
150 //}
151 
152 
153 //------------------------------------------------------------------------------
154 // analyze
155 //------------------------------------------------------------------------------
157 {
158 
159  // Get the Jet collection
160  //----------------------------------------------------------------------------
161 
162  std::vector<Jet> recoJet1;
163  recoJet1.clear();
164  std::vector<Jet> recoJet2;
165  recoJet2.clear();
166 
172 
173  if(std::string("VsCalo") == JetType1) {
174  mEvent.getByToken(caloJet1Token_, caloJet1);
175  for (unsigned ijet=0; ijet<caloJet1->size(); ++ijet) recoJet1.push_back((*caloJet1)[ijet]);
176  }
177  if(std::string("PuCalo") == JetType1) {
178  mEvent.getByToken(caloJet2Token_, caloJet1);
179  for (unsigned ijet=0; ijet<caloJet1->size(); ++ijet) recoJet1.push_back((*caloJet1)[ijet]);
180  }
181  if(std::string("VsPF") == JetType1) {
182  mEvent.getByToken(pfJetsToken_, pfJets);
183  for (unsigned ijet=0; ijet<pfJets->size(); ++ijet) recoJet1.push_back((*pfJets)[ijet]);
184  }
185  if(std::string("PuPF") == JetType1) {
186  mEvent.getByToken(basicJetsToken_, basicJets);
187  for (unsigned ijet=0; ijet<basicJets->size(); ++ijet) recoJet1.push_back((*basicJets)[ijet]);
188  }
189 
190  if(std::string("VsCalo") == JetType2) {
191  mEvent.getByToken(caloJet1Token_, caloJet2);
192  for (unsigned ijet=0; ijet<caloJet2->size(); ++ijet) recoJet2.push_back((*caloJet2)[ijet]);
193  }
194  if(std::string("PuCalo") == JetType2) {
195  mEvent.getByToken(caloJet2Token_, caloJet2);
196  for (unsigned ijet=0; ijet<caloJet2->size(); ++ijet) recoJet2.push_back((*caloJet2)[ijet]);
197  }
198  if(std::string("VsPF") == JetType2) {
199  mEvent.getByToken(pfJetsToken_, pfJets);
200  for (unsigned ijet=0; ijet<pfJets->size(); ++ijet) recoJet2.push_back((*pfJets)[ijet]);
201  }
202  if(std::string("PuPF") == JetType2) {
203  mEvent.getByToken(basicJetsToken_, basicJets);
204  for (unsigned ijet=0; ijet<basicJets->size(); ++ijet) recoJet2.push_back((*basicJets)[ijet]);
205  }
206 
207  // start to perform the matching - between recoJet1 and recoJet2.
208 
209  Int_t Jet1_nref = recoJet1.size();
210  Int_t Jet2_nref = recoJet2.size();
211 
212  int jet1 = 0;
213  int jet2 = 0;
214 
215  std::vector <MyJet> vJet1, vJet2;
216  std::vector <int> Jet1_ID(Jet1_nref), Jet2_ID(Jet2_nref);
217 
218  if(Jet1_nref == 0 || Jet2_nref == 0) return;
219 
220  for(unsigned ijet1 = 0; ijet1 < recoJet1.size(); ++ijet1){
221 
222 
223  if(recoJet1[ijet1].pt() < mRecoJetPtThreshold) continue;
224  if(fabs(recoJet1[ijet1].eta()) < mRecoJetEtaCut) continue;
225 
226  MyJet JET1;
227  JET1.eta = recoJet1[ijet1].eta();
228  JET1.phi = recoJet1[ijet1].phi();
229  JET1.pt = recoJet1[ijet1].pt();
230  JET1.id = ijet1;
231 
232  vJet1.push_back(JET1);
233  jet1++;
234 
235  }// first jet loop
236 
237  for(unsigned ijet2 = 0; ijet2 < recoJet2.size(); ++ijet2){
238 
239  if(recoJet2[ijet2].pt() < mRecoJetPtThreshold) continue;
240  if(fabs(recoJet2[ijet2].eta()) < mRecoJetEtaCut) continue;
241 
242  MyJet JET2;
243  JET2.eta = recoJet1[ijet2].eta();
244  JET2.phi = recoJet1[ijet2].phi();
245  JET2.pt = recoJet1[ijet2].pt();
246  JET2.id = ijet2;
247 
248  vJet2.push_back(JET2);
249  jet2++;
250 
251  }// second jet loop
252 
253  bool onlyJet1 = (jet1>0 && jet2==0) ? true : false;
254  bool onlyJet2 = (jet1==0 && jet2 >0) ? true : false;
255  bool bothJet1Jet2 = (jet1>0 && jet2 >0) ? true : false;
256 
257  int matchedJets = 0;
258  int unmatchedJet1 = 0;
259  int unmatchedJet2 = 0;
260 
261  std::vector < MyJet >::const_iterator iJet, jJet;
262 
263  if(onlyJet1) {
264 
265  for(iJet = vJet1.begin(); iJet != vJet1.end(); ++iJet){
266 
267  int pj = (*iJet).id;
268 
269  mpT_Jet1_unmatched->Fill(recoJet1[pj].pt());
270 
271  }
272 
273  }else if(onlyJet2) {
274 
275  for(iJet = vJet2.begin(); iJet != vJet2.end(); ++iJet){
276 
277  int cj = (*iJet).id;
278 
279  mpT_Jet2_unmatched->Fill(recoJet2[cj].pt());
280  }
281 
282  }else if (bothJet1Jet2){
283 
284  ABMatchedJets mABMatchedJets;
285 
286  for(iJet = vJet1.begin(); iJet != vJet1.end(); ++iJet){
287  for(jJet = vJet2.begin(); jJet != vJet2.end(); ++jJet){
288  mABMatchedJets.insert(std::make_pair(*iJet, *jJet));
289  }
290  }
291 
292  ABItr itr;
293  // matched Jets matching Jet 1 to Jet 2
294  for(itr = mABMatchedJets.begin(); itr != mABMatchedJets.end(); ++itr){
295 
296  ABJetPair jetpair = (*itr);
297  MyJet Aj = jetpair.first;
298  MyJet Bj = jetpair.second;
299 
300  float delr = JetAnalyzer_HeavyIons_matching::deltaRR(Bj.eta, Bj.phi, Aj.eta, Aj.phi);
301 
302  if( delr < mRecoDelRMatch && Jet1_ID[Aj.id] == 0){
303 
304  mpT_ratio_Jet1Jet2->Fill((Float_t) recoJet2[Bj.id].pt()/recoJet1[Aj.id].pt());
305 
306  mpT_Jet1_matched->Fill(recoJet1[Aj.id].pt());
307  mpT_Jet2_matched->Fill(recoJet2[Bj.id].pt());
308 
309  Jet1_ID[Aj.id] = 1;
310  Jet2_ID[Bj.id] = 1;
311 
312  matchedJets++;
313 
314 
315  }
316 
317  }
318 
319  // for unmatched Jets
320  for(itr = mABMatchedJets.begin(); itr != mABMatchedJets.end(); ++itr){
321 
322  ABJetPair jetpair = (*itr);
323 
324  MyJet Aj = jetpair.first;
325  MyJet Bj = jetpair.second;
326 
327  if(Jet1_ID[Aj.id] == 0) {
328 
329  mpT_Jet1_unmatched->Fill(recoJet1[Aj.id].pt());
330  unmatchedJet1++;
331  Jet1_ID[Aj.id] = 1;
332 
333  if(std::string("VsCalo") == JetType1 || std::string("PuCalo") == JetType1){
334  mHadEnergy_Jet1_unmatched->Fill((*caloJet1)[Aj.id].hadEnergyInHO()
335  + (*caloJet1)[Aj.id].hadEnergyInHB()
336  + (*caloJet1)[Aj.id].hadEnergyInHF()
337  + (*caloJet1)[Aj.id].hadEnergyInHE()
338  );
339  mEmEnergy_Jet1_unmatched->Fill((*caloJet1)[Aj.id].emEnergyInEB()
340  + (*caloJet1)[Aj.id].emEnergyInEE()
341  + (*caloJet1)[Aj.id].emEnergyInHF()
342  );
343  }
344 
345  if(std::string("VsPF") == JetType1){
346  mChargedHadronEnergy_Jet1_unmatched->Fill((*pfJets)[Aj.id].chargedHadronEnergy());
347  mNeutralHadronEnergy_Jet1_unmatched->Fill((*pfJets)[Aj.id].neutralHadronEnergy());
348  mChargedEmEnergy_Jet1_unmatched->Fill((*pfJets)[Aj.id].chargedEmEnergy());
349  mNeutralEmEnergy_Jet1_unmatched->Fill((*pfJets)[Aj.id].neutralEmEnergy());
350  mChargedMuEnergy_Jet1_unmatched->Fill((*pfJets)[Aj.id].chargedMuEnergy());
351  }
352 
353  }
354 
355  if(Jet2_ID[Bj.id] == 0) {
356 
357  mpT_Jet2_unmatched->Fill(recoJet2[Bj.id].pt());
358  unmatchedJet2++;
359  Jet2_ID[Bj.id] = 2;
360  if(std::string("VsCalo") == JetType2 || std::string("PuCalo") == JetType2){
361  mHadEnergy_Jet2_unmatched->Fill((*caloJet2)[Bj.id].hadEnergyInHO()
362  + (*caloJet2)[Bj.id].hadEnergyInHB()
363  + (*caloJet2)[Bj.id].hadEnergyInHF()
364  + (*caloJet2)[Bj.id].hadEnergyInHE()
365  );
366  mEmEnergy_Jet2_unmatched->Fill((*caloJet2)[Bj.id].emEnergyInEB()
367  + (*caloJet2)[Bj.id].emEnergyInEE()
368  + (*caloJet2)[Bj.id].emEnergyInHF()
369  );
370  }
371 
372 
373  if(std::string("VsPF") == JetType2){
374  mChargedHadronEnergy_Jet2_unmatched->Fill((*pfJets)[Bj.id].chargedHadronEnergy());
375  mNeutralHadronEnergy_Jet2_unmatched->Fill((*pfJets)[Bj.id].neutralHadronEnergy());
376  mChargedEmEnergy_Jet2_unmatched->Fill((*pfJets)[Bj.id].chargedEmEnergy());
377  mNeutralEmEnergy_Jet2_unmatched->Fill((*pfJets)[Bj.id].neutralEmEnergy());
378  mChargedMuEnergy_Jet2_unmatched->Fill((*pfJets)[Bj.id].chargedMuEnergy());
379  }
380 
381  }
382 
383  }
384 
385  }// both Jet1 and Jet2 in the event.
386 
387 
388 }
389 
390 
391 
edm::EDGetTokenT< reco::PFJetCollection > pfJetsToken_
std::multiset< ABJetPair, CompareMatchedJets > ABMatchedJets
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:464
void analyze(const edm::Event &, const edm::EventSetup &) override
void Fill(long long x)
MonitorElement * book1D(Args &&...args)
Definition: DQMStore.h:115
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:274
edm::EDGetTokenT< reco::CaloJetCollection > caloJet1Token_
static float deltaRR(float eta1, float phi1, float eta2, float phi2)
std::string const & label() const
Definition: InputTag.h:43
JetAnalyzer_HeavyIons_matching(const edm::ParameterSet &)
std::multiset< ABJetPair >::iterator ABItr
edm::EDGetTokenT< reco::CaloJetCollection > caloJet2Token_
tuple pfJets
Definition: pfJets_cff.py:8
Definition: Run.h:43
edm::EDGetTokenT< reco::BasicJetCollection > basicJetsToken_