ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/CmsHi/JetAnalysis/src/HiInclusiveJetAnalyzer.cc
Revision: 1.2
Committed: Tue Sep 20 20:06:06 2011 UTC (13 years, 7 months ago) by yilmaz
Content type: text/plain
Branch: MAIN
Changes since 1.1: +34 -29 lines
Log Message:
made independent of jet format

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    
182     int bin = -1;
183     double hf = 0.;
184     double b = 999.;
185    
186     if(useCentrality_){
187     //if(!isMC_){
188     if(!centrality_) centrality_ = new CentralityProvider(iSetup);
189     centrality_->newEvent(iEvent,iSetup); // make sure you do this first in every event
190     //double c = centrality_->centralityValue();
191     const reco::Centrality *cent = centrality_->raw();
192    
193     hf = cent->EtHFhitSum();
194    
195     bin = centrality_->getBin();
196     b = centrality_->bMean();
197     //}
198     /*
199     else{
200    
201     TFile * centFile = new TFile("/net/hidsk0001/d00/scratch/mnguyen/CMSSW_3_9_1_patch1/src/macros/Hydjet_CentTable.root");
202    
203     edm::Handle<reco::Centrality> cent;
204     iEvent.getByLabel(edm::InputTag("hiCentrality"),cent);
205     //cout<<" grabbed centrality "<<endl;
206     CentralityBins::RunMap cmap = getCentralityFromFile(centFile, "makeCentralityTableTFile", "HFhitsHydjet_2760GeV", 1, 1);
207    
208     // Still getting cent from hits. When tower tables become available, we need to switch
209     hf = cent->EtHFhitSum();
210     //cout<<" hf "<<hf<<endl;
211     bin = cmap[run]->getBin(hf);
212     b = cmap[run]->bMeanOfBin(bin);
213     }
214     */
215     }
216    
217    
218    
219    
220     //double jetPtMin = 35.;
221    
222    
223     // loop the events
224    
225     jets_.bin = bin;
226     jets_.hf = hf;
227    
228    
229     if (useVtx_) {
230     edm::Handle<vector<reco::Vertex> >vertex;
231     iEvent.getByLabel(vtxTag_, vertex);
232    
233     if(vertex->size()>0) {
234     jets_.vx=vertex->begin()->x();
235     jets_.vy=vertex->begin()->y();
236     jets_.vz=vertex->begin()->z();
237     }
238     }
239    
240 yilmaz 1.2 edm::Handle<pat::JetCollection> patjets;
241     if(usePat_)iEvent.getByLabel(jetTag_, patjets);
242    
243     edm::Handle<reco::JetView> jets;
244 yilmaz 1.1 iEvent.getByLabel(jetTag_, jets);
245    
246     // FILL JRA TREE
247    
248     jets_.b = b;
249     jets_.nref = 0;
250    
251     if(!isMC_){
252     fillL1Bits(iEvent);
253     fillHLTBits(iEvent);
254     }
255    
256     for(unsigned int j = 0 ; j < jets->size(); ++j){
257 yilmaz 1.2 const reco::Jet& jet = (*jets)[j];
258 yilmaz 1.1
259     //cout<<" jet pt "<<jet.pt()<<endl;
260     //if(jet.pt() < jetPtMin) continue;
261 yilmaz 1.2 if (useJEC_ && usePat_){
262     jets_.rawpt[jets_.nref]=(*patjets)[j].correctedJet("Uncorrected").pt();
263     }
264 yilmaz 1.1 jets_.jtpt[jets_.nref] = jet.pt();
265     jets_.jteta[jets_.nref] = jet.eta();
266     jets_.jtphi[jets_.nref] = jet.phi();
267     jets_.jty[jets_.nref] = jet.eta();
268 yilmaz 1.2 jets_.jtpu[jets_.nref] = jet.pileup();
269 yilmaz 1.1
270 yilmaz 1.2 if(isMC_ && usePat_){
271     const reco::GenJet* genjet = (*patjets)[j].genJet();
272     if(genjet){
273     jets_.refpt[jets_.nref] = genjet->pt();
274     jets_.refeta[jets_.nref] = genjet->eta();
275     jets_.refphi[jets_.nref] = genjet->phi();
276     jets_.refy[jets_.nref] = genjet->eta();
277     jets_.refdphijt[jets_.nref] = reco::deltaPhi(jet.phi(), genjet->phi());
278     jets_.refdrjt[jets_.nref] = reco::deltaR(jet.eta(),jet.phi(),genjet->eta(),genjet->phi());
279     }else{
280     jets_.refpt[jets_.nref] = -999.;
281     jets_.refeta[jets_.nref] = -999.;
282     jets_.refphi[jets_.nref] = -999.;
283     jets_.refy[jets_.nref] = -999.;
284     jets_.refdphijt[jets_.nref] = -999.;
285     jets_.refdrjt[jets_.nref] = -999.;
286     }
287 yilmaz 1.1
288 yilmaz 1.2 // matched partons
289     const reco::GenParticle * parton = (*patjets)[j].genParton();
290     if(parton){
291     jets_.refparton_pt[jets_.nref] = parton->pt();
292     jets_.refparton_flavor[jets_.nref] = parton->pdgId();
293 yilmaz 1.1 } else {
294     jets_.refparton_pt[jets_.nref] = -999;
295     jets_.refparton_flavor[jets_.nref] = -999;
296     }
297 yilmaz 1.2 }
298 yilmaz 1.1
299     jets_.nref++;
300    
301    
302     }
303    
304    
305     if(isMC_){
306    
307     edm::Handle<GenEventInfoProduct> hEventInfo;
308     iEvent.getByLabel(eventInfoTag_,hEventInfo);
309     //jets_.pthat = hEventInfo->binningValues()[0];
310    
311     // binning values and qscale appear to be equivalent, but binning values not always present
312     jets_.pthat = hEventInfo->qScale();
313    
314     edm::Handle<vector<reco::GenJet> >genjets;
315     iEvent.getByLabel(genjetTag_, genjets);
316    
317     jets_.ngen = 0;
318     for(unsigned int igen = 0 ; igen < genjets->size(); ++igen){
319     const reco::GenJet & genjet = (*genjets)[igen];
320    
321     float genjet_pt = genjet.pt();
322    
323     // threshold to reduce size of output in minbias PbPb
324     if(genjet_pt>20.){
325    
326     jets_.genpt[jets_.ngen] = genjet_pt;
327     jets_.geneta[jets_.ngen] = genjet.eta();
328     jets_.genphi[jets_.ngen] = genjet.phi();
329     jets_.geny[jets_.ngen] = genjet.eta();
330    
331    
332     // find matching patJet if there is one
333    
334     jets_.gendrjt[jets_.ngen] = -1.0;
335     jets_.genmatchindex[jets_.ngen] = -1;
336    
337    
338     for(unsigned int ijet = 0 ; ijet < jets->size(); ++ijet){
339     const pat::Jet& jet = (*jets)[ijet];
340    
341     if(jet.genJet()){
342     if(fabs(genjet.pt()-jet.genJet()->pt())<0.0001 &&
343     fabs(genjet.eta()-jet.genJet()->eta())<0.0001 &&
344     (fabs(genjet.phi()-jet.genJet()->phi())<0.0001 || fabs(genjet.phi()-jet.genJet()->phi() - 2.0*TMath::Pi()) < 0.0001 )){
345    
346     jets_.genmatchindex[jets_.ngen] = (int)ijet;
347     jets_.gendphijt[jets_.ngen] = reco::deltaPhi(jet.phi(),genjet.phi());
348     jets_.gendrjt[jets_.ngen] = reco::deltaR(jet,genjet);
349    
350     break;
351     }
352     }
353     }
354     jets_.ngen++;
355     }
356    
357     }
358    
359     }
360     t->Fill();
361     }
362    
363    
364    
365    
366     //--------------------------------------------------------------------------------------------------
367     void HiInclusiveJetAnalyzer::fillL1Bits(const edm::Event &iEvent)
368     {
369     edm::Handle< L1GlobalTriggerReadoutRecord > L1GlobalTrigger;
370    
371     iEvent.getByLabel(L1gtReadout_, L1GlobalTrigger);
372     const TechnicalTriggerWord& technicalTriggerWordBeforeMask = L1GlobalTrigger->technicalTriggerWord();
373    
374     for (int i=0; i<64;i++)
375     {
376     jets_.l1TBit[i] = technicalTriggerWordBeforeMask.at(i);
377     }
378     jets_.nL1TBit = 64;
379    
380     int ntrigs = L1GlobalTrigger->decisionWord().size();
381     jets_.nL1ABit = ntrigs;
382    
383     for (int i=0; i != ntrigs; i++) {
384     bool accept = L1GlobalTrigger->decisionWord()[i];
385     //jets_.l1ABit[i] = (accept == true)? 1:0;
386     if(accept== true){
387     jets_.l1ABit[i] = 1;
388     }
389     else{
390     jets_.l1ABit[i] = 0;
391     }
392    
393     }
394     }
395    
396     //--------------------------------------------------------------------------------------------------
397     void HiInclusiveJetAnalyzer::fillHLTBits(const edm::Event &iEvent)
398     {
399     // Fill HLT trigger bits.
400     Handle<TriggerResults> triggerResultsHLT;
401     getProduct(hltResName_, triggerResultsHLT, iEvent);
402    
403     const TriggerResults *hltResults = triggerResultsHLT.product();
404     const TriggerNames & triggerNames = iEvent.triggerNames(*hltResults);
405    
406     jets_.nHLTBit = hltTrgNames_.size();
407    
408     for(size_t i=0;i<hltTrgNames_.size();i++){
409    
410     for(size_t j=0;j<triggerNames.size();++j) {
411    
412     if(triggerNames.triggerName(j) == hltTrgNames_[i]){
413    
414     //cout <<"hltTrgNames_(i) "<<hltTrgNames_[i]<<endl;
415     //cout <<"triggerName(j) "<<triggerNames.triggerName(j)<<endl;
416     //cout<<" result "<<triggerResultsHLT->accept(j)<<endl;
417     jets_.hltBit[i] = triggerResultsHLT->accept(j);
418     }
419    
420     }
421     }
422     }
423    
424     //--------------------------------------------------------------------------------------------------
425     template <typename TYPE>
426     inline void HiInclusiveJetAnalyzer::getProduct(const std::string name, edm::Handle<TYPE> &prod,
427     const edm::Event &event) const
428     {
429     // Try to access data collection from EDM file. We check if we really get just one
430     // product with the given name. If not we throw an exception.
431    
432     event.getByLabel(edm::InputTag(name),prod);
433     if (!prod.isValid())
434     throw edm::Exception(edm::errors::Configuration, "HiInclusiveJetAnalyzer::GetProduct()\n")
435     << "Collection with label '" << name << "' is not valid" << std::endl;
436     }
437    
438     //--------------------------------------------------------------------------------------------------
439     template <typename TYPE>
440     inline bool HiInclusiveJetAnalyzer::getProductSafe(const std::string name, edm::Handle<TYPE> &prod,
441     const edm::Event &event) const
442     {
443     // Try to safely access data collection from EDM file. We check if we really get just one
444     // product with the given name. If not, we return false.
445    
446     if (name.size()==0)
447     return false;
448    
449     try {
450     event.getByLabel(edm::InputTag(name),prod);
451     if (!prod.isValid())
452     return false;
453     } catch (...) {
454     return false;
455     }
456     return true;
457     }
458    
459    
460    
461    
462     DEFINE_FWK_MODULE(HiInclusiveJetAnalyzer);