{{lastupdated_at}} by {{lastupdated_by}}

{{>toc}}

h1. Fermi/LAT interface

h2. Pass 7 IRF

Pass 7 IRF have been implement in GammaLib version 00-06-00. Please check [[Validation of Pass 7 IRFs]] for more information.


h2. Livetime cubes

The Fermi/LAT Science Tools implement the livetime cube information as @map_tools::Exposure@ objects. @map_tools::Exposure@ derives from @SkyExposure@ which is defined as
<pre>
typedef BasicExposure<SkyBinner, healpix::CosineBinner> SkyExposure;
</pre>
The @BasicExposure@ class is a template class which has the following access operator that is used to retrieve exposure values:
<pre>
    template<class F>
        double operator()(const astro::SkyDir& dir, const F& fun)const
    {
        const C& binner = m_sky[dir];
        return binner(fun);
    }
</pre>
Here @C@ is the @healpix::CosineBinner@ class. The @class F@ is typically the effective area class. The @binner(fun)@ operator is implemented in @healpix/healpix/CosineBinner@ by
<pre>
    template<class F>
    double operator()(const F& f)const
    {   
        double sum=0;
        int n(0);
        for(const_iterator it=begin(); it!=end_costh(); ++it, ++n){
            sum += (*it)*f(costheta(it));
        }
        return sum; 

    }
</pre>
The energy does not appear explicitly in these code fragments as it is stored as a member of the effective area helper class:
<pre>
   class Aeff : public AeffBase {
   public:
      Aeff(double energy, int evtType, const Observation & observation); 
      virtual ~Aeff() {}
   protected:
      double m_energy;
      int m_evtType;
      const Observation & m_observation;
      virtual double value(double cosTheta, double phi=0) const;
   };
</pre>


h2. Exposure computation using efficiency factors

The Fermi/LAT Science Tools compute the exposure with trigger rate- and energy-dependent efficiency corrections in @Likelihood/Likelihood/ExposureCube.h@ where the following template class is implemented:
<pre>
   template<class T>
   double value(const astro::SkyDir & dir, const T & aeff, double energy) const {
      double factor1(1), factor2(0);
      if (m_efficiencyFactor) {
         double met((m_tstart + m_tstop)/2.);
         m_efficiencyFactor->getLivetimeFactors(energy, factor1, factor2, met);
      }
      double value1(value(dir, aeff));
      double value2(0);
      double exposure(factor1*value1);
      if (factor2 != 0) {
         value2 = value(dir, aeff, true);
         exposure += factor2*value2;
      }
      if (exposure < 0) {
         throw std::runtime_error("ExposureCube::value: exposure < 0");
      }
      return exposure;
   }
</pre>
An approximation is made in the computation of the efficiency factors which are the same for the front and the back section (see @irfs/latResponse/src/EfficiencyFactor.cxx@):
<pre>
void EfficiencyFactor::getLivetimeFactors(double energy, double & factor1, double & factor2, double met) const {
   (void)(met);
   if (!m_havePars) {
      factor1 = 1;
      factor2 = 0;
      return;
   }
   double logE(std::log10(energy));

   // Since we do not have access to the front and back effective
   // areas separately, just return the averages.  We would really
   // want to return the averages weighted by effective area (as a function
   // of energy and theta).
   factor1 = (m_p1_front(logE) + m_p1_back(logE))/2.;
   factor2 = (m_p0_front(logE) + m_p0_back(logE))/2.;
}
</pre>

In GammaLib we do not make this approximation as the response is fully separated in front and back response functions. The efficiency factors are returned for front and back separately:
<pre>
double GLATLtCube::operator()(const GSkyDir& dir, const GEnergy& energy, const GLATAeff& aeff) const
{
    // Compute exposure
    double exposure = m_exposure(dir, energy, aeff);

    // Optionally compute livetime factors for trigger rate- and
    // energy-dependent efficiency corrections
    if (aeff.has_efficiency()) {
    
        // Compute correction factors
        double f1 = aeff.efficiency_factor1(energy);
        double f2 = aeff.efficiency_factor2(energy);

        // Compute correction
        double correction = m_weighted_exposure(dir, energy, aeff);

        // Set exposure
        exposure = f1 * exposure + f2 * correction;

    } // endif: corrections requested

    // Return exposure
    return exposure;
}
</pre>
Below an example for the efficiency factors for @P7REP_SOURCE_V15@ at 200 MeV:

| *Parameter* | *Science Tools* | *GammaLib* |
| factor1::front | -1.06527 | -1.05736 |
| factor2::front | 2.28282 | 2.27408 |
| factor1::back | -1.06527 | -1.07317 |
| factor2::back | 2.28282 | 2.29157 |


h2. Effective area

The effective area is accessed in the Fermi/LAT Science Tools using the @ExposureCube::Aeff::value@ method (recall that the energy is stored in the member @m_energy@):
<pre>
double ExposureCube::Aeff::value(double cosTheta, double phi) const {
   double inclination = acos(cosTheta)*180./M_PI;
   double epoch((m_observation.expCube().tstart() + m_observation.expCube().tstop())/2.);
   std::map<unsigned int, irfInterface::Irfs *>::const_iterator respIt = m_observation.respFuncs().begin();
   for ( ; respIt != m_observation.respFuncs().end(); ++respIt) {
      if (respIt->second->irfID() == m_evtType) {
         irfInterface::IAeff * aeff = respIt->second->aeff();
         double aeff_val = aeff->value(m_energy, inclination, phi, epoch);
         return aeff_val;
      }
   }
   return 0;
}
</pre>The @inclination@ is given in degrees.

The effective area is implemented by the @irfs/latResponse/Aeff.cxx@ class. The @value@ method is used to access the effective area values:
<pre>
double Aeff::value(double energy, double theta, double phi, double time) const {
   (void)(time);
   double costheta(std::cos(theta*M_PI/180.));
   if (costheta < m_aeffTable.minCosTheta()) {
      return 0;
   }
   if (costheta > 0.99999) {
      costheta = 0.99999;  // blech.
   }
   bool interpolate;
   double logE(std::log10(energy));
   double aeff_value(m_aeffTable.value(logE, costheta, interpolate=true)*1e4);
   double phi_mod(phi_modulation(logE, costheta, phi, interpolate=false));
   return aeff_value*phi_mod;
}
</pre>Note the unnecessary conversion of @inclination@ to @costheta@ (why is not simply @costheta@ passed in this interface?). Effective area values are determined by bilinear interpolation in logE and costheta. The code also multiplies-in a Phi modulation:
<pre>
double Aeff::phi_modulation(double logE, double costheta, double phi, bool interpolate) const {
   double par[2];
   if (!m_phiDepPars || !m_usePhiDependence) {
      return 1.;
   }
   m_phiDepPars->getPars(logE, costheta, par, interpolate);
   return phi_modulation(par[0], par[1], phi);
}
double Aeff::phi_modulation(double par0, double par1, double phi) const {
   if (phi < 0) {
      phi += 360.;
   }
   double norm(1./(1. + par0/(1. + par1)));
   double phi_pv(std::fmod(phi*M_PI/180., M_PI) - M_PI/2.);
   double xx(2.*std::fabs(2./M_PI*std::fabs(phi_pv) - 0.5));
   return norm*(1. + par0*std::pow(xx, par1));
}
</pre>

The equivalent GammaLib code is
<pre>
double GLATAeff::operator()(const double& logE, const double& ctheta)
{
    // Get effective area value
    double aeff = (ctheta >= m_min_ctheta) 
                  ? m_aeff_bins.interpolate(logE, ctheta, m_aeff) : 0.0;

    // Make sure that effective area is not negative
    if (aeff < 0.0) {
        aeff = 0.0;
    }
    
    // Return effective area value
    return aeff;
}
</pre>A variant of this operator takes also a phi value, yet *the effective area phi dependence is so far not implemented in GammaLib*.

Here a number of effective area values (in cm2) determined using the ScienceTools and GammaLib for the @P7REP_SOURCE_V15@ response for 200 MeV as function of costheta (no phi dependence was used in the Fermi/LAT Science Tools):
| *logE* | *costheta* | *ScienceTools* | *GammaLib* |
| 2.30103 | 0.999844 | 2334.4 | 2426.55 |
| 2.30103 | 0.998594 | 2334.4 | 2417.21 |
| 2.30103 | 0.996094 | 2334.4 | 2398.55 |
| 2.30103 | 0.992344 | 2334.4 | 2370.56 |
| 2.30103 | 0.987344 | 2333.23 | 2333.23 |
| 2.30103 | 0.981094 | 2286.58 | 2286.58 |
| 2.30103 | 0.418594 | 119.825 | 119.825 |
| 2.30103 | 0.379844 | 60.7843 | 60.7843 |
| 2.30103 | 0.339844 | 26.3115 | 26.3115 |
| 2.30103 | 0.298594 | 8.99221 | 8.99221 |
| 2.30103 | 0.256094 | 2.41828 | 2.41828 |
| 2.30103 | 0.212344 | 0.571005 | 0.566958 |

For costheta > 0.992344 (about 7°) the effective area in the Fermi/LAT Science Tools saturates while for GammaLib it still increases (maybe the ScienceTools bi-linear interpolator does not extrapolate?). The Fermi/LAT Science Tools determine the minimum costheta angle from the lower boundary of the front response (note that formerly this was fixed to 70°, but I now dropped this since Fermi/LAT Science Tools do no limit the response).

h2. Mean PSF

The Fermi/LAT Science Tools and GammaLib use a mean PSF for point sources in binned analysis. The mean PSF is implemented by @Likelihood/src/MeanPsf.cxx@:
<pre>
void MeanPsf::init() {
   computeExposure();
   if (s_separations.size() == 0) {
      createLogArray(1e-4, 70., 200, s_separations);
   }
   m_psfValues.reserve(m_energies.size()*s_separations.size());
   for (unsigned int k = 0; k < m_energies.size(); k++) {
      for (unsigned int j = 0; j < s_separations.size(); j++) {
         double value(0);
         std::map<unsigned int, irfInterface::Irfs *>::const_iterator 
            resp = m_observation.respFuncs().begin();
         for (; resp != m_observation.respFuncs().end(); ++resp) {
            int evtType = resp->second->irfID();
            Psf psf(s_separations[j], m_energies[k], evtType, m_observation);
            value += m_observation.expCube().value(m_srcDir, psf,
                                                   m_energies[k]);
         }
         if (m_exposure[k] > 0) {
            value /= m_exposure[k];
         } else {
            value = 0;
         }
         m_psfValues.push_back(value);
       }
   }
}
void MeanPsf::computeExposure() {
   m_exposure.reserve(m_energies.size());
   for (size_t k(0); k < m_energies.size(); k++) {
      double value(0);
      std::map<unsigned int, irfInterface::Irfs *>::const_iterator
         resp = m_observation.respFuncs().begin();
      for (; resp != m_observation.respFuncs().end(); ++resp) {
         int evtType = resp->second->irfID();
         ExposureCube::Aeff aeff(m_energies[k], evtType, m_observation);
         value += m_observation.expCube().value(m_srcDir, aeff, m_energies[k]);
      }
      m_exposure.push_back(value);
   }
}
</pre>