/*
 * TrkDedxHeCut.cpp
 *
 *  Created on: 28/gen/2010
 *      Author: Nicola Mori
 */

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

#include "TrkDedxHeCut.h"

int TrkDedxHeCut::Check(PamLevel2 *event) {

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

  float R = track->GetRigidity();
  //float dEdx = track->GetDEDX();

  // New smart dE/dx computation
  // This procedure removes high releases and computes dE/dx
  float totReleaseX = 0., maxReleaseX = 0.;
  float totReleaseY = 0., maxReleaseY = 0.;
  int nX = 0, nY = 0;
  float dEdxView = 0;
  for (int ip = 0; ip < 6; ip++) {
    // X view
    dEdxView = track->GetDEDX(ip, 0);
    if (dEdxView > 0 && track->XGood(ip)) {
      totReleaseX += dEdxView;
      nX++;
      if (dEdxView > maxReleaseX)
        maxReleaseX = dEdxView;
    }
    // Y view
    dEdxView = track->GetDEDX(ip, 1);
    if (dEdxView > 0 && track->YGood(ip)) {
      totReleaseY += dEdxView;
      nY++;
      if (dEdxView > maxReleaseY)
        maxReleaseY = dEdxView;
    }
  }

  // Discard highest release, eventually
  if (maxReleaseX > -9. + 4. * totReleaseX / (float) nX) {
    totReleaseX -= maxReleaseX;
    nX--;
  }
  if (maxReleaseY > -8. + 4. * totReleaseY / (float) nY) {
    totReleaseY -= maxReleaseY;
    nY--;
  }

  // Compute dE/dx
  //float dEdx = (totReleaseX + totReleaseY) / (float) (nX + nY);
  float dEdxX = totReleaseX / nX;
  float dEdxY = totReleaseY / nY;

  // Analyze the event
  // X
  float denLow = pow(R, 1.8); // The power of the denominator is the same for X and Y (with current calibration)
  float denHigh = pow(R, 1.5);
  if (dEdxX < 3.7 + 4.6 / denLow)
    return 0;
  if (dEdxX > 8.9 + 17. / denHigh)
    return 0;
  // Y
  if (dEdxY < 3.3 + 4.9 / denLow)
    return 0;
  if (dEdxY > 8.0 + 17. / denHigh)
    return 0;
  //X+Y
  /*if (dEdx < 3.7 + 4.8 / (R * R))
    return 0;
  if (dEdx > 8.9 + 17. / (R * R))
    return 0;*/

  return CUTOK;
}