my_star/StRoot/StFemtoDstMaker/StFemtoDstMaker.cxx

#include "StFemtoDstMaker.h"
#include "StBTofHeader.h"
#include "TBranch.h"
#include "StPhysicalHelixD.hh"
#include "StDcaGeometry.h"

ClassImp(StFemtoDstMaker)

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

//_________________
StFemtoDstMaker::StFemtoDstMaker(StMuDstMaker *muMaker,
                                 const char *calibPath,                                 
                                 const Char_t *oFileName) {

  std::cout << "StFemtoDstMaker::StFemtoDstMaker - Creating an instance...";

  mAcceptTriggerFail = 0;
  mRefMultFail = 0;
  mAntiPileupFail = 0;
  mRankingFail = 0;
  mVtxNullFail = 0;
  mNPrimVtxFail = 0;
  mAcceptPrimVtxFail = 0;
  matrix = new TMatrixTSym<double>(2);
  mPtCut = 0.1;

  mOutFileName = oFileName;
  mMuDstMaker = muMaker;
  mRefMultCorrUtil = NULL;
  mCompression = 9;

  mCollisionType = false;     // 0-pp, 1-AuAu

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

  mAuAuCorrectZdc = true;

  mIsGoodTrack = false;
  mCurrentTrigger = 0;

  mCalibPath = calibPath;

  //
  // Default cut values
  //

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

  // initialize single-particle cut variables
  mTrackMomentum[0] = 0.1;
  mTrackMomentum[1] = 2.;
  mTrackDca[0] = 0.;
  mTrackDca[1] = 5.;
  mTrackDcaGlobal[0] = 0.;
  mTrackDcaGlobal[1] = 5.;
  mTrackNHits[0] = 15;
  mTrackNHits[1] = 50;
  mTrackNHitsFit[0] = 15;
  mTrackNHitsFit[1] = 50;
  mTrackEta[0] = -1.1;
  mTrackEta[1] = 1.1;
  mTrackFlag[0] = 0;
  mTrackFlag[1] = 1000;

  // flags and cuts for exclusion
  mRemovePions = false;
  mRemoveKaons = false;
  mRemoveProtons = false;
  mPthresh   = 0.4;
  mTpcPionNSigma[0] = -1.5;
  mTpcPionNSigma[1] = 1.5;
  mTofPionMassSqr[0] = -0.1;
  mTofPionMassSqr[1] = 0.15;
  mTpcKaonNSigma[0] = -1.5;
  mTpcKaonNSigma[1] = 1.5;
  mTofKaonMassSqr[0] = 0.2;
  mTofKaonMassSqr[1] = 0.35;
  mTpcProtonNSigma[0] = -1.5;
  mTpcProtonNSigma[1] = 1.5;
  mTofProtonMassSqr[0] = 0.7;
  mTofProtonMassSqr[1] = 1.1;

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

//_________________
StFemtoDstMaker::~StFemtoDstMaker() {
  delete mQVMaker;
  delete matrix;
}

//_________________
Int_t StFemtoDstMaker::Init() {

  std::cout << "StFemtoDstMaker::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(mCollisionType == true) { // for AuAu only
    if(!mRefMultCorrUtil) {
      std::cout << "StFemtoDstMaker::Init - Initializing StRefMultCorr...";
      mRefMultCorrUtil = new StRefMultCorr("refmult");
      std::cout << "\t[DONE]" << endl;
    }

    // initialize reaction plane maker
    mQVMaker = new QVMaker(4 /* do not fill hEP histograms just calculate Psi */);
  }

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

//________________
void StFemtoDstMaker::Clear(Option_t *option) {
  StMaker::Clear();
}

//________________
Int_t StFemtoDstMaker::Make() {

  mNEventsIn++;

  mMuDst = NULL;
  mMuEvent = NULL;

  //
  // MuDstMaker initialization
  //
  if(!mMuDstMaker) {
    mMuDstMaker = (StMuDstMaker*)GetMaker("MuDst");
    if (!mMuDstMaker) {
      LOG_ERROR << "StFemtoDstMaker::Make [ERROR] - Cannot find StMuDstMaker"
                << std::endl;
      return kStOk;
    }
  }

  //
  // Obtaining MuDst
  //
  mMuDst = mMuDstMaker->muDst();
  if (!mMuDst) {
    mMuDst = (StMuDst*)GetInputDS("MuDst");
    if (!mMuDst) {
      gMessMgr->Warning() << "StFemtoDstMaker::Make [WARNING] - No MuDst has been found" 
                          << endm;
      return kStOk;
    }
  }

  //
  // obtaining MuEvent
  //
  mMuEvent = (StMuEvent*)mMuDst->event();

  if(mMuEvent->refMult() < 0) { // multiplicity cannot be negative
    mRefMultFail++;
    return kStOk;
  }

  Int_t mNVertices = mMuDst->numberOfPrimaryVertices();
  Int_t mNPrimTracks = mMuDst->numberOfPrimaryTracks();
  Int_t mNGlobTracks = mMuDst->numberOfGlobalTracks();

  int iVert = 0;
  mPrimVertex = mMuDst->primaryVertex(iVert);
  StBTofHeader* mTofHeader = mMuDst->btofHeader();

  if (!mTofHeader || !mPrimVertex)   return kStOk;

  if (mNVertices == 0)   return kStOk;

  if (mPrimVertex->ranking() <= 0) // positive ranking only
  {
    mRankingFail++;
    return kStOk;
  }

  //
  // Check that position of the primary vertex if not (0,0,0)
  //
  if(mPrimVertex->position().x() == 0 &&
     mPrimVertex->position().y() == 0 &&
     mPrimVertex->position().z() == 0)
  {
    mVtxNullFail++;
    return kStOk;
  }

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

  //
  // Some initializations of local variables
  //
  StThreeVectorF mVertPosition;
  Float_t mVpdVz = 0.;
  Int_t mCent16;
  Double_t mCentWeight;
  Double_t mRefMultCorrVal;
  Double_t mBeta;
  Double_t mMassSqr;
  Int_t   mPrimVertIndex = -999;
  Float_t mRanking = -999.;
  Bool_t mGoodRPevent = false;
  

  if(mCollisionType == false) { // pp collisions means false state
    mCent16 = -1;
    mCentWeight = -999.;
  }

  //Calculate the amount of the pile-up vertices
  int mNVtxInBunchCross = 0;
  for (unsigned int i = 0; i < mMuDst->numberOfPrimaryVertices(); i++)
  {
    StMuPrimaryVertex *primVtx = mMuDst->primaryVertex(i);
    if(!primVtx || primVtx->ranking() <= 0.)   continue;

    if(TMath::Abs(primVtx->position().x()) < 1e-5 &&
       TMath::Abs(primVtx->position().y()) < 1e-5 &&
       TMath::Abs(primVtx->position().z()) < 1e-5)     continue;

    if(TMath::Abs( primVtx->position().z() - mTofHeader->vpdVz()) < 5. )   mNVtxInBunchCross++;
  }

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

  mEventIsGood = false;

  mPrimVertIndex = iVert;
  mRanking = mPrimVertex->ranking();
  mVertPosition = mPrimVertex->position();
  mVpdVz = mTofHeader->vpdVz();

  if(!AcceptTrigger(mMuEvent)) // if trigger is not found
  {
    mAcceptTriggerFail++;
    return kStOk;
  }

  if(!AcceptPrimVtx(mVertPosition, mVpdVz))
  {
    mAcceptPrimVtxFail++;
    return kStOk;
  }

  mEventIsGood = true;
  //
  // Set StRefMultCorr
  //
  if(mCollisionType == true) // AuAu collisions means true state
  {
    int mass1 = mMuEvent->runInfo().beamMassNumber(StBeamDirection::east);
    int mass2 = mMuEvent->runInfo().beamMassNumber(StBeamDirection::west);
    
    // if Au+Au
    if (mass1 == 197 && mass2 == 197) {

      mRefMultCorrUtil->init(mMuEvent->runId());
      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) // some AuAu collisions have -1 value
        return kStOk;
      mCentWeight = mRefMultCorrUtil->getWeight();
      mRefMultCorrVal = mRefMultCorrUtil->getRefMultCorr();
    }
    else { // Assume that it is Cu+Au case
      mCent16 = GetCentralityBin(mMuEvent->refMult(), mVertPosition.z());
    }

    if (mCalibPath && LoadEPCalibParam(mMuEvent->runId()))
      mGoodRPevent = true;
  } //if(mCollisionType == true)

  mNEventsPassed++;
  mFemtoEvent->SetEventId(mMuEvent->eventId());
  mFemtoEvent->SetRunId(mMuEvent->runId());
  mFemtoEvent->SetCollisionType(mCollisionType);
  mFemtoEvent->SetRefMult(mMuEvent->refMult());
  mFemtoEvent->SetRefMultCorr(mRefMultCorrVal);
  mFemtoEvent->SetRefMultCorrWeight(mCentWeight);
  mFemtoEvent->SetRefMultPos(mMuEvent->refMultPos());
  mFemtoEvent->SetCent16(mCent16);
  //mFemtoEvent->SetZDCe(mMuEvent->zdcTriggerDetector().adc(4));
  //mFemtoEvent->SetZDCw(mMuEvent->zdcTriggerDetector().adc(0));
  mFemtoEvent->SetNumberOfBTofHit(mMuDst->numberOfBTofHit());
  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;

  StFemtoTrack *mFTrack = NULL;

  //
  // Calculate transverse sphericity
  //
  float sph          = -999.;
  float pt_full      = 0.0;

  matrix->Zero();
  for (int i = 0; i < mNPrimTracks; i++)
  {
    mPrimTrack = mMuDst->primaryTracks(i);
    if (mPrimTrack->vertexIndex() != 0)   return kStOk;

    float pt = mPrimTrack->pt();
    const StThreeVectorF &p = mPrimTrack->p();
    if (fabs(mPrimTrack->eta()) < 1.0 && pt > mPtCut)
    {
      float px = p.x();
      float py = p.y();

      (*matrix)(0, 0) += px*px/pt;
      (*matrix)(1, 1) += py*py/pt;
      (*matrix)(0, 1) += px*py/pt;
      (*matrix)(1, 0) += px*py/pt;

      pt_full += pt;
    }
  }

  *matrix *= 1./pt_full;
  TMatrixDSymEigen death_machine(*matrix);
  TVectorD eigen = death_machine.GetEigenValues();
  sph = 2.*eigen.Min()/(eigen[0] + eigen[1]);
  mFemtoEvent->SetSphericity(sph);

  //
  // Calculate reaction plane
  //
  float rpE = -999., rpW = -999.;
  if (mGoodRPevent) {
    int qZVtx = int((mPrimVertex->position().z() + 30.)/60. * NqvZvtx /* see ConstVar.h */);

    if (mCent16 >= 0 && mCent16 < NqvCent &&
        qZVtx >= 0 && qZVtx < NqvZvtx) {
      
      // SMD
      QV Q_SMD[NsubSMD];
      CalcQvSMD(mMuEvent, Q_SMD, mCent16, qZVtx);
      rpE = Q_SMD[0].get_Psi();
      rpW = Q_SMD[1].get_Psi();
    }
  }

  mFemtoEvent->SetRPeast(rpE);
  mFemtoEvent->SetRPwest(rpW);
  
  //
  // Loop over primary tracks
  //
  for(Int_t iTrk=0; iTrk<mNPrimTracks; iTrk++)
  {
    mPrimTrack = mMuDst->primaryTracks(iTrk);
    int idxGlob = mPrimTrack->index2Global();
    if (idxGlob < 0)   continue;
    mGlobTrack = mMuDst->globalTracks(idxGlob);

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

    if (!AcceptTrack(mPrimTrack, iVert))   continue;

    Bool_t mIsTofTrack = IsTofTrack(mGlobTrack);
    if(mIsTofTrack) {
      mBeta = mGlobTrack->btofPidTraits().beta();
      mMassSqr = mPrimTrack->momentum().mag2() * ( 1./(mBeta * mBeta) - 1.);
      //
      // Next lines are cheat. But we should remove garbage
      //
      if(mMassSqr < -0.02) continue;
      if(mPrimTrack->momentum().mag() > 0.5) {
        if(mMassSqr > 0.06 && mMassSqr < 0.18)   continue;
        if(mMassSqr > 0.35 && mMassSqr < 0.8)    continue;
        if(mMassSqr > 1.)                        continue;
      }
    } //if(mIsTofTrack)
    else {
      mBeta = -999.;
      mMassSqr = -999.;
    }

    Bool_t mIsTpcPion   = false;
    Bool_t mIsTpcKaon   = false;
    Bool_t mIsTpcProton = false;
    Bool_t mIsTofPion   = false;
    Bool_t mIsTofKaon   = false;
    Bool_t mIsTofProton = false;

    if(mRemovePions) {
      mIsTpcPion = IsTpcPion(mPrimTrack);
      mIsTofPion = IsTofPion(mGlobTrack);

      if(mIsTpcPion==true || mIsTofPion==true)
        continue;
    }

    if(mRemoveKaons) {
      mIsTpcKaon = IsTpcKaon(mPrimTrack);
      mIsTofKaon = IsTofKaon(mGlobTrack);

      if(mIsTpcKaon==true || mIsTofKaon==true)
        continue;
    }

    if(mRemoveProtons) {
      mIsTpcProton = IsTpcProton(mPrimTrack);
      mIsTofProton = IsTofProton(mGlobTrack);

      if(mIsTpcProton==true || mIsTofProton==true)
        continue;
    }

    mFTrack = mFemtoEvent->AddFemtoTrack();

    mFTrack->SetId( mPrimTrack->id() );
    mFTrack->SetNHits( (Char_t)(mPrimTrack->charge()*mPrimTrack->nHits()) );
    mFTrack->SetNHitsPoss( mPrimTrack->nHitsPoss() );
    //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->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++)

  mFemtoEvent->SetRefMult2(refMult2);
  mFemtoEvent->SetRefMult2Pos(refMult2Pos);

  if(mEventIsGood) {
    mNBytes += mTree->Fill();
    mFemtoEvent->Clear();
  }


  return kStOk;
}

//_________________
Bool_t StFemtoDstMaker::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 StFemtoDstMaker::AcceptPrimVtx(StThreeVectorF vtxPos, 
                                      Float_t vpdVz) {
  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;

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

  return mIsGoodVtx;
}

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

  Bool_t mIsGoodTrack = 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->nHitsFit() >= mTrackNHitsFit[0] && trk->nHitsFit() <= mTrackNHitsFit[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];

  return mIsGoodTrack;
}

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

  Bool_t mIsGoodTrack = 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;

  return mIsGoodTrack;
}

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

  Bool_t mIsTofHit = trk->btofPidTraits().beta() > 0 &&
    trk->btofPidTraits().timeOfFlight() > 0;

  return mIsTofHit;
}

//_________________
Int_t StFemtoDstMaker::Finish() {

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

  std::cout
    << "*************************************" << "\n"
    << "StFemtoDstMaker has been finished" << "\n"
    << "\t nEventsPassed   : " << mNEventsPassed  << "\n"
    << "\t nEventsProcessed: " << mNEventsIn << "\n"
    << "Accept trigger fail = " << mAcceptTriggerFail << "\n"
    << "RefMult fail = " << mRefMultFail << "\n"
    << "Anti-pileup fail = " << mAntiPileupFail << "\n"
    << "Ranking fail = " << mRankingFail << "\n"
    << "VtxNull fail = " << mVtxNullFail << "\n"
    << "NPrimVtx fail = " << mNPrimVtxFail << "\n"
    << "AcceptPrimVtx fail = " << mAcceptPrimVtxFail << "\n"
    << "*************************************" << "\n";

  return kStOk;
}

//_________________
void StFemtoDstMaker::CleanVariables() {

  mEventIsGood = false;
}

//_________________
Bool_t StFemtoDstMaker::IsTpcPion(StMuTrack *pTrk) {

  Bool_t isPion = false;

  if( pTrk->momentum().mag() <= mPthresh && 
      pTrk->nSigmaPion() >= mTpcPionNSigma[0] &&
      pTrk->nSigmaPion() <= mTpcPionNSigma[1] ) {
    isPion = true;
  }

  return isPion;
}

//_________________
Bool_t StFemtoDstMaker::IsTpcKaon(StMuTrack *pTrk) {

  Bool_t isKaon = false;

  if( pTrk->momentum().mag() <= mPthresh && 
      pTrk->nSigmaKaon() >= mTpcKaonNSigma[0] &&
      pTrk->nSigmaKaon() <= mTpcKaonNSigma[1] ) {
    isKaon = true;
  }

  return isKaon;
}

//_________________
Bool_t StFemtoDstMaker::IsTpcProton(StMuTrack *pTrk) {

  Bool_t isProton = false;

  if( pTrk->momentum().mag() <= mPthresh && 
      pTrk->nSigmaProton() >= mTpcProtonNSigma[0] &&
      pTrk->nSigmaProton() <= mTpcProtonNSigma[1] ) {
    isProton = true;
  }

  return isProton;
}

//_________________
Bool_t StFemtoDstMaker::IsTofPion(StMuTrack *gTrk) {

  Bool_t isPion = false;
  Float_t beta = -999.;
  Float_t mSqr = -999.;
  if( gTrk->btofPidTraits().beta() > 0 &&
      gTrk->momentum().mag() > mPthresh ) {
    beta = gTrk->btofPidTraits().beta();
    mSqr = gTrk->momentum().mag2() * (1./(beta*beta) - 1.);

    if(mSqr >= mTofPionMassSqr[0] &&
       mSqr <= mTofPionMassSqr[1]) {
      isPion = true;
    }
  }

  return isPion;
}

//_________________
Bool_t StFemtoDstMaker::IsTofKaon(StMuTrack *gTrk) {

  Bool_t isKaon = false;
  Float_t beta = -999.;
  Float_t mSqr = -999.;
  if( gTrk->btofPidTraits().beta() > 0 &&
      gTrk->momentum().mag() > mPthresh ) {
    beta = gTrk->btofPidTraits().beta();
    mSqr = gTrk->momentum().mag2() * (1./(beta*beta) - 1.);

    if(mSqr >= mTofKaonMassSqr[0] &&
       mSqr <= mTofKaonMassSqr[1]) {
      isKaon = true;
    }
  }

  return isKaon;
}

//_________________
Bool_t StFemtoDstMaker::IsTofProton(StMuTrack *gTrk) {

  Bool_t isProton = false;
  Float_t beta = -999.;
  Float_t mSqr = -999.;
  if( gTrk->btofPidTraits().beta() > 0 &&
      gTrk->momentum().mag() > 0.5 ) {

    //gTrk->momentum().mag() > mPthresh ) {
    beta = gTrk->btofPidTraits().beta();
    mSqr = gTrk->momentum().mag2() * (1./(beta*beta) - 1.);

    if(mSqr >= mTofProtonMassSqr[0] &&
       mSqr <= mTofProtonMassSqr[1]) {
      isProton = true;
    }
  }

  return isProton;
}

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

//_________________
int StFemtoDstMaker::LoadEPCalibParam(int runid) {
  if (!mCalibPath)   return 0;
  
  TFile *in = new TFile(Form("%s/OutMuAna_%d.root", mCalibPath, runid));
  if (!in)   return 0;
	
  mQVMaker->LoadParameters(in);

  in->Close();
  delete in;
  return 1;
}

//_________________
int StFemtoDstMaker::GetCentralityBin(float refMult, int fVz) {
  const int multMax = 1000;
  const int mMultiplicityCut[2][17] = {{multMax, /* 0      */   // vpd-zdce-tac-protected
                                        273,     /* 0   1  */   //-30<vz<-25
                                        235,     /* 2   1  */
                                        201,     /* 2   3  */
                                        270,     /* 4   3  */
                                        142,     /* 4   5  */
                                        118,     /* 6   5  */
                                        97,      /* 6   7  */
                                        78,      /* 8   7  */
                                        62,      /* 8   9  */
                                        49,      /* 10  9  */
                                        38,      /* 10  11 */
                                        29,      /* 12  11 */
                                        22,      /* 12  13 */
                                        16,      /* 14  13 */
                                        12,      /* 14  15 */
                                        8},      /*     15 */
                                       {multMax, // vz>-25
                                        273, 
                                        235, 
                                        201, 
                                        170, 
                                        142, 
                                        118, 
                                        97, 
                                        78, 
                                        62, 
                                        49, 
                                        38, 
                                        29, 
                                        22, 
                                        16, 
                                        12, 
                                        8}}; 

  // NqvCent = 16 5% step up to 80%
  int i_vz = fVz > 0 ? 1 : 0;
  
  for (int iCent = 1; iCent <= NqvCent; iCent++) {
    if (mMultiplicityCut[i_vz][iCent] <= refMult &&
        refMult < mMultiplicityCut[i_vz][iCent - 1]) {
      
      return iCent - 1;
    }
  }
  return -1;
}

//_________________
void StFemtoDstMaker::CalcQvSMD(StMuEvent *mEv, QV fQv[NsubSMD], int fCent, int fVz) {

  // Calculate raw flow vectors
  float fqv[NsubSMD][4] = {{0}}; //east-west-combined, xy+weight
  for (int iSub = 0; iSub < NsubSMD; iSub++) { // east-west

    int nStrip;
    for (int iXY = 0;iXY < 2; iXY++) {
      if (iXY == 0)   nStrip = 8;  // vertical strips (x-direction)
      else            nStrip = 9;  // horizontal strips (y-direction)

      for (int iStrip = 1; iStrip < nStrip; iStrip++) {
        // fqv[iSub][1-2] = position(cm) * adcThisTile
        fqv[iSub][iXY]   += ZDCSMD_GetPosition(iSub, iXY, iStrip)*ZDCSMD(mEv, iSub, iXY, iStrip);
        // fqv[iSub][3-4] = adcThisTile
        fqv[iSub][iXY+2] += ZDCSMD(mEv, iSub, iXY, iStrip);
      }

      // normalization
      if(fqv[iSub][iXY+2] > 0)
        fqv[iSub][iXY] /= fqv[iSub][iXY+2];
      else
        fqv[iSub][iXY] = -9999;
      
      if (iSub == 0 && mEv->zdcTriggerDetector().adcSum(east) < 1)
        fqv[iSub][iXY] = -9999;
      if (iSub == 1 && mEv->zdcTriggerDetector().adcSum(west) < 1)
        fqv[iSub][iXY] = -9999;
      if (iSub == 2 && (fqv[0][iXY] < -9000 || fqv[1][iXY] < -9000))
        fqv[iSub][iXY] = -9999;
    } // iXY
  } // iSub
  
  // combined
  for (int iXY = 0; iXY < 2; iXY++) {
    fqv[2][iXY] = fqv[0][iXY] - fqv[1][iXY];
    if (fqv[0][iXY] == -9999 || fqv[1][iXY] == -9999)
      fqv[2][iXY] = -9999;
  }

  
  // Set raw Q-vectors and Calibration 
  // -----------------------------------------------
  int fth   = 0; // fth = fOrd in QV class
  float fQw = 1.0;
  for (int iSub = 0; iSub < NsubSMD; iSub++) {
    if (fqv[iSub][0] == -9999 || fqv[iSub][1] == -9999 )   fQw = -1;
    fQv[iSub].set_Qv(fqv[iSub][0], fqv[iSub][1], fQw, dZDC, iSub, fth, fCent, fVz);

    mQVMaker->doCalibration(&(fQv[iSub]));
  } // iSub
}

//_________________
Float_t StFemtoDstMaker::ZDCSMD(StMuEvent *mEv, int eastwest, int verthori, int strip) {
  
  // return adc for each tile
  return mEv->zdcTriggerDetector().zdcSmd((StBeamDirection)eastwest, verthori, strip);
}

//_________________
Float_t StFemtoDstMaker::ZDCSMD_GetPosition(int eastwest, int verthori, int strip) {

  // Get position of each slat;strip starts from 1
  Float_t zdcsmd_x[7] = {0.5, 2, 3.5, 5, 6.5, 8, 9.5};
  Float_t zdcsmd_y[8] = {1.25, 3.25, 5.25, 7.25, 9.25, 11.25, 13.25, 15.25};

  if (eastwest == 0 && verthori == 0) { return zdcsmd_x[strip-1] - mZDCSMDCenterex; }
  // east && vertical strips (x-direction)
  if (eastwest == 1 && verthori == 0) { return mZDCSMDCenterwx - zdcsmd_x[strip-1]; }
  // west && vertical strips (x-direction)
  if (eastwest == 0 && verthori == 1) { return zdcsmd_y[strip-1]/sqrt(2.) - mZDCSMDCenterey; }
  // east && horisontal strips (y-direction)
  if (eastwest == 1 && verthori == 1) { return zdcsmd_y[strip-1]/sqrt(2.) - mZDCSMDCenterwy; }
  // west && horisontal strips (y-direction)
  
  return 0;
}