GCOMTools.cpp

Go to the documentation of this file.

/***************************************************************************
 *                     GCOMTools.cpp - COMPTEL tools                       *
 * ----------------------------------------------------------------------- *
 *  copyright (C) 2021-2025 by Juergen Knoedlseder                         *
 * ----------------------------------------------------------------------- *
 *                                                                         *
 *  This program is free software: you can redistribute it and/or modify   *
 *  it under the terms of the GNU General Public License as published by   *
 *  the Free Software Foundation, either version 3 of the License, or      *
 *  (at your option) any later version.                                    *
 *                                                                         *
 *  This program is distributed in the hope that it will be useful,        *
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
 *  GNU General Public License for more details.                           *
 *                                                                         *
 *  You should have received a copy of the GNU General Public License      *
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
 *                                                                         *
 ***************************************************************************/
/**
 * @file GCOMTools.hpp
 * @brief Implementation of COMPTEL tools
 * @author Juergen Knoedlseder
 */
 
/* __ Includes ___________________________________________________________ */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include "GTime.hpp"
#include "GEnergy.hpp"
#include "GNodeArray.hpp"
#include "GTools.hpp"
#include "GCOMTools.hpp"
 
/* __ Method name definitions ____________________________________________ */
#define G_COM_TOFCOR      "com_tofcor(GEnergy&, GEnergy&, int&, int& tofmax)"
 
/* __ Coding definitions _________________________________________________ */
 
/* __ Debug definitions __________________________________________________ */
 
 
/***********************************************************************//**
 * @brief Convert TJD and COMPTEL tics in GTime object
 *
 * @param[in] tjd Truncated Julian Days (days).
 * @param[in] tics COMPTEL tics (1/8 ms).
 * @return Time.
 *
 * Converts TJD and COMPTEL tics into a GTime object. COMPTEL times are
 * given in UTC, i.e. 8393:0 converts into 1991-05-17T00:00:00 UT
 * (see COM-RP-UNH-DRG-037).
 *
 * The method applies the CGRO clock correction. The CGRO clock was too
 * fast by 2.042144 seconds before 8798:28800000 (1992-06-25T01:00:00).
 * After that date the clock was corrected.
 ***************************************************************************/
GTime gammalib::com_time(const int& tjd, const int& tics)
{
    // Compute tics to days conversion factor
    const double tics2days     = 0.000125 * gammalib::sec2day;
    const double clockcor_days = 2.042144 * gammalib::sec2day;
 
    // Compute MJD
    double mjd = double(tjd) + 40000.0 + double(tics) * tics2days;
 
    // Apply CGRO clock correction as the clock was too fast by 2.042144 sec
    // before 8798:28800000 (1992-06-25T01:00:00).
    if ((tjd < 8798) || ((tjd == 8798) && (tics < 28800000))) {
        mjd -= clockcor_days;
    }
 
    // Set time in MJD in the UTC time system
    GTime time;
    time.mjd(mjd, "UTC");
 
    // Return time
    return time;
}
 
 
/***********************************************************************//**
 * @brief Convert GTime in COMPTEL TJD
 *
 * @param[in] time Time.
 * @return Truncated Julian Days (days).
 *
 * Converts GTime object in TJD.
 *
 * The method applies the CGRO clock correction. The CGRO clock was too
 * fast by 2.042144 seconds before 8798:28800000 (1992-06-25T01:00:00).
 * After that date the clock was corrected.
 ***************************************************************************/
int gammalib::com_tjd(const GTime& time)
{
    // Set MJD in UTC of clock correction
    const double mjd_of_clockcor = 48798.04166666666424134746;
    const double clockcor_days   = 2.042144 * gammalib::sec2day;
    const double half_tics2days  = 0.5 * 0.000125 * gammalib::sec2day;
 
    // Compute MJD, applying the CGRO clock correction before 8798:28800000
    // (1992-06-25T01:00:00)
    double mjd = time.mjd("UTC");
    if (mjd < mjd_of_clockcor) {
        mjd += clockcor_days;
    }
 
    // Compute TJD, rounding to the nearest tics
    int tjd = int(mjd + half_tics2days) - 40000;
 
    // Return TJD
    return tjd;
}
 
 
/***********************************************************************//**
 * @brief Convert GTime in COMPTEL tics
 *
 * @param[in] time Time.
 * @return COMPTEL tics (1/8 ms).
 *
 * Converts GTime object in COMPTEL tics (1/8 ms).
 *
 * The method applies the CGRO clock correction. The CGRO clock was too
 * fast by 2.042144 seconds before 8798:28800000 (1992-06-25T01:00:00).
 * After that date the clock was corrected.
 ***************************************************************************/
int gammalib::com_tics(const GTime& time)
{
    // Set MJD in UTC of clock correction
    const double mjd_of_clockcor = 48798.04166666666424134746;
    const double clockcor_days   = 2.042144 * gammalib::sec2day;
    const int    tics_in_day     = 691200000;
 
    // Compute MJD, applying the CGRO clock correction before 8798:28800000
    // (1992-06-25T01:00:00)
    double mjd = time.mjd("UTC");
    if (mjd < mjd_of_clockcor) {
        mjd += clockcor_days;
    }
 
    // Compute number of tics, rounding to nearest int
    int tics = int((mjd - double(int(mjd))) * tics_in_day + 0.5);
 
    // If number of tics exceeds number of tics in one day then decrement
    // number of tics
    while (tics >= tics_in_day) {
        tics -= tics_in_day;
    }
 
    // Return tics
    return tics;
}
 
 
/***********************************************************************//**
 * @brief Compute ToF correction
 *
 * @param[in] emin Minimum energy.
 * @param[in] emax Maximum energy.
 * @param[in] tofmin Minimum ToF channel.
 * @param[in] tofmax Maximum ToF channel.
 * @return Time of flight correction.
 *
 * Compute the ToF correction according to COM-RP-ROL-DRG-057.
 ***************************************************************************/
double gammalib::com_tofcor(const GEnergy& emin,   const GEnergy& emax,
                            const int&     tofmin, const int&     tofmax)
{
    // ToF correction factors according to Table 1 of COM-RP-ROL-DRG-057
    const double  tofcoreng[] = {0.8660, 1.7321, 5.4772, 17.3205};
    const double  tofcor110[] = {1.14,   1.07,   1.02,    1.01};
    const double  tofcor111[] = {1.17,   1.09,   1.03,    1.01};
    const double  tofcor112[] = {1.21,   1.11,   1.05,    1.02};
    const double  tofcor113[] = {1.26,   1.15,   1.07,    1.04};
    const double  tofcor114[] = {1.32,   1.20,   1.11,    1.06};
    const double  tofcor115[] = {1.40,   1.27,   1.17,    1.11};
    const double  tofcor116[] = {1.50,   1.36,   1.24,    1.17};
    const double  tofcor117[] = {1.63,   1.47,   1.35,    1.28};
    const double  tofcor118[] = {1.79,   1.63,   1.51,    1.43};
    const double  tofcor119[] = {2.01,   1.85,   1.73,    1.67};
    const double* tofcor[]    = {tofcor110, tofcor111, tofcor112,
                                 tofcor113, tofcor114, tofcor115,
                                 tofcor116, tofcor117, tofcor118,
                                 tofcor119};
 
    // Initialise ToF correction
    double tofcorr = 1.0;
 
    // Determine ToF correction
    if (tofmin < 110) {
        std::string msg = "Minimum of ToF selection window "+
                          gammalib::str(tofmin)+
                          " is smaller than 110. No ToF correction is "
                          "available for this value.";
        gammalib::warning(G_COM_TOFCOR, msg);
    }
    else if (tofmin > 119) {
        std::string msg = "Minimum of ToF selection window "+
                          gammalib::str(tofmin)+
                          " is larger than 119. No ToF correction is "
                          "available for this value.";
        gammalib::warning(G_COM_TOFCOR, msg);
    }
    else if (tofmax != 130) {
        std::string msg = "Maximum of ToF selection window "+
                          gammalib::str(tofmax)+
                          " is not 130. No ToF correction is available for "
                          "this value.";
        gammalib::warning(G_COM_TOFCOR, msg);
    }
    else {
 
        // Compute mean energy
        double energy = std::sqrt(emin.MeV() * emax.MeV());
 
        // Set node array
        GNodeArray nodes(4, tofcoreng);
 
        // Set interpolation value
        nodes.set_value(energy);
 
        // Get ToF correction index
        int i = tofmin - 110;
 
        // Interpolate
        tofcorr = tofcor[i][nodes.inx_left()]  * nodes.wgt_left() +
                  tofcor[i][nodes.inx_right()] * nodes.wgt_right();
 
    }
 
    // Return ToF correction
    return tofcorr;
}