[Rivet] Minor analysis changes

[Gavin Hesketh]

Hi,
A request: there are several analyses the plot differential cross 
sections like sigma(Z+jet) / sigma (Z). It would be very useful to also 
store sigma(Z) for these analyses (this is actually essential 
information to uses these analyses with alpgen).

The cases I know about:
D0_2010_S8671338.aida
D0_2009_S8349509.aida
D0_2008_S7863608.aida
D0_2009_S8202443.aida

I hacked these myself by adding a "normalisation" plot with one bin, 
just to hold the number of Z's. But maybe there is a better way to do this?
Below is an example, D0_2008_S7863608, which also needed updating to 
include the (Z+j)/(Z) plots as well as the absolute Z+jet plots (the 
stored aida file needs updating too).
thanks,
	Gavin

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Projections/ZFinder.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/RivetAIDA.hh"

namespace Rivet {


   /// @brief D0 differential Z/\f$ \gamma^* \f$ + jet + \f$ X \f$ cross 
sections
   /// @author Gavin Hesketh, Andy Buckley, Frank Siegert
   class D0_2008_S7863608 : public Analysis {
   public:

     /// @name Construction
     //@{
     /// Constructor
     D0_2008_S7863608() : Analysis("D0_2008_S7863608"), 
_sum_of_weights(0.0)
     {
       setBeams(PROTON, ANTIPROTON);
       setNeedsCrossSection(true);
     }

     //@}


     /// @name Analysis methods
     //@{

     /// Book histograms
     void init() {
       ZFinder zfinder(-1.7, 1.7, 15.0*GeV, MUON, 65.0*GeV, 115.0*GeV, 0.2);
       addProjection(zfinder, "ZFinder");

       FastJets conefinder(zfinder.remainingFinalState(), 
FastJets::D0ILCONE, 0.5);
       addProjection(conefinder, "ConeFinder");

       _h_jet_pT_cross_section = bookHistogram1D(1, 1, 1);
       _h_jet_y_cross_section = bookHistogram1D(2, 1, 1);
       _h_Z_pT_cross_section = bookHistogram1D(3, 1, 1);
       _h_Z_y_cross_section = bookHistogram1D(4, 1, 1);
       _h_total_cross_section = bookHistogram1D(5, 1, 1);

       _h_jet_pT_cross_section_norm = bookHistogram1D(1, 1, 2);
       _h_jet_y_cross_section_norm = bookHistogram1D(2, 1, 2);
       _h_Z_pT_cross_section_norm = bookHistogram1D(3, 1, 2);
       _h_Z_y_cross_section_norm = bookHistogram1D(4, 1, 2);

       _normalisation = bookHistogram1D("normalisation", 1, 0.0, 1.0);

     }



     // Do the analysis
     void analyze(const Event& e) {
       const double weight = e.weight();

       const ZFinder& zfinder = applyProjection<ZFinder>(e, "ZFinder");
       if (zfinder.particles().size()==1) {
	_sum_of_weights += weight;
	_normalisation->fill(0.5,weight);

         const JetAlg& jetpro = applyProjection<JetAlg>(e, "ConeFinder");
         const Jets& jets = jetpro.jetsByPt(20.0*GeV);
         Jets jets_cut;
         foreach (const Jet& j, jets) {
           if (fabs(j.momentum().pseudorapidity()) < 2.8) {
             jets_cut.push_back(j);
           }
         }

         // Return if there are no jets:
         if(jets_cut.size()<1) {
           getLog() << Log::DEBUG << "Skipping event " << 
e.genEvent().event_number()
                    << " because no jets pass cuts " << endl;
           vetoEvent;
         }

         // cut on Delta R between jet and muons
         foreach (const Jet& j, jets_cut) {
           foreach (const Particle& mu, 
zfinder.constituentsFinalState().particles()) {
             if (deltaR(mu.momentum().pseudorapidity(), 
mu.momentum().azimuthalAngle(),
                        j.momentum().pseudorapidity(), 
j.momentum().azimuthalAngle()) < 0.5) {
               vetoEvent;
             }
           }
         }

         const FourMomentum Zmom = zfinder.particles()[0].momentum();

         // In jet pT
         _h_jet_pT_cross_section->fill( jets_cut[0].momentum().pT(), 
weight);
         _h_jet_y_cross_section->fill( 
fabs(jets_cut[0].momentum().rapidity()), weight);
         _h_jet_pT_cross_section_norm->fill( 
jets_cut[0].momentum().pT(), weight);
         _h_jet_y_cross_section_norm->fill( 
fabs(jets_cut[0].momentum().rapidity()), weight);

         // In Z pT
         _h_Z_pT_cross_section->fill(Zmom.pT(), weight);
         _h_Z_y_cross_section->fill(fabs(Zmom.rapidity()), weight);
         _h_Z_pT_cross_section_norm->fill(Zmom.pT(), weight);
         _h_Z_y_cross_section_norm->fill(fabs(Zmom.rapidity()), weight);

         _h_total_cross_section->fill(1960.0, weight);
       }
     }



     /// Finalize
     void finalize() {
       const double invlumi = crossSection()/sumOfWeights();
       scale(_h_total_cross_section, invlumi);
       scale(_h_jet_pT_cross_section, invlumi);
       scale(_h_jet_y_cross_section, invlumi);
       scale(_h_Z_pT_cross_section, invlumi);
       scale(_h_Z_y_cross_section, invlumi);

       if(_sum_of_weights>0){
	scale(_h_jet_pT_cross_section_norm, 1.0/_sum_of_weights);
	scale(_h_jet_y_cross_section_norm, 1.0/_sum_of_weights);
	scale(_h_Z_pT_cross_section_norm, 1.0/_sum_of_weights);
	scale(_h_Z_y_cross_section_norm, 1.0/_sum_of_weights);
       }
     }

     //@}


   private:

     /// @name Histograms
     //@{
     AIDA::IHistogram1D * _h_jet_pT_cross_section;
     AIDA::IHistogram1D * _h_jet_y_cross_section;
     AIDA::IHistogram1D * _h_Z_pT_cross_section;
     AIDA::IHistogram1D * _h_Z_y_cross_section;
     AIDA::IHistogram1D * _h_total_cross_section;

     AIDA::IHistogram1D * _h_jet_pT_cross_section_norm;
     AIDA::IHistogram1D * _h_jet_y_cross_section_norm;
     AIDA::IHistogram1D * _h_Z_pT_cross_section_norm;
     AIDA::IHistogram1D * _h_Z_y_cross_section_norm;

     AIDA::IHistogram1D *_normalisation;

     //@}
     double _sum_of_weights;
   };



   // This global object acts as a hook for the plugin system
   AnalysisBuilder<D0_2008_S7863608> plugin_D0_2008_S7863608;

}