JetFitAnalyzer/test/JetFinderAnalyzer.cc

/* A JetFitAnalyze_hat makes histograms with smearing */

#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 *);

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