DetectorResponse.cc

Go to the documentation of this file.

#include "DetectorResponse.h"
#include "DetectorResponseConfig.h"

#include <sstream>

#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <algorithm>
#include <stdexcept>


using namespace std;
using namespace TabulatedTankSimulatorNS;

namespace TabulatedTankSimulatorNS {

  static const char* electronName   = "electron";
  static const char* gammaName      = "gamma";
  static const char* muonName       = "muon";
  static const char* muonNameLow   = "muon-low";              // The area over which particles are injected is larger
  static const char* michelName     = "michel-electron";
  static const char* michelNameLow = "michel-electron-low";   // The area over which particles are injected is larger

  //----- Zenith grid  [0] WCD [1] ASCII [2] AMIGA
  const unsigned int nTo_Grid[3] = { 19 , 19 , 22  };
  const double Theta_Grid[3][50] = {
    { 0, 12, 25, 36, 45, 53, 60, 63, 67, 70, 72, 74, 76, 78, 80, 82, 84, 86, 88},
    { 0, 12, 25, 36, 45, 53, 60, 63, 67, 70, 72, 74, 76, 78, 80, 82, 84, 86, 88},
    { 0, 8, 12, 20, 25, 27, 30, 36, 38, 40, 42, 45, 47, 50, 53, 57, 60, 62, 64, 66, 68, 70 }
  };

  //------Energy grid  [0] WCD [1] ASCII [2] AMIGA

  const unsigned int  nKE_Electron[3] = { 35 , 35 , 35 };
  const double lKE_Electron[3][50] = {
    {
      -4.000, -3.800, -3.400, -3.200, -3.000, -2.900, -2.800, -2.700, -2.600, -2.500,
      -2.400, -2.300, -2.200, -2.100, -2.000, -1.800, -1.600, -1.400, -1.200, -1.000,
      -0.800, -0.600, -0.400, -0.200,  0.000,  0.200,  0.400,  0.600,  0.800,  1.000,
      1.200,  1.400,  1.600,  1.800,  2.000
    },

    {
      -4.000, -3.800, -3.400, -3.200, -3.000, -2.900, -2.800, -2.700, -2.600, -2.500,
      -2.400, -2.300, -2.200, -2.100, -2.000, -1.800, -1.600, -1.400, -1.200, -1.000,
      -0.800, -0.600, -0.400, -0.200,  0.000,  0.200,  0.400,  0.600,  0.800,  1.000,
      1.200,  1.400,  1.600,  1.800,  2.000
    },

    {
      -4.000, -3.800, -3.400, -3.200, -3.000, -2.900, -2.800, -2.700, -2.600, -2.500,
      -2.400, -2.300, -2.200, -2.100, -2.000, -1.800, -1.600, -1.400, -1.200, -1.000,
      -0.800, -0.600, -0.400, -0.200,  0.000,  0.200,  0.400,  0.600,  0.800,  1.000,
      1.200,  1.400,  1.600,  1.800,  2.000
    }
  };


  const unsigned int  nKE_Gamma[3] = { 31 , 31 , 31 };
  const double lKE_Gamma[3][50] = {
    {
      -4.000, -3.800, -3.600, -3.400, -3.200, -3.000, -2.800, -2.600, -2.400, -2.200,
      -2.000, -1.800, -1.600, -1.400, -1.200, -1.000, -0.800, -0.600, -0.400, -0.200,
      0.000,  0.200,  0.400,  0.600,  0.800,  1.000,  1.200,  1.400,  1.600, 1.800,
      2.000
    },

    {
      -4.000, -3.800, -3.600, -3.400, -3.200, -3.000, -2.800, -2.600, -2.400, -2.200,
      -2.000, -1.800, -1.600, -1.400, -1.200, -1.000, -0.800, -0.600, -0.400, -0.200,
      0.000,  0.200,  0.400,  0.600,  0.800,  1.000,  1.200,  1.400,  1.600, 1.800,
      2.000
    },


    {
      -4.000, -3.800, -3.600, -3.400, -3.200, -3.000, -2.800, -2.600, -2.400, -2.200,
      -2.000, -1.800, -1.600, -1.400, -1.200, -1.000, -0.800, -0.600, -0.400, -0.200,
      0.000,  0.200,  0.400,  0.600,  0.800,  1.000,  1.200,  1.400,  1.600, 1.800,
      2.000
    }
  };



  const unsigned int nKE_Muon[3] = { 34 , 34 , 34 };
  const double lKE_Muon[3][50] = {

    {
      -3.000, -2.800, -2.600, -2.400, -2.200, -2.000, -1.850, -1.700, -1.550, -1.400,
      -1.250, -1.100, -0.950, -0.800, -0.650, -0.500, -0.350, -0.200, -0.050,  0.100,
      0.250,  0.400,  0.550,  0.700,  0.850,  1.000,  1.300,  1.600,  1.900,  2.200,
      2.500,  2.800,  3.100, 3.400
    },

    {
      -3.000, -2.800, -2.600, -2.400, -2.200, -2.000, -1.850, -1.700, -1.550, -1.400,
      -1.250, -1.100, -0.950, -0.800, -0.650, -0.500, -0.350, -0.200, -0.050,  0.100,
      0.250,  0.400,  0.550,  0.700,  0.850,  1.000,  1.300,  1.600,  1.900,  2.200,
      2.500,  2.800,  3.100, 3.400
    },

    {
      -0.600, -0.500, -0.400, -0.300, -0.200, -0.150, -0.100, -0.075, -0.050, -0.025,
      0.000, 0.025,  0.050,  0.075,  0.100,  0.125,  0.150,  0.175,  0.200,  0.225,
      0.250, 0.275,  0.300,  0.325,  0.350,  0.375,  0.400,  0.425,  0.500,  0.600,
      0.800, 1.000, 1.200, 1.400
    }
  };


  const unsigned int nKE_MuonLow[3] = { 19 , 19 , 34 };
  const double lKE_MuonLow[3][50] = {

    {
      -3.000, -2.800, -2.600, -2.400, -2.200, -2.000, -1.850, -1.700, -1.550, -1.400,
      -1.250, -1.100, -0.950, -0.800, -0.650, -0.500, -0.350, -0.200, -0.050
    },

    {
      -3.000, -2.800, -2.600, -2.400, -2.200, -2.000, -1.850, -1.700, -1.550, -1.400,
      -1.250, -1.100, -0.950, -0.800, -0.650, -0.500, -0.350, -0.200, -0.050
    },

    {
      -0.600, -0.500, -0.400, -0.300, -0.200, -0.150, -0.100, -0.075, -0.050, -0.025,
      0.000, 0.025,  0.050,  0.075,  0.100,  0.125,  0.150,  0.175,  0.200,  0.225,
      0.250, 0.275,  0.300,  0.325,  0.350,  0.375,  0.400,  0.425,  0.500,  0.600,
      0.800, 1.000, 1.200, 1.400
    }
  };



  inline bool peordering(const pair<int, double>& v1, const pair<int, double>& v2)
  {
    return (v1.first < v2.first);
  }

}


double
DetectorResponse::GetAreaProj(int itype, double theta, double ke)
{
  // NOTE : For ASCII the projected area has a modulation with azimuth.
  //        At 88 [deg] there is a modulation of +-7% with azimuth.
  //        The mean over all azimuths is used, which corresponds approximately to the case of azimuth=90 [deg]
  //        ( projected particle direction perpendicular to the long side of the box)

  if (itype < 0 || itype > 2) return 0.;

  double f = 1.;
  if (fUseEnlargeArea && itype == 0 && log10(ke) < lKE_MuonLow[fDetectorType][nKE_MuonLow[fDetectorType] - 1]) f = 2.;

  if (fDetectorType == 0)         return M_PI * tankr * tankr * f * f + 2.*tankh * tankr * f * tan(theta);
  else if (fDetectorType == 1) {
    double L = scinlength * f, W = scinwidth * f, H = scinheight;
    double /*Max=sqrt(L*L+W*W)*H*tan(theta), Min=W*H*tan(theta),*/ Mean = L * H * tan(theta);
    return  L * W + Mean;
  } else                          return MDlength * MDwidth * nMD;
}

double
DetectorResponse::GetAreaPerp(int itype, double theta, double ke)
{
  return GetAreaProj(itype, theta, ke) * cos(theta);
}

double
DetectorResponse::GetRelativeTrack(double theta)
{
  if (fDetectorType == 0) return 1. / (cos(theta) + 2.*tankh * sin(theta) / M_PI / tankr);
  else                  return 1. / cos(theta);
}


//-----------------------------------------------------------------------------
//------------Constructor
//-----------------------------------------------------------------------------



DetectorResponse::DetectorResponse(bool flag, int DetectorType, bool UseEnlargeArea)
{
  fUseEnlargeArea = UseEnlargeArea;
  fDetectorType = DetectorType;

  //Check valid detector
  if (fDetectorType < 0 || fDetectorType > 2) {
    printf("Incorrect detector type %d \n", fDetectorType);
    exit(1);
  }

  //Set iGrid
  iGrid = fDetectorType;

  //Set NPECUT
  if (fDetectorType == 2) NPECUT = UNDEF;
  else                    NPECUT = 10;

  //Set CDF/G4 Dir
  fDataDir = DATA_DIR;
  fDataDir += "/";

  if (fDetectorType == 0 && flag) fDataDir += "g4files-wcd";
  else if (fDetectorType == 0 && !flag) fDataDir += "cdffiles-wcd";
  else if (fDetectorType == 1 && flag) fDataDir += "g4files-scin";
  else if (fDetectorType == 1 && !flag) fDataDir += "cdffiles-scin";
  else if (fDetectorType == 2 && flag) fDataDir += "g4files-md";
  else if (fDetectorType == 2 && !flag) fDataDir += "cdffiles-md";


  for (int itype = 0; itype < NTYPES - 2; itype++) { // NTYPES particle species implemented so far
    //const char *pName;
    unsigned int nx = nTo_Grid[iGrid], ny;
    const double* xvalue = Theta_Grid[iGrid], *yvalue;
    vector <PeCDF>* fCDFs, *fCDFs_dec;
    switch (itype) {
      case 0: {
        /*pName=muonName;*/     ny = nKE_Muon[iGrid] ;
        yvalue = lKE_Muon[iGrid];
        fCDFs = &fMuonCDFs;
        fCDFs_dec = &fMichelElectronCDFs;
        break;
      }
      case 1: {
        /*pName=electronName;*/ ny = nKE_Electron[iGrid];
        yvalue = lKE_Electron[iGrid];
        fCDFs = &fElectronCDFs;
        fCDFs_dec = NULL;
        break;
      }
      case 2: {
        /*pName=gammaName;*/    ny = nKE_Gamma[iGrid];
        yvalue = lKE_Gamma[iGrid];
        fCDFs = &fGammaCDFs;
        fCDFs_dec = NULL;
        break;
      }
      case 3: {
        /*pName=muonNameLow;*/  ny = nKE_MuonLow[iGrid];
        yvalue = lKE_MuonLow[iGrid];
        fCDFs = &fMuonLowCDFs;
        fCDFs_dec = &fMichelElectronLowCDFs;
        break;
      }
      default:
        continue;
    }

    if (fDetectorType == 2 && itype != 0) continue;

    //--Create the CDFs
    for (unsigned int j = 0; j < ny; j++) {
      for (unsigned int i = 0; i < nx; i++) {
        PeCDF thisCDF, thisCDF_MichelElectron;
        thisCDF.itype = itype;

        thisCDF.theta = xvalue[i] * M_PI / 180.;
        thisCDF.theta_deg = xvalue[i];
        thisCDF.lke = yvalue[j];

        thisCDF.mean = 0.;
        thisCDF.rms = 0.;
        thisCDF.nentries = 0;
        thisCDF.nentriesCut = 0;
        thisCDF.probZero = 0.;
        thisCDF.probCut = 0.;
        thisCDF.probDecay = 0.;
        thisCDF.probHit = 0.;
        thisCDF.PeCut = NPECUT;
        thisCDF.median = 0.;

        thisCDF_MichelElectron = thisCDF;
        thisCDF_MichelElectron.itype = (itype == 0 ? 4 : (itype == 3 ? 5 : UNDEF));

        if (flag) ReadG4File(thisCDF, thisCDF_MichelElectron);

        fCDFs->push_back(thisCDF);
        //---Signal from Michel electron
        if (fCDFs_dec != NULL && DetectorType != 2) fCDFs_dec->push_back(thisCDF_MichelElectron);
      }
    }
    bool ok = true;
    if (!flag) {
      ok = ReadCDFFile(fCDFs, nx * ny);
      if (fCDFs_dec != NULL && DetectorType != 2) ReadCDFFile(fCDFs_dec, nx * ny);
    } else {
      ok = WriteCDFFile(fCDFs);
      if (fCDFs_dec != NULL && DetectorType != 2) WriteCDFFile(fCDFs_dec);
    }

    if (!ok) {
      printf("Severe problem initializing DetectorResponse | Failed reading/writing CDF Files  \n");
      exit(1);
    }
  }


  // Set Unit per PE
  if (fDetectorType == 0) UnitPerPE = UnitPerPE_wcd;
  else if (fDetectorType == 1) UnitPerPE = UnitPerPE_scin;
  else                         UnitPerPE = UnitPerPE_md;
}


//-----------------------------------------------------------------------------
//------------Get Name, Get bin number, Get CDF of a bin
//-----------------------------------------------------------------------------

const char* DetectorResponse::GetParticleName(int itype)
{
  if (itype == 0) return muonName;
  else if (itype == 1) return electronName;
  else if (itype == 2) return gammaName;
  else if (itype == 3) return muonNameLow;
  else if (itype == 4) return michelName;
  else if (itype == 5) return michelNameLow;

  return "unknownParticle";
}


//lke log10(ke) GeV , itype (internal coding), to in radians
// ibins[4] -->  [0] (to1,ke1) [1] (to1,ke2) [2] (to2,ke1) [3] (to2,ke2)
// cbin --> closest bin
int
DetectorResponse::GetBin(double lke, double to, int itype, int* ibins , bool& flag)
{
  int cbin;
  unsigned int nx = nTo_Grid[iGrid], ny;
  const double* xvalue = Theta_Grid[iGrid], *yvalue;
  switch (itype) {
    case 0: {
      ny = nKE_Muon[iGrid] ;
      yvalue = lKE_Muon[iGrid];
      break;
    }
    case 1: {
      ny = nKE_Electron[iGrid];
      yvalue = lKE_Electron[iGrid];
      break;
    }
    case 2: {
      ny = nKE_Gamma[iGrid];
      yvalue = lKE_Gamma[iGrid];
      break;
    }
    case 3: {
      ny = nKE_MuonLow[iGrid];
      yvalue = lKE_MuonLow[iGrid];
      break;
    }
    case 4: {
      ny = nKE_Muon[iGrid] ;
      yvalue = lKE_Muon[iGrid];
      break;
    }
    case 5: {
      ny = nKE_MuonLow[iGrid] ;
      yvalue = lKE_MuonLow[iGrid];
      break;
    }
    default:
      ostringstream what;
      what << "Unknown particle type (" << itype << ") in " << __func__ << " at " << __FILE__ << ":" << __LINE__;
      throw std::logic_error(what.str());
  }

  double lkemin = yvalue[0], lkemax = yvalue[ny - 1];
  double thetamin = xvalue[0] * M_PI / 180., thetamax = xvalue[nx - 1] * M_PI / 180.;
  int ix_bin = -1, iy_bin = -1;
  bool extrapol_flag_to = false;
  bool extrapol_flag_ke = false;
  flag = false;

  //Extrapolation check
  if (to < thetamin) {
    ix_bin = 0;
    extrapol_flag_to = true;
    flag = true;
  } else if (to >= thetamax) {
    ix_bin = nx - 2;
    extrapol_flag_to = true;
    flag = true;
  }

  if (lke < lkemin) {
    iy_bin = 0;
    extrapol_flag_ke = true;
    flag = true;
  } else if (lke >= lkemax) {
    iy_bin = ny - 2;
    extrapol_flag_ke = true;
    flag = true;
  }

  //Interpolation bins
  if (extrapol_flag_to) {
    if (DEBUG_RESP == 1)
      printf("Warning: Particle angle outside range, an extrapolation will be used, %10.2lf  deg !\n", to * 180. / M_PI);
  } else {
    for (unsigned int ix = 0; ix < (nx - 1); ix++)
      if (to >= (xvalue[ix]*M_PI / 180.) && to < (xvalue[ix + 1]*M_PI / 180.)) {
        ix_bin = ix;
        break;
      }
  }

  if (extrapol_flag_ke) {
    if (DEBUG_RESP == 1) printf("Warning: Particle Energy outside range, an extrapolation will be used, log(ke)=%10.7lf GeV \n", lke);
  } else {
    for (unsigned int iy = 0; iy < (ny - 1); iy++)
      if (lke >=  yvalue[iy]  && lke < yvalue[iy + 1]) {
        iy_bin = iy;
        break;
      }
  }

  if (nx == 1) {
    ibins[0] = iy_bin * nx;
    ibins[1] = (iy_bin + 1) * nx;
    ibins[2] = iy_bin * nx;
    ibins[3] = (iy_bin + 1) * nx;
  } else {
    ibins[0] = iy_bin * nx + ix_bin;
    ibins[1] = (iy_bin + 1) * nx + ix_bin;
    ibins[2] = iy_bin * nx + ix_bin + 1;
    ibins[3] = (iy_bin + 1) * nx + ix_bin + 1;
  }

  //Closest bin
  double min;
  int cbin_x = -1, cbin_y = -1;
  min = 1.e10;
  for (unsigned int ix = 0; ix < nx; ix++)
    if (fabs(to * 180. / M_PI - xvalue[ix]) < min) {
      cbin_x = ix;
      min = fabs(to * 180. / M_PI - xvalue[ix]);
    }
  min = 1.e10;
  for (unsigned int iy = 0; iy < ny; iy++)
    if (fabs(lke - yvalue[iy]) < min) {
      cbin_y = iy;
      min = fabs(lke - yvalue[iy]);
    }
  if (cbin_x < 0 || cbin_y < 0) return -1;

  cbin = cbin_y * nx + cbin_x;

  return cbin;
}


PeCDF
DetectorResponse::GetCDF(int itype, int ibin)
{
  vector <PeCDF>* fCDFs;
  switch (itype) {
    case 0: {
      fCDFs = &fMuonCDFs;
      break;
    }
    case 1: {
      fCDFs = &fElectronCDFs;
      break;
    }
    case 2: {
      fCDFs = &fGammaCDFs;
      break;
    }
    case 3: {
      fCDFs = &fMuonLowCDFs;
      break;
    }
    case 4: {
      fCDFs = &fMichelElectronCDFs;
      break;
    }
    case 5: {
      fCDFs = &fMichelElectronLowCDFs;
      break;
    }

    default:
      ostringstream what;
      what << "Unknown particle type (" << itype << ") in " << __func__ << " at " << __FILE__ << ":" << __LINE__;
      throw std::logic_error(what.str());
  }

  PeCDF pCDF = fCDFs->at(ibin);

  printf("%s (%d) log(ke)=%4.2lf  theta=%.0lf  mean=%e rms=%e  nentries=%d pzero=%e pcut=%e | (-npe,npeCut)(%d,%d) ",
         GetParticleName(pCDF.itype), pCDF.itype, pCDF.lke, pCDF.theta_deg,
         pCDF.mean, pCDF.rms, pCDF.nentries, pCDF.probZero, pCDF.probCut, pCDF.cdf[0].first, NPECUT);
  if (itype == 0) printf(" pdecay=%e \n", pCDF.probDecay);
  else printf("\n");
  return pCDF;
}



//------------------------------------------------------------
//--Read the G4  File and set up the CDF
//------------------------------------------------------------

bool
DetectorResponse::ReadG4File(PeCDF& pCDF, PeCDF& pCDF_dec)
{

  //G4File name
  char G4FileName[200];
  sprintf(G4FileName, "%s/%s/ke%5.3lf_to%d.dat", fDataDir.c_str()
          , GetParticleName(pCDF.itype), pCDF.lke, int(pCDF.theta_deg));

  //--Open G4 file
  FILE* fp = fopen(G4FileName, "r");
  if (fp == NULL) {
    printf("Could not open %s file \n", G4FileName);
    return false;
  }

  printf("Reading %s file \n", G4FileName);

  while (true) {
    int partId;
    double keo;
    double to = 0;
    double x, y, z;

    //int decayFlag=0;
    double petot = 0., petot_decay = 0.;

    if (fDetectorType < 2) {
      // Length units [mm]
      double ASCII_Track_p, ASCII_Edep_p;
      double ASCII_Track, ASCII_Edep, ASCII_Edep_Michel;

      double ASCII_petot, ASCII_petot_Michel;
      double ASCII_FADCSum_LG, ASCII_FADCSum_HG;

      double WCD_Track_p;
      double WCD_petot, WCD_petot_Michel;
      double WCD_FADCSum_LG, WCD_FADCSum_HG;

      int decayFlag_Hall;
      double MichelKE, MichelZ, MichelTheta;

      double dummy1, dummy2;
      double fTrack_Ref, fTrack_ASCII, fTrack_WCD;

      int nread = fscanf(fp, "%d %lf %lf %lf %lf %lf     %lf %lf %lf %lf   %lf %lf %lf %lf %lf   %lf %lf  %lf %lf %lf %lf  %d %lf %lf %lf   %lf %lf %lf %lf ",
                         &partId, &to, &keo, &x, &y, &z,
                         &fTrack_Ref, &fTrack_ASCII, &fTrack_WCD, &WCD_Track_p,
                         &ASCII_Track_p, &ASCII_Edep_p, &ASCII_Track, &ASCII_Edep, &ASCII_Edep_Michel,
                         &ASCII_petot, &WCD_petot,
                         &ASCII_FADCSum_LG, &ASCII_FADCSum_HG, &WCD_FADCSum_LG, &WCD_FADCSum_HG,
                         &decayFlag_Hall, &MichelKE, &MichelZ, &MichelTheta,
                         &ASCII_petot_Michel, &WCD_petot_Michel , &dummy1, &dummy2);
      if (nread != 29) break;

      petot = (fDetectorType == 0 ? WCD_petot - WCD_petot_Michel : ASCII_petot - ASCII_petot_Michel);
      petot_decay = (fDetectorType == 0 ?      WCD_petot_Michel       :    ASCII_petot_Michel);

      for (int icase = 0; icase < 2; icase++) {
        double pe = (icase == 0 ? petot :  petot_decay);
        PeCDF* thisCDF = (icase == 0 ? &pCDF : &pCDF_dec);

        thisCDF->mean += pe;
        thisCDF->rms += pow(pe, 2.);
        thisCDF->nentries++;

        if (pe  == 0.)                     thisCDF->probZero += 1.;
        if (icase == 0 && petot_decay > 0.)     thisCDF->probDecay += 1.;
        if (pe < NPECUT && NPECUT != UNDEF) {
          thisCDF->probCut += 1.;
          thisCDF->nentriesCut++;
        }

        int size = thisCDF->prob.size();
        bool flag_e = false;
        for (int ipe = 0; ipe < size; ipe++)
          if (int(pe) == thisCDF->prob[ipe].first) {
            flag_e = true;
            thisCDF->prob[ipe].second += 1.;
          }
        if (!flag_e) thisCDF->prob.push_back(pair<int, double>((int)pe, 1));
      }
    }
    //----
    else if (fDetectorType == 2) {
      double /*theta,ke,xinj,yinj,zinj,*/track_p, edep_p, track, edep, meanX, npe;
      //unsigned int id;
      double xabove    , yabove    , zabove;
      double xbelow    , ybelow    ,/*zbelow,*/ebelow, tbelow;
      double xabove_sec, yabove_sec,/*zabove_sec,*/eabove_sec, tabove_sec;
      double xbelow_sec, ybelow_sec,/*zbelow_sec,*/ebelow_sec, tbelow_sec;

      unsigned int nread = fscanf(fp, "%d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf ",
                                  &partId, &to, &keo, &x, &y, &z, &track_p, &edep_p, &track, &edep, &npe, &meanX,
                                  &xabove,    &yabove,    &zabove,
                                  &ebelow,    &xbelow,    &ybelow,    &tbelow,
                                  &eabove_sec, &xabove_sec, &yabove_sec, &tabove_sec,
                                  &ebelow_sec, &xbelow_sec, &ybelow_sec, &tbelow_sec);
      if (nread != 27) break;
      if (xbelow > 0. && xbelow / 10. < TabulatedTankSimulatorNS::MDlength && fabs(ybelow / 10.) < TabulatedTankSimulatorNS::MDwidth / 2.) pCDF.probHit += 1.;
      pCDF.nentries++;
    }

    // Checks
    if (fabs(to - pCDF.theta_deg) > 0.1) {
      printf("Reading wrong file \n");
      continue;
    }
    if (fabs(log10(keo) - pCDF.lke) > 0.01) {
      printf("Reading wrong file \n");
      continue;
    }

  }
  fclose(fp);

  //-----------------------------------
  //  Set (probZero,probCut,probDecay,probHit,mean,rms,CDF)
  //-----------------------------------
  if (fDetectorType == 2) {
    if (pCDF.nentries > 0) pCDF.probHit = pCDF.probHit / double(pCDF.nentries) * MDAreaRatio;
    else                   pCDF.probHit = 0;
  } else {
    for (int icase = 0; icase < 2; icase++) {
      PeCDF* thisCDF = (icase == 0 ? &pCDF : &pCDF_dec);
      if (thisCDF->nentries > 0) {
        thisCDF->probZero = thisCDF->probZero / double(thisCDF->nentries);
        thisCDF->probDecay = thisCDF->probDecay / double(thisCDF->nentries);
        thisCDF->probCut  = thisCDF->probCut / double(thisCDF->nentries);
        thisCDF->mean /= double(thisCDF->nentries);
        thisCDF->rms = sqrt(double(thisCDF->nentries) / double(thisCDF->nentries - 1) * (thisCDF->rms / double(thisCDF->nentries) - thisCDF->mean * thisCDF->mean));
        SetUpCDFs(thisCDF);
      } else {
        thisCDF->probZero = 1.;
        thisCDF->probDecay = 0.;
        thisCDF->probCut = 1.;
        thisCDF->mean = 0.;
        thisCDF->rms = 0.;
      }
    }
  }

  if (fDetectorType == 2) return true;

  return true;
}


void
DetectorResponse::SetUpCDFs(PeCDF* pCDF)
{

  /*Order CDFs in increasing of order of pes*/
  sort(pCDF->prob.begin(),    pCDF->prob.end(), peordering);

  /*Set Mean/RMS/Nentries */
  double nentries = double(pCDF->nentries);
  double nentriesCut = double(pCDF->nentriesCut);
  pCDF->median = 0;

  /*CDF when n_pe < NPECUT */
  if (NPECUT == UNDEF)
    pCDF->PeCut = (pCDF->prob.size() > 0 ? pCDF->prob[0].first : 0);
  else {
    pCDF->PeCut = NPECUT;
    for (int ipe = 0; ipe < NPECUT; ipe++) {
      double acc = 0.;
      for (unsigned int ii = 0; ii < pCDF->prob.size(); ii++) {
        if (pCDF->prob[ii].first > (NPECUT - 1)) break;
        if (pCDF->prob[ii].first <= ipe) acc += pCDF->prob[ii].second;
      }
      if (nentriesCut > 0.) acc = acc / double(nentriesCut);
      else acc = 0.;

      pCDF->cdf.push_back(pair<int, double>(ipe, acc));
    }
    if (nentries == nentriesCut) {
      pCDF->prob.clear();
      return;
    }
  }

  /*CDF when n_pe >=NPECUT */
  int nbins;
  if ((nentries - nentriesCut) < 10.) nbins = 1;
  else nbins = int(NCDFMAX < (nentries - nentriesCut) / 10. ?  NCDFMAX : (nentries - nentriesCut) / 10.);
  int step = int(nentries - nentriesCut) / nbins;

  double cnt_acc = 0., cnt = 0.;
  for (unsigned int ipe = 0; ipe < pCDF->prob.size(); ipe++) {
    int npe = pCDF->prob[ipe].first;
    double nparticles = pCDF->prob[ipe].second;

    if (npe < pCDF->PeCut) continue;          /*Skip n_pe<PeCut*/
    cnt_acc += nparticles;
    cnt += nparticles;              /*Counter*/

    if (pCDF->median == 0. && (cnt_acc + nentriesCut) / nentries > 0.5) pCDF->median = npe; /*Set the Median*/

    if (cnt >= step || ipe == (pCDF->prob.size() - 1)) {
      double prob = cnt_acc / (nentries - nentriesCut);
      pCDF->cdf.push_back(pair<int, double>(npe, prob)) ;
      cnt = 0.;
    }
  }

  pCDF->prob.clear(); /*Clear the single entry CDFs, too much memory space*/
}



//-----------------------------------------------------------------------------
//------------Read/Write CDF files
//-----------------------------------------------------------------------------

bool
DetectorResponse::ReadCDFFile(vector<PeCDF>* fCDFs, int size)
{
  char CDFFileName[200];
  if (fDetectorType == 0)
    sprintf(CDFFileName, "%s/%s.wcd.dat", fDataDir.c_str(), GetParticleName(fCDFs->at(0).itype));
  else if (fDetectorType == 1)
    sprintf(CDFFileName, "%s/%s.scin.dat", fDataDir.c_str(), GetParticleName(fCDFs->at(0).itype));
  else if (fDetectorType == 2)
    sprintf(CDFFileName, "%s/%s.md.dat", fDataDir.c_str(), GetParticleName(fCDFs->at(0).itype));


  FILE* fp;
  fp = fopen(CDFFileName, "r");

  if (!fp) {
    printf("Could not open %s file \n", CDFFileName);
    return false;
  }

  int cnt = 0;
  while (!feof(fp)) {
    if (cnt >= size) break;

    unsigned int nbins = 0;
    int PeCut = 0, itype = -1, nentries = 0, nentriesCut = 0;
    double lke, to_deg, mean, rms, median, probCut, probHit, probZero;
    int nread =
      fscanf(fp, "%d %lf %lf %lf %lf %lf %d %d %lf %lf %lf  %d %u ",
             &itype, &lke, &to_deg, &mean, &rms, &median, &nentries, &nentriesCut, &probZero, &probCut, &probHit, &PeCut, &nbins);
    if (nread != 13) return false;

    bool extrapol_flag;
    int bins[4];
    int ibin = GetBin(lke, to_deg * M_PI / 180., itype, bins, extrapol_flag);

    if (ibin < 0 || ibin > size) return false;

    if (itype !=  fCDFs->at(ibin).itype) return false;
    if (fabs(lke - fCDFs->at(ibin).lke)) return false;
    if (fabs(to_deg - fCDFs->at(ibin).theta_deg)) return false;
    if (fDetectorType < 0 || fDetectorType > 2) return false;
    if (PeCut != NPECUT  && NPECUT != UNDEF)  return false;

    fCDFs->at(ibin).mean = mean;
    fCDFs->at(ibin).rms = rms;
    fCDFs->at(ibin).median = median;
    fCDFs->at(ibin).nentries = nentries;
    fCDFs->at(ibin).nentriesCut = nentriesCut;
    fCDFs->at(ibin).PeCut = PeCut;
    fCDFs->at(ibin).probCut = probCut;
    fCDFs->at(ibin).probZero = probCut;
    fCDFs->at(ibin).probHit = probHit;

    for (unsigned int i = 0; i < nbins; i++) {
      int npe = 0;
      double prob = 0.;
      if (fscanf(fp, " %d %lf ", &npe, &prob) != 2) return false;
      fCDFs->at(ibin).cdf.push_back(pair<int, double>(npe, prob));
    }
    if (fCDFs->at(ibin).cdf.size() != nbins) return false;
    cnt++;
  }

  if (cnt != size) return false;

  return true;
}


bool
DetectorResponse::WriteCDFFile(vector<PeCDF>* fCDFs)
{
  //---Write CDF file
  char CDFFileName[200];
  if (fDetectorType == 0)
    sprintf(CDFFileName, "%s/%s.wcd.dat", fDataDir.c_str(), GetParticleName(fCDFs->at(0).itype));
  else if (fDetectorType == 1)
    sprintf(CDFFileName, "%s/%s.scin.dat", fDataDir.c_str(), GetParticleName(fCDFs->at(0).itype));
  else if (fDetectorType == 2)
    sprintf(CDFFileName, "%s/%s.md.dat", fDataDir.c_str(), GetParticleName(fCDFs->at(0).itype));

  printf("Writing %s CDFFileName \n", CDFFileName);

  FILE* fp = fopen(CDFFileName, "w");
  if (!fp) {
    printf("Could not open %s file \n", CDFFileName);
    return false;
  }

  for (unsigned int ibin = 0; ibin < fCDFs->size(); ibin++) {

    fprintf(fp, "%d %e %4.2lf %e %e %e %d %d %e %e %e %d %zu \n ",
            fCDFs->at(ibin).itype, fCDFs->at(ibin).lke,
            fCDFs->at(ibin).theta_deg, fCDFs->at(ibin).mean, fCDFs->at(ibin).rms,
            fCDFs->at(ibin).median, fCDFs->at(ibin).nentries, fCDFs->at(ibin).nentriesCut,
            fCDFs->at(ibin).probZero, fCDFs->at(ibin).probCut, fCDFs->at(ibin).probHit, fCDFs->at(ibin).PeCut, fCDFs->at(ibin).cdf.size());

    for (unsigned int i = 0; i < fCDFs->at(ibin).cdf.size(); i++)
      fprintf(fp, " %d  %e \n ", fCDFs->at(ibin).cdf[i].first, fCDFs->at(ibin).cdf[i].second);
  }
  fclose(fp);

  return true;
}


//-----------------------------------------------------------------------------
//------------Signal Interpolation utilities
//-----------------------------------------------------------------------------


//--Linear Interpolation
double
DetectorResponse::Interpolate(double* x, double* y, double x0)
{
  if (x[0] == x[1] && y[0] == y[1]) return y[0];

  double y0 = (y[1] - y[0]) / (x[1] - x[0]) * (x0 - x[0]) + y[0];
  return y0;
}

//---Calculate the signal given a value of the CDF
double
DetectorResponse::GetSignalBin(double f, PeCDF& fPeCDF)
{
  //---Give back the Mean
  if (f == -100) return fPeCDF.mean;

  //-Get the bins
  int start = (NPECUT == UNDEF ? 0 : NPECUT);
  int end = fPeCDF.cdf.size() ;
  int size = end - start;

  if (size == 0) return fPeCDF.PeCut;

  int icdf = 0;
  for (int i = start; i < end ; i++)
    if (fPeCDF.cdf[i].second > f) {
      icdf = i;
      break;
    }


  double x[2], y[2];
  if (icdf == start) {
    x[0] = 0.;
    x[1] = fPeCDF.cdf[icdf].second;
    y[0] = fPeCDF.PeCut - 1;
    y[1] = fPeCDF.cdf[icdf].first;
  } else {
    if (icdf == (end - 1) && size > 2) icdf = icdf - 1;

    x[0] = fPeCDF.cdf[icdf - 1].second;
    x[1] = fPeCDF.cdf[icdf].second;
    y[0] = fPeCDF.cdf[icdf - 1].first;
    y[1] = fPeCDF.cdf[icdf].first;
  }
  double res = Interpolate(x, y, f);

  return round(res + 0.5);

}

//-----------------------------------------------------------------------------
// Hit probability ( only for muons in the Muon Detector )
//     Interpolation in sec(theta)
//     theta in radians
//     ke in GeV
//-----------------------------------------------------------------------------
bool DetectorResponse::IsHit(double ke, double theta, int itype , double f)
{
  if (itype != 0) return false;
  if (fDetectorType != 2) return false;
  if (f <= 0.) return false;
  if (ke > 10.) return true;

  bool outsideRange = false;
  vector <PeCDF>* fCDFs = &fMuonCDFs;
  unsigned int nx = nTo_Grid[iGrid], ny = nKE_Muon[iGrid];

  //Check on f
  if (!(f == -100) &&  !(f > 0 && f < 1)) return false;

  //Check completeness of CDF
  if (fCDFs->size() < nx * ny) {
    printf("Warning: there are no enough CDFs for this particle specie %zu , %u %u || TERMINATE \n", fCDFs->size(), nx, ny);
    return false; //--Terminate the program
  }

  //x value to interpolate
  double lke0 = log10(ke);
  double track0 = GetRelativeTrack(theta);

  //--Bins to be used for interpolation
  int ibins[4];  // 0->(to1,ke1) 1->(to1,ke2) 2->(to2,ke1) 3->(to2,ke2)
  GetBin(lke0 , theta, itype, ibins, outsideRange);

  //Check entries
  bool entriesFlag = (fCDFs->at(ibins[0]).nentries == 0 || fCDFs->at(ibins[1]).nentries == 0 ||
                      fCDFs->at(ibins[2]).nentries == 0 || fCDFs->at(ibins[3]).nentries == 0);

  if (entriesFlag) {
    printf("Warning: one of the bins needed to interpolate have no entries, in both branches, Particle : %d  theta %5.2lf log(ke) %5.2lf \n", itype, theta * 180. / M_PI, lke0);
    return false;
  }

  //Set x values of interpolation
  double lke[2] = { fCDFs->at(ibins[0]).lke, fCDFs->at(ibins[1]).lke };
  double to[2] = { fCDFs->at(ibins[0]).theta, fCDFs->at(ibins[2]).theta};
  double track[2];
  track[0] = GetRelativeTrack(to[0]);
  track[1] = GetRelativeTrack(to[1]);

  //-----------------------------------------------------------------------------------------------------------------

  double ptheta[2];
  double hit[2];

  hit[0] = fCDFs->at(ibins[0]).probHit;
  hit[1] = fCDFs->at(ibins[1]).probHit;
  ptheta[0] = Interpolate(lke, hit, lke0);


  hit[0] = fCDFs->at(ibins[2]).probHit;
  hit[1] = fCDFs->at(ibins[3]).probHit;
  ptheta[1] = Interpolate(lke, hit, lke0);

  double probHit = Interpolate(track, ptheta, track0);

  if (f < probHit)  return true;
  else            return false;

}

//-----------------------------------------------------------------------------
// Signal
//-----------------------------------------------------------------------------

//It returns the signal in Total Number of PEs
//If the muon decays the electron signal is not added, unless the Mean Signal is requested
//Theta in radians, ke in GeV, particle Id
double DetectorResponse::GetPe(double ke, double theta, int itype , double f)
{
  //--Signal ~ E  for EM particles (calorimeter approach) and independent on Theta
  //--Signal ~ <track> for Muons and independent on energy

  //-- itype = 0  (Muons/AntiMuons)   1 (Electrons/Positrons)  2 (Gammas)  3 (Muons/AntiMuons enlarge area)  4 ( Michel electron )  5 ( Michel electron enlarge area )

  //--Therefore:
  //                       1) Interpolations in E are done in (log Signal,log E)
  //                       2) Interpolation in Theta are done in ( Signal, Track(Theta) )

  if (fDetectorType == 2) return 0.;

  if (fUseEnlargeArea && itype == 0 && log10(ke) < lKE_MuonLow[fDetectorType][nKE_MuonLow[fDetectorType] - 1]) itype = 3;
  if (fUseEnlargeArea && itype == 4 && log10(ke) < lKE_MuonLow[fDetectorType][nKE_MuonLow[fDetectorType] - 1]) itype = 5;

  vector <PeCDF>* fCDFs;
  unsigned int nx = nTo_Grid[iGrid], ny;
  bool outsideRange = false;

  //Check on f
  if (!(f == -100) &&  !(f >= 0. && f <= 1.)) return 0.;

  //Select CDF
  if (itype == 0)      {
    ny = nKE_Muon[iGrid];
    fCDFs = &fMuonCDFs;
  } else if (itype == 1) {
    ny = nKE_Electron[iGrid];
    fCDFs = &fElectronCDFs;
  } else if (itype == 2) {
    ny = nKE_Gamma[iGrid];
    fCDFs = &fGammaCDFs;
  } else if (itype == 3) {
    ny = nKE_MuonLow[iGrid];
    fCDFs = &fMuonLowCDFs;
  } else if (itype == 4) {
    ny = nKE_Muon[iGrid];
    fCDFs = &fMichelElectronCDFs;
  } else if (itype == 5) {
    ny = nKE_MuonLow[iGrid];
    fCDFs = &fMichelElectronLowCDFs;
  } else  return 0.;

  //Check completeness of CDF
  if (fCDFs->size() < nx * ny) {
    printf("Warning: there are no enough CDFs for this particle specie %zu , %u %u || TERMINATE \n", fCDFs->size(), nx, ny);
    return 0.; //--Terminate the program
  }

  //x value to interpolate
  double lke0 = log10(ke);
  double track0 = GetRelativeTrack(theta);

  //--Bins to be used for interpolation
  int ibins[4];  // 0->(to1,ke1) 1->(to1,ke2) 2->(to2,ke1) 3->(to2,ke2)
  GetBin(lke0 , theta, itype, ibins, outsideRange);

  //Check entries
  bool entriesFlag = (fCDFs->at(ibins[0]).nentries == 0 || fCDFs->at(ibins[1]).nentries == 0 ||
                      fCDFs->at(ibins[2]).nentries == 0 || fCDFs->at(ibins[3]).nentries == 0);
  if (entriesFlag) {
    printf("Warning: one of the bins needed to interpolate have no entries, Particle : %d  theta %5.2lf log(ke) %5.2lf \n", itype, theta * 180. / M_PI, lke0);
    return 0.;
  }

  //Set x values of interpolation
  double lke[2] = { fCDFs->at(ibins[0]).lke, fCDFs->at(ibins[1]).lke };
  double to[2] = { fCDFs->at(ibins[0]).theta, fCDFs->at(ibins[2]).theta};
  double track[2];
  track[0] = GetRelativeTrack(to[0]);
  track[1] = GetRelativeTrack(to[1]);


  //-----------------------------------------------------------------------------------------------------------------
  //----Interpolate the Probability for 0,1,...NPECUT-1 Photoelectrons
  //-----------------------------------------------------------------------------------------------------------------
  double fcorr = f;
  if (f != -100 && NPECUT != UNDEF) {
    double probCut = 0.;

    double zero[2], ptheta[2];
    //int entries[2];
    zero[0] = fCDFs->at(ibins[0]).probCut;
    zero[1] = fCDFs->at(ibins[1]).probCut;
    ptheta[0] = Interpolate(lke, zero, lke0);
    //entries[0]=fCDFs->at(ibins[0]).nentries; entries[1]=fCDFs->at(ibins[1]).nentries;

    zero[0] = fCDFs->at(ibins[2]).probCut;
    zero[1] = fCDFs->at(ibins[3]).probCut;
    //entries[0]=fCDFs->at(ibins[2]).nentries; entries[1]=fCDFs->at(ibins[3]).nentries;
    ptheta[1] = Interpolate(lke, zero, lke0);
    probCut = Interpolate(track, ptheta, track0);


    if (fcorr <  probCut) {

      //---Calculate the Interpolated probability <= 0...9 Pes
      double probLow = 0.;
      double npe = 0.;
      for (int jj = 0; jj < NPECUT; jj++) {
        zero[0] = fCDFs->at(ibins[0]).cdf[jj].second;
        zero[1] = fCDFs->at(ibins[1]).cdf[jj].second;
        ptheta[0] = Interpolate(lke, zero, lke0);

        zero[0] = fCDFs->at(ibins[2]).cdf[jj].second;
        zero[1] = fCDFs->at(ibins[3]).cdf[jj].second;
        ptheta[1] = Interpolate(lke, zero, lke0);
        probLow = Interpolate(track, ptheta, track0);

        npe = double(jj);
        if (probLow > fcorr / probCut) break;
      }
      return npe;
    }
    fcorr = (fcorr - probCut) / (1 - probCut);
  }

  //-----------------------------------------------------------------------------------------------------------
  //----Interpolate when  >NPECUT Photoelectrons
  //-----------------------------------------------------------------------------------------------------------
  {

    double signal[2], lsignal[2];

    //---Interpolation in Energy for Theta1
    signal[0] = double(GetSignalBin(fcorr, fCDFs->at(ibins[0])));
    signal[1] = double(GetSignalBin(fcorr, fCDFs->at(ibins[1])));

    lsignal[0] = log10(signal[0]);
    lsignal[1] = log10(signal[1]);
    double s1;
    if (signal[0] > 0. && signal[1] > 0.)
      s1 = pow(10., Interpolate(lke, lsignal, lke0));
    else
      s1 = Interpolate(lke, signal, lke0);

    //---Interpolation in Energy for Theta2
    signal[0] = double(GetSignalBin(fcorr, fCDFs->at(ibins[2])));
    signal[1] = double(GetSignalBin(fcorr, fCDFs->at(ibins[3])));

    lsignal[0] = log10(signal[0]);
    lsignal[1] = log10(signal[1]);

    double s2;
    if (signal[0] > 0. && signal[1] > 0.)
      s2 = pow(10., Interpolate(lke, lsignal, lke0));
    else
      s2 = Interpolate(lke, signal, lke0);


    double y[2] = {s1, s2};
    double s3 = Interpolate(track, y, track0);


    if (f == -100.) return s3;
    return         round(s3);
  }

}