common_tools.py

Go to the documentation of this file.

###########################
# CONSTANTS
STANDARD_POT = 6e20 
###########################
# print a failure message, from:
#  https://cdcvs.fnal.gov/redmine/projects/novaart/repository/entry/trunk/NovaGridUtils/lib/python/NovaGridUtils.py
def fail(msg):
    import sys
    print "\n\033[91mtools: --- ERROR"
    print msg
    print "\033[0m\n"
    sys.exit(1)
###########################
# print a warning message, from:
#  https://cdcvs.fnal.gov/redmine/projects/novaart/repository/entry/trunk/NovaGridUtils/lib/python/NovaGridUtils.py
def warn(msg):
    print "\n\033[93mtool: --- WARNING"
    print msg
    print "\033[0m\n"
########################### 
# Get jobsub client job id from submission log
def getJobID(log):
    return_id = False
    for l in log:
        if "JobsubJobId of first job" in l:
            return_id = l.split(":")[1].strip()
    if return_id == False:
        print "tools: getJobID failed. Given log:"
        for l in log: print l.strip()
    return return_id
########################### 
# Get SAM URL from submission log
def getSAMURL(log):
    return_id = False
    for l in log:
        if "Station monitor will be " in l:
            return_id = l.split("/")[-1].strip()
    if return_id == False:
        print "tools: getSAMURL failed. Given log:"
        for l in log: print l.strip()
    return return_id
########################### 
# Get SAM station monitor from submission log
def getStationMonitor(log):
    return_id = False
    for l in log:
        if "Station monitor will be " in l:
            return_id = l.split("Station monitor will be ")[-1].strip()
    if return_id == False:
        print "tools: getStationMonitor failed. Given log:"
        for l in log: print l.strip()
    return return_id
########################### 
# get creation time of a file from:
#   http://stackoverflow.com/questions/237079/how-to-get-file-creation-modification-date-times-in-python
def fileModificationDatetime(filename):
    import os
    import datetime
    t = os.path.getmtime(filename)
    return datetime.datetime.fromtimestamp(t)
########################### 
# fetch log files for a jobsub job
def fetchLogs(job_id, run_datetime=False, min_caching_age_days=10, verbose=False):
    import subprocess
    import datetime
    import getpass
    import time
    import glob
    import os
    if ("SETUP_JOBSUB_CLIENT" not in os.environ):
        print "tools: Jobsub client has not been setup. Cannot fetch logs without the tools to do this"
        print "tools: Aborting\n"
        exit(1)

    if job_id == False: return []
    if (not run_datetime): run_datetime = datetime.datetime.now()

    log_dir = "%s/jobsub_logs/%s"%(os.environ["NOVAANAVALID_GRID"],job_id)
    if verbose: print "tools: fetching logs to dir: %s"%log_dir
    fetch_logs = True
    # check for existing instance
    if os.path.exists(log_dir):
        if verbose: print "tools:     - log directory exists"
        # folder_datetime = fileModificationDatetime(log_dir)
        now = datetime.datetime.now()
        d_t = datetime.datetime.now() - run_datetime
        if verbose: print "tools:     run datetime: %s, now: %s, days: %i"%(run_datetime,now,d_t.days)
        if d_t.days < min_caching_age_days:
            if verbose: print "tools:     re-fetching logs"
            rm_command = "rm %s/*"%log_dir
            os.system(rm_command)
        else:
            if verbose: print "tools:     using existing logs"
            fetch_logs = False

    if fetch_logs:
        fetch_command = "jobsub_fetchlog -G nova --job=%s --unzipdir=%s"%(job_id,log_dir)
        if "novapro" in getpass.getuser(): fetch_command+=" --role=Production"

        process = subprocess.Popen(fetch_command.split(" "), stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
        while process.poll() == None: time.sleep(0.1)
        if process.returncode != 0:
            print "tools : --- Exception caught while executing fetch log, return code: %i"%process.returncode
            print "tools :     command run: %s"%fetch_command
            for line in process.stdout.readlines():
                print line.strip()
        # subprocess.check_call(fetch_command.split())
    
    file_list = glob.glob("%s/*"%log_dir)
    return file_list
########################### 
# find the root file output of a component
def findRootOutput(component, verbose=False, testing=False):
    if "glob" not in dir(): import glob
    import os
    output_dir = component.output_dir
    if testing: 
        output_dir = output_dir.replace(os.environ["NOVAANAVALID_GRID"],os.environ["NOVAANAVALID_GRID_TEST"])

    if hasattr(component,"sam") and component.sam:
        file_name = "%s/hadd_%s_%s.root"%(component.output_dir, component.name, component.version)
        if os.path.exists(file_name):
            root_files = [file_name]
        else:
            root_files = []
    else:
        root_files = glob.glob("%s/*.root"%output_dir)
        if verbose:
            print "tools: Found %i files"%len(root_files)
            for root_file in root_files:
                print "tools:     - %s"%root_file
    return root_files
########################### 
# open all histograms
def openPlots(root_file_name, skip_2d=False, directories=False, verbose=False):
    from ROOT import TFile
    master = {}
    if verbose: print "tools: Opening: %s"%root_file_name
    open_root_file = TFile.Open( root_file_name )
    assert open_root_file.IsOpen()
    if directories:
        folder_list = open_root_file.GetListOfKeys()
        folder_list.sort()
        if verbose: print "tools:    Found %i keys"%len(folder_list)
        for folder in folder_list:
            folder_name = folder.GetName()
            directory = open_root_file.GetDirectory(folder_name)
            histogram_list = directory.GetListOfKeys()
            histogram_list.sort()
            for histogram in histogram_list:
                name = histogram.GetName()
                master["%s_%s"%(folder_name,name)] = directory.Get(name)
                assert master["%s_%s"%(folder_name,name)], "tools: Failure opening file: %s, histogram: %s/%s"%(root_name,folder_name,name)
    else:
        histogram_list = open_root_file.GetListOfKeys()
        histogram_list.sort()
        print "run  :         Found %i keys"%len(histogram_list)
        for histogram in histogram_list:
            name = histogram.GetName()
            if (skip_2d and ("_vs_" in name)): continue
            master[name] = open_root_file.Get(name)
            assert master[name], "tools: Failure opening file: %s, histogram: %s"%(root_name,name)
    return master, open_root_file
########################### 
# check if a histogram name follows the logical naming convention
def isLogicalHistogramName(name, strict=False):
    import naming_schema
    if "/" in name:    name = name.split("/")[-1]
    if ".png" in name: name = name[:-4]
    # histograms should follow a logical naming convention:
    #   <category>/<observable>-<reconstruction level>-<truth level>
    if len(name.split("-")) < 2: return False
    if strict:
        category = getHistogramCategory(name)
        if category not in naming_schema.categories.keys():    return False

        observable = getHistogramObservable(name)
        if observable not in naming_schema.observables.keys(): return False

        reco_level = getHistogramRecoLevel(name)
        if reco_level not in naming_schema.reco_levels.keys(): return False

        true_level = getHistogramTrueLevel(name)
        if true_level not in naming_schema.true_levels.keys(): return False

    return True
########################### 
# return the logical category of a histogram
def getHistogramCategory(name):
    if "/" in name:    name = name.split("/")[-1]
    if ".png" in name: name = name[:-4]
    return name.split("-")[0]
########################### 
# return the logical observable of a histogram
def getHistogramObservable(name):
    if "/" in name:    name = name.split("/")[-1]
    if ".png" in name: name = name[:-4]
    return name.split("-")[1]
########################### 
# return the logical reco-level of a histogram
def getHistogramRecoLevel(name):
    import naming_schema
    if "/" in name:    name = name.split("/")[-1]
    if ".png" in name: name = name[:-4]
    name = name.split("-")
    if len(name) == 2:
        # no reco or truth level
        return ""
    if len(name) == 3:
        # could be reco or truth, disambiguate
        if name[2] in naming_schema.true_levels.keys(): return ""
        if not (name[2] in naming_schema.reco_levels.keys()):
            fail("tools: unknown reco level found:%s (from: %s)"%(name[2],name))
        return name[2]
    if len(name) == 4:
        # both reco and true, so return reco
        if not (name[2] in naming_schema.reco_levels.keys()):
            fail("tools: unknown reco level found: %s (from: %s)"%(name[2],name))
            assert False
        return name[2]
    return ""
########################### 
# return the logical true-level of a histogram
def getHistogramTrueLevel(name):
    import naming_schema
    if "/" in name:    name = name.split("/")[-1]
    if ".png" in name: name = name[:-4]
    name = name.split("-")
    if len(name) == 2:
        # no reco or truth level
        return ""
    if len(name) == 3:
        # could be reco or truth, disambiguate
        if name[2] in naming_schema.reco_levels.keys(): return ""
        if not (name[2] in naming_schema.true_levels.keys()):
            fail("tools: unknown true level found: %s (from: %s)"%(name[2],name))
        return name[2]
    if len(name) == 4:
        # both reco and true, so return true
        if not (name[3] in naming_schema.true_levels.keys()):
            fail("tools: unknown true level found: %s (from: %s)"%(name[3],name))
        return name[3]
    return ""
########################### 
# make a directory
def mkdir(directory, verbose=False):
    import os
    if verbose: print "tools: mkdir: %s"%directory
    if (not os.path.exists(directory)):
        os.system("mkdir -p %s"%directory)
########################### 
# Read a bunch of statistics from an array of plots
def getPlotStats(plots, area_norm=True, verbose=True):
    from PlotStatistics import PlotStatistics
    from ROOT import TH1F, TH1D, TGraph
    stats = PlotStatistics()    
    if (type(plots[0]) in [TH1F, TH1D]):
        if area_norm:
            stats.normalisations = [("area","area_norm_"),("none","")]
        else:
            stats.normalisations = [("none","")]
    else:
        stats.normalisations = []
    for i,plot in enumerate(plots):
        if type(plot) in [TGraph]:
            stats.entries.append(-9)
            stats.sum_of_weights.append(-9)
        else:
            stats.entries       .append(plot.GetEntries())
            stats.sum_of_weights.append(plot.GetSumOfWeights())
        if type(plot) in [TH1F,TH1D]:
            stats.means     .append(plot.GetMean())
            stats.medians   .append(getMedian(plot))
            stats.modes     .append(plot.GetXaxis().GetBinCenter(plot.GetMaximumBin()))
            if i == 0:
                stats.mean_diffs.append(0)
            else:
                try:
                    mean_diff = ((plot.GetMean() - stats.means[0]) / stats.means[0]) * 100.
                except:
                    mean_diff = 0.
                stats.mean_diffs.append(mean_diff)
        else:
            stats.means     .append(-9)
            stats.medians   .append(-9)
            stats.modes     .append(-9)
            stats.mean_diffs.append(-9)
    return stats
########################### 
# function for getting the median of a histogram
def getMedian(histogram):
    quantiles = getQuantiles(histogram,2)
    return quantiles[0]
########################### 
# function for dividing histograms into quantiles
def getQuantiles(histogram,n_bins):
    from array import array
    quantiles = array("d",[(a+1.)/n_bins for a in range(n_bins)])
    locations = array("d",[0 for a in range(n_bins)])
    histogram.GetQuantiles(n_bins,locations,quantiles)
    return locations
########################### 
# convert a directory to a URL
def directoryToURL(directory):
    return directory.replace("/nusoft/app/web/htdoc/","http://nusoft.fnal.gov/")
########################### 
# find the maximum and minimum of a list of histograms
def findMinMax(histograms,ignore_zero=False,include_errors=True,verbose=False):
    this_max = 0
    this_min = 10e15
    for i,d in enumerate(histograms):
        name = d.GetName()
        start = d.GetXaxis().GetFirst()
        for bin in range(start, d.GetXaxis().GetLast()+1):
            entries = d.GetBinContent(bin)
            error = 0
            if include_errors: error = d.GetBinError(bin)      

            if (ignore_zero and entries ==0): continue        
            if ((entries) and ((entries-error) < this_min)):
                this_min = entries-error
                if verbose: print 'tools: new minium minimum %.2f, from %s, bin[%i]'%(this_min, name, bin)
            if ((entries) and ((entries+error) > this_max)):
                this_max = entries+error
                if verbose: print 'tools: new maximum maximum %.2f, from %s, bin[%i]'%(this_max, name, bin)
        if ignore_zero and this_min == 0:
            this_min = 0.01
                    
    if this_min == 10e15:
        # print "no no-zero minimum found. Using zero instead, although these plots are probably rubbish"
        this_min = 0
    
    return this_min, this_max
########################### 
# plots can have different truth levels and different reco levels, group them here
def findPlotTypes(plots, verbose=True):
    plots.sort()
    plot_types  = []
    reco_levels = []
    true_levels = []
    for plot in plots:
        category   = getHistogramCategory(plot)
        observable = getHistogramObservable(plot)
        plot_type = "%s-%s"%(category,observable)
        if (plot_type not in plot_types): plot_types.append(plot_type)
        reco_level = getHistogramRecoLevel(plot)
        if (reco_level not in reco_levels) and (category not in ["meta"]): reco_levels.append(reco_level)
        true_level = getHistogramTrueLevel(plot)
        if (true_level not in true_levels) and (category not in ["meta"]): true_levels.append(true_level)
    if verbose:
        print "view :     found %i plot types, %i reco-levels and %i true-levels"%(len(plot_types),len(reco_levels),len(true_levels))
    return plot_types,reco_levels,true_levels
########################### 
# check if these levels are the standard ones
def isStandardTrueLevels(levels):
    import nameing_schema
    standard = set(naming_schema.standard_truth_components)
    return set(levels) == standard
########################### 
# get all of the different histogram observables
def getHistogramObservables(plots):
    observables = []
    for plot in plots:
        observable = getHistogramObservable(plot)
        if observable in observables: continue
        observables.append(observable)
    return observables
########################### 
# get all of the different histogram categories
def getHistogramCategoriesAndObservables(plots):
    plot_types = []
    for plot in plots:
        category   = getHistogramCategory(plot)
        observable = getHistogramObservable(plot)
        plot_type = "%s-%s"%(category,observable)
        if (plot_type not in plot_types): plot_types.append(plot_type)
    return plot_types
########################### 
# Find the maximum skew of a list of histograms
def findMaxSkew(plots,verbose=True):
    if verbose:
        print "vtool: --- find Max Skew"

    max_skew = plots[0].GetSkewness()
    for i,plot in enumerate(plots):
        this_skew = plot.GetSkewness()
        if abs(this_skew) > abs(max_skew):
            max_skew = this_skew
        if verbose:
            print "vtool:       [%i] %s: %f"%(i,plot.GetName(),this_skew)

    if verbose: print "vtool:     Max: %f"%max_skew
    return max_skew
########################### 
# Find the maximum asymmetry of a list of histograms
def findMaxAsymmetry(plots,verbose=False):
    if verbose:
        print "vtool: --- find Max Asymmetry"

    max_asym = findAsymmetry(plots[0])
    for i,plot in enumerate(plots):
        this_asym = findAsymmetry(plot)
        if abs(this_asym) > abs(max_asym):
            max_asym = this_asym
        if verbose:
            print "vtool:       [%i] %s: %f"%(i,plot.GetName(),this_asym)

    if verbose: print "vtool:     Max: %f"%max_asym
    return max_asym
########################### 
# Find the asymmetry of histogram
def findAsymmetry(plot):
    n_bins = plot.GetNbinsX()
    low_int  = plot.Integral(0,             int(n_bins/2))
    high_int = plot.Integral(int(n_bins/2), n_bins)
    try:
        asym = (high_int - low_int) / (low_int + high_int)
    except ZeroDivisionError:
        asym = 0.
    return asym
########################### 
# Sort a list of plots
def sortPlotList(plots):
    import naming_schema
    plots.sort(key=lambda plot: "%s-%s"%(getHistogramCategory(plot),naming_schema.observableName(getHistogramObservable(plot))))
    return plots
########################### 
# check if a validation matches a logical pattern
def filterValidation(filter_pattern, name, version):
    tokens = filter_pattern.split(":")
    if len(tokens) == 1: tokens.append("*")
    assert len(tokens) == 2, "Filter list doesn't match expected pattern (*:*): %s"%filter_pattern
    
    if (tokens[0] != "*") and (tokens[0] != name): return False
    if (tokens[1] != "*") and (tokens[1] != version): return False

    return True