BeamConstraint.hh

// -*- C++ -*-
#ifndef RIVET_BeamConstraint_HH
#define RIVET_BeamConstraint_HH
 
#include "Rivet/Tools/Beams.hh"
#include "Rivet/Tools/ParticleName.hh"
 
namespace Rivet {
 
 
 
 
  inline bool compatibleBeamID(PdgId p, PdgId allowed) {
    return (allowed == PID::ANY || p == allowed);
  }
 
  inline bool compatibleBeamIDs(const PdgIdPair& pair, const PdgIdPair& allowedpair) {
    const bool oneToOne = compatibleBeamID(pair.first, allowedpair.first);
    const bool twoToTwo = compatibleBeamID(pair.second, allowedpair.second);
    const bool oneToTwo = compatibleBeamID(pair.first, allowedpair.second);
    const bool twoToOne = compatibleBeamID(pair.second, allowedpair.first);
    return (oneToOne && twoToTwo) || (oneToTwo && twoToOne);
  }
 
  inline bool compatibleBeamIDs(const ParticlePair& ppair, const PdgIdPair& allowedpair) {
    return compatibleBeamIDs(PID::make_pdgid_pair(ppair.first.pid(), ppair.second.pid()), allowedpair);
  }
 
  inline bool compatibleBeamIDs(const ParticlePair& a, const ParticlePair& b) {
    return compatibleBeamIDs(PID::make_pdgid_pair(a.first.pid(), a.second.pid()),
                             PID::make_pdgid_pair(b.first.pid(), b.second.pid()) );
  }
 
 
  inline bool compatibleBeamEnergies(const pair<double,double>& energies,
                                     const pair<double,double>& allowedenergies, const double reltol=1e-3) {
    const bool oneToOne = fuzzyEquals(energies.first,  allowedenergies.first,  reltol);
    const bool twoToTwo = fuzzyEquals(energies.second, allowedenergies.second, reltol);
    const bool oneToTwo = fuzzyEquals(energies.first,  allowedenergies.second, reltol);
    const bool twoToOne = fuzzyEquals(energies.second, allowedenergies.first,  reltol);
    const bool absDiffOneToOne = abs(energies.first - allowedenergies.first)   < 1*GeV;
    const bool absDiffTwoToTwo = abs(energies.second - allowedenergies.second) < 1*GeV;
    const bool absDiffOneToTwo = abs(energies.first - allowedenergies.second)  < 1*GeV;
    const bool absDiffTwoToOne = abs(energies.second - allowedenergies.first)  < 1*GeV;
    return (oneToOne && twoToTwo) || (oneToTwo && twoToOne) ||
           (absDiffOneToOne && absDiffTwoToTwo) || (absDiffOneToTwo && absDiffTwoToOne);
  }
 
  inline bool compatibleBeamEnergies(const ParticlePair& ppair,
                                     const pair<double,double>& allowedenergies, double reltol=1e-3) {
    return compatibleBeamEnergies(make_pair(ppair.first.E(), ppair.second.E()), allowedenergies, reltol);
  }
 
  inline bool compatibleBeamEnergies(const ParticlePair& ppair,
                                     const ParticlePair& allowedppair, double reltol=1e-3) {
    return compatibleBeamEnergies(make_pair(ppair.first.E(), ppair.second.E()),
                                  make_pair(allowedppair.first.E(), allowedppair.second.E()), reltol);
  }
 
 
  inline bool compatibleBeamEnergy(const pair<double,double>& energies, double allowedsqrts, double reltol=1e-3) {
    return fuzzyEquals(sqrtS(energies.first, energies.second), allowedsqrts, reltol);
  }
 
  inline bool compatibleBeamEnergy(const ParticlePair& ppair,
                                   const pair<double,double>& allowedenergies, double reltol=1e-3) {
    return fuzzyEquals(sqrtS(ppair), sqrtS(allowedenergies), reltol);
  }
 
  inline bool compatibleBeamEnergy(const ParticlePair& ppair,
                                   const double allowedsqrts, double reltol=1e-3) {
    return fuzzyEquals(sqrtS(ppair), allowedsqrts, reltol);
  }
 
 
  inline bool compatibleBeams(const ParticlePair& ppair, const PdgIdPair& allowedpids, double allowedsqrts) {
    return compatibleBeamIDs(ppair, allowedpids) && compatibleBeamEnergy(ppair, allowedsqrts);
  }
 
  inline bool compatibleBeams(const ParticlePair& ppair, const PdgIdPair& allowedpids, const pair<double,double>& allowedenergies) {
    return compatibleBeamIDs(ppair, allowedpids) && compatibleBeamEnergies(ppair, allowedenergies);
  }
 
  inline bool compatibleBeams(const ParticlePair& a, const ParticlePair& b) {
    return compatibleBeamIDs(a, b) && compatibleBeamEnergies(a, b);
  }
 
 
 
 
  // /// Find whether PdgIdPair @a pair is compatible with at least one template beam pair in a set @a allowedpairs
  // inline bool anyCompatibleBeamIDs(const PdgIdPair& pair, const set<PdgIdPair>& allowedpairs) {
  //   for (const PdgIdPair& bp : allowedpairs)
  //     if (compatibleBeamIDs(pair, bp)) return true;
  //   return false;
  // }
 
  // /// Return the intersection of two sets of {PdgIdPair}s.
  // inline set<PdgIdPair> intersection(const set<PdgIdPair>& a, const set<PdgIdPair>& b) {
  //   set<PdgIdPair> ret;
  //   for (set<PdgIdPair>::const_iterator bp = a.begin(); bp != a.end(); ++bp) {
  //     if (compatibleBeamIDs(*bp, b)) ret.insert(*bp);
  //   }
  //   return ret;
  // }
 
 
 
}
 
#endif