| 1 |
#ifndef TUPLE_JET
|
| 2 |
#define TUPLE_JET
|
| 3 |
|
| 4 |
#include "boost/foreach.hpp"
|
| 5 |
#include "TopQuarkAnalysis/TopRefTuple/interface/fTypes.h"
|
| 6 |
#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 |
#include "CondFormats/JetMETObjects/interface/JetResolution.h"
|
| 15 |
#include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
|
| 16 |
#include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
|
| 17 |
#include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
|
| 18 |
#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 |
jCU->setJetPt(jet.pt());
|
| 30 |
try {return jCU->getUncertainty(true);}
|
| 31 |
catch (...) {
|
| 32 |
std::cout << "Tuple_Jet.h: JetCorrectionUncertainty::getUncertainty"
|
| 33 |
<< "threw exception on jet with pt( " << jet.pt() << " ) "
|
| 34 |
<< "and eta ( " << jet.eta() << " )." << std::endl;
|
| 35 |
return 0.0;
|
| 36 |
}
|
| 37 |
}
|
| 38 |
|
| 39 |
|
| 40 |
|
| 41 |
template< typename T >
|
| 42 |
class Tuple_Jet : public edm::EDProducer {
|
| 43 |
public:
|
| 44 |
explicit Tuple_Jet(const edm::ParameterSet&);
|
| 45 |
typedef fTypes::dPolarLorentzV LorentzV;
|
| 46 |
typedef edm::Handle<edm::View<T> > Handle_t;
|
| 47 |
private:
|
| 48 |
void produce(edm::Event&, const edm::EventSetup& );
|
| 49 |
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 |
|
| 52 |
const edm::InputTag jetsTag, allJetsTag;
|
| 53 |
const std::string prefix,jecRecord;
|
| 54 |
const std::vector<std::string> btagNames;
|
| 55 |
const bool pfInfo, genInfo;
|
| 56 |
JetResolution jer;
|
| 57 |
TH1D dataMcResRatio;
|
| 58 |
};
|
| 59 |
|
| 60 |
template<class T> Tuple_Jet<T>::
|
| 61 |
Tuple_Jet(const edm::ParameterSet& cfg) :
|
| 62 |
jetsTag(cfg.getParameter<edm::InputTag>("jetsTag")),
|
| 63 |
allJetsTag(cfg.getParameter<edm::InputTag>("allJetsTag")),
|
| 64 |
prefix(cfg.getParameter<std::string>("prefix")),
|
| 65 |
jecRecord(cfg.getParameter<std::string>("jecRecord")),
|
| 66 |
btagNames(cfg.getParameter<std::vector<std::string> >("bTags")),
|
| 67 |
pfInfo(cfg.getParameter<bool>("pfInfo")),
|
| 68 |
genInfo(cfg.getParameter<bool>("genInfo")),
|
| 69 |
jer(cfg.getParameter<std::string>("jetResolutionFile")),
|
| 70 |
dataMcResRatio("dataMcResRatio",
|
| 71 |
"",
|
| 72 |
cfg.getParameter<std::vector<double> >("resolutionRatioBins").size()-1,
|
| 73 |
&cfg.getParameter<std::vector<double> >("resolutionRatioBins").at(0) )
|
| 74 |
{
|
| 75 |
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 |
BOOST_FOREACH(const std::string& btag, btagNames)
|
| 84 |
produces <std::vector<float> > (prefix + btag);
|
| 85 |
|
| 86 |
if(pfInfo) initPF();
|
| 87 |
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 |
}
|
| 97 |
|
| 98 |
template< typename T >
|
| 99 |
void Tuple_Jet<T>::
|
| 100 |
produce(edm::Event& evt, const edm::EventSetup& setup) {
|
| 101 |
Handle_t jets; evt.getByLabel(jetsTag, jets);
|
| 102 |
Handle_t all; evt.getByLabel(allJetsTag, all);
|
| 103 |
|
| 104 |
std::auto_ptr<std::vector<LorentzV> > p4 ( new std::vector<LorentzV>() ) ;
|
| 105 |
std::auto_ptr<std::vector<float> > jecFactor( new std::vector<float>() ) ;
|
| 106 |
std::auto_ptr<std::vector<float> > jecUnc ( new std::vector<float>() ) ;
|
| 107 |
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 |
std::auto_ptr<std::vector<float> > eta2mom ( new std::vector<float>() ) ;
|
| 111 |
std::auto_ptr<std::vector<float> > phi2mom ( new std::vector<float>() ) ;
|
| 112 |
std::map<std::string, std::vector<float>* > btags;
|
| 113 |
BOOST_FOREACH(const std::string& btag, btagNames)
|
| 114 |
btags[btag] = new std::vector<float>();
|
| 115 |
|
| 116 |
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 |
p4->push_back(LorentzV(jet->pt(),jet->eta(),jet->phi(),jet->mass()));
|
| 120 |
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 |
eta2mom->push_back(jet->etaetaMoment());
|
| 125 |
phi2mom->push_back(jet->phiphiMoment());
|
| 126 |
|
| 127 |
BOOST_FOREACH(const std::string& btag, btagNames)
|
| 128 |
btags[btag]->push_back(jet->bDiscriminator(btag+"BJetTags"));
|
| 129 |
|
| 130 |
TF1 * jerEta = jer.parameterEta("sigma",jet->eta());
|
| 131 |
reso->push_back( jerEta->Eval( jet->pt() ) );
|
| 132 |
delete jerEta;
|
| 133 |
}
|
| 134 |
delete jCU;
|
| 135 |
}
|
| 136 |
|
| 137 |
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 |
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 |
}
|
| 151 |
|
| 152 |
template<class T> void Tuple_Jet<T>::
|
| 153 |
initPF() {
|
| 154 |
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 |
|
| 160 |
produces <std::vector<unsigned> > (prefix + "Ncharged" );
|
| 161 |
produces <std::vector<unsigned> > (prefix + "Nneutral" );
|
| 162 |
produces <std::vector<unsigned> > (prefix + "Nmuon" );
|
| 163 |
produces <std::vector<unsigned> > (prefix + "Nelectron" );
|
| 164 |
produces <std::vector<unsigned> > (prefix + "Ndaughters" );
|
| 165 |
|
| 166 |
produces <std::vector<bool> > ( prefix + "PFJetIDloose" );
|
| 167 |
produces <std::vector<bool> > ( prefix + "PFJetIDtight" );
|
| 168 |
produces <float> ( prefix + "FailedPtMax");
|
| 169 |
|
| 170 |
produces <std::vector<int> > (prefix + "MuonsChargeSum" );
|
| 171 |
produces <std::vector<int> > (prefix + "ElectronsChargeSum" );
|
| 172 |
}
|
| 173 |
|
| 174 |
template<class T> void Tuple_Jet<T>::
|
| 175 |
producePF(edm::Event& evt, const Handle_t& jets, const Handle_t& all) {
|
| 176 |
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 |
std::auto_ptr<std::vector<unsigned> > Nelectron( new std::vector<unsigned>() );
|
| 186 |
std::auto_ptr<std::vector<unsigned> > Ndaughters( new std::vector<unsigned>() );
|
| 187 |
|
| 188 |
std::auto_ptr<std::vector<bool> > pfjetidloose ( new std::vector<bool>() ) ;
|
| 189 |
std::auto_ptr<std::vector<bool> > pfjetidtight ( new std::vector<bool>() ) ;
|
| 190 |
std::auto_ptr<float> failedPt ( new float(-1) ) ;
|
| 191 |
|
| 192 |
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 |
PFJetIDSelectionFunctor
|
| 196 |
pfLooseJetID(PFJetIDSelectionFunctor::FIRSTDATA, PFJetIDSelectionFunctor::LOOSE),
|
| 197 |
pfTightJetID(PFJetIDSelectionFunctor::FIRSTDATA, PFJetIDSelectionFunctor::TIGHT);
|
| 198 |
|
| 199 |
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 |
Nelectron->push_back( (unsigned) jet->electronMultiplicity() );
|
| 214 |
Ndaughters->push_back( (unsigned) jet->numberOfDaughters() );
|
| 215 |
|
| 216 |
pfjetidloose->push_back(pfLooseJetID( *jet, passLooseCuts ));
|
| 217 |
pfjetidtight->push_back(pfTightJetID( *jet, passTightCuts ));
|
| 218 |
|
| 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 |
}
|
| 233 |
}
|
| 234 |
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 |
|
| 242 |
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 |
evt.put(Nelectron, prefix + "Nelectron" );
|
| 251 |
evt.put(Ndaughters, prefix + "Ndaughters" );
|
| 252 |
|
| 253 |
evt.put( pfjetidloose, prefix + "PFJetIDloose" );
|
| 254 |
evt.put( pfjetidtight, prefix + "PFJetIDtight" );
|
| 255 |
evt.put( failedPt, prefix + "FailedPtMax" );
|
| 256 |
|
| 257 |
evt.put( muonsChargeSum, prefix + "MuonsChargeSum" );
|
| 258 |
evt.put( electronsChargeSum, prefix + "ElectronsChargeSum" );
|
| 259 |
}
|
| 260 |
|
| 261 |
|
| 262 |
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 |
|
| 269 |
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 |
produces <std::vector<int> > (prefix + "GenPdgId" );
|
| 275 |
produces <std::vector<int> > (prefix + "GenMotherPdgId" );
|
| 276 |
}
|
| 277 |
|
| 278 |
|
| 279 |
template<class T>
|
| 280 |
void Tuple_Jet<T>::
|
| 281 |
produceGen(edm::Event& evt, const Handle_t& jets, const Handle_t& all){
|
| 282 |
|
| 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 |
if( !gen ) continue;
|
| 291 |
const int bin( dataMcResRatio.FindFixBin(fabs(gen->eta())) );
|
| 292 |
const double
|
| 293 |
c( dataMcResRatio.GetBinContent(bin) ),
|
| 294 |
cu( c + dataMcResRatio.GetBinError(bin) ),
|
| 295 |
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 |
std::auto_ptr<std::vector<int> > pdgId(new std::vector<int>() );
|
| 312 |
std::auto_ptr<std::vector<int> > motherPdgId(new std::vector<int>() );
|
| 313 |
|
| 314 |
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 |
const reco::GenParticle * genP = jet->genParton();
|
| 318 |
const int bin( gen ? dataMcResRatio.FindFixBin(fabs(gen->eta())) : -1 );
|
| 319 |
const double
|
| 320 |
c( dataMcResRatio.GetBinContent(bin) ),
|
| 321 |
cu( c + dataMcResRatio.GetBinError(bin) ),
|
| 322 |
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 |
pdgId->push_back( genP ? genP->pdgId() : 0.0 );
|
| 331 |
motherPdgId->push_back( genP ? genP->mother()->pdgId() : 0.0 );
|
| 332 |
}
|
| 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 |
evt.put( pdgId, prefix + "GenPdgId" );
|
| 340 |
evt.put( motherPdgId, prefix + "GenMotherPdgId" );
|
| 341 |
}
|
| 342 |
|
| 343 |
#endif
|
| 344 |
|
| 345 |
|
