ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Betchart/TopRefTuple/interface/Tuple_Jet.h
Revision: 1.8
Committed: Fri Jan 25 00:02:23 2013 UTC (12 years, 3 months ago) by bbetchar
Content type: text/plain
Branch: MAIN
CVS Tags: V00-03-02, V00-03-01, HEAD
Changes since 1.7: +25 -0 lines
Log Message:
keep jet- Nelectron, MuonsChargeSum, ElectronsChargeSum

File Contents

# User Rev Content
1 bbetchar 1.1 #ifndef TUPLE_JET
2     #define TUPLE_JET
3    
4 bbetchar 1.2 #include "boost/foreach.hpp"
5 bbetchar 1.3 #include "TopQuarkAnalysis/TopRefTuple/interface/fTypes.h"
6 bbetchar 1.1 #include "FWCore/Framework/interface/EDProducer.h"
7     #include "FWCore/Framework/interface/Frameworkfwd.h"
8     #include "FWCore/Framework/interface/Event.h"
9    
10     #include "DataFormats/PatCandidates/interface/Jet.h"
11    
12     #include "PhysicsTools/SelectorUtils/interface/PFJetIDSelectionFunctor.h"
13    
14 bbetchar 1.4 #include "CondFormats/JetMETObjects/interface/JetResolution.h"
15 bbetchar 1.1 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
16     #include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
17     #include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
18 bbetchar 1.4 #include "TH1D.h"
19    
20    
21     JetCorrectionUncertainty* jetCorrUnc(const edm::EventSetup& setup, const std::string& jecRecord) {
22     edm::ESHandle<JetCorrectorParametersCollection> JetCorParColl;
23     setup.get<JetCorrectionsRecord>().get(jecRecord,JetCorParColl);
24     return new JetCorrectionUncertainty((*JetCorParColl)["Uncertainty"]);
25     }
26    
27     float uncFunc(JetCorrectionUncertainty* jCU, const reco::Candidate::LorentzVector& jet) {
28     jCU->setJetEta(jet.eta());
29 bbetchar 1.5 jCU->setJetPt(jet.pt());
30     try {return jCU->getUncertainty(true);}
31 bbetchar 1.4 catch (...) {
32 bbetchar 1.5 std::cout << "Tuple_Jet.h: JetCorrectionUncertainty::getUncertainty"
33     << "threw exception on jet with pt( " << jet.pt() << " ) "
34     << "and eta ( " << jet.eta() << " )." << std::endl;
35 bbetchar 1.4 return 0.0;
36     }
37     }
38 bbetchar 1.1
39 bbetchar 1.5
40    
41 bbetchar 1.1 template< typename T >
42     class Tuple_Jet : public edm::EDProducer {
43     public:
44     explicit Tuple_Jet(const edm::ParameterSet&);
45 bbetchar 1.4 typedef fTypes::dPolarLorentzV LorentzV;
46 bbetchar 1.5 typedef edm::Handle<edm::View<T> > Handle_t;
47 bbetchar 1.1 private:
48     void produce(edm::Event&, const edm::EventSetup& );
49 bbetchar 1.5 void initPF(); void producePF(edm::Event&, const Handle_t&, const Handle_t&);
50     void initGen(); void produceGen(edm::Event&, const Handle_t&, const Handle_t&);
51 bbetchar 1.1
52 bbetchar 1.4 const edm::InputTag jetsTag, allJetsTag;
53 bbetchar 1.2 const std::string prefix,jecRecord;
54     const std::vector<std::string> btagNames;
55 bbetchar 1.4 const bool pfInfo, genInfo;
56     JetResolution jer;
57 bbetchar 1.5 TH1D dataMcResRatio;
58 bbetchar 1.1 };
59    
60     template<class T> Tuple_Jet<T>::
61     Tuple_Jet(const edm::ParameterSet& cfg) :
62 bbetchar 1.2 jetsTag(cfg.getParameter<edm::InputTag>("jetsTag")),
63 bbetchar 1.4 allJetsTag(cfg.getParameter<edm::InputTag>("allJetsTag")),
64 bbetchar 1.2 prefix(cfg.getParameter<std::string>("prefix")),
65     jecRecord(cfg.getParameter<std::string>("jecRecord")),
66     btagNames(cfg.getParameter<std::vector<std::string> >("bTags")),
67 bbetchar 1.4 pfInfo(cfg.getParameter<bool>("pfInfo")),
68     genInfo(cfg.getParameter<bool>("genInfo")),
69     jer(cfg.getParameter<std::string>("jetResolutionFile")),
70 bbetchar 1.5 dataMcResRatio("dataMcResRatio",
71     "",
72     cfg.getParameter<std::vector<double> >("resolutionRatioBins").size()-1,
73     &cfg.getParameter<std::vector<double> >("resolutionRatioBins").at(0) )
74 bbetchar 1.1 {
75 bbetchar 1.4 produces <std::vector<LorentzV> > ( prefix + "P4" );
76     produces <std::vector<float> > ( prefix + "JecFactor" );
77     produces <std::vector<float> > ( prefix + "JecUnc" );
78     produces <std::vector<float> > ( prefix + "Area" );
79     produces <std::vector<float> > ( prefix + "Resolution" );
80     produces <std::vector<float> > ( prefix + "Charge" );
81     produces <std::vector<float> > ( prefix + "Eta2Moment" );
82     produces <std::vector<float> > ( prefix + "Phi2Moment" );
83 bbetchar 1.5 BOOST_FOREACH(const std::string& btag, btagNames)
84     produces <std::vector<float> > (prefix + btag);
85 bbetchar 1.4
86     if(pfInfo) initPF();
87 bbetchar 1.5 if(genInfo) {
88     initGen();
89     for(int bin=1; bin<=dataMcResRatio.GetNbinsX(); ++bin) {
90     dataMcResRatio.SetBinContent( bin, cfg.getParameter<std::vector<double> >("resolutionRatio").at(bin-1) );
91     dataMcResRatio.SetBinError( bin, cfg.getParameter<std::vector<double> >("resolutionRatioErr").at(bin-1) );
92     }
93     dataMcResRatio.SetBinContent( 1 + dataMcResRatio.GetNbinsX(), 1.0);
94     dataMcResRatio.SetBinError( 1 + dataMcResRatio.GetNbinsX(), 0.0);
95     }
96 bbetchar 1.1 }
97    
98     template< typename T >
99     void Tuple_Jet<T>::
100     produce(edm::Event& evt, const edm::EventSetup& setup) {
101 bbetchar 1.5 Handle_t jets; evt.getByLabel(jetsTag, jets);
102     Handle_t all; evt.getByLabel(allJetsTag, all);
103 bbetchar 1.1
104     std::auto_ptr<std::vector<LorentzV> > p4 ( new std::vector<LorentzV>() ) ;
105 bbetchar 1.4 std::auto_ptr<std::vector<float> > jecFactor( new std::vector<float>() ) ;
106 bbetchar 1.1 std::auto_ptr<std::vector<float> > jecUnc ( new std::vector<float>() ) ;
107 bbetchar 1.4 std::auto_ptr<std::vector<float> > reso ( new std::vector<float>() ) ;
108     std::auto_ptr<std::vector<float> > area ( new std::vector<float>() ) ;
109     std::auto_ptr<std::vector<float> > charge ( new std::vector<float>() ) ;
110 bbetchar 1.1 std::auto_ptr<std::vector<float> > eta2mom ( new std::vector<float>() ) ;
111     std::auto_ptr<std::vector<float> > phi2mom ( new std::vector<float>() ) ;
112 bbetchar 1.5 std::map<std::string, std::vector<float>* > btags;
113     BOOST_FOREACH(const std::string& btag, btagNames)
114     btags[btag] = new std::vector<float>();
115 bbetchar 1.1
116 bbetchar 1.2 if(jets.isValid()) {
117     JetCorrectionUncertainty* jCU = jetCorrUnc(setup, jecRecord);
118     for(typename edm::View<T>::const_iterator jet = jets->begin(); jet!=jets->end(); jet++) {
119 bbetchar 1.3 p4->push_back(LorentzV(jet->pt(),jet->eta(),jet->phi(),jet->mass()));
120 bbetchar 1.4 jecFactor->push_back( jet->jecSetsAvailable() ? jet->energy() / jet->correctedJet("Uncorrected").energy() : 1.0 );
121     jecUnc->push_back(uncFunc(jCU, jet->p4()));
122     area->push_back(jet->jetArea());
123     charge->push_back(jet->jetCharge());
124 bbetchar 1.2 eta2mom->push_back(jet->etaetaMoment());
125     phi2mom->push_back(jet->phiphiMoment());
126 bbetchar 1.5
127     BOOST_FOREACH(const std::string& btag, btagNames)
128     btags[btag]->push_back(jet->bDiscriminator(btag+"BJetTags"));
129 bbetchar 1.6
130     TF1 * jerEta = jer.parameterEta("sigma",jet->eta());
131     reso->push_back( jerEta->Eval( jet->pt() ) );
132     delete jerEta;
133 bbetchar 1.2 }
134     delete jCU;
135 bbetchar 1.1 }
136    
137 bbetchar 1.4 evt.put( p4, prefix + "P4" );
138     evt.put(jecFactor, prefix + "JecFactor" );
139     evt.put( jecUnc, prefix + "JecUnc" );
140     evt.put( reso, prefix + "Resolution" );
141     evt.put( area, prefix + "Area" );
142     evt.put( charge, prefix + "Charge" );
143     evt.put( eta2mom, prefix + "Eta2Moment" );
144     evt.put( phi2mom, prefix + "Phi2Moment" );
145 bbetchar 1.5 BOOST_FOREACH(const std::string& btag, btagNames)
146     evt.put( std::auto_ptr<std::vector<float> >(btags[btag]), prefix + btag );
147    
148     if(pfInfo) producePF(evt,jets, all);
149     if(genInfo) produceGen(evt,jets, all);
150 bbetchar 1.1 }
151    
152 bbetchar 1.2 template<class T> void Tuple_Jet<T>::
153 bbetchar 1.4 initPF() {
154 bbetchar 1.2 produces <std::vector<float> > (prefix + "FchargedHad" );
155     produces <std::vector<float> > (prefix + "FneutralHad" );
156     produces <std::vector<float> > (prefix + "FchargedEm" );
157     produces <std::vector<float> > (prefix + "FneutralEm" );
158     produces <std::vector<float> > (prefix + "FchargedMu" );
159 bbetchar 1.1
160 bbetchar 1.2 produces <std::vector<unsigned> > (prefix + "Ncharged" );
161     produces <std::vector<unsigned> > (prefix + "Nneutral" );
162     produces <std::vector<unsigned> > (prefix + "Nmuon" );
163 bbetchar 1.8 produces <std::vector<unsigned> > (prefix + "Nelectron" );
164 bbetchar 1.4 produces <std::vector<unsigned> > (prefix + "Ndaughters" );
165 bbetchar 1.1
166 bbetchar 1.2 produces <std::vector<bool> > ( prefix + "PFJetIDloose" );
167     produces <std::vector<bool> > ( prefix + "PFJetIDtight" );
168 bbetchar 1.5 produces <float> ( prefix + "FailedPtMax");
169 bbetchar 1.1
170 bbetchar 1.8 produces <std::vector<int> > (prefix + "MuonsChargeSum" );
171     produces <std::vector<int> > (prefix + "ElectronsChargeSum" );
172 bbetchar 1.1 }
173    
174     template<class T> void Tuple_Jet<T>::
175 bbetchar 1.5 producePF(edm::Event& evt, const Handle_t& jets, const Handle_t& all) {
176 bbetchar 1.1 std::auto_ptr<std::vector<float> > FchargedHad( new std::vector<float>() );
177     std::auto_ptr<std::vector<float> > FneutralHad( new std::vector<float>() );
178     std::auto_ptr<std::vector<float> > FchargedEm( new std::vector<float>() );
179     std::auto_ptr<std::vector<float> > FneutralEm( new std::vector<float>() );
180     std::auto_ptr<std::vector<float> > FchargedMu( new std::vector<float>() );
181    
182     std::auto_ptr<std::vector<unsigned> > Ncharged( new std::vector<unsigned>() );
183     std::auto_ptr<std::vector<unsigned> > Nneutral( new std::vector<unsigned>() );
184     std::auto_ptr<std::vector<unsigned> > Nmuon( new std::vector<unsigned>() );
185 bbetchar 1.8 std::auto_ptr<std::vector<unsigned> > Nelectron( new std::vector<unsigned>() );
186 bbetchar 1.4 std::auto_ptr<std::vector<unsigned> > Ndaughters( new std::vector<unsigned>() );
187 bbetchar 1.1
188 bbetchar 1.2 std::auto_ptr<std::vector<bool> > pfjetidloose ( new std::vector<bool>() ) ;
189     std::auto_ptr<std::vector<bool> > pfjetidtight ( new std::vector<bool>() ) ;
190 bbetchar 1.5 std::auto_ptr<float> failedPt ( new float(-1) ) ;
191 bbetchar 1.1
192 bbetchar 1.8 std::auto_ptr<std::vector<int> > muonsChargeSum ( new std::vector<int>() ) ;
193     std::auto_ptr<std::vector<int> > electronsChargeSum ( new std::vector<int>() ) ;
194    
195 bbetchar 1.1 PFJetIDSelectionFunctor
196     pfLooseJetID(PFJetIDSelectionFunctor::FIRSTDATA, PFJetIDSelectionFunctor::LOOSE),
197     pfTightJetID(PFJetIDSelectionFunctor::FIRSTDATA, PFJetIDSelectionFunctor::TIGHT);
198    
199 bbetchar 1.2 if(jets.isValid()) {
200     for(typename edm::View<T>::const_iterator jet = jets->begin(); jet!=jets->end(); jet++) {
201     pat::strbitset
202     passLooseCuts( pfLooseJetID .getBitTemplate() ),
203     passTightCuts( pfTightJetID .getBitTemplate() );
204    
205     FchargedHad->push_back( jet->chargedHadronEnergyFraction() );
206     FneutralHad->push_back( jet->neutralHadronEnergyFraction() );
207     FchargedEm->push_back( jet->chargedEmEnergyFraction() );
208     FneutralEm->push_back( jet->neutralEmEnergyFraction() );
209     FchargedMu->push_back( jet->chargedMuEnergyFraction() );
210     Ncharged->push_back( (unsigned) jet->chargedMultiplicity() );
211     Nneutral->push_back( (unsigned) jet->neutralMultiplicity() );
212     Nmuon->push_back( (unsigned) jet->muonMultiplicity() );
213 bbetchar 1.8 Nelectron->push_back( (unsigned) jet->electronMultiplicity() );
214 bbetchar 1.4 Ndaughters->push_back( (unsigned) jet->numberOfDaughters() );
215 bbetchar 1.2
216     pfjetidloose->push_back(pfLooseJetID( *jet, passLooseCuts ));
217     pfjetidtight->push_back(pfTightJetID( *jet, passTightCuts ));
218 bbetchar 1.8
219     muonsChargeSum->push_back( 0 ) ;
220     electronsChargeSum->push_back( 0 ) ;
221     for(unsigned nMoreMu(Nmuon->back()),nMoreEl(Nelectron->back()), i(0); nMoreMu || nMoreEl; i++) {
222     int pdgId = jet->getPFConstituent(i)->pdgId();
223     if(abs(pdgId)==13) {
224     nMoreMu--;
225     muonsChargeSum->back() += (pdgId>0 ? 1 : -1) ;
226     }
227     if(abs(pdgId)==11) {
228     nMoreEl--;
229     electronsChargeSum->back() += (pdgId>0 ? 1 : -1) ;
230     }
231     }
232 bbetchar 1.2 }
233 bbetchar 1.1 }
234 bbetchar 1.5 if(all.isValid()) {
235     *failedPt=0;
236     for(typename edm::View<T>::const_iterator jet = all->begin(); jet!=all->end(); jet++) {
237     pat::strbitset passLooseCuts( pfLooseJetID .getBitTemplate() );
238     if( !pfLooseJetID( *jet, passLooseCuts ) && jet->pt() > *failedPt) *failedPt = jet->pt();
239     }
240     }
241 bbetchar 1.1
242 bbetchar 1.2 evt.put(FchargedHad, prefix + "FchargedHad" );
243     evt.put(FneutralHad, prefix + "FneutralHad" );
244     evt.put(FchargedEm, prefix + "FchargedEm" );
245     evt.put(FneutralEm, prefix + "FneutralEm" );
246     evt.put(FchargedMu, prefix + "FchargedMu" );
247     evt.put(Ncharged, prefix + "Ncharged" );
248     evt.put(Nneutral, prefix + "Nneutral" );
249     evt.put(Nmuon, prefix + "Nmuon" );
250 bbetchar 1.8 evt.put(Nelectron, prefix + "Nelectron" );
251 bbetchar 1.4 evt.put(Ndaughters, prefix + "Ndaughters" );
252 bbetchar 1.1
253 bbetchar 1.2 evt.put( pfjetidloose, prefix + "PFJetIDloose" );
254     evt.put( pfjetidtight, prefix + "PFJetIDtight" );
255 bbetchar 1.5 evt.put( failedPt, prefix + "FailedPtMax" );
256 bbetchar 1.1
257 bbetchar 1.8 evt.put( muonsChargeSum, prefix + "MuonsChargeSum" );
258     evt.put( electronsChargeSum, prefix + "ElectronsChargeSum" );
259 bbetchar 1.1 }
260    
261    
262 bbetchar 1.4 template<class T>
263     void Tuple_Jet<T>::
264     initGen() {
265     produces <LorentzV> (prefix + "DeltaMETSmear");
266     produces <LorentzV> (prefix + "DeltaMETSmearUp");
267     produces <LorentzV> (prefix + "DeltaMETSmearDown");
268 bbetchar 1.1
269 bbetchar 1.4 produces <std::vector<float> > (prefix + "Smear");
270     produces <std::vector<float> > (prefix + "SmearUp");
271     produces <std::vector<float> > (prefix + "SmearDown");
272     produces <std::vector<LorentzV> > (prefix + "GenP4" );
273     produces <std::vector<int> > (prefix + "GenFlavor" );
274 bbetchar 1.7 produces <std::vector<int> > (prefix + "GenPdgId" );
275     produces <std::vector<int> > (prefix + "GenMotherPdgId" );
276 bbetchar 1.1 }
277    
278    
279 bbetchar 1.4 template<class T>
280     void Tuple_Jet<T>::
281 bbetchar 1.5 produceGen(edm::Event& evt, const Handle_t& jets, const Handle_t& all){
282 bbetchar 1.4
283     std::auto_ptr<LorentzV> deltaMET( new LorentzV() );
284     std::auto_ptr<LorentzV> deltaMETu( new LorentzV() );
285     std::auto_ptr<LorentzV> deltaMETd( new LorentzV() );
286    
287     if(all.isValid()) {
288     for (typename edm::View<T>::const_iterator jet = all->begin(); jet!=all->end(); ++jet ) {
289     const reco::GenJet * gen = jet->genJet();
290 bbetchar 1.7 if( !gen ) continue;
291 bbetchar 1.5 const int bin( dataMcResRatio.FindFixBin(fabs(gen->eta())) );
292 bbetchar 1.4 const double
293 bbetchar 1.5 c( dataMcResRatio.GetBinContent(bin) ),
294     cu( c + dataMcResRatio.GetBinError(bin) ),
295 bbetchar 1.4 cd( 2*c - cu ),
296     pTgOverPt( gen->pt() / jet->pt() );
297     *deltaMET += (1 - std::max(0., c + (1-c ) * pTgOverPt ) ) * jet->p4() ;
298     *deltaMETu += (1 - std::max(0., cu + (1-cu) * pTgOverPt ) ) * jet->p4() ;
299     *deltaMETd += (1 - std::max(0., cd + (1-cd) * pTgOverPt ) ) * jet->p4() ;
300     }
301     }
302     evt.put( deltaMET, prefix + "DeltaMETSmear");
303     evt.put( deltaMETu, prefix + "DeltaMETSmearUp");
304     evt.put( deltaMETd, prefix + "DeltaMETSmearDown");
305    
306     std::auto_ptr<std::vector<LorentzV> > p4(new std::vector<LorentzV>() );
307     std::auto_ptr<std::vector<int> > flavor(new std::vector<int>() );
308     std::auto_ptr<std::vector<float> > smear(new std::vector<float>() );
309     std::auto_ptr<std::vector<float> > smearu(new std::vector<float>() );
310     std::auto_ptr<std::vector<float> > smeard(new std::vector<float>() );
311 bbetchar 1.7 std::auto_ptr<std::vector<int> > pdgId(new std::vector<int>() );
312     std::auto_ptr<std::vector<int> > motherPdgId(new std::vector<int>() );
313 bbetchar 1.1
314 bbetchar 1.4 if(jets.isValid()) {
315     for (typename edm::View<T>::const_iterator jet = jets->begin(); jet!=jets->end(); ++jet ) {
316     const reco::GenJet * gen = jet->genJet();
317 bbetchar 1.7 const reco::GenParticle * genP = jet->genParton();
318 bbetchar 1.5 const int bin( gen ? dataMcResRatio.FindFixBin(fabs(gen->eta())) : -1 );
319 bbetchar 1.4 const double
320 bbetchar 1.5 c( dataMcResRatio.GetBinContent(bin) ),
321     cu( c + dataMcResRatio.GetBinError(bin) ),
322 bbetchar 1.4 cd( 2*c - cu ),
323     pTgOverPt( gen ? gen->pt() / jet->pt() : 0 );
324    
325     p4->push_back( gen ? LorentzV(gen->pt(), gen->eta(), gen->phi(), gen->mass()) : LorentzV());
326     flavor->push_back(jet->partonFlavour());
327     smear->push_back( gen ? std::max(0., c + (1-c ) * pTgOverPt ) : 1.0);
328     smearu->push_back( gen ? std::max(0., cu + (1-cu) * pTgOverPt ) : 1.0);
329     smeard->push_back( gen ? std::max(0., cd + (1-cd) * pTgOverPt ) : 1.0);
330 bbetchar 1.7 pdgId->push_back( genP ? genP->pdgId() : 0.0 );
331     motherPdgId->push_back( genP ? genP->mother()->pdgId() : 0.0 );
332 bbetchar 1.4 }
333     }
334     evt.put( p4, prefix + "GenP4" );
335     evt.put( flavor, prefix + "GenFlavor" );
336     evt.put( smear , prefix + "Smear" );
337     evt.put( smearu, prefix + "SmearUp" );
338     evt.put( smeard, prefix + "SmearDown" );
339 bbetchar 1.7 evt.put( pdgId, prefix + "GenPdgId" );
340     evt.put( motherPdgId, prefix + "GenMotherPdgId" );
341 bbetchar 1.1 }
342    
343     #endif
344    
345