Long_term_correction_header.h

Go to the documentation of this file.

#ifndef Long_term_correction_header_h
#define Long_term_correction_header_h

/*  by Jakub Vicha
 *  slightly modified to be compiled within offline by Max Stadelmaier
 *
 *  long-term corrections for energ(FD)/energy(SD)
 *  probably this should be implemented by an FD expert somewhere else
 */

unsigned int DataRelease = 1, NumberOfBreaks = 2;

bool ReadNewLTparametrization = false;

double LongTermCorrectionPerEye( double EsdOverFd, double gpsmegatime, unsigned int eye )       // Esd / Efd, GPS time in Ms and zenith-angle bin
{
        const unsigned int nFDsites = 4; // LL, LM, LA, CO

        double GPSbreak[nFDsites], a1[nFDsites], a2[nFDsites], b[nFDsites];
        double GPSbreak1_2breaks[nFDsites], GPSbreak2_2breaks[nFDsites], a1_2breaks[nFDsites], a2_2breaks[nFDsites], a3_2breaks[nFDsites], b_2breaks[nFDsites];
        double CorrectedRatio, f; //, fRef;

        if( DataRelease == 0 )
        {
// above log Efd = 18.5 - ICRC 17

// GPS time [Ms] of long-term break
                GPSbreak[0] = 1059.67;
                GPSbreak[1] = 1004.65;
                GPSbreak[2] = 1063.43;
                GPSbreak[3] = 1001.87;

                a1[0] = 0.000301326;
                a1[1] = 0.00111561;
                a1[2] = 0.000521306;
                a1[3] = 0.000681541;

                a2[0] = 0.000240983;
                a2[1] = -0.000198824;
                a2[2] = -0.00170036;
                a2[3] = -0.00078575;

                b[0] = 0.707388;
                b[1] = -0.0298003;
                b[2] = 0.495244;
                b[3] = 0.336456;
        }

        if( DataRelease == 1 )
        {
                // above log Efd = 18.5 - ICRC 19 + 2018 + Core < 800m & angle < 5deg & LDF Cuts > 10^-10

                // GPS time [Ms] of long-term break
                GPSbreak[0] = 961.821;
                GPSbreak[1] = 985.511;
                GPSbreak[2] = 1047.96;
                GPSbreak[3] = 959.097;

                a1[0] = 0.000229442;
                a1[1] = 0.00110734;
                a1[2] = 0.00040724;
                a1[3] = 0.00068503;

                a2[0] = 0.000325319;
                a2[1] = 0.000192236;
                a2[2] = -0.000279451;
                a2[3] = -0.000187543;

                b[0] = 0.760063;
                b[1] = -0.0304074;
                b[2] = 0.60066;
                b[3] = 0.334217;


                // GPS time [Ms] of long-term break
                GPSbreak1_2breaks[0] = 1073.59;
                GPSbreak1_2breaks[1] = 1044.57;
                GPSbreak1_2breaks[2] = 1066.6;
                GPSbreak1_2breaks[3] = 1042.24;

                GPSbreak2_2breaks[0] = 1128.46;
                GPSbreak2_2breaks[1] = 1098.82;
                GPSbreak2_2breaks[2] = 1092.46;
                GPSbreak2_2breaks[3] = 1095;

                a1_2breaks[0] = 0.00040397;
                a1_2breaks[1] = 0.00103914;
                a1_2breaks[2] = 0.000437;
                a1_2breaks[3] = 0.000509163;

                a2_2breaks[0] = -0.00114323;
                a2_2breaks[1] = -0.00147673;
                a2_2breaks[2] = -0.00203909;
                a2_2breaks[3] = -0.00175069;

                a3_2breaks[0] = 0.00177153;
                a3_2breaks[1] = 0.00101904;
                a3_2breaks[2] = 3.87678e-05;
                a3_2breaks[3] = 0.000419606;

                b_2breaks[0] = -1.03302;
                b_2breaks[1] = -0.0856851;
                b_2breaks[2] = 0.944868;
                b_2breaks[3] = 0.466709;
        }

        if( NumberOfBreaks == 1 )
        {
                if( gpsmegatime < GPSbreak[eye] )
                {
                        f = a1[eye] * gpsmegatime + b[eye];

                        CorrectedRatio = EsdOverFd / f;
                }
                else
                {
                        f = a2[eye] * gpsmegatime + b[eye] + (a1[eye]-a2[eye]) * GPSbreak[eye];

                        CorrectedRatio = EsdOverFd / f;
                }
        }

        if( NumberOfBreaks == 2 )
        {
                if( gpsmegatime < GPSbreak1_2breaks[eye] )
                {
                        f = a1_2breaks[eye] * gpsmegatime + b_2breaks[eye] + ( a3_2breaks[eye] - a2_2breaks[eye] ) * GPSbreak2_2breaks[eye] + ( a2_2breaks[eye] - a1_2breaks[eye] ) * GPSbreak1_2breaks[eye]; // a1*x + b3 + (a3-a2)*X2 + (a2-a1)*X1

                        CorrectedRatio = EsdOverFd / f;
                }
                if( gpsmegatime > GPSbreak1_2breaks[eye] && gpsmegatime < GPSbreak2_2breaks[eye] )
                {
                        f = a2_2breaks[eye] * gpsmegatime + b_2breaks[eye] + (a3_2breaks[eye]-a2_2breaks[eye]) * GPSbreak2_2breaks[eye]; //y = a2*x + b3 + (a3-a2)*X2

                        CorrectedRatio = EsdOverFd / f;
                }
                if( gpsmegatime > GPSbreak2_2breaks[eye] )
                {
                        f = a3_2breaks[eye] * gpsmegatime + b_2breaks[eye]; //y = a3*x + b3

                        CorrectedRatio = EsdOverFd / f;
                }
        }

        return CorrectedRatio;
}

//=======================================================================================================================================


double LongTermCorrectionPerEyeAY( double EsdOverFd, double gpsmegatime, unsigned int eye )     // Esd / Efd, GPS time in Ms and zenith-angle bin
{
  const double MoonCycleLength  = 29.532;
  //Jan 1st 2004
  const double MoonCycleJan1st2004 = 11.45 - 344/MoonCycleLength;
  const int GPSJan1st2004 = 756950413;

  const double x = MoonCycleJan1st2004 + (gpsmegatime*1e6 - GPSJan1st2004)/(MoonCycleLength*24*3600);//Moon cycle

  double constant;
  double breaks[2];
  double slope[3];

  double LL_par[2] = {0, 0};    //linear part of fit >130 MoonCycles for LL
  double LL_stitch = 0;
  //fits
  switch (eye)
    {
    case 0:
      //Eye 1 - LL
      //special fit with 2 breaks in range < 130 && then linear
      constant = 0.99413;
      slope[0] = 0.000557776;
      slope[1] =  0.00536182;
      slope[2] = -0.0321092;
      breaks[0] = 100.198;
      breaks[1] = 124.123;

      //linear part > 130
      LL_par[0] = 9.33706e-01;
      LL_par[1] = 3.40356e-03;
      LL_stitch = 130;

      //default 2 break fit
      /*constant = 1.04467;
      slope[0] = 0.00106613;
      slope[1] = -0.00376464;
      slope[2] = 0.00295298;
      breaks[0] = 123.426;
      breaks[1] = 140.68;
      LL_stitch = 10000;
      */
      break;

    case 1:
      //Eye 2 - LM
      constant = 1.09387;
      slope[0] = 0.00205574;
      slope[1] = -0.000290713;
      slope[2] = 0.00323317;
      breaks[0] = 111.5;
      breaks[1] = 155.503;
      break;

    case 2:
      //Eye 3 - LA
      constant = 1.06842;
      slope[0] = 0.00222544;
      slope[1] = -0.00306826;
      slope[2] = 7.98857e-05;
      breaks[0] = 113.499;
      breaks[1] = 138.242;
      break;

    case 3:
      //Eye 4 - Co
      constant = 1.04585;
      slope[0] = 0.00134737;
      slope[1] = -0.0961635;
      slope[2] = 0.00152404;
      breaks[0] = 124.918;
      breaks[1] = 126.311;
      break;

    default:
    constant = 1.;
    slope[0] = 0.0;
    slope[1] = -0.0;
    slope[2] = 0.0;
    breaks[0] = 100.;
    breaks[1] = 100.;
    break;

    }

  double f =  x < breaks[1]? (x < breaks[0] ? constant + slope[0] * (x - breaks[0]) : constant + slope[1] * (x - breaks[0])) : constant + slope[1] * (breaks[1] - breaks[0]) + slope[2] * (x - breaks[1]);

  if (x > LL_stitch && eye == 0)
    f = LL_par[0] + LL_par[1] * (x - LL_stitch);

  return EsdOverFd / f;
}


#endif // #ifdef Long_term_correction_header_cxx