%s

Generated at %s

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

%s

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

%s

\n\n' % opts.TITLE) if opts.SINGLE: ## Write table of contents 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 opts.IGNORE_UNVALIDATED and ana.status().upper() != "VALIDATED": continue if opts.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"))) 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('

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

\n' % description) except UnicodeEncodeError as ue: print("Unicode error in analysis description for " + analysis + ": " + str(ue)) anapath = os.path.join(opts.OUTPUTDIR, analysis) if not opts.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(opts.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+"."+opts.VECTORFORMAT.lower() srcfile = obsname+".dat" if opts.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 opts.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 opts.NUMTHREADS: mp_cmd.append("--num-threads=%d" % opts.NUMTHREADS) if opts.NO_CLEANUP: mp_cmd.append("--no-cleanup") if opts.NO_SUBPROC: mp_cmd.append("--no-subproc") if opts.VECTORFORMAT == "PDF": mp_cmd.append("--pdfpng") elif opts.VECTORFORMAT == "PS": mp_cmd.append("--pspng") if opts.OUTPUT_FONT: mp_cmd.append("--font=%s" % opts.OUTPUT_FONT) # if opts.OUTPUT_FONT.upper() == "PALATINO": # mp_cmd.append("--palatino") # if opts.OUTPUT_FONT.upper() == "CM": # mp_cmd.append("--cm") # elif opts.OUTPUT_FONT.upper() == "TIMES": # mp_cmd.append("--times") # elif opts.OUTPUT_FONT.upper() == "HELVETICA": # mp_cmd.append("--helvetica") # elif opts.OUTPUT_FONT.upper() == "MINION": # mp_cmd.append("--minion") datfiles = [] for analysis in analyses: anapath = os.path.join(opts.OUTPUTDIR, analysis) #print(anapath) anadatfiles = glob.glob("%s/*.dat" % anapath) datfiles += sorted(anadatfiles) if datfiles: mp_cmd += datfiles if opts.VERBOSE: mp_cmd.append("--verbose") print("Calling make-plots with the following options:") print(" ".join(mp_cmd)) if not opts.DRY_RUN: Popen(mp_cmd).wait() if opts.BOOKLET and opts.VECTORFORMAT=="PDF": if which("pdftk") is not None: bookletcmd = ["pdftk"] for analysis in analyses: anapath = os.path.join(opts.OUTPUTDIR, analysis) bookletcmd += sorted(glob.glob("%s/*.pdf" % anapath)) bookletcmd += ["cat", "output", "%s/booklet.pdf" % opts.OUTPUTDIR] print(bookletcmd) Popen(bookletcmd).wait() elif which("pdfmerge") is not None: bookletcmd = ["pdfmerge"] for analysis in analyses: anapath = os.path.join(opts.OUTPUTDIR, analysis) bookletcmd += sorted(glob.glob("%s/*.pdf" % anapath)) bookletcmd += ["%s/booklet.pdf" % opts.OUTPUTDIR] print(bookletcmd) Popen(bookletcmd).wait() else: print("Neither pdftk nor pdfmerge available --- not booklet output possible") diff --git a/include/Rivet/Analyses/MC_Cent_pPb.hh b/include/Rivet/Analyses/MC_Cent_pPb.hh new file mode 100644 --- /dev/null +++ b/include/Rivet/Analyses/MC_Cent_pPb.hh @@ -0,0 +1,87 @@ +// -*- C++ -*- +#ifndef RIVET_MC_Cent_pPb_HH +#define RIVET_MC_Cent_pPb_HH + +#include "Rivet/Projections/FinalState.hh" +#include "Rivet/Projections/ChargedFinalState.hh" +#include "Rivet/Projections/SingleValueProjection.hh" +#include "Rivet/Projections/TriggerProjection.hh" + +namespace Rivet { + +/// Example of a centrality observable projection for pPb that uses +/// summed Et in the Pb direction. +class MC_SumETFwdPbCentrality: public SingleValueProjection { + +public: + + /// Constructor. + MC_SumETFwdPbCentrality() { + declare(FinalState(Cuts::eta < -3.2 && Cuts::eta > -4.9 && Cuts::pT > 0.1*GeV), + "FSSumETFwdCentrality"); + } + + /// Clone on the heap. + DEFAULT_RIVET_PROJ_CLONE(MC_SumETFwdPbCentrality); + +protected: + + /// Perform the projection on the Event + void project(const Event& e) { + clear(); + const FinalState & fsfwd = + apply(e, "FSSumETFwdCentrality"); + double estimate = 0.0; + for ( const Particle & p : fsfwd.particles() ) { + estimate += p.Et(); + } + set(estimate); + } + + /// Compare projections + int compare(const Projection& p) const { + return mkNamedPCmp(p, "FSSumETFwdCentrality"); + } + +}; + +/// Example of a trigger projection for minimum bias pPb requiring at +/// least one charged particle in both forward and backward direction. +class MC_pPbMinBiasTrigger: public TriggerProjection { + +public: + + /// Constructor. + MC_pPbMinBiasTrigger() { + declare(FinalState(Cuts::eta < -3.2 && Cuts::eta > -4.9 && Cuts::pT > 0.1*GeV), + "FSSumETFwdCentrality"); + declare(ChargedFinalState(Cuts::eta > 2.09 && Cuts::eta < 3.84 && + Cuts::pT > 0.1*GeV), "MBB"); + declare(ChargedFinalState(Cuts::eta < -2.09 && Cuts::eta > -3.84 && + Cuts::pT > 0.1*GeV), "MBF"); + } + + /// Clone on the heap. + DEFAULT_RIVET_PROJ_CLONE(MC_pPbMinBiasTrigger); + +protected: + + /// Perform the projection on the Event + void project(const Event& event) { + pass(); + if ( applyProjection(event,"MBF").particles().empty() || + applyProjection(event,"MBB").particles().empty() ) + fail(); + } + + /// Compare projections + int compare(const Projection& p) const { + return mkNamedPCmp(p, "MBF") || mkNamedPCmp(p, "MBB"); + } + +}; + + +} + +#endif diff --git a/include/Rivet/Analysis.hh b/include/Rivet/Analysis.hh --- a/include/Rivet/Analysis.hh +++ b/include/Rivet/Analysis.hh @@ -1,1150 +1,1226 @@ // -*- C++ -*- #ifndef RIVET_Analysis_HH #define RIVET_Analysis_HH #include "Rivet/Config/RivetCommon.hh" #include "Rivet/AnalysisInfo.hh" #include "Rivet/Event.hh" #include "Rivet/Projection.hh" #include "Rivet/ProjectionApplier.hh" #include "Rivet/ProjectionHandler.hh" #include "Rivet/AnalysisLoader.hh" #include "Rivet/Tools/Cuts.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Tools/ParticleUtils.hh" #include "Rivet/Tools/BinnedHistogram.hh" #include "Rivet/Tools/RivetMT2.hh" #include "Rivet/Tools/RivetYODA.hh" +#include "Rivet/Tools/Percentile.hh" #include "Rivet/Projections/CentralityProjection.hh" /// @def vetoEvent /// Preprocessor define for vetoing events, including the log message and return. #define vetoEvent \ do { MSG_DEBUG("Vetoing event on line " << __LINE__ << " of " << __FILE__); return; } while(0) namespace Rivet { // Forward declaration class AnalysisHandler; /// @brief This is the base class of all analysis classes in Rivet. /// /// There are /// three virtual functions which should be implemented in base classes: /// /// void init() is called by Rivet before a run is started. Here the /// analysis class should book necessary histograms. The needed /// projections should probably rather be constructed in the /// constructor. /// /// void analyze(const Event&) is called once for each event. Here the /// analysis class should apply the necessary Projections and fill the /// histograms. /// /// void finalize() is called after a run is finished. Here the analysis /// class should do whatever manipulations are necessary on the /// histograms. Writing the histograms to a file is, however, done by /// the Rivet class. class Analysis : public ProjectionApplier { /// The AnalysisHandler is a friend. friend class AnalysisHandler; public: /// @name Standard constructors and destructors. //@{ // /// The default constructor. // Analysis(); /// Constructor Analysis(const std::string& name); /// The destructor. virtual ~Analysis() {} //@} public: /// @name Main analysis methods //@{ /// Initialize this analysis object. A concrete class should here /// book all necessary histograms. An overridden function must make /// sure it first calls the base class function. virtual void init() { } /// Analyze one event. A concrete class should here apply the /// necessary projections on the \a event and fill the relevant /// histograms. An overridden function must make sure it first calls /// the base class function. virtual void analyze(const Event& event) = 0; /// Finalize this analysis object. A concrete class should here make /// all necessary operations on the histograms. Writing the /// histograms to a file is, however, done by the Rivet class. An /// overridden function must make sure it first calls the base class /// function. virtual void finalize() { } //@} public: /// @name Metadata /// Metadata is used for querying from the command line and also for /// building web pages and the analysis pages in the Rivet manual. //@{ /// Get the actual AnalysisInfo object in which all this metadata is stored. const AnalysisInfo& info() const { assert(_info && "No AnalysisInfo object :O"); return *_info; } /// @brief Get the name of the analysis. /// /// By default this is computed by combining the results of the /// experiment, year and Spires ID metadata methods and you should /// only override it if there's a good reason why those won't /// work. If options has been set for this instance, a /// corresponding string is appended at the end. virtual std::string name() const { return ( (info().name().empty()) ? _defaultname : info().name() ) + _optstring; } // get name of reference data file, which could be different from plugin name virtual std::string getRefDataName() const { return (info().getRefDataName().empty()) ? _defaultname : info().getRefDataName(); } // set name of reference data file, which could be different from plugin name virtual void setRefDataName(const std::string& ref_data="") { info().setRefDataName(!ref_data.empty() ? ref_data : name()); } /// Get the Inspire ID code for this analysis. virtual std::string inspireId() const { return info().inspireId(); } /// Get the SPIRES ID code for this analysis (~deprecated). virtual std::string spiresId() const { return info().spiresId(); } /// @brief Names & emails of paper/analysis authors. /// /// Names and email of authors in 'NAME \' format. The first /// name in the list should be the primary contact person. virtual std::vector authors() const { return info().authors(); } /// @brief Get a short description of the analysis. /// /// Short (one sentence) description used as an index entry. /// Use @a description() to provide full descriptive paragraphs /// of analysis details. virtual std::string summary() const { return info().summary(); } /// @brief Get a full description of the analysis. /// /// Full textual description of this analysis, what it is useful for, /// what experimental techniques are applied, etc. Should be treated /// as a chunk of restructuredText (http://docutils.sourceforge.net/rst.html), /// with equations to be rendered as LaTeX with amsmath operators. virtual std::string description() const { return info().description(); } /// @brief Information about the events needed as input for this analysis. /// /// Event types, energies, kinematic cuts, particles to be considered /// stable, etc. etc. Should be treated as a restructuredText bullet list /// (http://docutils.sourceforge.net/rst.html) virtual std::string runInfo() const { return info().runInfo(); } /// Experiment which performed and published this analysis. virtual std::string experiment() const { return info().experiment(); } /// Collider on which the experiment ran. virtual std::string collider() const { return info().collider(); } /// When the original experimental analysis was published. virtual std::string year() const { return info().year(); } /// The luminosity in inverse femtobarn virtual std::string luminosityfb() const { return info().luminosityfb(); } /// Journal, and preprint references. virtual std::vector references() const { return info().references(); } /// BibTeX citation key for this article. virtual std::string bibKey() const { return info().bibKey(); } /// BibTeX citation entry for this article. virtual std::string bibTeX() const { return info().bibTeX(); } /// Whether this analysis is trusted (in any way!) virtual std::string status() const { return (info().status().empty()) ? "UNVALIDATED" : info().status(); } /// Any work to be done on this analysis. virtual std::vector todos() const { return info().todos(); } /// Return the allowed pairs of incoming beams required by this analysis. virtual const std::vector& requiredBeams() const { return info().beams(); } /// Declare the allowed pairs of incoming beams required by this analysis. virtual Analysis& setRequiredBeams(const std::vector& requiredBeams) { info().setBeams(requiredBeams); return *this; } /// Sets of valid beam energy pairs, in GeV virtual const std::vector >& requiredEnergies() const { return info().energies(); } /// Get vector of analysis keywords virtual const std::vector & keywords() const { return info().keywords(); } /// Declare the list of valid beam energy pairs, in GeV virtual Analysis& setRequiredEnergies(const std::vector >& requiredEnergies) { info().setEnergies(requiredEnergies); return *this; } /// Return true if this analysis needs to know the process cross-section. /// @todo Remove this and require HepMC >= 2.06 bool needsCrossSection() const { return info().needsCrossSection(); } /// Declare whether this analysis needs to know the process cross-section from the generator. /// @todo Remove this and require HepMC >= 2.06 Analysis& setNeedsCrossSection(bool needed=true) { info().setNeedsCrossSection(needed); return *this; } //@} /// @name Internal metadata modifying methods //@{ /// Get the actual AnalysisInfo object in which all this metadata is stored (non-const). AnalysisInfo& info() { assert(_info && "No AnalysisInfo object :O"); return *_info; } //@} /// @name Run conditions //@{ /// Incoming beams for this run const ParticlePair& beams() const; /// Incoming beam IDs for this run const PdgIdPair beamIds() const; /// Centre of mass energy for this run double sqrtS() const; //@} /// @name Analysis / beam compatibility testing //@{ /// Check if analysis is compatible with the provided beam particle IDs and energies bool isCompatible(const ParticlePair& beams) const; /// Check if analysis is compatible with the provided beam particle IDs and energies bool isCompatible(PdgId beam1, PdgId beam2, double e1, double e2) const; /// Check if analysis is compatible with the provided beam particle IDs and energies bool isCompatible(const PdgIdPair& beams, const std::pair& energies) const; //@} /// Set the cross section from the generator Analysis& setCrossSection(double xs); //, double xserr=0.0); /// Access the controlling AnalysisHandler object. AnalysisHandler& handler() const { return *_analysishandler; } protected: /// Get a Log object based on the name() property of the calling analysis object. Log& getLog() const; /// Get the process cross-section in pb. Throws if this hasn't been set. double crossSection() const; /// Get the process cross-section per generated event in pb. Throws if this /// hasn't been set. double crossSectionPerEvent() const; /// @brief Get the number of events seen (via the analysis handler). /// /// @note Use in the finalize phase only. size_t numEvents() const; /// @brief Get the sum of event weights seen (via the analysis handler). /// /// @note Use in the finalize phase only. double sumW() const; /// Alias double sumOfWeights() const { return sumW(); } /// @brief Get the sum of squared event weights seen (via the analysis handler). /// /// @note Use in the finalize phase only. double sumW2() const; protected: /// @name Histogram paths //@{ /// Get the canonical histogram "directory" path for this analysis. const std::string histoDir() const; /// Get the canonical histogram path for the named histogram in this analysis. const std::string histoPath(const std::string& hname) const; /// Get the canonical histogram path for the numbered histogram in this analysis. const std::string histoPath(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const; /// Get the internal histogram name for given d, x and y (cf. HepData) const std::string mkAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const; /// Alias /// @deprecated Prefer the "mk" form, consistent with other "making function" names const std::string makeAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { return mkAxisCode(datasetId, xAxisId, yAxisId); } //@} /// @name Histogram reference data //@{ /// Get reference data for a named histo /// @todo SFINAE to ensure that the type inherits from YODA::AnalysisObject? template const T& refData(const string& hname) const { _cacheRefData(); MSG_TRACE("Using histo bin edges for " << name() << ":" << hname); if (!_refdata[hname]) { MSG_ERROR("Can't find reference histogram " << hname); throw Exception("Reference data " + hname + " not found."); } return dynamic_cast(*_refdata[hname]); } /// Get reference data for a numbered histo /// @todo SFINAE to ensure that the type inherits from YODA::AnalysisObject? template const T& refData(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { const string hname = makeAxisCode(datasetId, xAxisId, yAxisId); return refData(hname); } //@} /// @name Counter booking //@{ /// Book a counter. CounterPtr bookCounter(const std::string& name, const std::string& title=""); // const std::string& valtitle="" /// Book a counter, using a path generated from the dataset and axis ID codes /// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. CounterPtr bookCounter(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const std::string& title=""); // const std::string& valtitle="" //@} /// @name 1D histogram booking //@{ /// Book a 1D histogram with @a nbins uniformly distributed across the range @a lower - @a upper . Histo1DPtr bookHisto1D(const std::string& name, size_t nbins, double lower, double upper, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D histogram with non-uniform bins defined by the vector of bin edges @a binedges . Histo1DPtr bookHisto1D(const std::string& name, const std::vector& binedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D histogram with non-uniform bins defined by the vector of bin edges @a binedges . Histo1DPtr bookHisto1D(const std::string& name, const std::initializer_list& binedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D histogram with binning from a reference scatter. Histo1DPtr bookHisto1D(const std::string& name, const Scatter2D& refscatter, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D histogram, using the binnings in the reference data histogram. Histo1DPtr bookHisto1D(const std::string& name, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D histogram, using the binnings in the reference data histogram. /// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. Histo1DPtr bookHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); //@} /// @name 2D histogram booking //@{ /// Book a 2D histogram with @a nxbins and @a nybins uniformly /// distributed across the ranges @a xlower - @a xupper and @a /// ylower - @a yupper respectively along the x- and y-axis. Histo2DPtr bookHisto2D(const std::string& name, size_t nxbins, double xlower, double xupper, size_t nybins, double ylower, double yupper, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D histogram with non-uniform bins defined by the /// vectors of bin edges @a xbinedges and @a ybinedges. Histo2DPtr bookHisto2D(const std::string& name, const std::vector& xbinedges, const std::vector& ybinedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D histogram with non-uniform bins defined by the /// vectors of bin edges @a xbinedges and @a ybinedges. Histo2DPtr bookHisto2D(const std::string& name, const std::initializer_list& xbinedges, const std::initializer_list& ybinedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D histogram with binning from a reference scatter. Histo2DPtr bookHisto2D(const std::string& name, const Scatter3D& refscatter, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D histogram, using the binnings in the reference data histogram. Histo2DPtr bookHisto2D(const std::string& name, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D histogram, using the binnings in the reference data histogram. /// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. Histo2DPtr bookHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); //@} /// @name 1D profile histogram booking //@{ /// Book a 1D profile histogram with @a nbins uniformly distributed across the range @a lower - @a upper . Profile1DPtr bookProfile1D(const std::string& name, size_t nbins, double lower, double upper, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D profile histogram with non-uniform bins defined by the vector of bin edges @a binedges . Profile1DPtr bookProfile1D(const std::string& name, const std::vector& binedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D profile histogram with non-uniform bins defined by the vector of bin edges @a binedges . Profile1DPtr bookProfile1D(const std::string& name, const std::initializer_list& binedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D profile histogram with binning from a reference scatter. Profile1DPtr bookProfile1D(const std::string& name, const Scatter2D& refscatter, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D profile histogram, using the binnings in the reference data histogram. Profile1DPtr bookProfile1D(const std::string& name, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D profile histogram, using the binnings in the reference data histogram. /// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. Profile1DPtr bookProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); //@} /// @name 2D profile histogram booking //@{ /// Book a 2D profile histogram with @a nxbins and @a nybins uniformly /// distributed across the ranges @a xlower - @a xupper and @a ylower - @a /// yupper respectively along the x- and y-axis. Profile2DPtr bookProfile2D(const std::string& name, size_t nxbins, double xlower, double xupper, size_t nybins, double ylower, double yupper, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D profile histogram with non-uniform bins defined by the vectorx /// of bin edges @a xbinedges and @a ybinedges. Profile2DPtr bookProfile2D(const std::string& name, const std::vector& xbinedges, const std::vector& ybinedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D profile histogram with non-uniform bins defined by the vectorx /// of bin edges @a xbinedges and @a ybinedges. Profile2DPtr bookProfile2D(const std::string& name, const std::initializer_list& xbinedges, const std::initializer_list& ybinedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D profile histogram with binning from a reference scatter. Profile2DPtr bookProfile2D(const std::string& name, const Scatter3D& refscatter, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D profile histogram, using the binnings in the reference data histogram. Profile2DPtr bookProfile2D(const std::string& name, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D profile histogram, using the binnings in the reference data histogram. /// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. Profile2DPtr bookProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); //@} /// @name 2D scatter booking //@{ /// @brief Book a 2-dimensional data point set with the given name. /// /// @note Unlike histogram booking, scatter booking by default makes no /// attempt to use reference data to pre-fill the data object. If you want /// this, which is sometimes useful e.g. when the x-position is not really /// meaningful and can't be extracted from the data, then set the @a /// copy_pts parameter to true. This creates points to match the reference /// data's x values and errors, but with the y values and errors zeroed... /// assuming that there is a reference histo with the same name: if there /// isn't, an exception will be thrown. Scatter2DPtr bookScatter2D(const std::string& name, bool copy_pts=false, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// @brief Book a 2-dimensional data point set, using the binnings in the reference data histogram. /// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. /// /// @note Unlike histogram booking, scatter booking by default makes no /// attempt to use reference data to pre-fill the data object. If you want /// this, which is sometimes useful e.g. when the x-position is not really /// meaningful and can't be extracted from the data, then set the @a /// copy_pts parameter to true. This creates points to match the reference /// data's x values and errors, but with the y values and errors zeroed. Scatter2DPtr bookScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, bool copy_pts=false, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// @brief Book a 2-dimensional data point set with equally spaced x-points in a range. /// /// The y values and errors will be set to 0. Scatter2DPtr bookScatter2D(const std::string& name, size_t npts, double lower, double upper, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// @brief Book a 2-dimensional data point set based on provided contiguous "bin edges". /// /// The y values and errors will be set to 0. Scatter2DPtr bookScatter2D(const std::string& hname, const std::vector& binedges, const std::string& title, const std::string& xtitle, const std::string& ytitle); //@} public: /// @name accessing options for this Analysis instance. //@{ /// Get an option for this analysis instance as a string. std::string getOption(std::string optname) { if ( _options.find(optname) != _options.end() ) return _options.find(optname)->second; return ""; } /// Get an option for this analysis instance converted to a /// specific type (given by the specified @a def value). template T getOption(std::string optname, T def) { if (_options.find(optname) == _options.end()) return def; std::stringstream ss; ss << _options.find(optname)->second; T ret; ss >> ret; return ret; } //@} /// @brief Book a CentralityProjection /// /// Using a SingleValueProjection, @a proj, giving the value of an /// experimental observable to be used as a centrality estimator, /// book a CentralityProjection based on the experimentally /// measured pecentiles of this observable (as given by the /// reference data for the @a calHistName histogram in the @a /// calAnaName analysis. If a preloaded file with the output of a /// run using the @a calAnaName analysis contains a valid /// generated @a calHistName histogram, it will be used as an /// optional percentile binning. Also if this preloaded file /// contains a histogram with the name @a calHistName with an /// appended "_IMP" This histogram will be used to add an optional /// centrality percentile based on the generated impact /// parameter. If @increasing is true, a low (high) value of @proj /// is assumed to correspond to a more peripheral (central) event. const CentralityProjection& declareCentrality(const SingleValueProjection &proj, string calAnaName, string calHistName, const string projName, bool increasing = false); + /// @brief Book a Pecentile wrapper around AnalysisObjects. + /// + /// Based on a previously registered CentralityProjection named @a + /// projName book one AnalysisObject for each @a centralityBin and + /// name them according to the corresponding code in the @a ref + /// vector. + template + Percentile bookPercentile(string projName, + vector > centralityBins, + vector > ref) { + typedef typename ReferenceTraits::RefT RefT; + Percentile pctl(this, projName); + + const int nCent = centralityBins.size(); + for (int iCent = 0; iCent < nCent; ++iCent) { + const string axisCode = makeAxisCode(std::get<0>(ref[iCent]), + std::get<1>(ref[iCent]), + std::get<2>(ref[iCent])); + shared_ptr ao; + CounterPtr cnt; + try { + ao = getAnalysisObject(axisCode); + MSG_TRACE("Found old " << histoPath(axisCode)); + } + catch (Exception) { + const RefT & refscatter = refData(axisCode); + ao = make_shared(refscatter, histoPath(axisCode)); + addAnalysisObject(ao); + MSG_TRACE("Created new " << histoPath(axisCode)); + } + try { + cnt = getAnalysisObject("TMP/COUNTER/" + axisCode); + MSG_TRACE("Found old " << histoPath("TMP/COUNTER/" + axisCode)); + } + catch (Exception) { + cnt = make_shared(histoPath("TMP/COUNTER/" + axisCode)); + addAnalysisObject(cnt); + MSG_TRACE("Created new " << histoPath("TMP/COUNTER/" + axisCode)); + } + pctl.add(ao, cnt, centralityBins[iCent]); + } + return pctl; + } + + /// @brief Book Pecentile wrappers around AnalysisObjects. + /// + /// Based on a previously registered CentralityProjection named @a + /// projName book one (or several) AnalysisObject(s) named + /// according to @a ref where the x-axis will be filled according + /// to the percentile output(s) of the @projName. + template + PercentileXaxis bookPercentileXaxis(string projName, + tuple ref) { + typedef typename ReferenceTraits::RefT RefT; + PercentileXaxis pctl(this, projName); + + const string axisCode = makeAxisCode(std::get<0>(ref), + std::get<1>(ref), + std::get<2>(ref)); + shared_ptr ao; + CounterPtr cnt; + try { + ao = getAnalysisObject(histoPath(axisCode)); + } + catch (Exception) { + const RefT & refscatter = refData(axisCode); + ao = make_shared(refscatter, axisCode); + addAnalysisObject(ao); + } + pctl.add(proj, ao, make_shared()); + + return pctl; + } + + /// @name Analysis object manipulation /// @todo Should really be protected: only public to keep BinnedHistogram happy for now... //@{ /// Multiplicatively scale the given counter, @a cnt, by factor @s factor. void scale(CounterPtr cnt, double factor); /// Multiplicatively scale the given counters, @a cnts, by factor @s factor. /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh) /// @todo Use SFINAE for a generic iterable of CounterPtrs void scale(const std::vector& cnts, double factor) { for (auto& c : cnts) scale(c, factor); } /// @todo YUCK! template void scale(const CounterPtr (&cnts)[array_size], double factor) { // for (size_t i = 0; i < std::extent::value; ++i) scale(cnts[i], factor); for (auto& c : cnts) scale(c, factor); } /// Normalize the given histogram, @a histo, to area = @a norm. void normalize(Histo1DPtr histo, double norm=1.0, bool includeoverflows=true); /// Normalize the given histograms, @a histos, to area = @a norm. /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh) /// @todo Use SFINAE for a generic iterable of Histo1DPtrs void normalize(const std::vector& histos, double norm=1.0, bool includeoverflows=true) { for (auto& h : histos) normalize(h, norm, includeoverflows); } /// @todo YUCK! template void normalize(const Histo1DPtr (&histos)[array_size], double norm=1.0, bool includeoverflows=true) { for (auto& h : histos) normalize(h, norm, includeoverflows); } /// Multiplicatively scale the given histogram, @a histo, by factor @s factor. void scale(Histo1DPtr histo, double factor); /// Multiplicatively scale the given histograms, @a histos, by factor @s factor. /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh) /// @todo Use SFINAE for a generic iterable of Histo1DPtrs void scale(const std::vector& histos, double factor) { for (auto& h : histos) scale(h, factor); } /// @todo YUCK! template void scale(const Histo1DPtr (&histos)[array_size], double factor) { for (auto& h : histos) scale(h, factor); } /// Normalize the given histogram, @a histo, to area = @a norm. void normalize(Histo2DPtr histo, double norm=1.0, bool includeoverflows=true); /// Normalize the given histograms, @a histos, to area = @a norm. /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh) /// @todo Use SFINAE for a generic iterable of Histo2DPtrs void normalize(const std::vector& histos, double norm=1.0, bool includeoverflows=true) { for (auto& h : histos) normalize(h, norm, includeoverflows); } /// @todo YUCK! template void normalize(const Histo2DPtr (&histos)[array_size], double norm=1.0, bool includeoverflows=true) { for (auto& h : histos) normalize(h, norm, includeoverflows); } /// Multiplicatively scale the given histogram, @a histo, by factor @s factor. void scale(Histo2DPtr histo, double factor); /// Multiplicatively scale the given histograms, @a histos, by factor @s factor. /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh) /// @todo Use SFINAE for a generic iterable of Histo2DPtrs void scale(const std::vector& histos, double factor) { for (auto& h : histos) scale(h, factor); } /// @todo YUCK! template void scale(const Histo2DPtr (&histos)[array_size], double factor) { for (auto& h : histos) scale(h, factor); } /// Helper for counter division. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(CounterPtr c1, CounterPtr c2, Scatter1DPtr s) const; /// Helper for histogram division with raw YODA objects. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(const YODA::Counter& c1, const YODA::Counter& c2, Scatter1DPtr s) const; /// Helper for histogram division. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const; /// Helper for histogram division with raw YODA objects. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const; /// Helper for profile histogram division. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(Profile1DPtr p1, Profile1DPtr p2, Scatter2DPtr s) const; /// Helper for profile histogram division with raw YODA objects. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(const YODA::Profile1D& p1, const YODA::Profile1D& p2, Scatter2DPtr s) const; /// Helper for 2D histogram division. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(Histo2DPtr h1, Histo2DPtr h2, Scatter3DPtr s) const; /// Helper for 2D histogram division with raw YODA objects. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(const YODA::Histo2D& h1, const YODA::Histo2D& h2, Scatter3DPtr s) const; /// Helper for 2D profile histogram division. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(Profile2DPtr p1, Profile2DPtr p2, Scatter3DPtr s) const; /// Helper for 2D profile histogram division with raw YODA objects /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(const YODA::Profile2D& p1, const YODA::Profile2D& p2, Scatter3DPtr s) const; /// Helper for histogram efficiency calculation. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void efficiency(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const; /// Helper for histogram efficiency calculation. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void efficiency(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const; /// Helper for histogram asymmetry calculation. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void asymm(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const; /// Helper for histogram asymmetry calculation. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void asymm(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const; /// Helper for converting a differential histo to an integral one. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void integrate(Histo1DPtr h, Scatter2DPtr s) const; /// Helper for converting a differential histo to an integral one. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void integrate(const Histo1D& h, Scatter2DPtr s) const; //@} public: /// List of registered analysis data objects const vector& analysisObjects() const { return _analysisobjects; } protected: /// @name Data object registration, retrieval, and removal //@{ /// Register a data object in the histogram system void addAnalysisObject(AnalysisObjectPtr ao); /// Get a data object from the histogram system template const std::shared_ptr getAnalysisObject(const std::string& name) const { foreach (const AnalysisObjectPtr& ao, analysisObjects()) { if (ao->path() == histoPath(name)) return dynamic_pointer_cast(ao); } throw Exception("Data object " + histoPath(name) + " not found"); } /// Get a data object from the histogram system (non-const) template std::shared_ptr getAnalysisObject(const std::string& name) { foreach (const AnalysisObjectPtr& ao, analysisObjects()) { if (ao->path() == histoPath(name)) return dynamic_pointer_cast(ao); } throw Exception("Data object " + histoPath(name) + " not found"); } /// Unregister a data object from the histogram system (by name) void removeAnalysisObject(const std::string& path); /// Unregister a data object from the histogram system (by pointer) void removeAnalysisObject(AnalysisObjectPtr ao); /// Get all data object from the AnalysisHandler. vector getAllData(bool includeorphans) const; /// Get a data object from another analysis (e.g. preloaded /// calibration histogram). /// Get a data object from the histogram system (non-const) template std::shared_ptr getAnalysisObject(const std::string & ananame, const std::string& name) { std::string path = "/" + ananame + "/" + name; for ( AnalysisObjectPtr ao : getAllData(true) ) { if ( ao->path() == path ) return dynamic_pointer_cast(ao); } return std::shared_ptr(); } /// Get a named Histo1D object from the histogram system const Histo1DPtr getHisto1D(const std::string& name) const { return getAnalysisObject(name); } /// Get a named Histo1D object from the histogram system (non-const) Histo1DPtr getHisto1D(const std::string& name) { return getAnalysisObject(name); } /// Get a Histo1D object from the histogram system by axis ID codes (non-const) const Histo1DPtr getHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } /// Get a Histo1D object from the histogram system by axis ID codes (non-const) Histo1DPtr getHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } /// Get a named Histo2D object from the histogram system const Histo2DPtr getHisto2D(const std::string& name) const { return getAnalysisObject(name); } /// Get a named Histo2D object from the histogram system (non-const) Histo2DPtr getHisto2D(const std::string& name) { return getAnalysisObject(name); } /// Get a Histo2D object from the histogram system by axis ID codes (non-const) const Histo2DPtr getHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } /// Get a Histo2D object from the histogram system by axis ID codes (non-const) Histo2DPtr getHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } /// Get a named Profile1D object from the histogram system const Profile1DPtr getProfile1D(const std::string& name) const { return getAnalysisObject(name); } /// Get a named Profile1D object from the histogram system (non-const) Profile1DPtr getProfile1D(const std::string& name) { return getAnalysisObject(name); } /// Get a Profile1D object from the histogram system by axis ID codes (non-const) const Profile1DPtr getProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } /// Get a Profile1D object from the histogram system by axis ID codes (non-const) Profile1DPtr getProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } /// Get a named Profile2D object from the histogram system const Profile2DPtr getProfile2D(const std::string& name) const { return getAnalysisObject(name); } /// Get a named Profile2D object from the histogram system (non-const) Profile2DPtr getProfile2D(const std::string& name) { return getAnalysisObject(name); } /// Get a Profile2D object from the histogram system by axis ID codes (non-const) const Profile2DPtr getProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } /// Get a Profile2D object from the histogram system by axis ID codes (non-const) Profile2DPtr getProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } /// Get a named Scatter2D object from the histogram system const Scatter2DPtr getScatter2D(const std::string& name) const { return getAnalysisObject(name); } /// Get a named Scatter2D object from the histogram system (non-const) Scatter2DPtr getScatter2D(const std::string& name) { return getAnalysisObject(name); } /// Get a Scatter2D object from the histogram system by axis ID codes (non-const) const Scatter2DPtr getScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } /// Get a Scatter2D object from the histogram system by axis ID codes (non-const) Scatter2DPtr getScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) { return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); } //@} private: /// Name passed to constructor (used to find .info analysis data file, and as a fallback) string _defaultname; /// Pointer to analysis metadata object unique_ptr _info; /// Storage of all plot objects /// @todo Make this a map for fast lookup by path? vector _analysisobjects; /// @name Cross-section variables //@{ double _crossSection; bool _gotCrossSection; //@} /// The controlling AnalysisHandler object. AnalysisHandler* _analysishandler; /// Collection of cached refdata to speed up many autobookings: the /// reference data file should only be read once. mutable std::map _refdata; /// Options the (this instance of) the analysis map _options; /// The string of options. string _optstring; private: /// @name Utility functions //@{ /// Get the reference data for this paper and cache it. void _cacheRefData() const; //@} /// The assignment operator is private and must never be called. /// In fact, it should not even be implemented. Analysis& operator=(const Analysis&); }; } // Include definition of analysis plugin system so that analyses automatically see it when including Analysis.hh #include "Rivet/AnalysisBuilder.hh" /// @def DECLARE_RIVET_PLUGIN /// Preprocessor define to prettify the global-object plugin hook mechanism. #define DECLARE_RIVET_PLUGIN(clsname) Rivet::AnalysisBuilder plugin_ ## clsname /// @def DECLARE_ALIASED_RIVET_PLUGIN /// Preprocessor define to prettify the global-object plugin hook mechanism, with an extra alias name for this analysis. // #define DECLARE_ALIASED_RIVET_PLUGIN(clsname, alias) Rivet::AnalysisBuilder plugin_ ## clsname ## ( ## #alias ## ) #define DECLARE_ALIASED_RIVET_PLUGIN(clsname, alias) DECLARE_RIVET_PLUGIN(clsname)( #alias ) /// @def DEFAULT_RIVET_ANALYSIS_CONSTRUCTOR /// Preprocessor define to prettify the manky constructor with name string argument #define DEFAULT_RIVET_ANALYSIS_CONSTRUCTOR(clsname) clsname() : Analysis(# clsname) {} /// @def DEFAULT_RIVET_ANALYSIS_CTOR /// Slight abbreviation for DEFAULT_RIVET_ANALYSIS_CONSTRUCTOR #define DEFAULT_RIVET_ANALYSIS_CTOR(clsname) DEFAULT_RIVET_ANALYSIS_CONSTRUCTOR(clsname) #endif diff --git a/include/Rivet/AnalysisHandler.hh b/include/Rivet/AnalysisHandler.hh --- a/include/Rivet/AnalysisHandler.hh +++ b/include/Rivet/AnalysisHandler.hh @@ -1,272 +1,308 @@ // -*- C++ -*- #ifndef RIVET_RivetHandler_HH #define RIVET_RivetHandler_HH #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Particle.hh" #include "Rivet/AnalysisLoader.hh" #include "Rivet/Tools/RivetYODA.hh" namespace Rivet { // Forward declaration and smart pointer for Analysis class Analysis; typedef std::shared_ptr AnaHandle; // Needed to make smart pointers compare equivalent in the STL set struct CmpAnaHandle { bool operator() (const AnaHandle& a, const AnaHandle& b) const { return a.get() < b.get(); } }; /// A class which handles a number of analysis objects to be applied to /// generated events. An {@link Analysis}' AnalysisHandler is also responsible /// for handling the final writing-out of histograms. class AnalysisHandler { public: /// @name Constructors and destructors. */ //@{ /// Preferred constructor, with optional run name. AnalysisHandler(const string& runname=""); /// @brief Destructor /// The destructor is not virtual, as this class should not be inherited from. ~AnalysisHandler(); //@} private: /// Get a logger object. Log& getLog() const; public: /// @name Run properties //@{ /// Get the name of this run. string runName() const { return _runname; } /// Get the number of events seen. Should only really be used by external /// steering code or analyses in the finalize phase. size_t numEvents() const { return _eventcounter.numEntries(); } /// @brief Access the sum of the event weights seen /// /// This is the weighted equivalent of the number of events. It should only /// be used by external steering code or analyses in the finalize phase. double sumW() const { return _eventcounter.sumW(); } /// Access to the sum of squared-weights double sumW2() const { return _eventcounter.sumW2(); } /// @brief Compatibility alias for sumOfWeights /// /// @deprecated Prefer sumW double sumOfWeights() const { return sumW(); } /// @brief Set the sum of weights /// /// This is useful if Rivet is steered externally and /// the analyses are run for a sub-contribution of the events /// (but of course have to be normalised to the total sum of weights) /// /// @todo What about the sumW2 term? That needs to be set coherently. Need a /// new version, with all three N,sumW,sumW2 numbers (or a counter) /// supplied. /// /// @deprecated Weight sums are no longer tracked this way... void setSumOfWeights(const double& sum) { //_sumOfWeights = sum; throw Error("Can't set sum of weights independently, since it's now tracked by a Counter. " "Please contact the Rivet authors if you need this."); } /// Is cross-section information required by at least one child analysis? /// @deprecated Should no-longer be an issue: does any generator not write the cross-section? bool needCrossSection() const; /// Whether the handler knows about a cross-section. /// @deprecated Should no-longer be an issue: does any generator not write the cross-section? bool hasCrossSection() const; /// Get the cross-section known to the handler. double crossSection() const { return _xs; } /// Set the cross-section for the process being generated. /// @todo What about the xsec uncertainty? Add a second, optional argument? AnalysisHandler& setCrossSection(double xs); /// Set the beam particles for this run AnalysisHandler& setRunBeams(const ParticlePair& beams) { _beams = beams; MSG_DEBUG("Setting run beams = " << beams << " @ " << sqrtS()/GeV << " GeV"); return *this; } /// Get the beam particles for this run, usually determined from the first event. const ParticlePair& beams() const { return _beams; } /// Get beam IDs for this run, usually determined from the first event. /// @deprecated Use standalone beamIds(ah.beams()), to clean AH interface PdgIdPair beamIds() const; /// Get energy for this run, usually determined from the first event. /// @deprecated Use standalone sqrtS(ah.beams()), to clean AH interface double sqrtS() const; /// Setter for _ignoreBeams void setIgnoreBeams(bool ignore=true); //@} /// @name Handle analyses //@{ /// Get a list of the currently registered analyses' names. std::vector analysisNames() const; /// Get the collection of currently registered analyses. const std::set& analyses() const { return _analyses; } /// Get a registered analysis by name. const AnaHandle analysis(const std::string& analysisname) const; /// Add an analysis to the run list by object AnalysisHandler& addAnalysis(Analysis* analysis); /// @brief Add an analysis to the run list using its name. /// /// The actual Analysis to be used will be obtained via /// AnalysisLoader::getAnalysis(string). If no matching analysis is found, /// no analysis is added (i.e. the null pointer is checked and discarded. AnalysisHandler& addAnalysis(const std::string& analysisname); /// @brief Add an analysis with a map of analysis options. AnalysisHandler& addAnalysis(const std::string& analysisname, std::map pars); /// @brief Add analyses to the run list using their names. /// /// The actual {@link Analysis}' to be used will be obtained via /// AnalysisHandler::addAnalysis(string), which in turn uses /// AnalysisLoader::getAnalysis(string). If no matching analysis is found /// for a given name, no analysis is added, but also no error is thrown. AnalysisHandler& addAnalyses(const std::vector& analysisnames); /// Remove an analysis from the run list using its name. AnalysisHandler& removeAnalysis(const std::string& analysisname); /// Remove analyses from the run list using their names. AnalysisHandler& removeAnalyses(const std::vector& analysisnames); //@} /// @name Main init/execute/finalise //@{ /// Initialize a run, with the run beams taken from the example event. void init(const GenEvent& event); /// @brief Analyze the given \a event by reference. /// /// This function will call the AnalysisBase::analyze() function of all /// included analysis objects. void analyze(const GenEvent& event); /// @brief Analyze the given \a event by pointer. /// /// This function will call the AnalysisBase::analyze() function of all /// included analysis objects, after checking the event pointer validity. void analyze(const GenEvent* event); /// Finalize a run. This function calls the AnalysisBase::finalize() /// functions of all included analysis objects. void finalize(); //@} /// @name Histogram / data object access //@{ /// Add a vector of analysis objects to the current state. void addData(const std::vector& aos); /// Read analysis plots into the histo collection (via addData) from the named file. void readData(const std::string& filename); /// Get all analyses' plots as a vector of analysis objects. - std::vector getData(bool includeorphans = false) const; + std::vector getData(bool includeorphans = false, + bool includetmps = false) const; /// Write all analyses' plots (via getData) to the named file. void writeData(const std::string& filename) const; + /// Tell Rivet to dump intermediate result to a file named @a + /// dumpfile every @a period'th event. If @period is not positive, + /// no dumping will be done. + void dump(string dumpfile, int period) { + _dumpPeriod = period; + _dumpFile = dumpfile; + } + + /// Take the vector of yoda files and merge them together using + /// the cross section and weight information provided in each + /// file. Each file in @a aofiles is assumed to have been produced + /// by Rivet. By default the files are assumed to contain + /// different processes (or the same processs but mutually + /// exclusive cuts), but if @a equiv if ture, the files are + /// assumed to contain output of completely equivalent (but + /// statistically independent) Rivet runs. The corresponding + /// analyses will be loaded and their analysis objects will be + /// filled with the merged result. finalize() will be run on each + /// relevant anslysis. The resulting YODA file can then be rwitten + /// out by writeData(). + void mergeYodas(const vector & aofiles, bool equiv = false); + //@} - private: /// The collection of Analysis objects to be used. set _analyses; /// A vector of pre-loaded object which do not have a valid /// Analysis plugged in. vector _orphanedPreloads; + /// A vector containing copies of analysis objects after + /// finalize() has been run. + vector _finalizedAOs; /// @name Run properties //@{ /// Run name std::string _runname; /// Event counter Counter _eventcounter; /// Cross-section known to AH double _xs, _xserr; /// Beams used by this run. ParticlePair _beams; /// Flag to check if init has been called bool _initialised; /// Flag whether input event beams should be ignored in compatibility check bool _ignoreBeams; + /// Determines how often Rivet runs finalize() and writes the + /// result to a YODA file. + int _dumpPeriod; + + /// The name of a YODA file to which Rivet periodically dumps + /// results. + string _dumpFile; + + /// Flag to indicate periodic dumping is in progress + bool _dumping; + //@} private: /// The assignment operator is private and must never be called. /// In fact, it should not even be implemented. AnalysisHandler& operator=(const AnalysisHandler&); /// The copy constructor is private and must never be called. In /// fact, it should not even be implemented. AnalysisHandler(const AnalysisHandler&); }; } #endif diff --git a/include/Rivet/AnalysisInfo.hh b/include/Rivet/AnalysisInfo.hh --- a/include/Rivet/AnalysisInfo.hh +++ b/include/Rivet/AnalysisInfo.hh @@ -1,289 +1,298 @@ // -*- C++ -*- #ifndef RIVET_AnalysisInfo_HH #define RIVET_AnalysisInfo_HH #include "Rivet/Config/RivetCommon.hh" #include namespace Rivet { class AnalysisInfo { public: /// Static factory method: returns null pointer if no metadata found static unique_ptr make(const std::string& name); /// @name Standard constructors and destructors. //@{ /// The default constructor. AnalysisInfo() { clear(); } /// The destructor. ~AnalysisInfo() { } //@} public: /// @name Metadata /// Metadata is used for querying from the command line and also for /// building web pages and the analysis pages in the Rivet manual. //@{ /// Get the name of the analysis. By default this is computed using the /// experiment, year and Inspire/Spires ID metadata methods. std::string name() const { if (!_name.empty()) return _name; if (!experiment().empty() && !year().empty()) { if (!inspireId().empty()) { return experiment() + "_" + year() + "_I" + inspireId(); } else if (!spiresId().empty()) { return experiment() + "_" + year() + "_S" + spiresId(); } } return ""; } /// Set the name of the analysis. void setName(const std::string& name) { _name = name; } /// Get the reference data name of the analysis (if different from plugin name). std::string getRefDataName() const { if (!_refDataName.empty()) return _refDataName; return name(); } /// Set the reference data name of the analysis (if different from plugin name). void setRefDataName(const std::string& name) { _refDataName = name; } /// Get the Inspire (SPIRES replacement) ID code for this analysis. const std::string& inspireId() const { return _inspireId; } /// Set the Inspire (SPIRES replacement) ID code for this analysis. void setInspireId(const std::string& inspireId) { _inspireId = inspireId; } /// Get the SPIRES ID code for this analysis. const std::string& spiresId() const { return _spiresId; } /// Set the SPIRES ID code for this analysis. void setSpiresId(const std::string& spiresId) { _spiresId = spiresId; } /// @brief Names & emails of paper/analysis authors. /// Names and email of authors in 'NAME \' format. The first /// name in the list should be the primary contact person. const std::vector& authors() const { return _authors; } /// Set the author list. void setAuthors(const std::vector& authors) { _authors = authors; } /// @brief Get a short description of the analysis. /// Short (one sentence) description used as an index entry. /// Use @a description() to provide full descriptive paragraphs /// of analysis details. const std::string& summary() const { return _summary; } /// Set the short description for this analysis. void setSummary(const std::string& summary) { _summary = summary; } /// @brief Get a full description of the analysis. /// Full textual description of this analysis, what it is useful for, /// what experimental techniques are applied, etc. Should be treated /// as a chunk of restructuredText (http://docutils.sourceforge.net/rst.html), /// with equations to be rendered as LaTeX with amsmath operators. const std::string& description() const { return _description; } /// Set the full description for this analysis. void setDescription(const std::string& description) { _description = description; } /// @brief Information about the events needed as input for this analysis. /// Event types, energies, kinematic cuts, particles to be considered /// stable, etc. etc. Should be treated as a restructuredText bullet list /// (http://docutils.sourceforge.net/rst.html) const std::string& runInfo() const { return _runInfo; } /// Set the full description for this analysis. void setRunInfo(const std::string& runInfo) { _runInfo = runInfo; } /// Beam particle types const std::vector& beams() const { return _beams; } /// Set beam particle types void setBeams(const std::vector& beams) { _beams = beams; } /// Sets of valid beam energies const std::vector >& energies() const { return _energies; } /// Set the valid beam energies void setEnergies(const std::vector >& energies) { _energies = energies; } /// Experiment which performed and published this analysis. const std::string& experiment() const { return _experiment; } /// Set the experiment which performed and published this analysis. void setExperiment(const std::string& experiment) { _experiment = experiment; } /// Collider on which the experiment ran. const std::string& collider() const { return _collider; } /// Set the collider on which the experiment ran. void setCollider(const std::string& collider) { _collider = collider; } /// @brief When the original experimental analysis was published. /// When the refereed paper on which this is based was published, /// according to SPIRES. const std::string& year() const { return _year; } /// Set the year in which the original experimental analysis was published. void setYear(const std::string& year) { _year = year; } /// The integrated data luminosity of the data set const std::string& luminosityfb() const { return _luminosityfb; } /// Set the integrated data luminosity of the data set void setLuminosityfb(const std::string& luminosityfb) { _luminosityfb = luminosityfb; } /// Journal and preprint references. const std::vector& references() const { return _references; } /// Set the journal and preprint reference list. void setReferences(const std::vector& references) { _references = references; } /// Analysis Keywords for grouping etc const std::vector& keywords() const { return _keywords; } /// BibTeX citation key for this article. const std::string& bibKey() const { return _bibKey;} /// Set the BibTeX citation key for this article. void setBibKey(const std::string& bibKey) { _bibKey = bibKey; } /// BibTeX citation entry for this article. const std::string& bibTeX() const { return _bibTeX; } /// Set the BibTeX citation entry for this article. void setBibTeX(const std::string& bibTeX) { _bibTeX = bibTeX; } /// Whether this analysis is trusted (in any way!) const std::string& status() const { return _status; } /// Set the analysis code status. void setStatus(const std::string& status) { _status = status; } /// Any work to be done on this analysis. const std::vector& todos() const { return _todos; } /// Set the to-do list. void setTodos(const std::vector& todos) { _todos = todos; } /// Get the option list. const std::vector& options() const { return _options; } /// Check if the given option is valid. bool validOption(std::string key, std::string val) const; /// Set the option list. void setOptions(const std::vector& opts) { _options = opts; buildOptionMap(); } /// Build a map of options to facilitate checking. void buildOptionMap(); /// Return true if this analysis needs to know the process cross-section. bool needsCrossSection() const { return _needsCrossSection; } /// Return true if this analysis needs to know the process cross-section. void setNeedsCrossSection(bool needXsec) { _needsCrossSection = needXsec; } + /// Return true if finalize() can be run multiple times for this analysis. + bool reentrant() const { return _reentrant; } + + /// setReentrant + void setReentrant(bool ree = true) { _reentrant = ree; } + //@} private: std::string _name; std::string _refDataName; std::string _spiresId, _inspireId; std::vector _authors; std::string _summary; std::string _description; std::string _runInfo; std::string _experiment; std::string _collider; std::vector > _beams; std::vector > _energies; std::string _year; std::string _luminosityfb; std::vector _references; std::vector _keywords; std::string _bibKey; std::string _bibTeX; //std::string _bibTeXBody; ///< Was thinking of avoiding duplication of BibKey... std::string _status; std::vector _todos; bool _needsCrossSection; std::vector _options; std::map< std::string, std::set > _optionmap; + bool _reentrant; + void clear() { _name = ""; _refDataName = ""; _spiresId = ""; _inspireId = ""; _authors.clear(); _summary = ""; _description = ""; _runInfo = ""; _experiment = ""; _collider = ""; _beams.clear(); _energies.clear(); _year = ""; _luminosityfb = ""; _references.clear(); _keywords.clear(); _bibKey = ""; _bibTeX = ""; //_bibTeXBody = ""; _status = ""; _todos.clear(); _needsCrossSection = false; _options.clear(); _optionmap.clear(); + _reentrant = false; } }; /// String representation std::string toString(const AnalysisInfo& ai); /// Stream an AnalysisInfo as a text description inline std::ostream& operator<<(std::ostream& os, const AnalysisInfo& ai) { os << toString(ai); return os; } } #endif diff --git a/include/Rivet/Makefile.am b/include/Rivet/Makefile.am --- a/include/Rivet/Makefile.am +++ b/include/Rivet/Makefile.am @@ -1,181 +1,185 @@ ## Internal headers - not to be installed nobase_dist_noinst_HEADERS = ## Public headers - to be installed nobase_pkginclude_HEADERS = ## Rivet interface nobase_pkginclude_HEADERS += \ Rivet.hh \ Run.hh \ Event.hh \ ParticleBase.hh \ Particle.fhh Particle.hh \ Jet.fhh Jet.hh \ Projection.fhh Projection.hh \ ProjectionApplier.hh \ ProjectionHandler.hh \ Analysis.hh \ AnalysisHandler.hh \ AnalysisInfo.hh \ AnalysisBuilder.hh \ AnalysisLoader.hh ## Build config stuff nobase_pkginclude_HEADERS += \ Config/RivetCommon.hh \ Config/RivetConfig.hh \ Config/BuildOptions.hh ## Projections nobase_pkginclude_HEADERS += \ Projections/AliceCommon.hh \ Projections/AxesDefinition.hh \ Projections/Beam.hh \ Projections/BeamThrust.hh \ Projections/CentralEtHCM.hh \ Projections/CentralityProjection.hh \ Projections/ChargedFinalState.hh \ Projections/ChargedLeptons.hh \ Projections/ConstLossyFinalState.hh \ Projections/DirectFinalState.hh \ Projections/DISFinalState.hh \ Projections/DISKinematics.hh \ Projections/DISLepton.hh \ Projections/DressedLeptons.hh \ Projections/FastJets.hh \ Projections/FinalPartons.hh \ Projections/FinalState.hh \ Projections/FoxWolframMoments.hh \ Projections/FParameter.hh \ Projections/GeneratedPercentileProjection.hh \ Projections/HadronicFinalState.hh \ Projections/HeavyHadrons.hh \ Projections/Hemispheres.hh \ Projections/IdentifiedFinalState.hh \ Projections/ImpactParameterProjection.hh \ Projections/IndirectFinalState.hh \ Projections/InitialQuarks.hh \ Projections/InvMassFinalState.hh \ Projections/JetAlg.hh \ Projections/JetShape.hh \ Projections/LeadingParticlesFinalState.hh \ Projections/LossyFinalState.hh \ Projections/MergedFinalState.hh \ Projections/MissingMomentum.hh \ + Projections/MixedFinalState.hh \ Projections/NeutralFinalState.hh \ Projections/NonHadronicFinalState.hh \ Projections/NonPromptFinalState.hh \ Projections/ParisiTensor.hh \ Projections/ParticleFinder.hh \ Projections/PartonicTops.hh \ Projections/PercentileProjection.hh \ Projections/PrimaryHadrons.hh \ Projections/PrimaryParticles.hh \ Projections/PromptFinalState.hh \ Projections/SingleValueProjection.hh \ Projections/SmearedParticles.hh \ Projections/SmearedJets.hh \ Projections/SmearedMET.hh \ Projections/Sphericity.hh \ Projections/Spherocity.hh \ Projections/TauFinder.hh \ Projections/Thrust.hh \ Projections/TriggerCDFRun0Run1.hh \ Projections/TriggerCDFRun2.hh \ Projections/TriggerProjection.hh \ Projections/TriggerUA5.hh \ Projections/UnstableFinalState.hh \ Projections/UserCentEstimate.hh \ Projections/VetoedFinalState.hh \ Projections/VisibleFinalState.hh \ Projections/WFinder.hh \ Projections/ZFinder.hh ## Meta-projection convenience headers nobase_pkginclude_HEADERS += \ Projections/FinalStates.hh \ Projections/Smearing.hh ## Analysis base class headers # TODO: Move to Rivet/AnalysisTools header dir? nobase_pkginclude_HEADERS += \ + Analyses/MC_Cent_pPb.hh \ Analyses/MC_ParticleAnalysis.hh \ Analyses/MC_JetAnalysis.hh \ Analyses/MC_JetSplittings.hh ## Tools nobase_pkginclude_HEADERS += \ Tools/AliceCommon.hh \ + Tools/AtlasCommon.hh \ Tools/BeamConstraint.hh \ Tools/BinnedHistogram.hh \ Tools/CentralityBinner.hh \ Tools/Cmp.fhh \ Tools/Cmp.hh \ Tools/Cutflow.hh \ Tools/Cuts.fhh \ Tools/Cuts.hh \ Tools/Exceptions.hh \ Tools/JetUtils.hh \ Tools/Logging.hh \ Tools/Random.hh \ Tools/ParticleBaseUtils.hh \ Tools/ParticleIdUtils.hh \ Tools/ParticleUtils.hh \ Tools/ParticleName.hh \ + Tools/Percentile.hh \ Tools/PrettyPrint.hh \ Tools/RivetPaths.hh \ Tools/RivetSTL.hh \ Tools/RivetFastJet.hh \ Tools/RivetHepMC.hh \ Tools/RivetYODA.hh \ Tools/RivetMT2.hh \ Tools/SmearingFunctions.hh \ Tools/MomentumSmearingFunctions.hh \ Tools/ParticleSmearingFunctions.hh \ Tools/JetSmearingFunctions.hh \ Tools/TypeTraits.hh \ Tools/Utils.hh \ Tools/Nsubjettiness/AxesDefinition.hh \ Tools/Nsubjettiness/ExtraRecombiners.hh \ Tools/Nsubjettiness/MeasureDefinition.hh \ Tools/Nsubjettiness/Njettiness.hh \ Tools/Nsubjettiness/NjettinessPlugin.hh \ Tools/Nsubjettiness/Nsubjettiness.hh \ Tools/Nsubjettiness/TauComponents.hh \ Tools/Nsubjettiness/XConePlugin.hh nobase_dist_noinst_HEADERS += \ Tools/osdir.hh ## Maths nobase_pkginclude_HEADERS += \ Math/Matrices.hh \ Math/Vector3.hh \ Math/VectorN.hh \ Math/MatrixN.hh \ Math/MatrixDiag.hh \ Math/MathHeader.hh \ Math/Vectors.hh \ Math/LorentzTrans.hh \ Math/Matrix3.hh \ Math/MathUtils.hh \ Math/Vector4.hh \ Math/Math.hh \ Math/Units.hh \ Math/Constants.hh \ Math/eigen/util.h \ Math/eigen/regressioninternal.h \ Math/eigen/regression.h \ Math/eigen/vector.h \ Math/eigen/ludecompositionbase.h \ Math/eigen/ludecomposition.h \ Math/eigen/linearsolver.h \ Math/eigen/linearsolverbase.h \ Math/eigen/matrix.h \ Math/eigen/vectorbase.h \ Math/eigen/projective.h \ Math/eigen/matrixbase.h diff --git a/include/Rivet/Projections/MixedFinalState.hh b/include/Rivet/Projections/MixedFinalState.hh new file mode 100644 --- /dev/null +++ b/include/Rivet/Projections/MixedFinalState.hh @@ -0,0 +1,97 @@ +// -*- C++ -*- +#ifndef RIVET_MixedFinalState_HH +#define RIVET_MixedFinalState_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 MixedFinalState : public Projection { + public: + // Constructor + MixedFinalState(const ParticleFinder& fsp, const ParticleFinder& mix, size_t nMixIn, + double oMin, double oMax, double deltao ) : nMix(nMixIn){ + setName("MixedFinalState"); + addProjection(fsp,"FS"); + addProjection(mix,"MIX"); + + // 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(MixedFinalState); + + + // 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){ + 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 { + return mkNamedPCmp(p,"FS"); + } + + private: + // 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/include/Rivet/Tools/AtlasCommon.hh b/include/Rivet/Tools/AtlasCommon.hh new file mode 100644 --- /dev/null +++ b/include/Rivet/Tools/AtlasCommon.hh @@ -0,0 +1,85 @@ +// -*- C++ -*- +#ifndef RIVET_ATLAS_COMMON_HH +#define RIVET_ATLAS_COMMON_HH + +#include "Rivet/Projections/FinalState.hh" +#include "Rivet/Projections/ChargedFinalState.hh" +#include "Rivet/Projections/SingleValueProjection.hh" +#include "Rivet/Projections/TriggerProjection.hh" + +namespace Rivet { +/// Common projections for ATLAS trigger conditions and centrality. + + namespace ATLAS { +/// Centrality projection for pPb collisions (one sided) +class SumET_PB_Centrality: public SingleValueProjection { + +public: + + /// Constructor. + SumET_PB_Centrality() { + declare(FinalState(Cuts::eta < -3.2 && Cuts::eta > -4.9 && Cuts::pT > 0.1*GeV), + "SumET_PB_Centrality"); + } + + /// Clone on the heap. + DEFAULT_RIVET_PROJ_CLONE(SumET_PB_Centrality); + +protected: + + /// Perform the projection on the Event + void project(const Event& e) { + clear(); + const FinalState & fsfwd = + apply(e, "SumET_PB_Centrality"); + double estimate = 0.0; + for ( const Particle & p : fsfwd.particles() ) { + estimate += p.Et(); + } + set(estimate); + } + + /// Compare projections + int compare(const Projection& p) const { + return mkNamedPCmp(p, "SumET_PB_Centrality"); + } + +}; + +/// ATLAS min bias trigger conditions. +class MinBiasTrigger: public TriggerProjection { + +public: + + /// Constructor. + MinBiasTrigger() { + declare(ChargedFinalState(Cuts::eta > 2.09 && Cuts::eta < 3.84 && + Cuts::pT > 0.1*GeV), "MBB"); + declare(ChargedFinalState(Cuts::eta < -2.09 && Cuts::eta > -3.84 && + Cuts::pT > 0.1*GeV), "MBF"); + } + + /// Clone on the heap. + DEFAULT_RIVET_PROJ_CLONE(MinBiasTrigger); + +protected: + + /// Perform the projection on the Event + void project(const Event& event) { + pass(); + if ( applyProjection(event,"MBF").particles().empty() || + applyProjection(event,"MBB").particles().empty() ) + fail(); + } + + /// Compare projections + int compare(const Projection& p) const { + return mkNamedPCmp(p, "MBF") || mkNamedPCmp(p, "MBB"); + } + +}; + +} +} + +#endif diff --git a/include/Rivet/Tools/Percentile.hh b/include/Rivet/Tools/Percentile.hh new file mode 100644 --- /dev/null +++ b/include/Rivet/Tools/Percentile.hh @@ -0,0 +1,687 @@ +#ifndef PERCENTILE_HH +#define PERCENTILE_HH + +#include "Rivet/Event.hh" +#include "Rivet/Projections/CentralityProjection.hh" +#include "Rivet/ProjectionApplier.hh" + +namespace Rivet { + + +/// Forward declaration. +class Analysis; + +/// @brief PercentileBase is the base class of all Percentile classes. +/// +/// This base class contains all non-templated variables and +/// infrastructure needed. +class PercentileBase { + +public: + + /// @brief the main constructor + /// + /// requiring a pointer, @a ana, to the Analysis to which this + /// object belongs and the name of the CentralityProjection, @a + /// projname, to be used. + PercentileBase(Analysis * ana, string projName) + : _ana(ana), _projName(projName) {} + + /// @brief Default constructor. + PercentileBase() {} + + /// @initialize the PercentileBase for a new event. + /// + /// This will perform the assigned CentralityProjection and select + /// out the (indices) of the internal AnalysisObjects that are to be + /// active in this event. + void selectBins(const Event &); + + /// @brief Helper function to check if @a x is within @a range. + static bool inRange(double x, pair range) { + return x >= range.first && ( x < range.second || ( x == 100.0 && x == range.second ) ); + } + /// @brief Copy information from @a other PercentileBase + void copyFrom(const PercentileBase & other) { + _ana = other._ana; + _projName = other._projName; + _cent = other._cent; + } + + /// @brief check if @a other PercentileBase is compatible with this. + bool compatible(const PercentileBase & other) const { + return ( _ana == other._ana && + _projName == other._projName && + _cent == other._cent ); + } + + /// @breif return the list of centrality bins. + /// + /// The size of this vector is the same as number of internal + /// analysis objects in the sub class PercentileTBase. + const vector< pair > & centralities() const { + return _cent; + } + +protected: + + /// The Analysis object to which This object is assigned. + Analysis * _ana; + + /// The name of the CentralityProjection. + string _projName; + + /// The list of indices of the analysis objects that are to be + /// filled in the current event. + vector _activeBins; + + /// The list of centrality intervals, one for each included analysis + /// object. + vector > _cent; + +}; + +/// @brief PercentileTBase is the base class of all Percentile classes. +/// +/// This base class contains all template-dependent variables and +/// infrastructure needed for Percentile and PercentileXaxis. +template +class PercentileTBase : public PercentileBase { + +public: + + /// Convenient typedef. + typedef typename T::Ptr TPtr; + + /// @brief the main constructor + /// + /// requiring a pointer, @a ana, to the Analysis to which this + /// object belongs and the name of the CentralityProjection, @a + /// projname, to be used. + PercentileTBase(Analysis * ana, string projName) + : PercentileBase(ana, projName), _histos() {} + + /// @brief Default constructor. + PercentileTBase() {} + + /// @brief Empty destructor. + ~PercentileTBase() {} + + /// @brief add a new percentile bin. + /// + /// Add an analysis objects which are clones of @a temp that should + /// be active for events in the given centrality bin @a + /// cent. Several analysis objects may be added depending on the + /// number of alternative centrality definitions in the + /// CentralityProjection @a proj. This function is common for + /// Percentile and PecentileXaxis, but for the latter the @a cent + /// argument should be left to its default. + void add(shared_ptr ao, CounterPtr cnt, + pair cent = {0.0, 100.0} ) { + _cent.push_back(cent); + _histos.push_back( { ao, cnt } ); + } + + /// @brief Copy the information from an @a other Percentile object. + /// + /// This function differs from a simple assignement as the @a other + /// analysis objects are not copied, but supplied separately through + /// @a tv. + bool add(const PercentileBase & other, const vector & tv) { + copyFrom(other); + if ( tv.size() != _cent.size() ) return false; + for ( auto t : tv ) + _histos.push_back( { t, make_shared() } ); + return true; + } + + /// @brief initialize for a new event. Select which AnalysisObjects + /// should be filled for this event. Keeps track of the number of + /// events seen for each centrality bin and AnalysisAbject. + void init(const Event & event) { + selectBins(event); + for (const auto bin : _activeBins) + _histos[bin].second->fill(event.weight()); + } + + /// @brief Normalize each AnalysisObject + /// + /// by dividing by the sum of the events seen for each centrality + /// bin. + void normalizePerEvent() { + for (const auto &hist : _histos) + if ( hist.second->numEntries() > 0 && hist.first->numEntries() > 0) + hist.first->scaleW(1./hist.second->val()); + } + + /// Simple scaling of each AnalysisObject. + void scale(float scale) { + for (const auto hist : _histos) + hist.first->scaleW(scale); + } + + /// Execute a function for each AnalysisObject. + void exec(function f) { for ( auto hist : _histos) f(hist); } + + /// @brief Access the underlyng AnalysisObjects + /// + /// The returned vector contains a pair, where the first member is + /// the AnalysisObject and the second is a counter keeping track of + /// the sum of event weights for which the AnalysisObject has been + /// active. + const vector, shared_ptr > > & + analysisObjects() const{ + return _histos; + } + +protected: + + /// The returned vector contains a pair, where the first member is + /// the AnalysisObject and the second is a counter keeping track of + /// the sum of event weights for which the AnalysisObject has been + /// active. + vector, shared_ptr > > _histos; + +}; + +/// @brief The Percentile class for centrality binning. +/// +/// The Percentile class automatically handles the selection of which +/// AnalysisObject(s) should be filled depending on the centrality of +/// an event. It cointains a list of AnalysisObjects, one for each +/// centrality bin requested (note that these bins may be overlapping) +/// and each centrality definition is available in the assigned +/// CentralityProjection. +template +class Percentile : public PercentileTBase { + +public: + + /// @brief the main constructor + /// + /// requiring a pointer, @a ana, to the Analysis to which this + /// object belongs and the name of the CentralityProjection, @a + /// projname, to be used. + Percentile(Analysis * ana, string projName) + : PercentileTBase(ana, projName) {} + + /// @brief Default constructor. + Percentile() {} + + /// @brief Empty destructor. + ~Percentile() {} + + /// Needed to access members of the templated base class. + using PercentileTBase::_histos; + + /// Needed to access members of the templated base class. + using PercentileTBase::_activeBins; + + /// Fill each AnalysisObject selected in the last call to + /// PercentileTBaseinit + template + void fill(Args... args) { + for (const auto bin : _activeBins) { + _histos[bin].first->fill(args...); + } + } + + /// Subtract the contents fro another Pecentile. + Percentile &operator-=(const Percentile &rhs) { + const int nCent = _histos.size(); + for (int iCent = 0; iCent < nCent; ++iCent) { + *_histos[iCent].first -= *rhs._histos[iCent].first; + } + } + + /// Add the contents fro another Pecentile. + Percentile &operator+=(const Percentile &rhs) { + const int nCent = _histos.size(); + for (int iCent = 0; iCent < nCent; ++iCent) { + *_histos[iCent].first += *rhs._histos[iCent].first; + /// @todo should this also add the Counter? + } + } + + /// Make this object look like a pointer. + Percentile *operator->() { return this; } + + /// Pointer to member operator. + Percentile &operator->*(function f) { exec(f); return *this; } + +}; + +/// @brief The PercentileXaxis class for centrality binning. +/// +/// The PercentileXaxis class automatically handles the x-axis of an +/// AnalysisObject when the x-axis is to be the centrality of an +/// event. This could also be done by eg. filling directly a Histo1D +/// with the result of a CentralityProjection. However, since the +/// CentralityProjection may handle several centrality definitions at +/// the same time it is reasonable to instead use +/// PercentileXaxis which will fill one histogram for each +/// centrality definition. +/// +/// Operationally this class works like the Percentile class, but only +/// one centrality bin (0-100) is included. When fill()ed the first +/// argument is always given by the assigned CentralityProjection. +template +class PercentileXaxis : public PercentileTBase { + +public: + + /// @brief the main constructor + /// + /// requiring a pointer, @a ana, to the Analysis to which this + /// object belongs and the name of the CentralityProjection, @a + /// projname, to be used. + PercentileXaxis(Analysis * ana, string projName) + : PercentileTBase(ana, projName) {} + + /// @brief Default constructor. + PercentileXaxis() {} + + /// @brief Empty destructor. + ~PercentileXaxis() {} + + /// Needed to access members of the templated base class. + using PercentileTBase::_histos; + + /// Needed to access members of the templated base class. + using PercentileTBase::_activeBins; + + /// Fill each AnalysisObject selected in the last call to + /// PercentileTBaseinit + template + void fill(Args... args) { + for (const auto bin : _activeBins) { + _histos[bin].first->fill(bin, args...); + } + } + + /// Subtract the contents fro another PecentileXaxis. + PercentileXaxis &operator-=(const PercentileXaxis &rhs) { + const int nCent = _histos.size(); + for (int iCent = 0; iCent < nCent; ++iCent) { + *_histos[iCent].first -= *rhs._histos[iCent].first; + } + } + + /// Add the contents fro another PecentileXaxis. + PercentileXaxis &operator+=(const PercentileXaxis &rhs) { + const int nCent = this->_histos.size(); + for (int iCent = 0; iCent < nCent; ++iCent) { + *_histos[iCent].first += *rhs._histos[iCent].first; + } + } + + /// Make this object look like a pointer. + PercentileXaxis *operator->() { return this; } + + /// Pointer to member operator. + PercentileXaxis &operator->*(function f) { exec(f); return *this; } + +}; + +/// @name Combining Percentiles following the naming of functions for +// the underlying AnalysisObjects: global operators +// @{ + +template +Percentile::RefT> +divide(const Percentile numer, const Percentile denom) { + typedef typename ReferenceTraits::RefT ScatT; + Percentile ret; + vector scatters; + assert( numer.compatible(denom) ); + for ( int i = 0, N = numer.analysisObjects().size(); i < N; ++i ) + scatters.push_back(make_shared(divide(*numer.analysisObjects()[i].first, + *denom.analysisObjects()[i].first))); + ret.add(numer, scatters); + return ret; +} + +template +Percentile::RefT> +divide(const Percentile numer, + const Percentile::RefT> denom) { + typedef typename ReferenceTraits::RefT ScatT; + Percentile ret; + vector scatters; + assert( numer.compatible(denom) ); + for ( int i = 0, N = numer.analysisObjects().size(); i < N; ++i ) + scatters.push_back(make_shared(divide(*numer.analysisObjects()[i].first, + *denom.analysisObjects()[i].first))); + ret.add(numer, scatters); + return ret; +} + +template +Percentile::RefT> +divide(const Percentile::RefT> numer, + const Percentile denom) { + typedef typename ReferenceTraits::RefT ScatT; + Percentile::RefT> ret; + vector scatters; + assert( numer.compatible(denom) ); + for ( int i = 0, N = numer.analysisObjects().size(); i < N; ++i ) + scatters.push_back(make_shared(divide(*numer.analysisObjects()[i].first, + *denom.analysisObjects()[i].first))); + ret.add(numer, scatters); + return ret; +} + +template +Percentile add(const Percentile pctla, const Percentile pctlb) { + Percentile ret; + vector aos; + assert( pctla.compatible(pctlb) ); + for ( int i = 0, N = pctla.analysisObjects().size(); i < N; ++i ) + aos.push_back(make_shared(add(*pctla.analysisObjects()[i].first, + *pctlb.analysisObjects()[i].first))); + ret.add(pctla, aos); + return ret; +} + +template +Percentile::RefT> +add(const Percentile pctla, + const Percentile::RefT> pctlb) { + typedef typename ReferenceTraits::RefT ScatT; + Percentile ret; + vector scatters; + assert( pctla.compatible(pctlb) ); + for ( int i = 0, N = pctla.analysisObjects().size(); i < N; ++i ) + scatters.push_back(make_shared(add(*pctla.analysisObjects()[i].first, + *pctlb.analysisObjects()[i].first))); + ret.add(pctla, scatters); + return ret; +} + +template +Percentile::RefT> +add(const Percentile::RefT> pctla, + const Percentile pctlb) { + typedef typename ReferenceTraits::RefT ScatT; + Percentile ret; + vector scatters; + assert( pctla.compatible(pctlb) ); + for ( int i = 0, N = pctla.analysisObjects().size(); i < N; ++i ) + scatters.push_back(make_shared(add(*pctla.analysisObjects()[i].first, + *pctlb.analysisObjects()[i].first))); + ret.add(pctla, scatters); + return ret; +} + +template +Percentile subtract(const Percentile pctla, const Percentile pctlb) { + Percentile ret; + vector aos; + assert( pctla.compatible(pctlb) ); + for ( int i = 0, N = pctla.analysisObjects().size(); i < N; ++i ) + aos.push_back(make_shared(subtract(*pctla.analysisObjects()[i].first, + *pctlb.analysisObjects()[i].first))); + ret.add(pctla, aos); + return ret; +} + +template +Percentile::RefT> +subtract(const Percentile pctla, + const Percentile::RefT> pctlb) { + typedef typename ReferenceTraits::RefT ScatT; + Percentile ret; + vector scatters; + assert( pctla.compatible(pctlb) ); + for ( int i = 0, N = pctla.analysisObjects().size(); i < N; ++i ) + scatters.push_back(make_shared(subtract(*pctla.analysisObjects()[i].first, + *pctlb.analysisObjects()[i].first))); + ret.add(pctla, scatters); + return ret; +} + +template +Percentile::RefT> +subtract(const Percentile::RefT> pctla, + const Percentile pctlb) { + typedef typename ReferenceTraits::RefT ScatT; + Percentile ret; + vector scatters; + assert( pctla.compatible(pctlb) ); + for ( int i = 0, N = pctla.analysisObjects().size(); i < N; ++i ) + scatters.push_back(make_shared(subtract(*pctla.analysisObjects()[i].first, + *pctlb.analysisObjects()[i].first))); + ret.add(pctla, scatters); + return ret; +} + +template +Percentile::RefT> +multiply(const Percentile pctla, + const Percentile::RefT> pctlb) { + typedef typename ReferenceTraits::RefT ScatT; + Percentile ret; + vector scatters; + assert( pctla.compatible(pctlb) ); + for ( int i = 0, N = pctla.analysisObjects().size(); i < N; ++i ) + scatters.push_back(make_shared(multiply(*pctla.analysisObjects()[i].first, + *pctlb.analysisObjects()[i].first))); + ret.add(pctla, scatters); + return ret; +} + +template +Percentile::RefT> +multiply(const Percentile::RefT> pctla, + const Percentile pctlb) { + typedef typename ReferenceTraits::RefT ScatT; + Percentile ret; + vector scatters; + assert( pctla.compatible(pctlb) ); + for ( int i = 0, N = pctla.analysisObjects().size(); i < N; ++i ) + scatters.push_back(make_shared(multiply(*pctla.analysisObjects()[i].first, + *pctlb.analysisObjects()[i].first))); + ret.add(pctla, scatters); + return ret; +} + + + + +template +PercentileXaxis::RefT> +divide(const PercentileXaxis numer, const PercentileXaxis denom) { + typedef typename ReferenceTraits::RefT ScatT; + PercentileXaxis ret; + vector scatters; + assert( numer.compatible(denom) ); + for ( int i = 0, N = numer.analysisObjects().size(); i < N; ++i ) + scatters.push_back(make_shared(divide(*numer.analysisObjects()[i].first, + *denom.analysisObjects()[i].first))); + ret.add(numer, scatters); + return ret; +} + +template +PercentileXaxis::RefT> +divide(const PercentileXaxis numer, + const PercentileXaxis::RefT> denom) { + typedef typename ReferenceTraits::RefT ScatT; + PercentileXaxis ret; + vector scatters; + assert( numer.compatible(denom) ); + for ( int i = 0, N = numer.analysisObjects().size(); i < N; ++i ) + scatters.push_back(make_shared(divide(*numer.analysisObjects()[i].first, + *denom.analysisObjects()[i].first))); + ret.add(numer, scatters); + return ret; +} + +template +PercentileXaxis::RefT> +divide(const PercentileXaxis::RefT> numer, + const PercentileXaxis denom) { + typedef typename ReferenceTraits::RefT ScatT; + PercentileXaxis::RefT> ret; + vector scatters; + assert( numer.compatible(denom) ); + for ( int i = 0, N = numer.analysisObjects().size(); i < N; ++i ) + scatters.push_back(make_shared(divide(*numer.analysisObjects()[i].first, + *denom.analysisObjects()[i].first))); + ret.add(numer, scatters); + return ret; +} + +template +PercentileXaxis add(const PercentileXaxis pctla, const PercentileXaxis pctlb) { + PercentileXaxis ret; + vector aos; + assert( pctla.compatible(pctlb) ); + for ( int i = 0, N = pctla.analysisObjects().size(); i < N; ++i ) + aos.push_back(make_shared(add(*pctla.analysisObjects()[i].first, + *pctlb.analysisObjects()[i].first))); + ret.add(pctla, aos); + return ret; +} + +template +PercentileXaxis::RefT> +add(const PercentileXaxis pctla, + const PercentileXaxis::RefT> pctlb) { + typedef typename ReferenceTraits::RefT ScatT; + PercentileXaxis ret; + vector scatters; + assert( pctla.compatible(pctlb) ); + for ( int i = 0, N = pctla.analysisObjects().size(); i < N; ++i ) + scatters.push_back(make_shared(add(*pctla.analysisObjects()[i].first, + *pctlb.analysisObjects()[i].first))); + ret.add(pctla, scatters); + return ret; +} + +template +PercentileXaxis::RefT> +add(const PercentileXaxis::RefT> pctla, + const PercentileXaxis pctlb) { + typedef typename ReferenceTraits::RefT ScatT; + PercentileXaxis ret; + vector scatters; + assert( pctla.compatible(pctlb) ); + for ( int i = 0, N = pctla.analysisObjects().size(); i < N; ++i ) + scatters.push_back(make_shared(add(*pctla.analysisObjects()[i].first, + *pctlb.analysisObjects()[i].first))); + ret.add(pctla, scatters); + return ret; +} + +template +PercentileXaxis subtract(const PercentileXaxis pctla, const PercentileXaxis pctlb) { + PercentileXaxis ret; + vector aos; + assert( pctla.compatible(pctlb) ); + for ( int i = 0, N = pctla.analysisObjects().size(); i < N; ++i ) + aos.push_back(make_shared(subtract(*pctla.analysisObjects()[i].first, + *pctlb.analysisObjects()[i].first))); + ret.add(pctla, aos); + return ret; +} + +template +PercentileXaxis::RefT> +subtract(const PercentileXaxis pctla, + const PercentileXaxis::RefT> pctlb) { + typedef typename ReferenceTraits::RefT ScatT; + PercentileXaxis ret; + vector scatters; + assert( pctla.compatible(pctlb) ); + for ( int i = 0, N = pctla.analysisObjects().size(); i < N; ++i ) + scatters.push_back(make_shared(subtract(*pctla.analysisObjects()[i].first, + *pctlb.analysisObjects()[i].first))); + ret.add(pctla, scatters); + return ret; +} + +template +PercentileXaxis::RefT> +subtract(const PercentileXaxis::RefT> pctla, + const PercentileXaxis pctlb) { + typedef typename ReferenceTraits::RefT ScatT; + PercentileXaxis ret; + vector scatters; + assert( pctla.compatible(pctlb) ); + for ( int i = 0, N = pctla.analysisObjects().size(); i < N; ++i ) + scatters.push_back(make_shared(subtract(*pctla.analysisObjects()[i].first, + *pctlb.analysisObjects()[i].first))); + ret.add(pctla, scatters); + return ret; +} + +template +PercentileXaxis::RefT> +multiply(const PercentileXaxis pctla, + const PercentileXaxis::RefT> pctlb) { + typedef typename ReferenceTraits::RefT ScatT; + PercentileXaxis ret; + vector scatters; + assert( pctla.compatible(pctlb) ); + for ( int i = 0, N = pctla.analysisObjects().size(); i < N; ++i ) + scatters.push_back(make_shared(multiply(*pctla.analysisObjects()[i].first, + *pctlb.analysisObjects()[i].first))); + ret.add(pctla, scatters); + return ret; +} + +template +PercentileXaxis::RefT> +multiply(const PercentileXaxis::RefT> pctla, + const PercentileXaxis pctlb) { + typedef typename ReferenceTraits::RefT ScatT; + PercentileXaxis ret; + vector scatters; + assert( pctla.compatible(pctlb) ); + for ( int i = 0, N = pctla.analysisObjects().size(); i < N; ++i ) + scatters.push_back(make_shared(multiply(*pctla.analysisObjects()[i].first, + *pctlb.analysisObjects()[i].first))); + ret.add(pctla, scatters); + return ret; +} + +template +Percentile +operator+(const Percentile