diff --git a/bin/rivet-cmphistos b/bin/rivet-cmphistos --- a/bin/rivet-cmphistos +++ b/bin/rivet-cmphistos @@ -1,491 +1,495 @@ #! /usr/bin/env python """\ Generate histogram comparison plots USAGE: %(prog)s [options] yodafile1[:'PlotOption1=Value':'PlotOption2=Value':...] [path/to/yodafile2 ...] [PLOT:Key1=Val1:...] where the plot options are described in the make-plots manual in the HISTOGRAM section. ENVIRONMENT: * RIVET_ANALYSIS_PATH: list of paths to be searched for plugin analysis libraries at runtime * RIVET_DATA_PATH: list of paths to be searched for data files """ from __future__ import print_function import rivet, yoda, sys, os rivet.util.check_python_version() rivet.util.set_process_name(os.path.basename(__file__)) class Plot(dict): "A tiny Plot object to help writing out the head in the .dat file" def __repr__(self): return "# BEGIN PLOT\n" + "\n".join("%s=%s" % (k,v) for k,v in self.items()) + "\n# END PLOT\n\n" def sanitiseString(s): #s = s.replace('_','\\_') #s = s.replace('^','\\^{}') #s = s.replace('$','\\$') s = s.replace('#','\\#') s = s.replace('%','\\%') return s def getHistos(filelist): """Loop over all input files. Only use the first occurrence of any REF-histogram and the first occurrence in each MC file for every MC-histogram.""" refhistos, mchistos = {}, {} for infile in filelist: mchistos.setdefault(infile, {}) analysisobjects = yoda.read(infile, patterns=args.PATHPATTERNS, unpatterns=args.PATHUNPATTERNS) #print(analysisobjects) for path, ao in analysisobjects.items(): ## We can't plot non-histograms yet # TODO: support counter plotting with a faked x (or y) position and forced plot width/height if ao.type not in ("Counter", "Histo1D", "Histo2D", "Profile1D", "Profile2D", "Scatter1D", "Scatter2D", "Scatter3D"): continue ## Make a path object and ensure the path is in standard form try: aop = rivet.AOPath(path) except Exception as e: #print(e) print("Found analysis object with non-standard path structure:", path, "... skipping") continue ## We don't plot data objects with path components hidden by an underscore prefix if aop.istmp() or aop.israw(): continue ## Add it to the ref or mc paths, if this path isn't already known basepath = aop.basepath(keepref=False) if aop.isref() and basepath not in refhistos: ao.path = aop.varpath(keepref=False, defaultvarid=0) refhistos[basepath] = ao else: #if basepath not in mchistos[infile]: mchistos[infile].setdefault(basepath, {})[aop.varid(0)] = ao return refhistos, mchistos def getRivetRefData(anas=None): "Find all Rivet reference data files" refhistos = {} rivet_data_dirs = rivet.getAnalysisRefPaths() dirlist = [] for d in rivet_data_dirs: if anas is None: import glob dirlist.append(glob.glob(os.path.join(d, '*.yoda'))) else: dirlist.append([os.path.join(d, a+'.yoda') for a in anas]) for filelist in dirlist: # TODO: delegate to getHistos? for infile in filelist: analysisobjects = yoda.read(infile, patterns=args.PATHPATTERNS, unpatterns=args.PATHUNPATTERNS) for path, ao in analysisobjects.items(): aop = rivet.AOPath(ao.path) if aop.isref(): ao.path = aop.basepath(keepref=False) refhistos[ao.path] = ao return refhistos def parseArgs(args): """Look at the argument list and split it at colons, in order to separate the file names from the plotting options. Store the file names and file specific plotting options.""" filelist = [] plotoptions = {} for a in args: asplit = a.split(':') path = asplit[0] filelist.append(path) plotoptions[path] = [] has_title = False for i in range(1, len(asplit)): ## Add 'Title' if there is no = sign before math mode if '=' not in asplit[i] or ('$' in asplit[i] and asplit[i].index('$') < asplit[i].index('=')): asplit[i] = 'Title=%s' % asplit[i] if asplit[i].startswith('Title='): has_title = True plotoptions[path].append(asplit[i]) if path != "PLOT" and not has_title: plotoptions[path].append('Title=%s' % sanitiseString(os.path.basename( os.path.splitext(path)[0] )) ) return filelist, plotoptions def setStyle(ao, istyle, variation=False): """Set default plot styles (color and line width) colors borrowed from Google Ngrams""" # LINECOLORS = ['{[HTML]{EE3311}}', # red (Google uses 'DC3912') # '{[HTML]{3366FF}}', # blue # '{[HTML]{109618}}', # green # '{[HTML]{FF9900}}', # orange # '{[HTML]{990099}}'] # lilac LINECOLORS = ['red', 'blue', 'green', 'orange', 'lilac'] LINESTYLES = ['solid', 'dashed', 'dashdotted', 'dotted'] if args.STYLE == 'talk': ao.setAnnotation('LineWidth', '1pt') if args.STYLE == 'bw': LINECOLORS = ['black!90', 'black!50', 'black!30'] jc = istyle % len(LINECOLORS) c = LINECOLORS[jc] js = (istyle // len(LINECOLORS)) % len(LINESTYLES) s = LINESTYLES[js] ## If plotting a variation (i.e. band), fade the colour if variation: c += "!30" ao.setAnnotation('LineStyle', '%s' % s) ao.setAnnotation('LineColor', '%s' % c) def setOptions(ao, options): "Set arbitrary annotations" for opt in options: key, val = opt.split('=', 1) ao.setAnnotation(key, val) # TODO: move to rivet.utils def mkoutdir(outdir): "Function to make output directories" if not os.path.exists(outdir): try: os.makedirs(outdir) except: msg = "Can't make output directory '%s'" % outdir raise Exception(msg) if not os.access(outdir, os.W_OK): msg = "Can't write to output directory '%s'" % outdir raise Exception(msg) def mkOutput(hpath, aos, plot=None, special=None): """ Make the .dat file string. We can't use "yoda.writeFLAT(anaobjects, 'foobar.dat')" because the PLOT and SPECIAL blocks don't have a corresponding analysis object. """ output = '' if plot is not None: output += str(plot) if special is not None: output += "\n" output += "# BEGIN SPECIAL %s\n" % hpath output += special output += "# END SPECIAL\n\n" from io import StringIO sio = StringIO() yoda.writeFLAT(aos, sio) output += sio.getvalue() return output def writeOutput(output, h): "Choose output file name and dir" if args.HIER_OUTPUT: hparts = h.strip("/").split("/", 1) ana = "_".join(hparts[:-1]) if len(hparts) > 1 else "ANALYSIS" outdir = os.path.join(args.OUTDIR, ana) outfile = '%s.dat' % hparts[-1].replace("/", "_") else: hparts = h.strip("/").split("/") outdir = args.OUTDIR outfile = '%s.dat' % "_".join(hparts) mkoutdir(outdir) outfilepath = os.path.join(outdir, outfile) f = open(outfilepath, 'w') f.write(output) f.close() #-------------------------------------------------------------------------------------------- if __name__ == '__main__': ## Command line parsing import argparse parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("ARGS", nargs="*") parser.add_argument('-o', '--outdir', dest='OUTDIR', default='.', help='write data files into this directory') parser.add_argument("--hier-out", action="store_true", dest="HIER_OUTPUT", default=False, help="write output dat files into a directory hierarchy which matches the analysis paths") parser.add_argument('--plotinfodir', dest='PLOTINFODIRS', action='append', default=['.'], help='directory which may contain plot header information (in addition ' 'to standard Rivet search paths)') parser.add_argument("--no-rivet-refs", dest="RIVETREFS", action="store_false", default=True, help="don't use Rivet reference data files") # parser.add_argument("--refid", dest="REF_ID", # default="REF", help="ID of reference data set (file path for non-REF data)") parser.add_argument("--reftitle", dest="REFTITLE", default='Data', help="Reference data legend entry") parser.add_argument("--pwd", dest="PATH_PWD", action="store_true", default=False, help="append the current directory (pwd) to the analysis/data search paths (cf. $RIVET_ANALYSIS/DATA_PATH)") parser.add_argument("-v", "--verbose", dest="VERBOSE", action="store_true", default=False, help="produce debug output to the terminal") stygroup = parser.add_argument_group("Plot style") stygroup.add_argument("--linear", action="store_true", dest="LINEAR", default=False, help="plot with linear scale") stygroup.add_argument("--errs", "--mcerrs", "--mc-errs", action="store_true", dest="MC_ERRS", default=False, help="show vertical error bars on the MC lines") stygroup.add_argument("--no-ratio", action="store_false", dest="RATIO", default=True, help="disable the ratio plot") stygroup.add_argument("--rel-ratio", action="store_true", dest="RATIO_DEVIATION", default=False, help="show the ratio plots scaled to the ref error") stygroup.add_argument("--no-plottitle", action="store_true", dest="NOPLOTTITLE", default=False, help="don't show the plot title on the plot " "(useful when the plot description should only be given in a caption)") stygroup.add_argument("--style", dest="STYLE", default="default", help="change plotting style: default|bw|talk") stygroup.add_argument("-c", "--config", dest="CONFIGFILES", action="append", default=["~/.make-plots"], help="additional plot config file(s). Settings will be included in the output configuration.") + stygroup.add_argument("--remove-options", help="remove options label from legend", dest="REMOVE_OPTIONS", + action="store_true", default=False) selgroup = parser.add_argument_group("Selective plotting") # selgroup.add_argument("--show-single", dest="SHOW_SINGLE", choices=("no", "ref", "mc", "all"), # default="mc", help="control if a plot file is made if there is only one dataset to be plotted " # "[default=%(default)s]. If the value is 'no', single plots are always skipped, for 'ref' and 'mc', " # "the plot will be written only if the single plot is a reference plot or an MC " # "plot respectively, and 'all' will always create single plot files.\n The 'ref' and 'all' values " # "should be used with great care, as they will also write out plot files for all reference " # "histograms without MC traces: combined with the -R/--rivet-refs flag, this is a great way to " # "write out several thousand irrelevant reference data histograms!") # selgroup.add_argument("--show-mc-only", "--all", action="store_true", dest="SHOW_IF_MC_ONLY", # default=False, help="make a plot file even if there is only one dataset to be plotted and " # "it is an MC one. Deprecated and will be removed: use --show-single instead, which overrides this.") # # selgroup.add_argument("-l", "--histogram-list", dest="HISTOGRAMLIST", # # default=None, help="specify a file containing a list of histograms to plot, in the format " # # "/ANALYSIS_ID/histoname, one per line, e.g. '/DELPHI_1996_S3430090/d01-x01-y01'.") selgroup.add_argument("-m", "--match", action="append", help="only write out histograms whose $path/$name string matches these regexes. The argument " "may also be a text file.", dest="PATHPATTERNS") selgroup.add_argument("-M", "--unmatch", action="append", help="exclude histograms whose $path/$name string matches these regexes", dest="PATHUNPATTERNS") args = parser.parse_args() ## Add pwd to search paths if args.PATH_PWD: rivet.addAnalysisLibPath(os.path.abspath(".")) rivet.addAnalysisDataPath(os.path.abspath(".")) ## Split the input file names and the associated plotting options ## given on the command line into two separate lists filelist, plotoptions = parseArgs(args.ARGS) ## Remove the PLOT dummy file from the file list if "PLOT" in filelist: filelist.remove("PLOT") ## Check that the files exist for f in filelist: if not os.access(f, os.R_OK): print("Error: cannot read from %s" % f) sys.exit(1) ## Read the .plot files plotdirs = args.PLOTINFODIRS + [os.path.abspath(os.path.dirname(f)) for f in filelist] plotparser = rivet.mkStdPlotParser(plotdirs, args.CONFIGFILES) ## Create a list of all histograms to be plotted, and identify if they are 2D histos (which need special plotting) try: refhistos, mchistos = getHistos(filelist) except IOError as e: print("File reading error: ", e.strerror) exit(1) hpaths, h2ds = [], [] for aos in mchistos.values(): for p in aos.keys(): ps = rivet.stripOptions(p) if ps and ps not in hpaths: hpaths.append(ps) firstaop = aos[p][sorted(aos[p].keys())[0]] # TODO: Would be nicer to test via isHisto and dim or similar, or yoda.Scatter/Histo/Profile base classes if type(firstaop) in (yoda.Histo2D, yoda.Profile2D) and ps not in h2ds: h2ds.append(ps) ## Take reference data from the Rivet search paths, if there is not already if args.RIVETREFS: try: refhistos2 = getRivetRefData() except IOError as e: print("File reading error: ", e.strerror) exit(1) refhistos2.update(refhistos) refhistos = refhistos2 ## Purge unmatched ref data entries to save memory keylist = list(refhistos.keys()) # can't modify for-loop target for refhpath in keylist: if refhpath not in hpaths: del refhistos[refhpath] ## Now loop over all MC histograms and plot them # TODO: factorize much of this into a rivet.utils mkplotfile(mchists, refhist, kwargs, is2d=False) function for hpath in hpaths: #print('Currently looking at', hpath) ## The analysis objects to be plotted anaobjects = [] ## List of histos to be drawn, to sync the legend and plotted lines mainlines = [] varlines = [] ## Is this a 2D histo? is2d = (hpath in h2ds) ## Will we be drawing a ratio plot? showratio = args.RATIO and not is2d ## A Plot object to represent the PLOT section in the .dat file plot = Plot() if not is2d: plot['Legend'] = '1' plot['LogY'] = '1' headers = plotparser.getHeaders(hpath) if headers: plot.update(headers) # for key, val in headers.items(): # plot[key] = val if "PLOT" in plotoptions: for key_val in plotoptions["PLOT"]: key, val = [s.strip() for s in key_val.split("=", 1)] if 'ReplaceOption' in key_val: opt_group = key_val.split("=", 2) key, val = "=".join(opt_group[:2]), opt_group[2] plot[key] = val + if args.REMOVE_OPTIONS: + plot['RemoveOptions'] = '1' if args.LINEAR: plot['LogY'] = '0' if args.NOPLOTTITLE: plot['Title'] = '' if showratio and args.RATIO_DEVIATION: plot['RatioPlotMode'] = 'deviation' if args.STYLE == 'talk': plot['PlotSize'] = '8,6' elif args.STYLE == 'bw' and showratio: plot['RatioPlotErrorBandColor'] = 'black!10' ## Get a special object, if there is one for this path special = plotparser.getSpecial(hpath) ## Handle reference data histogram, if there is one ratioreference, hasdataref = None, False if hpath in refhistos: hasdataref = True refdata = refhistos[hpath] refdata.setAnnotation('Title', args.REFTITLE) if not is2d: refdata.setAnnotation('ErrorBars', '1') refdata.setAnnotation('PolyMarker', '*') refdata.setAnnotation('ConnectBins', '0') if showratio: ratioreference = hpath ## For 1D anaobjects.append(refdata) mainlines.append(hpath) ## For 2D if is2d: s = mkOutput(hpath, [refdata], plot, special) writeOutput(s, hpath) ## Loop over the MC files to plot all instances of the histogram styleidx = 0 for infile in filelist: if infile in mchistos: for xpath in sorted(mchistos[infile]): if rivet.stripOptions(xpath) != hpath: continue hmcs = mchistos[infile][xpath] ## For now, just plot all the different variation histograms (reversed, so [0] is on top) # TODO: calculate and plot an appropriate error band, somehow... for i in sorted(hmcs.keys(), reverse=True): iscanonical = (str(i) == "0") hmc = hmcs[i] ## Default linecolor, linestyle if not is2d: setStyle(hmc, styleidx, not iscanonical) if args.MC_ERRS: hmc.setAnnotation('ErrorBars', '1') ## Plot defaults from .plot files histopts = plotparser.getHistogramOptions(hpath) if histopts: for key, val in histargs.items(): hmc.setAnnotation(key, val) ## Command line plot options setOptions(hmc, plotoptions[infile]) ## Set path attribute fullpath = "/"+infile+xpath if not iscanonical: fullpath += "["+str(i)+"]" hmc.setAnnotation('Path', fullpath) ## Add object / path to appropriate lists #if hmc.hasAnnotation("Title"): # hmc.setAnnotation("Title", hmc.annotation("Title") + # rivet.extractOptionString(xpath)) anaobjects.append(hmc) if iscanonical: mainlines.append(fullpath) else: varlines.append(fullpath) if showratio and ratioreference is None and iscanonical: ratioreference = fullpath ## For 2D, plot each histo now (since overlay makes no sense) if is2d: s = mkOutput(hpath, [hmc], plot, special) writeOutput(s, fullpath) styleidx += 1 ## Finally render the combined plots; only show the first one if it's 2D # TODO: Only show the first *MC* one if 2D? if is2d: anaobjects = anaobjects[:1] ## Add final attrs to Plot if not plot.has_key('DrawOnly'): plot['DrawOnly'] = ' '.join(varlines + mainlines).strip() if not plot.has_key('LegendOnly'): plot['LegendOnly'] = ' '.join(mainlines).strip() if showratio and not plot.has_key('RatioPlot') and len(varlines + mainlines) > 1: plot['RatioPlot'] = '1' if showratio and plot.has_key('RatioPlot') and plot['RatioPlot'] == '1' and len(varlines + mainlines) > 0: if not plot.has_key('RatioPlotReference'): plot['RatioPlotReference'] = ratioreference if not hasdataref and "RatioPlotYLabel" not in plot: if plot.get('RatioPlotMode', '') == 'deviation': plot['RatioPlotYLabel'] = 'Deviation' #r'$\text{MC}-\text{MC}_\text{ref}$' else: plot['RatioPlotYLabel'] = 'Ratio' #r'$\text{MC}/\text{MC}_\text{ref}$' ## Make the output and write to file o = mkOutput(hpath, anaobjects, plot, special) writeOutput(o, hpath) diff --git a/bin/rivet-mkhtml b/bin/rivet-mkhtml --- a/bin/rivet-mkhtml +++ b/bin/rivet-mkhtml @@ -1,506 +1,510 @@ #! /usr/bin/env python """\ %(prog)s [options] [ ...] [PLOT:Key1=Val1:...] Make web pages from histogram files written out by Rivet. You can specify multiple Monte Carlo YODA files to be compared in the same syntax as for rivet-cmphistos, 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). Any existing output directory will be overwritten. ENVIRONMENT: * RIVET_ANALYSIS_PATH: list of paths to be searched for analysis plugin libraries * RIVET_DATA_PATH: list of paths to be searched for data files """ from __future__ import print_function import rivet, sys, os rivet.util.check_python_version() rivet.util.set_process_name(os.path.basename(__file__)) import glob, shutil from subprocess import Popen,PIPE import argparse parser = argparse.ArgumentParser(usage=__doc__) parser.add_argument("YODAFILES", nargs="+", help="data files to compare") parser.add_argument("-o", "--outputdir", dest="OUTPUTDIR", default="./rivet-plots", help="directory for Web page output") parser.add_argument("-c", "--config", dest="CONFIGFILES", action="append", default=["~/.make-plots"], help="plot config file(s) to be used with rivet-cmphistos") parser.add_argument("-n", "--num-threads", metavar="NUMTHREADS", dest="NUMTHREADS", type=int, default=None, help="request make-plots to use a specific number of threads") parser.add_argument("--ignore-missing", dest="IGNORE_MISSING", action="store_true", default=False, help="ignore missing YODA files") parser.add_argument("-i", "--ignore-unvalidated", dest="IGNORE_UNVALIDATED", action="store_true", default=False, help="ignore unvalidated analyses") # parser.add_argument("--ref", "--refid", dest="REF_ID", # default=None, help="ID of reference data set (file path for non-REF data)") parser.add_argument("--dry-run", help="don't actually do any plotting or HTML building", dest="DRY_RUN", action="store_true", default=False) parser.add_argument("--no-cleanup", dest="NO_CLEANUP", action="store_true", default=False, help="keep plotting temporary directory") parser.add_argument("--no-subproc", dest="NO_SUBPROC", action="store_true", default=False, help="don't use subprocesses to render the plots in parallel -- useful for debugging") parser.add_argument("--pwd", dest="PATH_PWD", action="store_true", default=False, help="append the current directory (pwd) to the analysis/data search paths (cf. $RIVET_ANALYSIS_PATH)") stygroup = parser.add_argument_group("Style options") stygroup.add_argument("-t", "--title", dest="TITLE", default="Plots from Rivet analyses", help="title to be displayed on the main web page") stygroup.add_argument("--reftitle", dest="REFTITLE", default="Data", help="legend entry for reference data") stygroup.add_argument("--no-plottitle", dest="NOPLOTTITLE", action="store_true", default=False, help="don't show the plot title on the plot " "(useful when the plot description should only be given in a caption)") stygroup.add_argument("-s", "--single", dest="SINGLE", action="store_true", default=False, help="display plots on single webpage.") stygroup.add_argument("--no-ratio", dest="SHOW_RATIO", action="store_false", default=True, help="don't draw a ratio plot under each main plot.") stygroup.add_argument("--errs", "--mcerrs", "--mc-errs", dest="MC_ERRS", action="store_true", default=False, help="plot error bars.") stygroup.add_argument("--offline", dest="OFFLINE", action="store_true", default=False, help="generate HTML that does not use external URLs.") stygroup.add_argument("--pdf", dest="VECTORFORMAT", action="store_const", const="PDF", default="PDF", help="use PDF as the vector plot format.") stygroup.add_argument("--ps", dest="VECTORFORMAT", action="store_const", const="PS", default="PDF", help="use PostScript as the vector plot format. DEPRECATED") stygroup.add_argument("--booklet", dest="BOOKLET", action="store_true", default=False, help="create booklet (currently only available for PDF with pdftk or pdfmerge).") stygroup.add_argument("--font", dest="OUTPUT_FONT", choices="palatino,cm,times,helvetica,minion".split(","), default="palatino", help="choose the font to be used in the plots") stygroup.add_argument("--palatino", dest="OUTPUT_FONT", action="store_const", const="palatino", default="palatino", help="use Palatino as font (default). DEPRECATED: Use --font") stygroup.add_argument("--cm", dest="OUTPUT_FONT", action="store_const", const="cm", default="palatino", help="use Computer Modern as font. DEPRECATED: Use --font") stygroup.add_argument("--times", dest="OUTPUT_FONT", action="store_const", const="times", default="palatino", help="use Times as font. DEPRECATED: Use --font") stygroup.add_argument("--helvetica", dest="OUTPUT_FONT", action="store_const", const="helvetica", default="palatino", help="use Helvetica as font. DEPRECATED: Use --font") stygroup.add_argument("--minion", dest="OUTPUT_FONT", action="store_const", const="minion", default="palatino", help="use Adobe Minion Pro as font. DEPRECATED: Use --font") +stygroup.add_argument("--remove-options", help="remove options label from legend", dest="REMOVE_OPTIONS", + action="store_true", default=False) selgroup = parser.add_argument_group("Selective plotting") selgroup.add_argument("-m", "--match", action="append", dest="PATHPATTERNS", default=[], help="only write out histograms whose $path/$name string matches any of these regexes") selgroup.add_argument("-M", "--unmatch", action="append", dest="PATHUNPATTERNS", default=[], help="exclude histograms whose $path/$name string matches any of these regexes") selgroup.add_argument("--ana-match", action="append", dest="ANAPATTERNS", default=[], help="only write out histograms from analyses whose name matches any of these regexes") selgroup.add_argument("--ana-unmatch", action="append", dest="ANAUNPATTERNS", default=[], help="exclude histograms from analyses whose name matches any of these regexes") vrbgroup = parser.add_argument_group("Verbosity control") vrbgroup.add_argument("-v", "--verbose", help="add extra debug messages", dest="VERBOSE", action="store_true", default=False) args = parser.parse_args() yodafiles = args.YODAFILES ## Add pwd to search paths if args.PATH_PWD: rivet.addAnalysisLibPath(os.path.abspath(".")) rivet.addAnalysisDataPath(os.path.abspath(".")) ## Check that there are some arguments! if not yodafiles: print("Error: You need to specify some YODA files to be plotted!") sys.exit(1) ## Make output directory if not args.DRY_RUN: if os.path.exists(args.OUTPUTDIR) and not os.path.realpath(args.OUTPUTDIR)==os.getcwd(): import shutil shutil.rmtree(args.OUTPUTDIR) try: os.makedirs(args.OUTPUTDIR) except: print("Error: failed to make new directory '%s'" % args.OUTPUTDIR) sys.exit(1) ## Get set of analyses involved in the runs plotarg = None analyses = set() blocked_analyses = set() import yoda for yodafile in yodafiles: if yodafile.startswith("PLOT:"): plotarg = yodafile continue yodafilepath = os.path.abspath(yodafile.split(":")[0]) if not os.access(yodafilepath, os.R_OK): print("Error: cannot read from %s" % yodafilepath) if args.IGNORE_MISSING: continue else: sys.exit(2) try: ## Note: we use -m/-M flags here as well as when calling rivet-cmphistos, to potentially speed this initial loading analysisobjects = yoda.read(yodafilepath, patterns=args.PATHPATTERNS, unpatterns=args.PATHUNPATTERNS) except IOError as e: print("File reading error: ", e.strerror) sys.exit(1) for path, ao in analysisobjects.items(): ## Make a path object and ensure the path is in standard form try: aop = rivet.AOPath(path) except Exception as e: #print(e) print("Found analysis object with non-standard path structure:", path, "... skipping") continue ## We don't plot data objects with path components hidden by an underscore prefix if aop.istmp() or aop.israw(): continue ## Identify analysis/histo name parts analysis = "ANALYSIS" if aop.basepathparts(keepref=False): analysis = rivet.stripOptions(aop.basepathparts(keepref=False)[0]) #< TODO: for compatibility with rivet-cmphistos... generalise? #analysis = "_".join(aop.dirnameparts(keepref=False)[:-1]) #< TODO: would this be nicer? Currently incompatible with rivet-cmphistos ## Optionally veto on analysis name pattern matching if analysis in analyses.union(blocked_analyses): continue import re matched = True if args.ANAPATTERNS: matched = False for patt in args.ANAPATTERNS: if re.match(patt, analysis) is not None: matched = True break if matched and args.ANAUNPATTERNS: for patt in args.ANAUNPATTERNS: if re.match(patt, analysis) is not None: matched = False break if matched: analyses.add(analysis) else: blocked_analyses.add(analysis) ## Sort analyses: group ascending by analysis name (could specialise grouping by collider), then ## descending by year, and finally descending by bibliographic archive ID code (INSPIRE first). def anasort(name): rtn = (1, name) if name.startswith("MC"): rtn = (99999999, name) else: stdparts = name.split("_") try: year = int(stdparts[1]) rtn = (0, stdparts[0], -year, 0) idcode = (0 if stdparts[2][0] == "I" else 1e10) - int(stdparts[2][1:]) rtn = (0, stdparts[0], -year, idcode) if len(stdparts) > 3: rtn += stdparts[3:] except: pass return rtn analyses = sorted(analyses, key=anasort) ## Uncomment to test analysis ordering on index page # print(analyses) # sys.exit(0) ## Run rivet-cmphistos to get plain .dat files from .yoda ## We do this here since it also makes the necessary directories ch_cmd = ["rivet-cmphistos"] if args.MC_ERRS: ch_cmd.append("--mc-errs") if not args.SHOW_RATIO: ch_cmd.append("--no-ratio") if args.NOPLOTTITLE: ch_cmd.append("--no-plottitle") +if args.REMOVE_OPTIONS: + ch_cmd.append("--remove-options") # if args.REF_ID is not None: # ch_cmd.append("--refid=%s" % os.path.abspath(args.REF_ID)) if args.REFTITLE: ch_cmd.append("--reftitle=%s" % args.REFTITLE ) if args.PATHPATTERNS: for patt in args.PATHPATTERNS: ch_cmd += ["-m", patt] #"'"+patt+"'"] if args.PATHUNPATTERNS: for patt in args.PATHUNPATTERNS: ch_cmd += ["-M", patt] #"'"+patt+"'"] ch_cmd.append("--hier-out") # TODO: Need to be able to override this: provide a --plotinfodir cmd line option? ch_cmd.append("--plotinfodir=%s" % os.path.abspath("../")) for af in yodafiles: yodafilepath = os.path.abspath(af.split(":")[0]) if af.startswith("PLOT:"): yodafilepath = "PLOT" elif not os.access(yodafilepath, os.R_OK): continue newarg = yodafilepath if ":" in af: newarg += ":" + af.split(":", 1)[1] # print(newarg) ch_cmd.append(newarg) ## Pass rivet-mkhtml -c args to rivet-cmphistos for configfile in args.CONFIGFILES: configfile = os.path.abspath(os.path.expanduser(configfile)) if os.access(configfile, os.R_OK): ch_cmd += ["-c", configfile] if args.VERBOSE: ch_cmd.append("--verbose") print("Calling rivet-cmphistos with the following command:") print(" ".join(ch_cmd)) ## Run rivet-cmphistos in a subdir, after fixing any relative paths in Rivet env vars if not args.DRY_RUN: for var in ("RIVET_ANALYSIS_PATH", "RIVET_DATA_PATH", "RIVET_REF_PATH", "RIVET_INFO_PATH", "RIVET_PLOT_PATH"): if var in os.environ: abspaths = [os.path.abspath(p) for p in os.environ[var].split(":")] os.environ[var] = ":".join(abspaths) subproc = Popen(ch_cmd, cwd=args.OUTPUTDIR, stdout=PIPE, stderr=PIPE) out, err = subproc.communicate() retcode = subproc.returncode if args.VERBOSE or retcode != 0: print('Output from rivet-cmphistos:\n', out) if err : print('Errors from rivet-cmphistos:\n', err) if retcode != 0: print('Crash in rivet-cmphistos code = ', retcode, ' exiting') exit(retcode) ## Write web page containing all (matched) plots ## Make web pages first so that we can load it locally in ## a browser to view the output before all plots are made if not args.DRY_RUN: style = ''' ''' ## Include MathJax configuration script = '' if not args.OFFLINE: script = '''\ ''' ## A helper function for metadata LaTeX -> HTML conversion from rivet.util import htmlify ## A timestamp HTML fragment to be used on each page: import datetime timestamp = '

Generated at %s

\n' % datetime.datetime.now().strftime("%A, %d. %B %Y %I:%M%p") index = open(os.path.join(args.OUTPUTDIR, "index.html"), "w") index.write('\n\n%s\n%s\n' % (args.TITLE, style + script)) if args.BOOKLET and args.VECTORFORMAT == "PDF": index.write('

%s

\n\n' % args.TITLE) else: index.write('

%s

\n\n' % args.TITLE) if args.SINGLE: ## Write table of contents index.write('
    \n') for analysis in analyses: summary = analysis ana = rivet.AnalysisLoader.getAnalysis(analysis) if ana: summary = "%s (%s)" % (ana.summary(), analysis) if args.IGNORE_UNVALIDATED and ana.status() != "VALIDATED": continue index.write('
  • %s\n' % (analysis, htmlify(summary)) ) index.write('
\n') for analysis in analyses: references = [] summary = htmlify(analysis) description, inspireid, spiresid = None, None, None if analysis.find("_I") > 0: inspireid = analysis[analysis.rfind('_I')+2:len(analysis)] elif analysis.find("_S") > 0: spiresid = analysis[analysis.rfind('_S')+2:len(analysis)] ana = rivet.AnalysisLoader.getAnalysis(analysis) if ana: if ana.summary(): summary = htmlify("%s (%s)" % (ana.summary(), analysis)) references = ana.references() description = htmlify(ana.description()) spiresid = ana.spiresId() if args.IGNORE_UNVALIDATED and ana.status().upper() != "VALIDATED": continue try: # a, s = (analysis.encode("utf-8"), summary.encode("utf-8")) if args.SINGLE: index.write('\n

%s

\n' % (analysis.encode("utf-8"), summary.encode("utf-8"))) else: index.write('\n

%s

\n' % (analysis.encode("utf-8"), summary.encode("utf-8"))) except UnicodeEncodeError as ue: print("Unicode error in analysis description for " + analysis + ": " + str(ue)) reflist = [] if inspireid: reflist.append('Inspire' % inspireid) reflist.append('HepData' % inspireid) elif spiresid: # elif ana.spiresId(): reflist.append('Inspire' % spiresid) reflist.append('HepData' % spiresid) reflist += references index.write('

%s

\n' % " | ".join(reflist)) if description: try: index.write('

%s

\n' % description) except UnicodeEncodeError as ue: print("Unicode error in analysis description for " + analysis + ": " + str(ue)) anapath = os.path.join(args.OUTPUTDIR, analysis) if not args.SINGLE: if not os.path.exists(anapath): os.makedirs(anapath) anaindex = open(os.path.join(anapath, "index.html"), 'w') anaindex.write('\n\n%s – %s\n%s\n\n' % (htmlify(args.TITLE), analysis, style + script)) anaindex.write('

%s

\n' % htmlify(analysis)) anaindex.write('

Back to index

\n') if description: try: anaindex.write('

\n %s\n

\n' % description) except UnicodeEncodeError as ue: print("Unicode error in analysis description for " + analysis + ": " + str(ue)) else: anaindex = index datfiles = glob.glob("%s/*.dat" % anapath) #print(datfiles) anaindex.write('
\n') for datfile in sorted(datfiles): obsname = os.path.basename(datfile).replace(".dat", "") pngfile = obsname+".png" vecfile = obsname+"."+args.VECTORFORMAT.lower() srcfile = obsname+".dat" if args.SINGLE: pngfile = os.path.join(analysis, pngfile) vecfile = os.path.join(analysis, vecfile) srcfile = os.path.join(analysis, srcfile) anaindex.write('
\n') anaindex.write(' ⚓⌘ %s:
\n' % (analysis, obsname, srcfile, os.path.splitext(vecfile)[0]) ) anaindex.write(' \n' % (analysis, obsname, vecfile) ) anaindex.write(' \n' % pngfile ) anaindex.write(' \n') anaindex.write('
\n') anaindex.write('
\n') if not args.SINGLE: anaindex.write('
%s\n
\n' % timestamp) anaindex.close() index.write('
%s\n' % timestamp) index.close() # http://stackoverflow.com/questions/377017/test-if-executable-exists-in-python def which(program): import os def is_exe(fpath): return os.path.isfile(fpath) and os.access(fpath, os.X_OK) fpath, fname = os.path.split(program) if fpath: if is_exe(program): return program else: for path in os.environ["PATH"].split(os.pathsep): path = path.strip('"') exe_file = os.path.join(path, program) if is_exe(exe_file): return exe_file return None ## Run make-plots on all generated .dat files # sys.exit(0) mp_cmd = ["make-plots"] if args.NUMTHREADS: mp_cmd.append("--num-threads=%d" % args.NUMTHREADS) if args.NO_CLEANUP: mp_cmd.append("--no-cleanup") if args.NO_SUBPROC: mp_cmd.append("--no-subproc") if args.VECTORFORMAT == "PDF": mp_cmd.append("--pdfpng") elif args.VECTORFORMAT == "PS": mp_cmd.append("--pspng") if args.OUTPUT_FONT: mp_cmd.append("--font=%s" % args.OUTPUT_FONT) # if args.OUTPUT_FONT.upper() == "PALATINO": # mp_cmd.append("--palatino") # if args.OUTPUT_FONT.upper() == "CM": # mp_cmd.append("--cm") # elif args.OUTPUT_FONT.upper() == "TIMES": # mp_cmd.append("--times") # elif args.OUTPUT_FONT.upper() == "HELVETICA": # mp_cmd.append("--helvetica") # elif args.OUTPUT_FONT.upper() == "MINION": # mp_cmd.append("--minion") datfiles = [] for analysis in analyses: anapath = os.path.join(args.OUTPUTDIR, analysis) #print(anapath) anadatfiles = glob.glob("%s/*.dat" % anapath) datfiles += sorted(anadatfiles) if datfiles: mp_cmd += datfiles if args.VERBOSE: mp_cmd.append("--verbose") print("Calling make-plots with the following options:") print(" ".join(mp_cmd)) if not args.DRY_RUN: Popen(mp_cmd).wait() if args.BOOKLET and args.VECTORFORMAT=="PDF": if which("pdftk") is not None: bookletcmd = ["pdftk"] for analysis in analyses: anapath = os.path.join(args.OUTPUTDIR, analysis) bookletcmd += sorted(glob.glob("%s/*.pdf" % anapath)) bookletcmd += ["cat", "output", "%s/booklet.pdf" % args.OUTPUTDIR] print(bookletcmd) Popen(bookletcmd).wait() elif which("pdfmerge") is not None: bookletcmd = ["pdfmerge"] for analysis in analyses: anapath = os.path.join(args.OUTPUTDIR, analysis) bookletcmd += sorted(glob.glob("%s/*.pdf" % anapath)) bookletcmd += ["%s/booklet.pdf" % args.OUTPUTDIR] print(bookletcmd) Popen(bookletcmd).wait() else: print("Neither pdftk nor pdfmerge available --- not booklet output possible") diff --git a/include/Rivet/Math/LorentzTrans.hh b/include/Rivet/Math/LorentzTrans.hh --- a/include/Rivet/Math/LorentzTrans.hh +++ b/include/Rivet/Math/LorentzTrans.hh @@ -1,296 +1,307 @@ #ifndef RIVET_MATH_LORENTZTRANS #define RIVET_MATH_LORENTZTRANS #include "Rivet/Math/MathHeader.hh" #include "Rivet/Math/MathUtils.hh" #include "Rivet/Math/MatrixN.hh" #include "Rivet/Math/Matrix3.hh" #include "Rivet/Math/Vector4.hh" #include namespace Rivet { /// @brief Object implementing Lorentz transform calculations and boosts. /// /// @note These boosts are defined actively, i.e. as modifications of vectors /// rather than frame transformations. So the boost vector is the opposite of /// what you might expect. /// /// @todo Review the active/passive convention choice. Seems counterintuitive now... class LorentzTransform { public: /// @name Simple Lorentz factor conversions //@{ /// Calculate the \f$ \gamma \f$ factor from \f$ \beta \f$ static double beta2gamma(double beta) { return 1.0 / sqrt(1 - sqr(beta)); } /// Calculate \f$ \beta \f$ from the \f$ \gamma \f$ factor static double gamma2beta(double gamma) { return sqrt(1 - sqr(1/gamma)); } // /// Calculate the \f$ \gamma \f$ factor from \f$ \bar{\beta}^2 = 1-\beta^2 \f$ // static double betabarsq2gamma(double betabarsq) { // return 1.0 / sqrt(betabarsq); // } // /// Calculate \f$ \bar{\beta}^2 = 1-\beta^2 \f$ from the \f$ \gamma \f$ factor // static double gamma2betabarsq(double gamma) { // return 1.0 / sqr(gamma); // } //@} /// @name Construction //@{ /// Default (identity) constructor LorentzTransform() { _boostMatrix = Matrix<4>::mkIdentity(); } + /// Constructor from a 4x4 matrix + LorentzTransform(const Matrix<4>& boostMatrix) { + _boostMatrix = boostMatrix; + } /// Make an LT for an active boost (i.e. object velocity += in boost direction) static LorentzTransform mkObjTransformFromBeta(const Vector3& vbeta) { LorentzTransform rtn; return rtn.setBetaVec(vbeta); } /// Make an LT for a passive boost (i.e. object velocity -= in boost direction) static LorentzTransform mkFrameTransformFromBeta(const Vector3& vbeta) { LorentzTransform rtn; return rtn.setBetaVec(-vbeta); } /// Make an LT for an active boost (i.e. object velocity += in boost direction) static LorentzTransform mkObjTransformFromGamma(const Vector3& vgamma) { LorentzTransform rtn; if (!vgamma.isZero()) rtn.setGammaVec(vgamma); return rtn; } /// Make an LT for a passive boost (i.e. object velocity -= in boost direction) static LorentzTransform mkFrameTransformFromGamma(const Vector3& vgamma) { LorentzTransform rtn; if (!vgamma.isZero()) rtn.setGammaVec(-vgamma); return rtn; } /// Make an LT for an active boost (i.e. object velocity += in boost direction) static LorentzTransform mkObjTransform(const FourMomentum& p4) { return mkObjTransformFromBeta(p4.betaVec()); } /// Make an LT for a passive boost (i.e. object velocity -= in boost direction) static LorentzTransform mkFrameTransform(const FourMomentum& p4) { return mkObjTransformFromBeta(-p4.betaVec()); } //@} /// @name Boost vector and beta/gamma factors //@{ + /// Set up an active Lorentz boost from a (unit) direction vector, and \f$ \beta \f$ & \f$ \gamma \f$ factors + LorentzTransform& _setBoost(const Vector3& vec, double beta, double gamma) { + // Set to identity for null boosts + _boostMatrix = Matrix<4>::mkIdentity(); + if (isZero(beta)) return *this; + // + // It's along the x, y, or z axis if 2 Cartesian components are zero + const bool alongxyz = (int(vec.x() == 0) + int(vec.y() == 0) + int(vec.z() == 0) == 2); + const int i = (!alongxyz || vec.x() != 0) ? 1 : (vec.y() != 0) ? 2 : 3; + const int isign = !alongxyz ? 1 : sign(vec[i-1]); + // + _boostMatrix.set(0, 0, gamma); + _boostMatrix.set(i, i, gamma); + _boostMatrix.set(0, i, +isign*beta*gamma); //< +ve coeff since active boost + _boostMatrix.set(i, 0, +isign*beta*gamma); //< +ve coeff since active boost + // + if (!alongxyz) _boostMatrix = rotate(Vector3::mkX(), vec)._boostMatrix; + return *this; + } + /// Set up an active Lorentz boost from the \f$ \vec\beta \f$ vector LorentzTransform& setBetaVec(const Vector3& vbeta) { // Set to identity for null boosts _boostMatrix = Matrix<4>::mkIdentity(); - if (vbeta.mod2() == 0) return *this; - assert(vbeta.mod2() < 1); - // - // It's along the x, y, or z axis if 2 Cartesian components are zero - const bool alongxyz = (int(isZero(vbeta.x())) + int(isZero(vbeta.y())) + int(isZero(vbeta.z())) == 2); - const int i = !isZero(vbeta.x()) ? 1 : !isZero(vbeta.y()) ? 2 : !isZero(vbeta.z()) ? 3 : 1; + if (isZero(vbeta.mod2())) return *this; const double beta = vbeta.mod(); const double gamma = beta2gamma(beta); - _boostMatrix.set(0, 0, gamma); - _boostMatrix.set(i, i, gamma); - _boostMatrix.set(0, i, +beta*gamma); //< +ve coeff since active boost - _boostMatrix.set(i, 0, +beta*gamma); //< +ve coeff since active boost - if (!alongxyz) _boostMatrix = rotate(Vector3::mkX(), vbeta)._boostMatrix; - return *this; + return _setBoost(vbeta.unit(), beta, gamma); } /// Get the \f$ \vec\beta \f$ vector for an active Lorentz boost Vector3 betaVec() const { FourMomentum boost(_boostMatrix.getColumn(0)); //< @todo WRONG?! //cout << "!!!" << boost << endl; if (boost.isZero()) return Vector3(); assert(boost.E() > 0); const double beta = boost.p3().mod() / boost.E(); return boost.p3().unit() * beta; } /// Get the \f$ \beta \f$ factor double beta() const { return betaVec().mod(); } /// Set up an active Lorentz boost from the \f$ \vec\gamma \f$ vector LorentzTransform& setGammaVec(const Vector3& vgamma) { + // Set to identity for null boosts + _boostMatrix = Matrix<4>::mkIdentity(); + if (isZero(vgamma.mod2() - 1)) return *this; const double gamma = vgamma.mod(); const double beta = gamma2beta(gamma); - _boostMatrix = Matrix<4>::mkIdentity(); - _boostMatrix.set(0, 0, gamma); - _boostMatrix.set(1, 1, gamma); - _boostMatrix.set(0, 1, +beta*gamma); //< +ve coeff since active boost - _boostMatrix.set(1, 0, +beta*gamma); //< +ve coeff since active boost - if (beta > 0) _boostMatrix = rotate(Vector3::mkX(), vgamma)._boostMatrix; - return *this; + return _setBoost(vgamma.unit(), beta, gamma); } /// Get the \f$ \vec\gamma \f$ vector for an active Lorentz boost Vector3 gammaVec() const { FourMomentum boost(_boostMatrix.getColumn(0)); //< @todo WRONG?! if (boost.isZero()) return Vector3(); assert(boost.E() > 0); const double beta = boost.p3().mod() / boost.E(); return boost.p3().unit() * beta; } /// Get the \f$ \gamma \f$ factor double gamma() const { return beta2gamma(beta()); } //@} /// Apply this transformation to the given 4-vector FourVector transform(const FourVector& v4) const { return multiply(_boostMatrix, v4); } /// Apply this transformation to the given 4-mometum FourMomentum transform(const FourMomentum& v4) const { return multiply(_boostMatrix, v4); } /// Apply this transformation to the given 4-vector FourVector operator () (const FourVector& v4) const { return transform(v4); } /// Apply this transformation to the given 4-mometum FourMomentum operator () (const FourMomentum& v4) const { return transform(v4); } /// @name Transform modifications //@{ /// Rotate the transformation cf. the difference between vectors @a from and @a to LorentzTransform rotate(const Vector3& from, const Vector3& to) const { return rotate(Matrix3(from, to)); } /// Rotate the transformation by @a angle radians about @a axis LorentzTransform rotate(const Vector3& axis, double angle) const { return rotate(Matrix3(axis, angle)); } /// Rotate the transformation by the 3D rotation matrix @a rot LorentzTransform rotate(const Matrix3& rot) const { LorentzTransform lt = *this; const Matrix4 rot4 = _mkMatrix4(rot); const Matrix4 newlt = rot4 * _boostMatrix * rot4.inverse(); lt._boostMatrix = newlt; return lt; } /// Calculate the inverse transform LorentzTransform inverse() const { LorentzTransform rtn; rtn._boostMatrix = _boostMatrix.inverse(); return rtn; } /// Combine LTs, treating @a this as the LH matrix. LorentzTransform combine(const LorentzTransform& lt) const { LorentzTransform rtn; rtn._boostMatrix = _boostMatrix * lt._boostMatrix; return rtn; } /// Operator combination of two LTs LorentzTransform operator*(const LorentzTransform& lt) const { return combine(lt); } /// Pre-multiply m3 by this LT LorentzTransform preMult(const Matrix3& m3) { _boostMatrix = multiply(_mkMatrix4(m3),_boostMatrix); return *this; } /// Post-multiply m3 by this LT LorentzTransform postMult(const Matrix3& m3) { _boostMatrix *= _mkMatrix4(m3); return *this; } //@} /// Return the matrix form Matrix4 toMatrix() const { return _boostMatrix; } private: Matrix4 _mkMatrix4(const Matrix3& m3) const { Matrix4 m4 = Matrix4::mkIdentity(); for (size_t i = 0; i < 3; ++i) { for (size_t j = 0; j < 3; ++j) { m4.set(i+1, j+1, m3.get(i, j)); } } return m4; } Matrix4 _boostMatrix; }; inline LorentzTransform inverse(const LorentzTransform& lt) { return lt.inverse(); } inline LorentzTransform combine(const LorentzTransform& a, const LorentzTransform& b) { return a.combine(b); } inline FourVector transform(const LorentzTransform& lt, const FourVector& v4) { return lt.transform(v4); } ////////////////////////// inline string toString(const LorentzTransform& lt) { return toString(lt.toMatrix()); } inline ostream& operator<<(std::ostream& out, const LorentzTransform& lt) { out << toString(lt); return out; } } #endif diff --git a/include/Rivet/Projections/EventMixingFinalState.hh b/include/Rivet/Projections/EventMixingFinalState.hh --- a/include/Rivet/Projections/EventMixingFinalState.hh +++ b/include/Rivet/Projections/EventMixingFinalState.hh @@ -1,99 +1,99 @@ // -*- C++ -*- #ifndef RIVET_EventMixingFinalState_HH #define RIVET_EventMixingFinalState_HH #include "Rivet/Projection.hh" #include "Rivet/Projections/ParticleFinder.hh" #include #include namespace Rivet { // @brief Projects out an event mixed of several events, given // a mixing observable (eg. number of final state particles), // defining what should qualify as a mixable event. // Binning in the mixing observable is defined in the constructor, // as is the number of events one wants to mix with. // The protected method calculateMixingObs() can be overloaded // in derived classes, to define other mixing observables, eg. // centrality or something even more elaborate. // // !!!DISCLAIMER!!! The projection makes no attempt at correct handling // of event weights - ie. what one should use as event weight for several // mixed events. Proceed with caution if you do not use an MC with // unit weights. // // @author Christian Bierlich typedef map > MixMap; class EventMixingFinalState : public Projection { public: // Constructor EventMixingFinalState(const ParticleFinder& fsp, const ParticleFinder& mix, size_t nMixIn, double oMin, double oMax, double deltao ) : nMix(nMixIn){ setName("EventMixingFinalState"); addProjection(fsp,"FS"); addProjection(mix,"MIX"); MSG_WARNING("EventMixingFinalState is a naive implementation, not currently " << "validated. Use with caution."); // Set up the map for mixing events for(double o = oMin; o < oMax; o+=deltao ) mixEvents[o] = deque(); } // Clone on the heap DEFAULT_RIVET_PROJ_CLONE(EventMixingFinalState); // Return a vector of mixing events. vector getMixingEvents() const { MixMap::const_iterator mixItr = mixEvents.lower_bound(mObs); if(mixItr == mixEvents.end() || mixItr->second.size() < nMix + 1) return vector(); return vector(mixItr->second.begin(), mixItr->second.end() - 1); } protected: // Calulate mixing observable. // Can be overloaded to define more elaborate observables. virtual void calculateMixingObs(const Particles& parts){ mObs = parts.size(); } /// Perform the projection on the Event - void project(const Event& e){ + virtual void project(const Event& e){ const Particles parts = applyProjection(e, "FS").particles(); calculateMixingObs(parts); MixMap::iterator mixItr = mixEvents.lower_bound(mObs); if(mixItr == mixEvents.end()){ // We are out of bounds. MSG_DEBUG("Mixing observable out of bounds."); return; } const Particles mix = applyProjection(e, "MIX").particles(); mixItr->second.push_back(mix); if(mixItr->second.size() > nMix + 1) mixItr->second.pop_front(); } /// Compare with other projections - int compare(const Projection& p) const { + virtual int compare(const Projection& p) const { return mkNamedPCmp(p,"FS"); } - private: + protected: // The number of event to mix with size_t nMix; // The mixing observable of the current event double mObs; // The event map; MixMap mixEvents; }; } #endif diff --git a/test/testMatVec.cc b/test/testMatVec.cc --- a/test/testMatVec.cc +++ b/test/testMatVec.cc @@ -1,155 +1,197 @@ #include "Rivet/Math/MathUtils.hh" #include "Rivet/Math/Vectors.hh" #include "Rivet/Math/Matrices.hh" -//#include "Rivet/Math/MatrixDiag.hh" +#include "Rivet/Math/MatrixDiag.hh" using namespace Rivet; #include #include #include using namespace std; + +/// Set up an active Lorentz boost in the z direction +LorentzTransform mkLTz(const double beta) { + assert(abs(beta) < 1); + const double gamma = LorentzTransform::beta2gamma(beta); + Matrix<4> rtn = Matrix<4>::mkIdentity(); + rtn.set(0, 0, gamma); + rtn.set(3, 3, gamma); + rtn.set(0, 3, +beta*gamma); //< +ve coeff since active boost + rtn.set(3, 0, +beta*gamma); //< +ve coeff since active boost + return LorentzTransform(rtn); +} + + int main() { FourVector a(1,0,0,0); cout << a << ": interval = " << a.invariant() << endl; assert(fuzzyEquals(a.invariant(), 1)); a.setZ(1); assert(isZero(a.invariant())); cout << a << ": interval = " << a.invariant() << endl; a.setY(2).setZ(3); cout << a << ": interval = " << a.invariant() << endl; assert(fuzzyEquals(a.invariant(), -12)); cout << a << ": vector = " << a.vector3() << endl << endl; FourMomentum b(1,0,0,0); cout << b << ": mass = " << b.mass() << endl; assert(fuzzyEquals(b.mass2(), 1)); b.setPz(1); cout << b << ": mass = " << b.mass() << endl; assert(isZero(b.mass2())); b.setPy(2).setPz(3).setE(6); cout << b << ": mass = " << b.mass2() << endl; assert(fuzzyEquals(b.mass2(), 23)); cout << b << ": vector = " << b.vector3() << endl << endl; Matrix3 m; m.set(0, 0, 7/4.0); m.set(0, 1, 3 * sqrt(3)/4.0); m.set(1, 0, 3 * sqrt(3)/4.0); m.set(1, 1, 13/4.0); m.set(2, 2, 9); cout << m << endl << endl; // EigenSystem<3> es = diagonalize(m); cout << "Matrices:" << endl; cout << Matrix3() << endl; cout << Matrix3::mkIdentity() << endl; const Matrix3 I3 = Matrix3::mkIdentity(); cout << Matrix3::mkIdentity() * m * I3 << endl; cout << "tr(0) & det(0): " << Matrix3().trace() << ", " << Matrix3().det() << endl; cout << "tr(I3) & det(I3): " << I3.trace() << ", " << I3.det() << endl; Matrix3 m1 = Matrix3::mkIdentity(); Matrix3 m2 = m1; m1.setRow(1, Vector3(1,2,3)); m2.setColumn(1, Vector3(3,2,1)); Matrix3 m3 = Matrix3::mkZero(); cout << m1 << " + " << m2 << " = " << m1 + m2 << endl; m3.setRow(0, Vector3(2,3,0)).setRow(1, Vector3(1,4,3)).setRow(2, Vector3(0,1,2)); cout << m1+m2 << " == " << m3 << ": " << (m1+m2 == m3 ? "true" : "false") << endl; cout << endl; Vector3 v3(1,2,3); cout << "Vector: " << v3 << endl; cout << "Invert: " << v3 << " --> " << -v3 << endl; const Matrix3 rot90(Vector3(0,0,1), PI/2.0); const Matrix3 rot90m(Vector3(0,0,1), -PI/2.0); const Matrix3 rot180(Vector3(0,0,1), PI); const Matrix3 rot180m(Vector3(0,0,1), -PI); const Vector3 v3_90 = rot90*v3; cout << "Rot 90: " << v3 << " ---> " << v3_90 << endl; const Vector3 v3_90m = rot90m*v3; cout << "Rot -90: " << v3 << " ---> " << v3_90m << endl; const Vector3 v3_180 = rot180*v3; cout << "Rot 180: " << v3 << " ---> " << v3_180 << endl; const Vector3 v3_180m = rot180m*v3; cout << "Rot -180: " << v3 << " ---> " << v3_180m << endl; assert(fuzzyEquals(v3_180, v3_180m)); const Vector3 v3_9090 = rot90*rot90*v3; cout << "Rot 2 x 90: " << v3 << " ---> " << v3_9090 << endl; assert(fuzzyEquals(v3_180, v3_9090)); const Vector3 v3_90m90m = rot90m*rot90m*v3; cout << "Rot 2 x -90: " << v3 << " ---> " << v3_90m90m << endl; assert(fuzzyEquals(v3_180, v3_90m90m)); const Vector3 v3_9090m = rot90*rot90m*v3; const Vector3 v3_90m90 = rot90m*rot90*v3; cout << "Rot 90*-90: "<< v3 << " ---> " << v3_9090m << endl; cout << "Rot -90*90: "<< v3 << " ---> " << v3_90m90 << endl; assert(fuzzyEquals(v3, v3_9090m)); assert(fuzzyEquals(v3, v3_90m90)); const Vector3 v3_90i = rot90.inverse()*v3; cout << "Rot (90)^-1: "<< v3 << " ---> " << v3_90i << endl; assert(fuzzyEquals(v3_90i, v3_90m)); const Vector3 v3_9090i = rot90*rot90.inverse()*v3; const Vector3 v3_90i90 = rot90.inverse()*rot90*v3; cout << "Rot 90*(90)^-1: "<< v3 << " ---> " << v3_9090i << endl; cout << "Rot (90)^-1*90: "<< v3 << " ---> " << v3_90i90 << endl; assert(fuzzyEquals(v3, v3_9090i)); assert(fuzzyEquals(v3, v3_90i90)); const Matrix3 rot1(Vector3(0,1,0), PI/180.0); cout << "Rot 0 x 45 x 1: " << v3 << endl; for (size_t i = 0; i < 8; ++i) { for (size_t j = 0; j < 45; ++j) { v3 = rot1*v3; } cout << "Rot " << i+1 << " x 45 x 1: " << v3 << endl; } assert(fuzzyEquals(v3, Vector3(1,2,3))); cout << endl; cout << "Boosts:" << endl; LorentzTransform ltX = LorentzTransform::mkObjTransformFromBeta(Vector3(0.5, 0, 0)); cout << "LTx: " << ltX << endl; cout << "I on LTx: " << ltX.rotate(Matrix3::mkIdentity()) << endl; cout << "Rot90 on LTx: " << ltX.rotate(rot90) << endl; cout << endl; + LorentzTransform ltY = LorentzTransform::mkObjTransformFromBeta(Vector3(0, 0.5, 0)); + cout << "LTy: " << ltY << endl; + cout << "I on LTy: " << ltY.rotate(Matrix3::mkIdentity()) << endl; + cout << "Rot90 on LTy: " << ltY.rotate(rot90) << endl; + cout << endl; + + LorentzTransform ltZ = LorentzTransform::mkObjTransformFromBeta(Vector3(0, 0, 0.5)); + cout << "LTz: " << ltZ << endl; + cout << "I on LTz: " << ltZ.rotate(Matrix3::mkIdentity()) << endl; + cout << "Rot90 on LTz: " << ltZ.rotate(rot90) << endl; + cout << endl; + + LorentzTransform ltZ2 = mkLTz(0.5); + cout << "LTz': " << ltZ2 << endl; + cout << endl; + + LorentzTransform ltXYZ = LorentzTransform::mkObjTransformFromBeta(Vector3(0.1, -0.2, 0.3)); + cout << "LTxyz: " << ltXYZ << endl; + cout << "I on LTxyz: " << ltXYZ.rotate(Matrix3::mkIdentity()) << endl; + cout << "Rot90 on LTxyz: " << ltXYZ.rotate(rot90) << endl; + cout << endl; + cout << "X-boosts:" << endl; const FourMomentum p1 = FourMomentum(10, 0, 0, 1); const FourMomentum p2 = ltX.transform(p1); cout << p1 << " -> " << p2 << endl; cout << p2 << " -> " << ltX.inverse().transform(p2) << endl; //cout << p1.boostVector() << endl; + cout << "LT(p1) = " << LorentzTransform::mkFrameTransformFromBeta(p1.boostVector()) << endl; const FourMomentum p3 = LorentzTransform::mkFrameTransformFromBeta(p1.boostVector()).transform(p1); cout << p1 << " -> " << p3 << endl; + assert(isZero(p3.x())); + assert(isZero(p3.y())); + assert(isZero(p3.z())); + assert(p3.E() < p1.E()); + assert(fuzzyEquals(p3.E(), p1.mass())); cout << endl; - LorentzTransform ltY = LorentzTransform::mkObjTransformFromGamma(Vector3(0, 1.2, 0)); + // LorentzTransform ltY = LorentzTransform::mkObjTransformFromGamma(Vector3(0, 1.2, 0)); cout << FourMomentum(1,0,0,1) << " -> " //<< "\n " << (ltX * ltY).transform(FourMomentum(1,0,0,1)) << endl; cout << FourMomentum(1,0,0,1) << " -> " //<< "\n " << (ltY * ltX).transform(FourMomentum(1,0,0,1)) << endl; cout << (ltX * ltY).betaVec() << endl; cout << (ltY * ltX).betaVec() << endl; cout << (ltX * ltX.inverse()).betaVec() << endl; // If we are already in the rest frame and there is no boost, then LT is trivial/identity LorentzTransform noBoost; cout << "Element 0,0 should be 1 and is " << noBoost.toMatrix().get(0,0) << endl; assert(noBoost.toMatrix().get(0,0)==1); cout << "Element 0,1 should be 0 and is " << noBoost.toMatrix().get(0,1) << endl; assert(noBoost.toMatrix().get(0,1)==0); cout << "Element 1,0 should be 0 and is " << noBoost.toMatrix().get(1,0) << endl; assert(noBoost.toMatrix().get(1,0)==0); cout << "Element 1,1 should be 1 and is " << noBoost.toMatrix().get(1,1) << endl; assert(noBoost.toMatrix().get(1,1)==1); return EXIT_SUCCESS; }