ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/CmsHi/JetAnalysis/src/HiInclusiveJetAnalyzer.cc
Revision: 1.7
Committed: Wed Nov 2 19:44:09 2011 UTC (13 years, 6 months ago) by yjlee
Content type: text/plain
Branch: MAIN
CVS Tags: hi44X_02
Changes since 1.6: +4 -3 lines
Log Message:
Remove L1 and HLT

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