ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/CmsHi/JetAnalysis/src/HiInclusiveJetAnalyzer.cc
Revision: 1.3
Committed: Tue Sep 20 22:24:42 2011 UTC (13 years, 7 months ago) by yilmaz
Content type: text/plain
Branch: MAIN
CVS Tags: hi413_03, hi413_02, hi441_1, hi441_0, hi413_11, hi413_10, hi413_09
Changes since 1.2: +0 -20 lines
Log Message:
option for neutral components only

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     useVtx_ = iConfig.getUntrackedParameter<bool>("useVtx",true);
61     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     // TTree* t= new TTree("t","Jet Response Analyzer");
107     t->Branch("run",&jets_.run,"run/I");
108     t->Branch("evt",&jets_.evt,"evt/I");
109     t->Branch("lumi",&jets_.lumi,"lumi/I");
110     t->Branch("b",&jets_.b,"b/F");
111     t->Branch("vx",&jets_.vx,"vx/F");
112     t->Branch("vy",&jets_.vy,"vy/F");
113     t->Branch("vz",&jets_.vz,"vz/F");
114     t->Branch("hf",&jets_.hf,"hf/F");
115     t->Branch("nref",&jets_.nref,"nref/I");
116     t->Branch("bin",&jets_.bin,"bin/I");
117     t->Branch("rawpt",jets_.rawpt,"rawpt[nref]/F");
118     t->Branch("jtpt",jets_.jtpt,"jtpt[nref]/F");
119     t->Branch("jteta",jets_.jteta,"jteta[nref]/F");
120     t->Branch("jty",jets_.jty,"jty[nref]/F");
121     t->Branch("jtphi",jets_.jtphi,"jtphi[nref]/F");
122 yilmaz 1.2 t->Branch("jtpu",jets_.jtpu,"jtpu[nref]/F");
123 yilmaz 1.1
124     if(isMC_){
125     t->Branch("pthat",&jets_.pthat,"pthat/F");
126    
127     // Only matched gen jets
128     t->Branch("refpt",jets_.refpt,"refpt[nref]/F");
129     t->Branch("refeta",jets_.refeta,"refeta[nref]/F");
130     t->Branch("refy",jets_.refy,"refy[nref]/F");
131     t->Branch("refphi",jets_.refphi,"refphi[nref]/F");
132     t->Branch("refdphijt",jets_.refdphijt,"refdphijt[nref]/F");
133     t->Branch("refdrjt",jets_.refdrjt,"refdrjt[nref]/F");
134     // matched parton
135     t->Branch("refparton_pt",jets_.refparton_pt,"refparton_pt[nref]/F");
136     t->Branch("refparton_flavor",jets_.refparton_flavor,"refparton_flavor[nref]/F");
137    
138     // For all gen jets, matched or unmatched
139     t->Branch("ngen",&jets_.ngen,"ngen/I");
140     t->Branch("genmatchindex",jets_.genmatchindex,"genmatchindex[ngen]/I");
141     t->Branch("genpt",jets_.genpt,"genpt[ngen]/F");
142     t->Branch("geneta",jets_.geneta,"geneta[ngen]/F");
143     t->Branch("geny",jets_.geny,"geny[ngen]/F");
144     t->Branch("genphi",jets_.genphi,"genphi[ngen]/F");
145     t->Branch("gendphijt",jets_.gendphijt,"gendphijt[ngen]/F");
146     t->Branch("gendrjt",jets_.gendrjt,"gendrjt[ngen]/F");
147     }
148    
149     if(!isMC_){
150     t->Branch("nL1TBit",&jets_.nL1TBit,"nL1TBit/I");
151     t->Branch("l1TBit",jets_.l1TBit,"l1TBit[nL1TBit]/O");
152    
153     t->Branch("nL1ABit",&jets_.nL1ABit,"nL1ABit/I");
154     t->Branch("l1ABit",jets_.l1ABit,"l1ABit[nL1ABit]/O");
155    
156     t->Branch("nHLTBit",&jets_.nHLTBit,"nHLTBit/I");
157     t->Branch("hltBit",jets_.hltBit,"hltBit[nHLTBit]/O");
158    
159     }
160    
161     TH1D::SetDefaultSumw2();
162    
163    
164     }
165    
166    
167     void
168     HiInclusiveJetAnalyzer::analyze(const Event& iEvent,
169     const EventSetup& iSetup) {
170    
171     int event = iEvent.id().event();
172     int run = iEvent.id().run();
173     int lumi = iEvent.id().luminosityBlock();
174    
175     jets_.run = run;
176     jets_.evt = event;
177     jets_.lumi = lumi;
178    
179     LogDebug("HiInclusiveJetAnalyzer")<<"START event: "<<event<<" in run "<<run<<endl;
180    
181     int bin = -1;
182     double hf = 0.;
183     double b = 999.;
184    
185     if(useCentrality_){
186     if(!centrality_) centrality_ = new CentralityProvider(iSetup);
187     centrality_->newEvent(iEvent,iSetup); // make sure you do this first in every event
188     //double c = centrality_->centralityValue();
189     const reco::Centrality *cent = centrality_->raw();
190    
191     hf = cent->EtHFhitSum();
192    
193     bin = centrality_->getBin();
194     b = centrality_->bMean();
195     }
196    
197    
198    
199    
200     //double jetPtMin = 35.;
201    
202    
203     // loop the events
204    
205     jets_.bin = bin;
206     jets_.hf = hf;
207    
208    
209     if (useVtx_) {
210     edm::Handle<vector<reco::Vertex> >vertex;
211     iEvent.getByLabel(vtxTag_, vertex);
212    
213     if(vertex->size()>0) {
214     jets_.vx=vertex->begin()->x();
215     jets_.vy=vertex->begin()->y();
216     jets_.vz=vertex->begin()->z();
217     }
218     }
219    
220 yilmaz 1.2 edm::Handle<pat::JetCollection> patjets;
221     if(usePat_)iEvent.getByLabel(jetTag_, patjets);
222    
223     edm::Handle<reco::JetView> jets;
224 yilmaz 1.1 iEvent.getByLabel(jetTag_, jets);
225    
226     // FILL JRA TREE
227    
228     jets_.b = b;
229     jets_.nref = 0;
230    
231     if(!isMC_){
232     fillL1Bits(iEvent);
233     fillHLTBits(iEvent);
234     }
235    
236     for(unsigned int j = 0 ; j < jets->size(); ++j){
237 yilmaz 1.2 const reco::Jet& jet = (*jets)[j];
238 yilmaz 1.1
239     //cout<<" jet pt "<<jet.pt()<<endl;
240     //if(jet.pt() < jetPtMin) continue;
241 yilmaz 1.2 if (useJEC_ && usePat_){
242     jets_.rawpt[jets_.nref]=(*patjets)[j].correctedJet("Uncorrected").pt();
243     }
244 yilmaz 1.1 jets_.jtpt[jets_.nref] = jet.pt();
245     jets_.jteta[jets_.nref] = jet.eta();
246     jets_.jtphi[jets_.nref] = jet.phi();
247     jets_.jty[jets_.nref] = jet.eta();
248 yilmaz 1.2 jets_.jtpu[jets_.nref] = jet.pileup();
249 yilmaz 1.1
250 yilmaz 1.2 if(isMC_ && usePat_){
251     const reco::GenJet* genjet = (*patjets)[j].genJet();
252     if(genjet){
253     jets_.refpt[jets_.nref] = genjet->pt();
254     jets_.refeta[jets_.nref] = genjet->eta();
255     jets_.refphi[jets_.nref] = genjet->phi();
256     jets_.refy[jets_.nref] = genjet->eta();
257     jets_.refdphijt[jets_.nref] = reco::deltaPhi(jet.phi(), genjet->phi());
258     jets_.refdrjt[jets_.nref] = reco::deltaR(jet.eta(),jet.phi(),genjet->eta(),genjet->phi());
259     }else{
260     jets_.refpt[jets_.nref] = -999.;
261     jets_.refeta[jets_.nref] = -999.;
262     jets_.refphi[jets_.nref] = -999.;
263     jets_.refy[jets_.nref] = -999.;
264     jets_.refdphijt[jets_.nref] = -999.;
265     jets_.refdrjt[jets_.nref] = -999.;
266     }
267 yilmaz 1.1
268 yilmaz 1.2 // matched partons
269     const reco::GenParticle * parton = (*patjets)[j].genParton();
270     if(parton){
271     jets_.refparton_pt[jets_.nref] = parton->pt();
272     jets_.refparton_flavor[jets_.nref] = parton->pdgId();
273 yilmaz 1.1 } else {
274     jets_.refparton_pt[jets_.nref] = -999;
275     jets_.refparton_flavor[jets_.nref] = -999;
276     }
277 yilmaz 1.2 }
278 yilmaz 1.1
279     jets_.nref++;
280    
281    
282     }
283    
284    
285     if(isMC_){
286    
287     edm::Handle<GenEventInfoProduct> hEventInfo;
288     iEvent.getByLabel(eventInfoTag_,hEventInfo);
289     //jets_.pthat = hEventInfo->binningValues()[0];
290    
291     // binning values and qscale appear to be equivalent, but binning values not always present
292     jets_.pthat = hEventInfo->qScale();
293    
294     edm::Handle<vector<reco::GenJet> >genjets;
295     iEvent.getByLabel(genjetTag_, genjets);
296    
297     jets_.ngen = 0;
298     for(unsigned int igen = 0 ; igen < genjets->size(); ++igen){
299     const reco::GenJet & genjet = (*genjets)[igen];
300    
301     float genjet_pt = genjet.pt();
302    
303     // threshold to reduce size of output in minbias PbPb
304     if(genjet_pt>20.){
305    
306     jets_.genpt[jets_.ngen] = genjet_pt;
307     jets_.geneta[jets_.ngen] = genjet.eta();
308     jets_.genphi[jets_.ngen] = genjet.phi();
309     jets_.geny[jets_.ngen] = genjet.eta();
310    
311    
312     // find matching patJet if there is one
313    
314     jets_.gendrjt[jets_.ngen] = -1.0;
315     jets_.genmatchindex[jets_.ngen] = -1;
316    
317    
318     for(unsigned int ijet = 0 ; ijet < jets->size(); ++ijet){
319     const pat::Jet& jet = (*jets)[ijet];
320    
321     if(jet.genJet()){
322     if(fabs(genjet.pt()-jet.genJet()->pt())<0.0001 &&
323     fabs(genjet.eta()-jet.genJet()->eta())<0.0001 &&
324     (fabs(genjet.phi()-jet.genJet()->phi())<0.0001 || fabs(genjet.phi()-jet.genJet()->phi() - 2.0*TMath::Pi()) < 0.0001 )){
325    
326     jets_.genmatchindex[jets_.ngen] = (int)ijet;
327     jets_.gendphijt[jets_.ngen] = reco::deltaPhi(jet.phi(),genjet.phi());
328     jets_.gendrjt[jets_.ngen] = reco::deltaR(jet,genjet);
329    
330     break;
331     }
332     }
333     }
334     jets_.ngen++;
335     }
336    
337     }
338    
339     }
340     t->Fill();
341     }
342    
343    
344    
345    
346     //--------------------------------------------------------------------------------------------------
347     void HiInclusiveJetAnalyzer::fillL1Bits(const edm::Event &iEvent)
348     {
349     edm::Handle< L1GlobalTriggerReadoutRecord > L1GlobalTrigger;
350    
351     iEvent.getByLabel(L1gtReadout_, L1GlobalTrigger);
352     const TechnicalTriggerWord& technicalTriggerWordBeforeMask = L1GlobalTrigger->technicalTriggerWord();
353    
354     for (int i=0; i<64;i++)
355     {
356     jets_.l1TBit[i] = technicalTriggerWordBeforeMask.at(i);
357     }
358     jets_.nL1TBit = 64;
359    
360     int ntrigs = L1GlobalTrigger->decisionWord().size();
361     jets_.nL1ABit = ntrigs;
362    
363     for (int i=0; i != ntrigs; i++) {
364     bool accept = L1GlobalTrigger->decisionWord()[i];
365     //jets_.l1ABit[i] = (accept == true)? 1:0;
366     if(accept== true){
367     jets_.l1ABit[i] = 1;
368     }
369     else{
370     jets_.l1ABit[i] = 0;
371     }
372    
373     }
374     }
375    
376     //--------------------------------------------------------------------------------------------------
377     void HiInclusiveJetAnalyzer::fillHLTBits(const edm::Event &iEvent)
378     {
379     // Fill HLT trigger bits.
380     Handle<TriggerResults> triggerResultsHLT;
381     getProduct(hltResName_, triggerResultsHLT, iEvent);
382    
383     const TriggerResults *hltResults = triggerResultsHLT.product();
384     const TriggerNames & triggerNames = iEvent.triggerNames(*hltResults);
385    
386     jets_.nHLTBit = hltTrgNames_.size();
387    
388     for(size_t i=0;i<hltTrgNames_.size();i++){
389    
390     for(size_t j=0;j<triggerNames.size();++j) {
391    
392     if(triggerNames.triggerName(j) == hltTrgNames_[i]){
393    
394     //cout <<"hltTrgNames_(i) "<<hltTrgNames_[i]<<endl;
395     //cout <<"triggerName(j) "<<triggerNames.triggerName(j)<<endl;
396     //cout<<" result "<<triggerResultsHLT->accept(j)<<endl;
397     jets_.hltBit[i] = triggerResultsHLT->accept(j);
398     }
399    
400     }
401     }
402     }
403    
404     //--------------------------------------------------------------------------------------------------
405     template <typename TYPE>
406     inline void HiInclusiveJetAnalyzer::getProduct(const std::string name, edm::Handle<TYPE> &prod,
407     const edm::Event &event) const
408     {
409     // Try to access data collection from EDM file. We check if we really get just one
410     // product with the given name. If not we throw an exception.
411    
412     event.getByLabel(edm::InputTag(name),prod);
413     if (!prod.isValid())
414     throw edm::Exception(edm::errors::Configuration, "HiInclusiveJetAnalyzer::GetProduct()\n")
415     << "Collection with label '" << name << "' is not valid" << std::endl;
416     }
417    
418     //--------------------------------------------------------------------------------------------------
419     template <typename TYPE>
420     inline bool HiInclusiveJetAnalyzer::getProductSafe(const std::string name, edm::Handle<TYPE> &prod,
421     const edm::Event &event) const
422     {
423     // Try to safely access data collection from EDM file. We check if we really get just one
424     // product with the given name. If not, we return false.
425    
426     if (name.size()==0)
427     return false;
428    
429     try {
430     event.getByLabel(edm::InputTag(name),prod);
431     if (!prod.isValid())
432     return false;
433     } catch (...) {
434     return false;
435     }
436     return true;
437     }
438    
439    
440    
441    
442     DEFINE_FWK_MODULE(HiInclusiveJetAnalyzer);