#include "StHbtMaker/Cut/gregsTrackCut.h"
#include <sstream>

#ifdef __ROOT__
  ClassImp(gregsTrackCut)
#endif

//#define STHBTDEBUG

//_________________
gregsTrackCut::gregsTrackCut() {
  mNTracksPassed = mNTracksFailed = 0;
  mNSigmaElectron[0] = -1e9;  mNSigmaElectron[1] = +1e9;
  mNSigmaPion[0]     = -1e9;  mNSigmaPion[1]     = +1e9;
  mNSigmaKaon[0]     = -1e9;  mNSigmaKaon[1]     = +1e9;
  mNSigmaProton[0]   = -1e9;  mNSigmaProton[1]   = +1e9;
  mEta[0]      = -1e9;   mEta[1]      = +1e9; 
  mRapidity[0] = -1e9;   mRapidity[1] = +1e9; 
  mP[0] = -1e9;  mP[1] = +1e9;
  mPt[0] = -1e9;  mPt[1] = +1e9;
  mPx[0] = -1e9;  mPx[1] = +1e9;
  mPy[0] = -1e9;  mPy[1] = +1e9;
  mPz[0] = -1e9;  mPz[1] = +1e9;
  mDCA[0] = -1e9;  mDCA[1] = +1e9;
  mDCAGlobal[0] = -1e9;  mDCAGlobal[1] = +1e9;
  mNHits[0] = 0; mNHits[1] = 60;
  mNHitsFit[0] = 0; mNHitsFit[1] = 60;
  mNdEdxHits[0] = 0; mNdEdxHits[1] = 60;
  mTofMassSqr[0] = -1e9; mTofMassSqr[1] = +1e9;
  mAntiSplit = 0.;
  mCharge = 0;
  mPIDSelection = 2;
  mType = 1;
  mPartId = 4; // kaon selection by befault
  mTpcTrackMom[0] = -1e9; mTpcTrackMom[1] = +1e9;
  mTofTrackMom[0] = -1e9; mTofTrackMom[1] = +1e9;
  mTpcTrackPt[0] = -1e9; mTpcTrackPt[1] = +1e9;
  mTofTrackPt[0] = -1e9; mTofTrackPt[1] = +1e9;
  mInverseBetaDifferce[0] = -1e9; mInverseBetaDifferce[1]=+1e9;
  mTofNSigmaElectron[0] = -3; mTofNSigmaElectron[1] = 3;
  mTofNSigmaPion[0] = -3; mTofNSigmaPion[1] = 3;
  mTofNSigmaKaon[0] = -3; mTofNSigmaKaon[1] = 3;
  mTofNSigmaProton[0] = -3; mTofNSigmaProton[1] = 3;
}

//_________________
gregsTrackCut::gregsTrackCut(gregsTrackCut& c) : StHbtTrackCut(c) {
  mNTracksPassed = mNTracksFailed = 0;
  mCharge = c.mCharge;
  mType = c.mType;
  mAntiSplit = c.mAntiSplit;
  mNSigmaElectron[0] = c.mNSigmaElectron[0];  mNSigmaElectron[1] = c.mNSigmaElectron[1];
  mNSigmaPion[0] = c.mNSigmaPion[0];  mNSigmaPion[1] = c.mNSigmaPion[1];
  mNSigmaKaon[0] = c.mNSigmaKaon[0];  mNSigmaKaon[1] = c.mNSigmaKaon[1];
  mNSigmaProton[0] = c.mNSigmaProton[0];  mNSigmaProton[1] = c.mNSigmaProton[1];
  mNHits[0] = c.mNHits[0];  mNHits[1] = c.mNHits[1];
  mNHitsFit[0] = c.mNHitsFit[0];  mNHitsFit[1] = c.mNHitsFit[1];
  mNdEdxHits[0] = c.mNdEdxHits[0];  mNdEdxHits[1] = c.mNdEdxHits[1];
  mP[0] = c.mP[0]; mP[1] = c.mP[1];
  mPt[0] = c.mPt[0]; mPt[1] = c.mPt[1];
  mPx[0] = c.mPx[0]; mPx[1] = c.mPx[1];
  mPy[0] = c.mPy[0]; mPy[1] = c.mPy[1];
  mPz[0] = c.mPz[0]; mPz[1] = c.mPz[1];
  mRapidity[0] = c.mRapidity[0];  mRapidity[1] = c.mRapidity[1];
  mEta[0] = c.mEta[0];  mEta[1] = c.mEta[1];
  mDCA[0] = c.mDCA[0];  mDCA[1] = c.mDCA[1];
  mDCAGlobal[0] = c.mDCAGlobal[0];  mDCAGlobal[1] = c.mDCAGlobal[1];
  mTofMassSqr[0] = c.mTofMassSqr[0];  mTofMassSqr[1] = c.mTofMassSqr[1];
  mPIDSelection = c.mPIDSelection;  
  mTpcTrackMom[0] = c.mTpcTrackMom[0];  mTpcTrackMom[1] = c.mTpcTrackMom[1];  
  mTofTrackMom[0] = c.mTofTrackMom[0];  mTofTrackMom[1] = c.mTofTrackMom[1];  
  mTpcTrackPt[0] = c.mTpcTrackPt[0];  mTpcTrackPt[1] = c.mTpcTrackPt[1];  
  mTofTrackPt[0] = c.mTofTrackPt[0];  mTofTrackPt[1] = c.mTofTrackPt[1];
  mPartId = c.mPartId;
  mInverseBetaDifferce[0] = c.mInverseBetaDifferce[0];  mInverseBetaDifferce[1] = c.mInverseBetaDifferce[1];
  mTofNSigmaElectron[0] = c.mTofNSigmaElectron[0];  mTofNSigmaElectron[1] = c.mTofNSigmaElectron[1];
  mTofNSigmaPion[0] = c.mTofNSigmaPion[0];  mTofNSigmaPion[1] = c.mTofNSigmaPion[1];
  mTofNSigmaKaon[0] = c.mTofNSigmaKaon[0];  mTofNSigmaKaon[1] = c.mTofNSigmaKaon[1];
  mTofNSigmaProton[0] = c.mTofNSigmaProton[0];  mTofNSigmaProton[1] = c.mTofNSigmaProton[1];
  mNTracksPassed=0;
  mNTracksFailed=0;
#ifdef STHBTDEBUG
  cout << " gregsTrackCut::gregsTrackCut(gregsTrackCut& c) " << endl;
#endif
}

//_________________
gregsTrackCut::~gregsTrackCut() {
  /* empty */
}

//_________________
bool gregsTrackCut::Pass(const StHbtTrack* track) {

#ifdef STHBTDEBUG
  cout << endl << " *** gregsTrackCut::Pass(const StHbtTrack* track) " << endl;
#endif

  float TEnergy = ::sqrt(track->P().mag2()+mMass*mMass);
  float TRapidity = 0.5*::log((TEnergy+track->P().z())/
			    (TEnergy-track->P().z()));
  float TPseudoRap = 0.5*::log( (track->P().mag()+track->P().z()) /
				(track->P().mag()-track->P().z()) );
  float aSplit = (float)track->NHitsFit()/track->NHitsPossible();

  //Checking the particle ID
  /*
  if(mPartId == 1) {
    cout << "gregsTrackCut::Pass() - mPartId=1 (electron)" << endl;
  }
  else if(mPartId == 3) {
    cout << "gregsTrackCut::Pass() - mPartId=3 (pion)" << endl;
  }
  else if(mPartId == 4) {
    cout << "gregsTrackCut::Pass() - mPartId=4 (kaon)" << endl;
  }
  else if(mPartId == 5) {
    cout << "gregsTrackCut::Pass() - mPartId=5 (proton)" << endl;
  }
  else {
    cout << "gregsTrackCut::Pass() - Unknown mPartId: " << mPartId 
	 << " Resetting to mPartId=4 (kaon)" << endl;
    mPartId=4;
  }
  */
  if(mPartId<1 || mPartId>5 || mPartId==2) {
    cout << "gregsTrackCut::Pass() - Unknown mPartId: " << mPartId
	 << " Resetting to mPartId=4 (kaon)" << endl;
    mPartId=4;
  }

  //Check 1/beta - 1/beta_{mass_hypothesis}
  float mInvBetaDifference = 0.;
  float mInvBetaPredicted = -999.;
  if(track->TofBeta() > 0) {
    mInvBetaPredicted = ::sqrt(mMass*mMass+track->P().mag2())/track->P().mag();
    mInvBetaDifference = 1./track->TofBeta() - mInvBetaPredicted;
  }

#ifdef STHBTDEBUG
  cout<< "NSigmaElec = " << track->NSigmaElectron() << endl
      << "NSigmaPion = " << track->NSigmaPion() << endl
      << "NSigmaKaon = " << track->NSigmaKaon() << endl
      << "NSigmaProt = " << track->NSigmaProton() << endl
      << "DCAxy      = " << track->DCAxy() << endl
      << "DCAGlobal  = " << track->DCAxyGlobal() << endl
      << "NHits      = " << track->NHits() << endl
      << "NHitsFit   = " << track->NHitsFit() << endl
      << "NHitsPoss  = " << track->NHitsPossible() << endl
      << "NHitsDedx  = " << track->NHitsDedx() << endl
      << "P          = " << track->P().mag() << endl
      << "Pt         = " << track->Pt() << endl
      << "Rapidity   = " << TRapidity << endl
      << "Charge     = " << track->Charge() << endl
      << "Type       = " << track->TrackType()<< endl
      << "TofMassSqr = " << track->TofMassSqr()<< endl
      << "TofBeta    = " << track->TofBeta()<< endl
      << "AntiSplit  = " << aSplit << endl
      << "PseudoRap  = " << TPseudoRap << endl
      << "PartId     = " << mPartId << endl
      << "SetMass    = " << mMass << endl
      << "InvBetaDiff= " << mInvBetaDifference << endl;
#endif

  bool mGoodCharge = (track->Charge() == mCharge);
  
  bool mGoodTrack=( true &&
		   (track->TrackType() == mType) &&
		   (aSplit >= mAntiSplit ) &&
		   (track->DCAxy() >= mDCA[0]) && (track->DCAxy() <= mDCA[1]) &&
		   (track->DCAxyGlobal() >= mDCAGlobal[0]) && (track->DCAxyGlobal() <= mDCAGlobal[1]) &&
		   (track->NHits() >= mNHits[0]) && (track->NHits() <= mNHits[1]) &&
		   (track->NHitsFit() >= mNHitsFit[0]) && (track->NHitsFit() <= mNHitsFit[1]) &&
		   (track->NHitsDedx() >= mNdEdxHits[0]) && (track->NHitsDedx() <= mNdEdxHits[1]) &&
		   (track->P().mag() >= mP[0]) && (track->P().mag() <= mP[1]) &&
		   (track->Pt() >= mPt[0]) && (track->Pt() <= mPt[1]) &&
		   (track->P().x() >= mPx[0]) && (track->P().x() <= mPx[1]) &&
		   (track->P().y() >= mPy[0]) && (track->P().y() <= mPy[1]) &&
		   (track->P().z() >= mPz[0]) && (track->P().z() <= mPz[1]) &&
		   (TPseudoRap >= mEta[0]) &&
		   (TPseudoRap <= mEta[1]) &&
		   (TRapidity >= mRapidity[0]) &&
		   (TRapidity <= mRapidity[1]) 
		   );

  //No need to check PID if track is parameters are bad
  if(!mGoodCharge || !mGoodTrack) {
    mNTracksFailed++;
    return false;
  }

  bool mGoodPID = false;
  bool mGoodTofTrack = false;
  bool mGoodTofPID = false;
  bool mGoodTpcTrack = false;
  bool mGoodTpcPID = false;
  bool mGoodTnTPID = false;
  bool mGoodToTPID = false;

#ifdef STHBTDEBUG
  cout << "mPIDSelection acceptable values: 0-TPC, 1-ToF, 2-TnT, 3-ToT" << endl;
  cout << "Chosen mPIDSelection: " << mPIDSelection << endl;
#endif

  if(mPIDSelection<0 || mPIDSelection>4) {
    cout << "gregsTrackCut::Pass() - Wrong mPIDSelection: " << mPIDSelection << endl;
    cout << "Resetting to ToT: 3" << endl;
    mPIDSelection = 3;
  }

  //Depends from selection flag: 0-TPC, 1-ToF, 2-TnT, 3-ToT

  /**********************************************/
  if(mPIDSelection == 0) {

#ifdef STHBTDEBUG
    cout << "TPC PID selection" << endl;
#endif

    if( (track->P().mag() >= mTpcTrackMom[0]) &&
	(track->P().mag() <= mTpcTrackMom[1]) &&
	(track->Pt() >= mTpcTrackPt[0]) &&
	(track->Pt() <= mTpcTrackPt[1]) )
      mGoodTpcTrack = true;


    if((mPartId == 1) && mGoodTpcTrack) {

      if( ((track->NSigmaElectron() >= mNSigmaElectron[0]) &&
	   (track->NSigmaElectron() <= mNSigmaElectron[1])) &&
	  ((track->NSigmaPion() <= mNSigmaPion[0]) ||
	   (track->NSigmaPion() >= mNSigmaPion[1])) &&
	  ((track->NSigmaKaon() <= mNSigmaKaon[0]) ||
	   (track->NSigmaKaon() >= mNSigmaKaon[1])) &&
	  ((track->NSigmaProton() <= mNSigmaProton[0]) ||
	   (track->NSigmaProton() >= mNSigmaProton[1])) )
	mGoodTpcPID = true;
      
#ifdef STHBTDEBUG
      cout << "TPC electron PID check: " << mGoodTpcPID << endl;
#endif
    }
    else if((mPartId == 3) && mGoodTpcTrack) {

      if( ((track->NSigmaPion() >= mNSigmaPion[0]) &&
	   (track->NSigmaPion() <= mNSigmaPion[1])) &&
	  ((track->NSigmaElectron() <= mNSigmaElectron[0]) ||
	   (track->NSigmaElectron() >= mNSigmaElectron[1])) &&
	  ((track->NSigmaKaon() <= mNSigmaKaon[0]) ||
	   (track->NSigmaKaon() >= mNSigmaKaon[1])) &&
	  ((track->NSigmaProton() <= mNSigmaProton[0]) ||
	   (track->NSigmaProton() >= mNSigmaProton[1])) )
	mGoodTpcPID = true;
      
#ifdef STHBTDEBUG
      cout << "TPC pion PID check: " << mGoodTpcPID << endl;
#endif
      
    }
    else if((mPartId == 4) && mGoodTpcTrack) {

      if( ((track->NSigmaKaon() >= mNSigmaKaon[0]) &&
	   (track->NSigmaKaon() <= mNSigmaKaon[1])) &&
	  ((track->NSigmaElectron() <= mNSigmaElectron[0]) ||
	   (track->NSigmaElectron() >= mNSigmaElectron[1])) &&
	  ((track->NSigmaPion() <= mNSigmaPion[0]) ||
	   (track->NSigmaPion() >= mNSigmaPion[1])) &&
	  ((track->NSigmaProton() <= mNSigmaProton[0]) ||
	   (track->NSigmaProton() >= mNSigmaProton[1])) )
	mGoodTpcPID = true;

#ifdef STHBTDEBUG
      cout << "TPC kaon PID check: " << mGoodTpcPID << endl;
#endif
   
    }
    else if((mPartId == 5) && mGoodTpcTrack) {

      if( ((track->NSigmaProton() >= mNSigmaProton[0]) &&
	   (track->NSigmaProton() <= mNSigmaProton[1])) &&
	  ((track->NSigmaElectron() <= mNSigmaElectron[0]) ||
	   (track->NSigmaElectron() >= mNSigmaElectron[1])) &&
	  ((track->NSigmaPion() <= mNSigmaPion[0]) ||
	   (track->NSigmaPion() >= mNSigmaPion[1])) &&
	  ((track->NSigmaKaon() <= mNSigmaKaon[0]) ||
	   (track->NSigmaKaon() >= mNSigmaKaon[1])) )
	mGoodTpcPID = true;

#ifdef STHBTDEBUG
      cout << "TPC proton PID check: " << mGoodTpcPID << endl;
#endif
      
    }
  } // if(mPIDSelection == 0)

  /*********************************************/
  else if(mPIDSelection == 1) {
    
#ifdef STHBTDEBUG
    cout << "ToF PID selection" << endl;
#endif

    //Track must have a ToF signal
    if(track->TofBeta() > 0) {
     
      if( (track->P().mag() >= mTofTrackMom[0]) &&
	  (track->P().mag() <= mTofTrackMom[1]) &&
	  (track->Pt() >= mTofTrackPt[0]) &&
	  (track->Pt() <= mTofTrackPt[1]) &&
	  mInvBetaDifference >= mInverseBetaDifferce[0] &&
	  mInvBetaDifference <= mInverseBetaDifferce[1])
	mGoodTofTrack = true;
      
      if( track->TofMassSqr() >= mTofMassSqr[0] &&
	  track->TofMassSqr() <= mTofMassSqr[1] &&
	  mGoodTofTrack)
	mGoodTofPID = true;

#ifdef STHBTDEBUG
      cout << "ToF PID check: " << mGoodTofPID << endl;
#endif
    } //if(track->TofBeta()>0)
  } //else if(mPIDSelection == 1)
  else if(mPIDSelection == 2) {

#ifdef STHBTDEBUG
    cout << "TnT PID selection" << endl;
#endif

    //Track must have a ToF signal
    if(track->TofBeta() > 0) {

      if( (track->P().mag() >= mTofTrackMom[0]) &&
	  (track->P().mag() <= mTofTrackMom[1]) &&
	  (track->Pt() >= mTofTrackPt[0]) &&
	  (track->Pt() <= mTofTrackPt[1]) &&
	  mInvBetaDifference >= mInverseBetaDifferce[0] &&
	  mInvBetaDifference <= mInverseBetaDifferce[1] )
	mGoodTofTrack = true;

      if(mGoodTofTrack && 
	 track->TofMassSqr() >= mTofMassSqr[0] &&
	 track->TofMassSqr() <= mTofMassSqr[1] )
	mGoodTofPID = true;

      if((mPartId==1) && mGoodTofPID) {
	if((track->NSigmaElectron() >= mNSigmaElectron[0]) &&
	   (track->NSigmaElectron() <= mNSigmaElectron[1]))
	  mGoodTnTPID = true;
      } //if((mPartId==1) && mGoodTofPID)
      
      if((mPartId==3) && mGoodTofPID) {
	if( (track->NSigmaPion() >= mNSigmaPion[0]) &&
	    (track->NSigmaPion() <= mNSigmaPion[1]) )
	  mGoodTnTPID = true;
      } //if((mPartId==3) && mGoodTofPID)
      
      if((mPartId==4) && mGoodTofPID) {
	if( (track->NSigmaKaon() >= mNSigmaKaon[0]) &&
	    (track->NSigmaKaon() <= mNSigmaKaon[1]) )
	  mGoodTnTPID = true;
      } //if((mPartId==4) && mGoodTofPID)

      if((mPartId==5) && mGoodTofPID) {
	if( (track->NSigmaProton() >= mNSigmaProton[0]) &&
	    (track->NSigmaProton() <= mNSigmaProton[1]) )
	  mGoodTnTPID = true;
      } //if((mPartId==5) && mGoodTofPID)
	
    } //if(track->TofBeta() > 0)

#ifdef STHBTDEBUG
    cout << "TnT PID check: " << mGoodTnTPID << endl;
#endif
  } //else if(mPIDSelection == 2)
  
  /****************************************************/
  else if(mPIDSelection == 3) {

#ifdef STHBTDEBUG
    cout << "ToT PID selection" << endl;
    cout << "Selecting particle type: " << mPartId << endl;
#endif

    //If track has a ToF signal then we use ToF PID only
    if(track->TofBeta() > 0) {
      
      if( (track->P().mag() >= mTofTrackMom[0]) &&
	  (track->P().mag() <= mTofTrackMom[1]) &&
	  (track->Pt() >= mTofTrackPt[0]) &&
	  (track->Pt() <= mTofTrackPt[1]) &&
	  mInvBetaDifference >= mInverseBetaDifferce[0] &&
	  mInvBetaDifference <= mInverseBetaDifferce[1] )
	mGoodTofTrack = true;
      
      if( track->TofMassSqr() >= mTofMassSqr[0] &&
	  track->TofMassSqr() <= mTofMassSqr[1] &&
	  mGoodTofTrack)
	mGoodToTPID = true;
#ifdef STHBTDEBUG
      cout << "TofCut passed: " << mGoodToTPID << endl;
#endif
    } //if(track->TofBeta() > 0)
    else {
      
      if( (track->P().mag() >= mTpcTrackMom[0]) &&
	  (track->P().mag() <= mTpcTrackMom[1]) &&
	  (track->Pt() >= mTpcTrackPt[0]) &&
	  (track->Pt() <= mTpcTrackPt[1]) )
	mGoodTpcTrack = true;


      if((mPartId == 1) && mGoodTpcTrack) {

	if( ((track->NSigmaElectron() >= mNSigmaElectron[0]) &&
	     (track->NSigmaElectron() <= mNSigmaElectron[1])) &&
	    ((track->NSigmaPion() <= mNSigmaPion[0]) ||
	     (track->NSigmaPion() >= mNSigmaPion[1])) &&
	    ((track->NSigmaKaon() <= mNSigmaKaon[0]) ||
	     (track->NSigmaKaon() >= mNSigmaKaon[1])) &&
	    ((track->NSigmaProton() <= mNSigmaProton[0]) ||
	     (track->NSigmaProton() >= mNSigmaProton[1])) )
	  mGoodToTPID = true;
	
#ifdef STHBTDEBUG
	cout << "ToT electron TPC PID check: " << mGoodTpcPID << endl;
#endif
      } //if((mPartId == 1) && mGoodTpcTrack)
      else if((mPartId == 3) && mGoodTpcTrack) {
	
	if( ((track->NSigmaPion() >= mNSigmaPion[0]) &&
	     (track->NSigmaPion() <= mNSigmaPion[1])) &&
	    ((track->NSigmaElectron() <= mNSigmaElectron[0]) ||
	     (track->NSigmaElectron() >= mNSigmaElectron[1])) &&
	    ((track->NSigmaKaon() <= mNSigmaKaon[0]) ||
	     (track->NSigmaKaon() >= mNSigmaKaon[1])) &&
	    ((track->NSigmaProton() <= mNSigmaProton[0]) ||
	     (track->NSigmaProton() >= mNSigmaProton[1])) )
	  mGoodToTPID = true;
      
#ifdef STHBTDEBUG
	cout << "ToT pion TPC PID check: " << mGoodTpcPID << endl;
#endif
      
      } //else if((mPartId == 3) && mGoodTpcTrack)
      else if((mPartId == 4) && mGoodTpcTrack) {

	if( ((track->NSigmaKaon() >= mNSigmaKaon[0]) &&
	     (track->NSigmaKaon() <= mNSigmaKaon[1])) &&
	    ((track->NSigmaElectron() <= mNSigmaElectron[0]) ||
	     (track->NSigmaElectron() >= mNSigmaElectron[1])) &&
	    ((track->NSigmaPion() <= mNSigmaPion[0]) ||
	     (track->NSigmaPion() >= mNSigmaPion[1])) &&
	    ((track->NSigmaProton() <= mNSigmaProton[0]) ||
	     (track->NSigmaProton() >= mNSigmaProton[1])) )
	  mGoodToTPID = true;
	
#ifdef STHBTDEBUG
	cout << "ToT kaon TPC PID check: " << mGoodTpcPID << endl;
#endif
   
      } //else if((mPartId == 4) && mGoodTpcTrack)
      else if((mPartId == 5) && mGoodTpcTrack) {
	
	if( ((track->NSigmaProton() >= mNSigmaProton[0]) &&
	     (track->NSigmaProton() <= mNSigmaProton[1])) &&
	    ((track->NSigmaElectron() <= mNSigmaElectron[0]) ||
	     (track->NSigmaElectron() >= mNSigmaElectron[1])) &&
	    ((track->NSigmaPion() <= mNSigmaPion[0]) ||
	     (track->NSigmaPion() >= mNSigmaPion[1])) &&
	    ((track->NSigmaKaon() <= mNSigmaKaon[0]) ||
	     (track->NSigmaKaon() >= mNSigmaKaon[1])) )
	  mGoodToTPID = true;

#ifdef STHBTDEBUG
	cout << "ToT proton TPC PID check: " << mGoodTpcPID << endl;
#endif
      } //else if((mPartId == 5) && mGoodTpcTrack)
    } //else to if(track->TofBeta() > 0)

#ifdef STHBTDEBUG
    cout << "ToT PID check: " << mGoodToTPID << endl;
#endif
  } //else if(mPIDSelection == 3)


  /***********************************************************/
  //If ToF signal, then ToF+TPC and TPC only if no ToF signal
  else if(mPIDSelection == 4) {

#ifdef STHBTDEBUG
    cout << "ToT+TPC PID selection" << endl;
    cout << "Selecting particle type: " << mPartId << endl;
#endif

    //If track has a ToF signal then we use ToF PID only
    if(track->TofBeta() > 0) {
      
      if( (track->P().mag() >= mTofTrackMom[0]) &&
	  (track->P().mag() <= mTofTrackMom[1]) &&
	  (track->Pt() >= mTofTrackPt[0]) &&
	  (track->Pt() <= mTofTrackPt[1]) &&
	  mInvBetaDifference >= mInverseBetaDifferce[0] &&
	  mInvBetaDifference <= mInverseBetaDifferce[1] )
	mGoodTofTrack = true;
      
      if( track->TofMassSqr() >= mTofMassSqr[0] &&
	  track->TofMassSqr() <= mTofMassSqr[1] &&
	  mGoodTofTrack)
	mGoodTofPID = true;

      if((mPartId == 1) && mGoodTofPID) {
	if(track->NSigmaElectron() >= mTofNSigmaElectron[0] &&
	   track->NSigmaElectron() <= mTofNSigmaElectron[1])
	  mGoodToTPID = true;
      }
      else if((mPartId == 3) && mGoodTofPID) {
	if(track->NSigmaPion() >= mTofNSigmaPion[0] &&
	   track->NSigmaPion() <= mTofNSigmaPion[1])
	  mGoodToTPID = true;
      }
      else if((mPartId == 4) && mGoodTofPID) {
	if(track->NSigmaKaon() >= mTofNSigmaKaon[0] &&
	   track->NSigmaKaon() <= mTofNSigmaKaon[1])
	  mGoodToTPID = true;
      }
      else if((mPartId == 5) && mGoodTofPID) {
	if(track->NSigmaProton() >= mTofNSigmaProton[0] &&
	   track->NSigmaProton() <= mTofNSigmaProton[1])
	  mGoodToTPID = true;
      }

#ifdef STHBTDEBUG
      cout << "TofCut passed: " << mGoodToTPID << endl;
#endif
    } //if(track->TofBeta() > 0)
    else {
      
      if( (track->P().mag() >= mTpcTrackMom[0]) &&
	  (track->P().mag() <= mTpcTrackMom[1]) &&
	  (track->Pt() >= mTpcTrackPt[0]) &&
	  (track->Pt() <= mTpcTrackPt[1]) )
	mGoodTpcTrack = true;


      if((mPartId == 1) && mGoodTpcTrack) {

	if( ((track->NSigmaElectron() >= mNSigmaElectron[0]) &&
	     (track->NSigmaElectron() <= mNSigmaElectron[1])) &&
	    ((track->NSigmaPion() <= mNSigmaPion[0]) ||
	     (track->NSigmaPion() >= mNSigmaPion[1])) &&
	    ((track->NSigmaKaon() <= mNSigmaKaon[0]) ||
	     (track->NSigmaKaon() >= mNSigmaKaon[1])) &&
	    ((track->NSigmaProton() <= mNSigmaProton[0]) ||
	     (track->NSigmaProton() >= mNSigmaProton[1])) )
	  mGoodToTPID = true;
	
#ifdef STHBTDEBUG
	cout << "ToT electron TPC PID check: " << mGoodTpcPID << endl;
#endif
      } //if((mPartId == 1) && mGoodTpcTrack)
      else if((mPartId == 3) && mGoodTpcTrack) {
	
	if( ((track->NSigmaPion() >= mNSigmaPion[0]) &&
	     (track->NSigmaPion() <= mNSigmaPion[1])) &&
	    ((track->NSigmaElectron() <= mNSigmaElectron[0]) ||
	     (track->NSigmaElectron() >= mNSigmaElectron[1])) &&
	    ((track->NSigmaKaon() <= mNSigmaKaon[0]) ||
	     (track->NSigmaKaon() >= mNSigmaKaon[1])) &&
	    ((track->NSigmaProton() <= mNSigmaProton[0]) ||
	     (track->NSigmaProton() >= mNSigmaProton[1])) )
	  mGoodToTPID = true;
      
#ifdef STHBTDEBUG
	cout << "ToT pion TPC PID check: " << mGoodTpcPID << endl;
#endif
      
      } //else if((mPartId == 3) && mGoodTpcTrack)
      else if((mPartId == 4) && mGoodTpcTrack) {

	if( ((track->NSigmaKaon() >= mNSigmaKaon[0]) &&
	     (track->NSigmaKaon() <= mNSigmaKaon[1])) &&
	    ((track->NSigmaElectron() <= mNSigmaElectron[0]) ||
	     (track->NSigmaElectron() >= mNSigmaElectron[1])) &&
	    ((track->NSigmaPion() <= mNSigmaPion[0]) ||
	     (track->NSigmaPion() >= mNSigmaPion[1])) &&
	    ((track->NSigmaProton() <= mNSigmaProton[0]) ||
	     (track->NSigmaProton() >= mNSigmaProton[1])) )
	  mGoodToTPID = true;
	
#ifdef STHBTDEBUG
	cout << "ToT kaon TPC PID check: " << mGoodTpcPID << endl;
#endif
   
      } //else if((mPartId == 4) && mGoodTpcTrack)
      else if((mPartId == 5) && mGoodTpcTrack) {
	
	if( ((track->NSigmaProton() >= mNSigmaProton[0]) &&
	     (track->NSigmaProton() <= mNSigmaProton[1])) &&
	    ((track->NSigmaElectron() <= mNSigmaElectron[0]) ||
	     (track->NSigmaElectron() >= mNSigmaElectron[1])) &&
	    ((track->NSigmaPion() <= mNSigmaPion[0]) ||
	     (track->NSigmaPion() >= mNSigmaPion[1])) &&
	    ((track->NSigmaKaon() <= mNSigmaKaon[0]) ||
	     (track->NSigmaKaon() >= mNSigmaKaon[1])) )
	  mGoodToTPID = true;

#ifdef STHBTDEBUG
	cout << "ToT proton TPC PID check: " << mGoodTpcPID << endl;
#endif
      } //else if((mPartId == 5) && mGoodTpcTrack)
    } //else to if(track->TofBeta() > 0)

#ifdef STHBTDEBUG
    cout << "ToT PID check: " << mGoodToTPID << endl;
#endif
  } //else if(mPIDSelection == 4)
  else {
    cout << "gregsTrackCut::Pass - If you see this message then you "
	 << "are in a big trouble with mPIDSelection" << endl;
  }

  if((mPIDSelection==0) && (mGoodTrack) && (mGoodTpcPID) && (mGoodCharge))
    mGoodPID = true;
  if((mPIDSelection==1) && (mGoodTrack) && (mGoodTofPID) && (mGoodCharge))
    mGoodPID = true;
  if((mPIDSelection==2) && (mGoodTrack) && (mGoodTnTPID) && (mGoodCharge))
    mGoodPID = true;
  if((mPIDSelection==3) && (mGoodTrack) && (mGoodToTPID) && (mGoodCharge))
    mGoodPID = true;
  if((mPIDSelection==4) && (mGoodTrack) && (mGoodToTPID) && (mGoodCharge))
    mGoodPID = true;
  

    
#ifdef STHBTDEBUG
  cout << " mGoodCharge: " << mGoodCharge << " mGoodTrack: " << mGoodTrack
       << " mGoodPID: " << mGoodPID << endl;
#endif

  (mGoodPID) ? mNTracksPassed++ : mNTracksFailed++;

#ifdef STHBTDEBUG
  if(mGoodPID)
    cout << "Single particle cut passed" << endl;
  else
    cout << "Single particle cut failed" << endl;
#endif

  return (mGoodPID);
}

//_________________
StHbtString gregsTrackCut::Report(){

  stringstream rep_stream;
  rep_stream << endl << "-- gregsTrackCut --" << endl
	     << "Charge: " << mCharge<< endl
	     << "mPartId: " << mPartId << endl
	     << mNSigmaPion[0] << " <= nSigma(pi) <= " << mNSigmaPion[1] << endl
	     << mNSigmaKaon[0] << " <= nSigma(K) <= " << mNSigmaKaon[1] << endl
	     << mNSigmaProton[0] << " <= nSigma(p) <= " << mNSigmaProton[1] << endl
	     << mTofNSigmaPion[0] << " <= ToF TPC nSigma(pi) <= " << mTofNSigmaPion[1] << endl
	     << mTofNSigmaKaon[0] << " <= ToF TPC nSigma(K) <= " << mTofNSigmaKaon[1] << endl
	     << mTofNSigmaProton[0] << " <= ToF TPC NSigma(p) <= " << mTofNSigmaProton[1] << endl
	     << mNHits[0] << " <= NHits <= " << mNHits[1] << endl
	     << mNHitsFit[0] << " <= NHitsFit <= " << mNHitsFit[1] << endl
	     << mNdEdxHits[0] << " <= NHits <= " << mNdEdxHits[1] << endl
	     << mPt[0] << " <= pT <= " << mPt[1] << endl
	     << mP[0] << " <= p <= " << mP[1] << endl
	     << mRapidity[0] << " <= rapidity <= " << mRapidity[1] << endl
	     << mEta[0] << " <= pseudorapidity <= " << mEta[1] << endl
	     << mDCA[0] << " <= dca <= " << mDCA[1] << endl
             << mTofMassSqr[0] << " <= ToF m^2 <= " << mTofMassSqr[1] << endl
	     << mInverseBetaDifferce[0] << " <= 1/beta-1/beta_hypo <= " << mInverseBetaDifferce[1] << endl
	     << "NHitsFit/NHitsPossible > " << mAntiSplit << endl
	     << " >>> N Tracks passed " << mNTracksPassed << endl
	     << " >>> N Tracks failed " << mNTracksFailed << endl;
  return (rep_stream.str());
}