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

[blackhole at projects.hepforge.org]

Author: buckley
Date: Thu Nov 25 15:25:34 2010
New Revision: 2767

Log:
Jet shape fixes and hacking in a temporary |y| restriction: we need a more powerful momentum filtering mechanism

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

Modified: trunk/include/Rivet/Projections/JetShape.hh
==============================================================================
--- trunk/include/Rivet/Projections/JetShape.hh	Thu Nov 25 14:11:16 2010	(r2766)
+++ trunk/include/Rivet/Projections/JetShape.hh	Thu Nov 25 15:25:34 2010	(r2767)
@@ -54,13 +54,13 @@
     JetShape(const JetAlg& jetalg,
              double rmin, double rmax, size_t nbins,
              double ptmin=0, double ptmax=MAXDOUBLE,
-             double rapmin=-MAXDOUBLE, double rapmax=-MAXDOUBLE,
+             double absrapmin=-MAXDOUBLE, double absrapmax=-MAXDOUBLE,
              RapScheme rapscheme=RAPIDITY);
 
     /// Constructor from vector of bin edges.
     JetShape(const JetAlg& jetalg, vector<double> binedges,
              double ptmin=0, double ptmax=MAXDOUBLE,
-             double rapmin=-MAXDOUBLE, double rapmax=-MAXDOUBLE,
+             double absrapmin=-MAXDOUBLE, double absrapmax=-MAXDOUBLE,
              RapScheme rapscheme=RAPIDITY);
 
     /// Clone on the heap.

Modified: trunk/src/Analyses/CDF_2005_S6217184.cc
==============================================================================
--- trunk/src/Analyses/CDF_2005_S6217184.cc	Thu Nov 25 14:11:16 2010	(r2766)
+++ trunk/src/Analyses/CDF_2005_S6217184.cc	Thu Nov 25 15:25:34 2010	(r2767)
@@ -11,9 +11,7 @@
 
 
   /// @brief CDF Run II jet shape analysis
-  /// @author Lars Sonnenschein
   /// @author Andy Buckley
-  /// @todo Convert to use the correct jet shape algorithm
   class CDF_2005_S6217184 : public Analysis {
   public:
 
@@ -24,6 +22,7 @@
       setBeams(PROTON, ANTIPROTON);
     }
 
+
     /// @name Analysis methods
     //@{
 
@@ -36,7 +35,7 @@
       addProjection(fj, "Jets");
 
       // Specify pT bins
-      _ptedges = { 37.0, 45.0, 55.0, 63.0, 73.0, 84.0, 97.0, 112.0, 128.0,
+      _ptedges = { 37.0, 45.0, 55.0, 63.0, 73.0, 84.0, 97.0, 112.0, 128.0, \
                    148.0, 166.0, 186.0, 208.0, 229.0, 250.0, 277.0, 304.0, 340.0, 380.0 };
 
       // Register a jet shape projection and histogram for each pT bin
@@ -60,30 +59,29 @@
 
 
     /// Do the analysis
-    void analyze(const Event& event) {
+    void analyze(const Event& evt) {
 
       // Get jets and require at least one to pass pT and y cuts
-      /// @todo Need a different rapidity cut -- discontinuous in |y|
-      const Jets jets = applyProjection<FastJets>(event, "Jets").jetsByPt(37, 380, 0.1, 0.7);
-      getLog() << Log::DEBUG << "Jet multiplicity before cuts = " << jets.size() << endl;
+      const Jets jets = applyProjection<FastJets>(evt, "Jets").jetsByPt(37, 380, 0.1, 0.7);
+      MSG_DEBUG("Jet multiplicity before cuts = " << jets.size());
       if (jets.size() == 0) {
         MSG_DEBUG("No jets found in required pT range");
-        vetoEvent(event);
+        vetoEvent(evt);
       }
 
       // Calculate and histogram jet shapes
-      const double weight = event.weight();
+      const double weight = evt.weight();
 
-      for (size_t ipT = 0; ipT < 18; ++ipT) {
-        const JetShape& jsipt = applyProjection<JetShape>(event, _jsnames_pT[ipT]);
+      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);
-          _profhistRho_pT[ipT]->fill(r_rho/0.7, (0.7/1.0)*jsipt.diffJetShape(rbin), weight);
+          _profhistRho_pT[ipt]->fill(r_rho/0.7, (0.7/1.0)*jsipt.diffJetShape(rbin), weight);
           const double r_Psi = jsipt.rBinMax(rbin);
-          _profhistPsi_pT[ipT]->fill(r_Psi/0.7, jsipt.intJetShape(rbin), weight);
+          _profhistPsi_pT[ipt]->fill(r_Psi/0.7, jsipt.intJetShape(rbin), weight);
         }
         // Final histo is 1 - Psi(0.3/R) as a function of jet pT bin
-        const double ptmid = (_ptedges[ipT] + _ptedges[ipT+1])/2.0;
+        const double ptmid = (_ptedges[ipt] + _ptedges[ipt+1])/2.0;
         _profhistPsi->fill(ptmid/GeV, jsipt.intJetShape(2), weight);
       }
 

Modified: trunk/src/Projections/JetShape.cc
==============================================================================
--- trunk/src/Projections/JetShape.cc	Thu Nov 25 14:11:16 2010	(r2766)
+++ trunk/src/Projections/JetShape.cc	Thu Nov 25 15:25:34 2010	(r2767)
@@ -9,13 +9,14 @@
   JetShape::JetShape(const JetAlg& jetalg,
                      double rmin, double rmax, size_t nbins,
                      double ptmin, double ptmax,
-                     double rapmin, double rapmax,
+                     double absrapmin, double absrapmax,
                      RapScheme rapscheme)
     : _rapscheme(rapscheme)
   {
     setName("JetShape");
     _binedges = linspace(rmin, rmax, nbins);
     _ptcuts = make_pair(ptmin, ptmax);
+    _rapcuts = make_pair(absrapmin, absrapmax);
     addProjection(jetalg, "Jets");
   }
 
@@ -24,12 +25,13 @@
   JetShape::JetShape(const JetAlg& jetalg,
                      vector<double> binedges,
                      double ptmin, double ptmax,
-                     double rapmin, double rapmax,
+                     double absrapmin, double absrapmax,
                      RapScheme rapscheme)
     : _binedges(binedges), _rapscheme(rapscheme)
   {
     setName("JetShape");
     _ptcuts = make_pair(ptmin, ptmax);
+    _rapcuts = make_pair(absrapmin, absrapmax);
     addProjection(jetalg, "Jets");
   }
 
@@ -60,10 +62,10 @@
 
     foreach (const Jet& j, jets) {
       FourMomentum pj = j.momentum();
-      // These cuts are already applied when getting the jets in ::project(), but we check again for direct calc use.
       if (!inRange(pj.pT(), _ptcuts)) continue;
-      if (_rapscheme == PSEUDORAPIDITY && !inRange(pj.eta(), _rapcuts)) continue;
-      if (_rapscheme == RAPIDITY && !inRange(pj.rapidity(), _rapcuts)) 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;
       foreach (const Particle& p, j.particles()) {
         const double dR = deltaR(pj, p.momentum(), _rapscheme);
         if (!inRange(dR, _binedges.front(), _binedges.back())) continue; //< Out of histo range
@@ -93,7 +95,7 @@
 
   void JetShape::project(const Event& e) {
     const Jets jets = applyProjection<JetAlg>(e, "Jets").jets(_ptcuts.first, _ptcuts.second,
-                                                              _rapcuts.first, _rapcuts.second, _rapscheme);
+                                                              -_rapcuts.second, _rapcuts.second, _rapscheme);
     calc(jets);
   }