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.16 2010/01/22 23:53:09 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");
}


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<TH2D *> histos = make_histo(iEvent, iSetup); // 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_)); // fit the histogram
       analyze_results(r, t, histo); // 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<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());

  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;

    /* save the true W decay products */
    /* Get the generator-level info */
    edm::Handle<edm::HepMCProduct> genProduct;
    evt.getByLabel("generator", genProduct);
    const HepMC::GenEvent *genEvt = genProduct->GetEvent();
    
    /* Save to file */
    ofstream ofs("true_W_decay.txt", ios_base::app);
    ofs << "Event "<<evt.id().event()<<", jet"<<i<<":" << endl;
    ofs << "PDG ID"<<'\t'<<"Eta"<<'\t'<<"Phi"<<'\t'<<"Energy"<<endl;
    for (HepMC::GenEvent::particle_const_iterator ppit = genEvt->particles_begin();
         ppit != genEvt->particles_end(); ppit++) {
      if ((*ppit)->pdg_id() == 24 || (*ppit)->pdg_id() == -24) {
        // found W-iterate over daughters
        HepMC::GenVertex *vtx = (*ppit)->end_vertex();
        if (vtx != 0)
          for (HepMC::GenVertex::particles_out_const_iterator pit = vtx->particles_out_const_begin();
               pit != vtx->particles_out_const_end(); pit++) {
            ofs << (*pit)->pdg_id() << '\t' << (*pit)->momentum().eta() << '\t'
                << (*pit)->momentum().phi() << '\t' << (*pit)->momentum().e() << endl;
          }
      }
    }
    ofs.close();
    
    /* 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);
    
    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) {
      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(4);

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