PRD2020/UnfoldUtilities.hh

Go to the documentation of this file.

/*  Utilities to  unfold the vertical spectrum such as done for ICRC 2015  */

/* Author: I. Valiño, based on functions provided by A. Schulz */

string CodeUsed;
map <string, Double_t> mapParameters;
string Verbose = "no";
string SDFile;
string OutPutFileName;
string OutPutFileNameRebinned;
string Model;
string MatrixType;
string ResolutionFlag;
string BiasFlag;
string EfficiencyFlag;
string ShowerToShowerFlag;
string VersionInput = "standard";
string UseNewTriggers = "no";
string PrintCompact = "no";
TVectorD GetSDCalPars() {
  //Calibration constants
  TVectorD pCal(2);
  pCal[0] = mapParameters["A"];
  pCal[1] = mapParameters["B"];
  return pCal;
}

double InverseSdCalibration(const double lgE, TVectorD pCal) {
  return (lgE - log10(pCal[0])) / pCal[1]; /*  returns lgS38 or lgS35 */
}

double ScaledErrorFunction(const double x, const std::vector<double>& pars) {
  return 0.5 * (1 + TMath::Erf((x - pars[0] ) / pars[1]));
}

// TRIGGER EFFICIENCY
double SDTriggerEfficiency_AveragedZenithAngle(const double lgE) {
//The functions can be used integrated in zenith angle
  const double thmax = mapParameters["MaxZenith"] * TMath::DegToRad();
  const double thmin = 0.;
  const double thStep = 0.0005; //0.0005 delivers a result consistent with root methods within 0.01%.
  double integralTheta = 0;
  const TVectorD pCal = GetSDCalPars();
  const bool oldTrig = ("no" == UseNewTriggers);

  for (double thIntegr = thmin; thIntegr < thmax; thIntegr += thStep) {
    if (EfficiencyFlag == "ICRC17_SD1500")
      integralTheta += SDTriggerEfficiencyICRC2017_SD1500(lgE, thIntegr, pCal) * thStep * TMath::Cos(thIntegr) * TMath::Sin(thIntegr);
    if (EfficiencyFlag == "ICRC17_SD750")
      integralTheta += SDTriggerEfficiencyICRC2017_SD750(lgE, thIntegr, pCal) * thStep * TMath::Cos(thIntegr) * TMath::Sin(thIntegr);
    if (EfficiencyFlag == "ICRC19_SD750")
      integralTheta += SDTriggerEfficiencyICRC2019_SD750(lgE, thIntegr, pCal, oldTrig) * thStep * TMath::Cos(thIntegr) * TMath::Sin(thIntegr);
    if (EfficiencyFlag == "ICRC19_SD1500_June2019")
      integralTheta += SD1500TriggerEfficiency_ICRC2019(lgE, thIntegr) * thStep * TMath::Cos(thIntegr) * TMath::Sin(thIntegr);
  }
  double fnorm = 0.5 * (pow(TMath::Sin(thmax), 2) - pow(TMath::Sin(thmin), 2));
  return integralTheta / fnorm;
}

double EfficiencyValues(double lgE, double zenith) {
  const TVectorD pCal = GetSDCalPars();
  double efficiency = 1;
  const bool oldTrig = ("no" == UseNewTriggers);

  if (EfficiencyFlag == "ICRC17_SD1500")
    efficiency = SDTriggerEfficiencyICRC2017_SD1500(lgE, zenith, pCal);
  else if (EfficiencyFlag == "ICRC17_SD750")
    efficiency = SDTriggerEfficiencyICRC2017_SD750(lgE, zenith, pCal);
  else if (EfficiencyFlag == "ICRC19_SD1500_May2019")
    efficiency = SD1500TriggerEfficiencyHybrids2019(lgE);
  else if (EfficiencyFlag == "ICRC19_SD1500_June2019")
    efficiency = SD1500TriggerEfficiency_ICRC2019(lgE, zenith);
  else  if (EfficiencyFlag == "ICRC19_SD750")
    efficiency = SDTriggerEfficiencyICRC2019_SD750(lgE, zenith, pCal, oldTrig);
  else  if ( EfficiencyFlag == "TriggerEfficiency_TestExample" )
    efficiency = TriggerEfficiency_TestExample(lgE);
  else {
    cerr << endl;
    cerr << EfficiencyFlag << ": no valid efficiency parameterization chosen " << endl;
    throw std::exception();
  }
  return efficiency;
}

double EfficiencyValuesNoDeclDep(double lgE) {
  //If values are zenith angle dependent average over the entire zenith angle range
  double efficiency = 1;
  if (EfficiencyFlag == "ICRC17_SD1500" ||
      EfficiencyFlag == "ICRC17_SD750" ||
      EfficiencyFlag == "ICRC19_SD1500_June2019" ||
      EfficiencyFlag == "ICRC19_SD750")
    efficiency = SDTriggerEfficiency_AveragedZenithAngle(lgE);
  else if (EfficiencyFlag == "ICRC19_SD1500_May2019")
    efficiency = SD1500TriggerEfficiencyHybrids2019(lgE);
  else  if ( EfficiencyFlag == "TriggerEfficiency_TestExample" )
    efficiency = TriggerEfficiency_TestExample(lgE);
  else {
    cerr << endl;
    cerr << EfficiencyFlag << ": no valid efficiency parameterization chosen " << endl;
    throw std::exception();
  }

  return efficiency;
}

// BIAS
TVectorD EnergyBias(TVectorD lgEs, TVectorD pCal) {
  //used by the non-declination dependent matrix calculation
  int ncol = lgEs.GetNoElements();
  TVectorD biasEnergy(ncol);

  for (int i = 0; i < ncol; ++i) {
    if (BiasFlag == "BiasICRC19_SD1500_preliminary") {
      //calculates the average over cos2theta of bias for each energy
      double integral = 0.;
      double deltacos2theta = 0.001;
      for (Double_t cos2theta = 0.25; cos2theta  < 1; cos2theta += deltacos2theta) {
        double costheta = sqrt(cos2theta);
        integral += (HybridBias_2019_SD1500(lgEs[i], costheta, mapParameters["MaxZenith"]) - 1.) * deltacos2theta;
      }
      integral /= 0.75; //normalization
      integral += 1.;

      if ( lgEs[i] > 18.4)
        integral = 1.;
      biasEnergy[i] = integral;
    } else if (BiasFlag == "BiasICRC19_SD1500_Offline") {
      //calculates the average over cos2theta of bias for each energy
      double integral = 0.;
      double deltacos2theta = 0.001;
      for (Double_t cos2theta = 0.25; cos2theta  < 1; cos2theta += deltacos2theta) {
        double costheta = sqrt(cos2theta);
        integral += (HybridBias_2019_11_06_Offline_SD1500(lgEs[i], costheta, mapParameters["MaxZenith"]) - 1.) * deltacos2theta;
      }
      integral /= 0.75; //normalization
      integral += 1.;

      if ( lgEs[i] > 18.4)
        integral = 1.;
      biasEnergy[i] = integral;
    } else if (BiasFlag == "BiasICRC19_SD1500_Herald") {
      //calculates the average over cos2theta of bias for each energy
      double integral = 0.;
      double deltacos2theta = 0.001;
      for (Double_t cos2theta = 0.25; cos2theta  < 1; cos2theta += deltacos2theta) {
        double costheta = sqrt(cos2theta);
        integral += (HybridBias_2019_11_06_Herald_SD1500(lgEs[i], costheta, mapParameters["MaxZenith"]) - 1.) * deltacos2theta;
      }
      integral /= 0.75; //normalization
      integral += 1.;

      if ( lgEs[i] > 18.4)
        integral = 1.;
      biasEnergy[i] = integral;
    } else if (BiasFlag == "BiasICRC19_SD750") {
      const double fdB = pCal[1];
      biasEnergy[i] = HybridBias_2019_SD750(lgEs[i], fdB);
    } else if (BiasFlag == "Bias_TestExample") {
      //something...
      //......
      biasEnergy[i] = 1.;
    } else if (BiasFlag == "TOZERO")
      biasEnergy[i] = 1.;
    else {
      cerr << endl;
      cerr << BiasFlag << ": no valid bias parameterization chosen " << endl;
      throw std::exception();
    }
  }

  return biasEnergy;
}

double EnergyBias_DD(double lgE, double cosTheta, TVectorD pCal) {
  double bias = 1;
  if (BiasFlag == "BiasICRC19_SD1500_preliminary")
    bias = HybridBias_2019_SD1500(lgE, cosTheta, mapParameters["MaxZenith"]);
  else if (BiasFlag == "BiasICRC19_SD1500_Offline")
    bias = HybridBias_2019_11_06_Offline_SD1500(lgE, cosTheta, mapParameters["MaxZenith"]);
  else if (BiasFlag == "BiasICRC19_SD1500_Herald")
    bias = HybridBias_2019_11_06_Herald_SD1500(lgE, cosTheta, mapParameters["MaxZenith"]);
  else if (BiasFlag == "BiasICRC19_SD750") 
   {
    const double fdB = pCal[1];
    bias = HybridBias_2019_SD750(lgE, fdB);
   }
  else if (BiasFlag == "Bias_TestExample")
    bias = Bias_TestExample(lgE, cosTheta);
  else if (BiasFlag == "TOZERO")
    bias = 1;
  else {
    cerr << endl;
    cerr << BiasFlag << ": no valid bias parameterization chosen " << endl;
    throw std::exception();
  }

  return bias;
}

//  SHOWER TO SHOWER FLUCTUATIONS
double ShowerToShowerFluctuation(const double lgE) {
  double shsh = 0;
  if (ShowerToShowerFlag == "ICRC17_SD1500Data") {
    const double xmin = 16.5;
    const double xmax = 20.5;
    const double z = (lgE - xmin) / (xmax - xmin);
    const double pars[2] = {0.2, 0.0783};
    shsh = (1 - z) * pars[0] + z * pars[1];
  } else if (ShowerToShowerFlag == "ICRC17_SD750Data")
    shsh = 0.13;
  else if (ShowerToShowerFlag == "TOZERO")
    shsh = 0;
  else {
    cerr << endl;
    cerr << ShowerToShowerFlag << ": no valid shower to shower parameterization chosen " << endl;
    throw std::exception();
  }
  return shsh;
}

// RESOLUTION
double ResolutionZenithDistrib(double lgE, double zenith) {
  TVectorD pCal = GetSDCalPars();//only for data resolution
  const double lgS38 = InverseSdCalibration(lgE, pCal);

  double res;
  if (ResolutionFlag == "ICRC17_SD750Data")
    res = DataSD750EnergyResolution_ICRC2017(lgE, zenith, pCal, lgS38);
  else if (ResolutionFlag == "ICRC17_SD1500Data")
    res = DataSD1500EnergyResolution(lgE, zenith, pCal, lgS38);
  else if (ResolutionFlag == "ICRC17_SD750Simu")
    res = SimulationSD750EnergyResolution_ICRC2017(lgE, zenith);
  else if (ResolutionFlag == "ICRC17_SD1500Simu") //Zenith angle independent... Not necessary here in the zenith integration
    res = SimulationSD1500EnergyResolution(lgE, zenith);
  else if (ResolutionFlag == "Valerio12_11_2018")
    res = Valerio12_11_2018(lgE);
  else if (ResolutionFlag == "Hybrids_ICRC_2019_preliminary") //Zenith angle independent... Not necessary here in the zenith integration
    res = Hybrids_ICRC_2019_SD1500(lgE);
  else if (ResolutionFlag == "Hybrids_ICRC_2019_SD750")
    res = Hybrids_ICRC_2019_SD750(lgE);
  else if (ResolutionFlag == "Hybrids_ICRC_2019_Offline")
    res = Hybrids_ICRC_2019_Offline_11_06(lgE);
  else if (ResolutionFlag == "Hybrids_ICRC_2019_Offline_05_07")
    res = Hybrids_ICRC_2019_Offline_05_07(lgE);
  else if (ResolutionFlag == "Hybrids_ICRC_2019_Herald")
    res = Hybrids_ICRC_2019_Herald_11_06(lgE);
  else if (ResolutionFlag == "Resolution_TestExample") //Zenith angle independent... Not necessary here in the zenith integration
    res = Resolution_TestExample(lgE);
  else {
    cerr << endl;
    cerr << ResolutionFlag << ": no valid resolution parameterization chosen " << endl;
    throw std::exception();
  }

  return res;
}

double IntegrateEnergyResolution(const double lgE) {
  const double thmax = mapParameters["MaxZenith"] * TMath::DegToRad();
  const double thmin = 0.;
  const double thStep = 0.0005;
  double integralTheta = 0;
  for (double thIntegr = thmin; thIntegr < thmax; thIntegr += thStep) {
    integralTheta += ResolutionZenithDistrib(lgE, thIntegr) * thStep * TMath::Cos(thIntegr) * TMath::Sin(thIntegr);
  }
  double fnorm = 0.5 * (pow(TMath::Sin(thmax), 2) - pow(TMath::Sin(thmin), 2));
  return integralTheta / fnorm;
}

TVectorD EnergyResolution(TVectorD lgEs) {
  int ncol = lgEs.GetNoElements();
  TVectorD EnergyRes(ncol);

  for (int i = 0; i < ncol; ++i) {
    double lgE = lgEs[i];
    double Esd = pow(10., lgE);

    double relErrs1 = IntegrateEnergyResolution(lgE);
    double relErrs2 = 0;
    relErrs2 = ShowerToShowerFluctuation(lgE);

    double relErrs = TMath::Sqrt(pow(relErrs1, 2)  + pow(relErrs2, 2));

    EnergyRes[i] = relErrs * Esd;
  }

  return EnergyRes;
}

double ResolutionValue(double lgE, double_t cosTheta) {
  TVectorD pCal = GetSDCalPars();//only for data resolution
  const double lgS38 = InverseSdCalibration(lgE, pCal);

  double resolution = 0;
  if (ResolutionFlag == "ICRC17_SD1500Data")
    resolution = DataSD1500EnergyResolution(lgE, acos(cosTheta), pCal, lgS38);
  else if (ResolutionFlag == "ICRC17_SD750Data")
    resolution = DataSD750EnergyResolution_ICRC2017(lgE, acos(cosTheta), pCal, lgS38);
  else if (ResolutionFlag == "ICRC17_SD1500Simu")
    resolution = SimulationSD1500EnergyResolution(lgE, acos(cosTheta));
  else if (ResolutionFlag == "ICRC17_SD750Simu")
    resolution = SimulationSD750EnergyResolution_ICRC2017(lgE, acos(cosTheta));
  else if (ResolutionFlag == "Valerio12_11_2018")
    resolution = Valerio12_11_2018(lgE);
  else if (ResolutionFlag == "Hybrids_ICRC_2019_preliminary")
    resolution = Hybrids_ICRC_2019_SD1500(lgE);
  else if (ResolutionFlag == "Hybrids_ICRC_2019_SD750")
    resolution = Hybrids_ICRC_2019_SD750(lgE);
  else if (ResolutionFlag == "Hybrids_ICRC_2019_Offline")
    resolution = Hybrids_ICRC_2019_Offline_11_06(lgE);
  else if (ResolutionFlag == "Hybrids_ICRC_2019_Offline_05_07")
    resolution = Hybrids_ICRC_2019_Offline_05_07(lgE);
  else if (ResolutionFlag == "Hybrids_ICRC_2019_Herald")
    resolution = Hybrids_ICRC_2019_Herald_11_06(lgE);
  else if (ResolutionFlag == "Resolution_TestExample")
    resolution = Resolution_TestExample(lgE);
  else {
    cerr << endl;
    cerr << ResolutionFlag << ": no valid resolution parameterization chosen " << endl;
    throw std::exception();
  }
  return resolution;
}


// Standard migration matrix
TMatrixD kResolutionMatrix(TVectorD pCal, TVectorD xos, double* resolutionOffset) {
  int N = xos.GetNoElements();
  TVectorD ses = EnergyResolution(xos);
  TVectorD Ebias = EnergyBias(xos, pCal);
  TMatrixD tmp_kmatrix(N, N);
  for (int j = 0; j < N; ++j) {
    double ex = pow(10, xos[j]);
    double efficiency = EfficiencyValuesNoDeclDep(xos[j]);

    for (int i = 0; i < N; ++i) {
      double ey = pow(10, xos[i]);
      tmp_kmatrix[i][j] = TMath::Gaus(ey, ex * Ebias[j], ses[j], true) * ey * log(10.) * efficiency;
    }
  }

  /* ofstream outMatrixFile("matrixFile_OLD.txt");
   for (int i = 0; i < N; ++i)
    {
     for (int j = 0; j < N; ++j)
      outMatrixFile<<tmp_kmatrix[i][j]<<" ";

     outMatrixFile<<endl;
    }
  */
  return tmp_kmatrix;
}

// Zenith dependent migration matrix
TMatrixD kResolutionMatrixZenith(TVectorD pCal, TVectorD xos, double* resolutionOffset){
  Int_t diagonalOfMatrix = 40;
  Double_t fac = 1;
  if (mapParameters["IntegrationSteps"]>10) {
    fac = mapParameters["IntegrationSteps"]/10.;
    diagonalOfMatrix = diagonalOfMatrix * fac;
  }
  cout << "Declination dependent migration matrix" << endl;
  cout << "Calculates only +-" << diagonalOfMatrix << " elements out of " << xos.GetNoElements() <<  "around diagonal to speed up calculation! " << endl;
  Double_t deltaTheta = 0.5;
  Double_t maxtheta = mapParameters["MaxZenith"];
  Double_t minTheta = mapParameters["MinZenith"];

  int N = xos.GetNoElements();
  TMatrixD tmp_kmatrixTH(N, N);
  TMatrixD tmp_kmatrixTH_norm(N, N);
  for (int i = 0; i < N; ++i)
    for (int j = 0; j < N; ++j) {
      tmp_kmatrixTH_norm[i][j] = 0.;
      tmp_kmatrixTH[i][j] = 0.;
    }  
  for (int i = 0; i < N; ++i) {
    double ey = pow(10, xos[i]);//reco
    for (int j = 0; j < N; ++j) {
      double ex = pow(10, xos[j]);//real
      if (fabs(i - j) < diagonalOfMatrix)
        for (Double_t th = minTheta; th < maxtheta; th += deltaTheta)
          {
            Double_t th_half = th + 0.5*deltaTheta;
            Double_t cosTheta = cos(th_half*TMath::DegToRad());
            double resolution = ResolutionValue(xos[j], cosTheta);
            double sh_shfluct = ShowerToShowerFluctuation(xos[j]);
            double tot_sigma = sqrt(pow(resolution, 2) + pow(sh_shfluct, 2));
            double efficiency = EfficiencyValues(xos[j], th_half*TMath::DegToRad());
            double bias = EnergyBias_DD(xos[j], cosTheta, pCal);
            tmp_kmatrixTH_norm[i][j] += sin(2.* th_half*TMath::DegToRad()) * 2. * TMath::Pi() * deltaTheta;
            tmp_kmatrixTH[i][j]      += sin(2.* th_half*TMath::DegToRad()) * 2. * TMath::Pi() * deltaTheta  * efficiency * TMath::Gaus(ey, ex * bias    , tot_sigma * ex, true) * ey * log(10.);
          }
      
      if (tmp_kmatrixTH_norm[i][j] > 0)
        tmp_kmatrixTH[i][j] /= tmp_kmatrixTH_norm[i][j]; 
    }
  }

  ofstream outMatrixFile("matrixFile_FULL_TH.txt");
  for (int i = 0; i < N; ++i)
   {
    for (int j = 0; j < N; ++j)
     outMatrixFile<<tmp_kmatrixTH[i][j]<<" ";

    outMatrixFile<<endl;
   }

  return tmp_kmatrixTH;

}



// Declination dependent migration matrix
TMatrixD kResolutionMatrixDD(TVectorD pCal, TVectorD xos, Double_t lat, Double_t deltaMin, Double_t deltaMax) {
  Int_t diagonalOfMatrix = 40;
// To speed up! verify always if there is a bias...
// The edges of the diagonal should be always very small. In first approximation they can be neglected..
  Double_t fac = 1;
  if (mapParameters["IntegrationSteps"]>10) {
    fac = mapParameters["IntegrationSteps"]/10.;
    diagonalOfMatrix = diagonalOfMatrix * fac;
  }
  cout << "Declination dependent migration matrix" << endl;
  cout << "Calculates only +-" << diagonalOfMatrix << " elements out of " << xos.GetNoElements() <<  "around diagonal to speed up calculation! " << endl;
  Double_t deltah = 0.1;
  Double_t deltadelta = 0.1;
  int N = xos.GetNoElements();
  TMatrixD tmp_kmatrixDD(N, N);
  TMatrixD tmp_kmatrixDD_norm(N, N);
  for (int i = 0; i < N; ++i)
    for (int j = 0; j < N; ++j) {
      tmp_kmatrixDD_norm[i][j] = 0.;
      tmp_kmatrixDD[i][j] = 0.;
    }

  for (int i = 0; i < N; ++i) {
    double ey = pow(10, xos[i]);//reco
    for (int j = 0; j < N; ++j) {
      double ex = pow(10, xos[j]);//real
      if (fabs(i - j) < diagonalOfMatrix)
        for (Double_t h = 0; h < 2.*TMath::Pi(); h += deltah) { //integrate over h
          for (Double_t delta = deltaMin; delta < deltaMax; delta += deltadelta) { // integrate over delta
            Double_t cosTheta = sin(lat) * sin(delta) + cos(lat) * cos(delta) * cos(h);

            double resolution = ResolutionValue(xos[j], cosTheta);
            double sh_shfluct = ShowerToShowerFluctuation(xos[j]);
            double tot_sigma = sqrt(pow(resolution, 2) + pow(sh_shfluct, 2));
            double efficiency = EfficiencyValues(xos[j], acos(cosTheta));
            double bias = EnergyBias_DD(xos[j], cosTheta, pCal);

            if (acos(cosTheta) < mapParameters["MaxZenith"] * TMath::DegToRad() && acos(cosTheta) > 0.) {
              tmp_kmatrixDD_norm[i][j] += deltah * deltadelta *  cos (delta) * cosTheta;
              tmp_kmatrixDD[i][j]      += deltah * deltadelta *  cos (delta) * cosTheta * efficiency  * TMath::Gaus(ey, ex * bias    , tot_sigma * ex, true) * ey * log(10.);
            }
          }
        }
      if (tmp_kmatrixDD_norm[i][j] > 0)
        tmp_kmatrixDD[i][j] /= tmp_kmatrixDD_norm[i][j];
    }
    if (i % 25 == 0)
      cout << i << endl;
  }
  /*
  ofstream outMatrixFile("matrixFile_FULL_NEW.txt");
  for (int i = 0; i < N; ++i)
   {
    for (int j = 0; j < N; ++j)
     outMatrixFile<<tmp_kmatrixDD[i][j]<<" ";

    outMatrixFile<<endl;
   }
  */
  return tmp_kmatrixDD;
}


// Zenith dependent migration matrix with J
// To produce a compact version of the matrix (Originally developed for the PRD paper - 2020)
void kResolutionMatrixZenith_J(TVectorD pCal, TVectorD xos, double* resolutionOffset,TF1 *FitFunctionFitted){
    
  double binSize = mapParameters["SpectrumBinSize"];//Bins in energy
  int N = (mapParameters["MaxIntegration"] - mapParameters["MinIntegration"]) / binSize ;//Bins in energy
  int binnumber = 10;//sub-bins in energy (Simpson calculation on the single energy bin)

  Double_t deltaTheta = 0.1;//bins in theta
  Double_t maxtheta = mapParameters["MaxZenith"];
  Double_t minTheta = mapParameters["MinZenith"];
  Int_t numBinTot = (maxtheta - minTheta) / deltaTheta;
  
  TMatrixD tmp_kmatrixTH(N, N);
  TMatrixD tmp_kmatrixTH_norm(N, N);
  for (int i = 0; i < N; ++i)
    for (int j = 0; j < N; ++j) {
      tmp_kmatrixTH_norm[i][j] = 0.;
      tmp_kmatrixTH[i][j] = 0.;
    }

  for (int i = 0; i < N; ++i) {
    double XoSi = mapParameters["MinIntegration"] + 0.5 * binSize  + binSize * i;//reco
    for (int j = 0; j < N; ++j) {
      double XoSj = mapParameters["MinIntegration"] + 0.5 * binSize + binSize * j;//real

      double norm =0.;
      double normE =0.;
      double ey_low = pow(10, XoSi-0.5 * binSize);
      double ey_high = pow(10, XoSi+0.5 * binSize);
      double dE_integr = (ey_high - ey_low)/binnumber;

      double ey_lowT = pow(10, XoSj-0.5 * binSize);
      double ey_highT = pow(10, XoSj+0.5 * binSize);
      double dE_integrT = (ey_highT - ey_lowT)/binnumber;

      double integralE=0;
      for(Int_t ww=0;ww<=binnumber;ww++)
        {
          double E_tmp = ey_low + ww * dE_integr;

          double integralEreal=0;
          normE =0.;
          for(Int_t w=0;w<=binnumber;w++)
            {
              double E_tmpT = ey_lowT + w * dE_integrT;
              double logEreal = log10(E_tmpT);
              double SpectrumJ = FitFunctionFitted->Eval(logEreal);
              norm =0.;
              double matrixTH =0.;
              Int_t countTh=0;   
                for (Double_t th = minTheta; th <= maxtheta; th += deltaTheta)
                  {
                    Double_t th_rad=th*TMath::DegToRad();
                    Double_t cosTheta = cos(th_rad);
                    double resolution = ResolutionValue(logEreal, cosTheta);
                    double sh_shfluct = ShowerToShowerFluctuation(logEreal);
                    double tot_sigma = sqrt(pow(resolution, 2) + pow(sh_shfluct, 2));
                    double efficiency = EfficiencyValues(logEreal, th_rad);
                    double bias = EnergyBias_DD(logEreal, cosTheta, pCal);
                    if(countTh==0 || countTh==numBinTot)
                      {
                        norm += sin(2.* th_rad) * 2. * TMath::Pi()* SpectrumJ;
                        matrixTH     += sin(2.* th_rad) * 2. * TMath::Pi() * efficiency * TMath::Gaus(E_tmp, E_tmpT * bias    , tot_sigma * E_tmpT, true) * SpectrumJ ;
                      }
                    else
                      {
                        if(countTh%2)
                          {
                            norm += 4 * sin(2.* th_rad) * 2. * TMath::Pi() * SpectrumJ;
                            matrixTH    += 4 * sin(2.* th_rad) * 2. * TMath::Pi() * efficiency * TMath::Gaus(E_tmp, E_tmpT * bias    , tot_sigma * E_tmpT, true) * SpectrumJ;
                          }
                        else
                          {
                            norm += 2 * sin(2.* th_rad) * 2. * TMath::Pi() * SpectrumJ;
                            matrixTH     += 2 * sin(2.* th_rad) * 2. * TMath::Pi() * efficiency * TMath::Gaus(E_tmp, E_tmpT * bias   , tot_sigma * E_tmpT, true) * SpectrumJ;
                          }
                      }  
                    countTh++;
                  }
              norm*=deltaTheta/3.;
              matrixTH*=deltaTheta/3.;

              if(w==0 || w==binnumber)
                {
                  integralEreal+=matrixTH;
                  normE+=norm;
                }
              else
                {
                  if(w%2)
                    {
                      integralEreal+=4.*matrixTH;
                      normE+=4*norm;
                    }
                  else
                    {
                      integralEreal+=2.*matrixTH;
                      normE+=2*norm;
                    }
                }

            }
          integralEreal*=dE_integrT/3.;
          normE*=dE_integrT/3.;


          if(ww==0 || ww==binnumber)
            integralE+=integralEreal;
          else
            {
              if(ww%2)
                integralE+=4.*integralEreal;
              else
                integralE+=2.*integralEreal;
            }     

        }
      integralE*=dE_integr/3.;
      
      if (normE > 0)
         integralE/=normE;

      tmp_kmatrixTH[i][j] = integralE;

      }
  }
  ofstream outMatrixFile("matrixFile_FULL_TH_J.txt");
  for (int i = 0; i < N; ++i)
   {
    for (int j = 0; j < N; ++j)
     outMatrixFile<<tmp_kmatrixTH[i][j]<<" ";

    outMatrixFile<<endl;
   }


}

Double_t DirectionalExposure(Double_t lat, Double_t deltaMin, Double_t deltaMax) {
  Double_t deltah = 0.001;
  Double_t deltadelta = 0.00001;//To make the exposure calculation more precise
  Double_t normalization=0;
  const double tilt = -30 * TMath::DegToRad();
  int countDelta=0;
  Double_t exposure = 0.;

  for (Double_t delta = mapParameters["ExposureMinRangeDecl"]*TMath::DegToRad(); delta < mapParameters["ExposureMaxRangeDecl"]*TMath::DegToRad(); delta += 0.0001) {
    for (Double_t h = 0; h < 2.*TMath::Pi(); h += 0.001) { //integrate over h
      Double_t cosTheta = sin(lat) * sin(delta) + cos(lat) * cos(delta) * cos(h);
      Double_t theta = acos(cosTheta);
      if (acos(cosTheta) < mapParameters["MaxZenith"]* TMath::DegToRad() && acos(cosTheta) > 0.) {
        double phi = atan2(sin(delta) * cos(lat) - cos(delta) * cos(h) * sin(lat), cos(delta) * sin(h));
        normalization += 0.001 * 0.0001 * cos(delta) * cosTheta * (1 + 0.00348 * tan(theta) * cos(phi - tilt));
        if(countDelta%10000000==0)
          cout<<"Normalization calculation step "<<countDelta<<endl;

        countDelta++;
        
      }
    }
  }
  cout<<"Normalization exposure integral "<<  normalization<<endl;

  countDelta=0;
  for (Double_t delta = deltaMin; delta < deltaMax; delta += deltadelta) {
    for (Double_t h = 0; h < 2.*TMath::Pi(); h += deltah) { //integrate over h
      Double_t cosTheta = sin(lat) * sin(delta) + cos(lat) * cos(delta) * cos(h);
      Double_t theta = acos(cosTheta);
      if (acos(cosTheta) < mapParameters["MaxZenith"]* TMath::DegToRad() && acos(cosTheta) > 0.) {
        double phi = atan2(sin(delta) * cos(lat) - cos(delta) * cos(h) * sin(lat), cos(delta) * sin(h));
        exposure += deltah * deltadelta * cos(delta) * cosTheta * (1 + 0.00348 * tan(theta) * cos(phi - tilt));
      }

     if(countDelta%10000000==0)
      cout<<"Exposure calculation step "<<countDelta<<endl;

     countDelta++;
    }
  }
 cout << " Exposure " <<  mapParameters["Exposure"] * exposure / normalization << "km2 sr yr" << endl;
 //return mapParameters["Exposure"] * exposure / normalization;
 return mapParameters["Exposure"] * exposure / 2.35718;
}


Double_t DirectionalExposureZenith(Double_t thetaMin, Double_t thetaMax) {
  Double_t deltatheta = 0.0001;
  Double_t normalization=0;
  for (Double_t theta = mapParameters["ExposureMinRangeZenith"]; theta < mapParameters["ExposureMaxRangeZenith"]; theta += deltatheta) {
    normalization+=sin(2.*theta*TMath::DegToRad())*deltatheta * TMath::DegToRad();
  }
  cout<<"Normalization exposure integral "<<  normalization<<endl;
                                              
  Double_t exposure = 0.;
  for (Double_t theta = thetaMin; theta < thetaMax; theta += deltatheta) {
    exposure+=sin(2.*theta*TMath::DegToRad())*deltatheta * TMath::DegToRad();
  }

  cout << " Exposure " <<  mapParameters["Exposure"] * exposure / normalization << "km2 sr yr" << endl;
//  return mapParameters["Exposure"] * exposure / normalization;
  return mapParameters["Exposure"] * exposure / 0.75;
}


void PrintMapParameters() {
  cerr << "\n---------------------------------------------------\n";
  cerr << "------------------Parameter Map--------------------\n";
  cerr << "---------------------------------------------------\n";

  std::map<string, Double_t>::iterator it = mapParameters.begin();
  while (it != mapParameters.end()) {
    std::string keyName = it->first;
    double val = it->second;
    cout << keyName << " :: " << val << endl;
    it++;
  }

  cout << "Model---------------------" << Model << endl;
  cout << "MatrixType---" << MatrixType << endl;
  cout << "ResolutionFlag------------" << ResolutionFlag << endl;
  cout << "BiasFlag------------------" << BiasFlag << endl;
  cout << "EfficiencyFlag------------" << EfficiencyFlag << endl;
  cout << "ShowerToShowerFlag--------" << ShowerToShowerFlag << endl;
  cout << "VersionInput--------------" << VersionInput << endl;
  cout << "UseNewTriggers------------" << UseNewTriggers << endl;
  cerr << "---------------------------------------------------\n";
  cerr << "---------------------------------------------------\n";
}

void LoadParam(string paramFileName) {
  ifstream in;
  in.open(paramFileName.c_str());
  string paramString, line;

  while (!in.eof()) {
    getline(in, line);
    Int_t pos = line.find('#');
    string lineNew = line.substr(0, pos);
    if (lineNew.length()) {
      Int_t pos2 = lineNew.find(' ');
      string key = lineNew.substr(0, pos2);

      string lineNew2 = lineNew.substr(pos2, 1000);
      Int_t posVal = lineNew2.find_first_not_of(' ');
      string lineNew3 = lineNew2.substr(posVal, 1000);
      Int_t posVal2 = lineNew3.find(' ');
      string val = lineNew3.substr(0, posVal2);

      Double_t tmpVal;
      if (key == "Code") { //if we read a string... Define a string variable
        CodeUsed = val;
      } else if (key == "DataFileName") {
        SDFile = val;
      } else if (key == "OutPutFileName") {
        OutPutFileName = val;
      } else if (key == "OutPutFileNameRebinned") {
        OutPutFileNameRebinned = val;
      } else if (key == "ModelFunction") {
        Model = val;
      } else if (key == "MatrixType") {
        MatrixType = val;
      } else if (key == "ResolutionFlag") {
        ResolutionFlag = val;
      } else if (key == "BiasFlag") {
        BiasFlag = val;
      } else if (key == "EfficiencyFlag") {
        EfficiencyFlag = val;
      } else if (key == "ShowerToShowerFlag") {
        ShowerToShowerFlag = val;
      } else if (key == "VersionInput") {
        VersionInput = val;
      } else if (key == "UseNewTriggers") {
        UseNewTriggers = val;
      } else if (key == "PrintCompact") {
        PrintCompact= val;
      } else if (key == "Verbose") {
        Verbose = val;
      } else { // if we read a number do nothing... Only add parameter into the config file.
        istringstream stringstream(val.c_str());
        stringstream >> tmpVal;
        mapParameters.insert(make_pair(key.c_str(), tmpVal));
      }
    }
  }

  if ("no" != Verbose)
    PrintMapParameters();

}


void ReadEventList(vector <double> &vLogESD, TH1D* th1, TH1D* th2, string FileName) {

  long int AugerId;
  long int EventId;
  int nst;
  double lgE, S1000, dS1000, theta, phi, dec, ra;
  Double_t A, B;
  string labelA, labelB;
  ifstream evtFile;

  vLogESD.clear();

  // test format. Minimal check: verify that the number of columns is compatible with expectations...
  // The user is responsible for the content of the event file!
  if ("no" != Verbose)
    TestFormat(FileName, 0);


  evtFile.open(FileName.c_str(), ios::in);
  if (evtFile.fail()) {
    cerr << endl;
    cerr << FileName << ": no valid file name  chosen " << endl;
    throw std::exception();
  }

  string headerline;
  getline(evtFile, headerline);
  headerline = headerline.substr(0, 1);
  evtFile >> labelA >> A;
  evtFile >> labelB >> B;
  if (headerline != "#" || labelA != "A" || labelB != "B") {
    cerr << endl;
    cerr << FileName << " does not look to be formatted properly " << endl;
    cerr << "Add header line and calibration parameters" << endl;
    throw std::exception();
  }
  mapParameters.insert(make_pair(labelA.c_str(), A));
  mapParameters.insert(make_pair(labelB.c_str(), B));

  Double_t maxDecl = mapParameters["MaxDeclination"];
  Double_t minDecl = mapParameters["MinDeclination"];
  Double_t maxZenith = mapParameters["MaxZenith"];
  Double_t minZenith = mapParameters["MinZenith"];

  if (MatrixType == "Standard")
    { //TODO
      maxDecl = 9000.;
      minDecl = -9000;
      maxZenith = 9000.;
      minZenith = -9000;
    }
  else if(MatrixType == "Zenith")
    {
      maxDecl = 9000.;
      minDecl = -9000;
    }
  else if(MatrixType == "Declination")
    {
      maxZenith = 9000.;
      minZenith = -9000;
    }

  while (!evtFile.eof()) {

    evtFile >> AugerId >> EventId >> lgE >> S1000 >> dS1000 >> nst >> theta >> phi >> dec >> ra;

    if (HeraldAnalysis)
      lgE = log10(lgE) + 18;

    if (!evtFile.eof()) {
      if (dec < maxDecl && dec >= minDecl && theta<maxZenith && theta>=minZenith) {
        th1->Fill(lgE);
        th2->Fill(lgE);
        if ( lgE >= th1->GetBinLowEdge(1) )
          vLogESD.push_back(lgE);
      }
    }
  }

  evtFile.close();
}


void ReadEventListICRC2019(vector <double> &vLogESD, TH1D* th1, TH1D* th2, string FileName) {

  long int AugerId;
  long int EventId;
  long int GPSTime;
  int isHybrid;
  double lgE, S1000, dS1000, theta, dtheta, phi, dphi, dec, ra, XCore, YCore, S1000_corr, S38, Nst, Nst_sat;
  Double_t A, B;
  string labelA, labelB;
  ifstream evtFile;

  vLogESD.clear();

  // test format. Minimal check: verify that the number of columns is compatible with expectations...
  // The user is responsible for the content of the event file!
  if ("no" != Verbose)
    TestFormat(FileName, 1);

  evtFile.open(FileName.c_str(), ios::in);
  if (evtFile.fail()) {
    cerr << endl;
    cerr << FileName << ": no valid file name  chosen " << endl;
    throw std::exception();
  }
  string headerline;
  getline(evtFile, headerline);
  headerline = headerline.substr(0, 1);
  evtFile >> labelA >> A;
  evtFile >> labelB >> B;
  if (headerline != "#" || labelA != "A" || labelB != "B") {
    cerr << endl;
    cerr << FileName << " does not look to be formatted properly " << endl;
    cerr << "Add header line and calibration parameters" << endl;
    throw std::exception();
  }

  mapParameters.insert(make_pair(labelA.c_str(), A));
  mapParameters.insert(make_pair(labelB.c_str(), B));

  Double_t maxDecl = mapParameters["MaxDeclination"];
  Double_t minDecl = mapParameters["MinDeclination"];
  Double_t maxZenith = mapParameters["MaxZenith"];
  Double_t minZenith = mapParameters["MinZenith"];

  if (MatrixType == "Standard")
    { //TODO
      maxDecl = 9000.;
      minDecl = -9000;
      maxZenith = 9000.;
      minZenith = -9000;
    }
  else if(MatrixType == "Zenith")
    {
      maxDecl = 9000.;
      minDecl = -9000;
    }
  else if(MatrixType == "Declination")
    {
      maxZenith = 9000.;
      minZenith = -9000;
    }

  while (!evtFile.eof()) {

    evtFile >> AugerId >> EventId >> GPSTime >> isHybrid >> Nst >> Nst_sat >> lgE >> theta >> dtheta >> phi >> dphi >> XCore >> YCore >> ra >> dec >> S1000 >> dS1000 >> S1000_corr >> S38;
//   evtFile >> AugerId >> EventId >> GPSTime>> isHybrid >> lgE >> theta >> dtheta >> phi >> dphi >> XCore >> YCore >> ra>> dec >> S1000 >> dS1000 >> S1000_corr >> S38;//PRELIMINARY VERSION FOR ICRC2019

    if (!evtFile.eof()) {
      if (dec < maxDecl && dec >= minDecl  && theta<maxZenith && theta>=minZenith) {
        th1->Fill(lgE);
        th2->Fill(lgE);
        if ( lgE >= th1->GetBinLowEdge(1) )
          vLogESD.push_back(lgE);
      }
    }
  }
  evtFile.close();
}


void TestFormat(string inputName, int flagversion) {
  ifstream evtFile;
  string headerline;
  evtFile.open(inputName.c_str(), ios::in);
  stringstream testString;
  getline(evtFile, headerline);
  getline(evtFile, headerline);
  getline(evtFile, headerline);
  getline(evtFile, headerline);
  testString << headerline;
  double x;
  int countCol = 0;
  cout << endl;

  while (testString >> x) {
    cout << x << " ";
    countCol++;
  }
  cout << endl;
  cout << countCol << " columns in the input event file" << endl;

  cout << "Check consistency with... " << endl;
  if (flagversion == 1)
    cout << "evtFile >> AugerId >> EventId >> GPSTime>> isHybrid>> Nst>> Nst_sat >> lgE >> theta >> dtheta >> phi >> dphi >> XCore >> YCore >> ra>> dec >> S1000 >> dS1000 >> S1000_corr >> S38;" << endl;
  if (flagversion == 0)
    cout << "evtFile >> AugerId >> EventId >> lgE >> S1000 >> dS1000 >> nst >> theta >> phi >> dec >> ra;" << endl;
  evtFile.close();
}


void PrintUsage(int nArgs) {
  cerr << "The correct number of inputs were not found " << endl;
  cerr << "Expected 2 and got " << nArgs << endl;
  cerr << "Run this program as ./UnfoldSpectrum <parameter file>" << endl;
}