UnitConst.cc

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
 
#include <framework/gearbox/Unit.h>
#include <framework/gearbox/Const.h>
#include <framework/logging/Logger.h>
#include <framework/particledb/EvtGenDatabasePDG.h>
 
#include <TMath.h>
 
#include <algorithm>
#include <fstream>
 
using namespace Belle2;
using namespace std;
 
/*** The implementation of the Unit class defined in Unit.h starts here ***/
 
namespace Belle2 {
  map<string, double> Unit::s_conversionFactors;
 
#define DEFINE_UNIT_NAME(var,value,name) const double Unit::var = Unit::registerConversion(name,value)
#define DEFINE_UNIT(var,value) DEFINE_UNIT_NAME(var,value,#var)
 
  DEFINE_UNIT(cm, 1.);     
  DEFINE_UNIT(ns, 1.);     
  DEFINE_UNIT(GHz, 1.);    
  DEFINE_UNIT(rad, 1.);    
  DEFINE_UNIT(GeV, 1.);    
  DEFINE_UNIT(K, 1.);      
  DEFINE_UNIT(e, 1.);      
  DEFINE_UNIT_NAME(g_cm3, 1., "g/cm3");  
// length units
  DEFINE_UNIT(km, Unit::cm * 1e5);     
  DEFINE_UNIT(m, Unit::cm * 1e2);      
  DEFINE_UNIT(mm, Unit::m  * 1e-3);    
  DEFINE_UNIT(um, Unit::m  * 1e-6);    
  DEFINE_UNIT(nm, Unit::m  * 1e-9);    
  DEFINE_UNIT(pm, Unit::m  * 1e-12);   
  DEFINE_UNIT(fm, Unit::m  * 1e-15);   
// area units
  DEFINE_UNIT(m2, Unit::m * Unit::m);        
  DEFINE_UNIT(cm2, Unit::cm * Unit::cm);     
  DEFINE_UNIT(mm2, Unit::mm * Unit::mm);     
  DEFINE_UNIT(b, Unit::m2 * 1e-28);    
  DEFINE_UNIT(mb, Unit::b  * 1e-3);    
  DEFINE_UNIT(ub, Unit::b  * 1e-6);    
  DEFINE_UNIT(nb, Unit::b  * 1e-9);    
  DEFINE_UNIT(pb, Unit::b  * 1e-12);   
  DEFINE_UNIT(fb, Unit::b  * 1e-15);   
  DEFINE_UNIT(ab, Unit::b  * 1e-18);   
// volume units
  DEFINE_UNIT(m3, Unit::m * Unit::m * Unit::m);      
  DEFINE_UNIT(cm3, Unit::cm * Unit::cm * Unit::cm);  
  DEFINE_UNIT(mm3, Unit::mm * Unit::mm * Unit::mm);  
// time units
  DEFINE_UNIT(s, Unit::ns * 1e9);      
  DEFINE_UNIT(ms, Unit::s  * 1e-3);    
  DEFINE_UNIT(us, Unit::s  * 1e-6);    
  DEFINE_UNIT(ps, Unit::s  * 1e-12);   
  DEFINE_UNIT(fs, Unit::s  * 1e-15);   
// frequency units
  DEFINE_UNIT(Hz,   Unit::GHz * 1e-9);        
  DEFINE_UNIT(kHz,  Unit::GHz * 1e-6);        
  DEFINE_UNIT(MHz,  Unit::GHz * 1e-3);        
  DEFINE_UNIT(mHz,  Unit::GHz * 1e-12);       
// angle units
  DEFINE_UNIT(mrad, Unit::rad * 1e-3);   
  DEFINE_UNIT(deg, TMath::DegToRad());   
// energy units
  DEFINE_UNIT(eV, Unit::GeV * 1e-9);   
  DEFINE_UNIT(keV, Unit::eV  * 1e3);   
  DEFINE_UNIT(MeV, Unit::eV  * 1e6);   
  DEFINE_UNIT(TeV, Unit::eV  * 1e12);   
  DEFINE_UNIT(J, 1.0 / TMath::Qe() * Unit::eV);    
  DEFINE_UNIT(V, Unit::eV / Unit::e);    
// charge units
  DEFINE_UNIT(C, 1.0 / TMath::Qe() * Unit::e);   
  DEFINE_UNIT(fC, Unit::C / 1.0e15);             
// magnetic field units
  DEFINE_UNIT(T, Unit::V * Unit::s / Unit::m2);   
  DEFINE_UNIT(mT, Unit::T * 1e-3);                
  DEFINE_UNIT(Gauss, Unit::T * 1e-4);             
  DEFINE_UNIT(kGauss, Unit::Gauss * 1e3);         
// density units
  DEFINE_UNIT_NAME(mg_cm3, Unit::g_cm3 * 1e-3, "mg/cm3");   
  DEFINE_UNIT_NAME(kg_cm3, Unit::g_cm3 * 1e3, "kg/cm3");    
  DEFINE_UNIT_NAME(g_mm3, Unit::g_cm3 / Unit::mm3, "g/mm3");     
  DEFINE_UNIT_NAME(mg_mm3, Unit::mg_cm3 / Unit::mm3, "mg/mm3");   
  DEFINE_UNIT_NAME(kg_mm3, Unit::kg_cm3 / Unit::mm3, "kg/mm3");   
  double Unit::convertValue(double value, const std::string& unitString)
  {
    auto it = s_conversionFactors.find(unitString);
    if (it == s_conversionFactors.end()) {
      B2ERROR("Could not find conversion factor for unit " << unitString << ", assuming 1.0");
      return value;
    }
    return it->second * value;
  }
 
  double Unit::convertValueToUnit(double value, const std::string& unitString)
  {
    auto it = s_conversionFactors.find(unitString);
    if (it == s_conversionFactors.end()) {
      B2ERROR("Could not find conversion factor for unit " << unitString << ", assuming 1.0");
      return value;
    }
    return value / it->second;
  }
 
  double Unit::registerConversion(const string& name, double value)
  {
    if (!s_conversionFactors.insert(make_pair(name, value)).second) {
      B2ERROR("Conversion factor for unit " << name << " already exists, ignoring");
    }
    return value;
  }
 
}
 
/*** The implementation of the Const class defined in Const.h starts here ***/
 
std::string Const::parseDetectors(EDetector det)
{
  if (det == Const::PXD)   return "PXD";
  else if (det == Const::SVD)   return "SVD";
  else if (det == Const::CDC)   return "CDC";
  else if (det == Const::TOP)   return "TOP";
  else if (det == Const::ARICH) return "ARICH";
  else if (det == Const::ECL)   return "ECL";
  else if (det == Const::KLM)   return "KLM";
  else if (det == Const::BKLM)  return "BKLM";
  else if (det == Const::EKLM)  return "EKLM";
  else if (det == Const::IR)    return "IR";
  else if (det == Const::TRG)   return "TRG";
  else if (det == Const::DAQ)   return "DAQ";
  else if (det == Const::BEAST) return "BEAST";
  else if (det == Const::TEST)  return "TEST";
  else B2FATAL("Unknown detector component: " << det);
  return "INVALID";
}
 
Const::DetectorSet::Iterator& Const::DetectorSet::Iterator::operator++()
{
  while (1) {
    m_SetBit = m_SetBit << 1;
    if (m_SetBit >= 0x1000) {
      m_SetBit = invalidDetector;
      return *this;
    }
    if ((m_DetectorSetBits & m_SetBit) != 0)
      break;
  }
  m_Index++;
  return *this;
}
 
bool Const::DetectorSet::Iterator::operator==(
  const Const::DetectorSet::Iterator& iterator)
{
  return m_SetBit == iterator.m_SetBit;
}
 
bool Const::DetectorSet::Iterator::operator!=(
  const Const::DetectorSet::Iterator& iterator)
{
  return m_SetBit != iterator.m_SetBit;
}
 
Const::DetectorSet::Iterator Const::DetectorSet::begin() const
{
  uint16_t setBit = 1;
  while ((m_bits & setBit) == 0) {
    setBit = setBit << 1;
    if (setBit >= 0x1000)
      return Const::DetectorSet::Iterator(0, m_bits, invalidDetector);
  }
  return Const::DetectorSet::Iterator(0, m_bits, setBit);
}
 
Const::DetectorSet::Iterator Const::DetectorSet::end() const
{
  return Const::DetectorSet::Iterator(0, m_bits, invalidDetector);
}
 
Const::DetectorSet operator + (const Const::DetectorSet& firstSet, const Const::DetectorSet& secondSet)
{
  Const::DetectorSet set(firstSet);
  set += secondSet;
  return set;
}
 
Const::DetectorSet operator - (const Const::DetectorSet& firstSet, const Const::DetectorSet& secondSet)
{
  Const::DetectorSet set(firstSet);
  set -= secondSet;
  return set;
}
 
Const::DetectorSet operator + (Const::EDetector firstDet, Const::EDetector secondDet)
{
  Const::DetectorSet set(firstDet);
  set += secondDet;
  return set;
}
 
unsigned short Const::DetectorSet::getBit(Const::EDetector det)
{
  if (det == 0 || det > 0x100) {
    B2ERROR("Const::DetectorSet::getBit() called for an invalid detector "
            "identifier (possibly a subdetector):" << det);
    return 0;
  }
  if (det > TEST) {
    B2ERROR("Const::DetectorSet::getBit(): Invalid detector ID");
  }
  return (1 << (det - 1));
}
 
Const::EDetector Const::DetectorSet::getDetector(unsigned short bit)
{
  switch (bit) {
    case 0x0001: return PXD;
    case 0x0002: return SVD;
    case 0x0004: return CDC;
    case 0x0008: return TOP;
    case 0x0010: return ARICH;
    case 0x0020: return ECL;
    case 0x0040: return KLM;
    case 0x0080: return IR;
    case 0x0100: return TRG;
    case 0x0200: return DAQ;
    case 0x0400: return BEAST;
    case 0x0800: return TEST;
    default: return invalidDetector;
  }
}
 
int Const::DetectorSet::getIndex(EDetector det) const
{
  unsigned short bit = getBit(det);
  if ((m_bits & bit) == 0) return -1;
  int index = 0;
  for (unsigned short setBit = 1; setBit < bit; setBit *= 2) {
    if ((m_bits & setBit) != 0) ++index;
  }
  return index;
}
 
size_t Const::DetectorSet::size() const
{
  int size = 0;
  for (unsigned short setBit = 1; setBit < 0x1000; setBit *= 2) {
    if ((m_bits & setBit) != 0) ++size;
  }
  return size;
}
 
std::string Const::DetectorSet::__str__() const
{
  std::string result = "<set: ";
  const std::string detectorNames[] = {"invalidDetector", "PXD", "SVD", "CDC", "TOP", "ARICH", "ECL", "KLM", "IR", "TRG", "DAQ", "BEAST", "TEST"};
  for (size_t index = 1; index <= Const::TEST; index++) {
    if (contains(EDetector(index))) {
      if (result.size() > 6) {
        result += ",";
      }
      result += detectorNames[index];
    }
  }
  result += ">";
  return result;
}
 
 
const Const::DetectorSet Const::VXDDetectors::c_set = Const::PXD + Const::SVD;
 
const Const::DetectorSet Const::TrackingDetectors::c_set = Const::PXD + Const::SVD + Const::CDC;
 
const Const::DetectorSet Const::PIDDetectors::c_set = Const::SVD + Const::CDC + Const::TOP + Const::ARICH + Const::ECL + Const::KLM;
 
 
const Const::DetectorSet Const::allDetectors = Const::IR + Const::PXD + Const::SVD + Const::CDC +
                                               Const::TOP + Const::ARICH + Const::ECL + Const::KLM +
                                               Const::TRG + Const::DAQ + Const::BEAST + Const::TEST;
 
const Const::DetectorSet Const::ClusterDetectors::c_set = Const::ECL + Const::KLM;
 
bool Const::ParticleType::operator < (const Const::ParticleType& other) const
{
  return m_pdgCode < other.m_pdgCode;
}
 
Const::ParticleType& Const::ParticleType::operator++()
{
  if (!m_set) {
    *this = invalidParticle;
  } else {
    m_pdgCode = m_set->at(++m_index).getPDGCode();
  }
 
  return *this;
}
 
Const::ParticleType Const::ParticleType::operator++(int)
{
  Const::ParticleType p = *this;
  ++(*this);
  return p;
}
 
const TParticlePDG* Const::ParticleType::getParticlePDG() const
{
  return EvtGenDatabasePDG::Instance()->GetParticle(m_pdgCode);
}
 
double Const::ParticleType::getMass() const
{
  return getParticlePDG()->Mass();
}
 
std::string Const::ParticleType::__repr__() const
{
  std::string result = "<type: ";
  result += getParticlePDG()->GetName();
  result += ">";
  return result;
}
 
const Const::ParticleSet Const::chargedStableSet =
  Const::ParticleType(11) + Const::ParticleType(13) + Const::ParticleType(211)
  + Const::ParticleType(321) + Const::ParticleType(2212) + Const::ParticleType(1000010020);
 
const Const::ParticleSet Const::clusterSet =
  Const::ParticleType(22) + Const::ParticleType(130) + Const::ParticleType(11)
  + Const::ParticleType(211) + Const::ParticleType(13) + Const::ParticleType(9900001);
 
 
const Const::ChargedStable Const::electron = Const::chargedStableSet.find(11);
const Const::ChargedStable Const::muon = Const::chargedStableSet.find(13);
const Const::ChargedStable Const::pion = Const::chargedStableSet.find(211);
const Const::ChargedStable Const::kaon = Const::chargedStableSet.find(321);
const Const::ChargedStable Const::proton = Const::chargedStableSet.find(2212);
const Const::ChargedStable Const::deuteron = Const::chargedStableSet.find(1000010020);
 
const Const::Cluster Const::clusterphoton = Const::clusterSet.find(22);
const Const::Cluster Const::clusterKlong = Const::clusterSet.find(130);
const Const::Cluster Const::clusterelectron = Const::clusterSet.find(11);
const Const::Cluster Const::clusterpion = Const::clusterSet.find(211);
const Const::Cluster Const::clustermuon = Const::clusterSet.find(13);
const Const::Cluster Const::clusterjunk = Const::clusterSet.find(9900001);
 
const Const::ParticleType Const::photon = Const::ParticleType(22);
const Const::ParticleType Const::pi0 = Const::ParticleType(111);
const Const::ParticleType Const::neutron = Const::ParticleType(2112);
const Const::ParticleType Const::antiNeutron = Const::ParticleType(-2112);
const Const::ParticleType Const::Kshort = Const::ParticleType(310);
const Const::ParticleType Const::Klong = Const::ParticleType(130);
const Const::ParticleType Const::Lambda = Const::ParticleType(3122);
const Const::ParticleType Const::antiLambda = Const::ParticleType(-3122);
const Const::ParticleType Const::invalidParticle = Const::ParticleType(9900000); // codes beginning with 99... are reserved
const Const::ParticleType Const::unspecifiedParticle = Const::ParticleType(9900001); // codes beginning with 99... are reserved
 
const Const::ParticleSet Const::finalStateParticlesSet =
  Const::chargedStableSet + Const::photon + Const::Kshort + Const::Klong
  + Const::Lambda + Const::antiLambda + Const::neutron;
 
const double Const::electronMass = Const::electron.getMass();
const double Const::muonMass = Const::muon.getMass();
const double Const::pionMass = Const::pion.getMass();
const double Const::kaonMass = Const::kaon.getMass();
const double Const::protonMass = Const::proton.getMass();
const double Const::deuteronMass = Const::deuteron.getMass();
const double Const::pi0Mass = Const::pi0.getMass();
const double Const::neutronMass = Const::neutron.getMass();
const double Const::K0Mass = Const::Kshort.getMass();
 
void Const::ParticleSet::add(const Const::ParticleType& p)
{
  if (contains(p))
    return;
  m_particles.emplace_back(p.getPDGCode(), this, m_particles.size());
}
 
bool Const::ParticleSet::contains(const Const::ParticleType& p) const
{
  return (std::find(m_particles.begin(), m_particles.end(), p) != m_particles.end());
}
 
Const::ParticleSet operator + (const Const::ParticleSet& firstSet, const Const::ParticleSet& secondSet)
{
  Const::ParticleSet result(firstSet);
  for (Const::ParticleType pdgIter = secondSet.begin(); pdgIter != secondSet.end(); ++pdgIter) {
    result.add(pdgIter);
  }
  return result;
}
 
 
const double Const::speedOfLight   = 29.9792458;
const double Const::kBoltzmann     = 8.617343 * 1.0e-5 * Unit::eV / Unit::K;
const double Const::ehEnergy       = 3.65 * Unit::eV;
const double Const::fineStrConst   = 1.0 / 137.036;
const double Const::permSi         = 11.9 * 8.8542 * 1e-18 * Unit::C / Unit::V / Unit::um;
const double Const::uTherm         = 0.026 * Unit::V;
const double Const::eMobilitySi    = 1415 * Unit::cm2 / Unit::V / Unit::s;
 
const double Const::doubleNaN = std::numeric_limits<double>::quiet_NaN();