/*
 * TrkDedxHCut.cpp
 *
 *  Created on: 26-mar-2009
 *      Author: Nico De Simone, S. Ricciarini
 */

/*! @file TrkDedxHCut.cpp The TrkDedxHCut class implementation file */

#include "TrkDedxHCut.h"

const Double_t TrkDedxHCut::_xLow[] = { 0.126414, 0.140142, 0.152714, 0.165854, 0.178357, 0.194244, 0.211669, 0.242859,
    0.264646, 0.293385, 0.314258, 0.360565, 0.392911, 0.450807, 0.508423, 0.603735, 0.680897, 0.754838, 0.851312,
    1.00032, 1.36201, 1.37, 1000. }; // _nLow=23 elements

const Double_t TrkDedxHCut::_yLow[] = { 18.9221, 16.0752, 14.2214, 12.6324, 11.2421, 9.71928, 8.1965, 6.93856, 6.14407,
    5.15095, 4.6875, 3.8268, 3.49576, 2.89989, 2.56886, 2.17161, 1.90678, 1.70816, 1.50953, 1.27754, 0.715042, 0, 0 }; // _nLow=23 elements

const Double_t TrkDedxHCut::_xHigh[] = { 0.126414, 0.392911, 0.474654, 0.508423, 0.563635, 0.680897, 0.76792, 0.851312,
    0.960115, 1.12068, 1.22122, 1.50085, 1.69267, 1.87649, 2.1903, 2.65934, 3.75055, 6.21589, 11.0139, 24.2773, 1000. }; // _nHigh=21 elements

const Double_t TrkDedxHCut::_yHigh[] = { 24.3512, 22.696, 21.107, 20.0477, 18.5249, 15.8766, 13.6917, 11.5731, 10.1827,
    8.26271, 7.2696, 5.5482, 5.01854, 4.6875, 4.15784, 3.83697, 3.59789, 3.42119, 3.36335, 3.36335, 3.36335 }; // _nHigh=21 elements


int TrkDedxHCut::Check(PamLevel2 *event) {

  TrkTrack *trkTrack = event->GetTrack(0)->GetTrkTrack();

  Float_t rigMod = trkTrack->GetRigidity();
  Float_t dedx = trkTrack->GetDEDX();

#ifdef DEBUGPAMCUT

  h_trk_he_dedx_rigmod[0]->Fill(rigMod, dedx);
  h_trk_le_dedx_rigmod[0]->Fill(rigMod, dedx);

#endif

  if (rigMod < _xHigh[0])
    return 0; // a proton below 126.414 MV is rejected

  // interpolation of _xLow
  for (int i = 0; i < _nLow - 1; i++) {
    if (_xLow[i] <= rigMod && rigMod < _xLow[i + 1]) {
      double slope = (_yLow[i + 1] - _yLow[i]) / (_xLow[i + 1] - _xLow[i]);
      double y = _yLow[i] + slope * (rigMod - _xLow[i]);

      if (dedx < y)
        return 0;
    }
  }

  // interpolation of _xHigh
  for (int i = 0; i < _nHigh - 1; i++) {
    if (_xHigh[i] <= rigMod && rigMod < _xHigh[i + 1]) {
      double slope = (_yHigh[i + 1] - _yHigh[i]) / (_xHigh[i + 1] - _xHigh[i]);
      double y = _yHigh[i] + slope * (rigMod - _xHigh[i]);

      if (dedx > y)
        return 0;
    }
  }

  // extrapolation above 1 TV
  if (rigMod >= _xHigh[_nHigh - 1]) {
    double slope = (_yHigh[_nHigh - 1] - _yHigh[_nHigh - 2]) / (_xHigh[_nHigh - 1] - _xHigh[_nHigh - 2]);
    double y = _yHigh[_nHigh - 1] + slope * (rigMod - _xHigh[_nHigh - 1]);

    if (dedx > y)
      return 0;
  }

#ifdef DEBUGPAMCUT

  h_trk_he_dedx_rigmod[1]->Fill(rigMod, dedx);
  h_trk_le_dedx_rigmod[1]->Fill(rigMod, dedx);

#endif

  return CUTOK;

}