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