my_star/StRoot/StHbtMaker/Cut/gregsFemtoDstEventCut.cxx

#include "StHbtMaker/Cut/gregsFemtoDstEventCut.h"
#include <TMath.h>
#include <sstream>

#ifdef __ROOT__
ClassImp(gregsFemtoDstEventCut)
#endif

//_________________
gregsFemtoDstEventCut::gregsFemtoDstEventCut() {

  mRefMult[0] = 0.;
  mRefMult[1] = 5000;
  mVtxX[0] = -1e9;
  mVtxX[1] = +1e9;
  mVtxY[0] = -1e9;
  mVtxY[1] = +1e9;
  mVtxZ[0] = -1e9;
  mVtxZ[1] = +1e9;
  mVtxR = 1e9;
  mVzVpdVzDiff[0] = -50;
  mVzVpdVzDiff[1] = +50;
  mCollisionType = 0; //0-pp; 1-AuAu
  mNumberOfCentralityBins = 9;
  mCentralityBinRange[0] = 0;
  mCentralityBinRange[1] = 8;

  mVtxXShift = 0.;
  mVtxYShift = 0.;
  mNEventsPassed = mNEventsFailed = 0;
}

//_________________
gregsFemtoDstEventCut::gregsFemtoDstEventCut(const gregsFemtoDstEventCut& c) :StHbtEventCut(c) {
  
  mRefMult[0] = c.mRefMult[0];
  mRefMult[1] = c.mRefMult[1];
  mVtxX[0] = c.mVtxX[0];
  mVtxX[1] = c.mVtxX[1];
  mVtxY[0] = c.mVtxY[0];
  mVtxY[1] = c.mVtxY[1];
  mVtxZ[0] = c.mVtxZ[0];
  mVtxZ[1] = c.mVtxZ[1];
  mVtxR = c.mVtxR;
  mVzVpdVzDiff[0] = c.mVzVpdVzDiff[0];
  mVzVpdVzDiff[1] = c.mVzVpdVzDiff[1];
  mCollisionType = c.mCollisionType; //0-pp; 1-AuAu
  mNumberOfCentralityBins = c.mNumberOfCentralityBins;
  mCentralityBinRange[0] = c.mCentralityBinRange[0];
  mCentralityBinRange[1] = c.mCentralityBinRange[1];

  mVtxXShift = c.mVtxXShift;
  mVtxYShift = c.mVtxYShift;

  mNEventsPassed = 0;
  mNEventsFailed = 0;
}

//_________________
gregsFemtoDstEventCut::~gregsFemtoDstEventCut() {
  /* no-op */
}

//_________________
bool gregsFemtoDstEventCut::Pass(const StHbtEvent *event) {

  //added by gnigmat
#ifdef STHBTDEBUG
  std::cout << "gregsFemtoDstEventCut::Pass" << std::endl;
#endif

  int    mult = event->RefMult();
  double vtxX = event->PrimVertPos().x();
  double vtxY = event->PrimVertPos().y();
  double vtxZ = event->PrimVertPos().z();
  double vtxXCorrected = vtxX - mVtxXShift;
  double vtxYCorrected = vtxY - mVtxYShift;
  double vtxR = TMath::Sqrt(vtxXCorrected * vtxXCorrected + 
			    vtxYCorrected * vtxYCorrected);
  float  vzVpdVzDiff = vtxZ - event->VpdVz();
  int mCentBin;
  if(mCollisionType == 1) { //If Au+Au collisions
    if(mNumberOfCentralityBins == 9) {
      mCentBin = event->Cent9();
    }
    else if(mNumberOfCentralityBins == 16) {
      mCentBin = event->Cent16();
    }
    else {
      std::cout << "Wrong centrality definition: " 
                << mNumberOfCentralityBins << std::endl;
      mCentBin = -1;
    }
  } //if(mCollisionType == 1)
  

  bool mGoodEvent = false;
  bool mGoodPrimVtx = false;
  bool mGoodRefMult = false;

  if(mCollisionType == 0) { // p+p collisions

    mGoodPrimVtx = ( (vtxXCorrected >= mVtxX[0]) && 
      (vtxXCorrected <= mVtxX[1]) &&
      (vtxYCorrected >= mVtxY[0]) && 
      (vtxYCorrected <= mVtxY[1]) &&
      (vtxZ >= mVtxZ[0]) && (vtxZ <= mVtxZ[1]) &&
      (vtxR <= mVtxR) &&
      (vzVpdVzDiff >= mVzVpdVzDiff[0]) && 
      (vzVpdVzDiff <= mVzVpdVzDiff[1]) );
    mGoodRefMult = ( (mult >= mRefMult[0]) && (mult<= mRefMult[1]) );
  }
  else { //Au+Au collisions
    mGoodPrimVtx = ( (vtxZ >= mVtxZ[0]) && (vtxZ <= mVtxZ[1]) &&
		     (vtxR <= mVtxR) &&
		     (vzVpdVzDiff >= mVzVpdVzDiff[0]) && 
                     (vzVpdVzDiff <= mVzVpdVzDiff[1]) );

    mGoodRefMult = ( (mCentBin >= mCentralityBinRange[0]) &&
                     (mCentBin <= mCentralityBinRange[1]) );
  }

  //Main decision
  if(mGoodPrimVtx && mGoodRefMult) {
    mGoodEvent = true;
  }
  else { //This is not needed but it's easy to check
    mGoodEvent = false;
  }

  mGoodEvent ? mNEventsPassed++ : mNEventsFailed++;
  return (mGoodEvent);
}

//_________________
string gregsFemtoDstEventCut::Report() {

  stringstream rep_stream;
  if(mCollisionType==0) {
    rep_stream << mRefMult[0] << " <= Multiplicity <= " << mRefMult[1] 
	       << std::endl;
  }
  else {
    rep_stream << mCentralityBinRange[0] << " <= CentralityBin <= " 
	       << mCentralityBinRange[1] << std::endl;
  }
  rep_stream << mVtxX[0] << " <= VertexX <= " << mVtxX[1] << std::endl
	     << mVtxY[0] << " <= VertexY <= " << mVtxY[1] << std::endl
	     << mVtxZ[0] << " <= VertexZ <= " << mVtxZ[1] << std::endl
	     << "VtxXShift: " << mVtxXShift << "\tVtxYShift: " 
	     << mVtxYShift << std::endl
	     << "Vertex R <= " << mVtxR << std::endl
	     << mVzVpdVzDiff[0] << " <= VpdVzDifference <= " 
	     << mVzVpdVzDiff[1] << std::endl
	     << "Number of events passed: " << mNEventsPassed << std::endl
	     << "Number of events failed: " << mNEventsFailed << std::endl;

  return ( rep_stream.str() );
}