#include "StFemtoDstInclusiveSelector.h"
#include "StBTofHeader.h"
#include "TBranch.h"
#include "StPhysicalHelixD.hh"

ClassImp(StFemtoDstInclusiveSelector)

// Set maximum file size to 1.9 GB (Root has a 2GB limit)
#define MAXFILESIZE 1900000000

//_________________
StFemtoDstInclusiveSelector::StFemtoDstInclusiveSelector(StMuDstMaker *muMaker,
                                                         const Char_t *oFileName) {

  std::cout << "StFemtoDstInclusiveSelector::Constructor - Creating an instance...";
  mOutFileName = oFileName;
  mMuDstMaker = muMaker;
  mMuDst = mMuDstMaker->muDst();
  mRefMultCorrUtil = NULL;
  mCompression = 9;

  mCollisionTypeAuAu = false; // pp: false; AuAu: true
  mLoopAllVert = false;       // Not first primary vertex does not contain tracks with fast detectors (TOF)

  mEventIsGood = false;
  mNEventsIn = 0;
  mNEventsPassed = 0;

  mAuAuCorrectZdc = true;

  mIsGoodTrack = false;
  mCurrentTrigger = 0;

  //Initialize event cut variables
  mPrimVtxZ[0] = -70.;
  mPrimVtxZ[1] = 70.;
  mPrimVtxR[0] = -1.;
  mPrimVtxR[1] = 10.;
  mPrimVtxVpdVzDiff[0] = -5.;
  mPrimVtxVpdVzDiff[1]= 5.;
  mPrimVtxXShift = 0.;
  mPrimVtxYShift = 0.;

  //Initialize single-particle cut variables
  mTrackMomentum[0] = 0.1;
  mTrackMomentum[1] = 1.7;
  mTrackDca[0] = 0.;
  mTrackDca[1] = 5.;
  mTrackDcaGlobal[0] = 0.;
  mTrackDcaGlobal[1] = 5.;
  mTrackNHits[0] = 14;
  mTrackNHits[1] = 50;

  mTrackEta[0] = -1.;
  mTrackEta[1] = 1.;
  mTrackFlag[0] = 0;
  mTrackFlag[1] = 1000;

  mSelectPions = true;
  mSelectTpcPions = true;
  mPionTpcMomentum[0] = 0.15;      mPionTpcMomentum[1] = 0.65;
  mPionTpcNSigmaElectron[0]= -2.; mPionTpcNSigmaElectron[1]= 2.;
  mPionTpcNSigmaPion[0]= -3.;      mPionTpcNSigmaPion[1]= 3.;
  mPionTpcNSigmaKaon[0]= -1.5;     mPionTpcNSigmaKaon[1]= 1.5;
  mPionTpcNSigmaProton[0]= -3.;   mPionTpcNSigmaProton[1]= 3.;
  mSelectTofPions = true;
  mPionTofMomentum[0] = 0.15;      mPionTofMomentum[1] = 1.75;
  mPionMassSqr[0] = -0.01;         mPionMassSqr[1] = 0.09;

  mSelectKaons = true;
  mSelectTpcKaons = true;
  mKaonTpcMomentum[0] = 0.15;      mKaonTpcMomentum[1] = 0.65;
  mKaonTpcNSigmaElectron[0]= -2.; mKaonTpcNSigmaElectron[1]= 2.;
  mKaonTpcNSigmaPion[0]= -1.5;     mKaonTpcNSigmaPion[1]= 1.5;
  mKaonTpcNSigmaKaon[0]= -3.;      mKaonTpcNSigmaKaon[1]= 3.;
  mKaonTpcNSigmaProton[0]= -3.;   mKaonTpcNSigmaProton[1]= 3.;
  mSelectTofKaons = true;
  mKaonTofMomentum[0] = 0.15;      mKaonTofMomentum[1] = 1.75;
  mKaonMassSqr[0] = 0.16;          mKaonMassSqr[1] = 0.35;

  mSelectProtons = false;
  mSelectTpcProtons = true;
  mProtonTpcMomentum[0] = 0.15;      mProtonTpcMomentum[1] = 1.2;
  mProtonTpcNSigmaElectron[0]= -3.; mProtonTpcNSigmaElectron[1]= 3.;
  mProtonTpcNSigmaPion[0]= -3.;     mProtonTpcNSigmaPion[1]= 3.;
  mProtonTpcNSigmaKaon[0]= -3.;     mProtonTpcNSigmaKaon[1]= 3.;
  mProtonTpcNSigmaProton[0]= -3.;    mProtonTpcNSigmaProton[1]= 3.;
  mSelectTofProtons = true;
  mProtonTofMomentum[0] = 0.15; mProtonTofMomentum[1] = 1.75;
  mProtonMassSqr[0] = 0.65; mProtonMassSqr[1] = 1.05;

  std::cout << "\t[DONE]" << std::endl;
}

//_________________
StFemtoDstInclusiveSelector::~StFemtoDstInclusiveSelector() {
  /* do nothing */
}

//_________________
Int_t StFemtoDstInclusiveSelector::Init() {

  std::cout << "StFemtoDstInclusiveSelector::Init - Initializing the maker"
      << std::endl 
      << "Creating output file: " << mOutFileName;

  mFemtoEvent = new StFemtoEvent();
  mOutFile = new TFile(mOutFileName, "RECREATE");
  mOutFile->SetCompressionLevel(mCompression);
  std::cout << "\t[DONE]" << std::endl;

  mTree = new TTree("StFemtoDst","StFemtoDst");
  mTree->SetMaxTreeSize(MAXFILESIZE);
  //mTree->SetAutoSave(1000000);
  mTree->Branch("StFemtoEvent","StFemtoEvent", &mFemtoEvent);

  if(mCollisionTypeAuAu) { //For AuAu only
    if(!mRefMultCorrUtil) {
      std::cout << "StFemtoDstInclusiveSelector::Init - Initializing StRefMultCorr...";
      mRefMultCorrUtil = new StRefMultCorr("refmult");
      std::cout << "\t[DONE]" << endl;
    }
  }

  mNBytes = 0;
  std::cout << "StFemtoDstInclusiveSelector::Init - Initialization has been finished"
      << std::endl;
  return StMaker::Init();
}

//_________________
void StFemtoDstInclusiveSelector::Clear(Option_t *option) {
  StMaker::Clear(option);
}

//_________________
Int_t StFemtoDstInclusiveSelector::Make() {

  // Get MuEvent from the MuDst
  mNEventsIn++;
  StMuEvent *mMuEvent = mMuDst->event();
  if (!AcceptTrigger(mMuEvent))  { // If trigger the trigger was not found
    return kStOk;
  }
  if (mMuEvent->refMult() < 0) { // Multiplicity can't be negative
    return kStOk;
  }

  UShort_t mNVertices = mLoopAllVert ? mMuDst->numberOfPrimaryVertices() : 1;
  Short_t mNPrimTracks = mMuDst->numberOfPrimaryTracks();
  Short_t mNGlobTracks = mMuDst->numberOfGlobalTracks();

  //Reasonable amount of tracks
  if (mNPrimTracks < 0 || mNPrimTracks > 10000) {
    return kStOk;
  }

  Bool_t mEventIsGood = false;
  Float_t mVpdVz = 0.;
  StBTofHeader *mTofHeader = mMuDst->btofHeader();

  // Centrality variables
  Int_t mCent16;
  Double_t mCentWeight, mRefMultCorrVal;

  // If pp collisions
  if (!mCollisionTypeAuAu) {
    mCent16 = -1;
    mCentWeight = -999.;
  }

  Float_t mRanking = 0;
  StMuPrimaryVertex *mPrimVertex = NULL;
  StMuTrack *mPrimTrack = NULL;
  StMuTrack *mGlobTrack = NULL;
  Float_t mMassSqr = -999.;
  Float_t mBeta = -999.;
  Bool_t mIsTofTrack = false;
  Bool_t mIsTpcPion = false;
  Bool_t mIsTofPion = false;
  Bool_t mIsTpcKaon = false;
  Bool_t mIsTofKaon = false;
  Bool_t mIsTpcProton = false;
  Bool_t mIsTofProton = false;

  //Calculate the amount of the pile-up vertices
  UShort_t mNVtxInBunchCross = 0;
  for(UShort_t iVert=0; iVert<mMuDst->numberOfPrimaryVertices(); iVert++) {

    mPrimVertex = mMuDst->primaryVertex(iVert);
    if(!mPrimVertex || !mTofHeader) continue;
    if(mPrimVertex->ranking() <= 0.) continue;

    if(TMath::Abs(mPrimVertex->position().x()) < 1e-5 &&
       TMath::Abs(mPrimVertex->position().y()) < 1e-5 &&
       TMath::Abs(mPrimVertex->position().z()) < 1e-5) {
      continue;
    }
    if(TMath::Abs( mPrimVertex->position().z() - mTofHeader->vpdVz()) < 5. ) {
      mNVtxInBunchCross++;
    }
  } //for(unsigned short iVert=0; iVert<mNVertices; iVert++)

  //There should be only one primary vertex that caused the trigger
  if(mNVtxInBunchCross > 1) {
    return kStOk;
  }

  // Vertex loop
  for (UShort_t iVert = 0; iVert < mNVertices; iVert++) {

    mEventIsGood = false;
    mPrimVertex = mMuDst->primaryVertex(iVert);
    if (!mPrimVertex) {
      //std::cout << "No primary vertices in the event. Skip..." << std::endl;
      continue;
    }

    mRanking = mPrimVertex->ranking();
    //Positive ranking only
    if (mPrimVertex->ranking() <= 0) continue;

    //Not (0,0,0) position of the primary vertex
    if (TMath::Abs(mPrimVertex->position().x()) < 1e-5 &&
        TMath::Abs(mPrimVertex->position().y()) < 1e-5 &&
        TMath::Abs(mPrimVertex->position().z()) < 1e-5 ) continue;

    StThreeVectorF mVertPosition = mPrimVertex->position();    
    if(mTofHeader) { //In case of the absence of the vpd information
      mVpdVz = mTofHeader->vpdVz();
    }
    else {
      mVpdVz = mVertPosition.z();
    }

    // Event cut
    if (!AcceptPrimVtx(mVertPosition, mVpdVz)) continue;
    mEventIsGood = true;

    // Set StRefMultCorr
    if (mCollisionTypeAuAu) { // Only for AuAu collisions
      mRefMultCorrUtil->init(mMuEvent->runId());
      //Check for the bad run
      if (mRefMultCorrUtil->isBadRun(mMuEvent->runId())) {
        std::cout << "StRefMultCorrUtil::IsBadRun" << std::endl;
        continue;
      }
      if (mAuAuCorrectZdc == true) {
        mRefMultCorrUtil->initEvent(mMuEvent->refMult(),
                                    mVertPosition.z(),
                                    mMuEvent->runInfo().zdcCoincidenceRate());
      }
      else {
        mRefMultCorrUtil->initEvent(mMuEvent->refMult(),
                                    mVertPosition.z());
      }
      mCent16 = mRefMultCorrUtil->getCentralityBin16();
      if (mCent16 < 0) continue;   //Some AuAu collisions have -1 value
      mCentWeight = mRefMultCorrUtil->getWeight();
      mRefMultCorrVal = mRefMultCorrUtil->getRefMultCorr();
    }

    // Write event information to the FemtoDst
    mNEventsPassed++;
    mFemtoEvent->SetEventId( mMuEvent->eventId() );
    mFemtoEvent->SetRunId( mMuEvent->runId() );
    mFemtoEvent->SetCollisionType( mCollisionTypeAuAu );
    mFemtoEvent->SetRefMult( mMuEvent->refMult() );
    mFemtoEvent->SetRefMultCorr( mRefMultCorrVal );
    mFemtoEvent->SetRefMultCorrWeight( mCentWeight );
    mFemtoEvent->SetCent16( mCent16 );
    mFemtoEvent->SetRefMultPos( mMuEvent->refMultPos() );
    mFemtoEvent->SetNumberOfBTofHit( mMuDst->numberOfBTofHit() ); //mMuEvent->btofTrayMultiplicity() );
    //mFemtoEvent->SetZDCe( mMuEvent->zdcTriggerDetector().adc(4) );
    //mFemtoEvent->SetZDCw( mMuEvent->zdcTriggerDetector().adc(0) );
    mFemtoEvent->SetNumberOfPrimaryTracks( mNPrimTracks );
    mFemtoEvent->SetNumberOfGlobalTracks( mNGlobTracks );
    mFemtoEvent->SetMagneticField( mMuEvent->magneticField() );
    mFemtoEvent->SetPrimaryVertexPosition( mVertPosition.x(),
                                          mVertPosition.y(),
                                          mVertPosition.z() );
    mFemtoEvent->SetVpdVz(mVpdVz);
    mFemtoEvent->SetTriggerId(mCurrentTrigger);
    mFemtoEvent->SetPrimaryVertexRanking(mRanking);

    UInt_t refMult2 = 0;
    UInt_t refMult2Pos = 0;

    //Primary track loop
    for(UShort_t iTrk = 0; iTrk < mNPrimTracks; iTrk++) {

      mPrimTrack = mMuDst->primaryTracks(iTrk);
      mGlobTrack = mMuDst->globalTracks(mPrimTrack->index2Global());

      //RefMult2 calculation
      if ( AcceptRefMult2(mPrimTrack) ) {
        refMult2++;

        if(mPrimTrack->charge() > 0) refMult2Pos++;
      }

      //Track cut
      if ( !AcceptTrack(mPrimTrack, iVert) ) continue;

      mIsTofTrack = false;
      mIsTofTrack = IsTofTrack(mGlobTrack);
      mIsTpcPion = false;
      mIsTofPion = false;
      mIsTpcKaon = false;
      mIsTofKaon = false;
      mIsTpcProton = false;
      mIsTofProton = false;

      mBeta = -999.;
      mMassSqr = -999.;

      // Check if track has ToF signal
      if(mIsTofTrack) {
        mBeta = mGlobTrack->btofPidTraits().beta();
        mMassSqr = mPrimTrack->momentum().mag2() * ( 1./(mBeta * mBeta) - 1.);
      }

      // Pion selection
      if (mSelectPions) {
        if (mSelectTpcPions) 
          mIsTpcPion = IsTpcPion(mPrimTrack);
        if (mSelectTofPions && mIsTofTrack)
          mIsTofPion = IsTofPion(mMassSqr, mGlobTrack->p().mag());
      }

      // Kaon selection
      if (mSelectKaons) {
        if (mSelectTpcKaons)
          mIsTpcKaon = IsTpcKaon(mPrimTrack);
        if (mSelectTofKaons && mIsTofTrack)
          mIsTofKaon = IsTofKaon(mMassSqr, mGlobTrack->p().mag());
      }

      // Proton selection
      if (mSelectProtons) {
        if (mSelectTpcProtons)
          mIsTpcProton = IsTpcProton(mPrimTrack);
        if (mSelectTofProtons && mIsTofTrack)
          mIsTofProton = IsTofProton(mMassSqr, mGlobTrack->p().mag());
      }

      // Write StFemtoTrack
      if (mIsTpcPion || mIsTofPion || mIsTpcKaon || mIsTofKaon || 
          mIsTpcProton || mIsTofProton) {

        mFTrack = mFemtoEvent->AddFemtoTrack();

        mFTrack->SetId( mPrimTrack->id() );
        mFTrack->SetNHits( (Char_t)(mPrimTrack->charge()*mPrimTrack->nHits()) );
        mFTrack->SetNHitsPoss( mPrimTrack->nHitsPoss(kTpcId) );
        //mFTrack->SetNHitsDedx( mPrimTrack->nHitsDedx() );
        mFTrack->SetNSigmaElectron( mPrimTrack->nSigmaElectron() );
        mFTrack->SetNSigmaPion( mPrimTrack->nSigmaPion() );
        mFTrack->SetNSigmaKaon( mPrimTrack->nSigmaKaon() );
        mFTrack->SetNSigmaProton( mPrimTrack->nSigmaProton() );
        mFTrack->SetDedx( mPrimTrack->dEdx() );
        //mFTrack->SetChi2( mPrimTrack->chi2() );
        mFTrack->SetMap0( mGlobTrack->topologyMap().data(0) );
        mFTrack->SetMap1( mGlobTrack->topologyMap().data(1) );
        mFTrack->SetP( mPrimTrack->momentum().x(), 
                      mPrimTrack->momentum().y(), 
                      mPrimTrack->momentum().z() );
        mFTrack->SetDCAxGlobal( mPrimTrack->dcaGlobal().x() );
        mFTrack->SetDCAyGlobal( mPrimTrack->dcaGlobal().y() );
        mFTrack->SetDCAzGlobal( mPrimTrack->dcaGlobal().z() );
        mFTrack->SetGlobalP( mGlobTrack->p().x(),
                            mGlobTrack->p().y(),
                            mGlobTrack->p().z() );

        if(mIsTofTrack) {
          mFTrack->SetBeta( mBeta );
        }
        else {
          mFTrack->SetBeta( -999. );
        }
      }
    } //for(Int_t iTrk = 0; iTrk < mNPrimTracks; iTrk++)

    //Store RefMult2
    mFemtoEvent->SetRefMult2(refMult2);
    mFemtoEvent->SetRefMult2Pos(refMult2Pos);

    //
    // If event is good then write it
    //
    if(mEventIsGood) {
      mNBytes += mTree->Fill();
      mFemtoEvent->Clear();
    }
  } //(Int_t iVert = 0; iVert < mNVertices; iVert++)

  return kStOk;
}

//_________________
Bool_t StFemtoDstInclusiveSelector::AcceptRefMult2(StMuTrack *trk) {

  Bool_t mIsGoodTrack = false;
  if ( trk->flag()> 0 && 
      trk->charge() != 0 && 
      trk->nHitsFit() >= 10 &&
      trk->dca().mag() < 3.0 &&
      trk->momentum().mag() >= 1e-10 &&
      TMath::Abs(trk->eta()) <= 1.0) {

    mIsGoodTrack = true;
  }

  return mIsGoodTrack;
}

//_________________
Bool_t StFemtoDstInclusiveSelector::AcceptTrigger(StMuEvent *muEvent) {
  //We will store the bit mask for the list of triggers.
  //The list of triggers should be the same and in the same
  //order, in order to read it later
  Bool_t mIsGoodTrigger = false;
  mCurrentTrigger = 0;
  if(mTriggerIdCollection.empty()) {
    mIsGoodTrigger = true;
  }
  else {
    for (UInt_t iTrg = 0; iTrg < mTriggerIdCollection.size(); iTrg++) {
      if(muEvent->triggerIdCollection().nominal().isTrigger(mTriggerIdCollection.at(iTrg))) {
        mIsGoodTrigger = true;
        mCurrentTrigger |= (1<<iTrg);
        break;
      }
    } //for(UInt_t iTrg=0; iTrg<mTriggerIdCollection.size(); iTrg++)
  }

  return mIsGoodTrigger;
}

//_________________
Bool_t StFemtoDstInclusiveSelector::AcceptPrimVtx(StThreeVectorF vtxPos, 
                                                  Float_t vpdVz) {
  Bool_t mIsGoodVtx = false;
  Float_t mVtxX = vtxPos.x() - mPrimVtxXShift;
  Float_t mVtxY = vtxPos.y() - mPrimVtxYShift;
  Float_t vtxR = TMath::Sqrt(mVtxX*mVtxX + mVtxY*mVtxY);
  Float_t vpdDiff = vtxPos.z() - vpdVz;

  if( vtxPos.z() >= mPrimVtxZ[0] &&
     vtxPos.z() <= mPrimVtxZ[1] &&
     vtxR >= mPrimVtxR[0] &&
     vtxR <= mPrimVtxR[1] &&
     vpdDiff >= mPrimVtxVpdVzDiff[0] &&
     vpdDiff <= mPrimVtxVpdVzDiff[1]) {
    mIsGoodVtx = true;
  }

  return mIsGoodVtx;
}

//_________________
Bool_t StFemtoDstInclusiveSelector::AcceptTrack(StMuTrack *trk, UShort_t vtxInd) {

  Bool_t mIsGoodTrack = false;
  if( trk->vertexIndex() == vtxInd &&
     trk->momentum().mag() >= mTrackMomentum[0] && trk->momentum().mag() <= mTrackMomentum[1] &&
     trk->type() == 1 &&
     trk->eta() >= mTrackEta[0] && trk->eta() <= mTrackEta[1] &&
     trk->nHits() >= mTrackNHits[0] && trk->nHits() <= mTrackNHits[1] &&
     trk->flag() >= mTrackFlag[0] && trk->flag() <= mTrackFlag[1] &&
     trk->dca(vtxInd).perp() >= mTrackDca[0] && trk->dca(vtxInd).perp() <= mTrackDca[1] &&
     trk->dcaGlobal(vtxInd).perp() >= mTrackDcaGlobal[0] &&
     trk->dcaGlobal(vtxInd).perp() <= mTrackDcaGlobal[1]) {

    mIsGoodTrack = true;
  }

  return mIsGoodTrack;
}

//_________________
Bool_t StFemtoDstInclusiveSelector::IsTofTrack(StMuTrack *trk) { //Only for globals

  Bool_t mIsTofHit = false;
  if(trk->btofPidTraits().beta() > 0 &&
     trk->btofPidTraits().timeOfFlight() > 0) {
    mIsTofHit = true;
  }
  return mIsTofHit;
}

//_________________
Int_t StFemtoDstInclusiveSelector::Finish() {

  mOutFile->cd();
  mOutFile->Write();
  //mTree->Write();
  mOutFile->Close();
  delete mTree;
  delete mOutFile;

  std::cout << "*********************************************" << std::endl
      << "StFemtoDstInclusiveSelector has been finished" << std::endl
      << "\t nEventsPassed   : " << mNEventsPassed  << std::endl
      << "\t nEventsProcessed: " << mNEventsIn << std::endl
      << "**********************************************" << std::endl;

  return kStOk;
}

//_________________
Bool_t StFemtoDstInclusiveSelector::IsTpcPion(StMuTrack *pTrk) {
  Float_t p = pTrk->p().mag();
  Float_t nSElectron = pTrk->nSigmaElectron();
  Float_t nSPion = pTrk->nSigmaPion();
  Float_t nSKaon = pTrk->nSigmaKaon();
  Float_t nSProton = pTrk->nSigmaProton();

  return ( (p >= mPionTpcMomentum[0]) && (p <= mPionTpcMomentum[1]) &&
          ( (nSPion >= mPionTpcNSigmaPion[0]) &&
           (nSPion <= mPionTpcNSigmaPion[1]) ) &&
          ( (nSElectron <= mPionTpcNSigmaElectron[0]) ||
           (nSElectron >= mPionTpcNSigmaElectron[1]))  &&
          ( (nSKaon <= mPionTpcNSigmaKaon[0]) ||
           (nSKaon >= mPionTpcNSigmaKaon[1]) ) &&
          ( (nSProton <= mPionTpcNSigmaProton[0]) ||
           (nSProton >= mPionTpcNSigmaProton[1]) ) );
}

//_________________
Bool_t StFemtoDstInclusiveSelector::IsTpcKaon(StMuTrack *pTrk) {
  Float_t p = pTrk->p().mag();
  Float_t nSElectron = pTrk->nSigmaElectron();
  Float_t nSPion = pTrk->nSigmaPion();
  Float_t nSKaon = pTrk->nSigmaKaon();
  Float_t nSProton = pTrk->nSigmaProton();

  return ( (p >= mKaonTpcMomentum[0]) && (p <= mKaonTpcMomentum[1]) &&
          ( (nSKaon >= mKaonTpcNSigmaKaon[0]) &&
           (nSKaon <= mKaonTpcNSigmaKaon[1]) ) &&
          ( (nSElectron <= mKaonTpcNSigmaElectron[0]) ||
           (nSElectron >= mKaonTpcNSigmaElectron[1]) ) &&
          ( (nSPion <= mKaonTpcNSigmaPion[0]) ||
           (nSPion >= mKaonTpcNSigmaPion[1]) ) &&
          ( (nSProton <= mKaonTpcNSigmaProton[0]) ||
           (nSProton >= mKaonTpcNSigmaProton[1]) ) );
}

//_________________
Bool_t StFemtoDstInclusiveSelector::IsTpcProton(StMuTrack *pTrk) {
  Float_t p = pTrk->p().mag();
  Float_t nSElectron = pTrk->nSigmaElectron();
  Float_t nSPion = pTrk->nSigmaPion();
  Float_t nSKaon = pTrk->nSigmaKaon();
  Float_t nSProton = pTrk->nSigmaProton();

  return ( (p >= mProtonTpcMomentum[0]) && (p <= mProtonTpcMomentum[1]) &&
          ( (nSProton >= mProtonTpcNSigmaProton[0]) &&
           (nSProton <= mProtonTpcNSigmaProton[1]) ) &&
          ( (nSElectron <= mProtonTpcNSigmaElectron[0]) ||
           (nSElectron >= mProtonTpcNSigmaElectron[1]) ) &&
          ( (nSPion <= mProtonTpcNSigmaPion[0]) ||
           (nSPion >= mProtonTpcNSigmaPion[1]) ) &&
          ( (nSKaon <= mProtonTpcNSigmaKaon[0]) ||
           (nSKaon >= mProtonTpcNSigmaKaon[1]) ) );
}

//_________________
Bool_t StFemtoDstInclusiveSelector::IsTofPion(Float_t mSqr, Float_t p) {

  return ( (p >= mPionTofMomentum[0]) && (p <= mPionTofMomentum[1]) &&
          (mSqr >= mPionMassSqr[0]) && (mSqr <= mPionMassSqr[1]) );
}

//_________________
Bool_t StFemtoDstInclusiveSelector::IsTofKaon(Float_t mSqr, Float_t p) {

  return ( (p >= mKaonTofMomentum[0]) && (p <= mKaonTofMomentum[1]) &&
          (mSqr >= mKaonMassSqr[0]) && (mSqr <= mKaonMassSqr[1]) );
}

//_________________
Bool_t StFemtoDstInclusiveSelector::IsTofProton(Float_t mSqr, Float_t p) {

  return ( (p >= mProtonTofMomentum[0]) && (p <= mProtonTofMomentum[1]) &&
          (mSqr >= mProtonMassSqr[0]) && (mSqr <= mProtonMassSqr[1]) );
}

//_________________
void StFemtoDstInclusiveSelector::AddTriggerId(const UInt_t& id) {
  mTriggerIdCollection.push_back(id);
}