ICRC2019/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 DeclinationDepentMatrix;
string ResolutionFlag;
string BiasFlag;
string EfficiencyFlag;
string ShowerToShowerFlag;
string VersionInput = "standard";
string UseNewTriggers = "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;
}

// 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..
  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;
}

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 deltadelta = 0.001;
  const double tilt = -30 * TMath::DegToRad();
  int countDelta=0;
  Double_t exposure = 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 / 2.35718 << "km2 sr yr" << endl;
 return mapParameters["Exposure"] * exposure / 2.35718;
}

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 << "DeclinationDepentMatrix---" << DeclinationDepentMatrix << 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 == "DeclinationDepentMatrix") {
        DeclinationDepentMatrix = 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 == "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(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;

  // 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"];
  if (DeclinationDepentMatrix == "no") { //TODO
    maxDecl = 9000.;
    minDecl = -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) {
        th1->Fill(lgE);
        th2->Fill(lgE);
      }
    }
  }

  evtFile.close();
}


void ReadEventListICRC2019(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;

  // 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"];
  if (DeclinationDepentMatrix == "no") {
    maxDecl = 90.;
    minDecl = -90;
  }

  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 <30.26)
        th1->Fill(lgE);
        th2->Fill(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;
}