|
|
[Rivet-svn] r1769 - in trunk: include/Rivet include/Rivet/Analyses src/Analyses
blackhole at projects.hepforge.org
blackhole at projects.hepforge.org
Wed Aug 19 17:52:55 BST 2009
Author: fsiegert
Date: Wed Aug 19 17:52:54 2009
New Revision: 1769
Log:
Introduce an analysis class MC_JetAnalysis which MC analyses can derive
from to look at jet observables, like differential/integrated jet
resolutions, pT, eta, dR, deta, multiplicities.
This is useful to avoid duplications in MC validation analyses looking
at e.g. Z+jets, W+jets, photon+jets, pure jets, ... (some of which I
will start adding/transforming in the next step).
Added:
trunk/include/Rivet/Analyses/MC_JetAnalysis.hh
trunk/src/Analyses/MC_JetAnalysis.cc
Modified:
trunk/include/Rivet/Makefile.am
trunk/src/Analyses/Makefile.am
Added: trunk/include/Rivet/Analyses/MC_JetAnalysis.hh
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/include/Rivet/Analyses/MC_JetAnalysis.hh Wed Aug 19 17:52:54 2009 (r1769)
@@ -0,0 +1,55 @@
+// -*- C++ -*-
+#ifndef RIVET_MC_JetAnalysis_HH
+#define RIVET_MC_JetAnalysis_HH
+
+#include "Rivet/Analysis.hh"
+#include "Rivet/Projections/FinalState.hh"
+
+namespace Rivet {
+
+
+ class MC_JetAnalysis : public Analysis {
+
+ public:
+
+ /// Default constructor.
+ MC_JetAnalysis(const std::string& name, const double& sqrts,
+ const size_t& njet, const std::string& jetpro_name);
+
+
+ /// @name Analysis methods
+ //@{
+ virtual void init();
+ virtual void analyze(const Event& event);
+ virtual void finalize();
+ //@}
+
+ protected:
+
+ /// The energy scale and number of jets for which histograms are to be
+ /// initialised
+ double m_sqrts;
+ size_t m_njet;
+
+ /// The name of the jet projection to be used for this analysis
+ /// (this projection has to be registered by the derived analysis!)
+ const std::string m_jetpro_name;
+
+ /// @name Histograms
+ //@{
+ std::vector<AIDA::IHistogram1D *> _h_log10_d;
+ std::vector<AIDA::IDataPointSet *> _h_log10_R;
+ std::vector<AIDA::IHistogram1D *> _h_pT_jet;
+ std::vector<AIDA::IHistogram1D *> _h_eta_jet;
+ std::map<std::pair<size_t, size_t>, AIDA::IHistogram1D*> _h_deta_jets;
+ std::map<std::pair<size_t, size_t>, AIDA::IHistogram1D*> _h_dR_jets;
+ AIDA::IHistogram1D * _h_jet_multi_exclusive;
+ AIDA::IHistogram1D * _h_jet_multi_inclusive;
+ AIDA::IDataPointSet * _h_jet_multi_ratio;
+ //@}
+
+ };
+
+}
+
+#endif
Modified: trunk/include/Rivet/Makefile.am
==============================================================================
--- trunk/include/Rivet/Makefile.am Wed Aug 19 17:46:01 2009 (r1768)
+++ trunk/include/Rivet/Makefile.am Wed Aug 19 17:52:54 2009 (r1769)
@@ -85,6 +85,7 @@
Analyses/PDG_Hadron_Multiplicities.hh \
Analyses/PDG_Hadron_Multiplicities_Ratios.hh \
Analyses/MC_TVT1960_ZJETS.hh \
+ Analyses/MC_JetAnalysis.hh \
Analyses/MC_LHC_LEADINGJETS.hh \
Analyses/MC_LHC_ZANALYSIS.hh \
Analyses/MC_LHC_DIJET.hh \
Added: trunk/src/Analyses/MC_JetAnalysis.cc
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/src/Analyses/MC_JetAnalysis.cc Wed Aug 19 17:52:54 2009 (r1769)
@@ -0,0 +1,172 @@
+// -*- C++ -*-
+#include "Rivet/Analyses/MC_JetAnalysis.hh"
+#include "Rivet/Tools/Logging.hh"
+#include "Rivet/Projections/FastJets.hh"
+#include "Rivet/RivetAIDA.hh"
+
+namespace Rivet {
+
+
+ MC_JetAnalysis::MC_JetAnalysis(const std::string& name, const double& sqrts,
+ const size_t& njet, const std::string& jetpro_name)
+ : Analysis(name), m_sqrts(sqrts), m_njet(njet), m_jetpro_name(jetpro_name),
+ _h_log10_d(njet, NULL), _h_log10_R(njet+1, NULL), _h_pT_jet(njet, NULL),
+ _h_eta_jet(njet, NULL)
+ {
+ setNeedsCrossSection(true);
+ }
+
+
+
+ // Book histograms
+ void MC_JetAnalysis::init() {
+
+ for (size_t i=0; i<m_njet; ++i) {
+ stringstream dname;
+ dname<<"log10_d_"<<i<<i+1;
+ _h_log10_d[i] = bookHistogram1D(dname.str(), 50, 0.2, 2.6);
+
+ stringstream Rname;
+ Rname<<"log10_R_"<<i;
+ _h_log10_R[i] = bookDataPointSet(Rname.str(), 50, 0.2, 2.6);
+
+ stringstream pTname;
+ pTname<<"pT_jet_"<<i;
+ double pTmax = 1.0/(double(i)+2.0)*m_sqrts/2.0;
+ int nbins = 100/(i+1);
+ _h_pT_jet[i] = bookHistogram1D(pTname.str(), nbins, 0.0, pTmax);
+
+ stringstream etaname;
+ etaname<<"eta_jet_"<<i;
+ _h_eta_jet[i] = bookHistogram1D(etaname.str(), 50, -5.0, 5.0);
+
+ for (size_t j=i+1; j<m_njet; ++j) {
+ std::pair<size_t, size_t> ij(std::make_pair(i, j));
+
+ stringstream detaname;
+ detaname<<"deta_jets_"<<i<<j;
+ _h_deta_jets.insert(make_pair(ij, bookHistogram1D(detaname.str(), 50, -5.0, 5.0)));
+
+ stringstream dRname;
+ dRname<<"dR_jets_"<<i<<j;
+ _h_dR_jets.insert(make_pair(ij, bookHistogram1D(dRname.str(), 25, 0.0, 5.0)));
+ }
+ }
+ stringstream Rname;
+ Rname<<"log10_R_"<<m_njet;
+ _h_log10_R[m_njet] = bookDataPointSet(Rname.str(), 50, 0.2, 2.6);
+
+ _h_jet_multi_exclusive = bookHistogram1D("jet_multi_exclusive", m_njet+3, -0.5, m_njet+3-0.5);
+ _h_jet_multi_inclusive = bookHistogram1D("jet_multi_inclusive", m_njet+3, -0.5, m_njet+3-0.5);
+ _h_jet_multi_ratio = bookDataPointSet("jet_multi_ratio", m_njet+2, 0.5, m_njet+3-0.5);
+ }
+
+
+
+ // Do the analysis
+ void MC_JetAnalysis::analyze(const Event & e) {
+ double weight = e.weight();
+
+ const FastJets& jetpro = applyProjection<FastJets>(e, m_jetpro_name);
+
+ // jet resolutions and integrated jet rates
+ const fastjet::ClusterSequence* seq = jetpro.clusterSeq();
+ if (seq!=NULL) {
+ double previous_dij = 10.0;
+ for (size_t i=0; i<m_njet; ++i) {
+ // jet resolution i -> j
+ double d_ij=log10(sqrt(seq->exclusive_dmerge_max(i)));
+
+ // fill differential jet resolution
+ _h_log10_d[i]->fill(d_ij, weight);
+
+ // fill integrated jet resolution
+ for (int ibin=0; ibin<_h_log10_R[i]->size(); ++ibin) {
+ IDataPoint* dp=_h_log10_R[i]->point(ibin);
+ double dcut=dp->coordinate(0)->value();
+ if (d_ij<dcut && previous_dij>dcut) {
+ dp->coordinate(1)->setValue(dp->coordinate(1)->value()+weight);
+ }
+ }
+ previous_dij = d_ij;
+ }
+ // one remaining integrated jet resolution
+ for (int ibin=0; ibin<_h_log10_R[m_njet]->size(); ++ibin) {
+ IDataPoint* dp=_h_log10_R[m_njet]->point(ibin);
+ double dcut=dp->coordinate(0)->value();
+ if (previous_dij>dcut) {
+ dp->coordinate(1)->setValue(dp->coordinate(1)->value()+weight);
+ }
+ }
+ }
+
+ const Jets& jets = jetpro.jetsByPt(20.0);
+
+ // the remaining direct jet observables
+ for (size_t i=0; i<m_njet; ++i) {
+ if (jets.size()<i+1) continue;
+ _h_pT_jet[i]->fill(jets[i].momentum().pT(), weight);
+ _h_eta_jet[i]->fill(jets[i].momentum().eta(), weight);
+
+ for (size_t j=i+1; j<m_njet; ++j) {
+ if (jets.size()<j+1) continue;
+ std::pair<size_t, size_t> ij(std::make_pair(i, j));
+ double deta = jets[i].momentum().eta()-jets[j].momentum().eta();
+ double dR = deltaR(jets[i].momentum(), jets[j].momentum());
+ _h_deta_jets[ij]->fill(deta, weight);
+ _h_dR_jets[ij]->fill(dR, weight);
+ }
+ }
+ _h_jet_multi_exclusive->fill(jets.size(), weight);
+
+ for (size_t i=0; i<m_njet+2; ++i) {
+ if (jets.size()>=i) {
+ _h_jet_multi_inclusive->fill(i, weight);
+ }
+ }
+ }
+
+
+ // Finalize
+ void MC_JetAnalysis::finalize() {
+ for (size_t i=0; i<m_njet; ++i) {
+ scale(_h_log10_d[i], crossSection()/sumOfWeights());
+ for (int ibin=0; ibin<_h_log10_R[i]->size(); ++ibin) {
+ IDataPoint* dp=_h_log10_R[i]->point(ibin);
+ dp->coordinate(1)->setValue(dp->coordinate(1)->value()*crossSection()/sumOfWeights());
+ }
+
+ scale(_h_pT_jet[i], crossSection()/sumOfWeights());
+ scale(_h_eta_jet[i], crossSection()/sumOfWeights());
+
+ for (size_t j=i+1; j<m_njet; ++j) {
+ }
+ }
+ for (int ibin=0; ibin<_h_log10_R[m_njet]->size(); ++ibin) {
+ IDataPoint* dp=_h_log10_R[m_njet]->point(ibin);
+ dp->coordinate(1)->setValue(dp->coordinate(1)->value()*crossSection()/sumOfWeights());
+ }
+
+ // fill inclusive jet multi ratio
+ int Nbins=_h_jet_multi_inclusive->axis().bins();
+ std::vector<double> ratio(Nbins-1, 0.0);
+ std::vector<double> err(Nbins-1, 0.0);
+ std::cout<<"Nbins="<<Nbins<<std::endl;
+ std::cout<<"m_njet+2="<<m_njet+2<<std::endl;
+ for (int i=0; i<Nbins-1; ++i) {
+ std::cout<<"i="<<i<<std::endl;
+ if (_h_jet_multi_inclusive->binHeight(i)>0.0 && _h_jet_multi_inclusive->binHeight(i+1)>0.0) {
+ ratio[i]=_h_jet_multi_inclusive->binHeight(i+1)/_h_jet_multi_inclusive->binHeight(i);
+ double relerr_i=_h_jet_multi_inclusive->binError(i)/_h_jet_multi_inclusive->binHeight(i);
+ double relerr_j=_h_jet_multi_inclusive->binError(i+1)/_h_jet_multi_inclusive->binHeight(i+1);
+ err[i]=ratio[i]*(relerr_i+relerr_j);
+ std::cout<<"ratio[i]="<<ratio[i]<<std::endl;
+ }
+ }
+ _h_jet_multi_ratio->setCoordinate(1, ratio, err);
+
+ scale(_h_jet_multi_exclusive, crossSection()/sumOfWeights());
+ scale(_h_jet_multi_inclusive, crossSection()/sumOfWeights());
+ }
+
+}
Modified: trunk/src/Analyses/Makefile.am
==============================================================================
--- trunk/src/Analyses/Makefile.am Wed Aug 19 17:46:01 2009 (r1768)
+++ trunk/src/Analyses/Makefile.am Wed Aug 19 17:52:54 2009 (r1769)
@@ -53,6 +53,7 @@
ZEUS_2001_S4815815.cc \
PDG_Hadron_Multiplicities.cc \
PDG_Hadron_Multiplicities_Ratios.cc \
+ MC_JetAnalysis.cc \
MC_TVT1960_ZJETS.cc \
MC_LHC_LEADINGJETS.cc \
MC_LHC_DIJET.cc \
More information about the Rivet-svn
mailing list
|