diff --git a/analyses/pluginATLAS/ATLAS_2016_I1502620.cc b/analyses/pluginATLAS/ATLAS_2016_I1502620.cc --- a/analyses/pluginATLAS/ATLAS_2016_I1502620.cc +++ b/analyses/pluginATLAS/ATLAS_2016_I1502620.cc @@ -1,246 +1,239 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/WFinder.hh" #include "Rivet/Projections/ZFinder.hh" namespace Rivet { /// @brief inclusive W/Z cross sections at 7 TeV class ATLAS_2016_I1502620 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2016_I1502620); //@} /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // Get options from the new option system // run everything _mode = 0; _runZ = true; _runW = true; if ( getOption("LMODE") == "EL" || - getOption("LMODE") == "ZEL" || - getOption("LMODE") == "WEL" ) - _mode = 1; + getOption("LMODE") == "ZEL" || + getOption("LMODE") == "WEL" ) + _mode = 1; if ( getOption("LMODE") == "MU" || - getOption("LMODE") == "ZMU" || - getOption("LMODE") == "WMU" ) - _mode = 2; + getOption("LMODE") == "ZMU" || + getOption("LMODE") == "WMU" ) + _mode = 2; if ( getOption("LMODE") == "Z" || - getOption("LMODE") == "ZEL" || - getOption("LMODE") == "ZMU" ) - _runW = false; + getOption("LMODE") == "ZEL" || + getOption("LMODE") == "ZMU" ) + _runW = false; if ( getOption("LMODE") == "W" || - getOption("LMODE") == "WEL" || - getOption("LMODE") == "WMU" ) - _runZ = false; - - + getOption("LMODE") == "WEL" || + getOption("LMODE") == "WMU" ) + _runZ = false; ///Initialise and register projections here const FinalState fs; Cut Wcuts = Cuts::pT >= 25*GeV; // minimum lepton pT Cut Zcuts = Cuts::pT >= 20.0*GeV; WFinder wfinder_edressed(fs, Wcuts, PID::ELECTRON, 40*GeV, 13*TeV, 25*GeV, 0.1, WFinder::CLUSTERNODECAY, WFinder::NOTRACK, WFinder::TRANSMASS); declare(wfinder_edressed, "WFinder_edressed"); ZFinder zfindere(fs, Zcuts, PID::ELECTRON, 46.0*GeV, 150*GeV, 0.1, ZFinder::CLUSTERNODECAY, ZFinder::NOTRACK); declare(zfindere, "ZFindere"); WFinder wfinder_mdressed(fs, Wcuts, PID::MUON, 40*GeV, 13*TeV, 25*GeV, 0.1, WFinder::CLUSTERNODECAY, WFinder::NOTRACK, WFinder::TRANSMASS); declare(wfinder_mdressed, "WFinder_mdressed"); ZFinder zfinderm(fs, Zcuts, PID::MUON, 46.0*GeV, 150*GeV, 0.1, ZFinder::CLUSTERNODECAY, ZFinder::NOTRACK); declare(zfinderm, "ZFinderm"); /// Book histograms here if (_runW) { - _h_Wp_eta = bookHisto1D( 9, 1, 1); - _h_Wm_eta = bookHisto1D( 10, 1, 1); - _h_W_asym = bookScatter2D(35, 1, 1); + _h_Wp_eta = bookHisto1D( 9, 1, 1); + _h_Wm_eta = bookHisto1D( 10, 1, 1); + _h_W_asym = bookScatter2D(35, 1, 1); } if (_runZ) { - _h_Zcenlow_y_dressed = bookHisto1D(11, 1, 1); - _h_Zcenpeak_y_dressed = bookHisto1D(12, 1, 1); - _h_Zcenhigh_y_dressed = bookHisto1D(13, 1, 1); - _h_Zfwdpeak_y_dressed = bookHisto1D(14, 1, 1); - _h_Zfwdhigh_y_dressed = bookHisto1D(15, 1, 1); + _h_Zcenlow_y_dressed = bookHisto1D(11, 1, 1); + _h_Zcenpeak_y_dressed = bookHisto1D(12, 1, 1); + _h_Zcenhigh_y_dressed = bookHisto1D(13, 1, 1); + _h_Zfwdpeak_y_dressed = bookHisto1D(14, 1, 1); + _h_Zfwdhigh_y_dressed = bookHisto1D(15, 1, 1); } } /// Perform the per-event analysis void analyze(const Event& event) { // W stuff const WFinder& wfindere = apply(event, "WFinder_edressed"); const WFinder& wfinderm = apply(event, "WFinder_mdressed"); if (wfindere.bosons().size()+wfinderm.bosons().size() == 1 && _runW) { const double weight = event.weight(); - Particle lep; - if (_mode !=2 && wfindere.bosons().size() == 1 ) { - lep = wfindere.constituentLeptons()[0]; - } - else if (_mode !=1 && wfinderm.bosons().size() == 1 ) { - lep = wfinderm.constituentLeptons()[0]; - } - if (lep.charge3() == 3) { - _h_Wp_eta->fill(lep.abseta()); - } - else if (lep.charge3() == -3) { - _h_Wm_eta->fill(lep.abseta()); - } + Particle lep; + if (_mode !=2 && wfindere.bosons().size() == 1 ) { + lep = wfindere.constituentLeptons()[0]; + } + else if (_mode !=1 && wfinderm.bosons().size() == 1 ) { + lep = wfinderm.constituentLeptons()[0]; + } + if (lep.charge3() == 3) { + _h_Wp_eta->fill(lep.abseta(), weight); + } + else if (lep.charge3() == -3) { + _h_Wm_eta->fill(lep.abseta(), weight); + } } // now the Z stuff. const ZFinder& zfindere = apply(event, "ZFindere"); const ZFinder& zfinderm = apply(event, "ZFinderm"); // must be one and only one candidate. if (zfindere.bosons().size()+zfinderm.bosons().size() == 1 && _runZ) { - const double weight = event.weight(); - - Particle Zboson; - ParticleVector leptons; - - // candidate is e+e- - if (_mode != 2 && zfindere.bosons().size() == 1 ) { - - Zboson = zfindere.boson(); - leptons = zfindere.constituents(); - } + const double weight = event.weight(); + + Particle Zboson; + ParticleVector leptons; + + // candidate is e+e- + if (_mode != 2 && zfindere.bosons().size() == 1 ) { + + Zboson = zfindere.boson(); + leptons = zfindere.constituents(); + } - // candidate is mu+mu- + // candidate is mu+mu- else if (_mode !=1 && zfinderm.bosons().size() == 1 ) { - - Zboson = zfinderm.boson(); - leptons = zfinderm.constituents(); - - } - const double zrap = Zboson.absrap(); - const double zmass = Zboson.mass(); - const double eta1 = leptons[0].abseta(); - const double eta2 = leptons[1].abseta(); - - // separation into central/forward and three mass bins - if (eta1 < 2.5 && eta2 < 2.5) { - if (zmass < 66.0*GeV) _h_Zcenlow_y_dressed->fill(zrap, weight); - else if (zmass < 116.0*GeV) _h_Zcenpeak_y_dressed->fill(zrap, weight); - else _h_Zcenhigh_y_dressed->fill(zrap, weight); - } - else if ((eta1 < 2.5 && 2.5 < eta2 && eta2 < 4.9) || (eta2 < 2.5 && 2.5 < eta1 && eta1 < 4.9)) { - if (zmass < 66.0*GeV) vetoEvent; - if (zmass < 116.0*GeV) _h_Zfwdpeak_y_dressed->fill(zrap, weight); - else _h_Zfwdhigh_y_dressed->fill(zrap, weight); - } + + Zboson = zfinderm.boson(); + leptons = zfinderm.constituents(); + + } + const double zrap = Zboson.absrap(); + const double zmass = Zboson.mass(); + const double eta1 = leptons[0].abseta(); + const double eta2 = leptons[1].abseta(); + + // separation into central/forward and three mass bins + if (eta1 < 2.5 && eta2 < 2.5) { + if (zmass < 66.0*GeV) _h_Zcenlow_y_dressed->fill(zrap, weight); + else if (zmass < 116.0*GeV) _h_Zcenpeak_y_dressed->fill(zrap, weight); + else _h_Zcenhigh_y_dressed->fill(zrap, weight); + } + else if ((eta1 < 2.5 && 2.5 < eta2 && eta2 < 4.9) || (eta2 < 2.5 && 2.5 < eta1 && eta1 < 4.9)) { + if (zmass < 66.0*GeV) vetoEvent; + if (zmass < 116.0*GeV) _h_Zfwdpeak_y_dressed->fill(zrap, weight); + else _h_Zfwdhigh_y_dressed->fill(zrap, weight); + } } - + } /// Normalise histograms etc., after the run void finalize() { // Construct asymmetry: (dsig+/deta - dsig-/deta) / (dsig+/deta + dsig-/deta) //divide(*_h_Wp_eta - *_h_Wm_eta, *_h_Wp_eta + *_h_Wm_eta, _h_W_asym); if (_runW) { - for (size_t i = 0; i < _h_Wp_eta->numBins(); ++i) { - YODA::HistoBin1D& bp = _h_Wp_eta->bin(i); - YODA::HistoBin1D& bm = _h_Wm_eta->bin(i); - const double sum = bp.height() + bm.height(); - //const double xerr = 0.5 * bp.xWidth(); - double val = 0., yerr = 0.; + for (size_t i = 0; i < _h_Wp_eta->numBins(); ++i) { + YODA::HistoBin1D& bp = _h_Wp_eta->bin(i); + YODA::HistoBin1D& bm = _h_Wm_eta->bin(i); + const double sum = bp.height() + bm.height(); + //const double xerr = 0.5 * bp.xWidth(); + double val = 0., yerr = 0.; - if (sum) { - const double pos2 = bp.height() * bp.height(); - const double min2 = bm.height() * bm.height(); - const double errp2 = bp.heightErr() * bp.heightErr(); - const double errm2 = bm.heightErr() * bm.heightErr(); - val = (bp.height() - bm.height()) / sum; - yerr = 2. * sqrt(errm2 * pos2 + errp2 * min2) / (sum * sum); - } - _h_W_asym->addPoint(bp.midpoint(), val, 0.5*bp.xWidth(), yerr); - } + if (sum) { + const double pos2 = bp.height() * bp.height(); + const double min2 = bm.height() * bm.height(); + const double errp2 = bp.heightErr() * bp.heightErr(); + const double errm2 = bm.heightErr() * bm.heightErr(); + val = (bp.height() - bm.height()) / sum; + yerr = 2. * sqrt(errm2 * pos2 + errp2 * min2) / (sum * sum); + } + _h_W_asym->addPoint(bp.midpoint(), val, 0.5*bp.xWidth(), yerr); + } } // Print summary info const double xs_pb(crossSection() / picobarn); const double sumw(sumOfWeights()); MSG_DEBUG( "Cross-section/pb : " << xs_pb ); MSG_DEBUG( "Sum of weights : " << sumw ); MSG_DEBUG( "nEvents : " << numEvents() ); /// Normalise, scale and otherwise manipulate histograms here double lfac = 1.0; // If we have been running on both electrons and muons, need to divide by two to // get the xsec for one flavour if (_mode == 0) lfac = 0.5; const double sf = lfac * 0.5 * xs_pb / sumw; // 0.5 accounts for rapidity bin width if (_runW){ - scale(_h_Wp_eta, sf); - scale(_h_Wm_eta, sf); + scale(_h_Wp_eta, sf); + scale(_h_Wm_eta, sf); } if (_runZ){ - scale(_h_Zcenlow_y_dressed, sf); - scale(_h_Zcenpeak_y_dressed, sf); - scale(_h_Zcenhigh_y_dressed, sf); - scale(_h_Zfwdpeak_y_dressed, sf); - scale(_h_Zfwdhigh_y_dressed, sf); + scale(_h_Zcenlow_y_dressed, sf); + scale(_h_Zcenpeak_y_dressed, sf); + scale(_h_Zcenhigh_y_dressed, sf); + scale(_h_Zfwdpeak_y_dressed, sf); + scale(_h_Zfwdhigh_y_dressed, sf); } } //@} protected: size_t _mode; bool _runZ, _runW; private: /// @name Histograms //@{ Histo1DPtr _h_Wp_eta, _h_Wm_eta; Scatter2DPtr _h_W_asym; Histo1DPtr _h_Zcenlow_y_dressed; Histo1DPtr _h_Zcenpeak_y_dressed; Histo1DPtr _h_Zcenhigh_y_dressed; Histo1DPtr _h_Zfwdpeak_y_dressed; Histo1DPtr _h_Zfwdhigh_y_dressed; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2016_I1502620); } -// END END END - - - - diff --git a/analyses/pluginBES/BESIII_2012_I1113599.cc b/analyses/pluginBES/BESIII_2012_I1113599.cc --- a/analyses/pluginBES/BESIII_2012_I1113599.cc +++ b/analyses/pluginBES/BESIII_2012_I1113599.cc @@ -1,126 +1,127 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/Beam.hh" #include "Rivet/Projections/FinalState.hh" +#include "Rivet/Projections/UnstableParticles.hh" namespace Rivet { /// @brief J/psi to p pbar n nbar class BESIII_2012_I1113599 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(BESIII_2012_I1113599); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // Initialise and register projections declare(Beam(), "Beams"); declare(UnstableParticles(), "UFS"); declare(FinalState(), "FS"); // Book histograms _h_proton = bookHisto1D("ctheta_p",20,-1.,1.); _h_neutron = bookHisto1D("ctheta_n",20,-1.,1.); } /// Perform the per-event analysis void analyze(const Event& event) { // get the axis, direction of incoming electron const ParticlePair& beams = apply(event, "Beams").beams(); Vector3 axis; if(beams.first.pdgId()>0) axis = beams.first .momentum().p3().unit(); else axis = beams.second.momentum().p3().unit(); // types of final state particles const FinalState& fs = apply(event, "FS"); map nCount; int ntotal(0); Particle outgoing; foreach (const Particle& p, fs.particles()) { nCount[p.pdgId()] += 1; if(p.pdgId()==2212 || p.pdgId()==2112) outgoing = p; ++ntotal; } if(ntotal==2) { if(nCount[2212]==1 && nCount[-2212]==1) { _h_proton->fill(outgoing.momentum().p3().unit().dot(axis),event.weight()); } else if(nCount[2112]==1 && nCount[-2112]==1) { _h_neutron->fill(outgoing.momentum().p3().unit().dot(axis),event.weight()); } } } pair > calcAlpha(Histo1DPtr hist) { if(hist->numEntries()==0.) return make_pair(0.,make_pair(0.,0.)); double sum1(0.),sum2(0.),sum3(0.),sum4(0.),sum5(0.); for (auto bin : hist->bins() ) { double Oi = bin.area(); if(Oi==0.) continue; double a = 1.5*(bin.xMax() - bin.xMin()); double b = 0.5*(pow(bin.xMax(),3) - pow(bin.xMin(),3)); double Ei = bin.areaErr(); sum1 += a*Oi/sqr(Ei); sum2 += b*Oi/sqr(Ei); sum3 += sqr(a)/sqr(Ei); sum4 += sqr(b)/sqr(Ei); sum5 += a*b/sqr(Ei); } // calculate alpha double alpha = (-3*sum1 + 9*sum2 + sum3 - 3*sum5)/(sum1 - 3*sum2 + 3*sum4 - sum5); // and error double cc = -pow((sum3 + 9*sum4 - 6*sum5),3); double bb = -2*sqr(sum3 + 9*sum4 - 6*sum5)*(sum1 - 3*sum2 + 3*sum4 - sum5); double aa = sqr(sum1 - 3*sum2 + 3*sum4 - sum5)*(-sum3 - 9*sum4 + sqr(sum1 - 3*sum2 + 3*sum4 - sum5) + 6*sum5); double dis = sqr(bb)-4.*aa*cc; if(dis>0.) { dis = sqrt(dis); return make_pair(alpha,make_pair(0.5*(-bb+dis)/aa,-0.5*(-bb-dis)/aa)); } else { return make_pair(alpha,make_pair(0.,0.)); } } /// Normalise histograms etc., after the run void finalize() { // proton normalize(_h_proton ); pair > alpha = calcAlpha(_h_proton); Scatter2DPtr _h_alpha_proton = bookScatter2D(1,1,1); _h_alpha_proton->addPoint(0.5, alpha.first, make_pair(0.5,0.5), make_pair(alpha.second.first,alpha.second.second) ); // neutron normalize(_h_neutron); alpha = calcAlpha(_h_neutron); Scatter2DPtr _h_alpha_neutron = bookScatter2D(1,1,2); _h_alpha_neutron->addPoint(0.5, alpha.first, make_pair(0.5,0.5), make_pair(alpha.second.first,alpha.second.second) ); } //@} /// @name Histograms //@{ Histo1DPtr _h_proton,_h_neutron; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(BESIII_2012_I1113599); } diff --git a/analyses/pluginLEP/OPAL_2003_I611415.cc b/analyses/pluginLEP/OPAL_2003_I611415.cc --- a/analyses/pluginLEP/OPAL_2003_I611415.cc +++ b/analyses/pluginLEP/OPAL_2003_I611415.cc @@ -1,173 +1,173 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/GammaGammaFinalState.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { /// @brief Add a short analysis description here class OPAL_2003_I611415 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(OPAL_2003_I611415); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // get the hadronic final state const GammaGammaKinematics& diskin = declare(GammaGammaKinematics(), "Kinematics"); const FinalState & fs = declare(GammaGammaFinalState(diskin, GammaGammaFinalState::LAB), "FS"); declare(FastJets(fs, FastJets::KT,1.),"Jets"); _h_theta[0] = bookHisto1D( 1,1,1); _h_theta[1] = bookHisto1D( 2,1,1); _h_ET[0] = bookHisto1D( 3,1,1); _h_ET[1] = bookHisto1D( 4,1,1); _h_ET[2] = bookHisto1D( 5,1,1); _h_xg[0][0] = bookHisto1D( 6,1,1); _h_xg[0][1] = bookHisto1D( 7,1,1); _h_xg[1][0] = bookHisto1D( 9,1,1); _h_xg[1][1] = bookHisto1D(10,1,1); _h_xg[2][0] = bookHisto1D(11,1,1); _h_xg[2][1] = bookHisto1D(12,1,1); _h_xg_high = bookHisto1D( 8,1,1); _h_xlog[0] = bookHisto1D(13,1,1); _h_xlog[1] = bookHisto1D(14,1,1); _h_xlog[2] = bookHisto1D(15,1,1); _h_eta_diff[0] = bookHisto1D(16,1,1); _h_eta_diff[1] = bookHisto1D(17,1,1); _h_eta_min[0] = bookHisto1D(18,1,1); _h_eta_min[1] = bookHisto1D(19,1,1); - _h_eta_min[2] = bookHisto1D(20,1,1); - _h_eta_min[3] = bookHisto1D(21,1,1); + _h_eta_max[0] = bookHisto1D(20,1,1); + _h_eta_max[1] = bookHisto1D(21,1,1); } /// Perform the per-event analysis void analyze(const Event& event) { double weight = event.weight(); // need at least two jets with |eta|<2 and pT>3 Jets jets = apply(event, "Jets").jetsByPt(Cuts::Et > 3.*GeV and Cuts::abseta < 2.); if(jets.size()<2) vetoEvent; if(jets[0].Et().25) vetoEvent; // calculate x_gamma FourMomentum psum; for(const Particle & part : apply(event,"FS").particles()) { psum += part.momentum(); } FourMomentum pj = jets[0].momentum()+jets[1].momentum(); double xp = (pj.E()+pj.pz())/(psum.E()+psum.pz()); double xm = (pj.E()-pj.pz())/(psum.E()-psum.pz()); double cost = tanh(0.5*(jets[0].eta()-jets[1].eta())); // cost distributions if(pj.mass()>15.*GeV && etaBar<=1.) { if(xp>0.75 && xm>0.75) _h_theta[0]->fill(abs(cost),weight); else if(xp<0.75 && xm<0.75) _h_theta[1]->fill(abs(cost),weight); } // ET distributions _h_ET[0]->fill(Etbar,weight); if((xp<0.75 && xm>0.75)|| (xm<0.75&&xp>0.75)) _h_ET[1]->fill(Etbar,weight); else if(xp<0.75 && xm <0.75) _h_ET[2]->fill(Etbar,weight); if(Etbar>=5.&&Etbar<7.) { _h_xg[0][0]->fill(xp,weight); _h_xg[0][0]->fill(xm,weight); _h_xlog[0]->fill(log(xp),weight); _h_xlog[0]->fill(log(xm),weight); if((xp<0.75 && xm>0.75)|| (xm<0.75&&xp>0.75)) { _h_xg[1][0]->fill(xp,weight); _h_xg[1][0]->fill(xm,weight); _h_xlog[1]->fill(log(xp),weight); _h_xlog[1]->fill(log(xm),weight); } else if(xp<0.75 && xm <0.75) { _h_xg[2][0]->fill(xp,weight); _h_xg[2][0]->fill(xm,weight); _h_xlog[2]->fill(log(xp),weight); _h_xlog[2]->fill(log(xm),weight); } } else if(Etbar>=7.&& Etbar<11.) { _h_xg[0][1]->fill(xp,weight); _h_xg[0][1]->fill(xm,weight); if((xp<0.75 && xm>0.75)|| (xm<0.75&&xp>0.75)) { _h_xg[1][1]->fill(xp,weight); _h_xg[1][1]->fill(xm,weight); } else if(xp<0.75 && xm <0.75) { _h_xg[2][1]->fill(xp,weight); _h_xg[2][1]->fill(xm,weight); } } else if(Etbar>=11.&& Etbar<25.) { _h_xg_high->fill(xp,weight); _h_xg_high->fill(xm,weight); } // vs eta double etaMin = min(abs(jets[0].eta()),abs(jets[1].eta())); double etaMax = max(abs(jets[0].eta()),abs(jets[1].eta())); if((xp<0.75 && xm>0.75)|| (xm<0.75&&xp>0.75)) { _h_eta_diff[0]->fill(abs(jets[0].eta()-jets[1].eta()),weight); _h_eta_min[0]->fill(etaMin,weight); _h_eta_max[0]->fill(etaMax,weight); } else if(xp<0.75 && xm <0.75) { _h_eta_diff[1]->fill(abs(jets[0].eta()-jets[1].eta()),weight); _h_eta_min[1]->fill(etaMin,weight); _h_eta_max[1]->fill(etaMax,weight); } } /// Normalise histograms etc., after the run void finalize() { double fact = crossSection()/picobarn/sumOfWeights(); for(unsigned int ix=0;ix<2;++ix) { scale(_h_theta[ix], fact); scale(_h_eta_diff[ix], fact); scale(_h_eta_min[ix], fact); scale(_h_eta_max[ix], fact); for(unsigned int iy=0;iy<3;++iy) { scale(_h_xg[iy][ix],fact); } } for(unsigned int ix=0;ix<3;++ix) { scale(_h_ET[ix],fact); scale(_h_xlog[ix],fact); } scale(_h_xg_high,fact); } //@} /// @name Histograms //@{ Histo1DPtr _h_theta[2],_h_ET[3],_h_xg[3][2],_h_xg_high; Histo1DPtr _h_xlog[3],_h_eta_diff[2],_h_eta_min[2],_h_eta_max[2]; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(OPAL_2003_I611415); } diff --git a/bin/rivet-cmphistos b/bin/rivet-cmphistos --- a/bin/rivet-cmphistos +++ b/bin/rivet-cmphistos @@ -1,501 +1,501 @@ #! /usr/bin/env python """\ Generate histogram comparison plots USAGE: %(prog)s [options] yodafile1[:'PlotOption1=Value':'PlotOption2=Value':...] [path/to/yodafile2 ...] [PLOT:Key1=Val1:...] where the plot options are described in the make-plots manual in the HISTOGRAM section. ENVIRONMENT: * RIVET_ANALYSIS_PATH: list of paths to be searched for plugin analysis libraries at runtime * RIVET_DATA_PATH: list of paths to be searched for data files """ from __future__ import print_function import rivet, yoda, sys, os rivet.util.check_python_version() rivet.util.set_process_name(os.path.basename(__file__)) class Plot(dict): "A tiny Plot object to help writing out the head in the .dat file" def __repr__(self): return "# BEGIN PLOT\n" + "\n".join("%s=%s" % (k,v) for k,v in self.items()) + "\n# END PLOT\n\n" def sanitiseString(s): #s = s.replace('_','\\_') #s = s.replace('^','\\^{}') #s = s.replace('$','\\$') s = s.replace('#','\\#') s = s.replace('%','\\%') return s def getHistos(filelist): """Loop over all input files. Only use the first occurrence of any REF-histogram and the first occurrence in each MC file for every MC-histogram.""" refhistos, mchistos = {}, {} for infile in filelist: mchistos.setdefault(infile, {}) analysisobjects = yoda.read(infile, patterns=args.PATHPATTERNS, unpatterns=args.PATHUNPATTERNS) #print(analysisobjects) for path, ao in analysisobjects.items(): ## We can't plot non-histograms yet # TODO: support counter plotting with a faked x (or y) position and forced plot width/height - if ao.type not in ("Counter", + if ao.type() not in ("Counter", "Histo1D", "Histo2D", "Profile1D", "Profile2D", "Scatter1D", "Scatter2D", "Scatter3D"): continue ## Make a path object and ensure the path is in standard form try: aop = rivet.AOPath(path) except Exception as e: #print(e) print("Found analysis object with non-standard path structure:", path, "... skipping") continue ## We don't plot data objects with path components hidden by an underscore prefix if aop.istmp() or aop.israw(): continue ## Add it to the ref or mc paths, if this path isn't already known basepath = aop.basepath(keepref=False) if aop.isref() and basepath not in refhistos: ao.path = aop.varpath(keepref=False, defaultvarid=0) refhistos[basepath] = ao else: #if basepath not in mchistos[infile]: mchistos[infile].setdefault(basepath, {})[aop.varid(0)] = ao return refhistos, mchistos def getRivetRefData(anas=None): "Find all Rivet reference data files" refhistos = {} rivet_data_dirs = rivet.getAnalysisRefPaths() dirlist = [] for d in rivet_data_dirs: if anas is None: import glob dirlist.append(glob.glob(os.path.join(d, '*.yoda'))) else: dirlist.append([os.path.join(d, a+'.yoda') for a in anas]) for filelist in dirlist: # TODO: delegate to getHistos? for infile in filelist: analysisobjects = yoda.read(infile, patterns=args.PATHPATTERNS, unpatterns=args.PATHUNPATTERNS) for path, ao in analysisobjects.items(): aop = rivet.AOPath(ao.path) if aop.isref(): ao.path = aop.basepath(keepref=False) refhistos[ao.path] = ao return refhistos def parseArgs(args): """Look at the argument list and split it at colons, in order to separate the file names from the plotting options. Store the file names and file specific plotting options.""" filelist = [] plotoptions = {} for a in args: asplit = a.split(':') path = asplit[0] filelist.append(path) plotoptions[path] = [] has_title = False for i in range(1, len(asplit)): ## Add 'Title' if there is no = sign before math mode if '=' not in asplit[i] or ('$' in asplit[i] and asplit[i].index('$') < asplit[i].index('=')): asplit[i] = 'Title=%s' % asplit[i] if asplit[i].startswith('Title='): has_title = True plotoptions[path].append(asplit[i]) if path != "PLOT" and not has_title: plotoptions[path].append('Title=%s' % sanitiseString(os.path.basename( os.path.splitext(path)[0] )) ) return filelist, plotoptions def setStyle(ao, istyle, variation=False): """Set default plot styles (color and line width) colors borrowed from Google Ngrams""" # LINECOLORS = ['{[HTML]{EE3311}}', # red (Google uses 'DC3912') # '{[HTML]{3366FF}}', # blue # '{[HTML]{109618}}', # green # '{[HTML]{FF9900}}', # orange # '{[HTML]{990099}}'] # lilac LINECOLORS = ['red', 'blue', 'green', 'orange', 'lilac'] LINESTYLES = ['solid', 'dashed', 'dashdotted', 'dotted'] if args.STYLE == 'talk': ao.setAnnotation('LineWidth', '1pt') if args.STYLE == 'bw': LINECOLORS = ['black!90', 'black!50', 'black!30'] jc = istyle % len(LINECOLORS) c = LINECOLORS[jc] js = (istyle // len(LINECOLORS)) % len(LINESTYLES) s = LINESTYLES[js] ## If plotting a variation (i.e. band), fade the colour if variation: c += "!30" ao.setAnnotation('LineStyle', '%s' % s) ao.setAnnotation('LineColor', '%s' % c) def setOptions(ao, options): "Set arbitrary annotations" for opt in options: key, val = opt.split('=', 1) ao.setAnnotation(key, val) # TODO: move to rivet.utils def mkoutdir(outdir): "Function to make output directories" if not os.path.exists(outdir): try: os.makedirs(outdir) except: msg = "Can't make output directory '%s'" % outdir raise Exception(msg) if not os.access(outdir, os.W_OK): msg = "Can't write to output directory '%s'" % outdir raise Exception(msg) def mkOutput(hpath, aos, plot=None, special=None): """ Make the .dat file string. We can't use "yoda.writeFLAT(anaobjects, 'foobar.dat')" because the PLOT and SPECIAL blocks don't have a corresponding analysis object. """ output = '' if plot is not None: output += str(plot) if special is not None: output += "\n" output += "# BEGIN SPECIAL %s\n" % hpath output += special output += "# END SPECIAL\n\n" from io import StringIO sio = StringIO() yoda.writeFLAT(aos, sio) output += sio.getvalue() return output def writeOutput(output, h): "Choose output file name and dir" if args.HIER_OUTPUT: hparts = h.strip("/").split("/", 1) ana = "_".join(hparts[:-1]) if len(hparts) > 1 else "ANALYSIS" outdir = os.path.join(args.OUTDIR, ana) outfile = '%s.dat' % hparts[-1].replace("/", "_") else: hparts = h.strip("/").split("/") outdir = args.OUTDIR outfile = '%s.dat' % "_".join(hparts) mkoutdir(outdir) outfilepath = os.path.join(outdir, outfile) f = open(outfilepath, 'w') f.write(output) f.close() #-------------------------------------------------------------------------------------------- if __name__ == '__main__': ## Command line parsing import argparse parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("ARGS", nargs="*") parser.add_argument('-o', '--outdir', dest='OUTDIR', default='.', help='write data files into this directory') parser.add_argument("--hier-out", action="store_true", dest="HIER_OUTPUT", default=False, help="write output dat files into a directory hierarchy which matches the analysis paths") parser.add_argument('--plotinfodir', dest='PLOTINFODIRS', action='append', default=['.'], help='directory which may contain plot header information (in addition ' 'to standard Rivet search paths)') parser.add_argument("--no-rivet-refs", dest="RIVETREFS", action="store_false", default=True, help="don't use Rivet reference data files") # parser.add_argument("--refid", dest="REF_ID", # default="REF", help="ID of reference data set (file path for non-REF data)") parser.add_argument("--reftitle", dest="REFTITLE", default='Data', help="Reference data legend entry") parser.add_argument("--pwd", dest="PATH_PWD", action="store_true", default=False, help="append the current directory (pwd) to the analysis/data search paths (cf. $RIVET_ANALYSIS/DATA_PATH)") parser.add_argument("-v", "--verbose", dest="VERBOSE", action="store_true", default=False, help="produce debug output to the terminal") stygroup = parser.add_argument_group("Plot style") stygroup.add_argument("--linear", action="store_true", dest="LINEAR", default=False, help="plot with linear scale") stygroup.add_argument("--errs", "--mcerrs", "--mc-errs", action="store_true", dest="MC_ERRS", default=False, help="show vertical error bars on the MC lines") stygroup.add_argument("--no-ratio", action="store_false", dest="RATIO", default=True, help="disable the ratio plot") stygroup.add_argument("--rel-ratio", action="store_true", dest="RATIO_DEVIATION", default=False, help="show the ratio plots scaled to the ref error") stygroup.add_argument("--no-plottitle", action="store_true", dest="NOPLOTTITLE", default=False, help="don't show the plot title on the plot " "(useful when the plot description should only be given in a caption)") stygroup.add_argument("--style", dest="STYLE", default="default", help="change plotting style: default|bw|talk") stygroup.add_argument("-c", "--config", dest="CONFIGFILES", action="append", default=["~/.make-plots"], help="additional plot config file(s). Settings will be included in the output configuration.") stygroup.add_argument("--remove-options", help="remove options label from legend", dest="REMOVE_OPTIONS", action="store_true", default=False) selgroup = parser.add_argument_group("Selective plotting") # selgroup.add_argument("--show-single", dest="SHOW_SINGLE", choices=("no", "ref", "mc", "all"), # default="mc", help="control if a plot file is made if there is only one dataset to be plotted " # "[default=%(default)s]. If the value is 'no', single plots are always skipped, for 'ref' and 'mc', " # "the plot will be written only if the single plot is a reference plot or an MC " # "plot respectively, and 'all' will always create single plot files.\n The 'ref' and 'all' values " # "should be used with great care, as they will also write out plot files for all reference " # "histograms without MC traces: combined with the -R/--rivet-refs flag, this is a great way to " # "write out several thousand irrelevant reference data histograms!") # selgroup.add_argument("--show-mc-only", "--all", action="store_true", dest="SHOW_IF_MC_ONLY", # default=False, help="make a plot file even if there is only one dataset to be plotted and " # "it is an MC one. Deprecated and will be removed: use --show-single instead, which overrides this.") # # selgroup.add_argument("-l", "--histogram-list", dest="HISTOGRAMLIST", # # default=None, help="specify a file containing a list of histograms to plot, in the format " # # "/ANALYSIS_ID/histoname, one per line, e.g. '/DELPHI_1996_S3430090/d01-x01-y01'.") selgroup.add_argument("-m", "--match", action="append", help="only write out histograms whose $path/$name string matches these regexes. The argument " "may also be a text file.", dest="PATHPATTERNS") selgroup.add_argument("-M", "--unmatch", action="append", help="exclude histograms whose $path/$name string matches these regexes", dest="PATHUNPATTERNS") args = parser.parse_args() ## Add pwd to search paths if args.PATH_PWD: rivet.addAnalysisLibPath(os.path.abspath(".")) rivet.addAnalysisDataPath(os.path.abspath(".")) ## Split the input file names and the associated plotting options ## given on the command line into two separate lists filelist, plotoptions = parseArgs(args.ARGS) ## Remove the PLOT dummy file from the file list if "PLOT" in filelist: filelist.remove("PLOT") ## Check that the files exist for f in filelist: if not os.access(f, os.R_OK): print("Error: cannot read from %s" % f) sys.exit(1) ## Read the .plot files plotdirs = args.PLOTINFODIRS + [os.path.abspath(os.path.dirname(f)) for f in filelist] plotparser = rivet.mkStdPlotParser(plotdirs, args.CONFIGFILES) ## Create a list of all histograms to be plotted, and identify if they are 2D histos (which need special plotting) try: refhistos, mchistos = getHistos(filelist) except IOError as e: print("File reading error: ", e.strerror) exit(1) hpaths, h2ds = [], [] for aos in mchistos.values(): for p in aos.keys(): ps = rivet.stripOptions(p) if ps and ps not in hpaths: hpaths.append(ps) firstaop = aos[p][sorted(aos[p].keys())[0]] # TODO: Would be nicer to test via isHisto and dim or similar, or yoda.Scatter/Histo/Profile base classes if type(firstaop) in (yoda.Histo2D, yoda.Profile2D, yoda.Scatter3D) and ps not in h2ds: h2ds.append(ps) # Unique list of analyses anas = list(set([x.split("/")[1] for x in hpaths])) ## Take reference data from the Rivet search paths, if there is not already if args.RIVETREFS: try: refhistos2 = getRivetRefData(anas) except IOError as e: print("File reading error: ", e.strerror) exit(1) refhistos2.update(refhistos) refhistos = refhistos2 ## Purge unmatched ref data entries to save memory keylist = list(refhistos.keys()) # can't modify for-loop target for refhpath in keylist: if refhpath not in hpaths: del refhistos[refhpath] ## Now loop over all MC histograms and plot them # TODO: factorize much of this into a rivet.utils mkplotfile(mchists, refhist, kwargs, is2d=False) function for hpath in hpaths: #print('Currently looking at', hpath) ## The analysis objects to be plotted anaobjects = [] ## List of histos to be drawn, to sync the legend and plotted lines mainlines = [] varlines = [] ## Is this a 2D histo? is2d = (hpath in h2ds) ## Will we be drawing a ratio plot? showratio = args.RATIO and not is2d ## A Plot object to represent the PLOT section in the .dat file plot = Plot() if not is2d: plot['Legend'] = '1' plot['LogY'] = '1' headers = plotparser.getHeaders(hpath) if headers: plot.update(headers) # for key, val in headers.items(): # plot[key] = val if "PLOT" in plotoptions: for key_val in plotoptions["PLOT"]: key, val = [s.strip() for s in key_val.split("=", 1)] if 'ReplaceOption' in key_val: opt_group = key_val.split("=", 2) key, val = "=".join(opt_group[:2]), opt_group[2] plot[key] = val if args.REMOVE_OPTIONS: plot['RemoveOptions'] = '1' if args.LINEAR: plot['LogY'] = '0' if args.NOPLOTTITLE: plot['Title'] = '' if showratio and args.RATIO_DEVIATION: plot['RatioPlotMode'] = 'deviation' if args.STYLE == 'talk': plot['PlotSize'] = '8,6' elif args.STYLE == 'bw' and showratio: plot['RatioPlotErrorBandColor'] = 'black!10' ## Get a special object, if there is one for this path special = plotparser.getSpecial(hpath) ## Handle reference data histogram, if there is one ratioreference, hasdataref = None, False if hpath in refhistos: hasdataref = True refdata = refhistos[hpath] refdata.setAnnotation('Title', args.REFTITLE) if not is2d: refdata.setAnnotation('ErrorBars', '1') refdata.setAnnotation('PolyMarker', '*') refdata.setAnnotation('ConnectBins', '0') if showratio: ratioreference = hpath ## For 1D anaobjects.append(refdata) mainlines.append(hpath) ## For 2D if is2d: s = mkOutput(hpath, [refdata], plot, special) writeOutput(s, hpath) ## Loop over the MC files to plot all instances of the histogram styleidx = 0 for infile in filelist: if infile in mchistos: for xpath in sorted(mchistos[infile]): if rivet.stripOptions(xpath) != hpath: continue hmcs = mchistos[infile][xpath] ## For now, just plot all the different variation histograms (reversed, so [0] is on top) # TODO: calculate and plot an appropriate error band, somehow... for i in sorted(hmcs.keys(), reverse=True): iscanonical = (str(i) == "0") hmc = hmcs[i] ## Default linecolor, linestyle if not is2d: setStyle(hmc, styleidx, not iscanonical) if args.MC_ERRS: hmc.setAnnotation('ErrorBars', '1') ## Plot defaults from .plot files histopts = plotparser.getHistogramOptions(hpath) if histopts: for key, val in histopts.items(): hmc.setAnnotation(key, val) ## Command line plot options setOptions(hmc, plotoptions[infile]) ## Set path attribute fullpath = "/"+infile+xpath if not iscanonical: fullpath += "["+str(i)+"]" hmc.setAnnotation('Path', fullpath) ## Add object / path to appropriate lists #if hmc.hasAnnotation("Title"): # hmc.setAnnotation("Title", hmc.annotation("Title") + # rivet.extractOptionString(xpath)) if not is2d : anaobjects.append(hmc) if iscanonical: mainlines.append(fullpath) else: varlines.append(fullpath) if showratio and ratioreference is None and iscanonical: ratioreference = fullpath ## For 2D, plot each histo now (since overlay makes no sense) if is2d: suffix=(infile.split("/")[-1].replace(".yoda","")) s = mkOutput(hpath, [hmc], plot, special) writeOutput(s, "%s_%s" % (hpath,suffix)) styleidx += 1 ## Finally render the combined plots; only show the first one if it's 2D # TODO: Only show the first *MC* one if 2D? if is2d: anaobjects = anaobjects[:1] ## Add final attrs to Plot if 'DrawOnly' not in plot: plot['DrawOnly'] = ' '.join(varlines + mainlines).strip() if 'LegendOnly' not in plot: plot['LegendOnly'] = ' '.join(mainlines).strip() if showratio and 'RatioPlot' not in plot and len(varlines + mainlines) > 1: plot['RatioPlot'] = '1' if showratio and 'RatioPlot' in plot and plot['RatioPlot'] == '1' and len(varlines + mainlines) > 0: if 'RatioPlotReference' not in plot: plot['RatioPlotReference'] = ratioreference if not hasdataref and "RatioPlotYLabel" not in plot: if plot.get('RatioPlotMode', '') == 'deviation': plot['RatioPlotYLabel'] = 'Deviation' #r'$\text{MC}-\text{MC}_\text{ref}$' else: plot['RatioPlotYLabel'] = 'Ratio' #r'$\text{MC}/\text{MC}_\text{ref}$' ## Make the output and write to file if(len(anaobjects)>0) : o = mkOutput(hpath, anaobjects, plot, special) writeOutput(o, hpath)