ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/CmsHi/JetAnalysis/src/HiInclusiveJetAnalyzer.cc
Revision: 1.6
Committed: Wed Nov 2 19:41:13 2011 UTC (13 years, 6 months ago) by yjlee
Content type: text/plain
Branch: MAIN
Changes since 1.5: +1 -1 lines
Log Message:
Bug Fix

File Contents

# User Rev Content
1 yilmaz 1.1 /*
2     Based on the jet response analyzer
3     Modified by Matt Nguyen, November 2010
4    
5     */
6    
7     #include "CmsHi/JetAnalysis/interface/HiInclusiveJetAnalyzer.h"
8    
9     #include "FWCore/Framework/interface/EventSetup.h"
10     #include "FWCore/Framework/interface/ESHandle.h"
11     #include <Math/DistFunc.h>
12     #include "TMath.h"
13    
14     #include "FWCore/Framework/interface/ESHandle.h"
15    
16     #include "FWCore/MessageLogger/interface/MessageLogger.h"
17     #include "FWCore/Utilities/interface/Exception.h"
18     #include "FWCore/Framework/interface/EventSetup.h"
19    
20     #include "DataFormats/HeavyIonEvent/interface/CentralityBins.h"
21     #include "DataFormats/HeavyIonEvent/interface/Centrality.h"
22    
23     #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
24     #include "DataFormats/PatCandidates/interface/Jet.h"
25     #include "DataFormats/JetReco/interface/CaloJetCollection.h"
26     #include "SimDataFormats/HiGenData/interface/GenHIEvent.h"
27     #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
28     #include "DataFormats/JetReco/interface/GenJetCollection.h"
29     #include "DataFormats/VertexReco/interface/Vertex.h"
30     #include "DataFormats/Math/interface/deltaPhi.h"
31     #include "DataFormats/Math/interface/deltaR.h"
32    
33     #include "DataFormats/Common/interface/TriggerResults.h"
34     #include "FWCore/Common/interface/TriggerNames.h"
35     #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetup.h"
36     #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
37     #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h"
38     #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMap.h"
39     #include "L1Trigger/GlobalTrigger/interface/L1GlobalTrigger.h"
40    
41     using namespace std;
42     using namespace edm;
43     using namespace reco;
44    
45     HiInclusiveJetAnalyzer::HiInclusiveJetAnalyzer(const edm::ParameterSet& iConfig) {
46    
47    
48     jetTag_ = iConfig.getParameter<InputTag>("jetTag");
49     vtxTag_ = iConfig.getUntrackedParameter<edm::InputTag>("vtxTag",edm::InputTag("hiSelectedVertex"));
50    
51     isMC_ = iConfig.getUntrackedParameter<bool>("isMC",false);
52    
53     if(isMC_){
54     genjetTag_ = iConfig.getParameter<InputTag>("genjetTag");
55     eventInfoTag_ = iConfig.getParameter<InputTag>("eventInfoTag");
56     }
57     verbose_ = iConfig.getUntrackedParameter<bool>("verbose",false);
58    
59     useCentrality_ = iConfig.getUntrackedParameter<bool>("useCentrality",false);
60 yjlee 1.4 useVtx_ = iConfig.getUntrackedParameter<bool>("useVtx",false);
61 yilmaz 1.1 useJEC_ = iConfig.getUntrackedParameter<bool>("useJEC",true);
62 yilmaz 1.2 usePat_ = iConfig.getUntrackedParameter<bool>("usePAT",true);
63 yilmaz 1.1
64     if(!isMC_){
65     L1gtReadout_ = iConfig.getParameter<edm::InputTag>("L1gtReadout");
66     hltResName_ = iConfig.getUntrackedParameter<string>("hltTrgResults","TriggerResults::HLT");
67    
68    
69     if (iConfig.exists("hltTrgNames"))
70     hltTrgNames_ = iConfig.getUntrackedParameter<vector<string> >("hltTrgNames");
71    
72     if (iConfig.exists("hltProcNames"))
73     hltProcNames_ = iConfig.getUntrackedParameter<vector<string> >("hltProcNames");
74     else {
75     hltProcNames_.push_back("FU");
76     hltProcNames_.push_back("HLT");
77     }
78     }
79    
80     cout<<" jet collection : "<<jetTag_<<endl;
81     if(isMC_)cout<<" genjet collection : "<<genjetTag_<<endl;
82    
83    
84    
85     }
86    
87    
88    
89     HiInclusiveJetAnalyzer::~HiInclusiveJetAnalyzer() { }
90    
91    
92    
93     void
94     HiInclusiveJetAnalyzer::beginRun(const edm::Run& run,
95     const edm::EventSetup & es) {}
96    
97     void
98     HiInclusiveJetAnalyzer::beginJob() {
99    
100     centrality_ = 0;
101    
102     //string jetTagName = jetTag_.label()+"_tree";
103     string jetTagTitle = jetTag_.label()+" Jet Analysis Tree";
104     t = fs1->make<TTree>("t",jetTagTitle.c_str());
105    
106 yjlee 1.6 // TTree* t= new TTree("t","Jet Response Analyzer");
107 yjlee 1.4 //t->Branch("run",&jets_.run,"run/I");
108 yilmaz 1.1 t->Branch("evt",&jets_.evt,"evt/I");
109 yjlee 1.4 //t->Branch("lumi",&jets_.lumi,"lumi/I");
110 yilmaz 1.1 t->Branch("b",&jets_.b,"b/F");
111 yjlee 1.4 if (useVtx_) {
112     t->Branch("vx",&jets_.vx,"vx/F");
113     t->Branch("vy",&jets_.vy,"vy/F");
114     t->Branch("vz",&jets_.vz,"vz/F");
115     }
116     if (useCentrality_) {
117     t->Branch("hf",&jets_.hf,"hf/F");
118     t->Branch("nref",&jets_.nref,"nref/I");
119     t->Branch("bin",&jets_.bin,"bin/I");
120     }
121 yilmaz 1.1 t->Branch("rawpt",jets_.rawpt,"rawpt[nref]/F");
122     t->Branch("jtpt",jets_.jtpt,"jtpt[nref]/F");
123     t->Branch("jteta",jets_.jteta,"jteta[nref]/F");
124     t->Branch("jty",jets_.jty,"jty[nref]/F");
125     t->Branch("jtphi",jets_.jtphi,"jtphi[nref]/F");
126 yilmaz 1.2 t->Branch("jtpu",jets_.jtpu,"jtpu[nref]/F");
127 yilmaz 1.1
128     if(isMC_){
129     t->Branch("pthat",&jets_.pthat,"pthat/F");
130    
131     // Only matched gen jets
132     t->Branch("refpt",jets_.refpt,"refpt[nref]/F");
133     t->Branch("refeta",jets_.refeta,"refeta[nref]/F");
134     t->Branch("refy",jets_.refy,"refy[nref]/F");
135     t->Branch("refphi",jets_.refphi,"refphi[nref]/F");
136     t->Branch("refdphijt",jets_.refdphijt,"refdphijt[nref]/F");
137     t->Branch("refdrjt",jets_.refdrjt,"refdrjt[nref]/F");
138     // matched parton
139     t->Branch("refparton_pt",jets_.refparton_pt,"refparton_pt[nref]/F");
140     t->Branch("refparton_flavor",jets_.refparton_flavor,"refparton_flavor[nref]/F");
141    
142     // For all gen jets, matched or unmatched
143     t->Branch("ngen",&jets_.ngen,"ngen/I");
144     t->Branch("genmatchindex",jets_.genmatchindex,"genmatchindex[ngen]/I");
145     t->Branch("genpt",jets_.genpt,"genpt[ngen]/F");
146     t->Branch("geneta",jets_.geneta,"geneta[ngen]/F");
147     t->Branch("geny",jets_.geny,"geny[ngen]/F");
148     t->Branch("genphi",jets_.genphi,"genphi[ngen]/F");
149     t->Branch("gendphijt",jets_.gendphijt,"gendphijt[ngen]/F");
150     t->Branch("gendrjt",jets_.gendrjt,"gendrjt[ngen]/F");
151     }
152    
153     if(!isMC_){
154     t->Branch("nL1TBit",&jets_.nL1TBit,"nL1TBit/I");
155     t->Branch("l1TBit",jets_.l1TBit,"l1TBit[nL1TBit]/O");
156    
157     t->Branch("nL1ABit",&jets_.nL1ABit,"nL1ABit/I");
158     t->Branch("l1ABit",jets_.l1ABit,"l1ABit[nL1ABit]/O");
159    
160     t->Branch("nHLTBit",&jets_.nHLTBit,"nHLTBit/I");
161     t->Branch("hltBit",jets_.hltBit,"hltBit[nHLTBit]/O");
162    
163     }
164    
165     TH1D::SetDefaultSumw2();
166    
167    
168     }
169    
170    
171     void
172     HiInclusiveJetAnalyzer::analyze(const Event& iEvent,
173     const EventSetup& iSetup) {
174    
175     int event = iEvent.id().event();
176     int run = iEvent.id().run();
177     int lumi = iEvent.id().luminosityBlock();
178    
179     jets_.run = run;
180     jets_.evt = event;
181     jets_.lumi = lumi;
182    
183     LogDebug("HiInclusiveJetAnalyzer")<<"START event: "<<event<<" in run "<<run<<endl;
184    
185     int bin = -1;
186     double hf = 0.;
187     double b = 999.;
188    
189     if(useCentrality_){
190     if(!centrality_) centrality_ = new CentralityProvider(iSetup);
191     centrality_->newEvent(iEvent,iSetup); // make sure you do this first in every event
192     //double c = centrality_->centralityValue();
193     const reco::Centrality *cent = centrality_->raw();
194    
195     hf = cent->EtHFhitSum();
196    
197     bin = centrality_->getBin();
198     b = centrality_->bMean();
199     }
200    
201    
202    
203    
204     //double jetPtMin = 35.;
205    
206    
207     // loop the events
208    
209     jets_.bin = bin;
210     jets_.hf = hf;
211    
212    
213     if (useVtx_) {
214     edm::Handle<vector<reco::Vertex> >vertex;
215     iEvent.getByLabel(vtxTag_, vertex);
216    
217     if(vertex->size()>0) {
218     jets_.vx=vertex->begin()->x();
219     jets_.vy=vertex->begin()->y();
220     jets_.vz=vertex->begin()->z();
221     }
222     }
223    
224 yilmaz 1.2 edm::Handle<pat::JetCollection> patjets;
225     if(usePat_)iEvent.getByLabel(jetTag_, patjets);
226    
227     edm::Handle<reco::JetView> jets;
228 yilmaz 1.1 iEvent.getByLabel(jetTag_, jets);
229    
230     // FILL JRA TREE
231    
232     jets_.b = b;
233     jets_.nref = 0;
234    
235     if(!isMC_){
236     fillL1Bits(iEvent);
237     fillHLTBits(iEvent);
238     }
239    
240     for(unsigned int j = 0 ; j < jets->size(); ++j){
241 yilmaz 1.2 const reco::Jet& jet = (*jets)[j];
242 yilmaz 1.1
243     //cout<<" jet pt "<<jet.pt()<<endl;
244     //if(jet.pt() < jetPtMin) continue;
245 yilmaz 1.2 if (useJEC_ && usePat_){
246     jets_.rawpt[jets_.nref]=(*patjets)[j].correctedJet("Uncorrected").pt();
247     }
248 yilmaz 1.1 jets_.jtpt[jets_.nref] = jet.pt();
249     jets_.jteta[jets_.nref] = jet.eta();
250     jets_.jtphi[jets_.nref] = jet.phi();
251     jets_.jty[jets_.nref] = jet.eta();
252 yilmaz 1.2 jets_.jtpu[jets_.nref] = jet.pileup();
253 yilmaz 1.1
254 yilmaz 1.2 if(isMC_ && usePat_){
255     const reco::GenJet* genjet = (*patjets)[j].genJet();
256     if(genjet){
257     jets_.refpt[jets_.nref] = genjet->pt();
258     jets_.refeta[jets_.nref] = genjet->eta();
259     jets_.refphi[jets_.nref] = genjet->phi();
260     jets_.refy[jets_.nref] = genjet->eta();
261     jets_.refdphijt[jets_.nref] = reco::deltaPhi(jet.phi(), genjet->phi());
262     jets_.refdrjt[jets_.nref] = reco::deltaR(jet.eta(),jet.phi(),genjet->eta(),genjet->phi());
263     }else{
264     jets_.refpt[jets_.nref] = -999.;
265     jets_.refeta[jets_.nref] = -999.;
266     jets_.refphi[jets_.nref] = -999.;
267     jets_.refy[jets_.nref] = -999.;
268     jets_.refdphijt[jets_.nref] = -999.;
269     jets_.refdrjt[jets_.nref] = -999.;
270     }
271 yilmaz 1.1
272 yilmaz 1.2 // matched partons
273     const reco::GenParticle * parton = (*patjets)[j].genParton();
274     if(parton){
275     jets_.refparton_pt[jets_.nref] = parton->pt();
276     jets_.refparton_flavor[jets_.nref] = parton->pdgId();
277 yilmaz 1.1 } else {
278     jets_.refparton_pt[jets_.nref] = -999;
279     jets_.refparton_flavor[jets_.nref] = -999;
280     }
281 yilmaz 1.2 }
282 yilmaz 1.1
283     jets_.nref++;
284    
285    
286     }
287    
288    
289     if(isMC_){
290    
291     edm::Handle<GenEventInfoProduct> hEventInfo;
292     iEvent.getByLabel(eventInfoTag_,hEventInfo);
293     //jets_.pthat = hEventInfo->binningValues()[0];
294    
295     // binning values and qscale appear to be equivalent, but binning values not always present
296     jets_.pthat = hEventInfo->qScale();
297    
298     edm::Handle<vector<reco::GenJet> >genjets;
299     iEvent.getByLabel(genjetTag_, genjets);
300    
301     jets_.ngen = 0;
302     for(unsigned int igen = 0 ; igen < genjets->size(); ++igen){
303     const reco::GenJet & genjet = (*genjets)[igen];
304    
305     float genjet_pt = genjet.pt();
306    
307     // threshold to reduce size of output in minbias PbPb
308     if(genjet_pt>20.){
309    
310     jets_.genpt[jets_.ngen] = genjet_pt;
311     jets_.geneta[jets_.ngen] = genjet.eta();
312     jets_.genphi[jets_.ngen] = genjet.phi();
313     jets_.geny[jets_.ngen] = genjet.eta();
314    
315    
316     // find matching patJet if there is one
317    
318     jets_.gendrjt[jets_.ngen] = -1.0;
319     jets_.genmatchindex[jets_.ngen] = -1;
320    
321    
322     for(unsigned int ijet = 0 ; ijet < jets->size(); ++ijet){
323     const pat::Jet& jet = (*jets)[ijet];
324    
325     if(jet.genJet()){
326     if(fabs(genjet.pt()-jet.genJet()->pt())<0.0001 &&
327     fabs(genjet.eta()-jet.genJet()->eta())<0.0001 &&
328     (fabs(genjet.phi()-jet.genJet()->phi())<0.0001 || fabs(genjet.phi()-jet.genJet()->phi() - 2.0*TMath::Pi()) < 0.0001 )){
329    
330     jets_.genmatchindex[jets_.ngen] = (int)ijet;
331     jets_.gendphijt[jets_.ngen] = reco::deltaPhi(jet.phi(),genjet.phi());
332     jets_.gendrjt[jets_.ngen] = reco::deltaR(jet,genjet);
333    
334     break;
335     }
336     }
337     }
338     jets_.ngen++;
339     }
340    
341     }
342    
343     }
344     t->Fill();
345     }
346    
347    
348    
349    
350     //--------------------------------------------------------------------------------------------------
351     void HiInclusiveJetAnalyzer::fillL1Bits(const edm::Event &iEvent)
352     {
353     edm::Handle< L1GlobalTriggerReadoutRecord > L1GlobalTrigger;
354    
355     iEvent.getByLabel(L1gtReadout_, L1GlobalTrigger);
356     const TechnicalTriggerWord& technicalTriggerWordBeforeMask = L1GlobalTrigger->technicalTriggerWord();
357    
358     for (int i=0; i<64;i++)
359     {
360     jets_.l1TBit[i] = technicalTriggerWordBeforeMask.at(i);
361     }
362     jets_.nL1TBit = 64;
363    
364     int ntrigs = L1GlobalTrigger->decisionWord().size();
365     jets_.nL1ABit = ntrigs;
366    
367     for (int i=0; i != ntrigs; i++) {
368     bool accept = L1GlobalTrigger->decisionWord()[i];
369     //jets_.l1ABit[i] = (accept == true)? 1:0;
370     if(accept== true){
371     jets_.l1ABit[i] = 1;
372     }
373     else{
374     jets_.l1ABit[i] = 0;
375     }
376    
377     }
378     }
379    
380     //--------------------------------------------------------------------------------------------------
381     void HiInclusiveJetAnalyzer::fillHLTBits(const edm::Event &iEvent)
382     {
383     // Fill HLT trigger bits.
384     Handle<TriggerResults> triggerResultsHLT;
385     getProduct(hltResName_, triggerResultsHLT, iEvent);
386    
387     const TriggerResults *hltResults = triggerResultsHLT.product();
388     const TriggerNames & triggerNames = iEvent.triggerNames(*hltResults);
389    
390     jets_.nHLTBit = hltTrgNames_.size();
391    
392     for(size_t i=0;i<hltTrgNames_.size();i++){
393    
394     for(size_t j=0;j<triggerNames.size();++j) {
395    
396     if(triggerNames.triggerName(j) == hltTrgNames_[i]){
397    
398     //cout <<"hltTrgNames_(i) "<<hltTrgNames_[i]<<endl;
399     //cout <<"triggerName(j) "<<triggerNames.triggerName(j)<<endl;
400     //cout<<" result "<<triggerResultsHLT->accept(j)<<endl;
401     jets_.hltBit[i] = triggerResultsHLT->accept(j);
402     }
403    
404     }
405     }
406     }
407    
408     //--------------------------------------------------------------------------------------------------
409     template <typename TYPE>
410     inline void HiInclusiveJetAnalyzer::getProduct(const std::string name, edm::Handle<TYPE> &prod,
411     const edm::Event &event) const
412     {
413     // Try to access data collection from EDM file. We check if we really get just one
414     // product with the given name. If not we throw an exception.
415    
416     event.getByLabel(edm::InputTag(name),prod);
417     if (!prod.isValid())
418     throw edm::Exception(edm::errors::Configuration, "HiInclusiveJetAnalyzer::GetProduct()\n")
419     << "Collection with label '" << name << "' is not valid" << std::endl;
420     }
421    
422     //--------------------------------------------------------------------------------------------------
423     template <typename TYPE>
424     inline bool HiInclusiveJetAnalyzer::getProductSafe(const std::string name, edm::Handle<TYPE> &prod,
425     const edm::Event &event) const
426     {
427     // Try to safely access data collection from EDM file. We check if we really get just one
428     // product with the given name. If not, we return false.
429    
430     if (name.size()==0)
431     return false;
432    
433     try {
434     event.getByLabel(edm::InputTag(name),prod);
435     if (!prod.isValid())
436     return false;
437     } catch (...) {
438     return false;
439     }
440     return true;
441     }
442    
443    
444    
445    
446     DEFINE_FWK_MODULE(HiInclusiveJetAnalyzer);