ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/CmsHi/JetAnalysis/src/HiInclusiveJetAnalyzer.cc
Revision: 1.1
Committed: Tue Sep 20 19:45:06 2011 UTC (13 years, 7 months ago) by yilmaz
Content type: text/plain
Branch: MAIN
Log Message:
combining all analyzers

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