[Rivet-svn] r2671 - in trunk: include/Rivet/Analyses src/Analyses

[blackhole at projects.hepforge.org]

Author: buckley
Date: Sun Aug 22 16:15:01 2010
New Revision: 2671

Log:
Adding new jet eta/y ratio plots (for +/- rapidities) and making histo booking a bit more robust in MC_GENERIC by using the histoPath(name) function

Modified:
   trunk/include/Rivet/Analyses/MC_JetAnalysis.hh
   trunk/src/Analyses/MC_GENERIC.cc
   trunk/src/Analyses/MC_JetAnalysis.cc

Modified: trunk/include/Rivet/Analyses/MC_JetAnalysis.hh
==============================================================================
--- trunk/include/Rivet/Analyses/MC_JetAnalysis.hh	Sun Aug 22 14:26:25 2010	(r2670)
+++ trunk/include/Rivet/Analyses/MC_JetAnalysis.hh	Sun Aug 22 16:15:01 2010	(r2671)
@@ -4,6 +4,7 @@
 
 #include "Rivet/Analysis.hh"
 #include "Rivet/Projections/FinalState.hh"
+#include "LWH/Histogram1D.h"
 
 namespace Rivet {
 
@@ -47,7 +48,9 @@
     std::vector<AIDA::IDataPointSet *> _h_log10_R;
     std::vector<AIDA::IHistogram1D *> _h_pT_jet;
     std::vector<AIDA::IHistogram1D *> _h_eta_jet;
+    std::vector<shared_ptr<LWH::Histogram1D> > _h_eta_jet_plus, _h_eta_jet_minus;
     std::vector<AIDA::IHistogram1D *> _h_rap_jet;
+    std::vector<shared_ptr<LWH::Histogram1D> > _h_rap_jet_plus, _h_rap_jet_minus;
     std::vector<AIDA::IHistogram1D *> _h_mass_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;

Modified: trunk/src/Analyses/MC_GENERIC.cc
==============================================================================
--- trunk/src/Analyses/MC_GENERIC.cc	Sun Aug 22 14:26:25 2010	(r2670)
+++ trunk/src/Analyses/MC_GENERIC.cc	Sun Aug 22 16:15:01 2010	(r2671)
@@ -133,10 +133,10 @@
       scale(_histPhi, 1/sumOfWeights());
       scale(_histPhiCh, 1/sumOfWeights());
 
-      histogramFactory().divide("/MC_GENERIC/EtaPMRatio", *_tmphistEtaPlus, *_tmphistEtaMinus);
-      histogramFactory().divide("/MC_GENERIC/EtaChPMRatio", *_tmphistEtaChPlus, *_tmphistEtaChMinus);
-      histogramFactory().divide("/MC_GENERIC/RapidityPMRatio", *_tmphistRapPlus, *_tmphistRapMinus);
-      histogramFactory().divide("/MC_GENERIC/RapidityChPMRatio", *_tmphistRapChPlus, *_tmphistRapChMinus);
+      histogramFactory().divide(histoPath("EtaPMRatio"), *_tmphistEtaPlus, *_tmphistEtaMinus);
+      histogramFactory().divide(histoPath("EtaChPMRatio"), *_tmphistEtaChPlus, *_tmphistEtaChMinus);
+      histogramFactory().divide(histoPath("RapidityPMRatio"), *_tmphistRapPlus, *_tmphistRapMinus);
+      histogramFactory().divide(histoPath("RapidityChPMRatio"), *_tmphistRapChPlus, *_tmphistRapChMinus);
     }
 
     //@}

Modified: trunk/src/Analyses/MC_JetAnalysis.cc
==============================================================================
--- trunk/src/Analyses/MC_JetAnalysis.cc	Sun Aug 22 14:26:25 2010	(r2670)
+++ trunk/src/Analyses/MC_JetAnalysis.cc	Sun Aug 22 16:15:01 2010	(r2671)
@@ -13,7 +13,9 @@
                                  double jetptcut)
     : Analysis(name), m_njet(njet), m_jetpro_name(jetpro_name), m_jetptcut(jetptcut),
       _h_log10_d(njet, NULL), _h_log10_R(njet+1, NULL), _h_pT_jet(njet, NULL),
-      _h_eta_jet(njet, NULL), _h_rap_jet(njet, NULL), _h_mass_jet(njet, NULL)
+      _h_eta_jet(njet, NULL), _h_eta_jet_plus(njet), _h_eta_jet_minus(njet),
+      _h_rap_jet(njet, NULL), _h_rap_jet_plus(njet), _h_rap_jet_minus(njet),
+      _h_mass_jet(njet, NULL)
   {
     setNeedsCrossSection(true);
   }
@@ -35,7 +37,7 @@
 
       stringstream pTname;
       pTname << "jet_pT_" << i+1;
-      double pTmax = 1.0/(double(i)+2.0)*sqrtS()/GeV/2.0;
+      double pTmax = 1.0/(double(i)+2.0) * sqrtS()/GeV/2.0;
       int nbins_pT = 100/(i+1);
       _h_pT_jet[i] = bookHistogram1D(pTname.str(), logBinEdges(nbins_pT, 10.0, pTmax));
 
@@ -47,14 +49,18 @@
 
       stringstream etaname;
       etaname << "jet_eta_" << i+1;
-      _h_eta_jet[i] = bookHistogram1D(etaname.str(), i>1 ? 25 : 50, -5.0, 5.0);
+      _h_eta_jet[i] = bookHistogram1D(etaname.str(), i > 1 ? 25 : 50, -5.0, 5.0);
+      _h_eta_jet_plus[i].reset(new LWH::Histogram1D(i > 1 ? 15 : 25, 0, 5));
+      _h_eta_jet_minus[i].reset(new LWH::Histogram1D(i > 1 ? 15 : 25, 0, 5));
 
       stringstream rapname;
       rapname << "jet_y_" << i+1;
       _h_rap_jet[i] = bookHistogram1D(rapname.str(), i>1 ? 25 : 50, -5.0, 5.0);
+      _h_rap_jet_plus[i].reset(new LWH::Histogram1D(i > 1 ? 15 : 25, 0, 5));
+      _h_rap_jet_minus[i].reset(new LWH::Histogram1D(i > 1 ? 15 : 25, 0, 5));
 
-      for (size_t j=i+1; j<m_njet; ++j) {
-        std::pair<size_t, size_t> ij(std::make_pair(i, j));
+      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 << "jets_deta_" << i+1 << j+1;
@@ -117,7 +123,7 @@
 
     // The remaining direct jet observables
     for (size_t i = 0; i < m_njet; ++i) {
-      if (jets.size()<i+1) continue;
+      if (jets.size() < i+1) continue;
       _h_pT_jet[i]->fill(jets[i].momentum().pT()/GeV, weight);
       // Check for numerical precision issues with jet masses
       double m2_i = jets[i].momentum().mass2();
@@ -128,13 +134,30 @@
         }
         m2_i = 0.0;
       }
+
+      // Jet mass
       _h_mass_jet[i]->fill(sqrt(m2_i)/GeV, weight);
+
+      // Jet eta
       _h_eta_jet[i]->fill(jets[i].momentum().eta(), weight);
+      if (jets[i].momentum().eta() > 0.0) {
+        _h_eta_jet_plus[i]->fill(jets[i].momentum().eta(), weight);
+      } else {
+        _h_eta_jet_minus[i]->fill(jets[i].momentum().eta(), weight);
+      }
+
+      // Jet rapidity
       _h_rap_jet[i]->fill(jets[i].momentum().rapidity(), weight);
+      if (jets[i].momentum().rapidity() > 0.0) {
+        _h_rap_jet_plus[i]->fill(jets[i].momentum().rapidity(), weight);
+      } else {
+        _h_rap_jet_minus[i]->fill(jets[i].momentum().rapidity(), weight);
+      }
 
+      // Inter-jet properties
       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));
+        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);
@@ -165,7 +188,16 @@
       scale(_h_eta_jet[i], crossSection()/sumOfWeights());
       scale(_h_rap_jet[i], crossSection()/sumOfWeights());
 
+      // Create eta/rapidity ratio plots
+      stringstream etaname;
+      etaname << "jet_eta_pmratio_" << i;
+      histogramFactory().divide(histoPath(etaname.str()), *_h_eta_jet_plus[i], *_h_eta_jet_minus[i]);
+      stringstream rapname;
+      rapname << "jet_rap_pmratio_" << i;
+      histogramFactory().divide(histoPath(rapname.str()), *_h_rap_jet_plus[i], *_h_rap_jet_minus[i]);
     }
+
+
     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());