JetFitAnalyzer/src/JetFitAnalyzer.cc

// -*- C++ -*-
//
// Package:    JetFitAnalyzer
// Class:      JetFitAnalyzer
// 
/**\class JetFitAnalyzer JetFitAnalyzer.cc UserCode/JetFitAnalyzer/src/JetFitAnalyzer.cc

 Description: Base class for using jetfit with CMSSW

 Implementation:
     o Method make_histo is used to prepare a histogram for analysis;
       the user can read this from a file or generate from the actual event
     o Method make_model_def is used to define a formula and initial
       parameters for the model used
     o Method analyze_results performs user analysis on results;
       the user is given a trouble vector and the original histogram in
       addition.
     o The function pointer user_minuit defines what to do after each fit
       is done. The pointer should be obtained in the subclass constructor.
*/
//
// Original Author: David Nisson  
//         Created:  Wed Aug  5 16:38:38 PDT 2009
// $Id: JetFitAnalyzer.cc,v 1.14.2.1 2010/01/10 01:11:25 dnisson Exp $
//
//

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

#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 "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "PhysicsTools/UtilAlgos/interface/TFileService.h"

#include "TFormula.h"
#define PI 3.141592

#include <iostream>
#include <sstream>
#include <string>

using namespace std;

JetFitAnalyzer::JetFitAnalyzer(const edm::ParameterSet& iConfig)
{
  // get parameters from iConfig
  P_cutoff_val_ = iConfig.getUntrackedParameter("P_cutoff_val", 0.5);
  infoType_ = iConfig.getUntrackedParameter("infoType", 1);
  smear_ = iConfig.getUntrackedParameter("smear", 0.05);
  jet_algo_ = iConfig.getParameter<string>("jet_algo");
}


JetFitAnalyzer::~JetFitAnalyzer()
{
 
   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.)

}


//
// member functions
//

// ------------ method called to for each event  ------------
void
JetFitAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
   using namespace edm;

   HistoFitter histoFitter;
   std::vector<HistoFitter::Trouble> t; // record troubles fitting

   TH2D *histo = make_histo(iEvent, iSetup); // create the histogram
   if (histo != 0) {
     HistoFitter::ModelDefinition &_mdef = make_model_def(iEvent, iSetup, histo);
     histoFitter.set_ModelDefinition(&_mdef);
     HistoFitter::FitResults r(histoFitter.fit_histo(histo, t, 1, P_cutoff_val_)); // fit the histogram
     analyze_results(r, t, histo); // perform user analysis
   }
   else {
     std::cerr << "Fitting not performed" << std::endl;
   }
}

/* create a histogram of the particles in the highest-E jet found */
TH2D * JetFitAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&) {
  const reco::Jet * highest_e_jet = 0;

  /* Find the highest-energy fat jet */
  if (infoType_ == 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 (infoType_ == 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;
  }

  /* Get the particles in the fat jet */
  vector<const reco::Candidate *> particles = 
    highest_e_jet->getJetConstituentsQuick();
 
  /* Create the histogram name */
  ostringstream oss;
  oss << "eta_phi_energy"<<evt.id().event() << flush;

  edm::Service<TFileService> fs;

  /* Initialize the histogram */
  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) {
    // get the limits of the histogram
    double Xlo = histo->GetXaxis()->GetXmin();
    double Xhi = histo->GetXaxis()->GetXmax();
    double Ylo = histo->GetYaxis()->GetXmin();
    double Yhi = histo->GetYaxis()->GetXmax();

    // get the bin sizes
    double XbinSize = (Xhi - Xlo) / static_cast<double>(histo->GetNbinsX());
    double YbinSize = (Yhi - Ylo) / static_cast<double>(histo->GetNbinsY());

    for (int i = 0; i < particles.size(); i++) {
      // compute the parameters of the Gaussian to be added
      double N = particles[i]->energy();
      double x = particles[i]->eta();
      double y = particles[i]->phi();

      // Compute correct phi for this histogram if it is out of range
      if (y < Ylo) {
        y += 2.0*PI;
      }
      if (y > Yhi) {
        y -= 2.0*PI;
      }

      // loop over bins and add Gaussian in eta and phi to histogram
      for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) {
        for (vector< double >::size_type j2 = 0; j2 < 30; j2++) {
          // compute the distance of the bin from the particle eta, phi
          double eta = static_cast<double>((signed int)i2) * XbinSize + 
            Xlo - x;
          double phi = static_cast<double>((signed int)j2) * YbinSize +
            Ylo - y;
          
          /* Fill the histogram with a Gaussian centered around particle */
          histo->SetBinContent(i2+1, j2+1, histo->GetBinContent(i2+1, j2+1) 
                               + (N*XbinSize*YbinSize/(2.0*PI*smear_*smear_))
            * exp(-0.5*(eta*eta + phi*phi)/(smear_*smear_)));
        }
      }
    }
  }
  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;
}

/* Create a model definition based on the histogram and event information */
HistoFitter::ModelDefinition& JetFitAnalyzer::make_model_def(const edm::Event& evt,
                                                 const edm::EventSetup&,
                                                 TH2 *histo) {
  class DefaultModelDef : public HistoFitter::ModelDefinition {
  public:
    virtual double chisquareSigma(double E) {
      return 0.6*E + 1.0; // based on latest info
    }
  };

  DefaultModelDef *_mdef = new DefaultModelDef();

  // 2D Gaussian formula
  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);

  // set fixed Gaussian sigma equal to smearing
  _mdef->setIndivParameter(3, string("sig"), smear_, 0.001, 0.0, 0.0);

  _mdef->setMaxNumberOfGaussians(6);

  return *_mdef;
}

// ------------ method called once each job just before starting event loop  ------------
void 
JetFitAnalyzer::beginJob(const edm::EventSetup&)
{
}

// ------------ method called once each job just after ending the event loop  ------------
void 
JetFitAnalyzer::endJob() {
}

//define this as a plug-in
// this is a base class-we can't do this DEFINE_FWK_MODULE(JetFitAnalyzer);
Revision:	1.15
Committed:	Sun Jan 10 01:36:29 2010 UTC (15 years, 4 months ago) by dnisson
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V02-00-00
Changes since 1.14:	+161 -6 lines
Log Message:	Merging Spectrum branch back to trunk
#	User	Rev	Content
1	dnisson	1.1	// -- C++ --
2			//
3			// Package: JetFitAnalyzer
4			// Class: JetFitAnalyzer
5			//
6			/**\class JetFitAnalyzer JetFitAnalyzer.cc UserCode/JetFitAnalyzer/src/JetFitAnalyzer.cc
7
8			Description: Base class for using jetfit with CMSSW
9
10			Implementation:
11			o Method make_histo is used to prepare a histogram for analysis;
12			the user can read this from a file or generate from the actual event
13			o Method make_model_def is used to define a formula and initial
14			parameters for the model used
15			o Method analyze_results performs user analysis on results;
16			the user is given a trouble vector and the original histogram in
17			addition.
18			o The function pointer user_minuit defines what to do after each fit
19			is done. The pointer should be obtained in the subclass constructor.
20			*/
21			//
22			// Original Author: David Nisson
23			// Created: Wed Aug 5 16:38:38 PDT 2009
24	dnisson	1.15	// $Id: JetFitAnalyzer.cc,v 1.14.2.1 2010/01/10 01:11:25 dnisson Exp $
25	dnisson	1.1	//
26			//
27
28			#include "UserCode/JetFitAnalyzer/interface/JetFitAnalyzer.h"
29
30	dnisson	1.15	#include "DataFormats/JetReco/interface/PFJetCollection.h"
31			#include "DataFormats/JetReco/interface/GenJetCollection.h"
32			#include "DataFormats/JetReco/interface/PFJet.h"
33			#include "DataFormats/JetReco/interface/GenJet.h"
34
35			#include "FWCore/ServiceRegistry/interface/Service.h"
36			#include "FWCore/MessageLogger/interface/MessageLogger.h"
37			#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
38
39			#include "TFormula.h"
40			#define PI 3.141592
41
42			#include <iostream>
43			#include <sstream>
44			#include <string>
45
46			using namespace std;
47
48	dnisson	1.1	JetFitAnalyzer::JetFitAnalyzer(const edm::ParameterSet& iConfig)
49			{
50			// get parameters from iConfig
51			P_cutoff_val_ = iConfig.getUntrackedParameter("P_cutoff_val", 0.5);
52	dnisson	1.15	infoType_ = iConfig.getUntrackedParameter("infoType", 1);
53			smear_ = iConfig.getUntrackedParameter("smear", 0.05);
54			jet_algo_ = iConfig.getParameter<string>("jet_algo");
55	dnisson	1.1	}
56
57
58			JetFitAnalyzer::~JetFitAnalyzer()
59			{
60
61			// do anything here that needs to be done at desctruction time
62			// (e.g. close files, deallocate resources etc.)
63
64			}
65
66
67			//
68			// member functions
69			//
70
71			// ------------ method called to for each event ------------
72			void
73			JetFitAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
74			{
75			using namespace edm;
76
77	dnisson	1.13	HistoFitter histoFitter;
78	dnisson	1.15	std::vector<HistoFitter::Trouble> t; // record troubles fitting
79	dnisson	1.1
80	dnisson	1.15	TH2D *histo = make_histo(iEvent, iSetup); // create the histogram
81	dnisson	1.2	if (histo != 0) {
82	dnisson	1.11	HistoFitter::ModelDefinition &_mdef = make_model_def(iEvent, iSetup, histo);
83	dnisson	1.13	histoFitter.set_ModelDefinition(&_mdef);
84	dnisson	1.15	HistoFitter::FitResults r(histoFitter.fit_histo(histo, t, 1, P_cutoff_val_)); // fit the histogram
85			analyze_results(r, t, histo); // perform user analysis
86	dnisson	1.2	}
87			else {
88	dnisson	1.3	std::cerr << "Fitting not performed" << std::endl;
89	dnisson	1.2	}
90	dnisson	1.1	}
91
92	dnisson	1.15	/* create a histogram of the particles in the highest-E jet found */
93			TH2D * JetFitAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&) {
94			const reco::Jet * highest_e_jet = 0;
95
96			/* Find the highest-energy fat jet */
97			if (infoType_ == 0) {
98			edm::Handle< vector<reco::GenJet> > evtJets;
99			evt.getByLabel(jet_algo_, evtJets);
100			for (unsigned i = 0; i < evtJets->size(); i++) {
101			if (highest_e_jet == 0
102			\|\| (*evtJets)[i].energy() > highest_e_jet->energy()) {
103			highest_e_jet = &((*evtJets)[i]);
104			}
105			}
106			}
107			else if (infoType_ == 1) {
108			edm::Handle< vector<reco::PFJet> > evtJets;
109			evt.getByLabel(jet_algo_, evtJets);
110			for (unsigned i = 0; i < evtJets->size(); i++) {
111			if (highest_e_jet == 0
112			\|\| (*evtJets)[i].energy() > highest_e_jet->energy()) {
113			highest_e_jet = &((*evtJets)[i]);
114			}
115			}
116			}
117
118			if (highest_e_jet == 0) {
119			cerr << "No fat jets found!" << endl;
120			return 0;
121			}
122
123			/* Get the particles in the fat jet */
124			vector<const reco::Candidate *> particles =
125			highest_e_jet->getJetConstituentsQuick();
126
127			/* Create the histogram name */
128			ostringstream oss;
129			oss << "eta_phi_energy"<<evt.id().event() << flush;
130
131			edm::Service<TFileService> fs;
132
133			/* Initialize the histogram */
134			TH2D *histo = fs->make<TH2D>(oss.str().c_str(), oss.str().c_str(),
135			30,
136			highest_e_jet->eta()-0.5,
137			highest_e_jet->eta()+0.5,
138			30,
139			highest_e_jet->phi()-0.5,
140			highest_e_jet->phi()+0.5);
141
142			if (smear_ > 0.0) {
143			// get the limits of the histogram
144			double Xlo = histo->GetXaxis()->GetXmin();
145			double Xhi = histo->GetXaxis()->GetXmax();
146			double Ylo = histo->GetYaxis()->GetXmin();
147			double Yhi = histo->GetYaxis()->GetXmax();
148
149			// get the bin sizes
150			double XbinSize = (Xhi - Xlo) / static_cast<double>(histo->GetNbinsX());
151			double YbinSize = (Yhi - Ylo) / static_cast<double>(histo->GetNbinsY());
152
153			for (int i = 0; i < particles.size(); i++) {
154			// compute the parameters of the Gaussian to be added
155			double N = particles[i]->energy();
156			double x = particles[i]->eta();
157			double y = particles[i]->phi();
158
159			// Compute correct phi for this histogram if it is out of range
160			if (y < Ylo) {
161			y += 2.0*PI;
162			}
163			if (y > Yhi) {
164			y -= 2.0*PI;
165			}
166
167			// loop over bins and add Gaussian in eta and phi to histogram
168			for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) {
169			for (vector< double >::size_type j2 = 0; j2 < 30; j2++) {
170			// compute the distance of the bin from the particle eta, phi
171			double eta = static_cast<double>((signed int)i2) * XbinSize +
172			Xlo - x;
173			double phi = static_cast<double>((signed int)j2) * YbinSize +
174			Ylo - y;
175
176			/* Fill the histogram with a Gaussian centered around particle */
177			histo->SetBinContent(i2+1, j2+1, histo->GetBinContent(i2+1, j2+1)
178			+ (NXbinSizeYbinSize/(2.0PIsmear_*smear_))
179			* exp(-0.5(etaeta + phiphi)/(smear_smear_)));
180			}
181			}
182			}
183			}
184			else {
185			// don't smear--just fill with particles
186			for (int i = 0; i < particles.size(); i++) {
187			histo->Fill(particles[i]->eta(), particles[i]->phi(),
188			particles[i]->energy());
189			}
190			}
191
192			return histo;
193			}
194
195			/* Create a model definition based on the histogram and event information */
196			HistoFitter::ModelDefinition& JetFitAnalyzer::make_model_def(const edm::Event& evt,
197			const edm::EventSetup&,
198			TH2 *histo) {
199			class DefaultModelDef : public HistoFitter::ModelDefinition {
200			public:
201			virtual double chisquareSigma(double E) {
202			return 0.6*E + 1.0; // based on latest info
203			}
204			};
205
206			DefaultModelDef *_mdef = new DefaultModelDef();
207
208			// 2D Gaussian formula
209			TFormula *formula = new TFormula("gaus2d",
210			"[0]exp(-0.5((x-[1])2 + (y-[2])2)/([3]*2))/(2pi[3]*2)");
211			_mdef->setFormula(formula);
212
213			_mdef->setIndivEnergy(0);
214			_mdef->setIndivMeanX(1);
215			_mdef->setIndivMeanY(2);
216			_mdef->setIndivParameter(0, string("N"), 0.0, 0.0, 0.0, 1.0e6);
217			_mdef->setIndivParameter(1, string("mu_x"), 0.0, 0.0, 0.0, 0.0);
218			_mdef->setIndivParameter(2, string("mu_y"), 0.0, 0.0, 0.0, 0.0);
219
220			// set fixed Gaussian sigma equal to smearing
221			_mdef->setIndivParameter(3, string("sig"), smear_, 0.001, 0.0, 0.0);
222
223			_mdef->setMaxNumberOfGaussians(6);
224
225			return *_mdef;
226			}
227	dnisson	1.1
228			// ------------ method called once each job just before starting event loop ------------
229			void
230			JetFitAnalyzer::beginJob(const edm::EventSetup&)
231			{
232			}
233
234			// ------------ method called once each job just after ending the event loop ------------
235			void
236			JetFitAnalyzer::endJob() {
237			}
238
239			//define this as a plug-in
240			// this is a base class-we can't do this DEFINE_FWK_MODULE(JetFitAnalyzer);