JetFitAnalyzer/test/PartonResAnalyzer.cc

/* A JetFitAnalyzer that makes plots of parton resolution */
#include "UserCode/JetFitAnalyzer/interface/JetFitAnalyzer.h"

#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
#include "DataFormats/JetReco/interface/Jet.h"

#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"

#include <iostream>
#include <sstream>

#include "TF2.h"
#include "TNtuple.h"
#include "TH1.h"

#define PI 3.141593

using namespace std;

class PartonResAnalyzer : public JetFitAnalyzer {
public:
  explicit PartonResAnalyzer( const edm::ParameterSet&);
  ~PartonResAnalyzer() {}

private:
  virtual void beginJob(const edm::EventSetup &es);
  virtual void analyze_results(HistoFitter::FitResults, std::vector<HistoFitter::Trouble>,
                               TH2 *, const edm::Event &evt);

  // parton resolution
  TH1D *etaDev;
  TH1D *phiDev;
  TH1D *energyDev;

  // W Mmass resolution
  TH1D * WMassRes;
  TH1D *matchWMassRes;

  // Top mass resolution
  TH1D *topMassRes;
}; 
PartonResAnalyzer::PartonResAnalyzer(const edm::ParameterSet &pSet) 
  : JetFitAnalyzer(pSet) // this is important!
{
}

void PartonResAnalyzer::beginJob(const edm::EventSetup &es) {
  edm::Service<TFileService> fs;
  etaDev = fs->make<TH1D>("etaDev", "Deviation from True Parton Eta",
                          1000, -1.0, 1.0);
  phiDev = fs->make<TH1D>("phiDev", "Deviation from True Parton Phi",
                          1000, -1.0, 1.0);
  energyDev = fs->make<TH1D>("energyDev", "Deviation from True Parton Energy",
                             1000, -100.0, 100.0);

  WMassRes = fs->make<TH1D>("WMassRes", "Mass width of jet pairs",
                                1000, 0.0, 200.0);
  matchWMassRes = fs->make<TH1D>("matchWMassRes", "Mass of jet pairs matched to true W partons", 1000, 0.0, 200.0);


  topMassRes = fs->make<TH1D>("topMassRes", "Top Mass", 1000, 0.0, 200.0);
}

                                 
void PartonResAnalyzer::analyze_results(HistoFitter::FitResults r,
                                     std::vector<HistoFitter::Trouble> t,
                                     TH2 *hist_orig, const edm::Event &evt) {
  // perform analysis of fit results
  edm::Service<TFileService> fs;
  unsigned i = r.pars.size()-1;
  // save jets found
  vector<reco::Jet::LorentzVector> jets;
  unsigned nGauss = r.pval[i].size() / 4;
  for (unsigned j = 0; j < nGauss; j++) {
    double energy = r.pval[i][4*j];
    double eta = r.pval[i][4*j + 1];
    double phi = r.pval[i][4*j + 2];
    double width = r.pval[i][4*j + 3];
    //true jet width
    width = sqrt(width*width - 0.0025);
    
    reco::LeafCandidate::LorentzVector P4;
    reco::LeafCandidate::Point vertex(0.0, 0.0, 0.0);
    P4.SetE(energy);
    // SetEta and SetPhi don't seem to be working
    //P4.SetEta(eta);
    //P4.SetPhi(phi);
    double theta = 2.0*atan(exp(-eta));
    P4.SetPx(energy*sin(theta)*cos(phi));
    P4.SetPy(energy*sin(theta)*sin(phi));
    P4.SetPz(energy*cos(theta));
    jets.push_back(P4);
  }

  // get the generator-level info
  edm::Handle<edm::HepMCProduct> genProduct;
  evt.getByLabel("generator", genProduct);
  const HepMC::GenEvent *genEvt = genProduct->GetEvent();
  
  // obtain true W decay partons
  std::vector<reco::Jet::LorentzVector> trueJets;
  for (HepMC::GenEvent::particle_const_iterator ppit 
         = genEvt->particles_begin();
       ppit != genEvt->particles_end(); ppit++) {
    if ((*ppit)->pdg_id() == 24 || (*ppit)->pdg_id() == -24) {
      // found W --iteratr over daugheters
      HepMC::GenVertex *vtx = (*ppit)->end_vertex();
      if (vtx != 0) {
        for (HepMC::GenVertex::particles_out_const_iterator pit
               = vtx->particles_out_const_begin();
             pit != vtx->particles_out_const_end(); pit++) {
          reco::Jet::LorentzVector parton;
          parton.SetE((*pit)->momentum().e());
          parton.SetPx((*pit)->momentum().px());
          parton.SetPy((*pit)->momentum().py());
          parton.SetPz((*pit)->momentum().pz());
          // if parton falls within histo
          if (parton.Eta() > hist_orig->GetXaxis()->GetXmin()
              && parton.Eta() < hist_orig->GetXaxis()->GetXmax()
              && parton.Phi() > hist_orig->GetYaxis()->GetXmin()
              && parton.Phi() < hist_orig->GetYaxis()->GetXmax())
            trueJets.push_back(parton);
        }
      }
    }
  }
  cout << trueJets.size() << endl;

  // find closest fit jets to true jets and compute deviations
  std::vector<unsigned> closestJetIndex(trueJets.size());
  for (unsigned i = 0; i < trueJets.size(); i++) {
    // find the closest jet
    double smallestDR = 1.0e10;
    for (unsigned i2 = 0; i2 < jets.size(); i2++) {
      double DR = sqrt(pow(trueJets[i].Eta() - jets[i2].Eta(), 2.0)
                       + pow(trueJets[i].Phi() - jets[i2].Phi(), 2.0));
      if (DR < smallestDR) {
        smallestDR = DR;
        closestJetIndex[i] = i2;
      }
    }

    // compute the deviations in eta, phi, and E
    reco::Jet::LorentzVector dev = jets[closestJetIndex[i]] - trueJets[i];
    energyDev->Fill(static_cast<double>(dev.E()));
    etaDev->Fill(static_cast<double>(jets[closestJetIndex[i]].Eta()
                                     - trueJets[i].Eta()));
    phiDev->Fill(static_cast<double>(jets[closestJetIndex[i]].Phi()
                                     - trueJets[i].Phi()));


  }

  // compute the mass of each jet pair
  for (unsigned s1 = 0; s1 < jets.size() - 1; s1++) {
    for (unsigned s2 = s1+1; s2 < jets.size(); s2++) {
      reco::Jet::LorentzVector pair = jets[s1] + jets[s2];
      WMassRes->Fill(pair.M());
    }
  }

  // compute the mass of W from pairs matched to true W partons
  for (unsigned s1 = 0; s1 < trueJets.size() - 1; s1++) {
    for (unsigned s2 = s1+1; s2 < trueJets.size(); s2++) {
      reco::Jet::LorentzVector pair = jets[closestJetIndex[s1]] 
        + jets[closestJetIndex[s2]];
      matchWMassRes->Fill(pair.M());
    }
  }

  // compute the top mass
  reco::Jet::LorentzVector topVec;
  topVec.SetE(0.0);
  topVec.SetPx(0.0);
  topVec.SetPy(0.0);
  topVec.SetPz(0.0);
  for (unsigned s  = 0; s < jets.size(); s++) {
    topVec = topVec + jets[s];
  }
  if (jets.size() > 3)
    topMassRes->Fill(topVec.M());
}

DEFINE_FWK_MODULE(PartonResAnalyzer);
Revision:	1.7
Committed:	Tue Apr 20 15:40:19 2010 UTC (15 years ago) by dnisson
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V02-04-04, HEAD
Changes since 1.6:	+2 -2 lines
Log Message:	use correct true jets for matching
#	User	Rev	Content
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	dnisson	1.3	// parton resolution
33	dnisson	1.1	TH1D *etaDev;
34			TH1D *phiDev;
35			TH1D *energyDev;
36
37	dnisson	1.3	// W Mmass resolution
38			TH1D * WMassRes;
39			TH1D *matchWMassRes;
40	dnisson	1.4
41			// Top mass resolution
42			TH1D *topMassRes;
43	dnisson	1.3	};
44	dnisson	1.1	PartonResAnalyzer::PartonResAnalyzer(const edm::ParameterSet &pSet)
45			: JetFitAnalyzer(pSet) // this is important!
46			{
47			}
48
49			void PartonResAnalyzer::beginJob(const edm::EventSetup &es) {
50			edm::Service<TFileService> fs;
51			etaDev = fs->make<TH1D>("etaDev", "Deviation from True Parton Eta",
52			1000, -1.0, 1.0);
53			phiDev = fs->make<TH1D>("phiDev", "Deviation from True Parton Phi",
54			1000, -1.0, 1.0);
55			energyDev = fs->make<TH1D>("energyDev", "Deviation from True Parton Energy",
56			1000, -100.0, 100.0);
57	dnisson	1.3
58			WMassRes = fs->make<TH1D>("WMassRes", "Mass width of jet pairs",
59			1000, 0.0, 200.0);
60			matchWMassRes = fs->make<TH1D>("matchWMassRes", "Mass of jet pairs matched to true W partons", 1000, 0.0, 200.0);
61	dnisson	1.4
62
63			topMassRes = fs->make<TH1D>("topMassRes", "Top Mass", 1000, 0.0, 200.0);
64	dnisson	1.1	}
65
66	dnisson	1.3
67	dnisson	1.1	void PartonResAnalyzer::analyze_results(HistoFitter::FitResults r,
68			std::vector<HistoFitter::Trouble> t,
69			TH2 *hist_orig, const edm::Event &evt) {
70			// perform analysis of fit results
71			edm::Service<TFileService> fs;
72			unsigned i = r.pars.size()-1;
73			// save jets found
74			vector<reco::Jet::LorentzVector> jets;
75			unsigned nGauss = r.pval[i].size() / 4;
76			for (unsigned j = 0; j < nGauss; j++) {
77			double energy = r.pval[i][4*j];
78			double eta = r.pval[i][4*j + 1];
79			double phi = r.pval[i][4*j + 2];
80			double width = r.pval[i][4*j + 3];
81			//true jet width
82			width = sqrt(width*width - 0.0025);
83
84			reco::LeafCandidate::LorentzVector P4;
85			reco::LeafCandidate::Point vertex(0.0, 0.0, 0.0);
86			P4.SetE(energy);
87			// SetEta and SetPhi don't seem to be working
88			//P4.SetEta(eta);
89			//P4.SetPhi(phi);
90			double theta = 2.0*atan(exp(-eta));
91			P4.SetPx(energysin(theta)cos(phi));
92			P4.SetPy(energysin(theta)sin(phi));
93			P4.SetPz(energy*cos(theta));
94			jets.push_back(P4);
95			}
96
97			// get the generator-level info
98			edm::Handle<edm::HepMCProduct> genProduct;
99			evt.getByLabel("generator", genProduct);
100			const HepMC::GenEvent *genEvt = genProduct->GetEvent();
101
102			// obtain true W decay partons
103			std::vector<reco::Jet::LorentzVector> trueJets;
104			for (HepMC::GenEvent::particle_const_iterator ppit
105			= genEvt->particles_begin();
106			ppit != genEvt->particles_end(); ppit++) {
107			if ((ppit)->pdg_id() == 24 \|\| (ppit)->pdg_id() == -24) {
108			// found W --iteratr over daugheters
109			HepMC::GenVertex vtx = (ppit)->end_vertex();
110			if (vtx != 0) {
111			for (HepMC::GenVertex::particles_out_const_iterator pit
112			= vtx->particles_out_const_begin();
113			pit != vtx->particles_out_const_end(); pit++) {
114			reco::Jet::LorentzVector parton;
115			parton.SetE((*pit)->momentum().e());
116			parton.SetPx((*pit)->momentum().px());
117			parton.SetPy((*pit)->momentum().py());
118			parton.SetPz((*pit)->momentum().pz());
119			// if parton falls within histo
120			if (parton.Eta() > hist_orig->GetXaxis()->GetXmin()
121			&& parton.Eta() < hist_orig->GetXaxis()->GetXmax()
122			&& parton.Phi() > hist_orig->GetYaxis()->GetXmin()
123			&& parton.Phi() < hist_orig->GetYaxis()->GetXmax())
124			trueJets.push_back(parton);
125			}
126			}
127			}
128			}
129			cout << trueJets.size() << endl;
130
131	dnisson	1.3	// find closest fit jets to true jets and compute deviations
132			std::vector<unsigned> closestJetIndex(trueJets.size());
133	dnisson	1.1	for (unsigned i = 0; i < trueJets.size(); i++) {
134			// find the closest jet
135	dnisson	1.6	double smallestDR = 1.0e10;
136	dnisson	1.1	for (unsigned i2 = 0; i2 < jets.size(); i2++) {
137	dnisson	1.7	double DR = sqrt(pow(trueJets[i].Eta() - jets[i2].Eta(), 2.0)
138			+ pow(trueJets[i].Phi() - jets[i2].Phi(), 2.0));
139	dnisson	1.6	if (DR < smallestDR) {
140			smallestDR = DR;
141	dnisson	1.3	closestJetIndex[i] = i2;
142	dnisson	1.1	}
143			}
144
145			// compute the deviations in eta, phi, and E
146	dnisson	1.3	reco::Jet::LorentzVector dev = jets[closestJetIndex[i]] - trueJets[i];
147	dnisson	1.1	energyDev->Fill(static_cast<double>(dev.E()));
148	dnisson	1.3	etaDev->Fill(static_cast<double>(jets[closestJetIndex[i]].Eta()
149	dnisson	1.2	- trueJets[i].Eta()));
150	dnisson	1.3	phiDev->Fill(static_cast<double>(jets[closestJetIndex[i]].Phi()
151	dnisson	1.2	- trueJets[i].Phi()));
152	dnisson	1.3
153
154			}
155
156			// compute the mass of each jet pair
157			for (unsigned s1 = 0; s1 < jets.size() - 1; s1++) {
158			for (unsigned s2 = s1+1; s2 < jets.size(); s2++) {
159			reco::Jet::LorentzVector pair = jets[s1] + jets[s2];
160			WMassRes->Fill(pair.M());
161			}
162			}
163
164			// compute the mass of W from pairs matched to true W partons
165			for (unsigned s1 = 0; s1 < trueJets.size() - 1; s1++) {
166			for (unsigned s2 = s1+1; s2 < trueJets.size(); s2++) {
167			reco::Jet::LorentzVector pair = jets[closestJetIndex[s1]]
168			+ jets[closestJetIndex[s2]];
169			matchWMassRes->Fill(pair.M());
170			}
171	dnisson	1.1	}
172	dnisson	1.4
173			// compute the top mass
174			reco::Jet::LorentzVector topVec;
175			topVec.SetE(0.0);
176			topVec.SetPx(0.0);
177			topVec.SetPy(0.0);
178			topVec.SetPz(0.0);
179			for (unsigned s = 0; s < jets.size(); s++) {
180			topVec = topVec + jets[s];
181			}
182	dnisson	1.5	if (jets.size() > 3)
183			topMassRes->Fill(topVec.M());
184	dnisson	1.1	}
185
186			DEFINE_FWK_MODULE(PartonResAnalyzer);