my_star/StRoot/StMuDstQAMaker/StMuDstMinvMaker.cxx

#include "StMuDstMinvMaker.h"
#include "StBTofHeader.h"
#include "TBranch.h"

ClassImp(StMuDstMinvMaker)

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



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

  std::cout << "StMuDstMinvMaker::StMuDstMinvMaker - Creating an instance...";
  mMuDstMaker = muMaker;

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

  mIsGoodTrack = false;
  mCurrentTrigger = 0;

  mFileName = oFileName;
  //
  // 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
  //
  mTrackP[0] = 0.1;
  mTrackP[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;

  /*
  //
  // Initialize TPC and TOF PID
  //
  mPionPionNSigma[0] = -3.;   mPionPionNSigma[1] = 3.;
  mPionKaonNSigma[0] = 1.;    mPionKaonNSigma[1] = -1.;
  */
  mTTTPThreshold = 0.7;

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

//_________________
StMuDstMinvMaker::~StMuDstMinvMaker() {
  /* nothing to do */
}

//_________________
Int_t StMuDstMinvMaker::Init() {

  std::cout << "StMuDstMinvMaker::Init - Initializing the maker"
      << std::endl;

  //
  // Create histograms and output file
  //
  mOutFile = new TFile(mFileName, "RECREATE");
  //
  // General event distributions
  //
  hR = new TH1F("hR", "R=#sqrt{v_{x}^{2} + v_{y}^{2}};R (cm);#",
                200, 0., 2.);
  hVx = new TH1F("hVx", ";v_{x} (cm);#", 200, -1.5, 1.5);
  hVy = new TH1F("hVy", ";v_{y} (cm);#", 200, -1.5, 1.5);
  hVz = new TH1F("hVz", ";v_{z} (cm);#", 100, -70., 70.);
  hVxVsVy = new TH2F("hVxVsVy", ";v_{x} (cm); v_{y} (cm)",
                     600, -1.5, 1.5,
                     600, -1.5, 1.5);
  hVxVsVz = new TH2F("hVxVsVz", ";v_{x} (cm); v_{z} (cm)",
                     600, -1.5, 1.5,
                     600, -70., 70.);
  hVyVsVz = new TH2F("hVyVsVz", ";v_{y} (cm); v_{z} (cm)",
                     600, -1.5, 1.5,
                     600, -70., 70.);
  hNPrimTr = new TH1F("hNPrimTr", ";N_{primary tracks};#", 150, 0, 1500);
  hNGlobTr = new TH1F("hNGlobTr", ";N_{global tracks};#", 150, 0, 1500);
  hTofRefMult = new TH1F("hTofRefMult", ";TofRefMult;#", 100, 0, 1000);
  hTofRefMultVsRefMult = new TH2F("hTofRefMultVsRefMult", ";TofRefMult;RefMult",
                                  500, 0., 1000,
                                  500, 0., 1000);
  hVpd = new TH1F("hVpd", ";v_{z}^{vpd} (cm);#", 100, -70., 70.);
  hVpdVz = new TH1F("hVpdVz", ";v_{z}^{vpd}-v_{z}", 100, -10., 10.);
  hNPrimVtx = new TH1F("hNPrimVtx", ";N_{primary vertex};#", 10, 0, 10);
  //
  // General track distributions
  //
  hP = new TH1F("hP", ";p (GeV/c);#", 200, 0., 5.);
  hPt = new TH1F("hPt", ";p_{t} (GeV/c);#", 200, 0., 5.);
  hPx = new TH1F("hPx", ";p_{x} (GeV/c);#", 200, 0., 5.);
  hPy = new TH1F("hPy", ";p_{y} (GeV/c);#", 200, 0., 5.);
  hPz = new TH1F("hPz", ";p_{z} (GeV/c);#", 200, 0., 5.);
  hEta = new TH1F("hEta", ";#eta;#", 200, -2., 2.);
  hPAccept = new TH1F("hPAccept", "accepted tracks;p (GeV/c);#", 200, 0., 3.);
  hPtAccept = new TH1F("hPtAccept", "accepted tracks;p_{t} (GeV/c);#", 200, 0., 3.);
  hPxAccept = new TH1F("hPxAccept", "accepted tracks;p_{x} (GeV/c);#", 200, 0., 3.);
  hPyAccept = new TH1F("hPyAccept", "accepted tracks;p_{y} (GeV/c);#", 200, 0., 3.);
  hPzAccept = new TH1F("hPzAccept", "accepted tracks;p_{z} (GeV/c);#", 200, 0., 3.);
  hEtaAccept = new TH1F("hEtaAccept", "accepted tracks;#eta;#", 200, -2., 2.);
  hMassSqrVsPt = new TH2F("hMassSqrVsPt", ";p_{t}Q (GeV/c);m^{2} (GeV/c^{2})",
                          400, -2., 2.,
                          1000, -0.1, 1.5);
  hDedxVsPt = new TH2F("hDedxVsPt", ";p_{t}Q (GeV/c);dE/dx (keV/cm)",
                       400, -2., 2.,
                       400, 0., 20.);
  hInvBetaExpVsPt = new TH2F("hInvBetaExpVsPt", ";p_{t}Q (GeV/c);1/beta_{exp}",
                             400, -2., 2.,
                             200, 0., 2.);
  hInvBetaThVsPt = new TH2F("hInvBetaThVsPt", ";p_{t}Q (GeV/c);1/beta_{exp}",
                            400, -2., 2.,
                            200, 0., 2.);
  hTOF = new TH1F("hTOF", "time of flight;t ns;#", 100, 0., 50.);
  // Pion PID
  hNSigmaPionPion = new TH1F("hNSigmaPionPion", "#pi;n#sigma(#pi^{#pm});#",
                             100, -5., 5.);
  hNSigmaPionKaon = new TH1F("hNSigmaPionKaon", "#pi;n#sigma(#K^{#pm});#",
                             100, -5., 5.);
  hNSigmaPionProton = new TH1F("hNSigmaPionProton", "#pi;n#sigma(p);#",
                               100, -5., 5.);
  hNSigmaPionPionVsPt = new TH2F("hNSigmaPionPionVsPt",
                                 "#pi;p_{t}Q (GeV/c);n#sigma(#pi^{#pm})",
                                 400, -2., 2.,
                                 1000, -5., 5.);
  hNSigmaPionKaonVsPt = new TH2F("hNSigmaPionKaonVsPt",
                                 "#pi;p_{t}Q (GeV/c);n#sigma(K^{#pm})",
                                 400, -2., 2.,
                                 1000, -5., 5.);
  hNSigmaPionProtonVsPt = new TH2F("hNSigmaPionProtonVsPt",
                                   "#pi;p_{t}Q (GeV/c);n#sigma(p)",
                                   400, -2., 2.,
                                   1000, -5., 5.);
  hDMassSqrVsPtPionTPC = new TH2F("hDMassSqrVsPtPionTPC",
                                  "TPC PID;p_{t} (GeV/c);dM(#pi^{#pm}) (GeV/c^{2})",
                                  400, -2., 2.,
                                  500, -.5, .5);
  hDMassSqrVsPtPionTOF = new TH2F("hDMassSqrVsPtPionTOF",
                                  "TOF PID;p_{t} (GeV/c);dM(#pi^{#pm}) (GeV/c^{2})",
                                  400, -2., 2.,
                                  500, -.5, .5);
  hInvBetaExpThPionTNT = new TH2F("hInvBetaExpThPionTNT", "TPC and TOF #pi^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
                                  400, -2., 2.,
                                  500, -0.5, 0.5);
  hInvBetaExpThPionTPC = new TH2F("hInvBetaExpThPionTPC", "only TPC #pi^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
                                  400, -2., 2.,
                                  500, -0.5, 0.5);
  hInvBetaExpThPionTOF = new TH2F("hInvBetaExpThPionTOF", "only TOF #pi^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
                                  400, -2., 2.,
                                  500, -0.5, 0.5);
  hInvBetaExpThPionTTT = new TH2F("hInvBetaExpThPionTTT",
                                  Form("if p #leq %.2f then TPC, else TOF #pi^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
                                       mTTTPThreshold),
                                  400, -2., 2.,
                                  500, -0.5, 0.5);
  hInvBetaExpThPionINT = new TH2F("hInvBetaExpThPionINT",
                                  "#pi^{#pm};p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
                                  400, -2., 2.,
                                  500, -5., 5.);
  // Kaon PID
  hNSigmaKaonPion = new TH1F("hNSigmaKaonPion", "K;n#sigma(#pi^{#pm}};#",
                             100, -5., 5.);
  hNSigmaKaonKaon = new TH1F("hNSigmaKaonKaon", "K;n#sigma(#K^{#pm}};#",
                             100, -5., 5.);
  hNSigmaKaonProton = new TH1F("hNSigmaKaonProton", "K;n#sigma(p);#",
                               100, -5., 5.);
  hNSigmaKaonPionVsPt = new TH2F("hNSigmaKaonPionVsPt",
                                 "K;p_{t}Q (GeV/c);n#sigma(#pi^{#pm})",
                                 400, -2., 2.,
                                 1000, -5., 5.);
  hNSigmaKaonKaonVsPt = new TH2F("hNSigmaKaonKaonVsPt",
                                 "K;p_{t}Q (GeV/c);n#sigma(K^{#pm})",
                                 400, -2., 2.,
                                 1000, -5., 5.);
  hNSigmaKaonProtonVsPt = new TH2F("hNSigmaKaonProtonVsPt",
                                   "K;p_{t}Q (GeV/c);n#sigma(p)",
                                   400, -2., 2.,
                                   1000, -5., 5.);
  hDMassSqrVsPtKaonTPC = new TH2F("hDMassSqrVsPtKaonTPC",
                                  "TPC PID;p_{t} (GeV/c);dM(K^{#pm}) (GeV/c^{2})",
                                  400, -2., 2.,
                                  500, -.5, .5);
  hDMassSqrVsPtKaonTOF = new TH2F("hDMassSqrVsPtKaonTOF",
                                  "TOF PID;p_{t} (GeV/c);dM(K^{#pm}) (GeV/c^{2})",
                                  400, -2., 2.,
                                  500, -.5, .5);
  hInvBetaExpThKaonTNT = new TH2F("hInvBetaExpThKaonTNT", "TPC and TOF K^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
                                  400, -2., 2.,
                                  500, -0.5, 0.5);
  hInvBetaExpThKaonTPC = new TH2F("hInvBetaExpThKaonTPC", "only TPC K^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
                                  400, -2., 2.,
                                  500, -0.5, 0.5);
  hInvBetaExpThKaonTOF = new TH2F("hInvBetaExpThKaonTOF", "only TOF K^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
                                  400, -2., 2.,
                                  500, -0.5, 0.5);
  hInvBetaExpThKaonTTT = new TH2F("hInvBetaExpThKaonTTT",
                                  Form("if p #leq %.2f then TPC, else TOF K^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
                                       mTTTPThreshold),
                                  400, -2., 2.,
                                  500, -0.5, 0.5);
  hInvBetaExpThKaonINT = new TH2F("hInvBetaExpThKaonINT",
                                  "K^{#pm};p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
                                  400, -2., 2.,
                                  500, -5., 5.);
  // Proton PID
  hNSigmaProtonPion = new TH1F("hNSigmaProtonPion", "p;n#sigma(#pi^{#pm}};#",
                               100, -5., 5.);
  hNSigmaProtonKaon = new TH1F("hNSigmaProtonKaon", "p;n#sigma(#K^{#pm}};#",
                               100, -5., 5.);
  hNSigmaProtonProton = new TH1F("hNSigmaProtonProton", "p;n#sigma(p);#",
                                 100, -5., 5.);
  hNSigmaProtonPionVsPt = new TH2F("hNSigmaProtonPionVsPt",
                                   "p;p_{t}Q (GeV/c);n#sigma(#pi^{#pm})",
                                   400, -2., 2.,
                                   1000, -5., 5.);
  hNSigmaProtonKaonVsPt = new TH2F("hNSigmaProtonKaonVsPt",
                                   "p;p_{t}Q (GeV/c);n#sigma(K^{#pm})",
                                   400, -2., 2.,
                                   1000, -5., 5.);
  hNSigmaProtonProtonVsPt = new TH2F("hNSigmaProtonProtonVsPt",
                                     "p;p_{t}Q (GeV/c);n#sigma(p)",
                                     400, -2., 2.,
                                     1000, -5., 5.);
  hDMassSqrVsPtProtonTPC = new TH2F("hDMassSqrVsPtProtonTPC",
                                    "TPC PID;p_{t} (GeV/c);dM(p) (GeV/c^{2})",
                                    400, -2., 2.,
                                    500, -.5, .5);
  hDMassSqrVsPtProtonTOF = new TH2F("hDMassSqrVsPtProtonTOF",
                                    "TOF PID;p_{t} (GeV/c);dM(p) (GeV/c^{2})",
                                    400, -2., 2.,
                                    500, -.5, .5);
  hInvBetaExpThProtonTNT = new TH2F("hInvBetaExpThProtonTNT", "TPC and TOF p PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
                                    400, -2., 2.,
                                    500, -0.5, 0.5);
  hInvBetaExpThProtonTPC = new TH2F("hInvBetaExpThProtonTPC", "only TPC p PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
                                    400, -2., 2.,
                                    500, -0.5, 0.5);
  hInvBetaExpThProtonTOF = new TH2F("hInvBetaExpThProtonTOF", "only TOF p PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
                                    400, -2., 2.,
                                    500, -0.5, 0.5);
  hInvBetaExpThProtonTTT = new TH2F("hInvBetaExpThProtonTTT",
                                    Form("if p #leq %.2f then TPC, else TOF p PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
                                         mTTTPThreshold),
                                    400, -2., 2.,
                                    500, -0.5, 0.5);
  hInvBetaExpThProtonINT = new TH2F("hInvBetaExpThProtonINT",
                                    "p;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}",
                                    400, -2., 2.,
                                    500, -5., 5.);
  //
  // Any charge
  //
  hRefMult = new TH1F("hRefMult", ";RefMult;#", 100, 0., 1000.);
  hRefMultAccept = new TH1F("hRefMultAccept", ";RefMultAccept;#", 100, 0., 1000.);

  std::cout << "StMuDstMinvMaker::Init - Initialization has been finished"
      << std::endl;
  return StMaker::Init();
}

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

//________________
Int_t StMuDstMinvMaker::Make() {

  mNEventsIn++;

  mMuDst = NULL;
  mMuEvent = NULL;

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

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

  //Obtaining MuEvent
  mMuEvent = (StMuEvent*)mMuDst->event();
  if(!AcceptTrigger(mMuEvent)) { //If trigger is not found
    return kStOk;
  }

  //Multiplicity cannot be negative
  unsigned short refMult = mMuEvent->refMult();
  unsigned short refMultPos = 0;
  unsigned short refMultNeg = 0;

  if(refMult < 0) {
    return kStOk;
  }

  hRefMult->Fill(refMult);

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

  //Some initializations of local variables
  mPrimVertex = NULL;
  StThreeVectorF mVertPosition;
  Float_t mVpdVz = 0.;
  Int_t   mPrimVertIndex = -999;
  Float_t mRanking = -999.;

  //Clean index vectors
  CleanVariables();

  //Vertex loop
  unsigned short refMultAccept = 0;
  unsigned short tofRefMult = 0;
  unsigned short refMultAcceptPos = 0;
  unsigned short refMultAcceptNeg = 0;
  for(Int_t iVert=0; iVert<mNVertices; iVert++) {

    mEventIsGood = false;

    hNPrimVtx->Fill(mNVertices);

    if(iVert!=0) continue; // Not first primary vertex does not contain tracks with fast detectors (TOF)

    mPrimVertex = mMuDst->primaryVertex(iVert);
    //Positive ranking only
    if(mPrimVertex->ranking() <= 0) continue;
    //Not (0,0,0) position of the primary vertex
    if(mPrimVertex->position().x() == 0 &&
       mPrimVertex->position().y() == 0 &&
       mPrimVertex->position().z() == 0) continue;
    //Reasonable amount of tracks
    if(mNPrimTracks < 0 || mNPrimTracks > 10000) continue;

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

    if( !AcceptPrimVtx(mVertPosition, mVpdVz) ) continue;

    mEventIsGood = true;

    hNPrimTr->Fill(mNPrimTracks);
    hNGlobTr->Fill(mNGlobTracks);
    //
    //Loop over primary tracks
    //
    for(int iTrk = 0; iTrk < mNPrimTracks; iTrk++) {
      mPrimTrack = mMuDst->primaryTracks(iTrk);
      mGlobTrack = mMuDst->globalTracks(mPrimTrack->index2Global());
      short charge = mPrimTrack->charge() > 0 ? 1 : -1;
      float eta = mPrimTrack->eta();
      float pt  = mPrimTrack->pt();
      float p   = mPrimTrack->p().mag();
      float p2  = mPrimTrack->p().mag2();
      float px  = mPrimTrack->p().x();
      float py  = mPrimTrack->p().y();
      float pz  = mPrimTrack->p().z();

      charge > 0 ? refMultPos++ : refMultNeg++;

      hP->Fill(p);
      hPt->Fill(pt);
      hPx->Fill(px);
      hPy->Fill(py);
      hPz->Fill(pz);
      hEta->Fill(eta);

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

      hPAccept->Fill(p);
      hPtAccept->Fill(pt);
      hPxAccept->Fill(px);
      hPyAccept->Fill(py);
      hPzAccept->Fill(pz);
      hEtaAccept->Fill(eta);

      refMultAccept++;
      charge > 0 ? refMultAcceptPos++ : refMultAcceptNeg++;

      float nSigmaPion = mPrimTrack->nSigmaPion();
      float nSigmaKaon = mPrimTrack->nSigmaKaon();
      float nSigmaProton = mPrimTrack->nSigmaProton();
      float dedx = mPrimTrack->dEdx();
      bool  tofTrack = IsTofTrack(mGlobTrack);
      float betaExp;
      float massSqr;
      bool pionTPC   = ((nSigmaPion   >= mPionPionNSigma[0]     || nSigmaPion   <= mPionPionNSigma[1])      && 
                        (nSigmaKaon   <  mPionKaonNSigma[0]     || nSigmaKaon   >  mPionKaonNSigma[1])      &&
                        (nSigmaProton <  mPionProtonNSigma[0]   || nSigmaProton >  mPionProtonNSigma[1]));
      bool kaonTPC   = ((nSigmaPion   <  mKaonPionNSigma[0]     || nSigmaPion   >  mKaonPionNSigma[1])      &&
                        (nSigmaKaon   >= mKaonKaonNSigma[0]     || nSigmaKaon   <= mKaonKaonNSigma[1])      &&
                        (nSigmaProton <  mKaonProtonNSigma[0]   || nSigmaProton >  mKaonProtonNSigma[1]));
      bool protonTPC = ((nSigmaPion   <  mProtonPionNSigma[0]   || nSigmaPion   >  mProtonPionNSigma[1])    &&
                        (nSigmaKaon   <  mProtonKaonNSigma[0]   || nSigmaKaon   >  mProtonKaonNSigma[1])    &&
                        (nSigmaProton >= mProtonProtonNSigma[0] || nSigmaProton <= mProtonProtonNSigma[1]));
      bool pionTOF   = false;
      bool kaonTOF   = false;
      bool protonTOF = false;
      //
      // Pion TPC
      //
      if (pionTPC) {
        hNSigmaPionPion->Fill(nSigmaPion);
        hNSigmaPionKaon->Fill(nSigmaKaon);
        hNSigmaPionProton->Fill(nSigmaProton);
        hNSigmaPionPionVsPt->Fill(pt*charge, nSigmaPion);
        hNSigmaPionKaonVsPt->Fill(pt*charge, nSigmaKaon);
        hNSigmaPionProtonVsPt->Fill(pt*charge, nSigmaProton);
      }
      //
      // Kaon TPC
      //
      if (kaonTPC) {
        hNSigmaKaonPion->Fill(nSigmaPion);
        hNSigmaKaonKaon->Fill(nSigmaKaon);
        hNSigmaKaonProton->Fill(nSigmaProton);
        hNSigmaKaonPionVsPt->Fill(pt*charge, nSigmaPion);
        hNSigmaKaonKaonVsPt->Fill(pt*charge, nSigmaKaon);
        hNSigmaKaonProtonVsPt->Fill(pt*charge, nSigmaProton);
      }
      //
      // Proton TPC
      //
      if (protonTPC) {
        hNSigmaProtonPion->Fill(nSigmaPion);
        hNSigmaProtonKaon->Fill(nSigmaKaon);
        hNSigmaProtonProton->Fill(nSigmaProton);
        hNSigmaProtonPionVsPt->Fill(pt*charge, nSigmaPion);
        hNSigmaProtonKaonVsPt->Fill(pt*charge, nSigmaKaon);
        hNSigmaProtonProtonVsPt->Fill(pt*charge, nSigmaProton);
      }

      hDedxVsPt->Fill(pt*charge, dedx*1e6);

      if (tofTrack) {
        hTOF->Fill(mGlobTrack->btofPidTraits().timeOfFlight());

        float betaThPion = sqrt(p2/(PION_MASS*PION_MASS + p2));
        float betaThKaon = sqrt(p2/(KAON_MASS*KAON_MASS + p2));
        float betaThProton = sqrt(p2/(PROTON_MASS*PROTON_MASS + p2));
        tofRefMult++;
        betaExp = mGlobTrack->btofPidTraits().beta();
        massSqr = p2*(1./(betaExp*betaExp) - 1.);
        hMassSqrVsPt->Fill(pt*charge, massSqr);
        pionTOF   = massSqr >= mPionMass[0]   && massSqr <= mPionMass[1];
        kaonTOF   = massSqr >= mKaonMass[0]   && massSqr <= mKaonMass[1];
        protonTOF = massSqr >= mProtonMass[0] && massSqr <= mProtonMass[1];

        float dMpion   = massSqr - PION_MASS*PION_MASS;
        float dMkaon   = massSqr - KAON_MASS*KAON_MASS;
        float dMproton = massSqr - PROTON_MASS*PROTON_MASS;
        hInvBetaExpVsPt->Fill(pt*charge, 1./betaExp);
        //
        // TPC PID
        //
        if (pionTPC) {
          hInvBetaExpThPionTPC->Fill(pt*charge, 1./betaExp - 1./betaThPion);
          hDMassSqrVsPtPionTPC->Fill(pt*charge, dMpion);
        }
        else if (kaonTPC) {
          hInvBetaExpThKaonTPC->Fill(pt*charge, 1./betaExp - 1./betaThKaon);
          hDMassSqrVsPtKaonTPC->Fill(pt, dMkaon);
        }
        else if (protonTPC) {
          hInvBetaExpThProtonTPC->Fill(pt*charge, 1./betaExp - 1./betaThProton);
          hDMassSqrVsPtProtonTPC->Fill(pt, dMproton);
        }
        //
        // TOF PID
        //
        if (pionTOF) {
          hInvBetaExpThPionTOF->Fill(pt*charge, 1./betaExp - 1./betaThPion);
          hDMassSqrVsPtPionTOF->Fill(pt, dMpion);
          hInvBetaThVsPt->Fill(pt*charge, 1./betaThPion);
          hInvBetaExpThPionINT->Fill(pt*charge, 1./betaExp - 1./betaThPion);
          hInvBetaExpThPionINT->Fill(pt*charge, 1./betaExp - 1./betaThKaon);
          hInvBetaExpThPionINT->Fill(pt*charge, 1./betaExp - 1./betaThProton);
        }
        else if (kaonTOF) {
          hInvBetaExpThKaonTOF->Fill(pt*charge, 1./betaExp - 1./betaThKaon);
          hDMassSqrVsPtKaonTOF->Fill(pt, dMkaon);
          hInvBetaThVsPt->Fill(pt*charge, 1./betaThKaon);
          hInvBetaExpThKaonINT->Fill(pt*charge, 1./betaExp - 1./betaThPion);
          hInvBetaExpThKaonINT->Fill(pt*charge, 1./betaExp - 1./betaThKaon);
          hInvBetaExpThKaonINT->Fill(pt*charge, 1./betaExp - 1./betaThProton);
        }
        else if (protonTOF) {
          hInvBetaExpThProtonTOF->Fill(pt*charge, 1./betaExp - 1./betaThProton);
          hDMassSqrVsPtProtonTOF->Fill(pt, dMproton);
          hInvBetaThVsPt->Fill(pt*charge, 1./betaThProton);
          hInvBetaExpThProtonINT->Fill(pt*charge, 1./betaExp - 1./betaThPion);
          hInvBetaExpThProtonINT->Fill(pt*charge, 1./betaExp - 1./betaThKaon);
          hInvBetaExpThProtonINT->Fill(pt*charge, 1./betaExp - 1./betaThProton);	  
        }
        //
        // TNT PID
        //
        if (pionTOF && pionTPC) {
          hInvBetaExpThPionTNT->Fill(pt*charge, 1./betaExp - 1./betaThPion);	  
        }
        else if (kaonTOF && kaonTPC) {
          hInvBetaExpThKaonTNT->Fill(pt*charge, 1./betaExp - 1./betaThKaon);	  
        }
        else if (protonTOF && protonTPC) {
          hInvBetaExpThProtonTNT->Fill(pt*charge, 1./betaExp - 1./betaThProton);	  
        }
        //
        // TTT PID
        //
        if (p >= mTTTPThreshold) {
          if (pionTOF) {
            hInvBetaExpThPionTTT->Fill(pt*charge, 1./betaExp - 1./betaThPion);
          }
          else if (kaonTOF) {
            hInvBetaExpThKaonTTT->Fill(pt*charge, 1./betaExp - 1./betaThKaon);
          }
          else if (protonTOF) {
            hInvBetaExpThProtonTTT->Fill(pt*charge, 1./betaExp - 1./betaThProton);
          }
        } else {
          if (pionTPC) {
            hInvBetaExpThPionTTT->Fill(pt*charge, 1./betaExp - 1./betaThPion);
          }
          else if (kaonTPC) {
            hInvBetaExpThKaonTTT->Fill(pt*charge, 1./betaExp - 1./betaThKaon);
          }
          else if (protonTPC) {
            hInvBetaExpThProtonTTT->Fill(pt*charge, 1./betaExp - 1./betaThProton);
          }
        }
      }
    } //for(Int_t iTrk=0; iTrk<mNPrimTracks; iTrk++)
  } //for(Int_t iVert=0; iVert<mNVertices; iVert++)

  if (mEventIsGood) {
    hRefMultAccept->Fill(refMultAccept);
    hTofRefMult->Fill(tofRefMult);
    hTofRefMultVsRefMult->Fill(tofRefMult, refMultAccept);
  }

  return kStOk;
}

//_________________
Bool_t StMuDstMinvMaker::AcceptTrigger(StMuEvent *muEvent) {

  Bool_t mIsGoodTrigger = false;
  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 = mTriggerIdCollection.at(iTrg);
        break;
      }
    } //for(UInt_t iTrg=0; iTrg<mTriggerIdCollection.size(); iTrg++)
  }

  return mIsGoodTrigger;
}

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

  hR->Fill(vtxR);
  hVx->Fill(mVtxX);
  hVy->Fill(mVtxY);
  hVz->Fill(mVtxZ);
  hVxVsVy->Fill(mVtxX, mVtxY);
  hVxVsVz->Fill(mVtxX, mVtxZ);
  hVyVsVz->Fill(mVtxY, mVtxZ);
  hVpd->Fill(vpdVz);
  hVpdVz->Fill(vpdDiff);

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

//_________________
Bool_t StMuDstMinvMaker::AcceptTrack(StMuTrack *trk, UShort_t vtxInd) {
  float          mom      = trk->momentum().mag();
  float          eta      = trk->eta();
  unsigned short nHitsFit = trk->nHitsFit();
  float          ratioH   = (float)trk->nHitsFit()/(float)trk->nHitsPoss();
  short          flag     = trk->flag();
  float          dca      = trk->dcaGlobal(vtxInd).perp();

  bool isMomOk      = ( mom >= mTrackP[0] &&
                       mom <= mTrackP[1] );
  bool isEtaOk      = ( eta >= mTrackEta[0] &&
                       eta <= mTrackEta[1] );
  bool isNHitsFitOk = ( nHitsFit >= mTrackNHitsFit[0] &&
                       nHitsFit <= mTrackNHitsFit[1] &&
                       ratioH >= 0.51 );
  bool isFlagOk     = ( flag >= mTrackFlag[0] &&
                       flag <= mTrackFlag[1] );
  //  bool isDcaOk      = ( dca >= mTrackDcaGlobal[0] &&
  //                        dca <= mTrackDcaGlobal[1] );

  return ( isMomOk && isEtaOk && isNHitsFitOk && isFlagOk );
}


//_________________
Int_t StMuDstMinvMaker::Finish() {

  std::cout << "*************************************" << std::endl
      << "StMuDstMinvMaker has been finished" << std::endl
      << "\t nEventsPassed   : " << mNEventsPassed  << std::endl
      << "\t nEventsProcessed: " << mNEventsIn << std::endl
      << "*************************************" << std::endl;
  mOutFile->cd();
  mOutFile->Write();
  mOutFile->Close();

  return kStOk;
}

//_________________
void StMuDstMinvMaker::CleanVariables() {

  mEventIsGood = false;
}

//_________________
void StMuDstMinvMaker::SetTriggerId(unsigned int id) {
  mTriggerIdCollection.push_back(id);
}


//_________________
void StMuDstMinvMaker::SetMuDstMaker(StMuDstMaker *maker) {
  mMuDstMaker = maker;
}


//_________________
void StMuDstMinvMaker::SetVtxZCut(float lo, float hi) {
  mPrimVtxZ[0] = lo;
  mPrimVtxZ[1] = hi;
}


//_________________
void StMuDstMinvMaker::SetVtxRCut(float lo, float hi) {
  mPrimVtxR[0] = lo;
  mPrimVtxR[1] = hi;  
}


//_________________
void StMuDstMinvMaker::SetVtxShift(float xShift, float yShift) {
  mPrimVtxXShift = xShift;
  mPrimVtxYShift = yShift;
}


//_________________
void StMuDstMinvMaker::SetVtxVpdVzDiffCut(float lo, float hi) {
  mPrimVtxVpdVzDiff[0] = lo;
  mPrimVtxVpdVzDiff[1] = hi;
}


//_________________
void StMuDstMinvMaker::SetP(float lo, float hi) {
  mTrackP[0] = lo;
  mTrackP[1] = hi;
}


//_________________
void StMuDstMinvMaker::SetPt(float lo, float hi) {
  mTrackPt[0] = lo;
  mTrackPt[1] = hi;
}


//_________________
void StMuDstMinvMaker::SetTrackNHits(int lo, int hi) {
  mTrackNHits[0] = lo;
  mTrackNHits[1] = hi;  
}


//_________________
void StMuDstMinvMaker::SetTrackNHitsFit(int lo, int hi) {
  mTrackNHitsFit[0] = lo;
  mTrackNHitsFit[1] = hi;
}


//_________________
void StMuDstMinvMaker::SetTrackEta(float lo, float hi) {
  mTrackEta[0] = lo;
  mTrackEta[1] = hi;
}


//_________________
void StMuDstMinvMaker::SetTrackFlag(short lo, short hi) {
  mTrackFlag[0] = lo;
  mTrackFlag[1] = hi;
}


//_________________
void StMuDstMinvMaker::SetTrackDCA(float lo, float hi) {
  mTrackDcaGlobal[0] = lo;
  mTrackDcaGlobal[1] = hi;
}


//_________________
Bool_t StMuDstMinvMaker::IsTofTrack(StMuTrack *trk) { //Only for globals
  Bool_t isTofHit = false;
  if(trk->btofPidTraits().beta() > 0 &&
     trk->btofPidTraits().timeOfFlight() > 0) {
    isTofHit = true;
  }
  return isTofHit;
}


//_________________
void StMuDstMinvMaker::SetPionPionNSigma(float lo, float hi) {
  mPionPionNSigma[0] = lo;   mPionPionNSigma[1] = hi;
}


//_________________
void StMuDstMinvMaker::SetPionKaonNSigma(float lo, float hi) {
  mPionKaonNSigma[0] = lo;   mPionKaonNSigma[1] = hi;
}


//_________________
void StMuDstMinvMaker::SetPionProtonNSigma(float lo, float hi) {
  mPionProtonNSigma[0] = lo;   mPionProtonNSigma[1] = hi;
}


//_________________
void StMuDstMinvMaker::SetKaonPionNSigma(float lo, float hi) {
  mKaonPionNSigma[0] = lo;   mKaonPionNSigma[1] = hi;
}


//_________________
void StMuDstMinvMaker::SetKaonKaonNSigma(float lo, float hi) {
  mKaonKaonNSigma[0] = lo;   mKaonKaonNSigma[1] = hi;
}


//_________________
void StMuDstMinvMaker::SetKaonProtonNSigma(float lo, float hi) {
  mKaonProtonNSigma[0] = lo;   mKaonProtonNSigma[1] = hi;
}


//_________________
void StMuDstMinvMaker::SetProtonPionNSigma(float lo, float hi) {
  mProtonPionNSigma[0] = lo;   mProtonPionNSigma[1] = hi;
}


//_________________
void StMuDstMinvMaker::SetProtonKaonNSigma(float lo, float hi) {
  mProtonKaonNSigma[0] = lo;   mProtonKaonNSigma[1] = hi;
}


//_________________
void StMuDstMinvMaker::SetProtonProtonNSigma(float lo, float hi) {
  mProtonProtonNSigma[0] = lo;   mProtonProtonNSigma[1] = hi;
}


//_________________
void StMuDstMinvMaker::SetPionMassSqr(float lo, float hi) {
  mPionMass[0] = lo;   mPionMass[1] = hi;
}


//_________________
void StMuDstMinvMaker::SetKaonMassSqr(float lo, float hi) {
  mKaonMass[0] = lo;   mKaonMass[1] = hi;
}


//_________________
void StMuDstMinvMaker::SetProtonMassSqr(float lo, float hi) {
  mProtonMass[0] = lo;   mProtonMass[1] = hi;
}


//_________________
void StMuDstMinvMaker::SetTTTPThreshold(float pThresh) {
  mTTTPThreshold = pThresh;
}