JetFitAnalyzer/test/JetFinderAnalyzer.cc

/* A simple jet-finding analyzer */

#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"

#define PI 3.141593

using namespace std;
using namespace fastjet;

class JetFinderAnalyzer : public JetFitAnalyzer {
public:
  struct jet {
    double energy;
    double eta;
    double phi;
  };

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

private:
  static map<TH2 *, vector< vector<jet> > > unique_jets;

  static double phi_cutoff_;

  static double g2int(double xlo, double xhi, double ylo, double yhi,
               double *pval) {
    double sum1 = 0.0;
    double sum2 = 0.0;
    double xmid = 0.5 * (xlo + xhi);
    double ymid = 0.5 * (ylo + yhi);
    double xstep = (xhi - xlo) / 50.0;
    double ystep = (yhi - ylo) / 50.0;
    for (int i = 0; i < 50; i++) {
      double x = (static_cast<double>(i) + 0.5) * xstep + xlo;
      sum1 += xstep * jetfit::fit_fcn(x, ymid, pval);
    }
    for (int i = 0; i < 50; i++) {
      double y = (static_cast<double>(i) + 0.5) * ystep + ylo;
      sum2 += ystep * jetfit::fit_fcn(xmid, y, pval);
    }
    return sum1 * sum2;
  }

  virtual void beginJob(const edm::EventSetup &es);
  virtual TH2D * make_histo(const edm::Event &, const edm::EventSetup &);
  virtual jetfit::model_def& make_model_def(const edm::Event&,
                                           const edm::EventSetup&,
                                           TH2 *);
  virtual void analyze_results(jetfit::results r, 
                               std::vector<jetfit::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_;
};

map<TH2 *, vector< vector< JetFinderAnalyzer::jet > > >
JetFinderAnalyzer::unique_jets;

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");
  set_user_minuit(0); // we don't need this feature--it's deprecated
}

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 smeared histo
  // 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) );
  switch (smear_coord_) {
  case 1:
    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;

          // transform eta, phi to proper angle
          double theta = 2.0*atan(exp(-eta));
          double iota = asin(sin(theta)*sin(phi));
          
          smeared[i2][j2] += (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
            * exp(-0.5*(theta*theta + iota*iota)/(smear_*smear_));
        }
      }
    }
    break;
  case 0:
  default:
    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 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;

          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,
                  jetfit::model_def &_mdef, double phi_lo = -PI,
                       double phi_hi = PI) {
  // seed with jet collection
  int ijset = 0;
  for (unsigned ij = 0; ij < jets.size(); ij++) {
    double N = jets[ij].energy();
    double eta = jets[ij].eta();
    double phi = jets[ij].phi();
    if (N > 10.0) {
      _mdef.set_special_par(ijset, 0, N, _mdef.chisquare_error(N)*0.1,
                            0.0, 1.0e6);
      _mdef.set_special_par(ijset, 1, eta, 0.01,
                            0.0, 0.0);
      if (phi < phi_lo) {
        phi += 2.0*PI;
      }
      if (phi > phi_hi) {
        phi -= 2.0*PI;
      }
      _mdef.set_special_par(ijset, 2, phi, 0.01,
                            0.0, 0.0);
      _mdef.set_special_par(ijset, 3, 0.1, 0.001, 0.0, 0.0);
      ijset++;
    }
  }
}

jetfit::model_def& JetFinderAnalyzer::make_model_def(const edm::Event& evt,
                                                 const edm::EventSetup&,
                                                 TH2 *histo) {
  class jf_model_def : public jetfit::model_def {
  public:
    virtual double chisquare_error(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->set_formula(formula);
  _mdef->set_indiv_max_E(0);
  _mdef->set_indiv_max_x(1);
  _mdef->set_indiv_max_y(2);
  _mdef->set_indiv_par(0, string("N"), 0.0, 0.0, 0.0, 1.0e6);
  _mdef->set_indiv_par(1, string("mu_x"), 0.0, 0.0, 0.0, 0.0);
  _mdef->set_indiv_par(2, string("mu_y"), 0.0, 0.0, 0.0, 0.0);
  _mdef->set_indiv_par(3, string("sig"), 0.1, 0.001, 0.0, 0.0);

  // 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());
  }

  jetfit::set_model_def(_mdef);

  // generate initial fit histogram
  edm::Service<TFileService> fs;
  TH2D *init_fit_histo = fs->make<TH2D>(("init_fit_"+string(histo->GetName()))
                                        .c_str(),
                                        ("Initial fit for "
                                         +string(histo->GetName())).c_str(),
                                        histo->GetXaxis()->GetNbins(),
                                        histo->GetXaxis()->GetXmin(),
                                        histo->GetXaxis()->GetXmax(),
                                        histo->GetXaxis()->GetNbins(),
                                        histo->GetXaxis()->GetXmin(),
                                        histo->GetXaxis()->GetXmax());
  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();

  for (int i = 0; i < 30; i++) {
    for (int j = 0; j < 30; j++) {
      double x = (static_cast<double>(i) + 0.5)*XbinSize + Xlo;
      double y = (static_cast<double>(j) + 0.5)*YbinSize + Ylo;
      double pval[256];
      if (_mdef->get_n_special_par_sets() > 64) {
        cerr << "Parameter overload" << endl;
        return *_mdef;
      }
      else {
        for (int is = 0; is < _mdef->get_n_special_par_sets(); is++) {
          for (int ii = 0; ii < 4; ii++) {
            double spval, sperr, splo, sphi;
            _mdef->get_special_par(is, ii, spval, sperr, splo, sphi);
            pval[4*is + ii] = spval;
          }
        }
      }
      jetfit::set_ngauss(_mdef->get_n_special_par_sets());
      init_fit_histo->SetBinContent(i+1, j+1,
                                    jetfit::fit_fcn(x, y, pval)
                                    * XbinSize * YbinSize);
    }
  }

  return *_mdef;
}

void JetFinderAnalyzer::beginJob(const edm::EventSetup &es) {
  ofs.open("jetfindlog.txt", ios::out);
  if (ofs.fail()) {
    cerr << "Opening jetfindlog.txt FAILED" << endl;
  }
  ofs << "Jetfinder log" << endl
      << "=============" << endl << endl;
}

ostream& operator<<(ostream &out, jetfit::trouble t) {
  string action, error_string;
  
  if (t.istat != 3) {
    switch(t.occ) {
    case jetfit::T_NULL:
      action = "Program"; break;
    case jetfit::T_SIMPLEX:
      action = "SIMPLEX"; break;
    case jetfit::T_MIGRAD:
      action = "MIGRAD"; break;
    case jetfit::T_MINOS:
      action = "MINOS"; break;
    default:
      action = "Program"; break;
    }

    switch (t.istat) {
    case 0:
      error_string = "Unable to calculate error matrix"; break;
    case 1:
      error_string = "Error matrix a diagonal approximation"; break;
    case 2:
      error_string = "Error matrix not positive definite"; break;
    case 3:
      error_string = "Converged successfully"; break;
    default:
      ostringstream oss;
      oss<<"Unknown status code "<<t.istat << flush;
      error_string = oss.str(); break;
    }

    if (t.occ != jetfit::T_NULL)
      out << action<<" trouble: "<<error_string;
    else
      out << "Not calculated" << endl;
  }
  else {
    out << "Error matrix accurate" << endl;
  }

  return out;
}

void JetFinderAnalyzer::analyze_results(jetfit::results r,
                                     std::vector<jetfit::trouble> t,
                                     TH2 *hist_orig) {
  ofs << "Histogram "<<hist_orig->GetName() << endl;
  for (int i = unique_jets[hist_orig].size() - r.chisquare.size();
       i < unique_jets[hist_orig].size(); i++) {
    int ir = i - unique_jets[hist_orig].size() + r.chisquare.size();
    ofs << "For "<<i+1<<" gaussians: " << endl
        << t.at(i) << endl;
    ofs << "chisquare="<<r.chisquare.at(ir)
        << endl;
    ofs << unique_jets[hist_orig][i].size()<<" unique jets found" << endl;
    for (int j = 0; j < unique_jets[hist_orig][i].size(); j++) {
      jet _jet = unique_jets[hist_orig][i][j];
      ofs << "Jet "<<j<<": Energy = "<<_jet.energy<<", eta = "<<_jet.eta
          << ", phi = "<<_jet.phi << endl;
    }
    ofs << endl;
  }
  ofs << endl;

  // save fit function histograms to root file
  edm::Service<TFileService> fs;
  for (vector< vector<double> >::size_type i = 0;
       i < r.pval.size(); i++) {
    jetfit::set_ngauss(r.pval[i].size() / 4);
    TF2 *tf2 = new TF2("fit_func", jetfit::fit_fcn_TF2,
                       hist_orig->GetXaxis()->GetXmin(),
                       hist_orig->GetXaxis()->GetXmax(),
                       hist_orig->GetYaxis()->GetXmin(),
                       hist_orig->GetYaxis()->GetXmax(),
                       r.pval[i].size());
    for (vector<double>::size_type j = 0; j < r.pval[i].size(); j++) {
      tf2->SetParameter(j, r.pval[i][j]);
    }
    ostringstream fit_histo_oss;
    fit_histo_oss << hist_orig->GetName()<<"_fit_"<<i << flush;
    tf2->SetNpx(hist_orig->GetXaxis()->GetNbins());
    tf2->SetNpy(hist_orig->GetYaxis()->GetNbins());
    TH2D *fit_histo = fs->make<TH2D>(fit_histo_oss.str().c_str(),
                                     fit_histo_oss.str().c_str(),
                                     hist_orig->GetXaxis()->GetNbins(),
                                     hist_orig->GetXaxis()->GetXmin(),
                                     hist_orig->GetXaxis()->GetXmax(),
                                     hist_orig->GetYaxis()->GetNbins(),
                                     hist_orig->GetYaxis()->GetXmin(),
                                     hist_orig->GetYaxis()->GetXmax());
    double XbinSize = (hist_orig->GetXaxis()->GetXmax()
                       - hist_orig->GetXaxis()->GetXmin())
      / hist_orig->GetXaxis()->GetNbins();
    double YbinSize = (hist_orig->GetYaxis()->GetXmax()
                       - hist_orig->GetYaxis()->GetXmin())
      / hist_orig->GetYaxis()->GetNbins();
    for (int i = 1; i <= fit_histo->GetXaxis()->GetNbins(); i++) {
      for (int j = 1; j <= fit_histo->GetYaxis()->GetNbins(); j++) {
        double x = static_cast<double>(i - 1)
          * (fit_histo->GetXaxis()->GetXmax()
             - fit_histo->GetXaxis()->GetXmin())
          / static_cast<double>(fit_histo->GetXaxis()->GetNbins())
          + fit_histo->GetXaxis()->GetXmin();;
        double y = static_cast<double>(j - 1)
          * (fit_histo->GetYaxis()->GetXmax()
             - fit_histo->GetYaxis()->GetXmin())
          / static_cast<double>(fit_histo->GetYaxis()->GetNbins())
          + fit_histo->GetYaxis()->GetXmin();
        fit_histo->SetBinContent(i, j, tf2->Eval(x, y) * XbinSize * YbinSize);
      }
    }
  }

  // save results to file
  ostringstream res_tree_oss, rt_title_oss;
  res_tree_oss << hist_orig->GetName()<<"_results" << flush;
  rt_title_oss << "Fit results for "<<hist_orig->GetName() << flush;
}

DEFINE_FWK_MODULE(JetFinderAnalyzer);
Revision:	1.19
Committed:	Tue Nov 10 01:25:49 2009 UTC (15 years, 5 months ago) by dnisson
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V00-01-03, V00-01-02, V00-01-01
Branch point for:	V00-01-03-branch
Changes since 1.18:	+30 -23 lines
Log Message:	New sigmas
#	User	Rev	Content
1	dnisson	1.1	/* A simple jet-finding analyzer */
2
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
34			#define PI 3.141593
35
36			using namespace std;
37			using namespace fastjet;
38
39			class JetFinderAnalyzer : public JetFitAnalyzer {
40			public:
41			struct jet {
42			double energy;
43			double eta;
44			double phi;
45			};
46
47			explicit JetFinderAnalyzer( const edm::ParameterSet&);
48			~JetFinderAnalyzer() {}
49
50			private:
51			static map<TH2 *, vector< vector<jet> > > unique_jets;
52
53			static double phi_cutoff_;
54
55			static double g2int(double xlo, double xhi, double ylo, double yhi,
56			double *pval) {
57			double sum1 = 0.0;
58			double sum2 = 0.0;
59			double xmid = 0.5 * (xlo + xhi);
60			double ymid = 0.5 * (ylo + yhi);
61			double xstep = (xhi - xlo) / 50.0;
62			double ystep = (yhi - ylo) / 50.0;
63			for (int i = 0; i < 50; i++) {
64			double x = (static_cast<double>(i) + 0.5) * xstep + xlo;
65			sum1 += xstep * jetfit::fit_fcn(x, ymid, pval);
66			}
67			for (int i = 0; i < 50; i++) {
68			double y = (static_cast<double>(i) + 0.5) * ystep + ylo;
69			sum2 += ystep * jetfit::fit_fcn(xmid, y, pval);
70			}
71			return sum1 * sum2;
72			}
73
74	dnisson	1.15	virtual void beginJob(const edm::EventSetup &es);
75			virtual TH2D * make_histo(const edm::Event &, const edm::EventSetup &);
76	dnisson	1.1	virtual jetfit::model_def& make_model_def(const edm::Event&,
77			const edm::EventSetup&,
78			TH2 *);
79			virtual void analyze_results(jetfit::results r,
80			std::vector<jetfit::trouble> t, TH2 *);
81	dnisson	1.4	vector<reco::Candidate *> get_particles(const edm::Event&);
82	dnisson	1.1
83			fstream ofs;
84	dnisson	1.2	edm::InputTag inputTagPFCandidates_;
85	dnisson	1.4	edm::InputTag inputTagGenParticles_;
86	dnisson	1.2	int info_type_;
87	dnisson	1.1	double smear_;
88			int smear_coord_;
89	dnisson	1.11	string jet_algo_;
90	dnisson	1.15	string reclustered_jet_algo_;
91	dnisson	1.1	};
92
93			map<TH2 *, vector< vector< JetFinderAnalyzer::jet > > >
94			JetFinderAnalyzer::unique_jets;
95
96			JetFinderAnalyzer::JetFinderAnalyzer(const edm::ParameterSet &pSet)
97			: JetFitAnalyzer(pSet) // this is important!
98			{
99	dnisson	1.2	info_type_ = pSet.getUntrackedParameter("info_type", 0);
100
101	dnisson	1.4	if (info_type_ == 0) {
102			inputTagGenParticles_ = pSet.getParameter<edm::InputTag>("GenParticles");
103			}
104	dnisson	1.2	if (info_type_ == 1) {
105			inputTagPFCandidates_ = pSet.getParameter<edm::InputTag>("PFCandidates");
106			}
107
108	dnisson	1.1	smear_ = pSet.getUntrackedParameter("smear", 0.02);
109			smear_coord_ = pSet.getUntrackedParameter("smear_coord", 0);
110			// 0 = eta-phi smear
111			// 1 = proper angle smear
112	dnisson	1.11	jet_algo_ = pSet.getParameter<string>("jet_algo");
113	dnisson	1.15	set_user_minuit(0); // we don't need this feature--it's deprecated
114	dnisson	1.1	}
115
116			TH2D * JetFinderAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&) {
117	dnisson	1.16	const reco::Jet * highest_e_jet = 0;
118	dnisson	1.1
119	dnisson	1.13	if (info_type_ == 0) {
120			edm::Handle< vector<reco::GenJet> > evtJets;
121			evt.getByLabel(jet_algo_, evtJets);
122			for (unsigned i = 0; i < evtJets->size(); i++) {
123	dnisson	1.16	if (highest_e_jet == 0
124			\|\| (*evtJets)[i].energy() > highest_e_jet->energy()) {
125	dnisson	1.13	highest_e_jet = &((*evtJets)[i]);
126			}
127			}
128	dnisson	1.1	}
129	dnisson	1.13	else if (info_type_ == 1) {
130			edm::Handle< vector<reco::PFJet> > evtJets;
131			evt.getByLabel(jet_algo_, evtJets);
132			for (unsigned i = 0; i < evtJets->size(); i++) {
133	dnisson	1.16	if (highest_e_jet == 0
134			\|\| (*evtJets)[i].energy() > highest_e_jet->energy()) {
135	dnisson	1.13	highest_e_jet = &((*evtJets)[i]);
136	dnisson	1.1	}
137			}
138			}
139	dnisson	1.13
140			if (highest_e_jet == 0) {
141	dnisson	1.17	cerr << "No fat jets found!" << endl;
142			return 0;
143	dnisson	1.13	}
144
145			vector<const reco::Candidate *> particles =
146			highest_e_jet->getJetConstituentsQuick();
147	dnisson	1.17
148	dnisson	1.13	ostringstream oss;
149			oss << "eta_phi_energy"<<evt.id().event() << flush;
150	dnisson	1.15	edm::Service<TFileService> fs;
151			TH2D *histo = fs->make<TH2D>(oss.str().c_str(), oss.str().c_str(),
152	dnisson	1.14	30,
153	dnisson	1.13	highest_e_jet->eta()-0.5,
154			highest_e_jet->eta()+0.5,
155	dnisson	1.14	30,
156	dnisson	1.13	highest_e_jet->phi()-0.5,
157			highest_e_jet->phi()+0.5);
158
159	dnisson	1.17	if (smear_ > 0.0) {
160	dnisson	1.13	// create a smeared histo
161	dnisson	1.1	// create a temporary 2D vector for smeared energies
162	dnisson	1.13	double XbinSize = (histo->GetXaxis()->GetXmax()
163	dnisson	1.1	- histo->GetXaxis()->GetXmin()) /
164			static_cast<double>(histo->GetXaxis()->GetNbins());
165	dnisson	1.13	double YbinSize = (histo->GetYaxis()->GetXmax()
166	dnisson	1.1	- histo->GetYaxis()->GetXmin()) /
167			static_cast<double>(histo->GetYaxis()->GetNbins());
168	dnisson	1.13	double Xlo = histo->GetXaxis()->GetXmin();
169	dnisson	1.15	double Xhi = histo->GetXaxis()->GetXmax();
170	dnisson	1.13	double Ylo = histo->GetYaxis()->GetXmin();
171	dnisson	1.15	double Yhi = histo->GetYaxis()->GetXmax();
172	dnisson	1.14	vector< vector<double> > smeared(30, vector<double>(30, 0.0) );
173	dnisson	1.1	switch (smear_coord_) {
174			case 1:
175			for (int i = 0; i < particles.size(); i++) {
176	dnisson	1.4	double N = particles[i]->energy();
177	dnisson	1.2	double x = particles[i]->eta();
178			double y = particles[i]->phi();
179	dnisson	1.15	if (y < Ylo) {
180			y += 2.0*PI;
181			}
182			if (y > Yhi) {
183			y -= 2.0*PI;
184			}
185	dnisson	1.1	// loop over bins and add Gaussian in proper angle to smeared
186	dnisson	1.14	for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) {
187			for (vector< double >::size_type j2 = 0; j2 < 30; j2++) {
188	dnisson	1.1	double eta = static_cast<double>((signed int)i2) * XbinSize +
189	dnisson	1.15	Xlo - x;
190			double phi = static_cast<double>((signed int)j2) * YbinSize +
191			Ylo - y;
192
193	dnisson	1.1	// transform eta, phi to proper angle
194			double theta = 2.0*atan(exp(-eta));
195			double iota = asin(sin(theta)*sin(phi));
196
197			smeared[i2][j2] += (NXbinSizeYbinSize/(2.0PIsmear_*smear_))
198			* exp(-0.5(thetatheta + iotaiota)/(smear_smear_));
199			}
200			}
201			}
202			break;
203			case 0:
204			default:
205			for (int i = 0; i < particles.size(); i++) {
206	dnisson	1.4	double N = particles[i]->energy();
207	dnisson	1.2	double x = particles[i]->eta();
208			double y = particles[i]->phi();
209	dnisson	1.15	if (y < Ylo) {
210			y += 2.0*PI;
211			}
212			if (y > Yhi) {
213			y -= 2.0*PI;
214			}
215	dnisson	1.1	// loop over bins and add Gaussian to smeared
216	dnisson	1.14	for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) {
217			for (vector< double >::size_type j2 = 0; j2 < 30; j2++) {
218	dnisson	1.1	double eta = static_cast<double>((signed int)i2) * XbinSize
219	dnisson	1.15	+ Xlo - x;
220			double phi = static_cast<double>((signed int)j2) * YbinSize +
221			Ylo - y;
222
223	dnisson	1.1	smeared[i2][j2] += (NXbinSizeYbinSize/(2.0PIsmear_*smear_))
224			* exp(-0.5(etaeta + phiphi)/(smear_smear_));
225			}
226			}
227			}
228			}
229			// set histogram to match smear vector
230	dnisson	1.14	for (int i = 1; i <= 30; i++) {
231			for (int j = 1; j <= 30; j++) {
232	dnisson	1.1	histo->SetBinContent(i, j, smeared[i-1][j-1]);
233			}
234			}
235	dnisson	1.17	}
236			else {
237			// don't smear--just fill with particles
238			for (int i = 0; i < particles.size(); i++) {
239			histo->Fill(particles[i]->eta(), particles[i]->phi(),
240			particles[i]->energy());
241			}
242			}
243	dnisson	1.1
244			return histo;
245			}
246
247	dnisson	1.12	template <class Jet>
248			void seed_with_jetcoll(vector<Jet> jets,
249	dnisson	1.15	jetfit::model_def &_mdef, double phi_lo = -PI,
250			double phi_hi = PI) {
251	dnisson	1.11	// seed with jet collection
252	dnisson	1.1	int ijset = 0;
253			for (unsigned ij = 0; ij < jets.size(); ij++) {
254	dnisson	1.11	double N = jets[ij].energy();
255			double eta = jets[ij].eta();
256			double phi = jets[ij].phi();
257			if (N > 10.0) {
258	dnisson	1.19	_mdef.set_special_par(ijset, 0, N, _mdef.chisquare_error(N)*0.1,
259	dnisson	1.1	0.0, 1.0e6);
260	dnisson	1.11	_mdef.set_special_par(ijset, 1, eta, 0.01,
261	dnisson	1.1	0.0, 0.0);
262	dnisson	1.17	if (phi < phi_lo) {
263			phi += 2.0*PI;
264	dnisson	1.15	}
265	dnisson	1.17	if (phi > phi_hi) {
266			phi -= 2.0*PI;
267	dnisson	1.15	}
268	dnisson	1.17	_mdef.set_special_par(ijset, 2, phi, 0.01,
269	dnisson	1.1	0.0, 0.0);
270			_mdef.set_special_par(ijset, 3, 0.1, 0.001, 0.0, 0.0);
271			ijset++;
272			}
273			}
274			}
275
276			jetfit::model_def& JetFinderAnalyzer::make_model_def(const edm::Event& evt,
277			const edm::EventSetup&,
278			TH2 *histo) {
279			class jf_model_def : public jetfit::model_def {
280			public:
281			virtual double chisquare_error(double E) {
282	dnisson	1.19	return 0.298 - 0.9557*E; // study 11-09-09
283	dnisson	1.1	}
284			};
285
286			jf_model_def *_mdef = new jf_model_def();
287			TFormula *formula = new TFormula("gaus2d",
288			"[0]exp(-0.5((x-[1])2 + (y-[2])2)/([3]*2))/(2pi[3]*2)");
289			_mdef->set_formula(formula);
290			_mdef->set_indiv_max_E(0);
291			_mdef->set_indiv_max_x(1);
292			_mdef->set_indiv_max_y(2);
293			_mdef->set_indiv_par(0, string("N"), 0.0, 0.0, 0.0, 1.0e6);
294			_mdef->set_indiv_par(1, string("mu_x"), 0.0, 0.0, 0.0, 0.0);
295			_mdef->set_indiv_par(2, string("mu_y"), 0.0, 0.0, 0.0, 0.0);
296			_mdef->set_indiv_par(3, string("sig"), 0.1, 0.001, 0.0, 0.0);
297
298	dnisson	1.13	// get jetcoll from event file and select highest e jet
299	dnisson	1.12	if (info_type_ == 0) {
300			edm::Handle< vector<reco::GenJet> > jet_collection;
301			evt.getByLabel(jet_algo_, jet_collection);
302	dnisson	1.13	reco::GenJet highest_e_jet;
303			bool found_jet = false;
304			for (unsigned i = 0; i < jet_collection->size(); i++) {
305			if (!found_jet \|\| (*jet_collection)[i].energy() > highest_e_jet.energy())
306			{
307			highest_e_jet = (*jet_collection)[i];
308	dnisson	1.17	found_jet = true;
309	dnisson	1.13	}
310			}
311			vector<reco::GenJet> highest_e_jet_coll;
312			highest_e_jet_coll.push_back(highest_e_jet);
313	dnisson	1.15	seed_with_jetcoll(highest_e_jet_coll, *_mdef,
314			histo->GetYaxis()->GetXmin(),
315			histo->GetYaxis()->GetXmax());
316	dnisson	1.12	}
317			if (info_type_ == 1) {
318			edm::Handle< vector<reco::PFJet> > jet_collection;
319			evt.getByLabel(jet_algo_, jet_collection);
320	dnisson	1.13	reco::PFJet highest_e_jet;
321			bool found_jet = false;
322			for (unsigned i = 0; i < jet_collection->size(); i++) {
323			if (!found_jet \|\| (*jet_collection)[i].energy() > highest_e_jet.energy())
324			{
325			highest_e_jet = (*jet_collection)[i];
326	dnisson	1.17	found_jet = true;
327	dnisson	1.13	}
328			}
329			vector<reco::PFJet> highest_e_jet_coll;
330			highest_e_jet_coll.push_back(highest_e_jet);
331	dnisson	1.15	seed_with_jetcoll(highest_e_jet_coll, *_mdef,
332			histo->GetYaxis()->GetXmin(),
333			histo->GetYaxis()->GetXmax());
334	dnisson	1.12	}
335	dnisson	1.1
336			jetfit::set_model_def(_mdef);
337
338			// generate initial fit histogram
339			edm::Service<TFileService> fs;
340			TH2D *init_fit_histo = fs->make<TH2D>(("init_fit_"+string(histo->GetName()))
341			.c_str(),
342			("Initial fit for "
343			+string(histo->GetName())).c_str(),
344			histo->GetXaxis()->GetNbins(),
345			histo->GetXaxis()->GetXmin(),
346			histo->GetXaxis()->GetXmax(),
347			histo->GetXaxis()->GetNbins(),
348			histo->GetXaxis()->GetXmin(),
349			histo->GetXaxis()->GetXmax());
350			double XbinSize = (histo->GetXaxis()->GetXmax()
351			- histo->GetXaxis()->GetXmin()) /
352			static_cast<double>(histo->GetXaxis()->GetNbins());
353			double YbinSize = (histo->GetYaxis()->GetXmax()
354			- histo->GetYaxis()->GetXmin()) /
355			static_cast<double>(histo->GetYaxis()->GetNbins());
356			double Xlo = histo->GetXaxis()->GetXmin();
357			double Xhi = histo->GetXaxis()->GetXmax();
358			double Ylo = histo->GetYaxis()->GetXmin();
359			double Yhi = histo->GetYaxis()->GetXmax();
360
361	dnisson	1.19	for (int i = 0; i < 30; i++) {
362			for (int j = 0; j < 30; j++) {
363	dnisson	1.1	double x = (static_cast<double>(i) + 0.5)*XbinSize + Xlo;
364			double y = (static_cast<double>(j) + 0.5)*YbinSize + Ylo;
365			double pval[256];
366			if (_mdef->get_n_special_par_sets() > 64) {
367			cerr << "Parameter overload" << endl;
368			return *_mdef;
369			}
370			else {
371			for (int is = 0; is < _mdef->get_n_special_par_sets(); is++) {
372			for (int ii = 0; ii < 4; ii++) {
373			double spval, sperr, splo, sphi;
374			_mdef->get_special_par(is, ii, spval, sperr, splo, sphi);
375			pval[4*is + ii] = spval;
376			}
377			}
378			}
379			jetfit::set_ngauss(_mdef->get_n_special_par_sets());
380			init_fit_histo->SetBinContent(i+1, j+1,
381	dnisson	1.10	jetfit::fit_fcn(x, y, pval)
382			* XbinSize * YbinSize);
383	dnisson	1.1	}
384			}
385
386			return *_mdef;
387			}
388
389			void JetFinderAnalyzer::beginJob(const edm::EventSetup &es) {
390			ofs.open("jetfindlog.txt", ios::out);
391			if (ofs.fail()) {
392			cerr << "Opening jetfindlog.txt FAILED" << endl;
393			}
394			ofs << "Jetfinder log" << endl
395			<< "=============" << endl << endl;
396			}
397
398			ostream& operator<<(ostream &out, jetfit::trouble t) {
399			string action, error_string;
400
401			if (t.istat != 3) {
402			switch(t.occ) {
403			case jetfit::T_NULL:
404			action = "Program"; break;
405			case jetfit::T_SIMPLEX:
406			action = "SIMPLEX"; break;
407			case jetfit::T_MIGRAD:
408			action = "MIGRAD"; break;
409			case jetfit::T_MINOS:
410			action = "MINOS"; break;
411			default:
412			action = "Program"; break;
413			}
414
415			switch (t.istat) {
416			case 0:
417			error_string = "Unable to calculate error matrix"; break;
418			case 1:
419			error_string = "Error matrix a diagonal approximation"; break;
420			case 2:
421			error_string = "Error matrix not positive definite"; break;
422			case 3:
423			error_string = "Converged successfully"; break;
424			default:
425			ostringstream oss;
426			oss<<"Unknown status code "<<t.istat << flush;
427			error_string = oss.str(); break;
428			}
429
430			if (t.occ != jetfit::T_NULL)
431			out << action<<" trouble: "<<error_string;
432			else
433			out << "Not calculated" << endl;
434			}
435	dnisson	1.9	else {
436			out << "Error matrix accurate" << endl;
437			}
438	dnisson	1.1
439			return out;
440			}
441
442			void JetFinderAnalyzer::analyze_results(jetfit::results r,
443			std::vector<jetfit::trouble> t,
444			TH2 *hist_orig) {
445			ofs << "Histogram "<<hist_orig->GetName() << endl;
446	dnisson	1.9	for (int i = unique_jets[hist_orig].size() - r.chisquare.size();
447			i < unique_jets[hist_orig].size(); i++) {
448			int ir = i - unique_jets[hist_orig].size() + r.chisquare.size();
449	dnisson	1.1	ofs << "For "<<i+1<<" gaussians: " << endl
450	dnisson	1.9	<< t.at(i) << endl;
451			ofs << "chisquare="<<r.chisquare.at(ir)
452			<< endl;
453			ofs << unique_jets[hist_orig][i].size()<<" unique jets found" << endl;
454	dnisson	1.1	for (int j = 0; j < unique_jets[hist_orig][i].size(); j++) {
455			jet _jet = unique_jets[hist_orig][i][j];
456			ofs << "Jet "<<j<<": Energy = "<<_jet.energy<<", eta = "<<_jet.eta
457			<< ", phi = "<<_jet.phi << endl;
458			}
459			ofs << endl;
460			}
461			ofs << endl;
462
463			// save fit function histograms to root file
464			edm::Service<TFileService> fs;
465			for (vector< vector<double> >::size_type i = 0;
466			i < r.pval.size(); i++) {
467			jetfit::set_ngauss(r.pval[i].size() / 4);
468	dnisson	1.19	TF2 *tf2 = new TF2("fit_func", jetfit::fit_fcn_TF2,
469			hist_orig->GetXaxis()->GetXmin(),
470			hist_orig->GetXaxis()->GetXmax(),
471			hist_orig->GetYaxis()->GetXmin(),
472			hist_orig->GetYaxis()->GetXmax(),
473			r.pval[i].size());
474			for (vector<double>::size_type j = 0; j < r.pval[i].size(); j++) {
475			tf2->SetParameter(j, r.pval[i][j]);
476			}
477	dnisson	1.1	ostringstream fit_histo_oss;
478			fit_histo_oss << hist_orig->GetName()<<"_fit_"<<i << flush;
479	dnisson	1.19	tf2->SetNpx(hist_orig->GetXaxis()->GetNbins());
480			tf2->SetNpy(hist_orig->GetYaxis()->GetNbins());
481	dnisson	1.1	TH2D *fit_histo = fs->make<TH2D>(fit_histo_oss.str().c_str(),
482			fit_histo_oss.str().c_str(),
483			hist_orig->GetXaxis()->GetNbins(),
484			hist_orig->GetXaxis()->GetXmin(),
485			hist_orig->GetXaxis()->GetXmax(),
486			hist_orig->GetYaxis()->GetNbins(),
487			hist_orig->GetYaxis()->GetXmin(),
488			hist_orig->GetYaxis()->GetXmax());
489	dnisson	1.17	double XbinSize = (hist_orig->GetXaxis()->GetXmax()
490			- hist_orig->GetXaxis()->GetXmin())
491	dnisson	1.19	/ hist_orig->GetXaxis()->GetNbins();
492	dnisson	1.17	double YbinSize = (hist_orig->GetYaxis()->GetXmax()
493			- hist_orig->GetYaxis()->GetXmin())
494	dnisson	1.19	/ hist_orig->GetYaxis()->GetNbins();
495			for (int i = 1; i <= fit_histo->GetXaxis()->GetNbins(); i++) {
496			for (int j = 1; j <= fit_histo->GetYaxis()->GetNbins(); j++) {
497			double x = static_cast<double>(i - 1)
498			* (fit_histo->GetXaxis()->GetXmax()
499			- fit_histo->GetXaxis()->GetXmin())
500			/ static_cast<double>(fit_histo->GetXaxis()->GetNbins())
501			+ fit_histo->GetXaxis()->GetXmin();;
502			double y = static_cast<double>(j - 1)
503			* (fit_histo->GetYaxis()->GetXmax()
504			- fit_histo->GetYaxis()->GetXmin())
505			/ static_cast<double>(fit_histo->GetYaxis()->GetNbins())
506			+ fit_histo->GetYaxis()->GetXmin();
507			fit_histo->SetBinContent(i, j, tf2->Eval(x, y) * XbinSize * YbinSize);
508	dnisson	1.1	}
509			}
510			}
511
512			// save results to file
513			ostringstream res_tree_oss, rt_title_oss;
514			res_tree_oss << hist_orig->GetName()<<"_results" << flush;
515			rt_title_oss << "Fit results for "<<hist_orig->GetName() << flush;
516			}
517
518			DEFINE_FWK_MODULE(JetFinderAnalyzer);