diff --git a/bin/rivet-mkanalysis b/bin/rivet-mkanalysis --- a/bin/rivet-mkanalysis +++ b/bin/rivet-mkanalysis @@ -1,375 +1,381 @@ #! /usr/bin/env python """\ %(prog)s: make templates of analysis source files for Rivet Usage: %(prog)s [--help|-h] [--srcroot=] Without the --srcroot flag, the analysis files will be created in the current directory. """ import rivet, sys, os rivet.util.check_python_version() rivet.util.set_process_name(os.path.basename(__file__)) import logging ## Handle command line import argparse parser = argparse.ArgumentParser(usage=__doc__) parser.add_argument("ANANAMES", nargs="+", help="names of analyses to make") parser.add_argument("--srcroot", metavar="DIR", dest="SRCROOT", default=None, help="install the templates into the Rivet source tree (rooted " + "at directory DIR) rather than just creating all in the current dir") parser.add_argument("-q", "--quiet", dest="LOGLEVEL", default=logging.INFO, action="store_const", const=logging.WARNING, help="only write out warning and error messages") parser.add_argument("-v", "--verbose", dest="LOGLEVEL", default=logging.INFO, action="store_const", const=logging.DEBUG, help="provide extra debugging messages") parser.add_argument("-i", "--inline-info", dest="INLINE", action="store_true", default=False, help="Put analysis info into source file instead of separate data file.") args = parser.parse_args() logging.basicConfig(format="%(msg)s", level=args.LOGLEVEL) ANANAMES = args.ANANAMES ## Work out installation paths ANAROOT = os.path.abspath(args.SRCROOT or os.getcwd()) if not os.access(ANAROOT, os.W_OK): logging.error("Can't write to source root directory %s" % ANAROOT) sys.exit(1) ANASRCDIR = os.getcwd() ANAINFODIR = os.getcwd() ANAPLOTDIR = os.getcwd() if args.SRCROOT: ANASRCDIR = os.path.join(ANAROOT, "src/Analyses") ANAINFODIR = os.path.join(ANAROOT, "data/anainfo") ANAPLOTDIR = os.path.join(ANAROOT, "data/plotinfo") if not (os.path.exists(ANASRCDIR) and os.path.exists(ANAINFODIR) and os.path.exists(ANAPLOTDIR)): logging.error("Rivet analysis dirs do not exist under %s" % ANAROOT) sys.exit(1) if not (os.access(ANASRCDIR, os.W_OK) and os.access(ANAINFODIR, os.W_OK) and os.access(ANAPLOTDIR, os.W_OK)): logging.error("Can't write to Rivet analysis dirs under %s" % ANAROOT) sys.exit(1) ## Check for disallowed characters in analysis names import string allowedchars = string.ascii_letters + string.digits + "_" all_ok = True for ananame in ANANAMES: for c in ananame: if c not in allowedchars: logging.error("Analysis name '%s' contains disallowed character '%s'!" % (ananame, c)) all_ok = False break if not all_ok: logging.error("Exiting... please ensure that all analysis names are valid") sys.exit(1) ## Now make each analysis for ANANAME in ANANAMES: logging.info("Writing templates for %s to %s" % (ANANAME, ANAROOT)) ## Extract some metadata from the name if it matches the standard pattern import re re_stdana = re.compile(r"^(\w+)_(\d{4})_(I|S)(\d+)$") match = re_stdana.match(ANANAME) STDANA = False ANAEXPT = "" ANACOLLIDER = "" ANAYEAR = "" INSPIRE_SPIRES = 'I' ANAINSPIREID = "" if match: STDANA = True ANAEXPT = match.group(1) if ANAEXPT.upper() in ("ALICE", "ATLAS", "CMS", "LHCB"): ANACOLLIDER = "LHC" elif ANAEXPT.upper() in ("CDF", "D0"): ANACOLLIDER = "Tevatron" elif ANAEXPT.upper() == "BABAR": ANACOLLIDER = "PEP-II" elif ANAEXPT.upper() == "BELLE": ANACOLLIDER = "KEKB" ANAYEAR = match.group(2) INSPIRE_SPIRES = match.group(3) ANAINSPIREID = match.group(4) if INSPIRE_SPIRES == "I": ANAREFREPO = "Inspire" else: ANAREFREPO = "Spires" KEYWORDS = { "ANANAME" : ANANAME, "ANAEXPT" : ANAEXPT, "ANACOLLIDER" : ANACOLLIDER, "ANAYEAR" : ANAYEAR, "ANAREFREPO" : ANAREFREPO, "ANAINSPIREID" : ANAINSPIREID } ## Try to get bib info from SPIRES ANABIBKEY = "" ANABIBTEX = "" bibkey, bibtex = None, None if STDANA: try: logging.info("Getting Inspire/SPIRES biblio data for '%s'" % ANANAME) bibkey, bibtex = rivet.spiresbib.get_bibtex_from_repo(INSPIRE_SPIRES, ANAINSPIREID) except Exception as e: logging.error("Inspire/SPIRES oops: %s" % e) if bibkey and bibtex: ANABIBKEY = bibkey ANABIBTEX = bibtex KEYWORDS["ANABIBKEY"] = ANABIBKEY KEYWORDS["ANABIBTEX"] = ANABIBTEX ## Try to download YODA data file from HepData if STDANA: try: try: from urllib.request import urlretrieve except: from urllib import urlretrieve # hdurl = None if INSPIRE_SPIRES == "I": hdurl = "http://www.hepdata.net/record/ins%s?format=yoda&rivet=%s" % (ANAINSPIREID, ANANAME) if not hdurl: raise Exception("Couldn't identify a URL for getting reference data from HepData") logging.info("Getting data file from HepData at %s" % hdurl) tmpfile = urlretrieve(hdurl)[0] # import tarfile, shutil tar = tarfile.open(tmpfile, mode="r") fnames = tar.getnames() if len(fnames) > 1: logging.warning("Found more than one file in downloaded archive. Treating first as canonical") tar.extractall() shutil.move(fnames[0], "%s.yoda" % ANANAME) except Exception as e: logging.warning("Problem encountered retrieving from HepData: %s" % hdurl) logging.warning("No reference data file written") logging.debug("HepData oops: %s: %s" % (str(type(e)), e)) INLINEMETHODS = "" if args.INLINE: KEYWORDS["ANAREFREPO_LOWER"] = KEYWORDS["ANAREFREPO"].lower() INLINEMETHODS = """ public: string experiment() const { return "%(ANAEXPT)s"; } string year() const { return "%(ANAYEAR)s"; } string %(ANAREFREPO_LOWER)sId() const { return "%(ANAINSPIREID)s"; } string collider() const { return ""; } string summary() const { return ""; } string description() const { return ""; } string runInfo() const { return ""; } string bibKey() const { return "%(ANABIBKEY)s"; } string bibTeX() const { return "%(ANABIBTEX)s"; } string status() const { return "UNVALIDATED"; } vector authors() const { return vector(); } vector references() const { return vector(); } vector todos() const { return vector(); } """ % KEYWORDS del KEYWORDS["ANAREFREPO_LOWER"] KEYWORDS["INLINEMETHODS"] = INLINEMETHODS if ANANAME.startswith("MC_"): HISTOBOOKING = """\ book(_h["XXXX"], "myh1", 20, 0.0, 100.0); book(_h["YYYY"], "myh2", logspace(20, 1e-2, 1e3)); book(_h["ZZZZ"], "myh3", {0.0, 1.0, 2.0, 4.0, 8.0, 16.0}); book(_p["AAAA"], "myp", 20, 0.0, 100.0); book(_c["BBBB"], "myc");""" % KEYWORDS else: HISTOBOOKING = """\ // specify custom binning book(_h["XXXX"], "myh1", 20, 0.0, 100.0); book(_h["YYYY"], "myh2", logspace(20, 1e-2, 1e3)); book(_h["ZZZZ"], "myh3", {0.0, 1.0, 2.0, 4.0, 8.0, 16.0}); // take binning from reference data using HEPData ID (digits in "d01-x01-y01" etc.) book(_h["AAAA"], 1, 1, 1); book(_p["BBBB"], 2, 1, 1); book(_c["CCCC"], 3, 1, 1);""" % KEYWORDS KEYWORDS["HISTOBOOKING"] = HISTOBOOKING ANASRCFILE = os.path.join(ANASRCDIR, ANANAME+".cc") logging.debug("Writing implementation template to %s" % ANASRCFILE) f = open(ANASRCFILE, "w") src = '''\ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/MissingMomentum.hh" #include "Rivet/Projections/PromptFinalState.hh" namespace Rivet { /// @brief Add a short analysis description here class %(ANANAME)s : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(%(ANANAME)s); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // Initialise and register projections - // the basic final-state projection: - // all final-state particles within + // the basic final-state projection: + // all final-state particles within // the given eta acceptance const FinalState fs(Cuts::abseta < 4.9); // the final-state particles declared above are clustered using FastJet with // the anti-kT algorithm and a jet-radius parameter 0.4 // muons and neutrinos are excluded from the clustering FastJets jetfs(fs, FastJets::ANTIKT, 0.4, JetAlg::Muons::NONE, JetAlg::Invisibles::NONE); declare(jetfs, "jets"); // FinalState of prompt photons and bare muons and electrons in the event PromptFinalState photons(Cuts::abspid == PID::PHOTON); PromptFinalState bare_leps(Cuts::abspid == PID::MUON || Cuts::abspid == PID::ELECTRON); // dress the prompt bare leptons with prompt photons within dR < 0.1 // apply some fiducial cuts on the dressed leptons Cut lepton_cuts = Cuts::abseta < 2.5 && Cuts::pT > 20*GeV; DressedLeptons dressed_leps(photons, bare_leps, 0.1, lepton_cuts); declare(dressed_leps, "leptons"); // missing momentum declare(MissingMomentum(fs), "MET"); // Book histograms %(HISTOBOOKING)s } /// Perform the per-event analysis void analyze(const Event& event) { /// @todo Do the event by event analysis here // retrieve dressed leptons, sorted by pT vector leptons = apply(event, "leptons").dressedLeptons(); // retrieve clustered jets, sorted by pT, with a minimum pT cut Jets jets = apply(event, "jets").jetsByPt(Cuts::pT > 30*GeV); // remove all jets within dR < 0.2 of a dressed lepton idiscardIfAnyDeltaRLess(jets, leptons, 0.2); // select jets ghost-associated to B-hadrons with a certain fiducial selection Jets bjets = filter_select(jets, [](const Jet& jet) { return jet.bTagged(Cuts::pT > 5*GeV && Cuts::abseta < 2.5); }); // veto event if there are no b-jets if (bjets.empty()) vetoEvent; // apply a missing-momentum cut if (apply(event, "MET").missingPt() < 30*GeV) vetoEvent; // fill histogram with leading b-jet pT _h["XXXX"]->fill(bjets[0].pT()/GeV); } /// Normalise histograms etc., after the run void finalize() { normalize(_h["YYYY"]); // normalize to unity scale(_h["ZZZZ"], crossSection()/picobarn/sumOfWeights()); // norm to cross section } //@} /// @name Histograms //@{ map _h; map _p; map _c; //@} %(INLINEMETHODS)s }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(%(ANANAME)s); } ''' % KEYWORDS f.write(src) f.close() ANAPLOTFILE = os.path.join(ANAPLOTDIR, ANANAME+".plot") logging.debug("Writing plot template to %s" % ANAPLOTFILE) f = open(ANAPLOTFILE, "w") src = '''\ BEGIN PLOT /%(ANANAME)s/d01-x01-y01 Title=[Insert title for histogram d01-x01-y01 here] XLabel=[Insert $x$-axis label for histogram d01-x01-y01 here] YLabel=[Insert $y$-axis label for histogram d01-x01-y01 here] # + any additional plot settings you might like, see make-plots documentation END PLOT # ... add more histograms as you need them ... ''' % KEYWORDS f.write(src) f.close() if args.INLINE: sys.exit(0) ANAINFOFILE = os.path.join(ANAINFODIR, ANANAME+".info") logging.debug("Writing info template to %s" % ANAINFOFILE) f = open(ANAINFOFILE, "w") src = """\ Name: %(ANANAME)s Year: %(ANAYEAR)s -Summary: +Summary: Experiment: %(ANAEXPT)s Collider: %(ANACOLLIDER)s %(ANAREFREPO)sID: %(ANAINSPIREID)s -Status: UNVALIDATED +Status: UNVALIDATED NOTREENTRY NOHEPDATA SINGLEWEIGHT UNPHYSICAL Authors: - Your Name #References: #- '' #- '' #- '' RunInfo: #Beams: #Energies: #Luminosity_fb: Description: - ' 50\;\GeV$.>' + 'A brief description of what is measured and what it is useful for + in terms of MC testing, tuning, reinterpretation, etc. Use LaTeX + for maths, e.g. $\pT > 50\;\GeV$.' +ValidationInfo: + 'A description of the process used to validate the Rivet code against + the original experimental analysis. Cut-flow tables and similar information + are welcome' +#ReleaseTests: +# - $A my-hepmc-prefix :MODE=some_rivet_flag Keywords: [] BibKey: %(ANABIBKEY)s BibTeX: '%(ANABIBTEX)s' ToDo: - Implement the analysis, test it, remove this ToDo, and mark as VALIDATED :-) """ % KEYWORDS f.write(src) f.close() logging.info("Use e.g. 'rivet-buildplugin Rivet%s.so %s.cc' to compile the plugin" % (ANANAME, ANANAME))