diff --git a/analyses/pluginATLAS/ATLAS_2013_I1230812.cc b/analyses/pluginATLAS/ATLAS_2013_I1230812.cc --- a/analyses/pluginATLAS/ATLAS_2013_I1230812.cc +++ b/analyses/pluginATLAS/ATLAS_2013_I1230812.cc @@ -1,307 +1,303 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/ZFinder.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" namespace Rivet { /// Z + jets in pp at 7 TeV (combined channel / base class) /// @note This base class contains a "mode" variable for combined, e, and mu channel derived classes class ATLAS_2013_I1230812 : public Analysis { public: /// @name Constructors etc. //@{ /// Constructor - ATLAS_2013_I1230812(const string name="ATLAS_2013_I1230812", size_t mode = 0, - const string ref_data="ATLAS_2013_I1230812") : Analysis(name) { - _mode = mode; //default is combined e+mu mode - setRefDataName(ref_data); - } - + DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2013_I1230812); //@} /// Book histograms and initialise projections before the run void init() { // Get options from the new option system + _mode = 0; if ( getOption("ZMODE") == "EL" ) _mode = 1; if ( getOption("ZMODE") == "MU" ) _mode = 2; // Determine the e/mu decay channels used (NB Prompt leptons only). /// @todo Note that Zs are accepted with any rapidity: the cuts are on the e/mu: is this correct? Cut pt20 = Cuts::pT >= 20*GeV; Cut eta_e = _mode? Cuts::abseta < 1.37 || Cuts::absetaIn(1.52, 2.47) : Cuts::abseta < 2.5; Cut eta_m = _mode? Cuts::abseta < 2.4 : Cuts::abseta < 2.5; ZFinder zfinder_el(FinalState(eta_e), pt20, PID::ELECTRON, 66*GeV, 116*GeV); ZFinder zfinder_mu(FinalState(eta_m), pt20, PID::MUON, 66*GeV, 116*GeV); declare(zfinder_el, "zfinder_el"); declare(zfinder_mu, "zfinder_mu"); // Define veto FS in order to prevent Z-decay products entering the jet algorithm VetoedFinalState had_fs; had_fs.addVetoOnThisFinalState(getProjection("zfinder_el")); had_fs.addVetoOnThisFinalState(getProjection("zfinder_mu")); FastJets jets(had_fs, FastJets::ANTIKT, 0.4); jets.useInvisibles(true); declare(jets, "jets"); _h_njet_incl = bookHisto1D( 1, 1, _mode + 1); _h_njet_incl_ratio = bookScatter2D(2, 1, _mode + 1, true); _h_njet_excl = bookHisto1D( 3, 1, _mode + 1); _h_njet_excl_ratio = bookScatter2D(4, 1, _mode + 1, true); _h_njet_excl_pt150 = bookHisto1D( 5, 1, _mode + 1); _h_njet_excl_pt150_ratio = bookScatter2D(6, 1, _mode + 1, true); _h_njet_excl_vbf = bookHisto1D ( 7, 1, _mode + 1); _h_njet_excl_vbf_ratio = bookScatter2D(8, 1, _mode + 1, true); _h_ptlead = bookHisto1D( 9, 1, _mode + 1); _h_ptseclead = bookHisto1D( 10, 1, _mode + 1); _h_ptthirdlead = bookHisto1D( 11, 1, _mode + 1); _h_ptfourthlead = bookHisto1D( 12, 1, _mode + 1); _h_ptlead_excl = bookHisto1D( 13, 1, _mode + 1); _h_pt_ratio = bookHisto1D( 14, 1, _mode + 1); _h_pt_z = bookHisto1D( 15, 1, _mode + 1); _h_pt_z_excl = bookHisto1D( 16, 1, _mode + 1); _h_ylead = bookHisto1D( 17, 1, _mode + 1); _h_yseclead = bookHisto1D( 18, 1, _mode + 1); _h_ythirdlead = bookHisto1D( 19, 1, _mode + 1); _h_yfourthlead = bookHisto1D( 20, 1, _mode + 1); _h_deltay = bookHisto1D( 21, 1, _mode + 1); _h_mass = bookHisto1D( 22, 1, _mode + 1); _h_deltaphi = bookHisto1D( 23, 1, _mode + 1); _h_deltaR = bookHisto1D( 24, 1, _mode + 1); _h_ptthirdlead_vbf = bookHisto1D( 25, 1, _mode + 1); _h_ythirdlead_vbf = bookHisto1D( 26, 1, _mode + 1); _h_ht = bookHisto1D( 27, 1, _mode + 1); _h_st = bookHisto1D( 28, 1, _mode + 1); } /// Perform the per-event analysis void analyze(const Event& event) { FourMomentum z, lp, lm; const ZFinder& zfinder_el = apply(event, "zfinder_el"); const ZFinder& zfinder_mu = apply(event, "zfinder_mu"); if (_mode == 0 && (zfinder_el.constituents().size() + zfinder_mu.constituents().size()) != 2) vetoEvent; if (_mode == 1 && !(zfinder_el.constituents().size() == 2 && zfinder_mu.constituents().empty())) vetoEvent; if (_mode == 2 && !(zfinder_el.constituents().empty() && zfinder_mu.constituents().size() == 2)) vetoEvent; if (zfinder_el.constituents().size() == 2) { z = zfinder_el.boson().momentum(); lp = zfinder_el.constituents()[0].momentum(); lm = zfinder_el.constituents()[1].momentum(); } else if (zfinder_mu.constituents().size() == 2) { z = zfinder_mu.boson().momentum(); lp = zfinder_mu.constituents()[0].momentum(); lm = zfinder_mu.constituents()[1].momentum(); } else vetoEvent; if (deltaR(lp, lm) < 0.2) vetoEvent; Jets jets = apply(event, "jets").jetsByPt(Cuts::pT > 30*GeV && Cuts::absrap < 4.4); ifilter_discard(jets, deltaRLess(lp, 0.5)); ifilter_discard(jets, deltaRLess(lm, 0.5)); const double weight = event.weight(); // Fill jet multiplicities for (size_t ijet = 0; ijet <= jets.size(); ++ijet) { _h_njet_incl->fill(ijet, weight); } _h_njet_excl->fill(jets.size(), weight); // Require at least one jet if (jets.size() >= 1) { // Leading jet histos const double ptlead = jets[0].pT()/GeV; const double yabslead = fabs(jets[0].rapidity()); const double ptz = z.pT()/GeV; _h_ptlead->fill(ptlead, weight); _h_ylead ->fill(yabslead, weight); _h_pt_z ->fill(ptz, weight); // Fill jet multiplicities if (ptlead > 150) _h_njet_excl_pt150->fill(jets.size(), weight); // Loop over selected jets, fill inclusive distributions double st = 0; double ht = lp.pT()/GeV + lm.pT()/GeV; for (size_t ijet = 0; ijet < jets.size(); ++ijet) { ht += jets[ijet].pT()/GeV; st += jets[ijet].pT()/GeV; } _h_ht->fill(ht, weight); _h_st->fill(st, weight); // Require exactly one jet if (jets.size() == 1) { _h_ptlead_excl->fill(ptlead, weight); _h_pt_z_excl ->fill(ptz, weight); } } // Require at least two jets if (jets.size() >= 2) { // Second jet histos const double ptlead = jets[0].pT()/GeV; const double pt2ndlead = jets[1].pT()/GeV; const double ptratio = pt2ndlead/ptlead; const double yabs2ndlead = fabs(jets[1].rapidity()); _h_ptseclead->fill(pt2ndlead, weight); _h_yseclead->fill( yabs2ndlead, weight); _h_pt_ratio->fill( ptratio, weight); // Dijet histos const double deltaphi = fabs(deltaPhi(jets[1], jets[0])); const double deltarap = fabs(jets[0].rapidity() - jets[1].rapidity()) ; const double deltar = fabs(deltaR(jets[0], jets[1], RAPIDITY)); const double mass = (jets[0].momentum() + jets[1].momentum()).mass()/GeV; _h_mass->fill( mass, weight); _h_deltay->fill( deltarap, weight); _h_deltaphi->fill(deltaphi, weight); _h_deltaR->fill( deltar, weight); if (mass > 350 && deltarap > 3) _h_njet_excl_vbf->fill(jets.size(), weight); } // Require at least three jets if (jets.size() >= 3) { // Third jet histos const double pt3rdlead = jets[2].pT()/GeV; const double yabs3rdlead = fabs(jets[2].rapidity()); _h_ptthirdlead->fill(pt3rdlead, weight); _h_ythirdlead->fill( yabs3rdlead, weight); //Histos after VBF preselection const double deltarap = fabs(jets[0].rapidity() - jets[1].rapidity()) ; const double mass = (jets[0].momentum() + jets[1].momentum()).mass(); if (mass > 350 && deltarap > 3) { _h_ptthirdlead_vbf->fill(pt3rdlead, weight); _h_ythirdlead_vbf->fill( yabs3rdlead, weight); } } // Require at least four jets if (jets.size() >= 4) { // Fourth jet histos const double pt4thlead = jets[3].pT()/GeV; const double yabs4thlead = fabs(jets[3].rapidity()); _h_ptfourthlead->fill(pt4thlead, weight); _h_yfourthlead->fill( yabs4thlead, weight); } } /// @name Ratio calculator util functions //@{ /// Calculate the efficiency error, being careful about div-by-zero double err_incl(const HistoBin1D &M, const HistoBin1D &N, bool hasWeights) { double r = safediv(M.sumW(), N.sumW()); if (hasWeights) { // use F. James's approximation for weighted events return sqrt( safediv((1 - 2 * r) * M.sumW2() + r * r * N.sumW2(), N.sumW() * N.sumW()) ); } return sqrt( safediv(r * (1 - r), N.sumW()) ); } /// Calculate the ratio error, being careful about div-by-zero double err_excl(const HistoBin1D &A, const HistoBin1D &B) { double r = safediv(A.sumW(), B.sumW()); double dAsquared = safediv(A.sumW2(), A.sumW() * A.sumW()); // squared relative error of A double dBsquared = safediv(B.sumW2(), B.sumW() * B.sumW()); // squared relative error of B return r * sqrt(dAsquared + dBsquared); } //@} void finalize() { bool hasWeights = _h_njet_incl->effNumEntries() != _h_njet_incl->numEntries(); for (size_t i = 0; i < 6; ++i) { _h_njet_incl_ratio->point(i).setY(safediv(_h_njet_incl->bin(i + 1).sumW(), _h_njet_incl->bin(i).sumW()), err_incl(_h_njet_incl->bin(i + 1), _h_njet_incl->bin(i), hasWeights)); _h_njet_excl_ratio->point(i).setY(safediv(_h_njet_excl->bin(i + 1).sumW(), _h_njet_excl->bin(i).sumW()), err_excl(_h_njet_excl->bin(i + 1), _h_njet_excl->bin(i))); if (i >= 1) { _h_njet_excl_pt150_ratio->point(i - 1).setY(safediv(_h_njet_excl_pt150->bin(i).sumW(), _h_njet_excl_pt150->bin(i - 1).sumW()), err_excl(_h_njet_excl_pt150->bin(i), _h_njet_excl_pt150->bin(i - 1))); if (i >= 2) { _h_njet_excl_vbf_ratio->point(i - 2).setY(safediv(_h_njet_excl_vbf->bin(i).sumW(), _h_njet_excl_vbf->bin(i - 1).sumW()), err_excl(_h_njet_excl_vbf->bin(i), _h_njet_excl_vbf->bin(i - 1))); } } } double sf = _mode? 1.0 : 0.5; const double xs = sf * crossSectionPerEvent()/picobarn; scale({_h_njet_incl, _h_njet_excl, _h_njet_excl_pt150, _h_njet_excl_vbf}, xs); scale({_h_ptlead, _h_ptseclead, _h_ptthirdlead, _h_ptfourthlead, _h_ptlead_excl}, xs); scale({_h_pt_ratio, _h_pt_z, _h_pt_z_excl}, xs); scale({_h_ylead, _h_yseclead, _h_ythirdlead, _h_yfourthlead}, xs); scale({_h_deltay, _h_mass, _h_deltaphi, _h_deltaR}, xs); scale({_h_ptthirdlead_vbf, _h_ythirdlead_vbf}, xs); scale({_h_ht, _h_st}, xs); } //@} protected: size_t _mode; private: Scatter2DPtr _h_njet_incl_ratio; Scatter2DPtr _h_njet_excl_ratio; Scatter2DPtr _h_njet_excl_pt150_ratio; Scatter2DPtr _h_njet_excl_vbf_ratio; Histo1DPtr _h_njet_incl; Histo1DPtr _h_njet_excl; Histo1DPtr _h_njet_excl_pt150; Histo1DPtr _h_njet_excl_vbf; Histo1DPtr _h_ptlead; Histo1DPtr _h_ptseclead; Histo1DPtr _h_ptthirdlead; Histo1DPtr _h_ptfourthlead; Histo1DPtr _h_ptlead_excl; Histo1DPtr _h_pt_ratio; Histo1DPtr _h_pt_z; Histo1DPtr _h_pt_z_excl; Histo1DPtr _h_ylead; Histo1DPtr _h_yseclead; Histo1DPtr _h_ythirdlead; Histo1DPtr _h_yfourthlead; Histo1DPtr _h_deltay; Histo1DPtr _h_mass; Histo1DPtr _h_deltaphi; Histo1DPtr _h_deltaR; Histo1DPtr _h_ptthirdlead_vbf; Histo1DPtr _h_ythirdlead_vbf; Histo1DPtr _h_ht; Histo1DPtr _h_st; }; /*struct ATLAS_2013_I1230812_EL : public ATLAS_2013_I1230812 { ATLAS_2013_I1230812_EL() : ATLAS_2013_I1230812("ATLAS_2013_I1230812_EL", 1) { } }; struct ATLAS_2013_I1230812_MU : public ATLAS_2013_I1230812 { ATLAS_2013_I1230812_MU() : ATLAS_2013_I1230812("ATLAS_2013_I1230812_MU", 2) { } };*/ DECLARE_RIVET_PLUGIN(ATLAS_2013_I1230812); //DECLARE_RIVET_PLUGIN(ATLAS_2013_I1230812_EL); //DECLARE_RIVET_PLUGIN(ATLAS_2013_I1230812_MU); } diff --git a/src/Core/Analysis.cc b/src/Core/Analysis.cc --- a/src/Core/Analysis.cc +++ b/src/Core/Analysis.cc @@ -1,896 +1,895 @@ // -*- C++ -*- #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Analysis.hh" #include "Rivet/AnalysisHandler.hh" #include "Rivet/AnalysisInfo.hh" #include "Rivet/Tools/BeamConstraint.hh" namespace Rivet { Analysis::Analysis(const string& name) : _crossSection(-1.0), _gotCrossSection(false), _analysishandler(NULL) { ProjectionApplier::_allowProjReg = false; _defaultname = name; unique_ptr ai = AnalysisInfo::make(name); assert(ai); _info = move(ai); assert(_info); } double Analysis::sqrtS() const { return handler().sqrtS(); } const ParticlePair& Analysis::beams() const { return handler().beams(); } const PdgIdPair Analysis::beamIds() const { return handler().beamIds(); } const string Analysis::histoDir() const { /// @todo Cache in a member variable string _histoDir; if (_histoDir.empty()) { _histoDir = "/" + name(); if (handler().runName().length() > 0) { _histoDir = "/" + handler().runName() + _histoDir; } replace_all(_histoDir, "//", "/"); //< iterates until none } return _histoDir; } const string Analysis::histoPath(const string& hname) const { const string path = histoDir() + "/" + hname; return path; } const string Analysis::histoPath(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { return histoDir() + "/" + mkAxisCode(datasetId, xAxisId, yAxisId); } const string Analysis::mkAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { stringstream axisCode; axisCode << "d"; if (datasetId < 10) axisCode << 0; axisCode << datasetId; axisCode << "-x"; if (xAxisId < 10) axisCode << 0; axisCode << xAxisId; axisCode << "-y"; if (yAxisId < 10) axisCode << 0; axisCode << yAxisId; return axisCode.str(); } Log& Analysis::getLog() const { string logname = "Rivet.Analysis." + name(); return Log::getLog(logname); } /////////////////////////////////////////// size_t Analysis::numEvents() const { return handler().numEvents(); } double Analysis::sumW() const { return handler().sumW(); } double Analysis::sumW2() const { return handler().sumW(); } /////////////////////////////////////////// bool Analysis::isCompatible(const ParticlePair& beams) const { return isCompatible(beams.first.pid(), beams.second.pid(), beams.first.energy(), beams.second.energy()); } bool Analysis::isCompatible(PdgId beam1, PdgId beam2, double e1, double e2) const { PdgIdPair beams(beam1, beam2); pair energies(e1, e2); return isCompatible(beams, energies); } bool Analysis::isCompatible(const PdgIdPair& beams, const pair& energies) const { // First check the beam IDs bool beamIdsOk = false; foreach (const PdgIdPair& bp, requiredBeams()) { if (compatible(beams, bp)) { beamIdsOk = true; break; } } if (!beamIdsOk) return false; // Next check that the energies are compatible (within 1% or 1 GeV, whichever is larger, for a bit of UI forgiveness) /// @todo Use some sort of standard ordering to improve comparisons, esp. when the two beams are different particles bool beamEnergiesOk = requiredEnergies().size() > 0 ? false : true; typedef pair DoublePair; foreach (const DoublePair& ep, requiredEnergies()) { if ((fuzzyEquals(ep.first, energies.first, 0.01) && fuzzyEquals(ep.second, energies.second, 0.01)) || (fuzzyEquals(ep.first, energies.second, 0.01) && fuzzyEquals(ep.second, energies.first, 0.01)) || (abs(ep.first - energies.first) < 1*GeV && abs(ep.second - energies.second) < 1*GeV) || (abs(ep.first - energies.second) < 1*GeV && abs(ep.second - energies.first) < 1*GeV)) { beamEnergiesOk = true; break; } } return beamEnergiesOk; /// @todo Need to also check internal consistency of the analysis' /// beam requirements with those of the projections it uses. } /////////////////////////////////////////// Analysis& Analysis::setCrossSection(double xs) { _crossSection = xs; _gotCrossSection = true; return *this; } double Analysis::crossSection() const { if (!_gotCrossSection || std::isnan(_crossSection)) { string errMsg = "You did not set the cross section for the analysis " + name(); throw Error(errMsg); } return _crossSection; } double Analysis::crossSectionPerEvent() const { const double sumW = sumOfWeights(); assert(sumW != 0.0); return _crossSection / sumW; } //////////////////////////////////////////////////////////// // Histogramming void Analysis::_cacheRefData() const { if (_refdata.empty()) { MSG_TRACE("Getting refdata cache for paper " << name()); - //_refdata = getRefData(( (info().name().empty()) ? _defaultname : info().name() )); _refdata = getRefData(getRefDataName()); } } CounterPtr Analysis::bookCounter(const string& cname, const string& title) { // const string& xtitle, // const string& ytitle) { const string path = histoPath(cname); CounterPtr ctr = make_shared(path, title); addAnalysisObject(ctr); MSG_TRACE("Made counter " << cname << " for " << name()); // hist->setAnnotation("XLabel", xtitle); // hist->setAnnotation("YLabel", ytitle); return ctr; } CounterPtr Analysis::bookCounter(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const string& title) { // const string& xtitle, // const string& ytitle) { const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId); return bookCounter(axisCode, title); } Histo1DPtr Analysis::bookHisto1D(const string& hname, size_t nbins, double lower, double upper, const string& title, const string& xtitle, const string& ytitle) { Histo1DPtr hist; try { // try to bind to pre-existing // AnalysisObjectPtr ao = getAnalysisObject(path); // hist = dynamic_pointer_cast(ao); hist = getHisto1D(hname); /// @todo Test that cast worked /// @todo Also test that binning is as expected? MSG_TRACE("Bound pre-existing histogram " << hname << " for " << name()); } catch (...) { // binding failed; make it from scratch hist = make_shared(nbins, lower, upper, histoPath(hname), title); addAnalysisObject(hist); MSG_TRACE("Made histogram " << hname << " for " << name()); } hist->setTitle(title); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); return hist; } Histo1DPtr Analysis::bookHisto1D(const string& hname, const vector& binedges, const string& title, const string& xtitle, const string& ytitle) { Histo1DPtr hist; try { // try to bind to pre-existing // AnalysisObjectPtr ao = getAnalysisObject(path); // hist = dynamic_pointer_cast(ao); hist = getHisto1D(hname); /// @todo Test that cast worked /// @todo Also test that binning is as expected? MSG_TRACE("Bound pre-existing histogram " << hname << " for " << name()); } catch (...) { // binding failed; make it from scratch hist = make_shared(binedges, histoPath(hname), title); addAnalysisObject(hist); MSG_TRACE("Made histogram " << hname << " for " << name()); } hist->setTitle(title); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); return hist; } Histo1DPtr Analysis::bookHisto1D(const string& hname, const initializer_list& binedges, const string& title, const string& xtitle, const string& ytitle) { return bookHisto1D(hname, vector{binedges}, title, xtitle, ytitle); } Histo1DPtr Analysis::bookHisto1D(const string& hname, const Scatter2D& refscatter, const string& title, const string& xtitle, const string& ytitle) { Histo1DPtr hist; try { // try to bind to pre-existing // AnalysisObjectPtr ao = getAnalysisObject(path); // hist = dynamic_pointer_cast(ao); hist = getHisto1D(hname); /// @todo Test that cast worked /// @todo Also test that binning is as expected? MSG_TRACE("Bound pre-existing histogram " << hname << " for " << name()); } catch (...) { // binding failed; make it from scratch hist = make_shared(refscatter, histoPath(hname)); if (hist->hasAnnotation("IsRef")) hist->rmAnnotation("IsRef"); addAnalysisObject(hist); MSG_TRACE("Made histogram " << hname << " for " << name()); } hist->setTitle(title); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); return hist; } Histo1DPtr Analysis::bookHisto1D(const string& hname, const string& title, const string& xtitle, const string& ytitle) { const Scatter2D& refdata = refData(hname); return bookHisto1D(hname, refdata, title, xtitle, ytitle); } Histo1DPtr Analysis::bookHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const string& title, const string& xtitle, const string& ytitle) { const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId); return bookHisto1D(axisCode, title, xtitle, ytitle); } /// @todo Add booking methods which take a path, titles and *a reference Scatter from which to book* ///////////////// Histo2DPtr Analysis::bookHisto2D(const string& hname, size_t nxbins, double xlower, double xupper, size_t nybins, double ylower, double yupper, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string path = histoPath(hname); Histo2DPtr hist = make_shared(nxbins, xlower, xupper, nybins, ylower, yupper, path, title); addAnalysisObject(hist); MSG_TRACE("Made 2D histogram " << hname << " for " << name()); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); hist->setAnnotation("ZLabel", ztitle); return hist; } Histo2DPtr Analysis::bookHisto2D(const string& hname, const vector& xbinedges, const vector& ybinedges, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string path = histoPath(hname); Histo2DPtr hist = make_shared(xbinedges, ybinedges, path, title); addAnalysisObject(hist); MSG_TRACE("Made 2D histogram " << hname << " for " << name()); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); hist->setAnnotation("ZLabel", ztitle); return hist; } Histo2DPtr Analysis::bookHisto2D(const string& hname, const initializer_list& xbinedges, const initializer_list& ybinedges, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { return bookHisto2D(hname, vector{xbinedges}, vector{ybinedges}, title, xtitle, ytitle, ztitle); } Histo2DPtr Analysis::bookHisto2D(const string& hname, const Scatter3D& refscatter, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string path = histoPath(hname); Histo2DPtr hist( new Histo2D(refscatter, path) ); addAnalysisObject(hist); MSG_TRACE("Made 2D histogram " << hname << " for " << name()); if (hist->hasAnnotation("IsRef")) hist->rmAnnotation("IsRef"); hist->setTitle(title); hist->setAnnotation("XLabel", xtitle); hist->setAnnotation("YLabel", ytitle); hist->setAnnotation("ZLabel", ztitle); return hist; } Histo2DPtr Analysis::bookHisto2D(const string& hname, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const Scatter3D& refdata = refData(hname); return bookHisto2D(hname, refdata, title, xtitle, ytitle, ztitle); } Histo2DPtr Analysis::bookHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId); return bookHisto2D(axisCode, title, xtitle, ytitle, ztitle); } ///////////////// Profile1DPtr Analysis::bookProfile1D(const string& hname, size_t nbins, double lower, double upper, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Profile1DPtr prof = make_shared(nbins, lower, upper, path, title); addAnalysisObject(prof); MSG_TRACE("Made profile histogram " << hname << " for " << name()); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); return prof; } Profile1DPtr Analysis::bookProfile1D(const string& hname, const vector& binedges, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Profile1DPtr prof = make_shared(binedges, path, title); addAnalysisObject(prof); MSG_TRACE("Made profile histogram " << hname << " for " << name()); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); return prof; } Profile1DPtr Analysis::bookProfile1D(const string& hname, const initializer_list& binedges, const string& title, const string& xtitle, const string& ytitle) { return bookProfile1D(hname, vector{binedges}, title, xtitle, ytitle); } Profile1DPtr Analysis::bookProfile1D(const string& hname, const Scatter2D& refscatter, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Profile1DPtr prof = make_shared(refscatter, path); addAnalysisObject(prof); MSG_TRACE("Made profile histogram " << hname << " for " << name()); if (prof->hasAnnotation("IsRef")) prof->rmAnnotation("IsRef"); prof->setTitle(title); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); return prof; } Profile1DPtr Analysis::bookProfile1D(const string& hname, const string& title, const string& xtitle, const string& ytitle) { const Scatter2D& refdata = refData(hname); return bookProfile1D(hname, refdata, title, xtitle, ytitle); } Profile1DPtr Analysis::bookProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const string& title, const string& xtitle, const string& ytitle) { const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId); return bookProfile1D(axisCode, title, xtitle, ytitle); } /////////////////// Profile2DPtr Analysis::bookProfile2D(const string& hname, size_t nxbins, double xlower, double xupper, size_t nybins, double ylower, double yupper, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string path = histoPath(hname); Profile2DPtr prof = make_shared(nxbins, xlower, xupper, nybins, ylower, yupper, path, title); addAnalysisObject(prof); MSG_TRACE("Made 2D profile histogram " << hname << " for " << name()); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); prof->setAnnotation("ZLabel", ztitle); return prof; } Profile2DPtr Analysis::bookProfile2D(const string& hname, const vector& xbinedges, const vector& ybinedges, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string path = histoPath(hname); Profile2DPtr prof = make_shared(xbinedges, ybinedges, path, title); addAnalysisObject(prof); MSG_TRACE("Made 2D profile histogram " << hname << " for " << name()); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); prof->setAnnotation("ZLabel", ztitle); return prof; } Profile2DPtr Analysis::bookProfile2D(const string& hname, const initializer_list& xbinedges, const initializer_list& ybinedges, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { return bookProfile2D(hname, vector{xbinedges}, vector{ybinedges}, title, xtitle, ytitle, ztitle); } Profile2DPtr Analysis::bookProfile2D(const string& hname, const Scatter3D& refscatter, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string path = histoPath(hname); Profile2DPtr prof( new Profile2D(refscatter, path) ); addAnalysisObject(prof); MSG_TRACE("Made 2D profile histogram " << hname << " for " << name()); if (prof->hasAnnotation("IsRef")) prof->rmAnnotation("IsRef"); prof->setTitle(title); prof->setAnnotation("XLabel", xtitle); prof->setAnnotation("YLabel", ytitle); prof->setAnnotation("ZLabel", ztitle); return prof; } Profile2DPtr Analysis::bookProfile2D(const string& hname, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const Scatter3D& refdata = refData(hname); return bookProfile2D(hname, refdata, title, xtitle, ytitle, ztitle); } Profile2DPtr Analysis::bookProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const string& title, const string& xtitle, const string& ytitle, const string& ztitle) { const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId); return bookProfile2D(axisCode, title, xtitle, ytitle, ztitle); } ///////////////// Scatter2DPtr Analysis::bookScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, bool copy_pts, const string& title, const string& xtitle, const string& ytitle) { const string axisCode = mkAxisCode(datasetId, xAxisId, yAxisId); return bookScatter2D(axisCode, copy_pts, title, xtitle, ytitle); } Scatter2DPtr Analysis::bookScatter2D(const string& hname, bool copy_pts, const string& title, const string& xtitle, const string& ytitle) { Scatter2DPtr s; const string path = histoPath(hname); if (copy_pts) { const Scatter2D& refdata = refData(hname); s = make_shared(refdata, path); for (Point2D& p : s->points()) p.setY(0, 0); } else { s = make_shared(path); } addAnalysisObject(s); MSG_TRACE("Made scatter " << hname << " for " << name()); if (s->hasAnnotation("IsRef")) s->rmAnnotation("IsRef"); s->setTitle(title); s->setAnnotation("XLabel", xtitle); s->setAnnotation("YLabel", ytitle); return s; } Scatter2DPtr Analysis::bookScatter2D(const string& hname, size_t npts, double lower, double upper, const string& title, const string& xtitle, const string& ytitle) { Scatter2DPtr s; const string path = histoPath(hname); try { // try to bind to pre-existing s = getAnalysisObject(hname); /// @todo Also test that binning is as expected? MSG_TRACE("Bound pre-existing scatter " << path << " for " << name()); } catch (...) { // binding failed; make it from scratch s = make_shared(path); const double binwidth = (upper-lower)/npts; for (size_t pt = 0; pt < npts; ++pt) { const double bincentre = lower + (pt + 0.5) * binwidth; s->addPoint(bincentre, 0, binwidth/2.0, 0); } addAnalysisObject(s); MSG_TRACE("Made scatter " << hname << " for " << name()); } s->setTitle(title); s->setAnnotation("XLabel", xtitle); s->setAnnotation("YLabel", ytitle); return s; } Scatter2DPtr Analysis::bookScatter2D(const string& hname, const vector& binedges, const string& title, const string& xtitle, const string& ytitle) { const string path = histoPath(hname); Scatter2DPtr s = make_shared(path); for (size_t pt = 0; pt < binedges.size()-1; ++pt) { const double bincentre = (binedges[pt] + binedges[pt+1]) / 2.0; const double binwidth = binedges[pt+1] - binedges[pt]; s->addPoint(bincentre, 0, binwidth/2.0, 0); } addAnalysisObject(s); MSG_TRACE("Made scatter " << hname << " for " << name()); s->setTitle(title); s->setAnnotation("XLabel", xtitle); s->setAnnotation("YLabel", ytitle); return s; } ///////////////////// void Analysis::divide(CounterPtr c1, CounterPtr c2, Scatter1DPtr s) const { const string path = s->path(); *s = *c1 / *c2; s->setPath(path); } void Analysis::divide(const Counter& c1, const Counter& c2, Scatter1DPtr s) const { const string path = s->path(); *s = c1 / c2; s->setPath(path); } void Analysis::divide(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const { const string path = s->path(); *s = *h1 / *h2; s->setPath(path); } void Analysis::divide(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const { const string path = s->path(); *s = h1 / h2; s->setPath(path); } void Analysis::divide(Profile1DPtr p1, Profile1DPtr p2, Scatter2DPtr s) const { const string path = s->path(); *s = *p1 / *p2; s->setPath(path); } void Analysis::divide(const Profile1D& p1, const Profile1D& p2, Scatter2DPtr s) const { const string path = s->path(); *s = p1 / p2; s->setPath(path); } void Analysis::divide(Histo2DPtr h1, Histo2DPtr h2, Scatter3DPtr s) const { const string path = s->path(); *s = *h1 / *h2; s->setPath(path); } void Analysis::divide(const Histo2D& h1, const Histo2D& h2, Scatter3DPtr s) const { const string path = s->path(); *s = h1 / h2; s->setPath(path); } void Analysis::divide(Profile2DPtr p1, Profile2DPtr p2, Scatter3DPtr s) const { const string path = s->path(); *s = *p1 / *p2; s->setPath(path); } void Analysis::divide(const Profile2D& p1, const Profile2D& p2, Scatter3DPtr s) const { const string path = s->path(); *s = p1 / p2; s->setPath(path); } /// @todo Counter and Histo2D efficiencies and asymms void Analysis::efficiency(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const { const string path = s->path(); *s = YODA::efficiency(*h1, *h2); s->setPath(path); } void Analysis::efficiency(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const { const string path = s->path(); *s = YODA::efficiency(h1, h2); s->setPath(path); } void Analysis::asymm(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const { const string path = s->path(); *s = YODA::asymm(*h1, *h2); s->setPath(path); } void Analysis::asymm(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const { const string path = s->path(); *s = YODA::asymm(h1, h2); s->setPath(path); } void Analysis::scale(CounterPtr cnt, double factor) { if (!cnt) { MSG_WARNING("Failed to scale counter=NULL in analysis " << name() << " (scale=" << factor << ")"); return; } if (std::isnan(factor) || std::isinf(factor)) { MSG_WARNING("Failed to scale counter=" << cnt->path() << " in analysis: " << name() << " (invalid scale factor = " << factor << ")"); factor = 0; } MSG_TRACE("Scaling counter " << cnt->path() << " by factor " << factor); try { cnt->scaleW(factor); } catch (YODA::Exception& we) { MSG_WARNING("Could not scale counter " << cnt->path()); return; } } void Analysis::normalize(Histo1DPtr histo, double norm, bool includeoverflows) { if (!histo) { MSG_WARNING("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << norm << ")"); return; } MSG_TRACE("Normalizing histo " << histo->path() << " to " << norm); try { histo->normalize(norm, includeoverflows); } catch (YODA::Exception& we) { MSG_WARNING("Could not normalize histo " << histo->path()); return; } } void Analysis::scale(Histo1DPtr histo, double factor) { if (!histo) { MSG_WARNING("Failed to scale histo=NULL in analysis " << name() << " (scale=" << factor << ")"); return; } if (std::isnan(factor) || std::isinf(factor)) { MSG_WARNING("Failed to scale histo=" << histo->path() << " in analysis: " << name() << " (invalid scale factor = " << factor << ")"); factor = 0; } MSG_TRACE("Scaling histo " << histo->path() << " by factor " << factor); try { histo->scaleW(factor); } catch (YODA::Exception& we) { MSG_WARNING("Could not scale histo " << histo->path()); return; } } void Analysis::normalize(Histo2DPtr histo, double norm, bool includeoverflows) { if (!histo) { MSG_ERROR("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << norm << ")"); return; } MSG_TRACE("Normalizing histo " << histo->path() << " to " << norm); try { histo->normalize(norm, includeoverflows); } catch (YODA::Exception& we) { MSG_WARNING("Could not normalize histo " << histo->path()); return; } } void Analysis::scale(Histo2DPtr histo, double factor) { if (!histo) { MSG_ERROR("Failed to scale histo=NULL in analysis " << name() << " (scale=" << factor << ")"); return; } if (std::isnan(factor) || std::isinf(factor)) { MSG_ERROR("Failed to scale histo=" << histo->path() << " in analysis: " << name() << " (invalid scale factor = " << factor << ")"); factor = 0; } MSG_TRACE("Scaling histo " << histo->path() << " by factor " << factor); try { histo->scaleW(factor); } catch (YODA::Exception& we) { MSG_WARNING("Could not scale histo " << histo->path()); return; } } void Analysis::integrate(Histo1DPtr h, Scatter2DPtr s) const { // preserve the path info const string path = s->path(); *s = toIntegralHisto(*h); s->setPath(path); } void Analysis::integrate(const Histo1D& h, Scatter2DPtr s) const { // preserve the path info const string path = s->path(); *s = toIntegralHisto(h); s->setPath(path); } /// @todo 2D versions of integrate... defined how, exactly?!? ////////////////////////////////// void Analysis::addAnalysisObject(AnalysisObjectPtr ao) { _analysisobjects.push_back(ao); } void Analysis::removeAnalysisObject(const string& path) { for (vector::iterator it = _analysisobjects.begin(); it != _analysisobjects.end(); ++it) { if ((*it)->path() == path) { _analysisobjects.erase(it); break; } } } void Analysis::removeAnalysisObject(AnalysisObjectPtr ao) { for (vector::iterator it = _analysisobjects.begin(); it != _analysisobjects.end(); ++it) { if (*it == ao) { _analysisobjects.erase(it); break; } } } }