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