[Rivet-svn] r2773 - in trunk: . include/Rivet/Projections src/Analyses src/Projections

[blackhole at projects.hepforge.org]

Author: buckley
Date: Fri Nov 26 20:35:48 2010
New Revision: 2773

Log:
Making JetAlg return per-bin rather than per-event jet shape histo values, and updating the CDF analyses to match

Modified:
   trunk/ChangeLog
   trunk/include/Rivet/Projections/JetShape.hh
   trunk/src/Analyses/CDF_2005_S6217184.cc
   trunk/src/Analyses/CDF_2008_S7782535.cc
   trunk/src/Projections/JetShape.cc

Modified: trunk/ChangeLog
==============================================================================
--- trunk/ChangeLog	Fri Nov 26 18:47:27 2010	(r2772)
+++ trunk/ChangeLog	Fri Nov 26 20:35:48 2010	(r2773)
@@ -1,3 +1,14 @@
+2010-11-26  Andy Buckley  <andy at insectnation.org>
+
+	* CDF_2005_S6217184.cc, CDF_2008_S7782535.cc: Updates to use the
+	new per-jet JetAlg interface (and some other fixes).
+
+	* JetAlg.cc: Changed the interface on request to return per-jet
+	rather than per-event jet shapes, with an extra jet index argument.
+
+	* MathUtils.hh: Adding index_between(...) function, which is handy
+	for working out which bin a value falls in, given a set of bin edges.
+
 2010-11-25  Andy Buckley  <andy at insectnation.org>
 
 	* Cmp.hh: Adding ASC/DESC (and ANTISORTED) as preferred

Modified: trunk/include/Rivet/Projections/JetShape.hh
==============================================================================
--- trunk/include/Rivet/Projections/JetShape.hh	Fri Nov 26 18:47:27 2010	(r2772)
+++ trunk/include/Rivet/Projections/JetShape.hh	Fri Nov 26 20:35:48 2010	(r2773)
@@ -85,6 +85,11 @@
       return _binedges.size() - 1;
     }
 
+    /// Number of jets which passed cuts.
+    size_t numJets() const {
+      return _diffjetshapes.size();
+    }
+
     /// \f$ r_\text{min} \f$ value.
     double rMin() const {
       return _binedges.front();
@@ -106,40 +111,37 @@
     }
 
     /// Central \f$ r \f$ value for bin @a rbin.
-    /// @todo This external indexing thing is a bit nasty...
     double rBinMin(size_t rbin) const {
       assert(inRange(rbin, 0, numBins()));
       return _binedges[rbin];
     }
 
     /// Central \f$ r \f$ value for bin @a rbin.
-    /// @todo This external indexing thing is a bit nasty...
     double rBinMax(size_t rbin) const {
       assert(inRange(rbin, 0, numBins()));
       return _binedges[rbin+1];
     }
 
     /// Central \f$ r \f$ value for bin @a rbin.
-    /// @todo This external indexing thing is a bit nasty...
     double rBinMid(size_t rbin) const {
       assert(inRange(rbin, 0, numBins()));
       return (_binedges[rbin] + _binedges[rbin+1])/2.0;
     }
 
     /// Return value of differential jet shape profile histo bin.
-    /// @todo This external indexing thing is a bit nasty...
-    double diffJetShape(size_t rbin) const {
+    double diffJetShape(size_t ijet, size_t rbin) const {
+      assert(inRange(ijet, 0, numJets()));
       assert(inRange(rbin, 0, numBins()));
-      return _diffjetshapes[rbin];
+      return _diffjetshapes[ijet][rbin];
     }
 
     /// Return value of integrated jet shape profile histo bin.
-    /// @todo This external indexing thing is a bit nasty...
-    double intJetShape(size_t rbin) const {
+    double intJetShape(size_t ijet, size_t rbin) const {
+      assert(inRange(ijet, 0, numJets()));
       assert(inRange(rbin, 0, numBins()));
       double rtn  = 0;
       for (size_t i = 0; i <= rbin; ++i) {
-        rtn += _diffjetshapes[i];
+        rtn += _diffjetshapes[ijet][i];
       }
       return rtn;
     }
@@ -185,8 +187,8 @@
     /// @name The projected jet shapes
     //@{
 
-    /// Jet shape histo
-    vector<double> _diffjetshapes;
+    /// Jet shape histo -- first index is jet number, second is r bin
+    vector< vector<double> > _diffjetshapes;
 
     //@}
 

Modified: trunk/src/Analyses/CDF_2005_S6217184.cc
==============================================================================
--- trunk/src/Analyses/CDF_2005_S6217184.cc	Fri Nov 26 18:47:27 2010	(r2772)
+++ trunk/src/Analyses/CDF_2005_S6217184.cc	Fri Nov 26 20:35:48 2010	(r2773)
@@ -76,18 +76,15 @@
 
       for (size_t ipt = 0; ipt < 18; ++ipt) {
         const JetShape& jsipt = applyProjection<JetShape>(evt, _jsnames_pT[ipt]);
-        for (size_t rbin = 0; rbin < jsipt.numBins(); ++rbin) {
-          const double r_rho = jsipt.rBinMid(rbin);
-          cout << ipt << " " << rbin << " " << jsipt.diffJetShape(rbin) << endl;
-          _profhistRho_pT[ipt]->fill(r_rho/0.7, 0.7*jsipt.diffJetShape(rbin), weight);
-          const double r_Psi = jsipt.rBinMax(rbin);
-          _profhistPsi_pT[ipt]->fill(r_Psi/0.7, jsipt.intJetShape(rbin), weight);
+        for (size_t ijet = 0; ijet < jsipt.numJets(); ++ijet) {
+          for (size_t rbin = 0; rbin < jsipt.numBins(); ++rbin) {
+            const double r_rho = jsipt.rBinMid(rbin);
+            // cout << ipt << " " << rbin << " " << jsipt.diffJetShape(ijet, rbin) << endl;
+            _profhistRho_pT[ipt]->fill(r_rho/0.7, 0.7*jsipt.diffJetShape(ijet, rbin), weight);
+            const double r_Psi = jsipt.rBinMax(rbin);
+            _profhistPsi_pT[ipt]->fill(r_Psi/0.7, jsipt.intJetShape(ijet, rbin), weight);
+          }
         }
-
-        // Final histo is Psi(0.3/R) as a function of jet pT bin
-        /// @todo Can this actually be calculated event by event, or does it need to be assembled in a DPS in finalize? See CDF_2008_S7782535.
-        // const double ptmid = (_ptedges[ipt] + _ptedges[ipt+1])/2.0;
-        // _profhistPsi->fill(ptmid/GeV, jsipt.intJetShape(2), weight);
       }
 
     }

Modified: trunk/src/Analyses/CDF_2008_S7782535.cc
==============================================================================
--- trunk/src/Analyses/CDF_2008_S7782535.cc	Fri Nov 26 18:47:27 2010	(r2772)
+++ trunk/src/Analyses/CDF_2008_S7782535.cc	Fri Nov 26 20:35:48 2010	(r2773)
@@ -79,13 +79,14 @@
       const double weight = event.weight();
       for (size_t ipt = 0; ipt < 4; ++ipt) {
         if (bjets_ptbinned[ipt].empty()) continue;
-
         // Don't use the cached result: copy construct and calculate for provided b-jets only
         JetShape jsipt = applyProjection<JetShape>(event, _jsnames_pT[ipt]);
         jsipt.calc(bjets_ptbinned[ipt]);
-        for (size_t rbin = 0; rbin < jsipt.numBins(); ++rbin) {
-          const double r_Psi = jsipt.rBinMax(rbin);
-          _h_Psi_pT[ipt]->fill(r_Psi/0.7, jsipt.intJetShape(rbin), weight);
+        for (size_t ijet = 0; ijet < jsipt.numJets(); ++ijet) {
+          for (size_t rbin = 0; rbin < jsipt.numBins(); ++rbin) {
+            const double r_Psi = jsipt.rBinMax(rbin);
+            _h_Psi_pT[ipt]->fill(r_Psi/0.7, jsipt.intJetShape(ijet, rbin), weight);
+          }
         }
       }
 

Modified: trunk/src/Projections/JetShape.cc
==============================================================================
--- trunk/src/Projections/JetShape.cc	Fri Nov 26 18:47:27 2010	(r2772)
+++ trunk/src/Projections/JetShape.cc	Fri Nov 26 20:35:48 2010	(r2773)
@@ -55,7 +55,7 @@
 
 
   void JetShape::clear() {
-    _diffjetshapes = vector<double>(numBins(), 0.0);
+    _diffjetshapes.clear();
   }
 
 
@@ -63,27 +63,39 @@
     clear();
 
     foreach (const Jet& j, jets) {
+      // Apply jet cuts
       FourMomentum pj = j.momentum();
       if (!inRange(pj.pT(), _ptcuts)) continue;
       /// @todo Introduce a better (i.e. more safe and general) eta/y selection mechanism: MomentumFilter
       if (_rapscheme == PSEUDORAPIDITY && !inRange(fabs(pj.eta()), _rapcuts)) continue;
       if (_rapscheme == RAPIDITY && !inRange(fabs(pj.rapidity()), _rapcuts)) continue;
+
+      // Fill bins
+      vector<double> bins(numBins(), 0.0);
       foreach (const Particle& p, j.particles()) {
         const double dR = deltaR(pj, p.momentum(), _rapscheme);
         const int dRindex = index_between(dR, _binedges);
         if (dRindex == -1) continue; //< Out of histo range
-        _diffjetshapes[dRindex] += p.momentum().pT();
+        bins[dRindex] += p.momentum().pT();
       }
+
+      // Add bin vector for this jet to the diffjetshapes container
+      _diffjetshapes += bins;
     }
 
     // Normalize to total pT
-    double integral = 0.0;
-    for (size_t i = 0; i < numBins(); ++i) {
-      integral += _diffjetshapes[i];
-    }
-    if (integral > 0) {
+    foreach (vector<double>& binsref, _diffjetshapes) {
+      double integral = 0.0;
       for (size_t i = 0; i < numBins(); ++i) {
-        _diffjetshapes[i] /= integral;
+        integral += binsref[i];
+      }
+      if (integral > 0) {
+        for (size_t i = 0; i < numBins(); ++i) {
+          binsref[i] /= integral;
+        }
+      } else {
+        // It's just-about conceivable that a jet would have no particles in the given Delta(r) range...
+        MSG_DEBUG("No pT contributions in jet Delta(r) range: weird!");
       }
     }