[Rivet-svn] r2876 - in trunk: . bin include/Rivet/Tools src/Analyses src/Core src/Tools

blackhole at projects.hepforge.org blackhole at projects.hepforge.org
Wed Jan 12 22:55:31 GMT 2011


Author: buckley
Date: Wed Jan 12 22:55:30 2011
New Revision: 2876

Log:
Many improvements to compare-histos and rivet-mkhtml, and adding getAnalysisPlotPaths()

Modified:
   trunk/ChangeLog
   trunk/bin/compare-histos
   trunk/bin/rivet-mkhtml
   trunk/include/Rivet/Tools/BinnedHistogram.hh
   trunk/include/Rivet/Tools/RivetPaths.hh
   trunk/src/Analyses/ATLAS_2010_S8817804.cc
   trunk/src/Core/AnalysisInfo.cc
   trunk/src/Tools/RivetPaths.cc

Modified: trunk/ChangeLog
==============================================================================
--- trunk/ChangeLog	Tue Jan 11 14:59:39 2011	(r2875)
+++ trunk/ChangeLog	Wed Jan 12 22:55:30 2011	(r2876)
@@ -1,3 +1,18 @@
+2011-01-12  Andy Buckley  <andy at insectnation.org>
+
+	* Fix several rendering and comparison bugs in rivet-mkhtml.
+
+	* Allow make-plots to write into an existing directory, at the
+	user's own risk.
+
+	* Make rivet-mkhtml produce PDF-based output rather than PS by
+	default (most people want PDF these days). Can we do the same
+	change of default for make-plots?
+
+	* Add getAnalysisPlotPaths() function, and use it in compare-histos
+
+	* Use proper .info file search path function internally in AnalysisInfo::make.
+
 2011-01-11  Andy Buckley  <andy at insectnation.org>
 
 	* Clean up ATLAS dijet analysis.

Modified: trunk/bin/compare-histos
==============================================================================
--- trunk/bin/compare-histos	Tue Jan 11 14:59:39 2011	(r2875)
+++ trunk/bin/compare-histos	Wed Jan 12 22:55:30 2011	(r2876)
@@ -143,11 +143,8 @@
         opts.SHOW_IF_REF_ONLY = True
 
 
-    ## Add RIVET_PLOT_PATH and ref data dirs to PLOTINFO path
-    if os.environ.has_key("RIVET_PLOT_PATH"):
-        opts.PLOTINFODIR += [os.path.abspath(i) for i in os.environ["RIVET_PLOT_PATH"].split(":") if i]
-    opts.PLOTINFODIR += rivet_data_dirs
-    #print opts.PLOTINFODIR
+    ## Add standard locations and the input files' dirs to the PLOTINFO search paths
+    opts.PLOTINFODIR += rivet.getAnalysisPlotPaths()
     for a in args:
         adir = os.path.abspath(os.path.split(a)[0])
         if not adir in opts.PLOTINFODIR:

Modified: trunk/bin/rivet-mkhtml
==============================================================================
--- trunk/bin/rivet-mkhtml	Tue Jan 11 14:59:39 2011	(r2875)
+++ trunk/bin/rivet-mkhtml	Wed Jan 12 22:55:30 2011	(r2876)
@@ -3,15 +3,20 @@
 """\
 %prog [options] <aidafile1> [<aidafile2> <aidafile3>...]
 
-Make web pages from histogram files written out by Rivet.
-You can specify multiple Monte Carlo AIDA files to be compared in the same
-syntax as for compare-histos, i.e. including plotting options.
-Reference data should be found automatically.
+Make web pages from histogram files written out by Rivet.  You can specify
+multiple Monte Carlo AIDA files to be compared in the same syntax as for
+compare-histos, i.e. including plotting options.
+
+Reference data, analysis metadata, and plot style information should be found
+automatically (if not, set the RIVET_ANALYSIS_PATH or similar variables
+appropriately).
+
+You can overwrite an existing output directory.
 """
 
 import sys
 if sys.version_info[:3] < (2,4,0):
-    print >> sys.stderr, "rivet scripts require Python version >= 2.4.0... exiting"
+    sys.stderr.write("rivet scripts require Python version >= 2.4.0... exiting\n")
     sys.exit(1)
 
 import traceback
@@ -19,7 +24,7 @@
     import rivet
 except ImportError:
     traceback.print_exc(file=sys.stderr)
-    print >> sys.stderr, "rivet is broken... exiting"
+    sys.stderr.write("rivet is broken... exiting\n")
     sys.exit(1)
 
 import sys, os, glob, shutil
@@ -41,9 +46,9 @@
 parser.add_option("-n", "--num-threads", metavar="NUMTHREADS", dest="NUMTHREADS", type=int,
                   default=None, help="request make-plots to use a specific number of threads.")
 parser.add_option("--pdf", dest="VECTORFORMAT", action="store_const", const="PDF",
-                  default="PS", help="use PDF as the vector plot format.")
+                  default="PDF", help="use PDF as the vector plot format.")
 parser.add_option("--ps", dest="VECTORFORMAT", action="store_const", const="PS",
-                  default="PS", help="use PostScript as the vector plot format.")
+                  default="PDF", help="use PostScript as the vector plot format.")
 parser.add_option("-i", "--ignore-unvalidated", dest="IGNORE_UNVALIDATED", action="store_true",
                   default=False, help="ignore unvalidated analyses.")
 parser.add_option("-m", "--match", action="append", dest="PATHPATTERNS",
@@ -64,6 +69,9 @@
 
 
 ## Make output directory
+if os.path.exists(opts.OUTPUTDIR):
+    import shutil
+    shutil.rmtree(opts.OUTPUTDIR)
 try:
     os.makedirs(opts.OUTPUTDIR)
 except:
@@ -94,7 +102,11 @@
     if not os.access(aidafilepath, os.R_OK):
         print "Error: cannot read from %s" % aidafilepath
         sys.exit(2)
-    tree = ET.parse(aidafilepath)
+    try:
+        tree = ET.parse(aidafilepath)
+    except Exception, e:
+        print "Problem parsing AIDA XML file '%s': %s. Skipping this file" % (aidafilepath, e)
+        continue
     for dps in tree.findall("dataPointSet"):
         path = dps.get("path")
         analysis = path[path.rfind("/")+1:]
@@ -128,6 +140,7 @@
     ch_cmd.append("--refid=%s" % os.path.abspath(opts.REF_ID))
 ch_cmd.append("--hier-out")
 ch_cmd.append("--rivet-refs")
+# TODO: This isn't very sensible... what's the intention? Provide --plotinfodir cmd line option?
 ch_cmd.append("--plotinfodir=../")
 for af in aidafiles:
     ch_cmd.append("%s" % os.path.abspath(af))
@@ -135,6 +148,7 @@
     ch_cmd.append("--verbose")
     print "Calling compare-histos with the following options:"
     print ch_cmd
+    print " ".join(ch_cmd)
 Popen(ch_cmd, cwd=opts.OUTPUTDIR, stderr=PIPE).wait()
 
 
@@ -147,8 +161,6 @@
   a { text-decoration: none; font-weight: bold; }
 </style>"""
 
-import sys
-
 index = open(os.path.join(opts.OUTPUTDIR, "index.html"), "w")
 index.write('<html>\n<head>\n<title>%s</title>\n%s</head>\n<body>' % (opts.OUTPUTDIR, style))
 index.write('<h2>Plots from Rivet analyses</h2>\n\n')
@@ -170,28 +182,34 @@
     references = []
     summary = analysis
     description = "NONE"
-    spiresid = analysis[analysis.rfind('S')+1:len(analysis)]
+    spiresid = analysis[analysis.rfind('S')+1:len(analysis)] if analysis.startswith("S") else "NONE"
     ana = rivet.AnalysisLoader.getAnalysis(analysis)
     if ana:
-        if ana.summary() is not "NONE":
+        if ana.summary() and ana.summary() != "NONE":
             summary = "%s (%s)" % (ana.summary(), analysis)
         references = ana.references()
         description = ana.description()
         spiresid = ana.spiresId()
-        if opts.IGNORE_UNVALIDATED and ana.status() != "VALIDATED":
+        if opts.IGNORE_UNVALIDATED and ana.status().upper() != "VALIDATED":
             continue
     if opts.SINGLE:
         index.write('<h3 style="clear:left; padding-top:2em;"><a name="%s">%s</a></h3>\n' % (analysis, summary) )
     else:
         index.write('<h3><a href="%s/index.html" style="text-decoration:none;">%s</a></h3>\n' % (analysis, summary))
-    if spiresid and spiresid is not "NONE":
-        index.write('<p><a href="http://durpdg.dur.ac.uk/cgi-bin/spiface/hep/www?irn+%s">Spires</a>' % spiresid)
-    for ref in references:
-        index.write(' | %s' % ref)
-    index.write('</p>\n')
-    index.write('<p style="font-size:small;">%s</p>\n' % description)
+    reflist = []
+    if spiresid and spiresid != "NONE":
+        reflist.append('<a href="http://durpdg.dur.ac.uk/cgi-bin/spiface/hep/www?irn+%s">Spires</a>' % spiresid)
+    reflist += references
+    index.write('<p>%s</p>\n' % " | ".join(reflist))
+    index.write('<p style="font-size:smaller;">%s</p>\n' % description)
     anapath = os.path.join(opts.OUTPUTDIR, analysis)
     if not opts.SINGLE:
+        if not os.path.exists(anapath):
+            try:
+                os.makedirs(anapath)
+            except:
+                print "Error: failed to make new directory '%s'. Skipping analysis %s" % (anapath, analysis)
+                continue
         anaindex = open(os.path.join(anapath, "index.html"), 'w')
         anaindex.write('<html>\n<head>\n<title>%s - %s</title>\n%s</head>\n<body>\n' %
                        (opts.OUTPUTDIR, analysis, style))
@@ -207,8 +225,8 @@
         pngfile = obsname+".png"
         vecfile = obsname+"."+opts.VECTORFORMAT.lower()
         if opts.SINGLE:
-            pngfile = analysis+"/"+pngfile
-            vecfile = analysis+"/"+vecfile
+            pngfile = os.path.join(analysis, pngfile)
+            vecfile = os.path.join(analysis, vecfile)
 
         anaindex.write('  <div style="float:left; font-size:smaller; font-weight:bold;">\n')
         anaindex.write('    <a href="#%s-%s">⚓</a> %s:<br>\n' % (analysis, obsname, vecfile) )
@@ -219,13 +237,11 @@
 
     if not opts.SINGLE:
         anaindex.write("</body>\n</html>\n")
-# TODO: Why doesn't the index file work until late in the image processing?
 index.write('</body>\n</html>')
 
 
-## Now make the plots
-
 ## Run make-plots on all generated .dat files
+# sys.exit(0)
 mp_cmd = ["make-plots"]
 if opts.NUMTHREADS:
     mp_cmd.append("--num-threads=%d" % opts.NUMTHREADS)
@@ -238,12 +254,16 @@
     if os.access(os.path.expanduser(configfile), os.R_OK):
         mp_cmd.append("-c")
         mp_cmd.append(os.path.expanduser(configfile))
+datfiles = []
 for analysis in sorted(analyses, reverse=True):
     anapath = os.path.join(opts.OUTPUTDIR, analysis)
-    datfiles = glob.glob("%s/*.dat" % anapath)
-    for datfile in sorted(datfiles):
-        mp_cmd.append(datfile)
-if opts.VERBOSE:
-    print "Calling make-plots with the following options:"
-    print mp_cmd
-Popen(mp_cmd).wait()
+    #print anapath
+    anadatfiles = glob.glob("%s/*.dat" % anapath)
+    datfiles += sorted(anadatfiles)
+if datfiles:
+    mp_cmd += datfiles
+    if opts.VERBOSE:
+        mp_cmd.append("--verbose")
+        print "Calling make-plots with the following options:"
+        print mp_cmd
+    Popen(mp_cmd).wait()

Modified: trunk/include/Rivet/Tools/BinnedHistogram.hh
==============================================================================
--- trunk/include/Rivet/Tools/BinnedHistogram.hh	Tue Jan 11 14:59:39 2011	(r2875)
+++ trunk/include/Rivet/Tools/BinnedHistogram.hh	Wed Jan 12 22:55:30 2011	(r2876)
@@ -4,7 +4,7 @@
 #include "Rivet/Rivet.hh"
 
 namespace Rivet {
-  
+
   class Analysis;
 
 
@@ -12,44 +12,44 @@
    * BinnedHistogram contains a series of histograms of the same quantity
    * each in a different region of a second quantity.  For example, a
    * BinnedHistogram may contain histograms of the cross section differential
-   * in PT in different eta regions.
+   * in \f$ p_T \f$ in different \f$ \eta \f$  regions.
    **/
-  template<typename T> 
+  template<typename T>
   class BinnedHistogram {
   public:
-    
+
     /// Create a new empty BinnedHistogram
     BinnedHistogram() {
       return;
     }
- 
+
     ///  Add a histogram in the region between @a binMin and @a binMax to this
     ///  set of BinnedHistograms.
     const BinnedHistogram<T>& addHistogram(const T& binMin,
                                            const T& binMax,
                                            AIDA::IHistogram1D* histo);
-    
+
     /// Fill the histogram that lies in the same region as @a bin with the value
     /// @a val of weight @a weight.
     AIDA::IHistogram1D* fill(const T& bin,
                              const T& val,
                              double weight);
-    
+
     /// Scale histograms taking into account its "external" binwidth, i.e. by
     /// scale/binWidth
     void scale(const T& scale, Analysis* ana);
-    
+
     const vector<AIDA::IHistogram1D*>& getHistograms() const { return _histos; }
     vector<AIDA::IHistogram1D*>& getHistograms() { return _histos; }
- 
+
 
   private:
- 
+
     map<T, AIDA::IHistogram1D*> _histosByUpperBound;
     map<T, AIDA::IHistogram1D*> _histosByLowerBound;
     vector<AIDA::IHistogram1D*> _histos;
     map<AIDA::IHistogram1D*, T> _binWidths;
- 
+
   };
 
 

Modified: trunk/include/Rivet/Tools/RivetPaths.hh
==============================================================================
--- trunk/include/Rivet/Tools/RivetPaths.hh	Tue Jan 11 14:59:39 2011	(r2875)
+++ trunk/include/Rivet/Tools/RivetPaths.hh	Wed Jan 12 22:55:30 2011	(r2876)
@@ -30,6 +30,9 @@
   /// Get Rivet analysis info metadata search paths
   const std::vector<std::string> getAnalysisInfoPaths();
 
+  /// Get Rivet analysis plot style search paths
+  const std::vector<std::string> getAnalysisPlotPaths();
+
 
 }
 

Modified: trunk/src/Analyses/ATLAS_2010_S8817804.cc
==============================================================================
--- trunk/src/Analyses/ATLAS_2010_S8817804.cc	Tue Jan 11 14:59:39 2011	(r2875)
+++ trunk/src/Analyses/ATLAS_2010_S8817804.cc	Wed Jan 12 22:55:30 2011	(r2876)
@@ -10,14 +10,14 @@
 
 
   /// @brief ATLAS inclusive jet pT spectrum, di-jet Mass and di-jet chi
-  class ATLAS_2010_S8817804: public Analysis{
+  class ATLAS_2010_S8817804: public Analysis {
   public:
 
-    ATLAS_2010_S8817804()
-      : Analysis("ATLAS_2010_S8817804")
+    ATLAS_2010_S8817804() : Analysis("ATLAS_2010_S8817804")
     {
       setBeams(PROTON, PROTON);
       setNeedsCrossSection(true);
+      _sumW = 0;
     }
 
 
@@ -64,34 +64,47 @@
       jetAr[AKT4] = applyProjection<FastJets>(evt, "AntiKT06").jetsByPt(30*GeV);
       jetAr[AKT6] = applyProjection<FastJets>(evt, "AntiKT04").jetsByPt(30*GeV);
 
-      // Cut values
-      const double rapMax(2.8), leadPTCut(60*GeV), asymmetryCut(log(30.)), ySumCut(2.2);
+      cout << -1 << endl;
 
+      // Identify the dijets
       for (size_t alg = 0; alg < 2; ++alg) {
+        cout << 0 << endl;
+
         vector<FourMomentum> leadjets;
         foreach (const Jet& jet, jetAr[alg]) {
+          cout << 0.1 << endl;
           double pT = jet.momentum().pT();
-          double fabsy = fabs(jet.momentum().rapidity());
-          _pTHistos[alg].fill(fabsy, pT/GeV, evt.weight());
+          double absy = fabs(jet.momentum().rapidity());
+          cout << 0.2 << endl;
+          _pTHistos[alg].fill(absy, pT/GeV, evt.weight());
+
+          cout << 0.3 << endl;
 
-          if (fabsy < rapMax && leadjets.size() < 2) {
-            if (leadjets.empty() && pT < leadPTCut) continue;
+          if (absy < 2.8 && leadjets.size() < 2) {
+            if (leadjets.empty() && pT < 60*GeV) continue;
             leadjets.push_back(jet.momentum());
           }
+          cout << 0.4 << endl;
         }
 
+        cout << 1 << endl;
+
+        // Veto if no acceptable dijet found
         if (leadjets.size() < 2) {
           MSG_DEBUG("Could not find two suitable leading jets");
           continue;
         }
+        _sumW += evt.weight();
+
+        cout << 2 << endl;
 
         const double rap1 = leadjets[0].rapidity();
         const double rap2 = leadjets[1].rapidity();
         const double mass = (leadjets[0] + leadjets[1]).mass();
-        const double ymax = (fabs(rap1) > fabs(rap2)) ? fabs(rap1) : fabs(rap2);
+        //const double ymax = (fabs(rap1) > fabs(rap2)) ? fabs(rap1) : fabs(rap2);
+        const double ymax = max(fabs(rap1), fabs(rap2));
         const double chi = exp(fabs(rap1 - rap2));
-
-        if (fabs(rap1 + rap2) < ySumCut && fabs(rap1 - rap2) < asymmetryCut) {
+        if (fabs(rap1 + rap2) < 2.2) {
           _chiVsMass[alg].fill(mass/GeV, chi, evt.weight());
         }
         _massVsY[alg].fill(ymax, mass/GeV, evt.weight());
@@ -101,17 +114,19 @@
 
 
     void finalize() {
-      const double xs = crossSectionPerEvent() / picobarn;
       for (size_t alg = 0; alg < 2; ++alg) {
-        _pTHistos[alg].scale(xs, this);
-        _massVsY[alg].scale(xs, this);
-        _chiVsMass[alg].scale(xs, this);
+        _pTHistos[alg].scale(crossSectionPerEvent()/picobarn, this);
+        _massVsY[alg].scale(crossSection()/_sumW/picobarn, this);
+        _chiVsMass[alg].scale(crossSection()/_sumW/picobarn, this);
       }
     }
 
 
   private:
 
+    /// Counter for weights passing the dijet requirement cut
+    double _sumW;
+
     /// The inclusive pT spectrum for akt6 and akt4 jets (array index is jet type from enum above)
     BinnedHistogram<double> _pTHistos[2];
 

Modified: trunk/src/Core/AnalysisInfo.cc
==============================================================================
--- trunk/src/Core/AnalysisInfo.cc	Tue Jan 11 14:59:39 2011	(r2875)
+++ trunk/src/Core/AnalysisInfo.cc	Wed Jan 12 22:55:30 2011	(r2876)
@@ -2,6 +2,7 @@
 #include "Rivet/RivetBoost.hh"
 #include "Rivet/AnalysisInfo.hh"
 #include "Rivet/Tools/Utils.hh"
+#include "Rivet/Tools/RivetPaths.hh"
 #include "Rivet/Tools/Logging.hh"
 #include "yaml-cpp/yaml.h"
 #include <iostream>
@@ -12,7 +13,6 @@
 
 
   /// Ideas:
-  ///  * search RIVET_INFO_PATH etc. for <name>.info.yaml
   ///  * how to determine the name?
   ///  * only populate pointer on Analysis when requested
   ///  * use smart pointer: deletes automatically when Analysis
@@ -22,15 +22,9 @@
   /// Static factory method
   AnalysisInfo* AnalysisInfo::make(const std::string& ananame) {
     // Build the list of directories to search
-    vector<string> dirs;
-    char* env = 0;
-    // First try to use the Rivet data path variable
-    env = getenv("RIVET_INFO_PATH");
-    if (env) dirs += pathsplit(env);
-    // Then try to use the Rivet data install path
-    dirs += getRivetDataPath();
+    vector<string> dirs = getAnalysisInfoPaths();
     // And the current dir
-    dirs += ".";
+    // dirs += ".";
 
     bool found = false;
     string datapath = "";

Modified: trunk/src/Tools/RivetPaths.cc
==============================================================================
--- trunk/src/Tools/RivetPaths.cc	Tue Jan 11 14:59:39 2011	(r2875)
+++ trunk/src/Tools/RivetPaths.cc	Wed Jan 12 22:55:30 2011	(r2876)
@@ -92,4 +92,21 @@
   }
 
 
+  const vector<string> getAnalysisPlotPaths() {
+    vector<string> dirs;
+    char* env = 0;
+    env = getenv("RIVET_PLOT_PATH");
+    if (env) {
+      // Use the Rivet analysis path variable if set...
+      dirs += pathsplit(env);
+    } else {
+      // ... otherwise fall back to the Rivet data install path
+      dirs += getRivetDataPath();
+      // ... and also add any analysis plugin search dirs for convenience
+      dirs += getAnalysisLibPaths();
+    }
+    return dirs;
+  }
+
+
 }


More information about the Rivet-svn mailing list