/*
 * TrkDedxHCut.cpp
 *
 *  Created on: 15-mar-2010
 *      Author: 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

Int_t TrkDedxHCut::Check(PamLevel2 *event) {

//  TrkTrack *trkTrack = event->GetTrack(0)->GetTrkTrack();
    if(event->GetNTracks(_trkAlg)==0)return 0;    
    ExtTrack *trkTrack = event->GetTrack(0,_trkAlg)->GetExtTrack();

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

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

  // lower limit taken from N. De Simone (interpolation of _xLow)
  
  for (Int_t i = 0; i < _nLow - 1; i++) {
    if (_xLow[i] <= rigMod && rigMod < _xLow[i + 1]) {
      Double_t slope = (_yLow[i + 1] - _yLow[i]) / (_xLow[i + 1] - _xLow[i]);
      Double_t y = _yLow[i] + slope * (rigMod - _xLow[i]);

      if (dedx < y)
        return 0;
    }
  }

  // higher limit taken from N. Mori lower limit for Helium

  // Compute dE/dx for X and Y separately
  Float_t totReleaseX = 0.;
  Float_t totReleaseY = 0.;
  Int_t nX = 0, nY = 0;
  Float_t dEdxView = 0;

  for (Int_t ip = 0; ip < 6; ip++) {
    // X view
    dEdxView = trkTrack->GetDEDX(ip, 0);
    if (dEdxView > 0 && trkTrack->XGood(ip)) {
      totReleaseX += dEdxView;
      nX++;
    }
    // Y view
    dEdxView = trkTrack->GetDEDX(ip, 1);
    if (dEdxView > 0 && trkTrack->YGood(ip)) {
      totReleaseY += dEdxView;
      nY++;
    }
  }

  Float_t dEdxX = totReleaseX / nX; // nX assumed > 0
  Float_t dEdxY = totReleaseY / nY; // nY assumed > 0
  
  Float_t denHigh = pow(rigMod, 1.8); // The power of the denominator is the same for X and Y (with current calibration)
  if (dEdxX > 3.7 + 4.6 / denHigh) {
    return 0;
  }
  if (dEdxY > 3.3 + 4.9 / denHigh) {
    return 0;
  }
   
  return CUTOK;

}