JetFitAnalyzer/src/JetFitAnalyzer.cc

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

 Defaription: 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.22 2010/03/27 01:58:10 dnisson Exp $
//
//

#include <HepMC/GenEvent.h>
#include <HepMC/GenVertex.h>
#include <HepMC/GenParticle.h>

#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 "SimDataFormats/GeneratorProducts/interface/HepMCProduct.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>
#include <fstream>
// to save the W decay product energy, mass, eta, phi

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");
  maxGaussians_ = iConfig.getUntrackedParameter("maxGaussians", 4);
  gausSpread_ = iConfig.getUntrackedParameter("gausSpread", 0.1);
}


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

   std::vector<unsigned> numberOfParticles;
   std::vector<TH2D *> histos = make_histo(iEvent, iSetup, numberOfParticles);
   // create the histogrms
   for (unsigned i= 0; i < histos.size(); i++) {
     TH2D *histo = histos[i]; // obtain 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_,
                                                       numberOfParticles[i])); 
       // fit the histogram
       analyze_results(r, t, histo, iEvent); // perform user analysis
     }
     else {
       std::cerr << "Fitting not performed" << std::endl;
     }
   }
}

/* create histograms of the particles in all jets found */
std::vector<TH2D *> JetFitAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&, std::vector<unsigned> &numberOfParticles) {
  std::vector<TH2D*> histos;

  std::vector<const reco::Jet *> jets;
  if (infoType_ == 0) {
    edm::Handle< vector<reco::GenJet> > evtJets;
    evt.getByLabel(jet_algo_, evtJets);
    jets.resize(evtJets->size());
    for (unsigned i = 0; i < evtJets->size(); i++) {
      jets[i] = &((*evtJets)[i]);
    }
  }
  else if (infoType_ == 1) {
    edm::Handle< vector<reco::PFJet> > evtJets;
    evt.getByLabel(jet_algo_, evtJets);
    jets.resize(evtJets->size());
    for (unsigned i = 0; i < evtJets->size(); i++) {
      jets[i] = &((*evtJets)[i]);
    }
  }
    
  histos.resize(jets.size());
  numberOfParticles.resize(jets.size());

  for (unsigned i = 0; i < jets.size(); i++) {
    /* Get the particles in the fat jet */
    vector<const reco::Candidate *> particles = 
      jets[i]->getJetConstituentsQuick();
 
    /* Create the histogram name */
    ostringstream oss;
    oss << "eta_phi_energy"<<evt.id().event() <<"jet"<<i << flush;
    
    edm::Service<TFileService> fs;

    /* Initialize the histogram */
    TH2D *histo = fs->make<TH2D>(oss.str().c_str(), oss.str().c_str(),
                                 30,
                                 jets[i]->eta()-0.5,
                                 jets[i]->eta()+0.5,
                                 30,
                                 jets[i]->phi()-0.5,
                                 jets[i]->phi()+0.5);
    numberOfParticles[i] = particles.size(); 
    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());
      }
    }
    histos[i] = histo;
  }
 
  return histos;
}

/* 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, double Emax) {
      double sigma = 0.14*Emax*(1.0 - exp(-2.3*E/Emax));

      return sigma;
      // based on Nisson Elog entry 33
    }
  };

  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 gausSpread
  _mdef->setIndivParameter(3, string("sig"), gausSpread_, 0.001, 0.0, 0.0);

  _mdef->setMaxNumberOfGaussians(maxGaussians_);

  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.23
Committed:	Thu Apr 8 05:01:28 2010 UTC (15 years ago) by dnisson
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V02-04-01, V02-04-00, HEAD
Changes since 1.22:	+2 -27 lines
Log Message:	Fixed bug in HistoFitter::fcn
#	Content
1	//- C++ --
2	//
3	// Package: JetFitAnalyzer
4	// Class: JetFitAnalyzer
5	//
6	/**\class JetFitAnalyzer JetFitAnalyzer.cc UserCode/JetFitAnalyzer/src/JetFitAnalyzer.cc
7
8	Defaription: 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	// $Id: JetFitAnalyzer.cc,v 1.22 2010/03/27 01:58:10 dnisson Exp $
25	//
26	//
27
28	#include <HepMC/GenEvent.h>
29	#include <HepMC/GenVertex.h>
30	#include <HepMC/GenParticle.h>
31
32	#include "UserCode/JetFitAnalyzer/interface/JetFitAnalyzer.h"
33
34	#include "DataFormats/JetReco/interface/PFJetCollection.h"
35	#include "DataFormats/JetReco/interface/GenJetCollection.h"
36	#include "DataFormats/JetReco/interface/PFJet.h"
37	#include "DataFormats/JetReco/interface/GenJet.h"
38	#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
39
40	#include "FWCore/ServiceRegistry/interface/Service.h"
41	#include "FWCore/MessageLogger/interface/MessageLogger.h"
42	#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
43
44	#include "TFormula.h"
45	#define PI 3.141592
46
47	#include <iostream>
48	#include <sstream>
49	#include <string>
50	#include <fstream>
51	// to save the W decay product energy, mass, eta, phi
52
53	using namespace std;
54
55	JetFitAnalyzer::JetFitAnalyzer(const edm::ParameterSet& iConfig)
56	{
57	// get parameters from iConfig
58	P_cutoff_val_ = iConfig.getUntrackedParameter("P_cutoff_val", 0.5);
59	infoType_ = iConfig.getUntrackedParameter("infoType", 1);
60	smear_ = iConfig.getUntrackedParameter("smear", 0.05);
61	jet_algo_ = iConfig.getParameter<string>("jet_algo");
62	maxGaussians_ = iConfig.getUntrackedParameter("maxGaussians", 4);
63	gausSpread_ = iConfig.getUntrackedParameter("gausSpread", 0.1);
64	}
65
66
67	JetFitAnalyzer::~JetFitAnalyzer()
68	{
69
70	// do anything here that needs to be done at desctruction time
71	// (e.g. close files, deallocate resources etc.)
72
73	}
74
75
76	//
77	// member functions
78	//
79
80	// ------------ method called to for each event ------------
81	void
82	JetFitAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
83	{
84	using namespace edm;
85
86	HistoFitter histoFitter;
87	std::vector<HistoFitter::Trouble> t; // record troubles fitting
88
89	std::vector<unsigned> numberOfParticles;
90	std::vector<TH2D *> histos = make_histo(iEvent, iSetup, numberOfParticles);
91	// create the histogrms
92	for (unsigned i= 0; i < histos.size(); i++) {
93	TH2D *histo = histos[i]; // obtain the histogram
94	if (histo != 0) {
95	HistoFitter::ModelDefinition &_mdef = make_model_def(iEvent, iSetup, histo);
96	histoFitter.set_ModelDefinition(&_mdef);
97	HistoFitter::FitResults r(histoFitter.fit_histo(histo, t, 1,
98	P_cutoff_val_,
99	numberOfParticles[i]));
100	// fit the histogram
101	analyze_results(r, t, histo, iEvent); // perform user analysis
102	}
103	else {
104	std::cerr << "Fitting not performed" << std::endl;
105	}
106	}
107	}
108
109	/* create histograms of the particles in all jets found */
110	std::vector<TH2D *> JetFitAnalyzer::make_histo(const edm::Event &evt, const edm::EventSetup&, std::vector<unsigned> &numberOfParticles) {
111	std::vector<TH2D*> histos;
112
113	std::vector<const reco::Jet *> jets;
114	if (infoType_ == 0) {
115	edm::Handle< vector<reco::GenJet> > evtJets;
116	evt.getByLabel(jet_algo_, evtJets);
117	jets.resize(evtJets->size());
118	for (unsigned i = 0; i < evtJets->size(); i++) {
119	jets[i] = &((*evtJets)[i]);
120	}
121	}
122	else if (infoType_ == 1) {
123	edm::Handle< vector<reco::PFJet> > evtJets;
124	evt.getByLabel(jet_algo_, evtJets);
125	jets.resize(evtJets->size());
126	for (unsigned i = 0; i < evtJets->size(); i++) {
127	jets[i] = &((*evtJets)[i]);
128	}
129	}
130
131	histos.resize(jets.size());
132	numberOfParticles.resize(jets.size());
133
134	for (unsigned i = 0; i < jets.size(); i++) {
135	/* Get the particles in the fat jet */
136	vector<const reco::Candidate *> particles =
137	jets[i]->getJetConstituentsQuick();
138
139	/* Create the histogram name */
140	ostringstream oss;
141	oss << "eta_phi_energy"<<evt.id().event() <<"jet"<<i << flush;
142
143	edm::Service<TFileService> fs;
144
145	/* Initialize the histogram */
146	TH2D *histo = fs->make<TH2D>(oss.str().c_str(), oss.str().c_str(),
147	30,
148	jets[i]->eta()-0.5,
149	jets[i]->eta()+0.5,
150	30,
151	jets[i]->phi()-0.5,
152	jets[i]->phi()+0.5);
153	numberOfParticles[i] = particles.size();
154	if (smear_ > 0.0) {
155	// get the limits of the histogram
156	double Xlo = histo->GetXaxis()->GetXmin();
157	double Xhi = histo->GetXaxis()->GetXmax();
158	double Ylo = histo->GetYaxis()->GetXmin();
159	double Yhi = histo->GetYaxis()->GetXmax();
160
161	// get the bin sizes
162	double XbinSize = (Xhi - Xlo) / static_cast<double>(histo->GetNbinsX());
163	double YbinSize = (Yhi - Ylo) / static_cast<double>(histo->GetNbinsY());
164
165	for (int i = 0; i < particles.size(); i++) {
166	// compute the parameters of the Gaussian to be added
167	double N = particles[i]->energy();
168	double x = particles[i]->eta();
169	double y = particles[i]->phi();
170
171	// Compute correct phi for this histogram if it is out of range
172	if (y < Ylo) {
173	y += 2.0*PI;
174	}
175	if (y > Yhi) {
176	y -= 2.0*PI;
177	}
178
179	// loop over bins and add Gaussian in eta and phi to histogram
180	for (vector< vector<double> >::size_type i2 = 0; i2 < 30; i2++) {
181	for (vector< double >::size_type j2 = 0; j2 < 30; j2++) {
182	// compute the distance of the bin from the particle eta, phi
183	double eta = static_cast<double>((signed int)i2) * XbinSize +
184	Xlo - x;
185	double phi = static_cast<double>((signed int)j2) * YbinSize +
186	Ylo - y;
187
188	/* Fill the histogram with a Gaussian centered around particle */
189	histo->SetBinContent(i2+1, j2+1, histo->GetBinContent(i2+1, j2+1)
190	+ (NXbinSizeYbinSize/(2.0PIsmear_*smear_))
191	* exp(-0.5(etaeta + phiphi)/(smear_smear_)));
192	}
193	}
194	}
195	}
196	else {
197	// don't smear--just fill with particles
198	for (int i = 0; i < particles.size(); i++) {
199	histo->Fill(particles[i]->eta(), particles[i]->phi(),
200	particles[i]->energy());
201	}
202	}
203	histos[i] = histo;
204	}
205
206	return histos;
207	}
208
209	/* Create a model definition based on the histogram and event information */
210	HistoFitter::ModelDefinition& JetFitAnalyzer::make_model_def(const edm::Event& evt,
211	const edm::EventSetup&,
212	TH2 *histo) {
213	class DefaultModelDef : public HistoFitter::ModelDefinition {
214	public:
215	virtual double chisquareSigma(double E, double Emax) {
216	double sigma = 0.14Emax(1.0 - exp(-2.3*E/Emax));
217
218	return sigma;
219	// based on Nisson Elog entry 33
220	}
221	};
222
223	DefaultModelDef *_mdef = new DefaultModelDef();
224
225	// 2D Gaussian formula
226	TFormula *formula = new TFormula("gaus2d",
227	"[0]exp(-0.5((x-[1])2 + (y-[2])2)/([3]*2))/(2pi[3]*2)");
228	_mdef->setFormula(formula);
229
230	_mdef->setIndivEnergy(0);
231	_mdef->setIndivMeanX(1);
232	_mdef->setIndivMeanY(2);
233	_mdef->setIndivParameter(0, string("N"), 0.0, 0.0, 0.0, 1.0e6);
234	_mdef->setIndivParameter(1, string("mu_x"), 0.0, 0.0, 0.0, 0.0);
235	_mdef->setIndivParameter(2, string("mu_y"), 0.0, 0.0, 0.0, 0.0);
236
237	// set fixed Gaussian sigma equal to gausSpread
238	_mdef->setIndivParameter(3, string("sig"), gausSpread_, 0.001, 0.0, 0.0);
239
240	_mdef->setMaxNumberOfGaussians(maxGaussians_);
241
242	return *_mdef;
243	}
244
245	// ------------ method called once each job just before starting event loop ------------
246	void
247	JetFitAnalyzer::beginJob(const edm::EventSetup&)
248	{
249	}
250
251	// ------------ method called once each job just after ending the event loop ------------
252	void
253	JetFitAnalyzer::endJob() {
254	}
255
256	//define this as a plug-in
257	// this is a base class-we can't do this DEFINE_FWK_MODULE(JetFitAnalyzer);