/*
 * 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();
  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;
  /*
   ToFLevel2 *tofL2 = event->GetToFLevel2();
   int trkSeqNo = 0;
   for (int i = 0; i < event->GetToFLevel2()->ntrk(); i++) {
   if (tofL2->GetToFTrkVar(i)->trkseqno == track->GetSeqNo()) {
   trkSeqNo = i;
   break;
   }
   }
   if ((tofL2->GetdEdx(trkSeqNo, 0, 100) + tofL2->GetdEdx(trkSeqNo, 1, 100) + tofL2->GetdEdx(trkSeqNo, 2, 100)
   + tofL2->GetdEdx(trkSeqNo, 3, 100)) / 4. < 2) {
   cout << "Check" << endl;
   for (int ip = 0; ip < 6; ip++) {
   cout << "TOF" << ip << ": " << tofL2->GetdEdx(trkSeqNo, ip, 100) << endl;
   dEdxView = track->GetDEDX(ip, 0);
   if (dEdxView > 0 && track->XGood(ip))
   cout << "X" << ip << ": " << dEdxView << endl;
   dEdxView = track->GetDEDX(ip, 1);
   if (dEdxView > 0 && track->YGood(ip))
   cout << "Y" << ip << ": " << dEdxView << endl;
   }
   }
   */
  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;
}