JetFitAnalyzer/test/JetFinderAnalyzer.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 <iostream>
#include <sstream>

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

 #define PI 3.141593

using namespace std;

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

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

  TH1D *pairMassHist;
  TH1D *topMassHist;
};

JetFinderAnalyzer::JetFinderAnalyzer(const edm::ParameterSet &pSet) 
  : JetFitAnalyzer(pSet) // this is important!
{
}

void JetFinderAnalyzer::beginJob(const edm::EventSetup &es) {
  edm::Service<TFileService> fs;
  pairMassHist = fs->make<TH1D>("pairMass", "Mass of Jet Pairs",
                                100, 0.0, 250.0);
  
}

double evalFitFunction(HistoFitter::FitResults r, double x, double y,
                       int i) {
  unsigned nGauss = r.pars[i].size() / 4;
  double fitVal = 0.0;
  for (unsigned j = 0; j < nGauss; j++) {
    double N = r.pval[i][4*j];
    double mu_x = r.pval[i][4*j + 1];
    double mu_y = r.pval[i][4*j + 2];
    double sig = r.pval[i][4*j + 3];

    double rel_x = x - mu_x; double rel_y = y - mu_y;
    fitVal += (N / 2.0 / M_PI / sig / sig) 
      * exp(-(rel_x * rel_x + rel_y * rel_y)/2.0/sig/sig);
  }
  return fitVal;
}

void JetFinderAnalyzer::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;
  for (unsigned i = 0; i < r.pars.size(); i++) {
    ostringstream fitHistoOss;
    fitHistoOss << hist_orig->GetName() << "_fit" << i << flush;
    TH2D *fitHisto = fs->make<TH2D>(fitHistoOss.str().c_str(),
                                  ("Fitted distribution to "
                                     +std::string(hist_orig->GetName())).c_str(),
                                    hist_orig->GetNbinsX(),
                                    hist_orig->GetXaxis()->GetXmin(),
                                    hist_orig->GetXaxis()->GetXmax(),
                                    hist_orig->GetNbinsY(),
                                    hist_orig->GetYaxis()->GetXmin(),
                                    hist_orig->GetYaxis()->GetXmax());
    
    double Xlo = fitHisto->GetXaxis()->GetXmin();
    double Xhi = fitHisto->GetXaxis()->GetXmax();
    double Ylo = fitHisto->GetYaxis()->GetXmin();
    double Yhi = fitHisto->GetYaxis()->GetXmax();
    double XbinSize = (Xhi - Xlo) / static_cast<double>(fitHisto->GetNbinsX());
    double YbinSize = (Yhi - Ylo) / static_cast<double>(fitHisto->GetNbinsY());
    
    for (int ib = 1; ib <= fitHisto->GetNbinsX(); ib++) {
      for (int jb = 1; jb <= fitHisto->GetNbinsY(); jb++) {
        double x = (static_cast<double>(ib) - 0.5) * XbinSize + Xlo;
        double y = (static_cast<double>(jb) - 0.5) * YbinSize + Ylo;
        fitHisto->SetBinContent(ib, jb, evalFitFunction(r, x, y, i) * XbinSize * YbinSize);
      }
    }
    
    // save fit results to an ntuple
    ostringstream fitResultsOss;
    fitResultsOss << hist_orig->GetName() << "_results" << i << flush;
    TNtuple *rNtuple = fs->make<TNtuple>(fitResultsOss.str().c_str(),
                                         ("Fit results for "+std::string(hist_orig->GetName())).c_str(),
                                         "N:mu_x:mu_y:sigma");
    unsigned nGauss = r.pval[i].size() / 4;
    for (unsigned j = 0; j < nGauss; j++) {
      rNtuple->Fill(r.pval[i][4*j], r.pval[i][4*j+1], r.pval[i][4*j+2],
                    r.pval[i][4*j+3]);
    }
    
    // save jets found
    vector<reco::Jet::LorentzVector> jets;
    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);
    }
    // compute pair masses
    for (unsigned s1 = 0; s1 < jets.size() - 1; s1++) {
      for (unsigned s2 = s1 + 1; s2 < jets.size(); s2++) {
        // create TOTAL MOMENTUM VECTOR
        reco::Jet::LorentzVector P4 = jets[s1] + jets[s2];
        pairMassHist->Fill(static_cast<double>(P4.M()));
      }
    }
   
    // save chisquares to ntuple
    for (unsigned j = 0; j < r.chisquare.size(); j++) {
      ostringstream csNtupleName, csNtupleTitle;
      csNtupleName << hist_orig->GetName() << "_chi2_" << j << flush;
      csNtupleTitle << "Chisquare "<<j<<" for histo "<<hist_orig->GetName()
                    << flush;
      TNtuple *csNtuple = fs->make<TNtuple>(csNtupleName.str().c_str(),
                                            csNtupleTitle.str().c_str(),
                                            "chisq");
      csNtuple->Fill(r.chisquare[j]);
    }

    // save jets found
    
  }
}

DEFINE_FWK_MODULE(JetFinderAnalyzer);
Revision:	1.30
Committed:	Thu Apr 8 05:00:01 2010 UTC (15 years ago) by dnisson
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V02-04-04, V02-04-03, V02-04-02, V02-04-01, V02-04-00, HEAD
Changes since 1.29:	+3 -4 lines
Log Message:	Fixed bug in fcn, added parton resolution code
#	User	Rev	Content
1	dnisson	1.30	/* A JetFitAnalyzer that makes plots of parton resolution */
2	dnisson	1.1	#include "UserCode/JetFitAnalyzer/interface/JetFitAnalyzer.h"
3
4			#include "FWCore/ServiceRegistry/interface/Service.h"
5	dnisson	1.2	#include "FWCore/MessageLogger/interface/MessageLogger.h"
6	dnisson	1.1	#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
7	dnisson	1.29	#include "DataFormats/JetReco/interface/Jet.h"
8	dnisson	1.1
9	dnisson	1.18	#include <iostream>
10	dnisson	1.1	#include <sstream>
11
12			#include "TF2.h"
13	dnisson	1.26	#include "TNtuple.h"
14	dnisson	1.29	#include "TH1.h"
15	dnisson	1.1
16	dnisson	1.29	#define PI 3.141593
17	dnisson	1.1
18			using namespace std;
19
20			class JetFinderAnalyzer : public JetFitAnalyzer {
21			public:
22			explicit JetFinderAnalyzer( const edm::ParameterSet&);
23			~JetFinderAnalyzer() {}
24
25			private:
26	dnisson	1.15	virtual void beginJob(const edm::EventSetup &es);
27	dnisson	1.26	virtual void analyze_results(HistoFitter::FitResults, std::vector<HistoFitter::Trouble>,
28	dnisson	1.30	TH2 *,const edm::Event &evt);
29	dnisson	1.29
30			TH1D *pairMassHist;
31			TH1D *topMassHist;
32	dnisson	1.1	};
33
34			JetFinderAnalyzer::JetFinderAnalyzer(const edm::ParameterSet &pSet)
35			: JetFitAnalyzer(pSet) // this is important!
36			{
37	dnisson	1.26	}
38	dnisson	1.2
39	dnisson	1.26	void JetFinderAnalyzer::beginJob(const edm::EventSetup &es) {
40	dnisson	1.29	edm::Service<TFileService> fs;
41			pairMassHist = fs->make<TH1D>("pairMass", "Mass of Jet Pairs",
42			100, 0.0, 250.0);
43
44	dnisson	1.1	}
45
46	dnisson	1.28	double evalFitFunction(HistoFitter::FitResults r, double x, double y,
47			int i) {
48			unsigned nGauss = r.pars[i].size() / 4;
49	dnisson	1.26	double fitVal = 0.0;
50	dnisson	1.28	for (unsigned j = 0; j < nGauss; j++) {
51			double N = r.pval[i][4*j];
52			double mu_x = r.pval[i][4*j + 1];
53			double mu_y = r.pval[i][4*j + 2];
54			double sig = r.pval[i][4*j + 3];
55	dnisson	1.26
56			double rel_x = x - mu_x; double rel_y = y - mu_y;
57			fitVal += (N / 2.0 / M_PI / sig / sig)
58			* exp(-(rel_x * rel_x + rel_y * rel_y)/2.0/sig/sig);
59	dnisson	1.13	}
60	dnisson	1.26	return fitVal;
61			}
62	dnisson	1.13
63	dnisson	1.26	void JetFinderAnalyzer::analyze_results(HistoFitter::FitResults r,
64			std::vector<HistoFitter::Trouble> t,
65	dnisson	1.30	TH2 *hist_orig, const edm::Event &evt) {
66	dnisson	1.26	// perform analysis of fit results
67	dnisson	1.15	edm::Service<TFileService> fs;
68	dnisson	1.28	for (unsigned i = 0; i < r.pars.size(); i++) {
69			ostringstream fitHistoOss;
70			fitHistoOss << hist_orig->GetName() << "_fit" << i << flush;
71			TH2D *fitHisto = fs->make<TH2D>(fitHistoOss.str().c_str(),
72	dnisson	1.29	("Fitted distribution to "
73	dnisson	1.28	+std::string(hist_orig->GetName())).c_str(),
74			hist_orig->GetNbinsX(),
75			hist_orig->GetXaxis()->GetXmin(),
76			hist_orig->GetXaxis()->GetXmax(),
77			hist_orig->GetNbinsY(),
78			hist_orig->GetYaxis()->GetXmin(),
79			hist_orig->GetYaxis()->GetXmax());
80
81			double Xlo = fitHisto->GetXaxis()->GetXmin();
82			double Xhi = fitHisto->GetXaxis()->GetXmax();
83			double Ylo = fitHisto->GetYaxis()->GetXmin();
84			double Yhi = fitHisto->GetYaxis()->GetXmax();
85			double XbinSize = (Xhi - Xlo) / static_cast<double>(fitHisto->GetNbinsX());
86			double YbinSize = (Yhi - Ylo) / static_cast<double>(fitHisto->GetNbinsY());
87
88			for (int ib = 1; ib <= fitHisto->GetNbinsX(); ib++) {
89			for (int jb = 1; jb <= fitHisto->GetNbinsY(); jb++) {
90			double x = (static_cast<double>(ib) - 0.5) * XbinSize + Xlo;
91			double y = (static_cast<double>(jb) - 0.5) * YbinSize + Ylo;
92			fitHisto->SetBinContent(ib, jb, evalFitFunction(r, x, y, i) * XbinSize * YbinSize);
93			}
94			}
95
96			// save fit results to an ntuple
97			ostringstream fitResultsOss;
98			fitResultsOss << hist_orig->GetName() << "_results" << i << flush;
99			TNtuple *rNtuple = fs->make<TNtuple>(fitResultsOss.str().c_str(),
100			("Fit results for "+std::string(hist_orig->GetName())).c_str(),
101			"N:mu_x:mu_y:sigma");
102			unsigned nGauss = r.pval[i].size() / 4;
103			for (unsigned j = 0; j < nGauss; j++) {
104			rNtuple->Fill(r.pval[i][4j], r.pval[i][4j+1], r.pval[i][4*j+2],
105			r.pval[i][4*j+3]);
106			}
107
108	dnisson	1.29	// save jets found
109			vector<reco::Jet::LorentzVector> jets;
110			for (unsigned j = 0; j < nGauss; j++) {
111			double energy = r.pval[i][4*j];
112			double eta = r.pval[i][4*j + 1];
113			double phi = r.pval[i][4*j + 2];
114			double width = r.pval[i][4*j + 3];
115			//true jet width
116			width = sqrt(width*width - 0.0025);
117
118			reco::LeafCandidate::LorentzVector P4;
119			reco::LeafCandidate::Point vertex(0.0, 0.0, 0.0);
120			P4.SetE(energy);
121			// SetEta and SetPhi don't seem to be working
122			//P4.SetEta(eta);
123			//P4.SetPhi(phi);
124			double theta = 2.0*atan(exp(-eta));
125			P4.SetPx(energysin(theta)cos(phi));
126			P4.SetPy(energysin(theta)sin(phi));
127			P4.SetPz(energy*cos(theta));
128			jets.push_back(P4);
129			}
130			// compute pair masses
131			for (unsigned s1 = 0; s1 < jets.size() - 1; s1++) {
132			for (unsigned s2 = s1 + 1; s2 < jets.size(); s2++) {
133			// create TOTAL MOMENTUM VECTOR
134			reco::Jet::LorentzVector P4 = jets[s1] + jets[s2];
135			pairMassHist->Fill(static_cast<double>(P4.M()));
136			}
137			}
138
139	dnisson	1.28	// save chisquares to ntuple
140			for (unsigned j = 0; j < r.chisquare.size(); j++) {
141			ostringstream csNtupleName, csNtupleTitle;
142			csNtupleName << hist_orig->GetName() << "_chi2_" << j << flush;
143			csNtupleTitle << "Chisquare "<<j<<" for histo "<<hist_orig->GetName()
144			<< flush;
145			TNtuple *csNtuple = fs->make<TNtuple>(csNtupleName.str().c_str(),
146			csNtupleTitle.str().c_str(),
147			"chisq");
148			csNtuple->Fill(r.chisquare[j]);
149	dnisson	1.1	}
150	dnisson	1.29
151			// save jets found
152
153	dnisson	1.1	}
154			}
155
156			DEFINE_FWK_MODULE(JetFinderAnalyzer);