JetFitAnalyzer/test/JetFinderAnalyzer.cc

/* A JetFitAnalyzer that makes histograms with smearing */

#include "UserCode/JetFitAnalyzer/interface/JetFitAnalyzer.h"

#include "fastjet/ClusterSequence.hh"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "PhysicsTools/UtilAlgos/interface/TFileService.h"

#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
#include "DataFormats/ParticleFlowReco/interface/PFBlock.h"

/* Jet reco stuff */
#include "DataFormats/JetReco/interface/PFJetCollection.h"
#include "DataFormats/JetReco/interface/GenJetCollection.h"
#include "DataFormats/JetReco/interface/PFJet.h"
#include "DataFormats/JetReco/interface/GenJet.h"

#include <map>
#include <vector>
#include <limits>
#include <cmath>
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <sstream>

#include "TFormula.h"
#include "TF2.h"
#include "TNtuple.h"

#define PI 3.141593

using namespace std;
using namespace fastjet;

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

private:
  virtual void beginJob(const edm::EventSetup &es);
  virtual TH2D * make_histo(const edm::Event &, const edm::EventSetup &);
  virtual HistoFitter::ModelDefinition& make_model_def(const edm::Event&,
                                           const edm::EventSetup&,
                                           TH2 *);
  virtual void analyze_results(HistoFitter::FitResults r, 
                               std::vector<HistoFitter::Trouble> t, TH2 *);
  vector<reco::Candidate *> get_particles(const edm::Event&);

  fstream ofs;
  edm::InputTag inputTagPFCandidates_;
  edm::InputTag inputTagGenParticles_;
  int info_type_;
  double smear_;
  int smear_coord_;
  string jet_algo_;
  string reclustered_jet_algo_;
};

JetFinderAnalyzer::JetFinderAnalyzer(const edm::ParameterSet &pSet) 
  : JetFitAnalyzer(pSet) // this is important!
{
  info_type_ = pSet.getUntrackedParameter("info_type", 0);

  if (info_type_ == 0) {
    inputTagGenParticles_ = pSet.getParameter<edm::InputTag>("GenParticles");
  }
  if (info_type_ == 1) {
    inputTagPFCandidates_ = pSet.getParameter<edm::InputTag>("PFCandidates");
  }

  smear_ = pSet.getUntrackedParameter("smear", 0.02);
  smear_coord_ = pSet.getUntrackedParameter("smear_coord", 0);
  // 0 = eta-phi smear
  // 1 = proper angle smear
  jet_algo_ = pSet.getParameter<string>("jet_algo");
}

TH2D * JetFinderAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&) {
  const reco::Jet * highest_e_jet = 0;

  if (info_type_ == 0) {
    edm::Handle< vector<reco::GenJet> > evtJets;
    evt.getByLabel(jet_algo_, evtJets);
    for (unsigned i = 0; i < evtJets->size(); i++) {
      if (highest_e_jet == 0
          || (*evtJets)[i].energy() > highest_e_jet->energy()) {
        highest_e_jet = &((*evtJets)[i]);
      }
    }
  }
  else if (info_type_ == 1) {
    edm::Handle< vector<reco::PFJet> > evtJets;
    evt.getByLabel(jet_algo_, evtJets);
    for (unsigned i = 0; i < evtJets->size(); i++) {
      if (highest_e_jet == 0
          || (*evtJets)[i].energy() > highest_e_jet->energy()) {
        highest_e_jet = &((*evtJets)[i]);
      }
    }
  }
    
  if (highest_e_jet == 0) {
    cerr << "No fat jets found!" << endl;
    return 0;
  }

  vector<const reco::Candidate *> particles = 
    highest_e_jet->getJetConstituentsQuick();
 
  ostringstream oss;
  oss << "eta_phi_energy"<<evt.id().event() << flush;
  edm::Service<TFileService> fs;
  TH2D *histo = fs->make<TH2D>(oss.str().c_str(), oss.str().c_str(),
                         30,
                         highest_e_jet->eta()-0.5,
                         highest_e_jet->eta()+0.5,
                         30,
                         highest_e_jet->phi()-0.5,
                         highest_e_jet->phi()+0.5);

  if (smear_ > 0.0) {
    // create a temporary 2D vector for smeared energies
    double XbinSize = (histo->GetXaxis()->GetXmax()
                       - histo->GetXaxis()->GetXmin()) /
      static_cast<double>(histo->GetXaxis()->GetNbins());
    double YbinSize = (histo->GetYaxis()->GetXmax()
                       - histo->GetYaxis()->GetXmin()) /
      static_cast<double>(histo->GetYaxis()->GetNbins());
    double Xlo = histo->GetXaxis()->GetXmin();
    double Xhi = histo->GetXaxis()->GetXmax();
    double Ylo = histo->GetYaxis()->GetXmin();
    double Yhi = histo->GetYaxis()->GetXmax();
    vector< vector<double> > smeared(30, vector<double>(30, 0.0) );
    for (int i = 0; i < particles.size(); i++) {
      double N = particles[i]->energy();
      double x = particles[i]->eta();
      double y = particles[i]->phi();
      if (y < Ylo) {
        y += 2.0*PI;
      }
      if (y > Yhi) {
        y -= 2.0*PI;
      }
      // loop over bins and add Gaussian in proper angle to smeared
      for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) {
        for (vector< double >::size_type j2 = 0; j2 < 30; j2++) {
          double eta = static_cast<double>((signed int)i2) * XbinSize + 
            Xlo - x;
          double phi = static_cast<double>((signed int)j2) * YbinSize +
            Ylo - y;
          
          double theta, iota;
          switch (smear_coord_) {
          case 1:
            // transform eta, phi to proper angle
            theta = 2.0*atan(exp(-eta));
            iota = asin(sin(theta)*sin(phi));
            
            // add smeared Gaussian
            smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
              * exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_));
            break;
          case 0:
          default:
            // eta-phi smear--just add Gaussian centered around particle'
            // in eta and phi
            smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
              * exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_));
          }
        }
      }
      // set histogram to match smear vector
      for (int i = 1; i <= 30; i++) {
        for (int j = 1; j <= 30; j++) {
          histo->SetBinContent(i, j, smeared[i-1][j-1]);
        }
      }
    }
  }
  else {
    // don't smear--just fill with particles
    for (int i = 0; i < particles.size(); i++) {
      histo->Fill(particles[i]->eta(), particles[i]->phi(),
                  particles[i]->energy());
    }
  }
  
  return histo;
}
  
template <class Jet>
void seed_with_jetcoll(vector<Jet> jets,
                  HistoFitter::ModelDefinition &_mdef, double phi_lo = -PI,
                       double phi_hi = PI) {
  // seed with jet collection
  double N = jets[0].energy();
  double eta = jets[0].eta();
  double phi = jets[0].phi();
  _mdef.setIndivParameter(0, string("N"), N, 0.02*N,
                          0.0, 1.0e6);
  _mdef.setIndivParameter(1, string("mu_x"), eta, 0.01,
                          0.0, 0.0);
  if (phi < phi_lo) {
    phi += 2.0*PI;
  }
  if (phi > phi_hi) {
    phi -= 2.0*PI;
  }
  _mdef.setIndivParameter(2, string("mu_y"), phi, 0.01,
                          0.0, 0.0);
}

HistoFitter::ModelDefinition& JetFinderAnalyzer::make_model_def(const edm::Event& evt,
                                                 const edm::EventSetup&,
                                                 TH2 *histo) {
  class jf_model_def : public HistoFitter::ModelDefinition {
  public:
    virtual double chisquareSigma(double E) {
      return 0.298 - 0.9557*E; // study 11-09-09
    }
  };

  jf_model_def *_mdef = new jf_model_def();
  TFormula *formula = new TFormula("gaus2d", 
                                     "[0]*exp(-0.5*((x-[1])**2 + (y-[2])**2)/([3]**2))/(2*pi*[3]**2)");
  _mdef->setFormula(formula);
  _mdef->setIndivEnergy(0);
  _mdef->setIndivMeanX(1);
  _mdef->setIndivMeanY(2);
  _mdef->setIndivParameter(0, string("N"), 0.0, 0.0, 0.0, 1.0e6);
  _mdef->setIndivParameter(1, string("mu_x"), 0.0, 0.0, 0.0, 0.0);
  _mdef->setIndivParameter(2, string("mu_y"), 0.0, 0.0, 0.0, 0.0);
  _mdef->setIndivParameter(3, string("sig"), smear_, 0.001, 0.0, 0.0);

  _mdef->setMaxNumberOfGaussians(7);

  /*
  // get jetcoll from event file and select highest e jet
  if (info_type_ == 0) {
    edm::Handle< vector<reco::GenJet> > jet_collection;
    evt.getByLabel(jet_algo_, jet_collection);
    reco::GenJet highest_e_jet;
    bool found_jet = false;
    for (unsigned i = 0; i < jet_collection->size(); i++) {
      if (!found_jet || (*jet_collection)[i].energy() > highest_e_jet.energy())
        {
          highest_e_jet = (*jet_collection)[i];
          found_jet = true;
        }
    }
    vector<reco::GenJet> highest_e_jet_coll;
    highest_e_jet_coll.push_back(highest_e_jet);
    seed_with_jetcoll(highest_e_jet_coll, *_mdef,
                      histo->GetYaxis()->GetXmin(),
                      histo->GetYaxis()->GetXmax());
  }
  if (info_type_ == 1) {
    edm::Handle< vector<reco::PFJet> > jet_collection;
    evt.getByLabel(jet_algo_, jet_collection);
    reco::PFJet highest_e_jet;
    bool found_jet = false;
    for (unsigned i = 0; i < jet_collection->size(); i++) {
      if (!found_jet || (*jet_collection)[i].energy() > highest_e_jet.energy())
        {
          highest_e_jet = (*jet_collection)[i];
          found_jet = true;
        }
    }
    vector<reco::PFJet> highest_e_jet_coll;
    highest_e_jet_coll.push_back(highest_e_jet);
    seed_with_jetcoll(highest_e_jet_coll, *_mdef,
                      histo->GetYaxis()->GetXmin(),
                      histo->GetYaxis()->GetXmax());
  }
  */

  return *_mdef;
}

void JetFinderAnalyzer::beginJob(const edm::EventSetup &es) {

}

double evalFitFunction(HistoFitter::FitResults r, double x, double y) {
  unsigned nFits = r.pars.size();
  unsigned nGauss = r.pars[nFits-1].size() / 4;
  double fitVal = 0.0;
  for (unsigned i = 0; i < nGauss; i++) {
    double N = r.pval[nFits-1][4*i];
    double mu_x = r.pval[nFits-1][4*i + 1];
    double mu_y = r.pval[nFits-1][4*i + 2];
    double sig = r.pval[nFits-1][4*i + 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) {
  // TODO perform analysis of fit results
  edm::Service<TFileService> fs;
  TH2D *fitHisto = fs->make<TH2D>((std::string(hist_orig->GetName())+"_fit").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 i = 1; i <= fitHisto->GetNbinsX(); i++) {
    for (int j = 1; j <= fitHisto->GetNbinsY(); j++) {
      double x = (static_cast<double>(i) - 0.5) * XbinSize + Xlo;
      double y = (static_cast<double>(j) - 0.5) * YbinSize + Ylo;
      fitHisto->SetBinContent(i, j, evalFitFunction(r, x, y) * XbinSize * YbinSize);
    }
  }

  // fit result ntuple
  TNtuple *rNtuple = fs->make<TNtuple>((std::string(hist_orig->GetName())+"_results").c_str(),
                                       ("Fit results for "+std::string(hist_orig->GetName())).c_str(),
                                       "N:mu_x:mu_y:sigma");
  unsigned nFits = r.pval.size();
  unsigned nGauss = r.pval[nFits-1].size() / 4;
  for (unsigned i = 0; i < nGauss; i++) {
    rNtuple->Fill(r.pval[nFits-1][4*i], r.pval[nFits-1][4*i+1], r.pval[nFits-1][4*i+2],
                  r.pval[nFits-1][4*i+3]);
  }
}

DEFINE_FWK_MODULE(JetFinderAnalyzer);
Revision:	1.25.2.3
Committed:	Sun Dec 20 20:36:09 2009 UTC (15 years, 4 months ago) by dnisson
Content type:	text/plain
Branch:	V01-00-01-SPECTRUM
Changes since 1.25.2.2:	+12 -0 lines
Log Message:	Saving ntuple
#	User	Rev	Content
1	dnisson	1.24	/* A JetFitAnalyzer that makes histograms with smearing */
2	dnisson	1.1
3			#include "UserCode/JetFitAnalyzer/interface/JetFitAnalyzer.h"
4
5			#include "fastjet/ClusterSequence.hh"
6			#include "FWCore/ServiceRegistry/interface/Service.h"
7	dnisson	1.2	#include "FWCore/MessageLogger/interface/MessageLogger.h"
8	dnisson	1.1	#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
9
10	dnisson	1.4	#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
11			#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
12	dnisson	1.2	#include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
13			#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
14			#include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
15
16	dnisson	1.11	/* Jet reco stuff */
17			#include "DataFormats/JetReco/interface/PFJetCollection.h"
18			#include "DataFormats/JetReco/interface/GenJetCollection.h"
19			#include "DataFormats/JetReco/interface/PFJet.h"
20			#include "DataFormats/JetReco/interface/GenJet.h"
21
22	dnisson	1.1	#include <map>
23			#include <vector>
24			#include <limits>
25			#include <cmath>
26			#include <cstdlib>
27	dnisson	1.18	#include <iostream>
28	dnisson	1.1	#include <fstream>
29			#include <sstream>
30
31			#include "TFormula.h"
32			#include "TF2.h"
33	dnisson	1.25.2.3	#include "TNtuple.h"
34	dnisson	1.1
35			#define PI 3.141593
36
37			using namespace std;
38			using namespace fastjet;
39
40			class JetFinderAnalyzer : public JetFitAnalyzer {
41			public:
42			explicit JetFinderAnalyzer( const edm::ParameterSet&);
43			~JetFinderAnalyzer() {}
44
45			private:
46	dnisson	1.15	virtual void beginJob(const edm::EventSetup &es);
47			virtual TH2D * make_histo(const edm::Event &, const edm::EventSetup &);
48	dnisson	1.23	virtual HistoFitter::ModelDefinition& make_model_def(const edm::Event&,
49	dnisson	1.1	const edm::EventSetup&,
50			TH2 *);
51	dnisson	1.22	virtual void analyze_results(HistoFitter::FitResults r,
52			std::vector<HistoFitter::Trouble> t, TH2 *);
53	dnisson	1.4	vector<reco::Candidate *> get_particles(const edm::Event&);
54	dnisson	1.1
55			fstream ofs;
56	dnisson	1.2	edm::InputTag inputTagPFCandidates_;
57	dnisson	1.4	edm::InputTag inputTagGenParticles_;
58	dnisson	1.2	int info_type_;
59	dnisson	1.1	double smear_;
60			int smear_coord_;
61	dnisson	1.11	string jet_algo_;
62	dnisson	1.15	string reclustered_jet_algo_;
63	dnisson	1.1	};
64
65			JetFinderAnalyzer::JetFinderAnalyzer(const edm::ParameterSet &pSet)
66			: JetFitAnalyzer(pSet) // this is important!
67			{
68	dnisson	1.2	info_type_ = pSet.getUntrackedParameter("info_type", 0);
69
70	dnisson	1.4	if (info_type_ == 0) {
71			inputTagGenParticles_ = pSet.getParameter<edm::InputTag>("GenParticles");
72			}
73	dnisson	1.2	if (info_type_ == 1) {
74			inputTagPFCandidates_ = pSet.getParameter<edm::InputTag>("PFCandidates");
75			}
76
77	dnisson	1.1	smear_ = pSet.getUntrackedParameter("smear", 0.02);
78			smear_coord_ = pSet.getUntrackedParameter("smear_coord", 0);
79			// 0 = eta-phi smear
80			// 1 = proper angle smear
81	dnisson	1.11	jet_algo_ = pSet.getParameter<string>("jet_algo");
82	dnisson	1.1	}
83
84			TH2D * JetFinderAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&) {
85	dnisson	1.16	const reco::Jet * highest_e_jet = 0;
86	dnisson	1.1
87	dnisson	1.13	if (info_type_ == 0) {
88			edm::Handle< vector<reco::GenJet> > evtJets;
89			evt.getByLabel(jet_algo_, evtJets);
90			for (unsigned i = 0; i < evtJets->size(); i++) {
91	dnisson	1.16	if (highest_e_jet == 0
92			\|\| (*evtJets)[i].energy() > highest_e_jet->energy()) {
93	dnisson	1.13	highest_e_jet = &((*evtJets)[i]);
94			}
95			}
96	dnisson	1.1	}
97	dnisson	1.13	else if (info_type_ == 1) {
98			edm::Handle< vector<reco::PFJet> > evtJets;
99			evt.getByLabel(jet_algo_, evtJets);
100			for (unsigned i = 0; i < evtJets->size(); i++) {
101	dnisson	1.16	if (highest_e_jet == 0
102			\|\| (*evtJets)[i].energy() > highest_e_jet->energy()) {
103	dnisson	1.13	highest_e_jet = &((*evtJets)[i]);
104	dnisson	1.1	}
105			}
106			}
107	dnisson	1.13
108			if (highest_e_jet == 0) {
109	dnisson	1.17	cerr << "No fat jets found!" << endl;
110			return 0;
111	dnisson	1.13	}
112
113			vector<const reco::Candidate *> particles =
114			highest_e_jet->getJetConstituentsQuick();
115	dnisson	1.17
116	dnisson	1.13	ostringstream oss;
117			oss << "eta_phi_energy"<<evt.id().event() << flush;
118	dnisson	1.15	edm::Service<TFileService> fs;
119			TH2D *histo = fs->make<TH2D>(oss.str().c_str(), oss.str().c_str(),
120	dnisson	1.14	30,
121	dnisson	1.13	highest_e_jet->eta()-0.5,
122			highest_e_jet->eta()+0.5,
123	dnisson	1.14	30,
124	dnisson	1.13	highest_e_jet->phi()-0.5,
125			highest_e_jet->phi()+0.5);
126
127	dnisson	1.17	if (smear_ > 0.0) {
128	dnisson	1.24	// create a temporary 2D vector for smeared energies
129			double XbinSize = (histo->GetXaxis()->GetXmax()
130			- histo->GetXaxis()->GetXmin()) /
131			static_cast<double>(histo->GetXaxis()->GetNbins());
132			double YbinSize = (histo->GetYaxis()->GetXmax()
133			- histo->GetYaxis()->GetXmin()) /
134			static_cast<double>(histo->GetYaxis()->GetNbins());
135			double Xlo = histo->GetXaxis()->GetXmin();
136			double Xhi = histo->GetXaxis()->GetXmax();
137			double Ylo = histo->GetYaxis()->GetXmin();
138			double Yhi = histo->GetYaxis()->GetXmax();
139			vector< vector<double> > smeared(30, vector<double>(30, 0.0) );
140	dnisson	1.1	for (int i = 0; i < particles.size(); i++) {
141	dnisson	1.4	double N = particles[i]->energy();
142	dnisson	1.2	double x = particles[i]->eta();
143			double y = particles[i]->phi();
144	dnisson	1.15	if (y < Ylo) {
145			y += 2.0*PI;
146			}
147			if (y > Yhi) {
148			y -= 2.0*PI;
149			}
150	dnisson	1.1	// loop over bins and add Gaussian in proper angle to smeared
151	dnisson	1.14	for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) {
152			for (vector< double >::size_type j2 = 0; j2 < 30; j2++) {
153	dnisson	1.1	double eta = static_cast<double>((signed int)i2) * XbinSize +
154	dnisson	1.15	Xlo - x;
155			double phi = static_cast<double>((signed int)j2) * YbinSize +
156			Ylo - y;
157	dnisson	1.1
158	dnisson	1.24	double theta, iota;
159			switch (smear_coord_) {
160			case 1:
161			// transform eta, phi to proper angle
162			theta = 2.0*atan(exp(-eta));
163			iota = asin(sin(theta)*sin(phi));
164
165			// add smeared Gaussian
166			smeared[i2][j2] += (NXbinSizeYbinSize/(2.0PIsmear_*smear_))
167			* exp(-0.5(thetatheta + iotaiota)/(smear_smear_));
168			break;
169			case 0:
170			default:
171			// eta-phi smear--just add Gaussian centered around particle'
172			// in eta and phi
173			smeared[i2][j2] += (NXbinSizeYbinSize/(2.0PIsmear_*smear_))
174			* exp(-0.5(etaeta + phiphi)/(smear_smear_));
175			}
176	dnisson	1.1	}
177			}
178	dnisson	1.24	// set histogram to match smear vector
179			for (int i = 1; i <= 30; i++) {
180			for (int j = 1; j <= 30; j++) {
181			histo->SetBinContent(i, j, smeared[i-1][j-1]);
182	dnisson	1.1	}
183			}
184			}
185			}
186	dnisson	1.17	else {
187			// don't smear--just fill with particles
188			for (int i = 0; i < particles.size(); i++) {
189			histo->Fill(particles[i]->eta(), particles[i]->phi(),
190			particles[i]->energy());
191			}
192			}
193	dnisson	1.24
194	dnisson	1.1	return histo;
195			}
196	dnisson	1.24
197	dnisson	1.12	template <class Jet>
198			void seed_with_jetcoll(vector<Jet> jets,
199	dnisson	1.22	HistoFitter::ModelDefinition &_mdef, double phi_lo = -PI,
200	dnisson	1.15	double phi_hi = PI) {
201	dnisson	1.11	// seed with jet collection
202	dnisson	1.23	double N = jets[0].energy();
203			double eta = jets[0].eta();
204			double phi = jets[0].phi();
205	dnisson	1.25.2.1	_mdef.setIndivParameter(0, string("N"), N, 0.02*N,
206	dnisson	1.23	0.0, 1.0e6);
207	dnisson	1.25.2.1	_mdef.setIndivParameter(1, string("mu_x"), eta, 0.01,
208	dnisson	1.23	0.0, 0.0);
209			if (phi < phi_lo) {
210			phi += 2.0*PI;
211			}
212			if (phi > phi_hi) {
213			phi -= 2.0*PI;
214			}
215	dnisson	1.25.2.1	_mdef.setIndivParameter(2, string("mu_y"), phi, 0.01,
216	dnisson	1.23	0.0, 0.0);
217	dnisson	1.1	}
218
219	dnisson	1.22	HistoFitter::ModelDefinition& JetFinderAnalyzer::make_model_def(const edm::Event& evt,
220	dnisson	1.1	const edm::EventSetup&,
221			TH2 *histo) {
222	dnisson	1.22	class jf_model_def : public HistoFitter::ModelDefinition {
223	dnisson	1.1	public:
224	dnisson	1.22	virtual double chisquareSigma(double E) {
225	dnisson	1.19	return 0.298 - 0.9557*E; // study 11-09-09
226	dnisson	1.1	}
227			};
228
229			jf_model_def *_mdef = new jf_model_def();
230			TFormula *formula = new TFormula("gaus2d",
231			"[0]exp(-0.5((x-[1])2 + (y-[2])2)/([3]*2))/(2pi[3]*2)");
232	dnisson	1.23	_mdef->setFormula(formula);
233	dnisson	1.22	_mdef->setIndivEnergy(0);
234			_mdef->setIndivMeanX(1);
235			_mdef->setIndivMeanY(2);
236			_mdef->setIndivParameter(0, string("N"), 0.0, 0.0, 0.0, 1.0e6);
237			_mdef->setIndivParameter(1, string("mu_x"), 0.0, 0.0, 0.0, 0.0);
238			_mdef->setIndivParameter(2, string("mu_y"), 0.0, 0.0, 0.0, 0.0);
239	dnisson	1.25.2.1	_mdef->setIndivParameter(3, string("sig"), smear_, 0.001, 0.0, 0.0);
240	dnisson	1.22
241	dnisson	1.25.2.1	_mdef->setMaxNumberOfGaussians(7);
242	dnisson	1.1
243	dnisson	1.25.2.1	/*
244	dnisson	1.13	// get jetcoll from event file and select highest e jet
245	dnisson	1.12	if (info_type_ == 0) {
246			edm::Handle< vector<reco::GenJet> > jet_collection;
247			evt.getByLabel(jet_algo_, jet_collection);
248	dnisson	1.13	reco::GenJet highest_e_jet;
249			bool found_jet = false;
250			for (unsigned i = 0; i < jet_collection->size(); i++) {
251			if (!found_jet \|\| (*jet_collection)[i].energy() > highest_e_jet.energy())
252			{
253			highest_e_jet = (*jet_collection)[i];
254	dnisson	1.17	found_jet = true;
255	dnisson	1.13	}
256			}
257			vector<reco::GenJet> highest_e_jet_coll;
258			highest_e_jet_coll.push_back(highest_e_jet);
259	dnisson	1.15	seed_with_jetcoll(highest_e_jet_coll, *_mdef,
260			histo->GetYaxis()->GetXmin(),
261			histo->GetYaxis()->GetXmax());
262	dnisson	1.12	}
263			if (info_type_ == 1) {
264			edm::Handle< vector<reco::PFJet> > jet_collection;
265			evt.getByLabel(jet_algo_, jet_collection);
266	dnisson	1.13	reco::PFJet highest_e_jet;
267			bool found_jet = false;
268			for (unsigned i = 0; i < jet_collection->size(); i++) {
269			if (!found_jet \|\| (*jet_collection)[i].energy() > highest_e_jet.energy())
270			{
271			highest_e_jet = (*jet_collection)[i];
272	dnisson	1.17	found_jet = true;
273	dnisson	1.13	}
274			}
275			vector<reco::PFJet> highest_e_jet_coll;
276			highest_e_jet_coll.push_back(highest_e_jet);
277	dnisson	1.15	seed_with_jetcoll(highest_e_jet_coll, *_mdef,
278			histo->GetYaxis()->GetXmin(),
279			histo->GetYaxis()->GetXmax());
280	dnisson	1.12	}
281	dnisson	1.25.2.1	*/
282	dnisson	1.1
283			return *_mdef;
284			}
285
286			void JetFinderAnalyzer::beginJob(const edm::EventSetup &es) {
287
288			}
289
290	dnisson	1.25.2.2	double evalFitFunction(HistoFitter::FitResults r, double x, double y) {
291			unsigned nFits = r.pars.size();
292			unsigned nGauss = r.pars[nFits-1].size() / 4;
293			double fitVal = 0.0;
294			for (unsigned i = 0; i < nGauss; i++) {
295			double N = r.pval[nFits-1][4*i];
296			double mu_x = r.pval[nFits-1][4*i + 1];
297			double mu_y = r.pval[nFits-1][4*i + 2];
298			double sig = r.pval[nFits-1][4*i + 3];
299
300			double rel_x = x - mu_x; double rel_y = y - mu_y;
301			fitVal += (N / 2.0 / M_PI / sig / sig)
302			* exp(-(rel_x * rel_x + rel_y * rel_y)/2.0/sig/sig);
303			}
304			return fitVal;
305			}
306
307	dnisson	1.23	void JetFinderAnalyzer::analyze_results(HistoFitter::FitResults r,
308			std::vector<HistoFitter::Trouble> t,
309	dnisson	1.1	TH2 *hist_orig) {
310	dnisson	1.24	// TODO perform analysis of fit results
311	dnisson	1.25.2.2	edm::Service<TFileService> fs;
312			TH2D *fitHisto = fs->make<TH2D>((std::string(hist_orig->GetName())+"_fit").c_str(),
313			("Fitted distribution to "
314			+std::string(hist_orig->GetName())).c_str(),
315			hist_orig->GetNbinsX(),
316			hist_orig->GetXaxis()->GetXmin(),
317			hist_orig->GetXaxis()->GetXmax(),
318			hist_orig->GetNbinsY(),
319			hist_orig->GetYaxis()->GetXmin(),
320			hist_orig->GetYaxis()->GetXmax());
321
322			double Xlo = fitHisto->GetXaxis()->GetXmin();
323			double Xhi = fitHisto->GetXaxis()->GetXmax();
324			double Ylo = fitHisto->GetYaxis()->GetXmin();
325			double Yhi = fitHisto->GetYaxis()->GetXmax();
326			double XbinSize = (Xhi - Xlo) / static_cast<double>(fitHisto->GetNbinsX());
327			double YbinSize = (Yhi - Ylo) / static_cast<double>(fitHisto->GetNbinsY());
328
329			for (int i = 1; i <= fitHisto->GetNbinsX(); i++) {
330			for (int j = 1; j <= fitHisto->GetNbinsY(); j++) {
331			double x = (static_cast<double>(i) - 0.5) * XbinSize + Xlo;
332			double y = (static_cast<double>(j) - 0.5) * YbinSize + Ylo;
333			fitHisto->SetBinContent(i, j, evalFitFunction(r, x, y) * XbinSize * YbinSize);
334			}
335			}
336	dnisson	1.25.2.3
337			// fit result ntuple
338			TNtuple *rNtuple = fs->make<TNtuple>((std::string(hist_orig->GetName())+"_results").c_str(),
339			("Fit results for "+std::string(hist_orig->GetName())).c_str(),
340			"N:mu_x:mu_y:sigma");
341			unsigned nFits = r.pval.size();
342			unsigned nGauss = r.pval[nFits-1].size() / 4;
343			for (unsigned i = 0; i < nGauss; i++) {
344			rNtuple->Fill(r.pval[nFits-1][4i], r.pval[nFits-1][4i+1], r.pval[nFits-1][4*i+2],
345			r.pval[nFits-1][4*i+3]);
346			}
347	dnisson	1.1	}
348
349			DEFINE_FWK_MODULE(JetFinderAnalyzer);