eic-opticks/sysrap/sdomain.h

0001 #pragma once
0002 /**
0003 sdomain.h
0004 ===========
0005
0006
0007 Regarding hc_eVnm see U4PhysicalConstantsTest,
0008 there is slight difference from smath.h float value : could be arising from CLHEP version difference
0009
0010 Wavelength fine domain np.linspace(60,820,761)
0011
0012 **/
0013
0014
0015 #include <string>
0016 #include <sstream>
0017 #include <cassert>
0018 #include <iomanip>
0019
0020 #include "NPFold.h"
0021
0022
0023 struct sdomain
0024 {
0025     static constexpr double hc_eVnm = 1239.84198433200208455673 ;
0026     static constexpr const double    DOMAIN_LOW  = 60. ;
0027     static constexpr const double    DOMAIN_HIGH = 820. ;
0028     static constexpr const double    COARSE_DOMAIN_STEP = 20. ;
0029     static constexpr const int       COARSE_DOMAIN_LENGTH = 39 ;
0030     static constexpr const double    FINE_DOMAIN_STEP = 1. ;
0031     static constexpr const int       FINE_DOMAIN_LENGTH = 761 ;  // 820-60+1
0032
0033     static constexpr const char     DOMAIN_TYPE = 'F' ;   // 'C'
0034
0035     static constexpr double DomainLow(){    return DOMAIN_LOW ; }
0036     static constexpr double DomainHigh(){   return DOMAIN_HIGH ; }
0037     static constexpr double DomainStep(){   return DOMAIN_TYPE == 'F' ? FINE_DOMAIN_STEP    : COARSE_DOMAIN_STEP ; }
0038     static constexpr double DomainRange(){  return DOMAIN_HIGH - DOMAIN_LOW ; }
0039
0040     static constexpr int    DomainLength(){  return DOMAIN_TYPE == 'F' ? FINE_DOMAIN_LENGTH : COARSE_DOMAIN_LENGTH ; }
0041
0042     sdomain();
0043
0044     NP* get_wavelength_nm() const ;
0045     NP* get_energy_eV() const ;
0046     NPFold* get_fold() const ;
0047
0048     static std::string Desc(const double* vv, int length, int edge);
0049     std::string desc() const ;
0050
0051     int length ;
0052     double step ;
0053     double* wavelength_nm ;
0054     double* energy_eV ;      // HMM: energy_eV  is descending following ascending wavelength_nm
0055     double* spec4 ;
0056 };
0057
0058 inline sdomain::sdomain()
0059     :
0060     length(DomainLength()),
0061     step(DomainStep()),
0062     wavelength_nm(new double[length]),
0063     energy_eV(new double[length]),
0064     spec4(new double[4])
0065 {
0066     for(int i=0 ; i < length ; i++) wavelength_nm[i] = DOMAIN_LOW + step*double(i) ;
0067     assert( wavelength_nm[0] == DOMAIN_LOW );
0068     assert( wavelength_nm[length-1] == DOMAIN_HIGH );
0069     for(int i=0 ; i < length ; i++) energy_eV[i] = hc_eVnm/wavelength_nm[i] ;
0070
0071     spec4[0] = DomainLow();
0072     spec4[1] = DomainHigh();
0073     spec4[2] = DomainStep();
0074     spec4[3] = DomainRange();
0075 }
0076
0077
0078
0079 inline NP* sdomain::get_wavelength_nm() const
0080 {
0081     return NP::MakeFromValues<double>( wavelength_nm, length ) ;
0082 }
0083 inline NP* sdomain::get_energy_eV() const
0084 {
0085     return NP::MakeFromValues<double>( energy_eV, length ) ;
0086 }
0087 inline NPFold* sdomain::get_fold() const
0088 {
0089     NPFold* fold = new NPFold ;
0090     fold->add("wavelength_nm", get_wavelength_nm() );
0091     fold->add("energy_eV", get_energy_eV() );
0092     return fold ;
0093 }
0094
0095
0096
0097 inline std::string sdomain::Desc(const double* vv, int length, int edge)
0098 {
0099     std::stringstream ss ;
0100     ss << "(" ;
0101     for(int i=0 ; i < length ; i++)
0102     {
0103         if( i < edge || i > length - edge )
0104             ss << std::fixed << std::setw(10) << std::setprecision(5) << vv[i] << " " ;
0105         else if( i == edge )
0106             ss << "... " ;
0107     }
0108     ss << ")" ;
0109     std::string s = ss.str();
0110     return s ;
0111 }
0112 inline std::string sdomain::desc() const
0113 {
0114     int edge = 5 ;
0115     std::stringstream ss ;
0116     ss
0117         << "sdomain::desc"
0118         << " length " << length  << std::endl
0119         << " wavelength_nm " << Desc(wavelength_nm, length, edge) << std::endl
0120         << " energy_eV     " << Desc(energy_eV,     length, edge) << std::endl
0121         ;
0122     std::string s = ss.str();
0123     return s ;
0124 }
0125