diff --git a/analyses/pluginATLAS/ATLAS_2013_I1216670.cc b/analyses/pluginATLAS/ATLAS_2013_I1216670.cc
--- a/analyses/pluginATLAS/ATLAS_2013_I1216670.cc
+++ b/analyses/pluginATLAS/ATLAS_2013_I1216670.cc
@@ -1,120 +1,120 @@
 // -*- C++ -*-
 #include "Rivet/Analysis.hh"
 #include "Rivet/Projections/FinalState.hh"
 #include "Rivet/Projections/WFinder.hh"
 #include "Rivet/Projections/FastJets.hh"
 #include "Rivet/Projections/VetoedFinalState.hh"
 
 namespace Rivet {
 
   /// @brief MPI sensitive di-jet balance variables for W->ejj or W->mujj events.
   class ATLAS_2013_I1216670 : public Analysis {
   public:
 
     /// @name Constructor
     ATLAS_2013_I1216670()
       : Analysis("ATLAS_2013_I1216670")
     {    }
 
 
     /// @name Analysis methods
     //@{
 
     /// Book histograms, set up projections for W and jets
     void init() {
 
       book(_h_delta_jets_n ,1, 1, 1);
       book(_h_delta_jets   ,2, 1, 1);
 
       FinalState fs;
 
       Cut cuts = Cuts::abseta < 2.5 && Cuts::pT >= 20*GeV;
 
-      WFinder w_e_finder(fs, cuts, PID::ELECTRON, 40*GeV, MAXDOUBLE, 0.0*GeV, 0.0, WFinder::ClusterPhotons::NODECAY,
+      WFinder w_e_finder(fs, cuts, PID::ELECTRON, 40*GeV, DBL_MAX, 0.0*GeV, 0.0, WFinder::ClusterPhotons::NODECAY,
                          WFinder::AddPhotons::NO, WFinder::MassWindow::MT);
       declare(w_e_finder, "W_E_FINDER");
 
-      WFinder w_mu_finder(fs, cuts, PID::MUON, 40*GeV, MAXDOUBLE, 0.0*GeV, 0.0, WFinder::ClusterPhotons::NODECAY,
+      WFinder w_mu_finder(fs, cuts, PID::MUON, 40*GeV, DBL_MAX, 0.0*GeV, 0.0, WFinder::ClusterPhotons::NODECAY,
                           WFinder::AddPhotons::NO, WFinder::MassWindow::MT);
       declare(w_mu_finder, "W_MU_FINDER");
 
       VetoedFinalState jet_fs(fs);
       jet_fs.addVetoOnThisFinalState(getProjection<WFinder>("W_E_FINDER"));
       jet_fs.addVetoOnThisFinalState(getProjection<WFinder>("W_MU_FINDER"));
       FastJets jets(jet_fs, FastJets::ANTIKT, 0.4);
       declare(jets, "JETS");
 
     }
 
     /// Do the analysis
     void analyze(const Event &e) {
 
       const WFinder& w_e_finder  = apply<WFinder>(e, "W_E_FINDER" );
       const WFinder& w_mu_finder = apply<WFinder>(e, "W_MU_FINDER");
       Particle lepton, neutrino;
       Jets all_jets, jets;
 
       // Find exactly 1 W->e or W->mu boson
       if(w_e_finder.bosons().size() == 1 && w_mu_finder.bosons().size() == 0) {
         MSG_DEBUG(" Event identified as W->e nu.");
         if( !(w_e_finder.mT() > 40*GeV && w_e_finder.constituentNeutrino().Et() > 25.0*GeV) )  vetoEvent;
         lepton = w_e_finder.constituentLepton();
       } else if(w_mu_finder.bosons().size() == 1 && w_e_finder.bosons().size() == 0) {
         MSG_DEBUG(" Event identified as W->mu nu.");
         if( !(w_mu_finder.mT() > 40*GeV && w_mu_finder.constituentNeutrino().Et() > 25.0*GeV) )  vetoEvent;
         lepton = w_mu_finder.constituentLepton();
       } else {
         MSG_DEBUG(" No W found passing cuts.");
         vetoEvent;
       }
 
       all_jets = apply<FastJets>(e, "JETS").jetsByPt(Cuts::pt>20.0*GeV && Cuts::absrap<2.8);
 
       // Remove jets DeltaR < 0.5 from W lepton
       for(Jets::iterator it = all_jets.begin(); it != all_jets.end(); ++it) {
         double distance = deltaR( lepton, (*it) );
         if(distance < 0.5) {
           MSG_DEBUG("   Veto jet DeltaR " << distance << " from W lepton");
         } else {
           jets.push_back(*it);
         }
       }
 
       // Exactly two jets required
       if( jets.size() != 2 )  vetoEvent;
 
       // Calculate analysis quantities from the two jets
       double delta_jets = (jets.front().momentum() + jets.back().momentum()).pT();
       double total_pt = jets.front().momentum().pT() + jets.back().momentum().pT();
       double delta_jets_n = delta_jets / total_pt;
 
       _h_delta_jets->fill(   delta_jets); // Jet pT balance
       _h_delta_jets_n->fill( delta_jets_n); // Jet pT balance, normalised by scalar dijet pT
 
     }
 
     /// Finalize
     void finalize() {
 
       // Data is normalised to 0.03 and 3
       normalize(_h_delta_jets_n, 0.03);
       normalize(_h_delta_jets  , 3.0 );
 
     }
 
     //@}
 
   private:
 
     /// @name Histograms
     //@{
     Histo1DPtr _h_delta_jets_n;
     Histo1DPtr _h_delta_jets;
     //@}
 
   };
 
   // The hook for the plugin system
   DECLARE_RIVET_PLUGIN(ATLAS_2013_I1216670);
 
 }
diff --git a/analyses/pluginATLAS/ATLAS_2013_I1219109.cc b/analyses/pluginATLAS/ATLAS_2013_I1219109.cc
--- a/analyses/pluginATLAS/ATLAS_2013_I1219109.cc
+++ b/analyses/pluginATLAS/ATLAS_2013_I1219109.cc
@@ -1,157 +1,157 @@
 // -*- C++ -*-
 #include "Rivet/Analysis.hh"
 #include "Rivet/Projections/FinalState.hh"
 #include "Rivet/Projections/WFinder.hh"
 #include "Rivet/Projections/VetoedFinalState.hh"
 #include "Rivet/Projections/FastJets.hh"
 #include "Rivet/Projections/HeavyHadrons.hh"
 
 namespace Rivet {
 
 
   /// @brief ATLAS W+b measurement
   class ATLAS_2013_I1219109: public Analysis {
   public:
 
     ATLAS_2013_I1219109(string name = "ATLAS_2013_I1219109")
       : Analysis(name)
     {
       // the electron mode is used by default
       _mode = 1;
     }
 
 
     void init() {
       FinalState fs;
       declare(fs, "FinalState");
 
       Cut cuts = Cuts::abseta < 2.5 && Cuts::pT >= 25*GeV;
 
       // W finder for electrons and muons
-      WFinder wf(fs, cuts, _mode==3? PID::MUON : PID::ELECTRON, 0.0*GeV, MAXDOUBLE, 0.0, 0.1,
+      WFinder wf(fs, cuts, _mode==3? PID::MUON : PID::ELECTRON, 0.0*GeV, DBL_MAX, 0.0, 0.1,
                  WFinder::ClusterPhotons::NODECAY, WFinder::AddPhotons::NO, WFinder::MassWindow::MT);
       declare(wf, "WF");
 
       // jets
       VetoedFinalState jet_fs(fs);
       jet_fs.addVetoOnThisFinalState(getProjection<WFinder>("WF"));
       FastJets fj(jet_fs, FastJets::ANTIKT, 0.4);
       fj.useInvisibles();
       declare(fj, "Jets");
       declare(HeavyHadrons(Cuts::abseta < 2.5 && Cuts::pT > 5*GeV), "BHadrons");
 
 
       // book histograms
       book(_njet     ,1, 1, _mode); // dSigma / dNjet
       book(_jet1_bPt ,2, 1, _mode); // dSigma / dBjetPt for Njet = 1
       book(_jet2_bPt ,2, 2, _mode); // dSigma / dBjetPt for Njet = 2
 
     }
 
 
     void analyze(const Event& event) {
 
       //  retrieve W boson candidate
       const WFinder& wf = apply<WFinder>(event, "WF");
       if( wf.bosons().size() != 1 )  vetoEvent; // only one W boson candidate
       if( !(wf.mT() > 60.0*GeV) )    vetoEvent;
       //const Particle& Wboson  = wf.boson();
 
 
       // retrieve constituent neutrino
       const Particle& neutrino = wf.constituentNeutrino();
       if( !(neutrino.pT() > 25.0*GeV) )  vetoEvent;
 
       // retrieve constituent lepton
       const Particle& lepton = wf.constituentLepton();
 
       // count good jets, check if good jet contains B hadron
       const Particles& bHadrons = apply<HeavyHadrons>(event, "BHadrons").bHadrons();
       const Jets& jets = apply<JetAlg>(event, "Jets").jetsByPt(25*GeV);
       int goodjets = 0, bjets = 0;
       double bPt = 0.;
       for(const Jet& j : jets) {
         if( (j.abseta() < 2.1) && (deltaR(lepton, j) > 0.5) ) {
           // this jet passes the selection!
           ++goodjets;
           // j.bTagged() uses ghost association which is
           // more elegant, but not what has been used in
           // this analysis originally, will match B had-
           // rons in eta-phi space instead
           for(const Particle& b : bHadrons) {
             if( deltaR(j, b) < 0.3 ) {
               // jet matched to B hadron!
               if(!bPt)  bPt = j.pT() * GeV; // leading b-jet pT
               ++bjets; // count number of b-jets
               break;
             }
           }
         }
       }
       if( goodjets > 2 )  vetoEvent; // at most two jets
       if( !bjets )  vetoEvent; // at least one of them b-tagged
 
       double njets = double(goodjets);
       double ncomb = 3.0;
       _njet->fill(njets);
       _njet->fill(ncomb);
 
       if(     goodjets == 1)  _jet1_bPt->fill(bPt);
       else if(goodjets == 2)  _jet2_bPt->fill(bPt);
     }
 
 
     void finalize() {
 
       // Print summary info
       const double xs_pb(crossSection() / picobarn);
       const double sumw(sumOfWeights());
       MSG_INFO("Cross-Section/pb: " << xs_pb      );
       MSG_INFO("Sum of weights  : " << sumw       );
       MSG_INFO("nEvents         : " << numEvents());
 
       const double sf(xs_pb / sumw);
 
       scale(_njet,     sf);
       scale(_jet1_bPt, sf);
       scale(_jet2_bPt, sf);
     }
 
   protected:
 
     size_t _mode;
 
   private:
 
     Histo1DPtr _njet;
     Histo1DPtr _jet1_bPt;
     Histo1DPtr _jet2_bPt;
 
     //bool _isMuon;
 
   };
 
   class ATLAS_2013_I1219109_EL : public ATLAS_2013_I1219109 {
   public:
     ATLAS_2013_I1219109_EL()
       : ATLAS_2013_I1219109("ATLAS_2013_I1219109_EL")
     {
       _mode = 2;
     }
   };
 
   class ATLAS_2013_I1219109_MU : public ATLAS_2013_I1219109 {
   public:
     ATLAS_2013_I1219109_MU()
       : ATLAS_2013_I1219109("ATLAS_2013_I1219109_MU")
     {
       _mode = 3;
     }
   };
 
   // The hook for the plugin system
   DECLARE_RIVET_PLUGIN(ATLAS_2013_I1219109);
   DECLARE_RIVET_PLUGIN(ATLAS_2013_I1219109_EL);
   DECLARE_RIVET_PLUGIN(ATLAS_2013_I1219109_MU);
 
 }
diff --git a/analyses/pluginATLAS/ATLAS_2014_I1298023.cc b/analyses/pluginATLAS/ATLAS_2014_I1298023.cc
--- a/analyses/pluginATLAS/ATLAS_2014_I1298023.cc
+++ b/analyses/pluginATLAS/ATLAS_2014_I1298023.cc
@@ -1,126 +1,126 @@
 // -*- C++ -*-
 #include "Rivet/Analysis.hh"
 #include "Rivet/Projections/ChargedLeptons.hh"
 #include "Rivet/Projections/VetoedFinalState.hh"
 #include "Rivet/Projections/MissingMomentum.hh"
 #include "Rivet/Projections/FastJets.hh"
 #include "Rivet/Projections/DressedLeptons.hh"
 
 namespace Rivet {
 
   /// ATLAS same-sign WW production at 8 TeV (PRL)
   class ATLAS_2014_I1298023 : public Analysis {
   public:
 
     /// @name Constructors etc.
     //@{
     ATLAS_2014_I1298023()
       : Analysis("ATLAS_2014_I1298023")
     {
 
     }
     //@}
 
 
   public:
 
     /// @name Analysis methods
     //@{
 
     void init() {
 
       const FinalState fs;
 
       // bare leptons
       ChargedLeptons bare_leptons(fs);
 
       // dressed leptons
       Cut cuts = (Cuts::abseta < 2.5) & (Cuts::pT > 25*GeV);
       DressedLeptons leptons(fs, bare_leptons, 0.1, cuts);
       declare(leptons, "leptons");
 
       // MET
       declare(MissingMomentum(fs), "MissingET");
 
       // jets
       VetoedFinalState vfs(fs);
       vfs.addVetoPairId(PID::MUON);
       vfs.vetoNeutrinos();
       declare(FastJets(vfs, FastJets::ANTIKT, 0.4), "jets");
 
       // book histogram
       book(_hist ,1, 1, 1);
     }
 
 
     /// Perform the per-event analysis
     void analyze(const Event& event) {
       const vector<DressedLepton>& leptons = apply<DressedLeptons>(event, "leptons").dressedLeptons();
       if ( leptons.size() < 2 )  vetoEvent;
 
-      double minDR_ll = MAXDOUBLE, mll = -1.0;
+      double minDR_ll = DBL_MAX, mll = -1.0;
       for (unsigned int i = 0; i < leptons.size(); ++i) {
         DressedLepton lep1 = leptons[i];
         for (unsigned int j = i + 1; j < leptons.size(); ++j) {
           DressedLepton lep2 = leptons[j];
           double dr = deltaR(lep1, lep2);
           if ( dr < minDR_ll )  minDR_ll = dr;
           double m = (lep1.momentum() + lep2.momentum()).mass();
           if ( mll < 0. || m < mll )  mll = m;
         }
       }
       if ( minDR_ll <= 0.3 || mll <= 20*GeV )  vetoEvent;
 
       if ( leptons[0].charge() * leptons[1].charge() < 0.0 )  vetoEvent;
 
       const MissingMomentum& met = apply<MissingMomentum>(event, "MissingET");
       if ( met.vectorEt().mod() <= 40*GeV )  vetoEvent;
 
       const Jets& all_jets = apply<FastJets>(event, "jets").jetsByPt( (Cuts::abseta < 4.5) && (Cuts::pT > 30*GeV) );
       Jets jets;
-      double minDR_overall = MAXDOUBLE;
+      double minDR_overall = DBL_MAX;
       for (const Jet& jet : all_jets) {
-        double minDR_jet = MAXDOUBLE, minDR_electrons = MAXDOUBLE;
+        double minDR_jet = DBL_MAX, minDR_electrons = DBL_MAX;
         for( DressedLepton lep : leptons ) {
           double dr = deltaR(jet, lep);
           if ( dr < minDR_jet )  minDR_jet = dr;
           if ( lep.abspid() == 11 && dr < minDR_electrons )  minDR_electrons = dr;
         }
         if ( minDR_electrons < 0.05 )  continue; // veto jet if it overlaps with electron
         if ( minDR_jet < minDR_overall )  minDR_overall = minDR_jet;
         jets += jet;
       }
       if ( jets.size() < 2 || minDR_overall <= 0.3 )  vetoEvent;
 
       FourMomentum dijet = jets[0].momentum() + jets[1].momentum();
       if ( dijet.mass() <= 500*GeV )  vetoEvent;
 
       // inclusive region
       _hist->fill(0.5);
 
       // VBS region
       if ( deltaRap(jets[0], jets[1]) > 2.4 )  _hist->fill(1.5);
     }
 
     /// Normalise histograms etc., after the run
     void finalize() {
 
       const double scalefactor( crossSection() / sumOfWeights() );
       scale(_hist, scalefactor);
 
     }
     //@}
 
   private:
 
     /// @name Histograms
     //@{
     Histo1DPtr _hist;
     //@}
 
   };
 
   // The hook for the plugin system
   DECLARE_RIVET_PLUGIN(ATLAS_2014_I1298023);
 
 }
diff --git a/analyses/pluginATLAS/ATLAS_2014_I1306615.cc b/analyses/pluginATLAS/ATLAS_2014_I1306615.cc
--- a/analyses/pluginATLAS/ATLAS_2014_I1306615.cc
+++ b/analyses/pluginATLAS/ATLAS_2014_I1306615.cc
@@ -1,484 +1,484 @@
 // -*- C++ -*-
 #include "Rivet/Analysis.hh"
 #include "Rivet/Projections/DressedLeptons.hh"
 #include "Rivet/Projections/VetoedFinalState.hh"
 #include "Rivet/Projections/FastJets.hh"
 
 namespace Rivet {
 
 
   /// @brief ATLAS H->yy differential cross-sections measurement
   ///
   /// @author Michaela Queitsch-Maitland <michaela.queitsch-maitland@cern.ch>
   //
   // arXiv: http://arxiv.org/abs/ARXIV:1407.4222
   // HepData: http://hepdata.cedar.ac.uk/view/ins1306615
   class ATLAS_2014_I1306615 : public Analysis {
   public:
 
     // Constructor
     ATLAS_2014_I1306615()
       : Analysis("ATLAS_2014_I1306615")
     {    }
 
 
     // Book histograms and initialise projections before the run
     void init() {
 
       // Final state
       // All particles within |eta| < 5.0
       const FinalState FS(Cuts::abseta<5.0);
       declare(FS,"FS");
 
       // Project photons with pT > 25 GeV and |eta| < 2.37
       IdentifiedFinalState ph_FS(Cuts::abseta<2.37 && Cuts::pT>25.0*GeV);
       ph_FS.acceptIdPair(PID::PHOTON);
       declare(ph_FS, "PH_FS");
 
       // Project photons for dressing
       IdentifiedFinalState ph_dressing_FS(FS);
       ph_dressing_FS.acceptIdPair(PID::PHOTON);
 
       // Project bare electrons
       IdentifiedFinalState el_bare_FS(FS);
       el_bare_FS.acceptIdPair(PID::ELECTRON);
       declare(el_bare_FS,"el_bare_FS");
 
       // Project dressed electrons with pT > 15 GeV and |eta| < 2.47
       DressedLeptons el_dressed_FS(ph_dressing_FS, el_bare_FS, 0.1, Cuts::abseta < 2.47 && Cuts::pT > 15*GeV);
       declare(el_dressed_FS,"EL_DRESSED_FS");
 
       // Project bare muons
       IdentifiedFinalState mu_bare_FS(FS);
       mu_bare_FS.acceptIdPair(PID::MUON);
 
       // Project dressed muons with pT > 15 GeV and |eta| < 2.47
       //DressedLeptons mu_dressed_FS(ph_dressing_FS, mu_bare_FS, 0.1, true, -2.47, 2.47, 15.0*GeV, false);
       DressedLeptons mu_dressed_FS(ph_dressing_FS, mu_bare_FS, 0.1, Cuts::abseta < 2.47 && Cuts::pT > 15*GeV);
       declare(mu_dressed_FS,"MU_DRESSED_FS");
 
       // Final state excluding muons and neutrinos (for jet building and photon isolation)
       VetoedFinalState veto_mu_nu_FS(FS);
       veto_mu_nu_FS.vetoNeutrinos();
       veto_mu_nu_FS.addVetoPairId(PID::MUON);
       declare(veto_mu_nu_FS, "VETO_MU_NU_FS");
 
       // Build the anti-kT R=0.4 jets, using FinalState particles (vetoing muons and neutrinos)
       FastJets jets(veto_mu_nu_FS, FastJets::ANTIKT, 0.4);
       declare(jets, "JETS");
 
       // Book histograms
       // 1D distributions
       book(_h_pT_yy         ,1,1,1);
       book(_h_y_yy          ,2,1,1);
       book(_h_Njets30       ,3,1,1);
       book(_h_Njets50       ,4,1,1);
       book(_h_pT_j1         ,5,1,1);
       book(_h_y_j1          ,6,1,1);
       book(_h_HT            ,7,1,1);
       book(_h_pT_j2         ,8,1,1);
       book(_h_Dy_jj         ,9,1,1);
       book(_h_Dphi_yy_jj    ,10,1,1);
       book(_h_cosTS_CS      ,11,1,1);
       book(_h_cosTS_CS_5bin ,12,1,1);
       book(_h_Dphi_jj       ,13,1,1);
       book(_h_pTt_yy        ,14,1,1);
       book(_h_Dy_yy         ,15,1,1);
       book(_h_tau_jet       ,16,1,1);
       book(_h_sum_tau_jet   ,17,1,1);
       book(_h_y_j2          ,18,1,1);
       book(_h_pT_j3         ,19,1,1);
       book(_h_m_jj          ,20,1,1);
       book(_h_pT_yy_jj      ,21,1,1);
 
       // 2D distributions of cosTS_CS x pT_yy
       book(_h_cosTS_pTyy_low  ,22,1,1);
       book(_h_cosTS_pTyy_high ,22,1,2);
       book(_h_cosTS_pTyy_rest ,22,1,3);
 
       // 2D distributions of Njets x pT_yy
       book(_h_pTyy_Njets0 ,23,1,1);
       book(_h_pTyy_Njets1 ,23,1,2);
       book(_h_pTyy_Njets2 ,23,1,3);
 
       book(_h_pTj1_excl ,24,1,1);
 
       // Fiducial regions
       book(_h_fidXSecs ,29,1,1);
     }
 
     // Perform the per-event analysis
     void analyze(const Event& event) {
 
       // Get final state particles
       const ParticleVector& FS_ptcls          = apply<FinalState>(event, "FS").particles();
       const ParticleVector& ptcls_veto_mu_nu  = apply<VetoedFinalState>(event, "VETO_MU_NU_FS").particles();
       const ParticleVector& photons           = apply<IdentifiedFinalState>(event, "PH_FS").particlesByPt();
       const vector<DressedLepton>& el_dressed = apply<DressedLeptons>(event, "EL_DRESSED_FS").dressedLeptons();
       const vector<DressedLepton>& mu_dressed = apply<DressedLeptons>(event, "MU_DRESSED_FS").dressedLeptons();
 
       // For isolation calculation
       float dR_iso    = 0.4;
       float ETcut_iso = 14.0;
       FourMomentum ET_iso;
 
       // Fiducial selection: pT > 25 GeV, |eta| < 2.37 and isolation (in cone deltaR = 0.4) is < 14 GeV
       vector<const Particle*> fid_photons;
       for (const Particle& ph : photons) {
 
 	// Veto photons from hadron or tau decay
 	if ( fromHadronDecay(ph) ) continue;
 
 	// Calculate isolation
 	ET_iso = - ph.momentum();
 	// Loop over fs truth particles (excluding muons and neutrinos)
 	for (const Particle& p : ptcls_veto_mu_nu) {
 	  // Check if the truth particle is in a cone of 0.4
 	  if ( deltaR(ph.momentum(), p.momentum()) < dR_iso )
 	    ET_iso += p.momentum();
 	}
 
 	// Check isolation
 	if ( ET_iso.Et() > ETcut_iso ) continue;
 
 	// Fill vector of photons passing fiducial selection
 	fid_photons.push_back(&ph);
       }
 
       if(fid_photons.size() < 2) vetoEvent;
 
       const FourMomentum& y1 = fid_photons[0]->momentum();
       const FourMomentum& y2 = fid_photons[1]->momentum();
 
       double m_yy = (y1 + y2).mass();
 
       // Relative pT cuts
       if ( y1.pT() < 0.35 * m_yy || y2.pT() < 0.25 * m_yy ) vetoEvent;
 
       // Mass window cut
       if ( m_yy < 105 || m_yy > 160 ) vetoEvent;
 
       // -------------------------------------------- //
       // Passed diphoton baseline fiducial selection! //
       // -------------------------------------------- //
 
       // Electron selection
       vector<const Particle*> good_el;
       for(const DressedLepton& els : el_dressed)  {
 
 	const Particle& el = els.constituentLepton();
 	if ( el.momentum().pT()        < 15   ) continue;
 	if ( fabs(el.momentum().eta()) > 2.47 ) continue;
 	if ( deltaR(el.momentum(), y1) < 0.4  ) continue;
 	if ( deltaR(el.momentum(), y2) < 0.4  ) continue;
 	if ( fromHadronDecay(el)              ) continue; // Veto electrons from hadron or tau decay
 	good_el.push_back(&el);
       }
 
       // Muon selection
       vector<const Particle*> good_mu;
       for(const DressedLepton& mus : mu_dressed)  {
 
 	const Particle& mu = mus.constituentLepton();
 	if ( mu.momentum().pT()        < 15   ) continue;
 	if ( fabs(mu.momentum().eta()) > 2.47 ) continue;
 	if ( deltaR(mu.momentum(), y1) < 0.4  ) continue;
 	if ( deltaR(mu.momentum(), y2) < 0.4  ) continue;
 	if ( fromHadronDecay(mu)              ) continue; // Veto muons from hadron or tau decay
 	good_mu.push_back(&mu);
       }
 
       // Find prompt, invisible particles for missing ET calculation
       // Based on VisibleFinalState projection
       FourMomentum invisible(0,0,0,0);
       for (const Particle& p : FS_ptcls) {
 
 	// Veto non-prompt particles (from hadron or tau decay)
         if ( fromHadronDecay(p) ) continue;
 	// Charged particles are visible
 	if ( PID::charge3( p.pid() ) != 0 ) continue;
 	// Neutral hadrons are visible
 	if ( PID::isHadron( p.pid() ) ) continue;
         // Photons are visible
 	if ( p.pid() == PID::PHOTON ) continue;
         // Gluons are visible (for parton level analyses)
 	if ( p.pid() == PID::GLUON ) continue;
 	// Everything else is invisible
 	invisible += p.momentum();
       }
       double MET = invisible.Et();
 
       // Jet selection
       // Get jets with pT > 25 GeV and |rapidity| < 4.4
-      //const Jets& jets = apply<FastJets>(event, "JETS").jetsByPt(25.0*GeV, MAXDOUBLE, -4.4, 4.4, RAPIDITY);
+      //const Jets& jets = apply<FastJets>(event, "JETS").jetsByPt(25.0*GeV, DBL_MAX, -4.4, 4.4, RAPIDITY);
       const Jets& jets = apply<FastJets>(event, "JETS").jetsByPt(Cuts::pT>25.0*GeV && Cuts::absrap <4.4);
 
       vector<const Jet*> jets_25;
       vector<const Jet*> jets_30;
       vector<const Jet*> jets_50;
 
       for (const Jet& jet : jets) {
 
 	bool passOverlap = true;
 	// Overlap with leading photons
 	if ( deltaR(y1, jet.momentum()) < 0.4 ) passOverlap = false;
 	if ( deltaR(y2, jet.momentum()) < 0.4 ) passOverlap = false;
 
 	// Overlap with good electrons
 	for (const Particle* el : good_el)
 	  if ( deltaR(el->momentum(), jet.momentum()) < 0.2 ) passOverlap = false;
 
 	if ( ! passOverlap ) continue;
 
 	if ( fabs(jet.momentum().eta()) < 2.4 || ( fabs(jet.momentum().eta()) > 2.4 && jet.momentum().pT() > 30 ) ) jets_25.push_back(&jet);
 	if ( jet.momentum().pT() > 30 ) jets_30.push_back(&jet);
 	if ( jet.momentum().pT() > 50 ) jets_50.push_back(&jet);
       }
 
       // Fiducial regions
       _h_fidXSecs->fill(1);
       if ( jets_30.size() >= 1 ) _h_fidXSecs->fill(2);
       if ( jets_30.size() >= 2 ) _h_fidXSecs->fill(3);
       if ( jets_30.size() >= 3 ) _h_fidXSecs->fill(4);
       if ( jets_30.size() >= 2 && passVBFCuts(y1 + y2, jets_30.at(0)->momentum(), jets_30.at(1)->momentum()) ) _h_fidXSecs->fill(5);
       if ( (good_el.size() + good_mu.size()) > 0 ) _h_fidXSecs->fill(6);
       if ( MET > 80 ) _h_fidXSecs->fill(7);
 
       // Fill histograms
       // Inclusive variables
       _pT_yy    = (y1 + y2).pT();
       _y_yy     = fabs( (y1 + y2).rapidity() );
       _cosTS_CS = cosTS_CS(y1, y2);
       _pTt_yy   = pTt(y1, y2);
       _Dy_yy    = fabs( deltaRap(y1, y2) );
 
       _Njets30 = jets_30.size() > 3 ? 3 : jets_30.size();
       _Njets50 = jets_50.size() > 3 ? 3 : jets_50.size();
       _h_Njets30->fill(_Njets30);
       _h_Njets50->fill(_Njets50);
 
       _pT_j1 = jets_30.size() > 0 ? jets_30.at(0)->momentum().pT() : 0.;
       _pT_j2 = jets_30.size() > 1 ? jets_30.at(1)->momentum().pT() : 0.;
       _pT_j3 = jets_30.size() > 2 ? jets_30.at(2)->momentum().pT() : 0.;
 
       _HT = 0.0;
       for (const Jet* jet : jets_30)
 	_HT += jet->momentum().pT();
 
       _tau_jet     = tau_jet_max(y1 + y2, jets_25);
       _sum_tau_jet = sum_tau_jet(y1 + y2, jets_25);
 
       _h_pT_yy        ->fill(_pT_yy);
       _h_y_yy         ->fill(_y_yy);
       _h_pT_j1        ->fill(_pT_j1);
       _h_cosTS_CS     ->fill(_cosTS_CS);
       _h_cosTS_CS_5bin->fill(_cosTS_CS);
       _h_HT           ->fill(_HT);
       _h_pTt_yy       ->fill(_pTt_yy);
       _h_Dy_yy        ->fill(_Dy_yy);
       _h_tau_jet      ->fill(_tau_jet);
       _h_sum_tau_jet  ->fill(_sum_tau_jet);
 
       // >=1 jet variables
       if ( jets_30.size() >= 1 ) {
 	FourMomentum j1 = jets_30[0]->momentum();
 	_y_j1 = fabs( j1.rapidity() );
 
 	_h_pT_j2->fill(_pT_j2);
 	_h_y_j1 ->fill(_y_j1);
       }
 
       // >=2 jet variables
       if ( jets_30.size() >= 2 ) {
 	FourMomentum j1 = jets_30[0]->momentum();
 	FourMomentum j2 = jets_30[1]->momentum();
 
 	_Dy_jj      = fabs( deltaRap(j1, j2) );
 	_Dphi_jj    = fabs( deltaPhi(j1, j2) );
 	_Dphi_yy_jj = fabs( deltaPhi(y1 + y2, j1 + j2) );
 	_m_jj       = (j1 + j2).mass();
 	_pT_yy_jj   = (y1 + y2 + j1 + j2).pT();
 	_y_j2       = fabs( j2.rapidity() );
 
 	_h_Dy_jj      ->fill(_Dy_jj);
 	_h_Dphi_jj    ->fill(_Dphi_jj);
 	_h_Dphi_yy_jj ->fill(_Dphi_yy_jj);
 	_h_m_jj       ->fill(_m_jj);
 	_h_pT_yy_jj   ->fill(_pT_yy_jj);
 	_h_pT_j3      ->fill(_pT_j3);
 	_h_y_j2       ->fill(_y_j2);
       }
 
       // 2D distributions of cosTS_CS x pT_yy
       if ( _pT_yy < 80 )
 	_h_cosTS_pTyy_low->fill(_cosTS_CS);
       else if ( _pT_yy > 80 && _pT_yy < 200 )
 	_h_cosTS_pTyy_high->fill(_cosTS_CS);
       else if ( _pT_yy > 200 )
 	_h_cosTS_pTyy_rest->fill(_cosTS_CS);
 
       // 2D distributions of pT_yy x Njets
       if ( _Njets30 == 0 )
 	_h_pTyy_Njets0->fill(_pT_yy);
       else if ( _Njets30 == 1 )
 	_h_pTyy_Njets1->fill(_pT_yy);
       else if ( _Njets30 >= 2 )
 	_h_pTyy_Njets2->fill(_pT_yy);
 
       if ( _Njets30 == 1 ) _h_pTj1_excl->fill(_pT_j1);
 
     }
 
     // Normalise histograms after the run
     void finalize() {
 
       const double xs = crossSectionPerEvent()/femtobarn;
 
       scale(_h_pT_yy, xs);
       scale(_h_y_yy, xs);
       scale(_h_pT_j1, xs);
       scale(_h_y_j1, xs);
       scale(_h_HT, xs);
       scale(_h_pT_j2, xs);
       scale(_h_Dy_jj, xs);
       scale(_h_Dphi_yy_jj, xs);
       scale(_h_cosTS_CS, xs);
       scale(_h_cosTS_CS_5bin, xs);
       scale(_h_Dphi_jj, xs);
       scale(_h_pTt_yy, xs);
       scale(_h_Dy_yy, xs);
       scale(_h_tau_jet, xs);
       scale(_h_sum_tau_jet, xs);
       scale(_h_y_j2, xs);
       scale(_h_pT_j3, xs);
       scale(_h_m_jj, xs);
       scale(_h_pT_yy_jj, xs);
       scale(_h_cosTS_pTyy_low, xs);
       scale(_h_cosTS_pTyy_high, xs);
       scale(_h_cosTS_pTyy_rest, xs);
       scale(_h_pTyy_Njets0, xs);
       scale(_h_pTyy_Njets1, xs);
       scale(_h_pTyy_Njets2, xs);
       scale(_h_pTj1_excl, xs);
       scale(_h_Njets30, xs);
       scale(_h_Njets50, xs);
       scale(_h_fidXSecs, xs);
     }
 
     // Trace event record to see if particle came from a hadron (or a tau from a hadron decay)
     // Based on fromDecay() function
     bool fromHadronDecay(const Particle& p ) {
       const GenParticle * gp = p.genParticle();
       if (!gp) return false; /// TODO: something weird to make this necessary
       const GenVertex* prodVtx = p.genParticle()->production_vertex();
       if (prodVtx == NULL) return false;
       for (const GenParticle* ancestor : particles(prodVtx, HepMC::ancestors)) {
         const PdgId pid = ancestor->pdg_id();
         if (ancestor->status() == 2 && PID::isHadron(pid)) return true;
         if (ancestor->status() == 2 && (abs(pid) == PID::TAU && fromHadronDecay(ancestor))) return true;
       }
       return false;
     }
 
     // VBF-enhanced dijet topology selection cuts
     bool passVBFCuts(const FourMomentum &H, const FourMomentum &j1, const FourMomentum &j2) {
       return ( fabs(deltaRap(j1, j2)) > 2.8 && (j1 + j2).mass() > 400 && fabs(deltaPhi(H, j1 + j2)) > 2.6 );
     }
 
     // Cosine of the decay angle in the Collins-Soper frame
     double cosTS_CS(const FourMomentum &y1, const FourMomentum &y2) {
       return fabs( ( (y1.E() + y1.pz())* (y2.E() - y2.pz()) - (y1.E() - y1.pz()) * (y2.E() + y2.pz()) )
 		   / ((y1 + y2).mass() * sqrt(pow((y1 + y2).mass(), 2) + pow((y1 + y2).pt(), 2)) ) );
     }
 
     // Diphoton pT along thrust axis
     double pTt(const FourMomentum &y1, const FourMomentum &y2) {
       return fabs(y1.px() * y2.py() - y2.px() * y1.py()) / (y1 - y2).pT()*2;
     }
 
     // Tau of jet  (see paper for description)
     // tau_jet = mT/(2*cosh(y*)), where mT = pT (+) m, and y* = rapidty in Higgs rest frame
     double tau_jet( const FourMomentum &H, const FourMomentum &jet ) {
       return sqrt( pow(jet.pT(),2) + pow(jet.mass(),2) ) / (2.0 * cosh( jet.rapidity() - H.rapidity() ) );
     }
 
     // Maximal (leading) tau_jet (see paper for description)
     double tau_jet_max( const FourMomentum &H, const vector<const Jet*> jets, double tau_jet_cut = 8. ) {
       double max_tj = 0;
       for (size_t i=0; i < jets.size(); ++i) {
 	FourMomentum jet = jets[i]->momentum();
 	if ( tau_jet(H, jet) > tau_jet_cut )
 	  max_tj = max( tau_jet(H, jet), max_tj );
       }
       return max_tj;
     }
 
     // Scalar sum of tau for all jets (see paper for description)
     double sum_tau_jet( const FourMomentum &H, const vector<const Jet*> jets, double tau_jet_cut = 8. ) {
       double sum_tj = 0;
       for (size_t i=0; i < jets.size(); ++i) {
 	FourMomentum jet = jets[i]->momentum();
 	if ( tau_jet(H, jet) > tau_jet_cut )
 	  sum_tj += tau_jet(H, jet);
       }
       return sum_tj;
     }
 
   private:
 
     Histo1DPtr _h_pT_yy;
     Histo1DPtr _h_y_yy;
     Histo1DPtr _h_Njets30;
     Histo1DPtr _h_Njets50;
     Histo1DPtr _h_pT_j1;
     Histo1DPtr _h_y_j1;
     Histo1DPtr _h_HT;
     Histo1DPtr _h_pT_j2;
     Histo1DPtr _h_Dy_jj;
     Histo1DPtr _h_Dphi_yy_jj;
     Histo1DPtr _h_cosTS_CS;
     Histo1DPtr _h_cosTS_CS_5bin;
     Histo1DPtr _h_Dphi_jj;
     Histo1DPtr _h_pTt_yy;
     Histo1DPtr _h_Dy_yy;
     Histo1DPtr _h_tau_jet;
     Histo1DPtr _h_sum_tau_jet;
     Histo1DPtr _h_y_j2;
     Histo1DPtr _h_pT_j3;
     Histo1DPtr _h_m_jj;
     Histo1DPtr _h_pT_yy_jj;
     Histo1DPtr _h_cosTS_pTyy_low;
     Histo1DPtr _h_cosTS_pTyy_high;
     Histo1DPtr _h_cosTS_pTyy_rest;
     Histo1DPtr _h_pTyy_Njets0;
     Histo1DPtr _h_pTyy_Njets1;
     Histo1DPtr _h_pTyy_Njets2;
     Histo1DPtr _h_pTj1_excl;
     Histo1DPtr _h_fidXSecs;
 
     int _Njets30;
     int _Njets50;
     double _pT_yy;
     double _y_yy;
     double _cosTS_CS;
     double _pT_j1;
     double _m_jj;
     double _y_j1;
     double _HT;
     double _pT_j2;
     double _y_j2;
     double _Dphi_yy_jj;
     double _pT_yy_jj;
     double _Dphi_jj;
     double _Dy_jj;
     double _pT_j3;
     double _pTt_yy;
     double _Dy_yy;
     double _tau_jet;
     double _sum_tau_jet;
   };
 
   // The hook for the plugin system
   DECLARE_RIVET_PLUGIN(ATLAS_2014_I1306615);
 
 }
diff --git a/analyses/pluginATLAS/ATLAS_2014_I1319490.cc b/analyses/pluginATLAS/ATLAS_2014_I1319490.cc
--- a/analyses/pluginATLAS/ATLAS_2014_I1319490.cc
+++ b/analyses/pluginATLAS/ATLAS_2014_I1319490.cc
@@ -1,232 +1,232 @@
 #include "Rivet/Analysis.hh"
 #include "Rivet/Projections/FinalState.hh"
 #include "Rivet/Projections/WFinder.hh"
 #include "Rivet/Projections/FastJets.hh"
 #include "Rivet/Projections/VetoedFinalState.hh"
 
 namespace Rivet {
 
 
   class ATLAS_2014_I1319490 : public Analysis {
   public:
 
     ATLAS_2014_I1319490(string name = "ATLAS_2014_I1319490")
       : Analysis(name)
     {
       _mode = 0; // using electron channel for combined data by default
 
     }
 
 
     // Book histograms and initialise projections before the run
     void init() {
 
       FinalState fs;
 
       Cut cuts;
       if (_mode == 2) { // muon channel
         cuts = (Cuts::pT > 25.0*GeV) & Cuts::etaIn(-2.4, 2.4);
       } else if (_mode) { // electron channel
         cuts = (Cuts::pT > 25.0*GeV) & ( Cuts::etaIn(-2.47, -1.52) | Cuts::etaIn(-1.37, 1.37) | Cuts::etaIn(1.52, 2.47) );
       } else { // combined data extrapolated to common phase space
         cuts = (Cuts::pT > 25.0*GeV) & Cuts::etaIn(-2.5, 2.5);
       }
 
       // bosons
-      WFinder wfinder(fs, cuts, _mode > 1? PID::MUON : PID::ELECTRON, 40.0*GeV, MAXDOUBLE, 0.0*GeV, 0.1,
+      WFinder wfinder(fs, cuts, _mode > 1? PID::MUON : PID::ELECTRON, 40.0*GeV, DBL_MAX, 0.0*GeV, 0.1,
                       WFinder::ClusterPhotons::NODECAY, WFinder::AddPhotons::NO, WFinder::MassWindow::MT);
       declare(wfinder, "WF");
 
       // jets
       VetoedFinalState jet_fs(fs);
       jet_fs.addVetoOnThisFinalState(getProjection<WFinder>("WF"));
       FastJets jets(jet_fs, FastJets::ANTIKT, 0.4);
       jets.useInvisibles(true);
       declare(jets, "Jets");
 
       // book histograms
       book(histos["h_N_incl"]            ,1,1,_mode+1);
       book(histos["h_N"]                 ,4,1,_mode+1);
       book(histos["h_pt_jet1_1jet"]      ,5,1,_mode+1);
       book(histos["h_pt_jet1_1jet_excl"] ,6,1,_mode+1);
       book(histos["h_pt_jet1_2jet"]      ,7,1,_mode+1);
       book(histos["h_pt_jet1_3jet"]      ,8,1,_mode+1);
       book(histos["h_pt_jet2_2jet"]      ,9,1,_mode+1);
       book(histos["h_pt_jet3_3jet"]      ,10,1,_mode+1);
       book(histos["h_pt_jet4_4jet"]      ,11,1,_mode+1);
       book(histos["h_pt_jet5_5jet"]      ,12,1,_mode+1);
       book(histos["h_y_jet1_1jet"]       ,13,1,_mode+1);
       book(histos["h_y_jet2_2jet"]       ,14,1,_mode+1);
       book(histos["h_HT_1jet"]           ,15,1,_mode+1);
       book(histos["h_HT_1jet_excl"]      ,16,1,_mode+1);
       book(histos["h_HT_2jet"]           ,17,1,_mode+1);
       book(histos["h_HT_2jet_excl"]      ,18,1,_mode+1);
       book(histos["h_HT_3jet"]           ,19,1,_mode+1);
       book(histos["h_HT_3jet_excl"]      ,20,1,_mode+1);
       book(histos["h_HT_4jet"]           ,21,1,_mode+1);
       book(histos["h_HT_5jet"]           ,22,1,_mode+1);
       book(histos["h_deltaPhi_jet12"]    ,23,1,_mode+1);
       book(histos["h_deltaRap_jet12"]    ,24,1,_mode+1);
       book(histos["h_deltaR_jet12"]      ,25,1,_mode+1);
       book(histos["h_M_Jet12_2jet"]      ,26,1,_mode+1);
       book(histos["h_y_jet3_3jet"]       ,27,1,_mode+1);
       book(histos["h_y_jet4_4jet"]       ,28,1,_mode+1);
       book(histos["h_y_jet5_5jet"]       ,29,1,_mode+1);
       book(histos["h_ST_1jet"]           ,30,1,_mode+1);
       book(histos["h_ST_2jet"]           ,31,1,_mode+1);
       book(histos["h_ST_2jet_excl"]      ,32,1,_mode+1);
       book(histos["h_ST_3jet"]           ,33,1,_mode+1);
       book(histos["h_ST_3jet_excl"]      ,34,1,_mode+1);
       book(histos["h_ST_4jet"]           ,35,1,_mode+1);
       book(histos["h_ST_5jet"]           ,36,1,_mode+1);
     }
 
 
     void fillPlots(const Particle& lepton, const double& missET, Jets& all_jets) {
       // do jet-lepton overlap removal
       Jets jets;
       double ST = 0.0; // scalar pT sum of all selected jets
       for (const Jet &j : all_jets) {
         if (deltaR(j, lepton) > 0.5) {
           jets += j;
           ST += j.pT() / GeV;
         }
       }
 
       const size_t njets = jets.size();
 
       const double HT = ST + lepton.pT() / GeV + missET;
 
       histos["h_N"]->fill(njets + 0.5);
       for (size_t i = 0; i <= njets; ++i) {
         histos["h_N_incl"]->fill(i + 0.5);
       }
 
       if (njets) {
         const double pT1  = jets[0].pT() / GeV;
         const double rap1 = jets[0].absrap();
         histos["h_pt_jet1_1jet" ]->fill(pT1);
         histos["h_y_jet1_1jet"]->fill(rap1);
         histos["h_HT_1jet"]->fill(HT);
         histos["h_ST_1jet"]->fill(ST);
         if (njets == 1) {
           histos["h_pt_jet1_1jet_excl"]->fill(pT1);
           histos["h_HT_1jet_excl"]->fill(HT);
         } else {
           const double pT2  = jets[1].pT() / GeV;
           const double rap2 = jets[1].absrap();
           const double dR   = deltaR(jets[0], jets[1]);
           const double dRap = deltaRap(jets[0], jets[1]);
           const double dPhi = deltaPhi(jets[0], jets[1]);
           const double mjj  = (jets[0].momentum() + jets[1].momentum()).mass() / GeV;
           histos["h_pt_jet1_2jet"]->fill(pT1);
           histos["h_pt_jet2_2jet"]->fill(pT2);
           histos["h_y_jet2_2jet"]->fill(rap2);
           histos["h_M_Jet12_2jet"]->fill(mjj);
           histos["h_HT_2jet"]->fill(HT);
           histos["h_ST_2jet"]->fill(ST);
           histos["h_deltaPhi_jet12"]->fill(dPhi);
           histos["h_deltaRap_jet12"]->fill(dRap);
           histos["h_deltaR_jet12"]->fill(dR);
           if (njets == 2) {
             histos["h_ST_2jet_excl"]->fill(ST);
             histos["h_HT_2jet_excl"]->fill(HT);
           } else {
             const double pT3  = jets[2].pT() / GeV;
             const double rap3 = jets[2].absrap();
             histos["h_pt_jet1_3jet"]->fill(pT1);
             histos["h_pt_jet3_3jet"]->fill(pT3);
             histos["h_y_jet3_3jet"]->fill(rap3);
             histos["h_HT_3jet"]->fill(HT);
             histos["h_ST_3jet"]->fill(ST);
             if(njets == 3) {
               histos["h_ST_3jet_excl"]->fill(ST);
               histos["h_HT_3jet_excl"]->fill(HT);
             } else {
               const double pT4  = jets[3].pT() / GeV;
               const double rap4 = jets[3].absrap();
               histos["h_pt_jet4_4jet"]->fill(pT4);
               histos["h_y_jet4_4jet"]->fill(rap4);
               histos["h_HT_4jet"]->fill(HT);
               histos["h_ST_4jet"]->fill(ST);
               if (njets > 4) {
                 const double pT5  = jets[4].pT() / GeV;
                 const double rap5 = jets[4].absrap();
                 histos["h_pt_jet5_5jet"]->fill(pT5);
                 histos["h_y_jet5_5jet"]->fill(rap5);
                 histos["h_HT_5jet"]->fill(HT);
                 histos["h_ST_5jet"]->fill(ST);
               }
             }
           }
         }
       }
     }
 
 
     // Perform the per-event analysis
     void analyze(const Event& event) {
       // Retrieve boson candidate
       const WFinder& wf = apply<WFinder>(event, "WF");
       if (wf.empty()) vetoEvent;
 
       // Retrieve jets
       const JetAlg& jetfs = apply<JetAlg>(event, "Jets");
       Jets all_jets = jetfs.jetsByPt(Cuts::pT > 30.0*GeV && Cuts::absrap < 4.4);
 
       const Particles& leptons = wf.constituentLeptons();
       const double missET = wf.constituentNeutrino().pT() / GeV;
       if (leptons.size() == 1 && missET > 25.0 && wf.mT() > 40.0*GeV) {
         const Particle& lep = leptons[0];
         fillPlots(lep, missET, all_jets);
       }
     }
 
 
     void finalize() {
       const double scalefactor(crossSection() / sumOfWeights());
       /// @todo Update to use C++11 range-for
       for (map<string, Histo1DPtr>::iterator hit = histos.begin(); hit != histos.end(); ++hit) {
         scale(hit->second, scalefactor);
       }
     }
 
 
   protected:
 
     size_t _mode;
 
 
   private:
 
     map<string, Histo1DPtr> histos;
 
   };
 
 
   class ATLAS_2014_I1319490_EL : public ATLAS_2014_I1319490 {
   public:
     ATLAS_2014_I1319490_EL()
       : ATLAS_2014_I1319490("ATLAS_2014_I1319490_EL")
     {
       _mode = 1;
     }
   };
 
 
   class ATLAS_2014_I1319490_MU : public ATLAS_2014_I1319490 {
   public:
     ATLAS_2014_I1319490_MU()
       : ATLAS_2014_I1319490("ATLAS_2014_I1319490_MU")
     {
       _mode = 2;
     }
   };
 
 
   // The hooks for the plugin system
   DECLARE_RIVET_PLUGIN(ATLAS_2014_I1319490);
   DECLARE_RIVET_PLUGIN(ATLAS_2014_I1319490_EL);
   DECLARE_RIVET_PLUGIN(ATLAS_2014_I1319490_MU);
 
 }
diff --git a/analyses/pluginCMS/CMS_2013_I1218372.cc b/analyses/pluginCMS/CMS_2013_I1218372.cc
--- a/analyses/pluginCMS/CMS_2013_I1218372.cc
+++ b/analyses/pluginCMS/CMS_2013_I1218372.cc
@@ -1,162 +1,162 @@
 // Samantha Dooling DESY
 // February 2012
 //
 // -*- C++ -*-
 // =============================
 //
 // Ratio of the energy deposited in the pseudorapidity range
 // -6.6 < eta < -5.2 for events with a charged particle jet
 //
 // =============================
 #include "Rivet/Analysis.hh"
 #include "Rivet/Projections/FinalState.hh"
 #include "Rivet/Projections/ChargedFinalState.hh"
 #include "Rivet/Projections/FastJets.hh"
 #include "Rivet/Projections/Beam.hh"
 #include "Rivet/Projections/VetoedFinalState.hh"
 
 namespace Rivet {
 
 
   class CMS_2013_I1218372 : public Analysis {
   public:
 
   /// Constructor
   CMS_2013_I1218372()
     : Analysis("CMS_2013_I1218372")
     { }
 
     void init() {
 
       // gives the range of eta and min pT for the final state from which I get the jets
       FastJets jetpro (ChargedFinalState(-2.5, 2.5, 0.3*GeV), FastJets::ANTIKT, 0.5);
       declare(jetpro, "Jets");
 
       // skip Neutrinos and Muons
       VetoedFinalState fsv(FinalState(-7.0, -4.0, 0.*GeV));
       fsv.vetoNeutrinos();
       fsv.addVetoPairId(PID::MUON);
       declare(fsv, "fsv");
 
       // for the hadron level selection
-      VetoedFinalState sfsv(FinalState(-MAXDOUBLE, MAXDOUBLE, 0.*GeV));
+      VetoedFinalState sfsv(FinalState(-DBL_MAX, DBL_MAX, 0.*GeV));
       sfsv.vetoNeutrinos();
       sfsv.addVetoPairId(PID::MUON);
       declare(sfsv, "sfsv");
 
       //counters
       book(passedSumOfWeights, "passedSumOfWeights");
       book(inclEflow, "inclEflow");
 
       // Temporary histograms to fill the energy flow for leading jet events.
       // Ratios are calculated in finalyze().
       int id = 0;
       if (fuzzyEquals(sqrtS()/GeV,  900, 1e-3)) id=1;
       if (fuzzyEquals(sqrtS()/GeV, 2760, 1e-3)) id=2;
       if (fuzzyEquals(sqrtS()/GeV, 7000, 1e-3)) id=3;
       book(_h_ratio, id, 1, 1);
       book(_tmp_jet , "TMP/eflow_jet"  ,refData(id, 1, 1));  // Leading jet energy flow in pt
       book(_tmp_njet, "TMP/number_jet" ,refData(id, 1, 1)); // Number of events in pt
     }
 
 
     /// Perform the per-event analysis
     void analyze(const Event& event) {
 
       // Skip if the event is empty
       const FinalState& fsv = apply<FinalState>(event, "fsv");
       if (fsv.empty()) vetoEvent;
 
       // ====================== Minimum Bias selection
 
       const FinalState& sfsv = apply<FinalState>(event, "sfsv");
       Particles parts = sfsv.particles(cmpMomByRap);
       if (parts.empty()) vetoEvent;
 
       // find dymax
       double dymax = 0;
       int gap_pos  = -1;
       for (size_t i = 0; i < parts.size()-1; ++i) {
         double dy = parts[i+1].rapidity() - parts[i].rapidity();
         if (dy > dymax) {
           dymax = dy;
           gap_pos = i;
         }
       }
 
       // calculate mx2 and my2
       FourMomentum xmom;
       for (int i=0; i<=gap_pos; ++i) {
         xmom += parts[i].momentum();
       }
       double mx2 = xmom.mass2();
       if (mx2<0) vetoEvent;
 
       FourMomentum ymom;
       for (size_t i=gap_pos+1; i<parts.size(); ++i) {
         ymom += parts[i].momentum();
       }
       double my2 = ymom.mass2();
       if (my2<0) vetoEvent;
 
       // calculate xix and xiy and xidd
       double xix  = mx2 / sqr(sqrtS());
       double xiy  = my2 / sqr(sqrtS());
       double xidd = mx2*my2 / sqr(sqrtS()*0.938*GeV);
 
       // combine the selection: xi cuts
       bool passedHadronCuts = false;
       if (fuzzyEquals(sqrtS()/GeV,  900, 1e-3) && (xix > 0.1  || xiy > 0.4 || xidd > 0.5)) passedHadronCuts = true;
       if (fuzzyEquals(sqrtS()/GeV, 2760, 1e-3) && (xix > 0.07 || xiy > 0.2 || xidd > 0.5)) passedHadronCuts = true;
       if (fuzzyEquals(sqrtS()/GeV, 7000, 1e-3) && (xix > 0.04 || xiy > 0.1 || xidd > 0.5)) passedHadronCuts = true;
       if (!passedHadronCuts) vetoEvent;
 
       //  ============================== MINIMUM BIAS EVENTS
 
       // loop over particles to calculate the energy
       passedSumOfWeights->fill();
 
       for (const Particle& p : fsv.particles()) {
         if (-5.2 > p.eta() && p.eta() > -6.6) inclEflow->fill(p.E()/GeV);
       }
 
       //  ============================== JET EVENTS
 
       const FastJets& jetpro = apply<FastJets>(event, "Jets");
       const Jets& jets = jetpro.jetsByPt(1.0*GeV);
       if (jets.size()<1) vetoEvent;
 
       if (fabs(jets[0].eta()) < 2.0) {
         _tmp_njet->fill(jets[0].pT()/GeV);
 
         // energy flow
         for (const Particle& p : fsv.particles()) {
           if (p.eta() > -6.6 && p.eta() < -5.2) {  // ask for the CASTOR region
             _tmp_jet->fill(jets[0].pT()/GeV, p.E()/GeV);
           }
         }
       }
 
     }// analysis
 
     void finalize() {
       scale(_tmp_jet, *passedSumOfWeights / *inclEflow);
       divide(_tmp_jet, _tmp_njet, _h_ratio);
     }
 
   private:
     // counters
     CounterPtr passedSumOfWeights;
     CounterPtr inclEflow;
 
     // histograms
     Scatter2DPtr _h_ratio;
     Histo1DPtr   _tmp_jet;
     Histo1DPtr   _tmp_njet;
   };
 
 
   // The hook for the plugin system
   DECLARE_RIVET_PLUGIN(CMS_2013_I1218372);
 
 }
diff --git a/analyses/pluginCMS/CMS_2015_I1370682.cc b/analyses/pluginCMS/CMS_2015_I1370682.cc
--- a/analyses/pluginCMS/CMS_2015_I1370682.cc
+++ b/analyses/pluginCMS/CMS_2015_I1370682.cc
@@ -1,604 +1,604 @@
 #include "Rivet/Analysis.hh"
 #include "Rivet/Math/LorentzTrans.hh"
 #include "Rivet/Projections/FinalState.hh"
 #include "Rivet/Projections/FastJets.hh"
 
 namespace Rivet {
 
 
   namespace { //< only visible in this compilation unit
 
     /// @brief Pseudo top finder
     ///
     /// Find top quark in the particle level.
     /// The definition is based on the agreement at the LHC working group.
     class PseudoTop : public FinalState {
     public:
       /// @name Standard constructors and destructors.
       //@{
 
       /// The default constructor. May specify the minimum and maximum
       /// pseudorapidity \f$ \eta \f$ and the min \f$ p_T \f$ (in GeV).
       PseudoTop(double lepR = 0.1, double lepMinPt = 20, double lepMaxEta = 2.4,
                 double jetR = 0.4, double jetMinPt = 30, double jetMaxEta = 4.7)
-        : FinalState(-MAXDOUBLE, MAXDOUBLE, 0*GeV),
+        : FinalState(-DBL_MAX, DBL_MAX, 0*GeV),
           _lepR(lepR), _lepMinPt(lepMinPt), _lepMaxEta(lepMaxEta),
           _jetR(jetR), _jetMinPt(jetMinPt), _jetMaxEta(jetMaxEta)
       {
         setName("PseudoTop");
       }
 
       enum TTbarMode {CH_NONE=-1, CH_FULLHADRON = 0, CH_SEMILEPTON, CH_FULLLEPTON};
       enum DecayMode {CH_HADRON = 0, CH_MUON, CH_ELECTRON};
 
       TTbarMode mode() const {
         if (!_isValid) return CH_NONE;
         if (_mode1 == CH_HADRON && _mode2 == CH_HADRON) return CH_FULLHADRON;
         else if ( _mode1 != CH_HADRON && _mode2 != CH_HADRON) return CH_FULLLEPTON;
         else return CH_SEMILEPTON;
       }
       DecayMode mode1() const {return _mode1;}
       DecayMode mode2() const {return _mode2;}
 
       /// Clone on the heap.
       virtual unique_ptr<Projection> clone() const {
         return unique_ptr<Projection>(new PseudoTop(*this));
       }
 
       //@}
 
     public:
       Particle t1() const {return _t1;}
       Particle t2() const {return _t2;}
       Particle b1() const {return _b1;}
       Particle b2() const {return _b2;}
       ParticleVector wDecays1() const {return _wDecays1;}
       ParticleVector wDecays2() const {return _wDecays2;}
       Jets jets() const {return _jets;}
       Jets bjets() const {return _bjets;}
       Jets ljets() const {return _ljets;}
 
     protected:
       // Apply the projection to the event
       void project(const Event& e); // override; ///< @todo Re-enable when C++11 allowed
       void cleanup(std::map<double, std::pair<size_t, size_t> >& v, const bool doCrossCleanup=false) const;
 
     private:
       const double _lepR, _lepMinPt, _lepMaxEta;
       const double _jetR, _jetMinPt, _jetMaxEta;
 
       //constexpr ///< @todo Re-enable when C++11 allowed
       static double _tMass; // = 172.5*GeV; ///< @todo Re-enable when C++11 allowed
       //constexpr ///< @todo Re-enable when C++11 allowed
       static double _wMass; // = 80.4*GeV; ///< @todo Re-enable when C++11 allowed
 
     private:
       bool _isValid;
       DecayMode _mode1, _mode2;
 
       Particle _t1, _t2;
       Particle _b1, _b2;
       ParticleVector _wDecays1, _wDecays2;
       Jets _jets, _bjets, _ljets;
 
     };
 
     // More implementation below the analysis code
 
   }
 
 
 
   /// Pseudo-top analysis from CMS
   class CMS_2015_I1370682 : public Analysis {
   public:
 
     CMS_2015_I1370682()
       : Analysis("CMS_2015_I1370682"),
         _applyCorrection(true),
         _doShapeOnly(true)
     {    }
 
 
     void init() {
       declare(PseudoTop(0.1, 20, 2.4, 0.5, 30, 2.4), "ttbar");
 
       // Lepton + Jet channel
       book(_hSL_topPt         ,"d15-x01-y01"); // 1/sigma dsigma/dpt(top)
       book(_hSL_topPtTtbarSys ,"d16-x01-y01"); // 1/sigma dsigma/dpt*(top)
       book(_hSL_topY          ,"d17-x01-y01"); // 1/sigma dsigma/dy(top)
       book(_hSL_ttbarDelPhi   ,"d18-x01-y01"); // 1/sigma dsigma/ddeltaphi(t,tbar)
       book(_hSL_topPtLead     ,"d19-x01-y01"); // 1/sigma dsigma/dpt(t1)
       book(_hSL_topPtSubLead  ,"d20-x01-y01"); // 1/sigma dsigma/dpt(t2)
       book(_hSL_ttbarPt       ,"d21-x01-y01"); // 1/sigma dsigma/dpt(ttbar)
       book(_hSL_ttbarY        ,"d22-x01-y01"); // 1/sigma dsigma/dy(ttbar)
       book(_hSL_ttbarMass     ,"d23-x01-y01"); // 1/sigma dsigma/dm(ttbar)
 
       // Dilepton channel
       book(_hDL_topPt         ,"d24-x01-y01"); // 1/sigma dsigma/dpt(top)
       book(_hDL_topPtTtbarSys ,"d25-x01-y01"); // 1/sigma dsigma/dpt*(top)
       book(_hDL_topY          ,"d26-x01-y01"); // 1/sigma dsigma/dy(top)
       book(_hDL_ttbarDelPhi   ,"d27-x01-y01"); // 1/sigma dsigma/ddeltaphi(t,tbar)
       book(_hDL_topPtLead     ,"d28-x01-y01"); // 1/sigma dsigma/dpt(t1)
       book(_hDL_topPtSubLead  ,"d29-x01-y01"); // 1/sigma dsigma/dpt(t2)
       book(_hDL_ttbarPt       ,"d30-x01-y01"); // 1/sigma dsigma/dpt(ttbar)
       book(_hDL_ttbarY        ,"d31-x01-y01"); // 1/sigma dsigma/dy(ttbar)
       book(_hDL_ttbarMass     ,"d32-x01-y01"); // 1/sigma dsigma/dm(ttbar)
 
     }
 
 
     void analyze(const Event& event) {
 
       // Get the ttbar candidate
       const PseudoTop& ttbar = apply<PseudoTop>(event, "ttbar");
       if ( ttbar.mode() == PseudoTop::CH_NONE ) vetoEvent;
 
       const FourMomentum& t1P4 = ttbar.t1().momentum();
       const FourMomentum& t2P4 = ttbar.t2().momentum();
       const double pt1 = std::max(t1P4.pT(), t2P4.pT());
       const double pt2 = std::min(t1P4.pT(), t2P4.pT());
       const double dPhi = deltaPhi(t1P4, t2P4);
       const FourMomentum ttP4 = t1P4 + t2P4;
       const FourMomentum t1P4AtCM = LorentzTransform::mkFrameTransformFromBeta(ttP4.betaVec()).transform(t1P4);
 
       const double weight = 1.0;
 
       if ( ttbar.mode() == PseudoTop::CH_SEMILEPTON ) {
         const Particle lCand1 = ttbar.wDecays1()[0]; // w1 dau0 is the lepton in the PseudoTop
         if (lCand1.pT() < 33*GeV || lCand1.abseta() > 2.1) vetoEvent;
         _hSL_topPt->fill(t1P4.pT(), weight);
         _hSL_topPt->fill(t2P4.pT(), weight);
         _hSL_topPtTtbarSys->fill(t1P4AtCM.pT(), weight);
         _hSL_topY->fill(t1P4.rapidity(), weight);
         _hSL_topY->fill(t2P4.rapidity(), weight);
         _hSL_ttbarDelPhi->fill(dPhi, weight);
         _hSL_topPtLead->fill(pt1, weight);
         _hSL_topPtSubLead->fill(pt2, weight);
         _hSL_ttbarPt->fill(ttP4.pT(), weight);
         _hSL_ttbarY->fill(ttP4.rapidity(), weight);
         _hSL_ttbarMass->fill(ttP4.mass(), weight);
       }
       else if ( ttbar.mode() == PseudoTop::CH_FULLLEPTON ) {
         const Particle lCand1 = ttbar.wDecays1()[0]; // dau0 are the lepton in the PseudoTop
         const Particle lCand2 = ttbar.wDecays2()[0]; // dau0 are the lepton in the PseudoTop
         if (lCand1.pT() < 20*GeV || lCand1.abseta() > 2.4) vetoEvent;
         if (lCand2.pT() < 20*GeV || lCand2.abseta() > 2.4) vetoEvent;
         _hDL_topPt->fill(t1P4.pT(), weight);
         _hDL_topPt->fill(t2P4.pT(), weight);
         _hDL_topPtTtbarSys->fill(t1P4AtCM.pT(), weight);
         _hDL_topY->fill(t1P4.rapidity(), weight);
         _hDL_topY->fill(t2P4.rapidity(), weight);
         _hDL_ttbarDelPhi->fill(dPhi, weight);
         _hDL_topPtLead->fill(pt1, weight);
         _hDL_topPtSubLead->fill(pt2, weight);
         _hDL_ttbarPt->fill(ttP4.pT(), weight);
         _hDL_ttbarY->fill(ttP4.rapidity(), weight);
         _hDL_ttbarMass->fill(ttP4.mass(), weight);
       }
 
     }
 
 
     void finalize() {
       if ( _applyCorrection ) {
         // Correction functions for TOP-12-028 paper, (parton bin height)/(pseudotop bin height)
         const double ch15[] = { 5.473609, 4.941048, 4.173346, 3.391191, 2.785644, 2.371346, 2.194161, 2.197167, };
         const double ch16[] = { 5.470905, 4.948201, 4.081982, 3.225532, 2.617519, 2.239217, 2.127878, 2.185918, };
         const double ch17[] = { 10.003667, 4.546519, 3.828115, 3.601018, 3.522194, 3.524694, 3.600951, 3.808553, 4.531891, 9.995370, };
         const double ch18[] = { 4.406683, 4.054041, 3.885393, 4.213646, };
         const double ch19[] = { 6.182537, 5.257703, 4.422280, 3.568402, 2.889408, 2.415878, 2.189974, 2.173210, };
         const double ch20[] = { 5.199874, 4.693318, 3.902882, 3.143785, 2.607877, 2.280189, 2.204124, 2.260829, };
         const double ch21[] = { 6.053523, 3.777506, 3.562251, 3.601356, 3.569347, 3.410472, };
         const double ch22[] = { 11.932351, 4.803773, 3.782709, 3.390775, 3.226806, 3.218982, 3.382678, 3.773653, 4.788191, 11.905338, };
         const double ch23[] = { 7.145255, 5.637595, 4.049882, 3.025917, 2.326430, 1.773824, 1.235329, };
 
         const double ch24[] = { 2.268193, 2.372063, 2.323975, 2.034655, 1.736793, };
         const double ch25[] = { 2.231852, 2.383086, 2.341894, 2.031318, 1.729672, 1.486993, };
         const double ch26[] = { 3.993526, 2.308249, 2.075136, 2.038297, 2.036302, 2.078270, 2.295817, 4.017713, };
         const double ch27[] = { 2.205978, 2.175010, 2.215376, 2.473144, };
         const double ch28[] = { 2.321077, 2.371895, 2.338871, 2.057821, 1.755382, };
         const double ch29[] = { 2.222707, 2.372591, 2.301688, 1.991162, 1.695343, };
         const double ch30[] = { 2.599677, 2.026855, 2.138620, 2.229553, };
         const double ch31[] = { 5.791779, 2.636219, 2.103642, 1.967198, 1.962168, 2.096514, 2.641189, 5.780828, };
         const double ch32[] = { 2.006685, 2.545525, 2.477745, 2.335747, 2.194226, 2.076500, };
 
         applyCorrection(_hSL_topPt, ch15);
         applyCorrection(_hSL_topPtTtbarSys, ch16);
         applyCorrection(_hSL_topY, ch17);
         applyCorrection(_hSL_ttbarDelPhi, ch18);
         applyCorrection(_hSL_topPtLead, ch19);
         applyCorrection(_hSL_topPtSubLead, ch20);
         applyCorrection(_hSL_ttbarPt, ch21);
         applyCorrection(_hSL_ttbarY, ch22);
         applyCorrection(_hSL_ttbarMass, ch23);
 
         applyCorrection(_hDL_topPt, ch24);
         applyCorrection(_hDL_topPtTtbarSys, ch25);
         applyCorrection(_hDL_topY, ch26);
         applyCorrection(_hDL_ttbarDelPhi, ch27);
         applyCorrection(_hDL_topPtLead, ch28);
         applyCorrection(_hDL_topPtSubLead, ch29);
         applyCorrection(_hDL_ttbarPt, ch30);
         applyCorrection(_hDL_ttbarY, ch31);
         applyCorrection(_hDL_ttbarMass, ch32);
       }
 
       if ( _doShapeOnly ) {
         normalize(_hSL_topPt        );
         normalize(_hSL_topPtTtbarSys);
         normalize(_hSL_topY         );
         normalize(_hSL_ttbarDelPhi  );
         normalize(_hSL_topPtLead    );
         normalize(_hSL_topPtSubLead );
         normalize(_hSL_ttbarPt      );
         normalize(_hSL_ttbarY       );
         normalize(_hSL_ttbarMass    );
 
         normalize(_hDL_topPt        );
         normalize(_hDL_topPtTtbarSys);
         normalize(_hDL_topY         );
         normalize(_hDL_ttbarDelPhi  );
         normalize(_hDL_topPtLead    );
         normalize(_hDL_topPtSubLead );
         normalize(_hDL_ttbarPt      );
         normalize(_hDL_ttbarY       );
         normalize(_hDL_ttbarMass    );
       }
       else {
         const double s = 1./sumOfWeights();
         scale(_hSL_topPt        , s);
         scale(_hSL_topPtTtbarSys, s);
         scale(_hSL_topY         , s);
         scale(_hSL_ttbarDelPhi  , s);
         scale(_hSL_topPtLead    , s);
         scale(_hSL_topPtSubLead , s);
         scale(_hSL_ttbarPt      , s);
         scale(_hSL_ttbarY       , s);
         scale(_hSL_ttbarMass    , s);
         scale(_hDL_topPt        , s);
         scale(_hDL_topPtTtbarSys, s);
         scale(_hDL_topY         , s);
         scale(_hDL_ttbarDelPhi  , s);
         scale(_hDL_topPtLead    , s);
         scale(_hDL_topPtSubLead , s);
         scale(_hDL_ttbarPt      , s);
         scale(_hDL_ttbarY       , s);
         scale(_hDL_ttbarMass    , s);
       }
 
     }
 
 
     void applyCorrection(Histo1DPtr h, const double* cf) {
       vector<YODA::HistoBin1D>& bins = h->bins();
       for (size_t i=0, n=bins.size(); i<n; ++i ) {
         const double s = cf[i];
         YODA::HistoBin1D& bin = bins[i];
         bin.scaleW(s);
       }
     }
 
 
   private:
 
     const bool _applyCorrection, _doShapeOnly;
     Histo1DPtr _hSL_topPt, _hSL_topPtTtbarSys, _hSL_topY, _hSL_ttbarDelPhi, _hSL_topPtLead,
       _hSL_topPtSubLead, _hSL_ttbarPt, _hSL_ttbarY, _hSL_ttbarMass;
     Histo1DPtr _hDL_topPt, _hDL_topPtTtbarSys, _hDL_topY, _hDL_ttbarDelPhi, _hDL_topPtLead,
       _hDL_topPtSubLead, _hDL_ttbarPt, _hDL_ttbarY, _hDL_ttbarMass;
 
   };
 
 
 
   DECLARE_RIVET_PLUGIN(CMS_2015_I1370682);
 
 
   ///////////////
 
   // More PseudoTop implementation
   namespace {
 
 
     double PseudoTop::_tMass = 172.5*GeV;
     double PseudoTop::_wMass = 80.4*GeV;
 
 
     void PseudoTop::cleanup(map<double, pair<size_t, size_t> >& v, const bool doCrossCleanup) const {
       vector<map<double, pair<size_t, size_t> >::iterator> toErase;
       set<size_t> usedLeg1, usedLeg2;
       if ( !doCrossCleanup ) {
         /// @todo Reinstate when C++11 allowed: for (auto key = v.begin(); key != v.end(); ++key) {
         for (map<double, pair<size_t, size_t> >::iterator key = v.begin(); key != v.end(); ++key) {
           const size_t leg1 = key->second.first;
           const size_t leg2 = key->second.second;
           if (usedLeg1.find(leg1) == usedLeg1.end() and
               usedLeg2.find(leg2) == usedLeg2.end()) {
             usedLeg1.insert(leg1);
             usedLeg2.insert(leg2);
           } else {
             toErase.push_back(key);
           }
         }
       }
       else {
         /// @todo Reinstate when C++11 allowed: for (auto key = v.begin(); key != v.end(); ++key) {
         for (map<double, pair<size_t, size_t> >::iterator key = v.begin(); key != v.end(); ++key) {
           const size_t leg1 = key->second.first;
           const size_t leg2 = key->second.second;
           if (usedLeg1.find(leg1) == usedLeg1.end() and
               usedLeg1.find(leg2) == usedLeg1.end()) {
             usedLeg1.insert(leg1);
             usedLeg1.insert(leg2);
           } else {
             toErase.push_back(key);
           }
         }
       }
       /// @todo Reinstate when C++11 allowed:  for (auto& key : toErase) v.erase(key);
       for (size_t i = 0; i < toErase.size(); ++i) v.erase(toErase[i]);
     }
 
 
     void PseudoTop::project(const Event& e) {
       // Leptons : do the lepton clustering anti-kt R=0.1 using stable photons and leptons not from hadron decay
       // Neutrinos : neutrinos not from hadron decay
       // MET : vector sum of all invisible particles in x-y plane
       // Jets : anti-kt R=0.4 using all particles excluding neutrinos and particles used in lepton clustering
       //        add ghost B hadrons during the jet clustering to identify B jets.
 
       // W->lv : dressed lepton and neutrino pairs
       // W->jj : light flavored dijet
       // W candidate : select lv or jj pairs which minimise |mW1-80.4|+|mW2-80.4|
       //               lepton-neutrino pair will be selected with higher priority
 
       // t->Wb : W candidate + b jet
       // t candidate : select Wb pairs which minimise |mtop1-172.5|+|mtop2-172.5|
 
       _isValid = false;
       _theParticles.clear();
       _wDecays1.clear();
       _wDecays2.clear();
       _jets.clear();
       _bjets.clear();
       _ljets.clear();
       _mode1 = _mode2 = CH_HADRON;
 
       // Collect final state particles
       Particles pForLep, pForJet;
       Particles neutrinos; // Prompt neutrinos
       /// @todo Avoid this unsafe jump into HepMC -- all this can be done properly via VisibleFS and HeavyHadrons projections
       for (const GenParticle* p : Rivet::particles(e.genEvent())) {
         const int status = p->status();
         const int pdgId = p->pdg_id();
         if (status == 1) {
           Particle rp = *p;
           if (!PID::isHadron(pdgId) && !rp.fromHadron()) {
             // Collect particles not from hadron decay
             if (rp.isNeutrino()) {
               // Prompt neutrinos are kept in separate collection
               neutrinos.push_back(rp);
             } else if (pdgId == 22 || rp.isLepton()) {
               // Leptons and photons for the dressing
               pForLep.push_back(rp);
             }
           } else if (!rp.isNeutrino()) {
             // Use all particles from hadron decay
             pForJet.push_back(rp);
           }
         } else if (PID::isHadron(pdgId) && PID::hasBottom(pdgId)) {
           // NOTE: Consider B hadrons with pT > 5GeV - not in CMS proposal
           //if ( p->momentum().perp() < 5 ) continue;
 
           // Do unstable particles, to be used in the ghost B clustering
           // Use last B hadrons only
           bool isLast = true;
           for (const GenParticlePtr pp : Rivet::particles(p->end_vertex(), HepMC::children)) {
             if (PID::hasBottom(pp->pdg_id())) {
               isLast = false;
               break;
             }
           }
           if (!isLast) continue;
 
           // Rescale momentum by 10^-20
           Particle ghost(pdgId, FourMomentum(p->momentum())*1e-20/p->momentum().rho());
           pForJet.push_back(ghost);
         }
       }
 
       // Start object building from trivial thing - prompt neutrinos
       sortByPt(neutrinos);
 
       // Proceed to lepton dressing
       FastJets fjLep(FinalState(), FastJets::ANTIKT, _lepR);
       fjLep.calc(pForLep);
 
       Jets leptons;
       vector<int> leptonsId;
       set<int> dressedIdxs;
       for (const Jet& lep : fjLep.jetsByPt(_lepMinPt)) {
         if (lep.abseta() > _lepMaxEta) continue;
         double leadingPt = -1;
         int leptonId = 0;
         for (const Particle& p : lep.particles()) {
           /// @warning Barcodes aren't future-proof in HepMC
           dressedIdxs.insert(p.genParticle()->barcode());
           if (p.isLepton() && p.pT() > leadingPt) {
             leadingPt = p.pT();
             leptonId = p.pid();
           }
         }
         if (leptonId == 0) continue;
         leptons.push_back(lep);
         leptonsId.push_back(leptonId);
       }
 
       // Re-use particles not used in lepton dressing
       for (const Particle& rp : pForLep) {
         /// @warning Barcodes aren't future-proof in HepMC
         const int barcode = rp.genParticle()->barcode();
         // Skip if the particle is used in dressing
         if (dressedIdxs.find(barcode) != dressedIdxs.end()) continue;
         // Put back to be used in jet clustering
         pForJet.push_back(rp);
       }
 
       // Then do the jet clustering
       FastJets fjJet(FinalState(), FastJets::ANTIKT, _jetR);
       //fjJet.useInvisibles(); // NOTE: CMS proposal to remove neutrinos (AB: wouldn't work anyway, since they were excluded from clustering inputs)
       fjJet.calc(pForJet);
       for (const Jet& jet : fjJet.jetsByPt(_jetMinPt)) {
         if (jet.abseta() > _jetMaxEta) continue;
         _jets.push_back(jet);
         bool isBJet = false;
         for (const Particle& rp : jet.particles()) {
           if (PID::hasBottom(rp.pdgId())) {
             isBJet = true;
             break;
           }
         }
         if ( isBJet ) _bjets.push_back(jet);
         else _ljets.push_back(jet);
       }
 
       // Every building blocks are ready. Continue to pseudo-W and pseudo-top combination
 
       if (_bjets.size() < 2) return; // Ignore single top for now
       map<double, pair<size_t, size_t> > wLepCandIdxs;
       map<double, pair<size_t, size_t> > wHadCandIdxs;
 
       // Collect leptonic-decaying W's
       for (size_t iLep = 0, nLep = leptons.size(); iLep < nLep; ++iLep) {
         const Jet& lep = leptons.at(iLep);
         for (size_t iNu = 0, nNu = neutrinos.size(); iNu < nNu; ++iNu) {
           const Particle& nu = neutrinos.at(iNu);
           const double m = (lep.momentum()+nu.momentum()).mass();
           const double dm = std::abs(m-_wMass);
           wLepCandIdxs[dm] = make_pair(iLep, iNu);
         }
       }
 
       // Continue to hadronic decaying W's
       for (size_t i = 0, nLjet = _ljets.size(); i < nLjet; ++i) {
         const Jet& ljet1 = _ljets[i];
         for (size_t j = i+1; j < nLjet; ++j) {
           const Jet& ljet2 = _ljets[j];
           const double m = (ljet1.momentum()+ljet2.momentum()).mass();
           const double dm = std::abs(m-_wMass);
           wHadCandIdxs[dm] = make_pair(i, j);
         }
       }
 
       // Cleanup W candidate, choose pairs with minimum dm if they share decay products
       cleanup(wLepCandIdxs);
       cleanup(wHadCandIdxs, true);
       const size_t nWLepCand = wLepCandIdxs.size();
       const size_t nWHadCand = wHadCandIdxs.size();
 
       if (nWLepCand + nWHadCand < 2) return; // We skip single top
 
       int w1Q = 1, w2Q = -1;
       int w1dau1Id = 1, w2dau1Id = -1;
       FourMomentum w1dau1LVec, w1dau2LVec;
       FourMomentum w2dau1LVec, w2dau2LVec;
       if (nWLepCand == 0) { // Full hadronic case
         const pair<size_t, size_t>& idPair1 = wHadCandIdxs.begin()->second;
         const pair<size_t, size_t>& idPair2 = (++wHadCandIdxs.begin())->second;  ///< @todo Reinstate std::next
         const Jet& w1dau1 = _ljets[idPair1.first];
         const Jet& w1dau2 = _ljets[idPair1.second];
         const Jet& w2dau1 = _ljets[idPair2.first];
         const Jet& w2dau2 = _ljets[idPair2.second];
         w1dau1LVec = w1dau1.momentum();
         w1dau2LVec = w1dau2.momentum();
         w2dau1LVec = w2dau1.momentum();
         w2dau2LVec = w2dau2.momentum();
       } else if (nWLepCand == 1) { // Semi-leptonic case
         const pair<size_t, size_t>& idPair1 = wLepCandIdxs.begin()->second;
         const pair<size_t, size_t>& idPair2 = wHadCandIdxs.begin()->second;
         const Jet& w1dau1 = leptons[idPair1.first];
         const Particle& w1dau2 = neutrinos[idPair1.second];
         const Jet& w2dau1 = _ljets[idPair2.first];
         const Jet& w2dau2 = _ljets[idPair2.second];
         w1dau1LVec = w1dau1.momentum();
         w1dau2LVec = w1dau2.momentum();
         w2dau1LVec = w2dau1.momentum();
         w2dau2LVec = w2dau2.momentum();
         w1dau1Id = leptonsId[idPair1.first];
         w1Q = w1dau1Id > 0 ? -1 : 1;
         w2Q = -w1Q;
         switch (w1dau1Id) {
         case 13: case -13: _mode1 = CH_MUON; break;
         case 11: case -11: _mode1 = CH_ELECTRON; break;
         }
       } else { // Full leptonic case
         const pair<size_t, size_t>& idPair1 = wLepCandIdxs.begin()->second;
         const pair<size_t, size_t>& idPair2 = (++wLepCandIdxs.begin())->second;  ///< @todo Reinstate std::next
         const Jet& w1dau1 = leptons[idPair1.first];
         const Particle& w1dau2 = neutrinos[idPair1.second];
         const Jet& w2dau1 = leptons[idPair2.first];
         const Particle& w2dau2 = neutrinos[idPair2.second];
         w1dau1LVec = w1dau1.momentum();
         w1dau2LVec = w1dau2.momentum();
         w2dau1LVec = w2dau1.momentum();
         w2dau2LVec = w2dau2.momentum();
         w1dau1Id = leptonsId[idPair1.first];
         w2dau1Id = leptonsId[idPair2.first];
         w1Q = w1dau1Id > 0 ? -1 : 1;
         w2Q = w2dau1Id > 0 ? -1 : 1;
         switch (w1dau1Id) {
         case 13: case -13: _mode1 = CH_MUON; break;
         case 11: case -11: _mode1 = CH_ELECTRON; break;
         }
         switch (w2dau1Id) {
         case 13: case -13: _mode2 = CH_MUON; break;
         case 11: case -11: _mode2 = CH_ELECTRON; break;
         }
       }
       const FourMomentum w1LVec = w1dau1LVec+w1dau2LVec;
       const FourMomentum w2LVec = w2dau1LVec+w2dau2LVec;
 
       // Combine b jets
       double sumDm = 1e9;
       FourMomentum b1LVec, b2LVec;
       for (size_t i = 0, n = _bjets.size(); i < n; ++i) {
         const Jet& bjet1 = _bjets[i];
         const double mtop1 = (w1LVec+bjet1.momentum()).mass();
         const double dmtop1 = std::abs(mtop1-_tMass);
         for (size_t j=0; j<n; ++j) {
           if (i == j) continue;
           const Jet& bjet2 = _bjets[j];
           const double mtop2 = (w2LVec+bjet2.momentum()).mass();
           const double dmtop2 = std::abs(mtop2-_tMass);
 
           if (sumDm <= dmtop1+dmtop2) continue;
 
           sumDm = dmtop1+dmtop2;
           b1LVec = bjet1.momentum();
           b2LVec = bjet2.momentum();
         }
       }
       if (sumDm >= 1e9) return; // Failed to make top, but this should not happen.
 
       const FourMomentum t1LVec = w1LVec + b1LVec;
       const FourMomentum t2LVec = w2LVec + b2LVec;
 
       // Put all of them into candidate collection
       _t1 = Particle(w1Q*6, t1LVec);
       _b1 = Particle(w1Q*5, b1LVec);
       _wDecays1.push_back(Particle(w1dau1Id, w1dau1LVec));
       _wDecays1.push_back(Particle(-w1dau1Id+w1Q, w1dau2LVec));
 
       _t2 = Particle(w2Q*6, t2LVec);
       _b2 = Particle(w2Q*5, b2LVec);
       _wDecays2.push_back(Particle(w2dau1Id, w2dau1LVec));
       _wDecays2.push_back(Particle(-w2dau1Id+w2Q, w2dau2LVec));
 
       _isValid = true;
     }
 
   }
 
 
 }
diff --git a/analyses/pluginCMS/CMS_2016_I1413748.cc b/analyses/pluginCMS/CMS_2016_I1413748.cc
--- a/analyses/pluginCMS/CMS_2016_I1413748.cc
+++ b/analyses/pluginCMS/CMS_2016_I1413748.cc
@@ -1,328 +1,328 @@
 // -*- C++ -*-
 #include "Rivet/Analysis.hh"
 #include "Rivet/Projections/FinalState.hh"
 #include "Rivet/Projections/IdentifiedFinalState.hh"
 #include "Rivet/Projections/PromptFinalState.hh"
 #include "Rivet/Projections/DressedLeptons.hh"
 #include "Rivet/Projections/PartonicTops.hh"
 
 namespace Rivet {
 
 
   /// CMS 8 TeV dilepton channel ttbar spin correlations and polarisation analysis
   class CMS_2016_I1413748 : public Analysis {
   public:
 
     /// Constructor
     DEFAULT_RIVET_ANALYSIS_CTOR(CMS_2016_I1413748);
 
 
     /// Book histograms and initialise projections
     void init() {
 
       // Complete final state
-      FinalState fs(-MAXDOUBLE, MAXDOUBLE, 0*GeV);
+      FinalState fs(-DBL_MAX, DBL_MAX, 0*GeV);
 
       // Projection for dressed electrons and muons
       IdentifiedFinalState photons(fs);
       photons.acceptIdPair(PID::PHOTON);
 
       IdentifiedFinalState el_id(fs);
       el_id.acceptIdPair(PID::ELECTRON);
       PromptFinalState electrons(el_id);
       declare(electrons, "Electrons");
       DressedLeptons dressed_electrons(photons, electrons, 0.1);
       declare(dressed_electrons, "DressedElectrons");
 
       IdentifiedFinalState mu_id(fs);
       mu_id.acceptIdPair(PID::MUON);
       PromptFinalState muons(mu_id);
       declare(muons, "Muons");
       DressedLeptons dressed_muons(photons, muons, 0.1);
       declare(dressed_muons, "DressedMuons");
 
       // Parton-level top quarks
       declare(PartonicTops(PartonicTops::DecayMode::E_MU, false), "LeptonicPartonTops");
 
 
       // Booking of histograms
 
       // This histogram is independent of the parton-level information, and is an addition to the original analysis.
       // It is compared to the same data as the parton-level delta_phi histogram d02-x01-y01.
       book(_h_dphidressedleptons, "d00-x01-y01", _bins_dphi);
 
       // The remaining histos use parton-level information
       book(_h_dphi, "d02-x01-y01", _bins_dphi);
       book(_h_cos_opening_angle, "d05-x01-y01", _bins_cos_opening_angle);
       book(_h_c1c2, "d08-x01-y01", _bins_c1c2);
       book(_h_lep_costheta, "d11-x01-y01", _bins_lep_costheta);
       book(_h_lep_costheta_CPV, "d14-x01-y01", _bins_lep_costheta_CPV);
 
       // 2D histos
       book(_h_dphi_var[0], "d20-x01-y01", _bins_dphi, _bins_tt_mass);
       book(_h_cos_opening_angle_var[0], "d26-x01-y01", _bins_cos_opening_angle, _bins_tt_mass);
       book(_h_c1c2_var[0], "d32-x01-y01", _bins_c1c2, _bins_tt_mass);
       book(_h_lep_costheta_var[0], "d38-x01-y01", _bins_lep_costheta, _bins_tt_mass);
       book(_h_lep_costheta_CPV_var[0], "d44-x01-y01", _bins_lep_costheta_CPV, _bins_tt_mass);
 
       book(_h_dphi_var[1], "d50-x01-y01", _bins_dphi, _bins_tt_pT);
       book(_h_cos_opening_angle_var[1], "d56-x01-y01", _bins_cos_opening_angle, _bins_tt_pT);
       book(_h_c1c2_var[1], "d62-x01-y01", _bins_c1c2, _bins_tt_pT);
       book(_h_lep_costheta_var[1], "d68-x01-y01", _bins_lep_costheta, _bins_tt_pT);
       book(_h_lep_costheta_CPV_var[1], "d74-x01-y01", _bins_lep_costheta_CPV, _bins_tt_pT);
 
       book(_h_dphi_var[2], "d80-x01-y01", _bins_dphi, _bins_tt_absrapidity);
       book(_h_cos_opening_angle_var[2], "d86-x01-y01", _bins_cos_opening_angle, _bins_tt_absrapidity);
       book(_h_c1c2_var[2], "d92-x01-y01", _bins_c1c2, _bins_tt_absrapidity);
       book(_h_lep_costheta_var[2], "d98-x01-y01", _bins_lep_costheta, _bins_tt_absrapidity);
       book(_h_lep_costheta_CPV_var[2], "d104-x01-y01", _bins_lep_costheta_CPV, _bins_tt_absrapidity);
 
       // Profile histos for asymmetries
       book(_h_dphi_profile[0], "d17-x01-y01", _bins_tt_mass);
       book(_h_cos_opening_angle_profile[0], "d23-x01-y01", _bins_tt_mass);
       book(_h_c1c2_profile[0], "d29-x01-y01", _bins_tt_mass);
       book(_h_lep_costheta_profile[0], "d35-x01-y01", _bins_tt_mass);
       book(_h_lep_costheta_CPV_profile[0], "d41-x01-y01", _bins_tt_mass);
 
       book(_h_dphi_profile[1], "d47-x01-y01", _bins_tt_pT);
       book(_h_cos_opening_angle_profile[1], "d53-x01-y01", _bins_tt_pT);
       book(_h_c1c2_profile[1], "d59-x01-y01", _bins_tt_pT);
       book(_h_lep_costheta_profile[1], "d65-x01-y01", _bins_tt_pT);
       book(_h_lep_costheta_CPV_profile[1], "d71-x01-y01", _bins_tt_pT);
 
       book(_h_dphi_profile[2], "d77-x01-y01", _bins_tt_absrapidity);
       book(_h_cos_opening_angle_profile[2], "d83-x01-y01", _bins_tt_absrapidity);
       book(_h_c1c2_profile[2], "d89-x01-y01", _bins_tt_absrapidity);
       book(_h_lep_costheta_profile[2], "d95-x01-y01", _bins_tt_absrapidity);
       book(_h_lep_costheta_CPV_profile[2], "d101-x01-y01", _bins_tt_absrapidity);
 
     }
 
 
     /// Perform the per-event analysis
     void analyze(const Event& event) {
 
       const double weight = 1.0;
 
       // Use particle-level leptons for the first histogram
       const DressedLeptons& dressed_electrons = applyProjection<DressedLeptons>(event, "DressedElectrons");
       const DressedLeptons& dressed_muons = applyProjection<DressedLeptons>(event, "DressedMuons");
 
       const vector<DressedLepton> dressedels = dressed_electrons.dressedLeptons();
       const vector<DressedLepton> dressedmus = dressed_muons.dressedLeptons();
 
       const size_t ndressedel = dressedels.size();
       const size_t ndressedmu = dressedmus.size();
 
       // For the particle-level histogram, require exactly one electron and exactly one muon, to select
       // the ttbar->emu channel. Note this means ttbar->emu events with additional PromptFinalState
       // dilepton pairs from the shower are vetoed - for PYTHIA8, this affects ~0.5% of events, so the
       // effect is well below the level of sensitivity of the measured distribution.
       if ( ndressedel == 1 && ndressedmu == 1 ) {
 
         const int electrontouse = 0, muontouse = 0;
 
         // Opposite-charge leptons only
         if ( sameSign(dressedels[electrontouse],dressedmus[muontouse]) ) {
           MSG_INFO("Error, e and mu have same charge, skipping event");
         }
         else {
           //Get the four-momenta of the positively- and negatively-charged leptons
           FourMomentum lepPlus = dressedels[electrontouse].charge() > 0 ? dressedels[electrontouse] : dressedmus[muontouse];
           FourMomentum lepMinus = dressedels[electrontouse].charge() > 0 ? dressedmus[muontouse] : dressedels[electrontouse];
 
           // Now calculate the variable
           double dphi_temp = deltaPhi(lepPlus,lepMinus);
 
           fillWithUFOF( _h_dphidressedleptons, dphi_temp, weight );
         }
 
       }
 
 
       // The remaining variables use parton-level information.
 
       // Get the leptonically decaying tops
       const Particles& leptonicpartontops = apply<ParticleFinder>(event, "LeptonicPartonTops").particlesByPt();
       Particles chargedleptons;
       unsigned int ntrueleptonictops = 0;
       bool oppositesign = false;
 
       if ( leptonicpartontops.size() == 2 ) {
         for (size_t k = 0; k < leptonicpartontops.size(); ++k) {
 
           // Get the lepton
           const Particle lepTop = leptonicpartontops[k];
           const auto isPromptChargedLepton = [](const Particle& p){return (isChargedLepton(p) && isPrompt(p, false, false));};
           Particles lepton_candidates = lepTop.allDescendants(firstParticleWith(isPromptChargedLepton), false);
           if ( lepton_candidates.size() < 1 ) MSG_WARNING("error, PartonicTops::DecayMode::E_MU top quark had no daughter lepton candidate, skipping event.");
 
           // In some cases there is no lepton from the W decay but only leptons from the decay of a radiated gamma.
           // These hadronic PartonicTops are currently being mistakenly selected by PartonicTops::DecayMode::E_MU (as of April 2017), and need to be rejected.
           // PartonicTops::DecayMode::E_MU is being fixed in Rivet, and when it is the veto below should do nothing.
           /// @todo Should no longer be necessary -- remove
           bool istrueleptonictop = false;
           for (size_t i = 0; i < lepton_candidates.size(); ++i) {
             const Particle& lepton_candidate = lepton_candidates[i];
             if ( lepton_candidate.hasParent(PID::PHOTON) ) {
               MSG_DEBUG("Found gamma parent, top: " << k+1 << " of " << leptonicpartontops.size() << " , lepton: " << i+1 << " of " << lepton_candidates.size());
               continue;
             }
             if ( !istrueleptonictop && sameSign(lepTop,lepton_candidate) ) {
               chargedleptons.push_back(lepton_candidate);
               istrueleptonictop = true;
             }
             else MSG_WARNING("Found extra prompt charged lepton from top decay (and without gamma parent), ignoring it.");
           }
           if ( istrueleptonictop ) ++ntrueleptonictops;
         }
       }
 
       if ( ntrueleptonictops == 2 ) {
         oppositesign = !( sameSign(chargedleptons[0],chargedleptons[1]) );
         if ( !oppositesign ) MSG_WARNING("error, same charge tops, skipping event.");
       }
 
       if ( ntrueleptonictops == 2 && oppositesign ) {
 
         // Get the four-momenta of the positively- and negatively-charged leptons
         FourMomentum lepPlus = chargedleptons[0].charge() > 0 ? chargedleptons[0] : chargedleptons[1];
         FourMomentum lepMinus = chargedleptons[0].charge() > 0 ? chargedleptons[1] : chargedleptons[0];
 
         const double dphi_temp = deltaPhi(lepPlus,lepMinus);
 
         // Get the four-momenta of the positively- and negatively-charged tops
         FourMomentum topPlus_p4 = leptonicpartontops[0].pdgId() > 0 ? leptonicpartontops[0] : leptonicpartontops[1];
         FourMomentum topMinus_p4 = leptonicpartontops[0].pdgId() > 0 ? leptonicpartontops[1] : leptonicpartontops[0];
         const FourMomentum ttbar_p4 = topPlus_p4 + topMinus_p4;
 
         const double tt_mass_temp = ttbar_p4.mass();
         const double tt_absrapidity_temp = ttbar_p4.absrapidity();
         const double tt_pT_temp = ttbar_p4.pT();
 
         // Lorentz transformations to calculate the spin observables in the helicity basis
 
         // Transform everything to the ttbar CM frame
         LorentzTransform ttCM;
         ttCM.setBetaVec(-ttbar_p4.betaVec());
 
         topPlus_p4 = ttCM.transform(topPlus_p4);
         topMinus_p4 = ttCM.transform(topMinus_p4);
 
         lepPlus = ttCM.transform(lepPlus);
         lepMinus = ttCM.transform(lepMinus);
 
         // Now boost the leptons to their parent top CM frames
         LorentzTransform topPlus, topMinus;
         topPlus.setBetaVec(-topPlus_p4.betaVec());
         topMinus.setBetaVec(-topMinus_p4.betaVec());
 
         lepPlus = topPlus.transform(lepPlus);
         lepMinus = topMinus.transform(lepMinus);
 
         const double lepPlus_costheta_temp = lepPlus.vector3().dot(topPlus_p4.vector3()) / (lepPlus.vector3().mod() * topPlus_p4.vector3().mod());
         const double lepMinus_costheta_temp = lepMinus.vector3().dot(topMinus_p4.vector3()) / (lepMinus.vector3().mod() * topMinus_p4.vector3().mod());
         const double c1c2_temp = lepPlus_costheta_temp * lepMinus_costheta_temp;
         const double cos_opening_angle_temp = lepPlus.vector3().dot(lepMinus.vector3()) / (lepPlus.vector3().mod() * lepMinus.vector3().mod());
 
         // Fill parton-level histos
         fillWithUFOF( _h_dphi, dphi_temp, weight );
         fillWithUFOF( _h_cos_opening_angle, cos_opening_angle_temp, weight );
         fillWithUFOF( _h_c1c2, c1c2_temp, weight );
         fillWithUFOF( _h_lep_costheta, lepPlus_costheta_temp, weight );
         fillWithUFOF( _h_lep_costheta, lepMinus_costheta_temp, weight );
         fillWithUFOF( _h_lep_costheta_CPV, lepPlus_costheta_temp, weight );
         fillWithUFOF( _h_lep_costheta_CPV, -lepMinus_costheta_temp, weight );
 
         // Now fill the same variables in the 2D and profile histos vs ttbar invariant mass, pT, and absolute rapidity
         for (int i_var = 0; i_var < 3; ++i_var) {
           double var;
           if ( i_var == 0 ) {
             var = tt_mass_temp;
           } else if ( i_var == 1 ) {
             var = tt_pT_temp;
           } else {
             var = tt_absrapidity_temp;
           }
 
           fillWithUFOF( _h_dphi_var[i_var], dphi_temp, var, weight );
           fillWithUFOF( _h_cos_opening_angle_var[i_var], cos_opening_angle_temp, var, weight );
           fillWithUFOF( _h_c1c2_var[i_var], c1c2_temp, var, weight );
           fillWithUFOF( _h_lep_costheta_var[i_var], lepPlus_costheta_temp, var, weight );
           fillWithUFOF( _h_lep_costheta_var[i_var], lepMinus_costheta_temp, var, weight );
           fillWithUFOF( _h_lep_costheta_CPV_var[i_var], lepPlus_costheta_temp, var, weight );
           fillWithUFOF( _h_lep_costheta_CPV_var[i_var], -lepMinus_costheta_temp, var, weight );
 
           fillWithUFOF( _h_dphi_profile[i_var], dphi_temp, var, weight, (_h_dphi->xMax() + _h_dphi->xMin())/2. );
           fillWithUFOF( _h_cos_opening_angle_profile[i_var], cos_opening_angle_temp, var, weight, (_h_cos_opening_angle->xMax() + _h_cos_opening_angle->xMin())/2. );
           fillWithUFOF( _h_c1c2_profile[i_var], c1c2_temp, var, weight, (_h_c1c2->xMax() + _h_c1c2->xMin())/2. );
           fillWithUFOF( _h_lep_costheta_profile[i_var], lepPlus_costheta_temp, var, weight, (_h_lep_costheta->xMax() + _h_lep_costheta->xMin())/2. );
           fillWithUFOF( _h_lep_costheta_profile[i_var], lepMinus_costheta_temp, var, weight, (_h_lep_costheta->xMax() + _h_lep_costheta->xMin())/2. );
           fillWithUFOF( _h_lep_costheta_CPV_profile[i_var], lepPlus_costheta_temp, var, weight, (_h_lep_costheta_CPV->xMax() + _h_lep_costheta_CPV->xMin())/2. );
           fillWithUFOF( _h_lep_costheta_CPV_profile[i_var], -lepMinus_costheta_temp, var, weight, (_h_lep_costheta_CPV->xMax() + _h_lep_costheta_CPV->xMin())/2. );
 
         }
 
       }
 
     }
 
 
     /// Normalise histograms to unit area
     void finalize() {
 
       normalize(_h_dphidressedleptons);
 
       normalize(_h_dphi);
       normalize(_h_cos_opening_angle);
       normalize(_h_c1c2);
       normalize(_h_lep_costheta);
       normalize(_h_lep_costheta_CPV);
 
       for (int i_var = 0; i_var < 3; ++i_var) {
         normalize(_h_dphi_var[i_var]);
         normalize(_h_cos_opening_angle_var[i_var]);
         normalize(_h_c1c2_var[i_var]);
         normalize(_h_lep_costheta_var[i_var]);
         normalize(_h_lep_costheta_CPV_var[i_var]);
       }
 
     }
 
 
   private:
 
     Histo1DPtr _h_dphidressedleptons, _h_dphi, _h_lep_costheta, _h_lep_costheta_CPV, _h_c1c2, _h_cos_opening_angle;
     Histo2DPtr _h_dphi_var[3], _h_lep_costheta_var[3], _h_lep_costheta_CPV_var[3], _h_c1c2_var[3], _h_cos_opening_angle_var[3];
     Profile1DPtr _h_dphi_profile[3], _h_lep_costheta_profile[3], _h_lep_costheta_CPV_profile[3], _h_c1c2_profile[3], _h_cos_opening_angle_profile[3];
 
     const vector<double> _bins_tt_mass = {300., 430., 530., 1200.};
     const vector<double> _bins_tt_pT = {0., 41., 92., 300.};
     const vector<double> _bins_tt_absrapidity = {0., 0.34, 0.75, 1.5};
     const vector<double> _bins_dphi = {0., 5.*M_PI/60., 10.*M_PI/60., 15.*M_PI/60., 20.*M_PI/60., 25.*M_PI/60., 30.*M_PI/60., 35.*M_PI/60., 40.*M_PI/60., 45.*M_PI/60., 50.*M_PI/60., 55.*M_PI/60., M_PI};
     const vector<double> _bins_lep_costheta = {-1., -2./3., -1./3., 0., 1./3., 2./3., 1.};
     const vector<double> _bins_lep_costheta_CPV = {-1., -2./3., -1./3., 0., 1./3., 2./3., 1.};
     const vector<double> _bins_c1c2 = {-1., -0.4, -10./60., 0., 10./60., 0.4, 1.};
     const vector<double> _bins_cos_opening_angle = {-1., -2./3., -1./3., 0., 1./3., 2./3., 1.};
 
     void fillWithUFOF(Histo1DPtr h, double x, double w) {
       h->fill(std::max(std::min(x, h->xMax()-1e-9),h->xMin()+1e-9), w);
     }
 
     void fillWithUFOF(Histo2DPtr h, double x, double y, double w) {
       h->fill(std::max(std::min(x, h->xMax()-1e-9),h->xMin()+1e-9), std::max(std::min(y, h->yMax()-1e-9),h->yMin()+1e-9), w);
     }
 
     void fillWithUFOF(Profile1DPtr h, double x, double y, double w, double c) {
       h->fill(std::max(std::min(y, h->xMax()-1e-9),h->xMin()+1e-9), float(x > c) - float(x < c), w);
     }
 
 
   };
 
 
   // The hook for the plugin system
   DECLARE_RIVET_PLUGIN(CMS_2016_I1413748);
 
 
 }
diff --git a/analyses/pluginCMS/CMS_2016_I1430892.cc b/analyses/pluginCMS/CMS_2016_I1430892.cc
--- a/analyses/pluginCMS/CMS_2016_I1430892.cc
+++ b/analyses/pluginCMS/CMS_2016_I1430892.cc
@@ -1,259 +1,259 @@
 // -*- C++ -*-
 #include "Rivet/Analysis.hh"
 #include "Rivet/Projections/FinalState.hh"
 #include "Rivet/Projections/IdentifiedFinalState.hh"
 #include "Rivet/Projections/PromptFinalState.hh"
 #include "Rivet/Projections/DressedLeptons.hh"
 #include "Rivet/Projections/PartonicTops.hh"
 
 namespace Rivet {
 
 
   /// CMS 8 TeV dilepton channel ttbar charge asymmetry analysis
   class CMS_2016_I1430892 : public Analysis {
   public:
 
     /// Constructor
     DEFAULT_RIVET_ANALYSIS_CTOR(CMS_2016_I1430892);
 
 
     /// Book histograms and initialise projections
     void init() {
 
       // Complete final state
-      FinalState fs(-MAXDOUBLE, MAXDOUBLE, 0*GeV);
+      FinalState fs(-DBL_MAX, DBL_MAX, 0*GeV);
 
       // Projection for dressed electrons and muons
       IdentifiedFinalState photons(fs);
       photons.acceptIdPair(PID::PHOTON);
 
       IdentifiedFinalState el_id(fs);
       el_id.acceptIdPair(PID::ELECTRON);
       PromptFinalState electrons(el_id);
       declare(electrons, "Electrons");
       DressedLeptons dressed_electrons(photons, electrons, 0.1);
       declare(dressed_electrons, "DressedElectrons");
 
       IdentifiedFinalState mu_id(fs);
       mu_id.acceptIdPair(PID::MUON);
       PromptFinalState muons(mu_id);
       declare(muons, "Muons");
       DressedLeptons dressed_muons(photons, muons, 0.1);
       declare(dressed_muons, "DressedMuons");
 
       // Parton-level top quarks
       declare(PartonicTops(PartonicTops::DecayMode::E_MU, false), "LeptonicPartonTops");
 
 
       // Booking of histograms
 
       // This histogram is independent of the parton-level information, and is an
       // addition to the original analysis. It is compared to the same data as
       // the parton-level delta_abseta histogram d05-x01-y01.
       book(_h_dabsetadressedleptons, "d00-x01-y01", _bins_dabseta);
 
       // The remaining histos use parton-level information
       book(_h_dabseta, "d05-x01-y01", _bins_dabseta);
       book(_h_dabsrapidity, "d02-x01-y01", _bins_dabsrapidity);
 
       // 2D histos
       book(_h_dabsrapidity_var[0], "d11-x01-y01", _bins_dabsrapidity, _bins_tt_mass);
       book(_h_dabseta_var[0], "d17-x01-y01", _bins_dabseta, _bins_tt_mass);
 
       book(_h_dabsrapidity_var[1], "d23-x01-y01", _bins_dabsrapidity, _bins_tt_pT);
       book(_h_dabseta_var[1], "d29-x01-y01", _bins_dabseta, _bins_tt_pT);
 
       book(_h_dabsrapidity_var[2], "d35-x01-y01", _bins_dabsrapidity, _bins_tt_absrapidity);
       book(_h_dabseta_var[2], "d41-x01-y01", _bins_dabseta, _bins_tt_absrapidity);
 
       // Profile histos for asymmetries
       book(_h_dabsrapidity_profile[0], "d08-x01-y01", _bins_tt_mass);
       book(_h_dabseta_profile[0], "d14-x01-y01", _bins_tt_mass);
 
       book(_h_dabsrapidity_profile[1], "d20-x01-y01", _bins_tt_pT);
       book(_h_dabseta_profile[1], "d26-x01-y01", _bins_tt_pT);
 
       book(_h_dabsrapidity_profile[2], "d32-x01-y01", _bins_tt_absrapidity);
       book(_h_dabseta_profile[2], "d38-x01-y01", _bins_tt_absrapidity);
 
     }
 
 
     /// Perform the per-event analysis
     void analyze(const Event& event) {
 
       const double weight = 1.0;
 
       // Use particle-level leptons for the first histogram
       const DressedLeptons& dressed_electrons = applyProjection<DressedLeptons>(event, "DressedElectrons");
       const DressedLeptons& dressed_muons = applyProjection<DressedLeptons>(event, "DressedMuons");
 
       const vector<DressedLepton> dressedels = dressed_electrons.dressedLeptons();
       const vector<DressedLepton> dressedmus = dressed_muons.dressedLeptons();
 
       const size_t ndressedel = dressedels.size();
       const size_t ndressedmu = dressedmus.size();
 
       // For the particle-level histogram, require exactly one electron and exactly one muon, to select the ttbar->emu channel.
       // Note this means ttbar->emu events with additional PromptFinalState dilepton pairs from the shower are vetoed - for PYTHIA8,
       // this affects ~0.5% of events, so the effect is well below the level of sensitivity of the measured distribution.
       if ( ndressedel == 1 && ndressedmu == 1 ) {
 
         const int electrontouse = 0, muontouse = 0;
 
         // Opposite-charge leptons only
         if ( sameSign(dressedels[electrontouse], dressedmus[muontouse]) ) {
           MSG_INFO("Error, e and mu have same charge, skipping event");
         } else {
           // Get the four-momenta of the positively- and negatively-charged leptons
           FourMomentum lepPlus = dressedels[electrontouse].charge() > 0 ? dressedels[electrontouse] : dressedmus[muontouse];
           FourMomentum lepMinus = dressedels[electrontouse].charge() > 0 ? dressedmus[muontouse] : dressedels[electrontouse];
 
           // Now calculate the variable
           double dabseta_temp = lepPlus.abseta() - lepMinus.abseta();
 
           fillWithUFOF( _h_dabsetadressedleptons, dabseta_temp, weight );
         }
 
       }
 
 
       // The remaining variables use parton-level information.
 
       // Get the leptonically decaying tops
       const Particles leptonicpartontops = apply<ParticleFinder>(event, "LeptonicPartonTops").particlesByPt();
       Particles chargedleptons;
 
       unsigned int ntrueleptonictops = 0;
       bool oppositesign = false;
 
       if ( leptonicpartontops.size() == 2 ) {
         for (size_t k = 0; k < leptonicpartontops.size(); ++k) {
 
           // Get the lepton
           const Particle lepTop = leptonicpartontops[k];
           const auto isPromptChargedLepton = [](const Particle& p){return (isChargedLepton(p) && isPrompt(p, false, false));};
           Particles lepton_candidates = lepTop.allDescendants(firstParticleWith(isPromptChargedLepton), false);
           if ( lepton_candidates.size() < 1 ) MSG_WARNING("error, PartonicTops::DecayMode::E_MU top quark had no daughter lepton candidate, skipping event.");
 
           // In some cases there is no lepton from the W decay but only leptons from the decay of a radiated gamma.
           // These hadronic PartonicTops are currently being mistakenly selected by PartonicTops::DecayMode::E_MU (as of April 2017), and need to be rejected.
           // PartonicTops::DecayMode::E_MU is being fixed in Rivet, and when it is the veto below should do nothing.
           /// @todo Should no longer be necessary -- remove
           bool istrueleptonictop = false;
           for (size_t i = 0; i < lepton_candidates.size(); ++i) {
             const Particle& lepton_candidate = lepton_candidates[i];
             if ( lepton_candidate.hasParent(PID::PHOTON) ) {
               MSG_DEBUG("Found gamma parent, top: " << k+1 << " of " << leptonicpartontops.size() << " , lepton: " << i+1 << " of " << lepton_candidates.size());
               continue;
             }
             if ( !istrueleptonictop && sameSign(lepTop,lepton_candidate) ) {
               chargedleptons.push_back(lepton_candidate);
               istrueleptonictop = true;
             }
             else MSG_WARNING("Error, found extra prompt charged lepton from top decay (and without gamma parent), ignoring it.");
           }
           if ( istrueleptonictop ) ++ntrueleptonictops;
         }
       }
 
       if ( ntrueleptonictops == 2 ) {
         oppositesign = !( sameSign(chargedleptons[0],chargedleptons[1]) );
         if ( !oppositesign ) MSG_WARNING("error, same charge tops, skipping event.");
       }
 
       if ( ntrueleptonictops == 2 && oppositesign ) {
 
         // Get the four-momenta of the positively- and negatively-charged leptons
         const FourMomentum lepPlus = chargedleptons[0].charge() > 0 ? chargedleptons[0] : chargedleptons[1];
         const FourMomentum lepMinus = chargedleptons[0].charge() > 0 ? chargedleptons[1] : chargedleptons[0];
 
         const double dabseta_temp = lepPlus.abseta() - lepMinus.abseta();
 
         // Get the four-momenta of the positively- and negatively-charged tops
         const FourMomentum topPlus_p4 = leptonicpartontops[0].pdgId() > 0 ? leptonicpartontops[0] : leptonicpartontops[1];
         const FourMomentum topMinus_p4 = leptonicpartontops[0].pdgId() > 0 ? leptonicpartontops[1] : leptonicpartontops[0];
 
         const FourMomentum ttbar_p4 = topPlus_p4 + topMinus_p4;
 
         const double tt_mass_temp = ttbar_p4.mass();
         const double tt_absrapidity_temp = ttbar_p4.absrapidity();
         const double tt_pT_temp = ttbar_p4.pT();
         const double dabsrapidity_temp = topPlus_p4.absrapidity() - topMinus_p4.absrapidity();
 
         // Fill parton-level histos
         fillWithUFOF( _h_dabseta, dabseta_temp, weight );
         fillWithUFOF( _h_dabsrapidity, dabsrapidity_temp, weight );
 
         // Now fill the same variables in the 2D and profile histos vs ttbar invariant mass, pT, and absolute rapidity
         for (int i_var = 0; i_var < 3; ++i_var) {
           double var;
           if ( i_var == 0 ) {
             var = tt_mass_temp;
           } else if ( i_var == 1 ) {
             var = tt_pT_temp;
           } else {
             var = tt_absrapidity_temp;
           }
 
           fillWithUFOF( _h_dabsrapidity_var[i_var], dabsrapidity_temp, var, weight );
           fillWithUFOF( _h_dabseta_var[i_var], dabseta_temp, var, weight );
 
           fillWithUFOF( _h_dabsrapidity_profile[i_var], dabsrapidity_temp, var, weight, (_h_dabsrapidity->xMax() + _h_dabsrapidity->xMin())/2. );
           fillWithUFOF( _h_dabseta_profile[i_var], dabseta_temp, var, weight, (_h_dabseta->xMax() + _h_dabseta->xMin())/2. );
         }
 
       }
 
     }
 
 
     /// Normalise histograms to unit area
     void finalize() {
 
       normalize(_h_dabsetadressedleptons);
 
       normalize(_h_dabseta);
       normalize(_h_dabsrapidity);
 
       for (int i_var = 0; i_var < 3; ++i_var) {
         normalize(_h_dabsrapidity_var[i_var]);
         normalize(_h_dabseta_var[i_var]);
       }
 
     }
 
 
   private:
 
     Histo1DPtr _h_dabsetadressedleptons, _h_dabseta, _h_dabsrapidity;
     Histo2DPtr _h_dabseta_var[3], _h_dabsrapidity_var[3];
     Profile1DPtr _h_dabseta_profile[3], _h_dabsrapidity_profile[3];
 
     const vector<double> _bins_tt_mass = {300., 430., 530., 1200.};
     const vector<double> _bins_tt_pT = {0., 41., 92., 300.};
     const vector<double> _bins_tt_absrapidity = {0., 0.34, 0.75, 1.5};
     const vector<double> _bins_dabseta = { -2., -68./60., -48./60., -32./60., -20./60., -8./60., 0., 8./60., 20./60., 32./60., 48./60., 68./60., 2.};
     const vector<double> _bins_dabsrapidity = {-2., -44./60., -20./60., 0., 20./60., 44./60., 2.};
 
     void fillWithUFOF(Histo1DPtr h, double x, double w) {
       h->fill(std::max(std::min(x, h->xMax()-1e-9),h->xMin()+1e-9), w);
     }
 
     void fillWithUFOF(Histo2DPtr h, double x, double y, double w) {
       h->fill(std::max(std::min(x, h->xMax()-1e-9),h->xMin()+1e-9), std::max(std::min(y, h->yMax()-1e-9),h->yMin()+1e-9), w);
     }
 
     void fillWithUFOF(Profile1DPtr h, double x, double y, double w, double c) {
       h->fill(std::max(std::min(y, h->xMax()-1e-9),h->xMin()+1e-9), float(x > c) - float(x < c), w);
     }
 
 
   };
 
 
   // The hook for the plugin system
   DECLARE_RIVET_PLUGIN(CMS_2016_I1430892);
 
 
 }
diff --git a/analyses/pluginCMS/CMS_2016_I1473674.cc b/analyses/pluginCMS/CMS_2016_I1473674.cc
--- a/analyses/pluginCMS/CMS_2016_I1473674.cc
+++ b/analyses/pluginCMS/CMS_2016_I1473674.cc
@@ -1,124 +1,124 @@
 #include "Rivet/Analysis.hh"
 #include "Rivet/Projections/FinalState.hh"
 #include "Rivet/Projections/FastJets.hh"
 #include "Rivet/Projections/PartonicTops.hh"
 #include "Rivet/Projections/DressedLeptons.hh"
 #include "Rivet/Projections/IdentifiedFinalState.hh"
 #include "Rivet/Projections/PromptFinalState.hh"
 #include "Rivet/Projections/VetoedFinalState.hh"
 #include "Rivet/Projections/InvMassFinalState.hh"
 #include "Rivet/Projections/MissingMomentum.hh"
 
 namespace Rivet {
 
 
   class CMS_2016_I1473674 : public Analysis {
   public:
 
     // Minimal constructor
     DEFAULT_RIVET_ANALYSIS_CTOR(CMS_2016_I1473674);
 
 
     // Set up projections and book histograms
     void init() {
 
       // Complete final state
       FinalState fs;
 
       // Parton level top quarks
       declare(PartonicTops(PartonicTops::DecayMode::E_MU, false), "LeptonicPartonTops");
       declare(PartonicTops(PartonicTops::DecayMode::HADRONIC),    "HadronicPartonTops");
 
       // Projections for dressed electrons and muons
       IdentifiedFinalState photons(fs);
       photons.acceptIdPair(PID::PHOTON);
       //
       IdentifiedFinalState el_id(fs);
       el_id.acceptIdPair(PID::ELECTRON);
       PromptFinalState electrons(el_id);
       declare(electrons, "Electrons");
       DressedLeptons dressed_electrons(photons, electrons, 0.1);
       declare(dressed_electrons, "DressedElectrons");
       //
       IdentifiedFinalState mu_id(fs);
       mu_id.acceptIdPair(PID::MUON);
       PromptFinalState muons(mu_id);
       declare(muons, "Muons");
       DressedLeptons dressed_muons(photons, muons, 0.1);
       declare(dressed_muons, "DressedMuons");
 
       // Projection for jets
-      VetoedFinalState fs_jets(FinalState(-MAXDOUBLE, MAXDOUBLE, 0*GeV));
+      VetoedFinalState fs_jets(FinalState(-DBL_MAX, DBL_MAX, 0*GeV));
       fs_jets.addVetoOnThisFinalState(dressed_muons);
       declare(FastJets(fs_jets, FastJets::ANTIKT, 0.5), "Jets");
 
       // Projections for MET
       declare(MissingMomentum(), "MET");
 
 
       // Booking of histograms
       book(_hist_met ,5, 1, 1);
       book(_hist_ht  ,6, 1, 1);
       book(_hist_st  ,7, 1, 1);
       book(_hist_wpt ,8, 1, 1);
     }
 
 
     /// Per-event analysis
     void analyze(const Event& event) {
       const double weight = 1.0;
 
       // Select ttbar -> lepton+jets at parton level, removing tau decays
       const Particles leptonicpartontops = apply<ParticleFinder>(event, "LeptonicPartonTops").particlesByPt();
       if (leptonicpartontops.size() != 1) vetoEvent;
       const Particles hadronicpartontops = apply<ParticleFinder>(event, "HadronicPartonTops").particlesByPt();
       if (hadronicpartontops.size() != 1) vetoEvent;
 
       // Select ttbar -> lepton+jets at particle level
       const DressedLeptons& dressed_electrons = applyProjection<DressedLeptons>(event, "DressedElectrons");
       const DressedLeptons& dressed_muons = applyProjection<DressedLeptons>(event, "DressedMuons");
       if (dressed_electrons.dressedLeptons().size() + dressed_muons.dressedLeptons().size() != 1) vetoEvent;
       const FourMomentum lepton = (dressed_electrons.dressedLeptons().empty() ? dressed_muons : dressed_electrons).dressedLeptons()[0];
 
       // MET
       const MissingMomentum& met = applyProjection<MissingMomentum>(event, "MET");
       _hist_met->fill(met.visibleMomentum().pT()/GeV, weight);
 
       // HT and ST
       const FastJets& jetpro = applyProjection<FastJets>(event, "Jets");
       const Jets jets = jetpro.jetsByPt(20*GeV);
 
       double ht = 0.0;
       for (const Jet& j : jets) {
         if (deltaR(j.momentum(), lepton) > 0.3) {
           ht += j.pT();
         }
       }
 
       double st = ht + lepton.pT() + met.visibleMomentum().pT();
       _hist_ht->fill(ht/GeV, weight);
       _hist_st->fill(st/GeV, weight);
 
       // WPT
       const FourMomentum w = lepton - met.visibleMomentum();
       _hist_wpt->fill(w.pT()/GeV, weight);
     }
 
 
     /// Normalize histograms
     void finalize() {
       normalize(_hist_met);
       normalize(_hist_ht);
       normalize(_hist_st);
       normalize(_hist_wpt);
     }
 
   private:
     Histo1DPtr _hist_met, _hist_ht, _hist_st, _hist_wpt;
 
   };
 
 
   DECLARE_RIVET_PLUGIN(CMS_2016_I1473674);
 
 }
diff --git a/analyses/pluginD0/D0_2004_S5992206.cc b/analyses/pluginD0/D0_2004_S5992206.cc
--- a/analyses/pluginD0/D0_2004_S5992206.cc
+++ b/analyses/pluginD0/D0_2004_S5992206.cc
@@ -1,137 +1,137 @@
 // -*- C++ -*-
 #include "Rivet/Analysis.hh"
 #include "Rivet/Projections/FastJets.hh"
 #include "Rivet/Projections/VetoedFinalState.hh"
 #include "Rivet/Projections/VisibleFinalState.hh"
 #include "Rivet/Projections/MissingMomentum.hh"
 
 namespace Rivet {
 
 
   /* @brief D0 Run II angular correlations in di-jet events
    * @author Lars Sonnenschein
    *
    * Measurement of angular correlations in di-jet events.
    *
    * @par Run conditions
    *
    * @arg \f$ \sqrt{s} = \f$ 1960 GeV
    * @arg Run with generic QCD events.
    * @arg Several \f$ p_\perp^\text{min} \f$ cutoffs are probably required to fill the histograms:
    *   @arg \f$ p_\perp^\text{min} = \f$ 50, 75, 100, 150 GeV for the four pT ranges respecively
    *
    */
   class D0_2004_S5992206 : public Analysis {
 
   public:
 
     /// @name Constructors etc.
     //@{
 
     /// Constructor.
     D0_2004_S5992206()
       : Analysis("D0_2004_S5992206")
     {  }
 
     //@}
 
 
     /// @name Analysis methods
     //@{
 
     void init() {
       // Final state for jets, mET etc.
       const FinalState fs(-3.0, 3.0);
       declare(fs, "FS");
       // Veto neutrinos, and muons with pT above 1.0 GeV
       VetoedFinalState vfs(fs);
       vfs.vetoNeutrinos();
-      vfs.addVetoPairDetail(PID::MUON, 1.0*GeV, MAXDOUBLE);
+      vfs.addVetoPairDetail(PID::MUON, 1.0*GeV, DBL_MAX);
       declare(vfs, "VFS");
       declare(FastJets(vfs, FastJets::D0ILCONE, 0.7), "Jets");
       declare(MissingMomentum(vfs), "CalMET");
 
       // Book histograms
       book(_histJetAzimuth_pTmax75_100  ,1, 2, 1);
       book(_histJetAzimuth_pTmax100_130 ,2, 2, 1);
       book(_histJetAzimuth_pTmax130_180 ,3, 2, 1);
       book(_histJetAzimuth_pTmax180_    ,4, 2, 1);
     }
 
 
     /// Do the analysis
     void analyze(const Event& event) {
 
       // Analyse and print some info
       const JetAlg& jetpro = apply<JetAlg>(event, "Jets");
       MSG_DEBUG("Jet multiplicity before any pT cut = " << jetpro.size());
 
       const Jets jets  = jetpro.jetsByPt(40.0*GeV);
       if (jets.size() >= 2) {
         MSG_DEBUG("Jet multiplicity after pT > 40 GeV cut = " << jets.size());
       } else {
         vetoEvent;
       }
       const double rap1 = jets[0].rapidity();
       const double rap2 = jets[1].rapidity();
       if (fabs(rap1) > 0.5 || fabs(rap2) > 0.5) {
         vetoEvent;
       }
       MSG_DEBUG("Jet eta and pT requirements fulfilled");
       const double pT1 = jets[0].pT();
 
       const MissingMomentum& caloMissEt = apply<MissingMomentum>(event, "CalMET");
       MSG_DEBUG("Missing vector Et = " << caloMissEt.vectorEt()/GeV << " GeV");
       if (caloMissEt.vectorEt().mod() > 0.7*pT1) {
         MSG_DEBUG("Vetoing event with too much missing ET: "
                   << caloMissEt.vectorEt()/GeV << " GeV > "
                   << 0.7*pT1/GeV << " GeV");
         vetoEvent;
       }
 
       if (pT1/GeV >= 75.0) {
         const double dphi = deltaPhi(jets[0].phi(), jets[1].phi());
         if (inRange(pT1/GeV, 75.0, 100.0)) {
           _histJetAzimuth_pTmax75_100->fill(dphi);
         } else if (inRange(pT1/GeV, 100.0, 130.0)) {
           _histJetAzimuth_pTmax100_130->fill(dphi);
         } else if (inRange(pT1/GeV, 130.0, 180.0)) {
           _histJetAzimuth_pTmax130_180->fill(dphi);
         } else if (pT1/GeV > 180.0) {
           _histJetAzimuth_pTmax180_->fill(dphi);
         }
       }
 
     }
 
 
     // Finalize
     void finalize() {
       // Normalize histograms to unit area
       normalize(_histJetAzimuth_pTmax75_100);
       normalize(_histJetAzimuth_pTmax100_130);
       normalize(_histJetAzimuth_pTmax130_180);
       normalize(_histJetAzimuth_pTmax180_);
     }
 
     //@}
 
 
   private:
 
     /// @name Histograms
     //@{
     Histo1DPtr _histJetAzimuth_pTmax75_100;
     Histo1DPtr _histJetAzimuth_pTmax100_130;
     Histo1DPtr _histJetAzimuth_pTmax130_180;
     Histo1DPtr _histJetAzimuth_pTmax180_;
     //@}
 
   };
 
 
 
   // The hook for the plugin system
   DECLARE_RIVET_PLUGIN(D0_2004_S5992206);
 
 }
diff --git a/analyses/pluginPetra/TASSO_1990_S2148048.cc b/analyses/pluginPetra/TASSO_1990_S2148048.cc
--- a/analyses/pluginPetra/TASSO_1990_S2148048.cc
+++ b/analyses/pluginPetra/TASSO_1990_S2148048.cc
@@ -1,136 +1,136 @@
 // -*- C++ -*-
 #include "Rivet/Analysis.hh"
 #include "Rivet/Projections/Thrust.hh"
 #include "Rivet/Projections/Sphericity.hh"
 #include "Rivet/Projections/ChargedFinalState.hh"
 
 namespace Rivet {
 
 
   class TASSO_1990_S2148048 : public Analysis {
   public:
 
     /// Constructor
     DEFAULT_RIVET_ANALYSIS_CTOR(TASSO_1990_S2148048);
 
 
     /// @name Analysis methods
     //@{
 
     /// Book histograms and initialise projections before the run
     void init() {
-      const ChargedFinalState cfs(-MAXDOUBLE, MAXDOUBLE, 0.1/GeV);
+      const ChargedFinalState cfs(-DBL_MAX, DBL_MAX, 0.1/GeV);
       declare(cfs, "CFS");
 
       // Thrust and sphericity
       declare(Thrust(cfs), "Thrust");
       declare(Sphericity(cfs), "Sphericity");
 
       // Histos
       int offset = 0;
       switch (int(sqrtS()/GeV)) {
         case 14:
           offset = 0;
           break;
         case 22:
           offset = 1;
           break;
         case 35:
           offset = 2;
           break;
         case 44:
           offset = 3;
           break;
       }
       //book(_h_xp         , 2, 1, 1+offset);
       book(_h_sphericity , 6, 1, 1+offset);
       book(_h_aplanarity , 7, 1, 1+offset);
       book(_h_thrust     , 8, 1, 1+offset);
     }
 
 
     /// Perform the per-event analysis
     void analyze(const Event& event) {
       const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
 
       //// Get beams and average beam momentum
       //const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
       //const double meanBeamMom = ( beams.first.p3().mod() +
                                    //beams.second.p3().mod() ) / 2.0;
 
       // TASSO hadronic event selection TODO: move this into a trigger definition
       // See page 2 in publication
       // Condition 1)  --- require at least 5 (4) 'good' tracks
       int nch = cfs.particles().size();
       if ( (int(sqrtS()/GeV) > 27 && nch < 5) || (int(sqrtS()/GeV) <= 27 && nch < 4 ) ) {
         MSG_DEBUG("Failed # good tracks cut: " << nch);
         vetoEvent;
       }
       // Condition 2) ---
       // Condition 5) --- scalar momentum (not pT!!!) sum >= 0.265*s
       double momsum = 0.0;
       for (const Particle& p : cfs.particles()) {
         const double mom = p.p3().mod();
         momsum += mom;
       }
       if (momsum <=0.265 * sqrtS()/GeV) {
         MSG_DEBUG("Failed pTsum cut: " << momsum << " < " << 0.265 * sqrtS()/GeV);
         vetoEvent;
       }
 
       // Raise counter for events that pass trigger conditions
       //_sumWPassed += 1.0;
 
       const Thrust& thrust = apply<Thrust>(event, "Thrust");
       //const Vector3 & thrustAxis = thrust.thrustAxis ();
       //double theta = thrustAxis.theta();
       //if ( fabs(cos(theta)) >= 0.8 ) {
         //MSG_DEBUG("Failed thrust angle cut: " << fabs(cos(theta)));
         //vetoEvent;
       //}
 
       const Sphericity& sphericity = apply<Sphericity>(event, "Sphericity");
 
       //// Fill histograms in order of appearance in paper
       //for (const Particle& p : cfs.particles()) {
         //// Get momentum and energy of each particle.
         //const Vector3 mom3 = p.p3();
         //// Scaled momenta.
         //const double mom = mom3.mod();
         //const double scaledMom = mom/meanBeamMom;
         //_h_xp->fill(scaledMom);
       //}
       //
       _h_sphericity->fill(sphericity.sphericity());
       _h_aplanarity->fill(sphericity.aplanarity());
       _h_thrust->fill(thrust.thrust());
     }
 
 
     /// Normalise histograms etc., after the run
     void finalize() {
       //scale(_h_xp, _sumWPassed/(crossSection()*sumOfWeights()));
       normalize(_h_sphericity);
       normalize(_h_aplanarity);
       normalize(_h_thrust    );
     }
 
     //@}
 
 
   private:
 
     /// @name Histograms
     //@{
     Histo1DPtr _h_xp, _h_sphericity, _h_aplanarity, _h_thrust;
     //@}
 
 
   };
 
 
 
   // The hook for the plugin system
   DECLARE_RIVET_PLUGIN(TASSO_1990_S2148048);
 
 }
diff --git a/include/Rivet/Math/MathHeader.hh b/include/Rivet/Math/MathHeader.hh
--- a/include/Rivet/Math/MathHeader.hh
+++ b/include/Rivet/Math/MathHeader.hh
@@ -1,35 +1,31 @@
 #ifndef RIVET_Math_MathHeader
 #define RIVET_Math_MathHeader
 
 #include "Rivet/Tools/Exceptions.hh"
 #include "Rivet/Tools/Utils.hh"
 
 namespace Rivet {
 
 
   /// Pre-defined numeric type limits
-  /// @deprecated Prefer the standard DBL/INT_MAX
-  static const double MAXDOUBLE = DBL_MAX; // was std::numeric_limits<double>::max(); -- warns in GCC5
-  static const double MAXINT = INT_MAX; // was std::numeric_limits<int>::max(); -- warns in GCC5
-
   /// A pre-defined value of \f$ \pi \f$.
   static const double PI = M_PI;
 
   /// A pre-defined value of \f$ 2\pi \f$.
   static const double TWOPI = 2*M_PI;
 
   /// A pre-defined value of \f$ \pi/2 \f$.
   static const double HALFPI = M_PI_2;
 
   /// Enum for signs of numbers.
   enum Sign { MINUS = -1, ZERO = 0, PLUS = 1 };
 
   /// Enum for rapidity variable to be used in calculating \f$ R \f$, applying rapidity cuts, etc.
   enum RapScheme { PSEUDORAPIDITY = 0, ETARAP = 0, RAPIDITY = 1, YRAP = 1 };
 
   /// Enum for range of \f$ \phi \f$ to be mapped into
   enum PhiMapping { MINUSPI_PLUSPI, ZERO_2PI, ZERO_PI };
 
 }
 
 #endif
diff --git a/include/Rivet/Math/MathUtils.hh b/include/Rivet/Math/MathUtils.hh
--- a/include/Rivet/Math/MathUtils.hh
+++ b/include/Rivet/Math/MathUtils.hh
@@ -1,626 +1,626 @@
 // -*- C++ -*-
 #ifndef RIVET_MathUtils_HH
 #define RIVET_MathUtils_HH
 
 #include "Rivet/Math/MathHeader.hh"
 #include <type_traits>
 #include <cassert>
 
 namespace Rivet {
 
 
   /// @name Comparison functions for safe (floating point) equality tests
   //@{
 
   /// @brief Compare a number to zero
   ///
   /// This version for floating point types has a degree of fuzziness expressed
   /// by the absolute @a tolerance parameter, for floating point safety.
   template <typename NUM>
   inline typename std::enable_if<std::is_floating_point<NUM>::value, bool>::type
   isZero(NUM val, double tolerance=1e-8) {
     return fabs(val) < tolerance;
   }
 
   /// @brief Compare a number to zero
   ///
   /// SFINAE template specialisation for integers, since there is no FP
   /// precision issue.
   template <typename NUM>
   inline typename std::enable_if<std::is_integral<NUM>::value, bool>::type
   isZero(NUM val, double=1e-5) { //< NB. unused tolerance parameter for ints, still needs a default value!
     return val == 0;
   }
 
 
   /// @brief Compare two numbers for equality with a degree of fuzziness
   ///
   /// This version for floating point types (if any argument is FP) has a degree
   /// of fuzziness expressed by the fractional @a tolerance parameter, for
   /// floating point safety.
   template <typename N1, typename N2>
   inline typename std::enable_if<
     std::is_arithmetic<N1>::value && std::is_arithmetic<N2>::value &&
    (std::is_floating_point<N1>::value || std::is_floating_point<N2>::value), bool>::type
   fuzzyEquals(N1 a, N2 b, double tolerance=1e-5) {
     const double absavg = (std::abs(a) + std::abs(b))/2.0;
     const double absdiff = std::abs(a - b);
     const bool rtn = (isZero(a) && isZero(b)) || absdiff < tolerance*absavg;
     return rtn;
   }
 
   /// @brief Compare two numbers for equality with a degree of fuzziness
   ///
   /// Simpler SFINAE template specialisation for integers, since there is no FP
   /// precision issue.
   template <typename N1, typename N2>
   inline typename std::enable_if<
     std::is_integral<N1>::value && std::is_integral<N2>::value, bool>::type
     fuzzyEquals(N1 a, N2 b, double) { //< NB. unused tolerance parameter for ints, still needs a default value!
     return a == b;
   }
 
 
   /// @brief Compare two numbers for >= with a degree of fuzziness
   ///
   /// The @a tolerance parameter on the equality test is as for @c fuzzyEquals.
   template <typename N1, typename N2>
   inline typename std::enable_if<
     std::is_arithmetic<N1>::value && std::is_arithmetic<N2>::value, bool>::type
   fuzzyGtrEquals(N1 a, N2 b, double tolerance=1e-5) {
     return a > b || fuzzyEquals(a, b, tolerance);
   }
 
 
   /// @brief Compare two floating point numbers for <= with a degree of fuzziness
   ///
   /// The @a tolerance parameter on the equality test is as for @c fuzzyEquals.
   template <typename N1, typename N2>
   inline typename std::enable_if<
     std::is_arithmetic<N1>::value && std::is_arithmetic<N2>::value, bool>::type
   fuzzyLessEquals(N1 a, N2 b, double tolerance=1e-5) {
     return a < b || fuzzyEquals(a, b, tolerance);
   }
 
   //@}
 
 
   /// @name Ranges and intervals
   //@{
 
   /// Represents whether an interval is open (non-inclusive) or closed (inclusive).
   ///
   /// For example, the interval \f$ [0, \pi) \f$ is closed (an inclusive
   /// boundary) at 0, and open (a non-inclusive boundary) at \f$ \pi \f$.
   enum RangeBoundary { OPEN=0, SOFT=0, CLOSED=1, HARD=1 };
 
   /// @brief Determine if @a value is in the range @a low to @a high, for floating point numbers
   ///
   /// Interval boundary types are defined by @a lowbound and @a highbound.
   template <typename N1, typename N2, typename N3>
   inline typename std::enable_if<
     std::is_arithmetic<N1>::value && std::is_arithmetic<N2>::value && std::is_arithmetic<N3>::value, bool>::type
   inRange(N1 value, N2 low, N3 high,
           RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) {
     if (lowbound == OPEN && highbound == OPEN) {
       return (value > low && value < high);
     } else if (lowbound == OPEN && highbound == CLOSED) {
       return (value > low && value <= high);
     } else if (lowbound == CLOSED && highbound == OPEN) {
       return (value >= low && value < high);
     } else { // if (lowbound == CLOSED && highbound == CLOSED) {
       return (value >= low && value <= high);
     }
   }
 
   /// @brief Determine if @a value is in the range @a low to @a high, for floating point numbers
   ///
   /// Interval boundary types are defined by @a lowbound and @a highbound.
   /// Closed intervals are compared fuzzily.
   template <typename N1, typename N2, typename N3>
   inline typename std::enable_if<
     std::is_arithmetic<N1>::value && std::is_arithmetic<N2>::value && std::is_arithmetic<N3>::value, bool>::type
   fuzzyInRange(N1 value, N2 low, N3 high,
                RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) {
     if (lowbound == OPEN && highbound == OPEN) {
       return (value > low && value < high);
     } else if (lowbound == OPEN && highbound == CLOSED) {
       return (value > low && fuzzyLessEquals(value, high));
     } else if (lowbound == CLOSED && highbound == OPEN) {
       return (fuzzyGtrEquals(value, low) && value < high);
     } else { // if (lowbound == CLOSED && highbound == CLOSED) {
       return (fuzzyGtrEquals(value, low) && fuzzyLessEquals(value, high));
     }
   }
 
   /// Alternative version of inRange which accepts a pair for the range arguments.
   template <typename N1, typename N2, typename N3>
   inline typename std::enable_if<
     std::is_arithmetic<N1>::value && std::is_arithmetic<N2>::value && std::is_arithmetic<N3>::value, bool>::type
   inRange(N1 value, pair<N2, N3> lowhigh,
           RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) {
     return inRange(value, lowhigh.first, lowhigh.second, lowbound, highbound);
   }
 
 
   // Alternative forms, with snake_case names and boundary types in names rather than as args -- from MCUtils
 
   /// @brief Boolean function to determine if @a value is within the given range
   ///
   /// @note The interval is closed (inclusive) at the low end, and open (exclusive) at the high end.
   template <typename N1, typename N2, typename N3>
   inline typename std::enable_if<
     std::is_arithmetic<N1>::value && std::is_arithmetic<N2>::value && std::is_arithmetic<N3>::value, bool>::type
   in_range(N1 val, N2 low, N3 high) {
     return inRange(val, low, high, CLOSED, OPEN);
   }
 
   /// @brief Boolean function to determine if @a value is within the given range
   ///
   /// @note The interval is closed at both ends.
   template <typename N1, typename N2, typename N3>
   inline typename std::enable_if<
     std::is_arithmetic<N1>::value && std::is_arithmetic<N2>::value && std::is_arithmetic<N3>::value, bool>::type
   in_closed_range(N1 val, N2 low, N3 high) {
     return inRange(val, low, high, CLOSED, CLOSED);
   }
 
   /// @brief Boolean function to determine if @a value is within the given range
   ///
   /// @note The interval is open at both ends.
   template <typename N1, typename N2, typename N3>
   inline typename std::enable_if<
     std::is_arithmetic<N1>::value && std::is_arithmetic<N2>::value && std::is_arithmetic<N3>::value, bool>::type
   in_open_range(N1 val, N2 low, N3 high) {
     return inRange(val, low, high, OPEN, OPEN);
   }
 
   /// @todo Add pair-based versions of the named range-boundary functions
 
   //@}
 
 
   /// @name Miscellaneous numerical helpers
   //@{
 
   /// Named number-type squaring operation.
   template <typename NUM>
   inline typename std::enable_if<std::is_arithmetic<NUM>::value, NUM>::type
   sqr(NUM a) {
     return a*a;
   }
 
   /// @brief Named number-type addition in quadrature operation.
   ///
   /// @note Result has the sqrt operation applied.
   /// @todo When std::common_type can be used, generalise to multiple numeric types with appropriate return type.
   // template <typename N1, typename N2>
   template <typename NUM>
   inline typename std::enable_if<std::is_arithmetic<NUM>::value, NUM>::type
   //std::common_type<N1, N2>::type
   add_quad(NUM a, NUM b) {
     return sqrt(a*a + b*b);
   }
 
   /// Named number-type addition in quadrature operation.
   ///
   /// @note Result has the sqrt operation applied.
   /// @todo When std::common_type can be used, generalise to multiple numeric types with appropriate return type.
   // template <typename N1, typename N2>
   template <typename NUM>
   inline typename std::enable_if<std::is_arithmetic<NUM>::value, NUM>::type
   //std::common_type<N1, N2, N3>::type
   add_quad(NUM a, NUM b, NUM c) {
     return sqrt(a*a + b*b + c*c);
   }
 
   /// Return a/b, or @a fail if b = 0
   /// @todo When std::common_type can be used, generalise to multiple numeric types with appropriate return type.
   inline double safediv(double num, double den, double fail=0.0) {
     return (!isZero(den)) ? num/den : fail;
   }
 
   /// A more efficient version of pow for raising numbers to integer powers.
   template <typename NUM>
   inline typename std::enable_if<std::is_arithmetic<NUM>::value, NUM>::type
   intpow(NUM val, unsigned int exp) {
     assert(exp >= 0);
     if (exp == 0) return (NUM) 1;
     else if (exp == 1) return val;
     return val * intpow(val, exp-1);
   }
 
   /// Find the sign of a number
   template <typename NUM>
   inline typename std::enable_if<std::is_arithmetic<NUM>::value, int>::type
   sign(NUM val) {
     if (isZero(val)) return ZERO;
     const int valsign = (val > 0) ? PLUS : MINUS;
     return valsign;
   }
 
   //@}
 
 
   /// @name Physics statistical distributions
   //@{
 
   /// @brief CDF for the Breit-Wigner distribution
   inline double cdfBW(double x, double mu, double gamma) {
     // normalize to (0;1) distribution
     const double xn = (x - mu)/gamma;
     return std::atan(xn)/M_PI + 0.5;
   }
 
   /// @brief Inverse CDF for the Breit-Wigner distribution
   inline double invcdfBW(double p, double mu, double gamma) {
     const double xn = std::tan(M_PI*(p-0.5));
     return gamma*xn + mu;
   }
 
   //@}
 
 
   /// @name Binning helper functions
   //@{
 
   /// @brief Make a list of @a nbins + 1 values equally spaced between @a start and @a end inclusive.
   ///
   /// NB. The arg ordering and the meaning of the nbins variable is "histogram-like",
   /// as opposed to the Numpy/Matlab version.
   inline vector<double> linspace(size_t nbins, double start, double end, bool include_end=true) {
     assert(end >= start);
     assert(nbins > 0);
     vector<double> rtn;
     const double interval = (end-start)/static_cast<double>(nbins);
     for (size_t i = 0; i < nbins; ++i) {
       rtn.push_back(start + i*interval);
     }
     assert(rtn.size() == nbins);
     if (include_end) rtn.push_back(end); //< exact end, not result of n * interval
     return rtn;
   }
 
 
   /// @brief Make a list of @a nbins + 1 values exponentially spaced between @a start and @a end inclusive.
   ///
   /// NB. The arg ordering and the meaning of the nbins variable is "histogram-like",
   /// as opposed to the Numpy/Matlab version, and the start and end arguments are expressed
   /// in "normal" space, rather than as the logarithms of the start/end values as in Numpy/Matlab.
   inline vector<double> logspace(size_t nbins, double start, double end, bool include_end=true) {
     assert(end >= start);
     assert(start > 0);
     assert(nbins > 0);
     const double logstart = std::log(start);
     const double logend = std::log(end);
     const vector<double> logvals = linspace(nbins, logstart, logend, false);
     assert(logvals.size() == nbins);
     vector<double> rtn; rtn.reserve(nbins+1);
     rtn.push_back(start); //< exact start, not exp(log(start))
     for (size_t i = 1; i < logvals.size(); ++i) {
       rtn.push_back(std::exp(logvals[i]));
     }
     assert(rtn.size() == nbins);
     if (include_end) rtn.push_back(end); //< exact end, not exp(n * loginterval)
     return rtn;
   }
 
 
   /// @todo geomspace
 
 
   /// @brief Make a list of @a nbins + 1 values spaced for equal area
   /// Breit-Wigner binning between @a start and @a end inclusive. @a
   /// mu and @a gamma are the Breit-Wigner parameters.
   ///
   /// @note The arg ordering and the meaning of the nbins variable is "histogram-like",
   /// as opposed to the Numpy/Matlab version, and the start and end arguments are expressed
   /// in "normal" space.
   inline vector<double> bwspace(size_t nbins, double start, double end, double mu, double gamma) {
     assert(end >= start);
     assert(nbins > 0);
     const double pmin = cdfBW(start, mu, gamma);
     const double pmax = cdfBW(end,   mu, gamma);
     const vector<double> edges = linspace(nbins, pmin, pmax);
     assert(edges.size() == nbins+1);
     vector<double> rtn;
     for (double edge : edges) {
       rtn.push_back(invcdfBW(edge, mu, gamma));
     }
     assert(rtn.size() == nbins+1);
     return rtn;
   }
 
 
   /// Actual helper implementation of binIndex (so generic and specific overloading can work)
   template <typename NUM, typename CONTAINER>
   inline typename std::enable_if<std::is_arithmetic<NUM>::value && std::is_arithmetic<typename CONTAINER::value_type>::value, int>::type
   _binIndex(NUM val, const CONTAINER& binedges, bool allow_overflow=false) {
     if (val < *begin(binedges)) return -1; ///< Below/out of histo range
     // CONTAINER::iterator_type itend =
     if (val >= *(end(binedges)-1)) return allow_overflow ? int(binedges.size())-1 : -1; ///< Above/out of histo range
     auto it = std::upper_bound(begin(binedges), end(binedges), val);
     return std::distance(begin(binedges), --it);
   }
 
   /// @brief Return the bin index of the given value, @a val, given a vector of bin edges
   ///
   /// An underflow always returns -1. If allow_overflow is false (default) an overflow
   /// also returns -1, otherwise it returns the Nedge-1, the index of an inclusive bin
   /// starting at the last edge.
   ///
   /// @note The @a binedges vector must be sorted
   /// @todo Use std::common_type<NUM1, NUM2>::type x = val; ?
   template <typename NUM1, typename NUM2>
   inline typename std::enable_if<std::is_arithmetic<NUM1>::value && std::is_arithmetic<NUM2>::value, int>::type
   binIndex(NUM1 val, std::initializer_list<NUM2> binedges, bool allow_overflow=false) {
     return _binIndex(val, binedges, allow_overflow);
   }
 
   /// @brief Return the bin index of the given value, @a val, given a vector of bin edges
   ///
   /// An underflow always returns -1. If allow_overflow is false (default) an overflow
   /// also returns -1, otherwise it returns the Nedge-1, the index of an inclusive bin
   /// starting at the last edge.
   ///
   /// @note The @a binedges vector must be sorted
   /// @todo Use std::common_type<NUM1, NUM2>::type x = val; ?
   template <typename NUM, typename CONTAINER>
   inline typename std::enable_if<std::is_arithmetic<NUM>::value && std::is_arithmetic<typename CONTAINER::value_type>::value, int>::type
   binIndex(NUM val, const CONTAINER& binedges, bool allow_overflow=false) {
     return _binIndex(val, binedges, allow_overflow);
   }
 
   //@}
 
 
   /// @name Discrete statistics functions
   //@{
 
   /// Calculate the median of a sample
   /// @todo Support multiple container types via SFINAE
   template <typename NUM>
   inline typename std::enable_if<std::is_arithmetic<NUM>::value, NUM>::type
   median(const vector<NUM>& sample) {
     if (sample.empty()) throw RangeError("Can't compute median of an empty set");
     vector<NUM> tmp = sample;
     std::sort(tmp.begin(), tmp.end());
     const size_t imid = tmp.size()/2; // len1->idx0, len2->idx1, len3->idx1, len4->idx2, ...
     if (sample.size() % 2 == 0) return (tmp.at(imid-1) + tmp.at(imid)) / 2.0;
     else return tmp.at(imid);
   }
 
 
   /// Calculate the mean of a sample
   /// @todo Support multiple container types via SFINAE
   template <typename NUM>
   inline typename std::enable_if<std::is_arithmetic<NUM>::value, double>::type
   mean(const vector<NUM>& sample) {
     if (sample.empty()) throw RangeError("Can't compute mean of an empty set");
     double mean = 0.0;
     for (size_t i = 0; i < sample.size(); ++i) {
       mean += sample[i];
     }
     return mean/sample.size();
   }
 
   // Calculate the error on the mean, assuming Poissonian errors
   /// @todo Support multiple container types via SFINAE
   template <typename NUM>
   inline typename std::enable_if<std::is_arithmetic<NUM>::value, double>::type
   mean_err(const vector<NUM>& sample) {
     if (sample.empty()) throw RangeError("Can't compute mean_err of an empty set");
     double mean_e = 0.0;
     for (size_t i = 0; i < sample.size(); ++i) {
       mean_e += sqrt(sample[i]);
     }
     return mean_e/sample.size();
   }
 
 
   /// Calculate the covariance (variance) between two samples
   /// @todo Support multiple container types via SFINAE
   template <typename NUM>
   inline typename std::enable_if<std::is_arithmetic<NUM>::value, double>::type
   covariance(const vector<NUM>& sample1, const vector<NUM>& sample2) {
     if (sample1.empty() || sample2.empty()) throw RangeError("Can't compute covariance of an empty set");
     if (sample1.size() != sample2.size()) throw RangeError("Sizes of samples must be equal for covariance calculation");
     const double mean1 = mean(sample1);
     const double mean2 = mean(sample2);
     const size_t N = sample1.size();
     double cov = 0.0;
     for (size_t i = 0; i < N; i++) {
       const double cov_i = (sample1[i] - mean1)*(sample2[i] - mean2);
       cov += cov_i;
     }
     if (N > 1) return cov/(N-1);
     else return 0.0;
   }
 
   /// Calculate the error on the covariance (variance) of two samples, assuming poissonian errors
   /// @todo Support multiple container types via SFINAE
   template <typename NUM>
   inline typename std::enable_if<std::is_arithmetic<NUM>::value, double>::type
   covariance_err(const vector<NUM>& sample1, const vector<NUM>& sample2) {
     if (sample1.empty() || sample2.empty()) throw RangeError("Can't compute covariance_err of an empty set");
     if (sample1.size() != sample2.size()) throw RangeError("Sizes of samples must be equal for covariance_err calculation");
     const double mean1 = mean(sample1);
     const double mean2 = mean(sample2);
     const double mean1_e = mean_err(sample1);
     const double mean2_e = mean_err(sample2);
     const size_t N = sample1.size();
     double cov_e = 0.0;
     for (size_t i = 0; i < N; i++) {
       const double cov_i = (sqrt(sample1[i]) - mean1_e)*(sample2[i] - mean2) +
         (sample1[i] - mean1)*(sqrt(sample2[i]) - mean2_e);
       cov_e += cov_i;
     }
     if (N > 1) return cov_e/(N-1);
     else return 0.0;
   }
 
 
   /// Calculate the correlation strength between two samples
   /// @todo Support multiple container types via SFINAE
   template <typename NUM>
   inline typename std::enable_if<std::is_arithmetic<NUM>::value, double>::type
   correlation(const vector<NUM>& sample1, const vector<NUM>& sample2) {
     const double cov = covariance(sample1, sample2);
     const double var1 = covariance(sample1, sample1);
     const double var2 = covariance(sample2, sample2);
     const double correlation = cov/sqrt(var1*var2);
     const double corr_strength = correlation*sqrt(var2/var1);
     return corr_strength;
   }
 
   /// Calculate the error of the correlation strength between two samples assuming Poissonian errors
   /// @todo Support multiple container types via SFINAE
   template <typename NUM>
   inline typename std::enable_if<std::is_arithmetic<NUM>::value, double>::type
   correlation_err(const vector<NUM>& sample1, const vector<NUM>& sample2) {
     const double cov = covariance(sample1, sample2);
     const double var1 = covariance(sample1, sample1);
     const double var2 = covariance(sample2, sample2);
     const double cov_e = covariance_err(sample1, sample2);
     const double var1_e = covariance_err(sample1, sample1);
     const double var2_e = covariance_err(sample2, sample2);
 
     // Calculate the correlation
     const double correlation = cov/sqrt(var1*var2);
     // Calculate the error on the correlation
     const double correlation_err = cov_e/sqrt(var1*var2) -
       cov/(2*pow(3./2., var1*var2)) * (var1_e * var2 + var1 * var2_e);
 
     // Calculate the error on the correlation strength
     const double corr_strength_err = correlation_err*sqrt(var2/var1) +
       correlation/(2*sqrt(var2/var1)) * (var2_e/var1 - var2*var1_e/pow(2, var2));
 
     return corr_strength_err;
   }
 
   //@}
 
 
   /// @name Angle range mappings
   //@{
 
   /// @brief Reduce any number to the range [-2PI, 2PI]
   ///
   /// Achieved by repeated addition or subtraction of 2PI as required. Used to
   /// normalise angular measures.
   inline double _mapAngleM2PITo2Pi(double angle) {
     double rtn = fmod(angle, TWOPI);
     if (isZero(rtn)) return 0;
     assert(rtn >= -TWOPI && rtn <= TWOPI);
     return rtn;
   }
 
   /// Map an angle into the range (-PI, PI].
   inline double mapAngleMPiToPi(double angle) {
     double rtn = _mapAngleM2PITo2Pi(angle);
     if (isZero(rtn)) return 0;
     if (rtn > PI) rtn -= TWOPI;
     if (rtn <= -PI) rtn += TWOPI;
     assert(rtn > -PI && rtn <= PI);
     return rtn;
   }
 
   /// Map an angle into the range [0, 2PI).
   inline double mapAngle0To2Pi(double angle) {
     double rtn = _mapAngleM2PITo2Pi(angle);
     if (isZero(rtn)) return 0;
     if (rtn < 0) rtn += TWOPI;
     if (rtn == TWOPI) rtn = 0;
     assert(rtn >= 0 && rtn < TWOPI);
     return rtn;
   }
 
   /// Map an angle into the range [0, PI].
   inline double mapAngle0ToPi(double angle) {
     double rtn = fabs(mapAngleMPiToPi(angle));
     if (isZero(rtn)) return 0;
     assert(rtn > 0 && rtn <= PI);
     return rtn;
   }
 
   /// Map an angle into the enum-specified range.
   inline double mapAngle(double angle, PhiMapping mapping) {
     switch (mapping) {
     case MINUSPI_PLUSPI:
       return mapAngleMPiToPi(angle);
     case ZERO_2PI:
       return mapAngle0To2Pi(angle);
     case ZERO_PI:
       return mapAngle0To2Pi(angle);
     default:
       throw Rivet::UserError("The specified phi mapping scheme is not implemented");
     }
   }
 
   //@}
 
 
   /// @name Phase-space measure helpers
   //@{
 
   /// @brief Calculate the difference between two angles in radians
   ///
   /// Returns in the range [0, PI].
   inline double deltaPhi(double phi1, double phi2) {
     return mapAngle0ToPi(phi1 - phi2);
   }
 
   /// Calculate the abs difference between two pseudorapidities
   ///
   /// @note Just a cosmetic name for analysis code clarity.
   inline double deltaEta(double eta1, double eta2) {
     return fabs(eta1 - eta2);
   }
 
   /// Calculate the abs difference between two rapidities
   ///
   /// @note Just a cosmetic name for analysis code clarity.
   inline double deltaRap(double y1, double y2) {
     return fabs(y1 - y2);
   }
 
   /// Calculate the squared distance between two points in 2D rapidity-azimuthal
   /// ("\f$ \eta-\phi \f$") space. The phi values are given in radians.
   inline double deltaR2(double rap1, double phi1, double rap2, double phi2) {
     const double dphi = deltaPhi(phi1, phi2);
     return sqr(rap1-rap2) + sqr(dphi);
   }
 
   /// Calculate the distance between two points in 2D rapidity-azimuthal
   /// ("\f$ \eta-\phi \f$") space. The phi values are given in radians.
   inline double deltaR(double rap1, double phi1, double rap2, double phi2) {
     return sqrt(deltaR2(rap1, phi1, rap2, phi2));
   }
 
   /// Calculate a rapidity value from the supplied energy @a E and longitudinal momentum @a pz.
   inline double rapidity(double E, double pz) {
     if (isZero(E - pz)) {
       throw std::runtime_error("Divergent positive rapidity");
-      return MAXDOUBLE;
+      return DBL_MAX;
     }
     if (isZero(E + pz)) {
       throw std::runtime_error("Divergent negative rapidity");
-      return -MAXDOUBLE;
+      return -DBL_MAX;
     }
     return 0.5*log((E+pz)/(E-pz));
   }
 
   //@}
 
 
   /// Calculate transverse mass of two vectors from provided pT and deltaPhi
   /// @note Several versions taking two vectors are found in Vector4.hh
   inline double mT(double pT1, double pT2, double dphi) {
     return sqrt(2*pT1*pT2 * (1 - cos(dphi)) );
   }
 
 
 }
 
 
 #endif
diff --git a/include/Rivet/Projections/ConstLossyFinalState.hh b/include/Rivet/Projections/ConstLossyFinalState.hh
--- a/include/Rivet/Projections/ConstLossyFinalState.hh
+++ b/include/Rivet/Projections/ConstLossyFinalState.hh
@@ -1,77 +1,77 @@
 // -*- C++ -*-
 #ifndef RIVET_ConstLossyFinalState_HH
 #define RIVET_ConstLossyFinalState_HH
 
 #include "Rivet/Tools/Logging.hh"
 #include "Rivet/Config/RivetCommon.hh"
 #include "Rivet/Particle.hh"
 #include "Rivet/Event.hh"
 #include "Rivet/Projection.hh"
 #include "Rivet/Projections/FinalState.hh"
 #include "Rivet/Projections/LossyFinalState.hh"
 
 namespace Rivet {
 
 
   /// Functor used to implement constant random lossiness.
   class ConstRandomFilter {
   public:
 
     ConstRandomFilter(double lossFraction)
       : _lossFraction(lossFraction)
     {
       assert(_lossFraction >= 0);
     }
 
     // If operator() returns true, particle is deleted ("lost")
     bool operator()(const Particle&) {
       return rand01() < _lossFraction;
     }
 
     CmpState compare(const ConstRandomFilter& other) const {
       return cmp(_lossFraction, other._lossFraction);
     }
 
   private:
 
     double _lossFraction;
 
   };
 
 
 
   /// @brief Randomly lose a constant fraction of particles.
   class ConstLossyFinalState : public LossyFinalState<ConstRandomFilter> {
   public:
 
     /// @name Constructors
     //@{
 
     /// Constructor from a FinalState.
     ConstLossyFinalState(const FinalState& fsp, double lossfraction)
       : LossyFinalState<ConstRandomFilter>(fsp, ConstRandomFilter(lossfraction))
     {
       setName("ConstLossyFinalState");
     }
 
     /// Stand-alone constructor. Initialises the base FinalState projection.
     ConstLossyFinalState(double lossfraction,
-                         double mineta = -MAXDOUBLE,
-                         double maxeta = MAXDOUBLE,
+                         double mineta = -DBL_MAX,
+                         double maxeta = DBL_MAX,
                          double minpt = 0.0)
       : LossyFinalState<ConstRandomFilter>(ConstRandomFilter(lossfraction), mineta, maxeta, minpt)
     {
       setName("ConstLossyFinalState");
     }
 
     /// Clone on the heap.
     DEFAULT_RIVET_PROJ_CLONE(ConstLossyFinalState);
 
     //@}
 
   };
 
 
 }
 
 #endif
diff --git a/include/Rivet/Projections/HadronicFinalState.hh b/include/Rivet/Projections/HadronicFinalState.hh
--- a/include/Rivet/Projections/HadronicFinalState.hh
+++ b/include/Rivet/Projections/HadronicFinalState.hh
@@ -1,51 +1,51 @@
 // -*- C++ -*-
 #ifndef RIVET_HadronicFinalState_HH
 #define RIVET_HadronicFinalState_HH
 
 #include "Rivet/Tools/Logging.hh"
 #include "Rivet/Config/RivetCommon.hh"
 #include "Rivet/Particle.hh"
 #include "Rivet/Event.hh"
 #include "Rivet/Projection.hh"
 #include "Rivet/Projections/FinalState.hh"
 
 namespace Rivet {
 
 
   /// @brief Project only hadronic final state particles.
   class HadronicFinalState : public FinalState {
   public:
 
     /// Constructor: the supplied FinalState projection is assumed to live through the run.
     HadronicFinalState(const FinalState& fsp)
     {
       setName("HadronicFinalState");
       declare(fsp, "FS");
     }
 
-    HadronicFinalState(double mineta = -MAXDOUBLE,
-                       double maxeta = MAXDOUBLE,
+    HadronicFinalState(double mineta = -DBL_MAX,
+                       double maxeta = DBL_MAX,
                        double minpt = 0.0*GeV)
     {
       setName("HadronicFinalState");
       declare(FinalState(mineta, maxeta, minpt), "FS");
     }
 
     /// Clone on the heap.
     DEFAULT_RIVET_PROJ_CLONE(HadronicFinalState);
 
   protected:
 
     /// Apply the projection on the supplied event.
     void project(const Event& e);
 
     /// Compare projections.
     CmpState compare(const Projection& p) const;
 
   };
 
 
 }
 
 
 #endif
diff --git a/include/Rivet/Projections/JetShape.hh b/include/Rivet/Projections/JetShape.hh
--- a/include/Rivet/Projections/JetShape.hh
+++ b/include/Rivet/Projections/JetShape.hh
@@ -1,199 +1,199 @@
 // -*- C++ -*-
 #ifndef RIVET_JetShape_HH
 #define RIVET_JetShape_HH
 
 #include "Rivet/Config/RivetCommon.hh"
 #include "Rivet/Projection.hh"
 #include "Rivet/Projections/JetAlg.hh"
 #include "Rivet/Particle.hh"
 #include "Rivet/Event.hh"
 #include "Rivet/Tools/Utils.hh"
 
 namespace Rivet {
 
 
   ///*
   /// @brief Calculate the jet shape.
   ///
   /// Calculate the differential and integral jet shapes in \f$P_{\perp}\f$ for a given
   /// set of jets. This particular jet shape projection calculates jet shapes relative
   /// to jet centroids, using only the particles associated to each jet, for the hardest
   /// \f$ n \f$ jets.
   ///
   /// The rapidity scheme (\f$ \eta \f$ or \f$ y \f$) has to be specified when
   /// invoking the constructor.
   ///
   /// The differential jet shape around a given jet axis at distance interval
   /// \f$ r \pm \delta{r}/2 \f$ is defined as
   /// \f[
   /// \rho(r) =
   ///   \frac{1}{\delta r} \frac{1}{N_\mathrm{jets}}
   ///   \sum_\mathrm{jets} \frac{P_\perp(r - \delta r/2, r+\delta r/2)}{p_\perp(0, R)}
   /// \f]
   /// with \f$ 0 \le r \le R \f$ and \f$ P_\perp(r_1, r_2) = \sum_{\in [r_1, r_2)} p_\perp \f$.
   ///
   /// The integral jet shape around a given jet axes until distance \f$ r \f$ is defined as
   /// \f[
   /// \Psi(r) =
   ///   \frac{1}{N_\mathrm{jets}}
   ///   \sum_\mathrm{jets} \frac{P_\perp(0, r)}{p_\perp(0, R)}
   /// \f]
   /// with \f$ 0 \le r \le R \f$ and \f$ P_\perp(r_1, r_2) = \sum_{\in [r_1, r_2)} p_\perp \f$.
   ///
   /// The constructor expects also the binning in radius \f$ r \f$ to be supplied.
   ///
   class JetShape : public Projection {
   public:
 
     /// @name Constructors etc.
     //@{
 
     /// Constructor from histo range and number of bins.
     JetShape(const JetAlg& jetalg,
              double rmin, double rmax, size_t nbins,
-             double ptmin=0, double ptmax=MAXDOUBLE,
-             double absrapmin=-MAXDOUBLE, double absrapmax=-MAXDOUBLE,
+             double ptmin=0, double ptmax=DBL_MAX,
+             double absrapmin=-DBL_MAX, double absrapmax=-DBL_MAX,
              RapScheme rapscheme=RAPIDITY);
 
     /// Constructor from vector of bin edges.
     JetShape(const JetAlg& jetalg, vector<double> binedges,
-             double ptmin=0, double ptmax=MAXDOUBLE,
-             double absrapmin=-MAXDOUBLE, double absrapmax=-MAXDOUBLE,
+             double ptmin=0, double ptmax=DBL_MAX,
+             double absrapmin=-DBL_MAX, double absrapmax=-DBL_MAX,
              RapScheme rapscheme=RAPIDITY);
 
     /// Clone on the heap.
     DEFAULT_RIVET_PROJ_CLONE(JetShape);
 
     //@}
 
 
     /// Reset projection between events.
     void clear();
 
 
     /// Do the calculation directly on a supplied collection of Jet objects.
     void calc(const Jets& jets);
 
 
   public:
 
 
     /// Number of equidistant radius bins.
     size_t numBins() const {
       return _binedges.size() - 1;
     }
 
     /// Number of jets which passed cuts.
     size_t numJets() const {
       return _diffjetshapes.size();
     }
 
     /// \f$ r_\text{min} \f$ value.
     double rMin() const {
       return _binedges.front();
     }
 
     /// \f$ r_\text{max} \f$ value.
     double rMax() const {
       return _binedges.back();
     }
 
     /// \f$ p_\perp^\text{min} \f$ value.
     double ptMin() const {
       return _ptcuts.first;
     }
 
     /// \f$ p_\perp^\text{max} \f$ value.
     double ptMax() const {
       return _ptcuts.second;
     }
 
     /// Central \f$ r \f$ value for bin @a rbin.
     double rBinMin(size_t rbin) const {
       assert(inRange(rbin, 0u, numBins()));
       return _binedges[rbin];
     }
 
     /// Central \f$ r \f$ value for bin @a rbin.
     double rBinMax(size_t rbin) const {
       assert(inRange(rbin, 0u, numBins()));
       return _binedges[rbin+1];
     }
 
     /// Central \f$ r \f$ value for bin @a rbin.
     double rBinMid(size_t rbin) const {
       assert(inRange(rbin, 0u, numBins()));
       //cout << _binedges << '\n';
       return (_binedges[rbin] + _binedges[rbin+1])/2.0;
     }
 
     /// Return value of differential jet shape profile histo bin.
     double diffJetShape(size_t ijet, size_t rbin) const {
       assert(inRange(ijet, 0u, numJets()));
       assert(inRange(rbin, 0u, numBins()));
       return _diffjetshapes[ijet][rbin];
     }
 
     /// Return value of integrated jet shape profile histo bin.
     double intJetShape(size_t ijet, size_t rbin) const {
       assert(inRange(ijet, 0u, numJets()));
       assert(inRange(rbin, 0u, numBins()));
       double rtn  = 0;
       for (size_t i = 0; i <= rbin; ++i) {
         rtn += _diffjetshapes[ijet][i];
       }
       return rtn;
     }
 
     /// @todo Provide int and diff jet shapes with some sort of area normalisation?
 
     // /// Return value of \f$ \Psi \f$ (integrated jet shape) at given radius for a \f$ p_T \f$ bin.
     // /// @todo Remove this external indexing thing
     // double psi(size_t pTbin) const {
     //   return _PsiSlot[pTbin];
     // }
 
 
   protected:
 
     /// Apply the projection to the event.
     void project(const Event& e);
 
     /// Compare projections.
     CmpState compare(const Projection& p) const;
 
 
   private:
 
     /// @name Jet shape parameters
     //@{
 
     /// Vector of radius bin edges
     vector<double> _binedges;
 
     /// Lower and upper cuts on contributing jet \f$ p_\perp \f$.
     pair<double, double> _ptcuts;
 
     /// Lower and upper cuts on contributing jet (pseudo)rapidity.
     pair<double, double> _rapcuts;
 
     /// Rapidity scheme
     RapScheme _rapscheme;
 
     //@}
 
 
     /// @name The projected jet shapes
     //@{
 
     /// Jet shape histo -- first index is jet number, second is r bin
     vector< vector<double> > _diffjetshapes;
 
     //@}
 
   };
 
 
 }
 
 #endif
diff --git a/include/Rivet/Projections/LossyFinalState.hh b/include/Rivet/Projections/LossyFinalState.hh
--- a/include/Rivet/Projections/LossyFinalState.hh
+++ b/include/Rivet/Projections/LossyFinalState.hh
@@ -1,81 +1,81 @@
 // -*- C++ -*-
 #ifndef RIVET_LossyFinalState_HH
 #define RIVET_LossyFinalState_HH
 
 #include "Rivet/Tools/Logging.hh"
 #include "Rivet/Config/RivetCommon.hh"
 #include "Rivet/Particle.hh"
 #include "Rivet/Event.hh"
 #include "Rivet/Projection.hh"
 #include "Rivet/Projections/FinalState.hh"
 
 namespace Rivet {
 
 
   /// @brief Templated FS projection which can lose some of the supplied particles.
   template <typename FILTER>
   class LossyFinalState : public FinalState {
   public:
 
     /// @name Constructors
     //@{
 
     /// Constructor from FinalState.
     LossyFinalState(const FinalState& fsp, FILTER filter)
       : _filter(filter)
     {
       setName("LossyFinalState");
       declare(fsp, "FS");
     }
 
     /// Stand-alone constructor. Initialises the base FinalState projection.
     LossyFinalState(FILTER filter,
-                    double mineta = -MAXDOUBLE,
-                    double maxeta = MAXDOUBLE,
+                    double mineta = -DBL_MAX,
+                    double maxeta = DBL_MAX,
                     double minpt = 0.0)
       : _filter(filter)
     {
       setName("LossyFinalState");
       declare(FinalState(mineta, maxeta, minpt), "FS");
     }
 
     /// Virtual destructor, to allow subclassing
     virtual ~LossyFinalState() { }
 
     /// Clone on the heap.
     DEFAULT_RIVET_PROJ_CLONE(LossyFinalState);
 
     //@}
 
 
     /// Apply the projection on the supplied event.
     void project(const Event& e) {
       const FinalState& fs = applyProjection<FinalState>(e, "FS");
       getLog() << Log::DEBUG << "Pre-loss number of FS particles = " << fs.particles().size() << '\n';
       _theParticles.clear();
       std::remove_copy_if(fs.particles().begin(), fs.particles().end(),
                           std::back_inserter(_theParticles), _filter);
       getLog() << Log::DEBUG << "Filtered number of FS particles = " << _theParticles.size() << '\n';
     }
 
 
     /// Compare projections.
     CmpState compare(const Projection& p) const {
       const LossyFinalState<FILTER>& other = pcast< LossyFinalState<FILTER> >(p);
       const CmpState fscmp = mkNamedPCmp(other, "FS");
       if (fscmp != CmpState::EQ) return fscmp;
       return _filter.compare(other._filter);
     }
 
 
   protected:
 
     /// Filtering object: must support operator(const Particle&) and compare(const Filter&)
     FILTER _filter;
 
   };
 
 
 }
 
 #endif
diff --git a/include/Rivet/Projections/NonHadronicFinalState.hh b/include/Rivet/Projections/NonHadronicFinalState.hh
--- a/include/Rivet/Projections/NonHadronicFinalState.hh
+++ b/include/Rivet/Projections/NonHadronicFinalState.hh
@@ -1,50 +1,50 @@
 // -*- C++ -*-
 #ifndef RIVET_NonHadronicFinalState_HH
 #define RIVET_NonHadronicFinalState_HH
 
 #include "Rivet/Tools/Logging.hh"
 #include "Rivet/Config/RivetCommon.hh"
 #include "Rivet/Particle.hh"
 #include "Rivet/Event.hh"
 #include "Rivet/Projection.hh"
 #include "Rivet/Projections/FinalState.hh"
 
 namespace Rivet {
 
 
   /// @brief Project only hadronic final state particles.
   class NonHadronicFinalState : public FinalState {
   public:
 
     /// Constructor: the supplied FinalState projection is assumed to live through the run.
     NonHadronicFinalState(FinalState& fsp)
     {
       setName("NonHadronicFinalState");
       declare(fsp, "FS");
     }
 
-    NonHadronicFinalState(double mineta = -MAXDOUBLE,
-			  double maxeta = MAXDOUBLE,
+    NonHadronicFinalState(double mineta = -DBL_MAX,
+			  double maxeta = DBL_MAX,
 			  double minpt = 0.0*GeV)
     {
       setName("NonHadronicFinalState");
       declare(FinalState(mineta, maxeta, minpt), "FS");
     }
 
     /// Clone on the heap.
     DEFAULT_RIVET_PROJ_CLONE(NonHadronicFinalState);
 
 
     /// Apply the projection on the supplied event.
     void project(const Event& e);
 
     /// Compare projections.
     CmpState compare(const Projection& p) const;
 
   };
 
 
 }
 
 
 #endif
diff --git a/include/Rivet/Projections/VisibleFinalState.hh b/include/Rivet/Projections/VisibleFinalState.hh
--- a/include/Rivet/Projections/VisibleFinalState.hh
+++ b/include/Rivet/Projections/VisibleFinalState.hh
@@ -1,55 +1,55 @@
 // -*- C++ -*-
 #ifndef RIVET_VisibleFinalState_HH
 #define RIVET_VisibleFinalState_HH
 
 #include "Rivet/Tools/Logging.hh"
 #include "Rivet/Config/RivetCommon.hh"
 #include "Rivet/Particle.hh"
 #include "Rivet/Event.hh"
 #include "Rivet/Projection.hh"
 #include "Rivet/Projections/FinalState.hh"
 
 namespace Rivet {
 
 
   /// @brief Final state modifier excluding particles which are not experimentally visible
   class VisibleFinalState : public FinalState {
   public:
 
     /// @name Constructors
     //@{
 
     /// Constructor with min and max pseudorapidity \f$ \eta \f$ and min \f$ p_T \f$ (in GeV).
-    VisibleFinalState(double mineta = -MAXDOUBLE,
-                      double maxeta =  MAXDOUBLE,
+    VisibleFinalState(double mineta = -DBL_MAX,
+                      double maxeta =  DBL_MAX,
                       double minpt  =  0.0*GeV)
     {
       setName("VisibleFinalState");
       declare(FinalState(mineta, maxeta, minpt), "FS");
     }
 
     /// Constructor with specific FinalState.
     VisibleFinalState(const FinalState& fsp)
     {
       setName("VisibleFinalState");
       declare(fsp, "FS");
     }
 
     /// Clone on the heap.
     DEFAULT_RIVET_PROJ_CLONE(VisibleFinalState);
 
     //@}
 
 
     /// Apply the projection on the supplied event.
     void project(const Event& e);
 
     /// Compare projections.
     CmpState compare(const Projection& p) const;
 
   };
 
 
 }
 
 #endif
diff --git a/src/Projections/FinalState.cc b/src/Projections/FinalState.cc
--- a/src/Projections/FinalState.cc
+++ b/src/Projections/FinalState.cc
@@ -1,101 +1,101 @@
 // -*- C++ -*-
 #include "Rivet/Projections/FinalState.hh"
 
 namespace Rivet {
 
 
   FinalState::FinalState(const Cut& c)
     : ParticleFinder(c)
   {
     setName("FinalState");
     const bool isopen = (c == Cuts::open());
     MSG_TRACE("Check for open FS conditions: " << std::boolalpha << isopen);
     if (!isopen) declare(FinalState(), "OpenFS");
   }
 
 
   FinalState::FinalState(const FinalState& fsp, const Cut& c)
     : ParticleFinder(c)
   {
     setName("FinalState");
     MSG_TRACE("Registering base FSP as 'PrevFS'");
     declare(fsp, "PrevFS");
   }
 
 
   /// @deprecated, keep for backwards compatibility for now.
   FinalState::FinalState(double mineta, double maxeta, double minpt) {
     setName("FinalState");
     const bool openpt = isZero(minpt);
-    const bool openeta = (mineta <= -MAXDOUBLE && maxeta >= MAXDOUBLE);
+    const bool openeta = (mineta <= -DBL_MAX && maxeta >= DBL_MAX);
     MSG_TRACE("Check for open FS conditions:" << std::boolalpha << " eta=" << openeta << ", pt=" << openpt);
     if (openpt && openeta) {
       _cuts = Cuts::open();
     } else {
       declare(FinalState(), "OpenFS");
       if (openeta)
         _cuts = (Cuts::pT >= minpt);
       else if ( openpt )
         _cuts = Cuts::etaIn(mineta, maxeta);
       else
         _cuts = (Cuts::etaIn(mineta, maxeta) && Cuts::pT >= minpt);
     }
   }
 
 
   CmpState FinalState::compare(const Projection& p) const {
     const FinalState& other = dynamic_cast<const FinalState&>(p);
     // First check if there is a PrevFS and it it matches
     if (hasProjection("PrevFS") != other.hasProjection("PrevFS")) return CmpState::UNDEF;
     if (hasProjection("PrevFS")) {
       const PCmp prevcmp = mkPCmp(other, "PrevFS");
       if (prevcmp != CmpState::EQ) return prevcmp;
     }
     // Then check the extra cuts
     const bool cutcmp = _cuts == other._cuts;
     MSG_TRACE(_cuts << " VS " << other._cuts << " -> EQ == " << std::boolalpha << cutcmp);
     if (!cutcmp) return CmpState::UNDEF;
     // Checks all passed: these FSes are equivalent
     return CmpState::EQ;
   }
 
 
   void FinalState::project(const Event& e) {
     _theParticles.clear();
 
     // Handle "open FS" special case, which should not/cannot recurse
     if (_cuts == Cuts::OPEN) {
       MSG_TRACE("Open FS processing: should only see this once per event (" << e.genEvent()->event_number() << ")");
       for (const GenParticle* p : Rivet::particles(e.genEvent())) {
         if (p->status() == 1) {
           MSG_TRACE("FS GV = " << p->production_vertex());
           _theParticles.push_back(Particle(*p));
         }
       }
       return;
     }
 
     // Base the calculation on PrevFS if available, otherwise OpenFS
     /// @todo In general, we'd like to calculate a restrictive FS based on the most restricted superset FS.
     const Particles& allstable = applyProjection<FinalState>(e, (hasProjection("PrevFS") ? "PrevFS" : "OpenFS")).particles();
     for (const Particle& p : allstable) {
       const bool passed = accept(p);
       MSG_TRACE("Choosing: ID = " << p.pid()
                 << ", pT = " << p.pT()/GeV << " GeV"
                 << ", eta = " << p.eta()
                 << ": result = " << std::boolalpha << passed);
       if (passed) _theParticles.push_back(p);
     }
     MSG_TRACE("Number of final-state particles = " << _theParticles.size());
   }
 
 
   /// Decide if a particle is to be accepted or not.
   bool FinalState::accept(const Particle& p) const {
     // Not having status == 1 should never happen!
     assert(p.genParticle() == NULL || p.genParticle()->status() == 1);
     return _cuts->accept(p);
   }
 
 
 }