| 1 |
/* A JetFitAnalyzer that makes plots of parton resolution */
|
| 2 |
#include "UserCode/JetFitAnalyzer/interface/JetFitAnalyzer.h"
|
| 3 |
|
| 4 |
#include "FWCore/ServiceRegistry/interface/Service.h"
|
| 5 |
#include "FWCore/MessageLogger/interface/MessageLogger.h"
|
| 6 |
#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
|
| 7 |
#include "DataFormats/JetReco/interface/Jet.h"
|
| 8 |
|
| 9 |
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
|
| 10 |
|
| 11 |
#include <iostream>
|
| 12 |
#include <sstream>
|
| 13 |
|
| 14 |
#include "TF2.h"
|
| 15 |
#include "TNtuple.h"
|
| 16 |
#include "TH1.h"
|
| 17 |
|
| 18 |
#define PI 3.141593
|
| 19 |
|
| 20 |
using namespace std;
|
| 21 |
|
| 22 |
class QCDResAnalyzer : public JetFitAnalyzer {
|
| 23 |
public:
|
| 24 |
explicit QCDResAnalyzer( const edm::ParameterSet&);
|
| 25 |
~QCDResAnalyzer() {}
|
| 26 |
|
| 27 |
private:
|
| 28 |
virtual void beginJob(const edm::EventSetup &es);
|
| 29 |
virtual void analyze_results(HistoFitter::FitResults, std::vector<HistoFitter::Trouble>,
|
| 30 |
TH2 *, const edm::Event &evt);
|
| 31 |
|
| 32 |
// parton resolution
|
| 33 |
TH1D *etaDev;
|
| 34 |
TH1D *phiDev;
|
| 35 |
TH1D *energyDev;
|
| 36 |
|
| 37 |
};
|
| 38 |
QCDResAnalyzer::QCDResAnalyzer(const edm::ParameterSet &pSet)
|
| 39 |
: JetFitAnalyzer(pSet) // this is important!
|
| 40 |
{
|
| 41 |
}
|
| 42 |
|
| 43 |
void QCDResAnalyzer::beginJob(const edm::EventSetup &es) {
|
| 44 |
edm::Service<TFileService> fs;
|
| 45 |
etaDev = fs->make<TH1D>("etaDev", "Deviation from True Parton Eta",
|
| 46 |
1000, -1.0, 1.0);
|
| 47 |
phiDev = fs->make<TH1D>("phiDev", "Deviation from True Parton Phi",
|
| 48 |
1000, -1.0, 1.0);
|
| 49 |
energyDev = fs->make<TH1D>("energyDev", "Deviation from True Parton Energy",
|
| 50 |
1000, -100.0, 100.0);
|
| 51 |
|
| 52 |
}
|
| 53 |
|
| 54 |
|
| 55 |
void QCDResAnalyzer::analyze_results(HistoFitter::FitResults r,
|
| 56 |
std::vector<HistoFitter::Trouble> t,
|
| 57 |
TH2 *hist_orig, const edm::Event &evt) {
|
| 58 |
// perform analysis of fit results
|
| 59 |
edm::Service<TFileService> fs;
|
| 60 |
unsigned i = r.pars.size()-1;
|
| 61 |
// save jets found
|
| 62 |
vector<reco::Jet::LorentzVector> jets;
|
| 63 |
unsigned nGauss = r.pval[i].size() / 4;
|
| 64 |
for (unsigned j = 0; j < nGauss; j++) {
|
| 65 |
double energy = r.pval[i][4*j];
|
| 66 |
double eta = r.pval[i][4*j + 1];
|
| 67 |
double phi = r.pval[i][4*j + 2];
|
| 68 |
double width = r.pval[i][4*j + 3];
|
| 69 |
//true jet width
|
| 70 |
width = sqrt(width*width - 0.0025);
|
| 71 |
|
| 72 |
reco::LeafCandidate::LorentzVector P4;
|
| 73 |
reco::LeafCandidate::Point vertex(0.0, 0.0, 0.0);
|
| 74 |
P4.SetE(energy);
|
| 75 |
// SetEta and SetPhi don't seem to be working
|
| 76 |
//P4.SetEta(eta);
|
| 77 |
//P4.SetPhi(phi);
|
| 78 |
double theta = 2.0*atan(exp(-eta));
|
| 79 |
P4.SetPx(energy*sin(theta)*cos(phi));
|
| 80 |
P4.SetPy(energy*sin(theta)*sin(phi));
|
| 81 |
P4.SetPz(energy*cos(theta));
|
| 82 |
jets.push_back(P4);
|
| 83 |
}
|
| 84 |
|
| 85 |
// get the generator-level info
|
| 86 |
edm::Handle<edm::HepMCProduct> genProduct;
|
| 87 |
evt.getByLabel("generator", genProduct);
|
| 88 |
const HepMC::GenEvent *genEvt = genProduct->GetEvent();
|
| 89 |
|
| 90 |
// obtain true dijet partons
|
| 91 |
std::vector<reco::Jet::LorentzVector> trueJets;
|
| 92 |
for (HepMC::GenEvent::particle_const_iterator ppit
|
| 93 |
= genEvt->particles_begin();
|
| 94 |
ppit != genEvt->particles_end(); ppit++) {
|
| 95 |
if ((*ppit)->pdg_id() == 21 || (*ppit)->pdg_id() == -21) {
|
| 96 |
// found gluon --iterate over daugheters
|
| 97 |
HepMC::GenVertex *vtx = (*ppit)->end_vertex();
|
| 98 |
if (vtx != 0) {
|
| 99 |
for (HepMC::GenVertex::particles_out_const_iterator pit
|
| 100 |
= vtx->particles_out_const_begin();
|
| 101 |
pit != vtx->particles_out_const_end(); pit++) {
|
| 102 |
reco::Jet::LorentzVector parton;
|
| 103 |
parton.SetE((*pit)->momentum().e());
|
| 104 |
parton.SetPx((*pit)->momentum().px());
|
| 105 |
parton.SetPy((*pit)->momentum().py());
|
| 106 |
parton.SetPz((*pit)->momentum().pz());
|
| 107 |
// if parton falls within histo
|
| 108 |
if (parton.Eta() > hist_orig->GetXaxis()->GetXmin()
|
| 109 |
&& parton.Eta() < hist_orig->GetXaxis()->GetXmax()
|
| 110 |
&& parton.Phi() > hist_orig->GetYaxis()->GetXmin()
|
| 111 |
&& parton.Phi() < hist_orig->GetYaxis()->GetXmax())
|
| 112 |
trueJets.push_back(parton);
|
| 113 |
}
|
| 114 |
}
|
| 115 |
}
|
| 116 |
}
|
| 117 |
cout << trueJets.size() << endl;
|
| 118 |
|
| 119 |
// find closest fit jets to true jets and compute deviations
|
| 120 |
std::vector<unsigned> closestJetIndex(trueJets.size());
|
| 121 |
for (unsigned i = 0; i < trueJets.size(); i++) {
|
| 122 |
// find the closest jet
|
| 123 |
double smallestDP = 1.0e10;
|
| 124 |
for (unsigned i2 = 0; i2 < jets.size(); i2++) {
|
| 125 |
double DP = (trueJets[i] - jets[i2]).P();
|
| 126 |
if (DP < smallestDP) {
|
| 127 |
smallestDP = DP;
|
| 128 |
closestJetIndex[i] = i2;
|
| 129 |
}
|
| 130 |
}
|
| 131 |
|
| 132 |
// compute the deviations in eta, phi, and E
|
| 133 |
reco::Jet::LorentzVector dev = jets[closestJetIndex[i]] - trueJets[i];
|
| 134 |
energyDev->Fill(static_cast<double>(dev.E()));
|
| 135 |
etaDev->Fill(static_cast<double>(jets[closestJetIndex[i]].Eta()
|
| 136 |
- trueJets[i].Eta()));
|
| 137 |
phiDev->Fill(static_cast<double>(jets[closestJetIndex[i]].Phi()
|
| 138 |
- trueJets[i].Phi()));
|
| 139 |
|
| 140 |
|
| 141 |
}
|
| 142 |
}
|
| 143 |
|
| 144 |
DEFINE_FWK_MODULE(QCDResAnalyzer);
|
