|
[Rivet-svn] r3614 - in branches/2011-07-aida2yoda: . bin data data/anainfo data/plotinfo doc include/Rivet/Projections pyext src/Analyses src/Core src/Projectionsblackhole at projects.hepforge.org blackhole at projects.hepforge.orgFri Mar 16 15:23:52 GMT 2012
Author: hoeth Date: Fri Mar 16 15:23:52 2012 New Revision: 3614 Log: merge c3549-3560 from trunk Modified: branches/2011-07-aida2yoda/ChangeLog branches/2011-07-aida2yoda/bin/compare-histos branches/2011-07-aida2yoda/bin/make-plots branches/2011-07-aida2yoda/bin/rivet branches/2011-07-aida2yoda/bin/rivet-mkhtml branches/2011-07-aida2yoda/data/anainfo/MC_TTBAR.info branches/2011-07-aida2yoda/data/plotinfo/ATLAS_2012_I1083318.plot branches/2011-07-aida2yoda/data/plotinfo/MC_TTBAR.plot branches/2011-07-aida2yoda/data/rivet-completion branches/2011-07-aida2yoda/doc/make-plots.txt branches/2011-07-aida2yoda/include/Rivet/Projections/ChargedLeptons.hh branches/2011-07-aida2yoda/include/Rivet/Projections/UnstableFinalState.hh branches/2011-07-aida2yoda/pyext/lighthisto.py branches/2011-07-aida2yoda/src/Analyses/ALEPH_1996_S3486095.cc branches/2011-07-aida2yoda/src/Analyses/DELPHI_1996_S3430090.cc branches/2011-07-aida2yoda/src/Analyses/MC_TTBAR.cc branches/2011-07-aida2yoda/src/Core/Run.cc branches/2011-07-aida2yoda/src/Projections/ChargedLeptons.cc branches/2011-07-aida2yoda/src/Projections/UnstableFinalState.cc Modified: branches/2011-07-aida2yoda/ChangeLog ============================================================================== --- branches/2011-07-aida2yoda/ChangeLog Fri Mar 16 15:12:00 2012 (r3613) +++ branches/2011-07-aida2yoda/ChangeLog Fri Mar 16 15:23:52 2012 (r3614) @@ -1,3 +1,48 @@ +2012-02-14 Hendrik Hoeth <hendrik.hoeth at cern.ch> + + * DELPHI_1996_S3430090 and ALEPH_1996_S3486095: + fix rapidity vs {Thrust,Sphericity}-axis. + +2012-02-14 Andy Buckley <andy.buckley at cern.ch> + + * bin/compare-histos: Don't attempt to remove bins from MC histos + where they aren't found in the ref file, if the ref file is not + expt data, or if the new --no-rmgapbins arg is given. + + * bin/rivet: Remove the conversion of requested analysis names to + upper-case: mixed-case analysis names will now work. + +2012-02-14 Frank Siegert <frank.siegert at cern.ch> + + * Bugfixes and improvements for MC_TTBAR: + - Avoid assert failure with logspace starting at 0.0 + - Ignore charged lepton in jet finding (otherwise jet multi is always + +1). + - Add some dR/deta/dphi distributions as noted in TODO + - Change pT plots to logspace as well (to avoid low-stat high pT bins) + +2012-02-10 Hendrik Hoeth <hendrik.hoeth at cern.ch> + + * rivet-mkhtml -c option now has the semantics of a .plot + file. The contents are appended to the dat output by + compare-histos. + +2012-02-09 David Grellscheid <david.grellscheid at durham.ac.uk> + + * Fixed broken UnstableFS behaviour + +2012-01-25 Frank Siegert <frank.siegert at cern.ch> + + * Improvements in make-plots: + - Add PlotTickLabels and RatioPlotTickLabels options (cf. + make-plots.txt) + - Make ErrorBars and ErrorBands non-exclusive (and change + their order, such that Bars are on top of Bands) + +2012-01-25 Holger Schulz <holger.schulz at physik.hu-berlin.de> + + * Add ATLAS diffractive gap analysis + 2012-01-23 Andy Buckley <andy.buckley at cern.ch> * bin/rivet: When using --list-analyses, the analysis summary is Modified: branches/2011-07-aida2yoda/bin/compare-histos ============================================================================== --- branches/2011-07-aida2yoda/bin/compare-histos Fri Mar 16 15:12:00 2012 (r3613) +++ branches/2011-07-aida2yoda/bin/compare-histos Fri Mar 16 15:23:52 2012 (r3614) @@ -93,7 +93,7 @@ h.expt = dpsname.split("_")[0][1:] ## Hard-coded cosmetic handling for the D0 experiment name! if h.expt == "D0": - h.expt = "D\O\ " + h.expt = "D\O{}" histos[dpsname] = h return histos, titles, xlabels, ylabels @@ -113,6 +113,9 @@ parser.add_option("--plotinfodir", dest="PLOTINFODIR", action="append", default=["."], help="directory which may contain plot header information (in addition " "to standard Rivet search paths)") + parser.add_option("--no-rmgapbins", dest="RMBINS", action="store_false", + default=True, help="disable attempting to remove 'gap' bins from MC histos when they " + "don't appear in the ref file") stygroup = OptionGroup(parser, "Plot style") stygroup.add_option("--refid", dest="REF_ID", @@ -134,6 +137,8 @@ "(useful when the plot description should only be given in a caption)") stygroup.add_option("--style", dest="STYLE", default="default", help="change plotting style: default|bw|talk") + stygroup.add_option("-c", "--config", dest="CONFIGFILES", action="append", default=["~/.make-plots"], + help="additional plot config file(s). Settings will be included in the output configuration.") parser.add_option_group(stygroup) selgroup = OptionGroup(parser, "Selective plotting") @@ -385,7 +390,7 @@ ## Write out histos num_written = 0 - plotparser = PlotParser(opts.PLOTINFODIR) + plotparser = PlotParser(opts.PLOTINFODIR, opts.CONFIGFILES) for name in sorted(activenames): logging.debug("Writing histos for plot '%s'" % name) @@ -465,7 +470,12 @@ if opts.RATIO_DEVIATION: paramdefaults["RatioPlotMode"] = "deviation" + + ## Behaviour changes if "ref" data is just another MC file if not HISTOS[ref][name].isdata: + ## Don't attempt to remove "gap" bins + opts.RMBINS = False + ## Different ratio plot label (not MC/Data) paramdefaults["RatioPlotYLabel"] = "Ratio" if opts.RATIO_DEVIATION: paramdefaults["RatioPlotYLabel"] = "Deviation" @@ -534,14 +544,17 @@ numskipped = 0 #print hfile, name, HISTOS[hfile][name].numBins(), HISTOS[ref][name].numBins() for ibin, bin in enumerate(HISTOS[hfile][name].getBins()): - xmin, xmax = bin.getXRange() ## Skip writing this MC bin if the bin edges don't match, and the MC histo has too many bins - if hfile != ref and HISTOS[hfile][name].numBins() > HISTOS[ref][name].numBins(): - rxmin, rxmax = HISTOS[ref][name].getBin(ibin-numskipped).getXRange() - if not eq(rxmin, xmin) or not eq(rxmax, xmax): - numskipped += 1 - assert(numskipped <= (HISTOS[hfile][name].numBins() - HISTOS[ref][name].numBins())) - continue + if opts.RMBINS: + xmin, xmax = bin.getXRange() + #print HISTOS[hfile][name].numBins(), HISTOS[ref][name].numBins() + if hfile != ref and HISTOS[hfile][name].numBins() > HISTOS[ref][name].numBins(): + rxmin, rxmax = HISTOS[ref][name].getBin(ibin-numskipped).getXRange() + if not eq(rxmin, xmin) or not eq(rxmax, xmax): + numskipped += 1 + #print numskipped, (HISTOS[hfile][name].numBins() - HISTOS[ref][name].numBins()) + assert(numskipped <= abs(HISTOS[hfile][name].numBins() - HISTOS[ref][name].numBins())) + continue histstr += '%s\n' % (bin.asFlat()) histstr += "# END HISTOGRAM\n" histstrs.append(histstr) Modified: branches/2011-07-aida2yoda/bin/make-plots ============================================================================== --- branches/2011-07-aida2yoda/bin/make-plots Fri Mar 16 15:12:00 2012 (r3613) +++ branches/2011-07-aida2yoda/bin/make-plots Fri Mar 16 15:23:52 2012 (r3614) @@ -532,7 +532,7 @@ for i in range(len(FOO)): xcustomminorticks.append({'Value': float(FOO[i])}) xticks = XTicks(inputdata.description, self.coors) - if inputdata.description.has_key('RatioPlot') and inputdata.description['RatioPlot']=='1': + if (inputdata.description.has_key('RatioPlot') and inputdata.description['RatioPlot']=='1') or (inputdata.description.has_key('PlotTickLabels') and inputdata.description['PlotTickLabels']=='0'): drawlabels=False else: drawlabels=True @@ -738,10 +738,15 @@ for i in range(len(FOO)): xcustomminorticks.append({'Value': float(FOO[i])}) xticks = XTicks(inputdata.description, self.coors) + if inputdata.description.has_key('RatioPlotTickLabels') and inputdata.description['RatioPlotTickLabels']=='0': + drawlabels=False + else: + drawlabels=True out += xticks.draw(custommajortickmarks=xcustommajortickmarks,\ customminortickmarks=xcustomminortickmarks,\ custommajorticks=xcustommajorticks,\ - customminorticks=xcustomminorticks) + customminorticks=xcustomminorticks, + drawlabels=drawlabels) if inputdata.description.has_key('YMajorTickMarks') and inputdata.description['YMajorTickMarks']!='': ycustommajortickmarks=int(inputdata.description['YMajorTickMarks']) @@ -1477,6 +1482,14 @@ + coors.strphys2frameX(self.data[i]['UpEdge'][0]) + ', ' \ + coors.strphys2frameY(self.data[i]['UpEdge'][1]) + ')\n') else: + if self.getErrorBands(): + self.description['SmoothLine']=0 + for i in range(len(self.data)): + out += ('\\psframe[dimen=inner,linewidth=0pt,linestyle=none,fillstyle=solid,fillcolor=%s,opacity=%s]' %(self.getErrorBandColor(),self.getErrorBandOpacity())) + out += ('(' + coors.strphys2frameX(self.data[i]['LowEdge']) + ', ' \ + + coors.strphys2frameY(self.data[i]['Content']-self.data[i]['Error'][0]) + ')(' \ + + coors.strphys2frameX(self.data[i]['UpEdge']) + ', ' \ + + coors.strphys2frameY(self.data[i]['Content']+self.data[i]['Error'][1]) + ')\n') if self.getErrorBars(): for i in range(len(self.data)): if self.data[i]['Content']==0. and self.data[i]['Error']==[0.,0.]: continue @@ -1491,41 +1504,32 @@ + coors.strphys2frameY(self.data[i]['Content']-self.data[i]['Error'][0]) + ')(' \ + bincenter + ', ' \ + coors.strphys2frameY(self.data[i]['Content']+self.data[i]['Error'][1]) + ')\n') + if self.getSmoothLine(): + out += ('\psbezier') else: - if self.getErrorBands(): - self.description['SmoothLine']=0 - for i in range(len(self.data)): - out += ('\\psframe[dimen=inner,linewidth=0pt,linestyle=none,fillstyle=solid,fillcolor=%s,opacity=%s]' %(self.getErrorBandColor(),self.getErrorBandOpacity())) - out += ('(' + coors.strphys2frameX(self.data[i]['LowEdge']) + ', ' \ - + coors.strphys2frameY(self.data[i]['Content']-self.data[i]['Error'][0]) + ')(' \ - + coors.strphys2frameX(self.data[i]['UpEdge']) + ', ' \ - + coors.strphys2frameY(self.data[i]['Content']+self.data[i]['Error'][1]) + ')\n') + out += ('\psline') + if (self.getFillStyle() != 'none'): # make sure that filled areas go all the way down to the x-axis + if (coors.phys2frameX(self.data[0]['LowEdge']) > 1e-4): + out += '(' + coors.strphys2frameX(self.data[0]['LowEdge']) + ', -0.1)\n' + else: + out += '(-0.1, -0.1)\n' + for i in range(len(self.data)): if self.getSmoothLine(): - out += ('\psbezier') + out += ('(' + coors.strphys2frameX(0.5*(self.data[i]['LowEdge']+self.data[i]['UpEdge'])) + ', ' \ + + coors.strphys2frameY(self.data[i]['Content']) + ')\n') else: - out += ('\psline') - if (self.getFillStyle() != 'none'): # make sure that filled areas go all the way down to the x-axis - if (coors.phys2frameX(self.data[0]['LowEdge']) > 1e-4): - out += '(' + coors.strphys2frameX(self.data[0]['LowEdge']) + ', -0.1)\n' - else: - out += '(-0.1, -0.1)\n' - for i in range(len(self.data)): - if self.getSmoothLine(): - out += ('(' + coors.strphys2frameX(0.5*(self.data[i]['LowEdge']+self.data[i]['UpEdge'])) + ', ' \ - + coors.strphys2frameY(self.data[i]['Content']) + ')\n') - else: - out += ('(' + coors.strphys2frameX(self.data[i]['LowEdge']) + ', ' \ - + coors.strphys2frameY(self.data[i]['Content']) + ')(' \ - + coors.strphys2frameX(self.data[i]['UpEdge']) + ', ' \ - + coors.strphys2frameY(self.data[i]['Content']) + ')\n') - if not (self.description.has_key('ConnectGaps') and self.description['ConnectGaps']=='1'): - if (i+1 < len(self.data)) and (abs(coors.phys2frameX(self.data[i]['UpEdge']) - coors.phys2frameX(self.data[i+1]['LowEdge'])) > 1e-4): - out += ('\\psline') - if (self.getFillStyle() != 'none'): # make sure that filled areas go all the way down to the x-axis - if (coors.phys2frameX(self.data[-1]['UpEdge']) < 1-1e-4): - out += '(' + coors.strphys2frameX(self.data[-1]['UpEdge']) + ', -0.1)\n' - else: - out += '(1.1, -0.1)\n' + out += ('(' + coors.strphys2frameX(self.data[i]['LowEdge']) + ', ' \ + + coors.strphys2frameY(self.data[i]['Content']) + ')(' \ + + coors.strphys2frameX(self.data[i]['UpEdge']) + ', ' \ + + coors.strphys2frameY(self.data[i]['Content']) + ')\n') + if not (self.description.has_key('ConnectGaps') and self.description['ConnectGaps']=='1'): + if (i+1 < len(self.data)) and (abs(coors.phys2frameX(self.data[i]['UpEdge']) - coors.phys2frameX(self.data[i+1]['LowEdge'])) > 1e-4): + out += ('\\psline') + if (self.getFillStyle() != 'none'): # make sure that filled areas go all the way down to the x-axis + if (coors.phys2frameX(self.data[-1]['UpEdge']) < 1-1e-4): + out += '(' + coors.strphys2frameX(self.data[-1]['UpEdge']) + ', -0.1)\n' + else: + out += '(1.1, -0.1)\n' # if self.getPolyMarker() != '': for i in range(len(self.data)): Modified: branches/2011-07-aida2yoda/bin/rivet ============================================================================== --- branches/2011-07-aida2yoda/bin/rivet Fri Mar 16 15:12:00 2012 (r3613) +++ branches/2011-07-aida2yoda/bin/rivet Fri Mar 16 15:23:52 2012 (r3614) @@ -414,16 +414,16 @@ # if opts.ALL_ANALYSES: # opts.ANALYSES = all_analyses for a in opts.ANALYSES: - a_up = a.upper() + #a_up = a.upper() ## Print warning message and exit if not a valid analysis name - if not a_up in all_analyses: - logging.warning("'%s' is not a valid analysis. Available analyses are:" % a_up) + if not a in all_analyses: + logging.warning("'%s' is not a valid analysis. Available analyses are:" % a) for aa in all_analyses: logging.warning(" %s" % aa) logging.warning("Exiting...") sys.exit(1) - logging.debug("Adding analysis '%s'" % a_up) - ah.addAnalysis(a_up) + logging.debug("Adding analysis '%s'" % a) + ah.addAnalysis(a) ## Read and process events run = rivet.Run(ah) Modified: branches/2011-07-aida2yoda/bin/rivet-mkhtml ============================================================================== --- branches/2011-07-aida2yoda/bin/rivet-mkhtml Fri Mar 16 15:12:00 2012 (r3613) +++ branches/2011-07-aida2yoda/bin/rivet-mkhtml Fri Mar 16 15:23:52 2012 (r3614) @@ -38,7 +38,7 @@ parser.add_option("-t", "--title", dest="TITLE", default="Plots from Rivet analyses", help="title to be displayed on the main web page") parser.add_option("-c", "--config", dest="CONFIGFILES", action="append", default=["~/.make-plots"], - help="plot config file(s) to be used with make-plots.") + help="plot config file(s) to be used with compare-histos.") parser.add_option("-s", "--single", dest="SINGLE", action="store_true", default=False, help="display plots on single webpage.") parser.add_option("--no-ratio", dest="SHOW_RATIO", action="store_false", @@ -182,7 +182,13 @@ newarg += ":%s" % opt # print newarg ch_cmd.append(newarg) +for configfile in opts.CONFIGFILES: + configfile = os.path.abspath(os.path.expanduser(configfile)) + if os.access(configfile, os.R_OK): + ch_cmd.append("-c") + ch_cmd.append(configfile) # TODO: Pass rivet-mkhtml -m and -M args to compare-histos + if opts.VERBOSE: ch_cmd.append("--verbose") print "Calling compare-histos with the following command:" @@ -313,10 +319,6 @@ elif opts.OUTPUT_FONT == "minion": mp_cmd.append("--minion") mp_cmd.append("--full-range") -for configfile in opts.CONFIGFILES: - 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 analyses: anapath = os.path.join(opts.OUTPUTDIR, analysis) Modified: branches/2011-07-aida2yoda/data/anainfo/MC_TTBAR.info ============================================================================== --- branches/2011-07-aida2yoda/data/anainfo/MC_TTBAR.info Fri Mar 16 15:12:00 2012 (r3613) +++ branches/2011-07-aida2yoda/data/anainfo/MC_TTBAR.info Fri Mar 16 15:23:52 2012 (r3614) @@ -13,5 +13,3 @@ PtCuts: [0] Description: This is a pure Monte Carlo study for semi-leptonic $t\bar{t}$ production. -ToDo: - - Add delta(R,phi,eta) plots between various b and light jets, and w.r.t the top and W candidates Modified: branches/2011-07-aida2yoda/data/plotinfo/ATLAS_2012_I1083318.plot ============================================================================== --- branches/2011-07-aida2yoda/data/plotinfo/ATLAS_2012_I1083318.plot Fri Mar 16 15:12:00 2012 (r3613) +++ branches/2011-07-aida2yoda/data/plotinfo/ATLAS_2012_I1083318.plot Fri Mar 16 15:23:52 2012 (r3614) @@ -206,7 +206,7 @@ # END PLOT # BEGIN PLOT /ATLAS_2012_I1083318/d24-x01-y0[12] -Title=$\Delta y$ Distance of Leading Jets +Title=Rapidity Distance of Leading Jets XLabel=$\Delta y$(First Jet, Second Jet) YLabel=d$\sigma$/d$\Delta y$ [pb] LogY=0 Modified: branches/2011-07-aida2yoda/data/plotinfo/MC_TTBAR.plot ============================================================================== --- branches/2011-07-aida2yoda/data/plotinfo/MC_TTBAR.plot Fri Mar 16 15:12:00 2012 (r3613) +++ branches/2011-07-aida2yoda/data/plotinfo/MC_TTBAR.plot Fri Mar 16 15:23:52 2012 (r3614) @@ -1,72 +1,89 @@ +# BEGIN PLOT /MC_TTBAR/jet.?_[1234]_pT +XLabel=$p_\perp$ [GeV] +YLabel=$1/\sigma \, \mathrm{d}\sigma/\mathrm{d}p_\perp$ [GeV$^{-1}$] +LogX=1 +# END PLOT + # BEGIN PLOT /MC_TTBAR/jet_1_pT Title=Transverse momentum distribution for jet 1 -XLabel=$p_\perp$ [GeV] -YLabel=$1/\sigma \, \mathrm{d}\sigma/\mathrm{d}p_\perp$ [GeV$^-1$] # END PLOT # BEGIN PLOT /MC_TTBAR/jet_2_pT Title=Transverse momentum distribution for jet 2 -XLabel=$p_\perp$ [GeV] -YLabel=$1/\sigma \, \mathrm{d}\sigma/\mathrm{d}p_\perp$ [GeV$^-1$] # END PLOT # BEGIN PLOT /MC_TTBAR/jet_3_pT Title=Transverse momentum distribution for jet 3 -XLabel=$p_\perp$ [GeV] -YLabel=$1/\sigma \, \mathrm{d}\sigma/\mathrm{d}p_\perp$ [GeV$^-1$] # END PLOT # BEGIN PLOT /MC_TTBAR/jet_4_pT Title=Transverse momentum distribution for jet 4 -XLabel=$p_\perp$ [GeV] -YLabel=$1/\sigma \, \mathrm{d}\sigma/\mathrm{d}p_\perp$ [GeV$^-1$] # END PLOT # BEGIN PLOT /MC_TTBAR/jet_HT Title=$H_T$ distribution for all jets XLabel=$H_T$ [GeV] -YLabel=$1/\sigma \, \mathrm{d}\sigma/\mathrm{d}H_T$ [GeV$^-1$] +YLabel=$1/\sigma \, \mathrm{d}\sigma/\mathrm{d}H_T$ [GeV$^{-1}$] LogX=1 # END PLOT # BEGIN PLOT /MC_TTBAR/jetb_1_pT Title=Transverse momentum distribution for $b$-jet 1 -XLabel=$p_\perp$ [GeV] -YLabel=$1/\sigma \, \mathrm{d}\sigma/\mathrm{d}p_\perp$ [GeV$^-1$] # END PLOT # BEGIN PLOT /MC_TTBAR/jetb_2_pT Title=Transverse momentum distribution for $b$-jet 2 -XLabel=$p_\perp$ [GeV] -YLabel=$1/\sigma \, \mathrm{d}\sigma/\mathrm{d}p_\perp$ [GeV$^-1$] # END PLOT -# + # BEGIN PLOT /MC_TTBAR/jetl_1_pT Title=Transverse momentum distribution for light jet 1 -XLabel=$p_\perp$ [GeV] -YLabel=$1/\sigma \, \mathrm{d}\sigma/\mathrm{d}p_\perp$ [GeV$^-1$] # END PLOT # BEGIN PLOT /MC_TTBAR/jetl_2_pT Title=Transverse momentum distribution for light jet 2 -XLabel=$p_\perp$ [GeV] -YLabel=$1/\sigma \, \mathrm{d}\sigma/\mathrm{d}p_\perp$ [GeV$^-1$] # END PLOT # BEGIN PLOT /MC_TTBAR/W_mass Title=Mass distribution for $W$ bosons XLabel=$m_{jj}$ [GeV] -YLabel=$1/\sigma \, \mathrm{d}\sigma/\mathrm{d}m_{jj}$ [GeV$^-1$] +YLabel=$1/\sigma \, \mathrm{d}\sigma/\mathrm{d}m_{jj}$ [GeV$^{-1}$] +LogY=0 # END PLOT # BEGIN PLOT /MC_TTBAR/t_mass Title=Mass distribution for reconstructed top -XLabel=$m_{q\bar{q}b$ [GeV] -YLabel=$1/\sigma \, \mathrm{d}\sigma/\mathrm{d}m_{q\bar{q}b}$ [GeV$^-1$] +XLabel=$m_{q\bar{q}b}$ [GeV] +YLabel=$1/\sigma \, \mathrm{d}\sigma/\mathrm{d}m_{q\bar{q}b}$ [GeV$^{-1}$] +LogY=0 # END PLOT # BEGIN PLOT /MC_TTBAR/t_mass_W_cut Title=Mass distribution for reconstructed top after $m_W$ cut -XLabel=$m_{q\bar{q}b$ [GeV] -YLabel=$1/\sigma \, \mathrm{d}\sigma/\mathrm{d}m_{q\bar{q}b}$ [GeV$^-1$] +XLabel=$m_{q\bar{q}b}$ [GeV] +YLabel=$1/\sigma \, \mathrm{d}\sigma/\mathrm{d}m_{q\bar{q}b}$ [GeV$^{-1}$] +LogY=0 +# END PLOT + +# BEGIN PLOT /MC_TTBAR/.*_mass +XLabel=$m$ [GeV] +YLabel=$1/\sigma \, \mathrm{d}\sigma/\mathrm{d}m$ [GeV$^{-1}$] +# END PLOT + +# BEGIN PLOT /MC_TTBAR/.*_dR +XLabel=$\Delta R$ +YLabel=$1/\sigma \, \mathrm{d}\sigma/\mathrm{d}\Delta R$ +LogY=0 # END PLOT + +# BEGIN PLOT /MC_TTBAR/.*_deta +XLabel=$\Delta \eta$ +YLabel=$1/\sigma \, \mathrm{d}\sigma/\mathrm{d}\Delta \eta$ +LogY=0 +# END PLOT + +# BEGIN PLOT /MC_TTBAR/.*_dphi +XLabel=$\Delta \phi$ +YLabel=$1/\sigma \, \mathrm{d}\sigma/\mathrm{d}\Delta \phi$ +LogY=0 +# END PLOT + Modified: branches/2011-07-aida2yoda/data/rivet-completion ============================================================================== --- branches/2011-07-aida2yoda/data/rivet-completion Fri Mar 16 15:12:00 2012 (r3613) +++ branches/2011-07-aida2yoda/data/rivet-completion Fri Mar 16 15:23:52 2012 (r3614) @@ -109,6 +109,7 @@ opts="$opts --all --show-mc-only --show-single --refid" opts="$opts --no-plottitle" opts="$opts --plotinfodir" + opts="$opts --no-rmgapbins" opts="$opts --quiet -q --verbose -v" if [[ ${cur} == -* ]] ; then COMPREPLY=( $(compgen -W "${opts}" -- ${cur}) ) Modified: branches/2011-07-aida2yoda/doc/make-plots.txt ============================================================================== --- branches/2011-07-aida2yoda/doc/make-plots.txt Fri Mar 16 15:12:00 2012 (r3613) +++ branches/2011-07-aida2yoda/doc/make-plots.txt Fri Mar 16 15:23:52 2012 (r3614) @@ -143,6 +143,13 @@ To specify minor ticks at arbitrary positions. `<list>` is a tab separated list of format `value1 <tab> value2 <tab> value3 ...`. +-------------------- +PlotTickLabels=<0|1> +RatioPlotTickLabels=<0|1> +-------------------- +Disable/enable plotting of the tick labels in the plot and ratio plot (useful +if multiple plots are to be combined manually later). + Axes ^^^^ Modified: branches/2011-07-aida2yoda/include/Rivet/Projections/ChargedLeptons.hh ============================================================================== --- branches/2011-07-aida2yoda/include/Rivet/Projections/ChargedLeptons.hh Fri Mar 16 15:12:00 2012 (r3613) +++ branches/2011-07-aida2yoda/include/Rivet/Projections/ChargedLeptons.hh Fri Mar 16 15:23:52 2012 (r3614) @@ -13,7 +13,7 @@ /// @brief Get charged final-state leptons /// /// NB. This is just electrons and muons, unless you set taus stable! - class ChargedLeptons : public Projection { + class ChargedLeptons : public FinalState { public: /// Constructor @@ -40,14 +40,9 @@ /// Access the projected leptons. const ParticleVector& chargedLeptons() const { - return _theChargedLeptons; + return _theParticles; } - private: - - /// The leptons - ParticleVector _theChargedLeptons; - }; Modified: branches/2011-07-aida2yoda/include/Rivet/Projections/UnstableFinalState.hh ============================================================================== --- branches/2011-07-aida2yoda/include/Rivet/Projections/UnstableFinalState.hh Fri Mar 16 15:12:00 2012 (r3613) +++ branches/2011-07-aida2yoda/include/Rivet/Projections/UnstableFinalState.hh Fri Mar 16 15:23:52 2012 (r3614) @@ -32,12 +32,9 @@ UnstableFinalState(double mineta = -MAXRAPIDITY, double maxeta = MAXRAPIDITY, double minpt = 0.0*GeV) - : _etamin(mineta), _etamax(maxeta), _ptmin(minpt) + : FinalState(mineta,maxeta,minpt) { setName("UnstableFinalState"); - // addCut("eta", MORE_EQ, mineta); - // addCut("eta", LESS_EQ, maxeta); - // addCut("pT", MORE_EQ, minpt); } @@ -48,45 +45,11 @@ //@} - - /// @name Accessors - //@{ - - /// Access the projected final-state particles. - virtual const ParticleVector& particles() const { return _theParticles; } - - /// Is this final state empty? - virtual bool empty() const { return _theParticles.empty(); } - - /// @deprecated Is this final state empty? - virtual bool isEmpty() const { return _theParticles.empty(); } - - //@} - - protected: /// Apply the projection to the event. virtual void project(const Event& e); - /// Compare projections. - virtual int compare(const Projection& p) const; - - - protected: - - /// The minimum allowed pseudorapidity. - double _etamin; - - /// The maximum allowed pseudorapidity. - double _etamax; - - /// The minimum allowed transverse momentum. - double _ptmin; - - /// The final-state particles. - ParticleVector _theParticles; - }; Modified: branches/2011-07-aida2yoda/pyext/lighthisto.py ============================================================================== --- branches/2011-07-aida2yoda/pyext/lighthisto.py Fri Mar 16 15:12:00 2012 (r3613) +++ branches/2011-07-aida2yoda/pyext/lighthisto.py Fri Mar 16 15:23:52 2012 (r3614) @@ -554,7 +554,7 @@ pat_property = re.compile('^(\w+?)=(.*)$') pat_path_property = re.compile('^(\S+?)::(\w+?)=(.*)$') - def __init__(self, plotpaths=None): + def __init__(self, plotpaths=None, addfiles=[]): """ Parameters ---------- @@ -573,6 +573,8 @@ plotpaths = [] self.plotpaths = plotpaths + self.addfiles = addfiles + if len(self.plotpaths) == 0: try: import rivet @@ -615,39 +617,43 @@ ret = {'PLOT': {}, 'SPECIAL': None, 'HISTOGRAM': {}} for pidir in self.plotpaths: plotfile = os.path.join(pidir, base) - if os.access(plotfile, os.R_OK): - #print plotfile - startreading = False - f = open(plotfile) - for line in f: - m = self.pat_begin_block.match(line) - if m: - tag, pathpat = m.group(1,2) - # pathpat could be a regex - if tag == section and re.match(pathpat,hpath): - startreading = True - if section in ['SPECIAL']: - ret[section] = '' - continue - if not startreading: - continue - if self.isEndMarker(line, section): - startreading = False - continue - elif self.isComment(line): - continue - if section in ['PLOT', 'HISTOGRAM']: - vm = self.pat_property.match(line) - if vm: - prop, value = vm.group(1,2) - #print prop, value - ret[section][prop] = value - elif section in ['SPECIAL']: - ret[section] += line - f.close() - # no break, as we can collect settings from multiple .plot files + self.readHeadersFromFile(plotfile, ret, section, hpath) + # no break, as we can collect settings from multiple .plot files + for extrafile in self.addfiles: + self.readHeadersFromFile(extrafile, ret, section, hpath) return ret[section] + def readHeadersFromFile(self, plotfile, ret, section, hpath): + """Get a section for a histogram from a .plot file.""" + if os.access(plotfile, os.R_OK): + startreading = False + f = open(plotfile) + for line in f: + m = self.pat_begin_block.match(line) + if m: + tag, pathpat = m.group(1,2) + # pathpat could be a regex + if tag == section and re.match(pathpat,hpath): + startreading = True + if section in ['SPECIAL']: + ret[section] = '' + continue + if not startreading: + continue + if self.isEndMarker(line, section): + startreading = False + continue + elif self.isComment(line): + continue + if section in ['PLOT', 'HISTOGRAM']: + vm = self.pat_property.match(line) + if vm: + prop, value = vm.group(1,2) + #print prop, value + ret[section][prop] = value + elif section in ['SPECIAL']: + ret[section] += line + f.close() def getHeaders(self, hpath): """Get the plot headers for histogram hpath. Modified: branches/2011-07-aida2yoda/src/Analyses/ALEPH_1996_S3486095.cc ============================================================================== --- branches/2011-07-aida2yoda/src/Analyses/ALEPH_1996_S3486095.cc Fri Mar 16 15:12:00 2012 (r3613) +++ branches/2011-07-aida2yoda/src/Analyses/ALEPH_1996_S3486095.cc Fri Mar 16 15:23:52 2012 (r3614) @@ -215,7 +215,7 @@ // Calculate rapidities w.r.t. thrust. const double rapidityT = 0.5 * std::log((energy + momT) / (energy - momT)); - _histRapidityT->fill(rapidityT, weight); + _histRapidityT->fill(fabs(rapidityT), weight); if (std::fabs(rapidityT) <= 0.5) { rapt05 += 1.0; } Modified: branches/2011-07-aida2yoda/src/Analyses/DELPHI_1996_S3430090.cc ============================================================================== --- branches/2011-07-aida2yoda/src/Analyses/DELPHI_1996_S3430090.cc Fri Mar 16 15:12:00 2012 (r3613) +++ branches/2011-07-aida2yoda/src/Analyses/DELPHI_1996_S3430090.cc Fri Mar 16 15:23:52 2012 (r3614) @@ -251,8 +251,8 @@ // Calculate rapidities w.r.t. thrust and sphericity. const double rapidityT = 0.5 * std::log((energy + momT) / (energy - momT)); const double rapidityS = 0.5 * std::log((energy + momS) / (energy - momS)); - _histRapidityT->fill(rapidityT, weight); - _histRapidityS->fill(rapidityS, weight); + _histRapidityT->fill(fabs(rapidityT), weight); + _histRapidityS->fill(fabs(rapidityS), weight); MSG_TRACE(fabs(rapidityT) << " " << scaledMom/GeV); } Evis2 = Evis*Evis; Modified: branches/2011-07-aida2yoda/src/Analyses/MC_TTBAR.cc ============================================================================== --- branches/2011-07-aida2yoda/src/Analyses/MC_TTBAR.cc Fri Mar 16 15:12:00 2012 (r3613) +++ branches/2011-07-aida2yoda/src/Analyses/MC_TTBAR.cc Fri Mar 16 15:23:52 2012 (r3614) @@ -1,5 +1,6 @@ #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" +#include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/ChargedLeptons.hh" #include "Rivet/Projections/MissingMomentum.hh" #include "Rivet/Projections/FastJets.hh" @@ -15,9 +16,6 @@ /// Minimal constructor MC_TTBAR() : Analysis("MC_TTBAR") { - _sumwPassedLepJetMET = 0; - _sumwPassedJetID = 0; - _sumwPassedWMass = 0; } @@ -30,32 +28,51 @@ // A FinalState is used to select particles within |eta| < 4.2 and with pT // > 30 GeV, out of which the ChargedLeptons projection picks only the // electrons and muons, to be accessed later as "LFS". - addProjection(ChargedLeptons(FinalState(-4.2, 4.2, 30*GeV)), "LFS"); + ChargedLeptons lfs(FinalState(-4.2, 4.2, 30*GeV)); + addProjection(lfs, "LFS"); // A second FinalState is used to select all particles in |eta| < 4.2, // with no pT cut. This is used to construct jets and measure missing // transverse energy. - FinalState fs(-4.2, 4.2, 0*GeV); + VetoedFinalState fs(FinalState(-4.2, 4.2, 0*GeV)); + fs.addVetoOnThisFinalState(lfs); addProjection(FastJets(fs, FastJets::ANTIKT, 0.6), "Jets"); addProjection(MissingMomentum(fs), "MissingET"); // Booking of histograms _h_njets = bookHisto1D("jet_mult", 11, -0.5, 10.5); // - _h_jet_1_pT = bookHisto1D("jet_1_pT", 50, 0, 500); - _h_jet_2_pT = bookHisto1D("jet_2_pT", 50, 0, 400); - _h_jet_3_pT = bookHisto1D("jet_3_pT", 50, 0, 300); - _h_jet_4_pT = bookHisto1D("jet_4_pT", 50, 0, 200); - _h_jet_HT = bookHisto1D("jet_HT", logspace(0, 2000, 50)); + _h_jet_1_pT = bookHisto1D("jet_1_pT", logspace(20.0, 500.0, 50)); + _h_jet_2_pT = bookHisto1D("jet_2_pT", logspace(20.0, 400.0, 50)); + _h_jet_3_pT = bookHisto1D("jet_3_pT", logspace(20.0, 300.0, 50)); + _h_jet_4_pT = bookHisto1D("jet_4_pT", logspace(20.0, 200.0, 50)); + _h_jet_HT = bookHisto1D("jet_HT", logspace(100.0, 2000.0, 50)); // - _h_bjet_1_pT = bookHisto1D("jetb_1_pT", 50, 0, 400); - _h_bjet_2_pT = bookHisto1D("jetb_2_pT", 50, 0, 300); + _h_bjet_1_pT = bookHisto1D("jetb_1_pT", logspace(20.0, 400.0, 50)); + _h_bjet_2_pT = bookHisto1D("jetb_2_pT", logspace(20.0, 300.0, 50)); // - _h_ljet_1_pT = bookHisto1D("jetl_1_pT", 50, 0, 400); - _h_ljet_2_pT = bookHisto1D("jetl_2_pT", 50, 0, 300); + _h_ljet_1_pT = bookHisto1D("jetl_1_pT", logspace(20.0, 400.0, 50)); + _h_ljet_2_pT = bookHisto1D("jetl_2_pT", logspace(20.0, 300.0, 50)); // _h_W_mass = bookHisto1D("W_mass", 75, 30, 180); _h_t_mass = bookHisto1D("t_mass", 150, 130, 430); _h_t_mass_W_cut = bookHisto1D("t_mass_W_cut", 150, 130, 430); + // + _h_jetb_1_jetb_2_dR = bookHisto1D("jetb_1_jetb_2_dR", 20, 0.0, 7.0); + _h_jetb_1_jetb_2_deta = bookHisto1D("jetb_1_jetb_2_deta", 20, 0.0, 7.0); + _h_jetb_1_jetb_2_dphi = bookHisto1D("jetb_1_jetb_2_dphi", 20, 0.0, M_PI); + _h_jetb_1_jetl_1_dR = bookHisto1D("jetb_1_jetl_1_dR", 20, 0.0, 7.0); + _h_jetb_1_jetl_1_deta = bookHisto1D("jetb_1_jetl_1_deta", 20, 0.0, 7.0); + _h_jetb_1_jetl_1_dphi = bookHisto1D("jetb_1_jetl_1_dphi", 20, 0.0, M_PI); + _h_jetl_1_jetl_2_dR = bookHisto1D("jetl_1_jetl_2_dR", 20, 0.0, 7.0); + _h_jetl_1_jetl_2_deta = bookHisto1D("jetl_1_jetl_2_deta", 20, 0.0, 7.0); + _h_jetl_1_jetl_2_dphi = bookHisto1D("jetl_1_jetl_2_dphi", 20, 0.0, M_PI); + _h_jetb_1_W_dR = bookHisto1D("jetb_1_W_dR", 20, 0.0, 7.0); + _h_jetb_1_W_deta = bookHisto1D("jetb_1_W_deta", 20, 0.0, 7.0); + _h_jetb_1_W_dphi = bookHisto1D("jetb_1_W_dphi", 20, 0.0, M_PI); + _h_jetb_1_l_dR = bookHisto1D("jetb_1_l_dR", 20, 0.0, 7.0); + _h_jetb_1_l_deta = bookHisto1D("jetb_1_l_deta", 20, 0.0, 7.0); + _h_jetb_1_l_dphi = bookHisto1D("jetb_1_l_dphi", 20, 0.0, M_PI); + _h_jetb_1_l_mass = bookHisto1D("jetb_1_l_mass", 40, 0.0, 500.0); } @@ -96,7 +113,6 @@ } // Update passed-cuts counter and fill all-jets histograms - _sumwPassedLepJetMET += weight; _h_jet_1_pT->fill(alljets[0].momentum().pT()/GeV, weight); _h_jet_2_pT->fill(alljets[1].momentum().pT()/GeV, weight); _h_jet_3_pT->fill(alljets[2].momentum().pT()/GeV, weight); @@ -142,6 +158,7 @@ } } MSG_DEBUG("Number of b-jets = " << bjets.size()); + MSG_DEBUG("Number of l-jets = " << ljets.size()); if (bjets.size() != 2) { MSG_DEBUG("Event failed post-lepton-isolation b-tagging cut"); vetoEvent; @@ -152,7 +169,6 @@ } // Plot the pTs of the identified jets. - _sumwPassedJetID += weight; _h_bjet_1_pT->fill(bjets[0].momentum().pT(), weight); _h_bjet_2_pT->fill(bjets[1].momentum().pT(), weight); _h_ljet_1_pT->fill(ljets[0].momentum().pT(), weight); @@ -186,30 +202,67 @@ _h_t_mass->fill(t2.mass(), weight); // Placing a cut on the well-known W mass helps to reduce backgrounds - if (inRange(W.mass()/GeV, 75, 85)) { + if (inRange(W.mass()/GeV, 75.0, 85.0)) { MSG_DEBUG("W found with mass " << W.mass()/GeV << " GeV"); - _sumwPassedWMass += weight; _h_t_mass_W_cut->fill(t1.mass(), weight); _h_t_mass_W_cut->fill(t2.mass(), weight); + + _h_jetb_1_jetb_2_dR->fill(deltaR(bjets[0].momentum(), bjets[1].momentum()),weight); + _h_jetb_1_jetb_2_deta->fill(fabs(bjets[0].momentum().eta()-bjets[1].momentum().eta()),weight); + _h_jetb_1_jetb_2_dphi->fill(deltaPhi(bjets[0].momentum(),bjets[1].momentum()),weight); + + _h_jetb_1_jetl_1_dR->fill(deltaR(bjets[0].momentum(), ljets[0].momentum()),weight); + _h_jetb_1_jetl_1_deta->fill(fabs(bjets[0].momentum().eta()-ljets[0].momentum().eta()),weight); + _h_jetb_1_jetl_1_dphi->fill(deltaPhi(bjets[0].momentum(),ljets[0].momentum()),weight); + + _h_jetl_1_jetl_2_dR->fill(deltaR(ljets[0].momentum(), ljets[1].momentum()),weight); + _h_jetl_1_jetl_2_deta->fill(fabs(ljets[0].momentum().eta()-ljets[1].momentum().eta()),weight); + _h_jetl_1_jetl_2_dphi->fill(deltaPhi(ljets[0].momentum(),ljets[1].momentum()),weight); + + _h_jetb_1_W_dR->fill(deltaR(bjets[0].momentum(), W),weight); + _h_jetb_1_W_deta->fill(fabs(bjets[0].momentum().eta()-W.eta()),weight); + _h_jetb_1_W_dphi->fill(deltaPhi(bjets[0].momentum(),W),weight); + + FourMomentum l=lfs.chargedLeptons()[0].momentum(); + _h_jetb_1_l_dR->fill(deltaR(bjets[0].momentum(), l),weight); + _h_jetb_1_l_deta->fill(fabs(bjets[0].momentum().eta()-l.eta()),weight); + _h_jetb_1_l_dphi->fill(deltaPhi(bjets[0].momentum(),l),weight); + _h_jetb_1_l_mass->fill(FourMomentum(bjets[0].momentum()+l).mass(), weight); } } void finalize() { - scale(_h_njets, 1/_sumwPassedLepJetMET); - scale(_h_jet_1_pT, 1/_sumwPassedLepJetMET); - scale(_h_jet_2_pT, 1/_sumwPassedLepJetMET); - scale(_h_jet_3_pT, 1/_sumwPassedLepJetMET); - scale(_h_jet_4_pT, 1/_sumwPassedLepJetMET); - scale(_h_jet_HT, 1/_sumwPassedLepJetMET); - scale(_h_bjet_1_pT, 1/_sumwPassedJetID); - scale(_h_bjet_2_pT, 1/_sumwPassedJetID); - scale(_h_ljet_1_pT, 1/_sumwPassedJetID); - scale(_h_ljet_2_pT, 1/_sumwPassedJetID); - scale(_h_W_mass, 1/_sumwPassedJetID); - scale(_h_t_mass, 1/_sumwPassedJetID); - scale(_h_t_mass_W_cut, 1/_sumwPassedWMass); + normalize(_h_njets); + normalize(_h_jet_1_pT); + normalize(_h_jet_2_pT); + normalize(_h_jet_3_pT); + normalize(_h_jet_4_pT); + normalize(_h_jet_HT); + normalize(_h_bjet_1_pT); + normalize(_h_bjet_2_pT); + normalize(_h_ljet_1_pT); + normalize(_h_ljet_2_pT); + normalize(_h_W_mass); + normalize(_h_t_mass); + normalize(_h_t_mass_W_cut); + normalize(_h_jetb_1_jetb_2_dR); + normalize(_h_jetb_1_jetb_2_deta); + normalize(_h_jetb_1_jetb_2_dphi); + normalize(_h_jetb_1_jetl_1_dR); + normalize(_h_jetb_1_jetl_1_deta); + normalize(_h_jetb_1_jetl_1_dphi); + normalize(_h_jetl_1_jetl_2_dR); + normalize(_h_jetl_1_jetl_2_deta); + normalize(_h_jetl_1_jetl_2_dphi); + normalize(_h_jetb_1_W_dR); + normalize(_h_jetb_1_W_deta); + normalize(_h_jetb_1_W_dphi); + normalize(_h_jetb_1_l_dR); + normalize(_h_jetb_1_l_deta); + normalize(_h_jetb_1_l_dphi); + normalize(_h_jetb_1_l_mass); } //@} @@ -217,9 +270,6 @@ private: - // Passed-cuts counters - double _sumwPassedLepJetMET, _sumwPassedJetID, _sumwPassedWMass; - // @name Histogram data members //@{ @@ -230,6 +280,12 @@ Histo1DPtr _h_ljet_1_pT, _h_ljet_2_pT; Histo1DPtr _h_W_mass; Histo1DPtr _h_t_mass, _h_t_mass_W_cut; + Histo1DPtr _h_jetb_1_jetb_2_dR, _h_jetb_1_jetb_2_deta, _h_jetb_1_jetb_2_dphi; + Histo1DPtr _h_jetb_1_jetl_1_dR, _h_jetb_1_jetl_1_deta, _h_jetb_1_jetl_1_dphi; + Histo1DPtr _h_jetl_1_jetl_2_dR, _h_jetl_1_jetl_2_deta, _h_jetl_1_jetl_2_dphi; + Histo1DPtr _h_jetb_1_W_dR, _h_jetb_1_W_deta, _h_jetb_1_W_dphi; + Histo1DPtr _h_jetb_1_l_dR, _h_jetb_1_l_deta, _h_jetb_1_l_dphi,_h_jetb_1_l_mass; + //@} Modified: branches/2011-07-aida2yoda/src/Core/Run.cc ============================================================================== --- branches/2011-07-aida2yoda/src/Core/Run.cc Fri Mar 16 15:12:00 2012 (r3613) +++ branches/2011-07-aida2yoda/src/Core/Run.cc Fri Mar 16 15:23:52 2012 (r3614) @@ -43,7 +43,7 @@ } // Rescale event weights by file-level weight, if scaling is non-trivial if (!fuzzyEquals(_fileweight, 1.0)) { - for (size_t i = 0; i < _evt->weights().size(); ++i) { + for (size_t i = 0; i < (size_t) _evt->weights().size(); ++i) { _evt->weights()[i] *= _fileweight; } } Modified: branches/2011-07-aida2yoda/src/Projections/ChargedLeptons.cc ============================================================================== --- branches/2011-07-aida2yoda/src/Projections/ChargedLeptons.cc Fri Mar 16 15:12:00 2012 (r3613) +++ branches/2011-07-aida2yoda/src/Projections/ChargedLeptons.cc Fri Mar 16 15:23:52 2012 (r3614) @@ -14,16 +14,16 @@ void ChargedLeptons::project(const Event& evt) { // Reset result - _theChargedLeptons.clear(); + _theParticles.clear(); // Loop over charged particles and fill vector with leptons const FinalState& fs = applyProjection<FinalState>(evt, "ChFS"); foreach (const Particle& p, fs.particles()) { if (PID::isLepton(p.pdgId())) { - _theChargedLeptons += Particle(p); + _theParticles += Particle(p); } } - std::sort(_theChargedLeptons.begin(), _theChargedLeptons.end(), cmpParticleByPt); + std::sort(_theParticles.begin(), _theParticles.end(), cmpParticleByPt); } Modified: branches/2011-07-aida2yoda/src/Projections/UnstableFinalState.cc ============================================================================== --- branches/2011-07-aida2yoda/src/Projections/UnstableFinalState.cc Fri Mar 16 15:12:00 2012 (r3613) +++ branches/2011-07-aida2yoda/src/Projections/UnstableFinalState.cc Fri Mar 16 15:23:52 2012 (r3614) @@ -8,25 +8,28 @@ namespace Rivet { - - int UnstableFinalState::compare(const Projection& p) const { - const UnstableFinalState& other = dynamic_cast<const UnstableFinalState&>(p); - return \ - cmp(_etamin, other._etamin) || - cmp(_etamax, other._etamax) || - cmp(_ptmin, other._ptmin); - } - - void UnstableFinalState::project(const Event& e) { _theParticles.clear(); + // \todo We should really implement the FinalState algorithm here instead + double etamin, etamax; + if ( _etaRanges.empty() ) { + etamin = -MAXRAPIDITY; + etamax = MAXRAPIDITY; + } + else { + // with our current constructor choice, we can only ever have one entry + assert( _etaRanges.size() == 1 ); + etamin = _etaRanges[0].first; + etamax = _etaRanges[0].second; + } + foreach (GenParticle* p, Rivet::particles(e.genEvent())) { const int st = p->status(); bool passed = \ ( st == 1 || (st == 2 && abs(p->pdg_id()) != 22) ) && !isZero(p->momentum().perp()) && p->momentum().perp() >= _ptmin && - p->momentum().eta() > _etamin && p->momentum().eta() < _etamax && + p->momentum().eta() > etamin && p->momentum().eta() < etamax && !IS_PARTON_PDGID(p->pdg_id()); // Avoid double counting by re-marking as unpassed if particle ID == parent ID @@ -68,7 +71,7 @@ } } } - MSG_DEBUG("Number of final-state particles = " << _theParticles.size()); + MSG_DEBUG("Number of unstable final-state particles = " << _theParticles.size()); }
More information about the Rivet-svn mailing list |