my_star/StRoot/StHbtMaker/Reader/gregsStHbtMuDstReader.cxx

#include "StChain.h"

#include "StEvent/StEventTypes.h"

#include "StEventUtilities/StuRefMult.hh"
#include "StEventUtilities/StuProbabilityPidAlgorithm.h"

#include "StarClassLibrary/StPhysicalHelixD.hh"

#include "StMuDSTMaker/COMMON/StMuException.hh"
#include "StMuDSTMaker/COMMON/StMuTrack.h"
#include "StMuDSTMaker/COMMON/StMuDebug.h"
#include "StMuDSTMaker/COMMON/StMuCut.h"
#include "StMuDSTMaker/COMMON/StMuDst.h"
#include "StBTofHeader.h"

#include "gregsStHbtMuDstReader.h"
#include "StHbtMaker/Infrastructure/StHbtEvent.hh"
#include "StHbtMaker/Base/StHbtEventCut.h"

#include "TFile.h"
#include "TTree.h"
#include "TClass.h"
#include "TChain.h"
#include "TStreamerInfo.h"
#include "TClonesArray.h"

ClassImp(gregsStHbtMuDstReader)

#if !(ST_NO_NAMESPACES)
  using namespace units;
#endif


//_________________
gregsStHbtMuDstReader::gregsStHbtMuDstReader(StMuDstMaker* maker) : 
  mMuDstMaker(maker), mTrackType(1), mReadTracks(1), 
  mReadV0s(0) {

  mEventCounter=0;
  mReaderStatus = 0;
  mProcessAllVertices = false;
  mFinish = false;
  mHbtEvent = 0;
}

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

//_________________
void gregsStHbtMuDstReader::clear(){
  DEBUGMESSAGE1("");
  //  if (mHbtEvent) { delete mHbtEvent; mHbtEvent=0; }
}

//_________________
int gregsStHbtMuDstReader::Init(){
  DEBUGMESSAGE1("");
//   if (!ioMaker) mIOMaker = (StIOMaker*)GetMaker("IOMaker");
//   if (!mStStrangeMuDstMaker) mStStrangeMuDstMaker = (StStrangeMuDstMaker*)GetMaker("StrangeMaker");
//   if (!mFlowMaker)  = (StFlowMaker*)GetMaker("FlowMaker");
  return 0;
}

//_________________
void gregsStHbtMuDstReader::Clear(){
  DEBUGMESSAGE1("");
  clear();
}

//_________________
StHbtEvent* gregsStHbtMuDstReader::ReturnHbtEvent(){

  DEBUGMESSAGE1("");

  clear();
  mMuDst=mMuDstMaker->muDst();

  if(!mMuDst) {
    std::cout << "gregsStHbtMuDstReader::ReturnHbtEvent [WARNING] - No MuDst file found"
	      << std::endl;
    return 0;
  }

  mHbtEvent = 0;
  mPrimVert = 0;

  //Event processing
  if (mMuDst && mMuDst->event()) {

    DEBUGVALUE3(mMuDst);
    mMuDst->fixTrackIndices();

    StMuEvent* mMuEvent = mMuDst->event();
    StBTofHeader *mTofHeader = mMuDst->btofHeader();
    float mVpdVz;

    //If trigger collection is not empty then check the trigger list
    bool mTriggerIsPassed = IsGoodTrigger(mMuEvent);
    if(mTriggerIsPassed == false) {
      return 0;
    }

    //Get number of primary vertices and tracks
    Int_t mPrimVertNum = mMuDst->numberOfPrimaryVertices();
    Int_t mPrimTracksNum = mMuDst->numberOfPrimaryTracks();
    Int_t mGlobTracksNum = mMuDst->numberOfGlobalTracks();

    //There should be some tracks and vertices in the event
    if(mPrimTracksNum<0 || mPrimTracksNum>10000 ||
       mGlobTracksNum<0 || mGlobTracksNum>50000 ||
       mPrimVertNum<0) {
      return 0;
    }

    //Primary vertex loop
    for(int iVert=0; iVert<mPrimVertNum; iVert++) {

      mHbtEvent = 0;
      //std::cout << "iVert: " << iVert << std::endl;

      //Process only the first primary vertex if it is set
      if(iVert!=0 && mProcessAllVertices==false) continue;

      //std::cout << "Processing vertex: " << iVert << std::endl;

      //Get the current primary vertex
      mPrimVert = mMuDst->primaryVertex(iVert);

      if(mTofHeader) {
	mVpdVz = mTofHeader->vpdVz();
      }
      else {
	mVpdVz = mPrimVert->position().z();
      }

      bool mVertexPassed = IsGoodPrimaryVertex(mPrimVert, mVpdVz);
      if(mVertexPassed == false) {
	//std::cout << "Is bad vertex" << std::endl;
	continue;
      }

      //Create new StHbtEvent
      mHbtEvent = new StHbtEvent();

      //Fill StHbtEvent
      mHbtEvent->SetEventNumber(mMuEvent->eventNumber());
      mHbtEvent->SetRunNumber(mMuEvent->runNumber());
      mHbtEvent->SetCtbMult( (unsigned int) mMuEvent->ctbMultiplicity() );
      mHbtEvent->SetZdcAdcEast( (unsigned int) mMuEvent->zdcAdcAttentuatedSumEast() );
      mHbtEvent->SetZdcAdcWest( (unsigned int) mMuEvent->zdcAdcAttentuatedSumWest() );
      mHbtEvent->SetNumberOfTpcHits( 0 );
      mHbtEvent->SetNumberOfTracks( mMuEvent->eventSummary().numberOfTracks() );
      mHbtEvent->SetNumberOfGoodTracks( mMuEvent->eventSummary().numberOfGoodTracks() );
      mHbtEvent->SetUncorrectedNumberOfPositivePrimaries( mMuEvent->refMultPos() );
      mHbtEvent->SetUncorrectedNumberOfNegativePrimaries( mMuEvent->refMultNeg() );
      mHbtEvent->SetUncorrectedNumberOfPrimaries( mMuEvent->refMult() );   //This is the refMult()
      mHbtEvent->SetReactionPlane(0., 0);
      mHbtEvent->SetReactionPlane(0., 1);
      mHbtEvent->SetReactionPlaneError(0., 0);
      mHbtEvent->SetReactionPlaneError(0., 1);
      mHbtEvent->SetPrimVertPos( mPrimVert->position() );
      mHbtEvent->SetMagneticField( mMuEvent->magneticField() );
      mHbtEvent->SetZdcCoincidenceRate( mMuEvent->runInfo().zdcCoincidenceRate() );
      mHbtEvent->SetBbcCoincidenceRate( mMuEvent->runInfo().bbcCoincidenceRate() );
      mHbtEvent->SetTriggerWord( mMuEvent->l0Trigger().triggerWord() );
      mHbtEvent->SetTriggerActionWord( mMuEvent->l0Trigger().triggerActionWord() );
      mHbtEvent->SetVpdVz( mVpdVz );

      StMuTrack *mPrimTrack = 0;
      StMuTrack *mGlobTrack = 0;
      
      //Fill tracks
      if(mReadTracks) {
	TObjArray* tracks=0;
	switch(mTrackType) {
	case 0: tracks = mMuDst->globalTracks();  break;
	case 1: tracks = mMuDst->primaryTracks(); break;
	default: DEBUGMESSAGE("don't know how to handle this track type");
	}
	
	if(!tracks) {
	  DEBUGMESSAGE("No tracks have been found");
	  delete mHbtEvent;
	  mHbtEvent = 0;
	  return 0;
	}

	int mNTracks = tracks->GetEntries();
	for(int iTrk=0; iTrk<mNTracks; iTrk++) {

	  bool mTrackPassed = false;
	  if(mTrackType == 0) {
	    mGlobTrack = mMuDst->globalTracks(iTrk);
	    mTrackPassed = IsGoodTrack(mGlobTrack, iVert);
	  }
	  else if(mTrackType == 1) {
	    mPrimTrack = mMuDst->primaryTracks(iTrk);
	    mTrackPassed = IsGoodTrack(mPrimTrack, iVert);
	  }

	  if( mTrackPassed == false ) continue;

	  StHbtTrack *trackCopy = new StHbtTrack(mMuDst, (StMuTrack*)tracks->UncheckedAt(iTrk));
	  mHbtEvent->TrackCollection()->push_back(trackCopy);
	} //for(int iTrk=0; iTrk<mNTracks; iTrk++)
      } //if(mReadTracks)

      //Fill V0s
      if(mReadV0s) {
	
	//Should be filled for V0s
      } //if(mReadV0s)

      // Apply event cut, if available
      if(mEventCut) {

	//std::cout << "Why am I here?" << std::endl;
	if (!(mEventCut->Pass(mHbtEvent))){
	  delete mHbtEvent;
	  mHbtEvent = 0;
	  return 0;
	}
      }

      return mHbtEvent;
    } //for(Int_t iVert=0; iVert<mPrimVertNum; iVert++)
    
  } //if (mMuDst && mMuDst->event() )
  else {
    // DTSMaker returned NULL, so we are not OK
    mReaderStatus = 1;
  }

  return mHbtEvent;
}

//_________________
void gregsStHbtMuDstReader::Finish() {
  if (mFinish) {
    for ( int i=0; i<10; i++) {
      cout << "why are you calling the Finish() again  ???????" << endl;
      cout << "are you the stupid chain destructor ???????????" << endl;
    }
  }
  else {
    mFinish = true;
  }
  return;
}

//_________________
void gregsStHbtMuDstReader::setProbabilityPidFile(const char* file) {
  if (mProbabilityPidAlgorithm)
    mProbabilityPidAlgorithm->readParametersFromFile(file);
}

//_________________
bool gregsStHbtMuDstReader::IsGoodTrigger(StMuEvent *event) {

  bool isGood = true;

  if(!mTriggerIdCollection.empty()) {

    isGood = false;
    for(unsigned int i=0; i<mTriggerIdCollection.size(); i++) {

      if(event->triggerIdCollection().nominal().isTrigger(mTriggerIdCollection[i])) {
	isGood = true;
	break;
      }
    }
  } //if(!mTriggerIdCollection.empty())

  return isGood;
}

//_________________
bool gregsStHbtMuDstReader::IsGoodPrimaryVertex(StMuPrimaryVertex *vertex, 
						Float_t vpdVz) {

  bool isGood = false;
  float x = vertex->position().x();
  float y = vertex->position().y();
  float z = vertex->position().z();
  float r = ::sqrt(x*x + y*y);
  float vpdVzDiff = z - vpdVz;
  if( (z >= -80.) && (z <= 80.) && (r <= 3.) &&
      (vpdVzDiff >= -10.) && (vpdVzDiff <= 10.) &&
      (vertex->ranking() > 0) ) {
    isGood = true;
  }

  /*
  if( (::abs(x)<1e-6) && 
      (::abs(y)<1e-6) && 
      (::abs(z)<1e-6) ) {
    isGood = false;
  }
  */

  return isGood;
}
 
//_________________
bool gregsStHbtMuDstReader::IsGoodTrack(StMuTrack *trk, unsigned short vertInd) {

  bool isGood = false;
  if( trk->vertexIndex() == vertInd &&
      trk->momentum().mag() >= 0.1 &&
      trk->momentum().mag() <= 1.7 &&
      trk->eta() >= -1.1 &&
      trk->eta() <= 1.1 &&
      trk->nHits() >= 5 &&
      trk->nHits() <= 50 &&
      trk->flag() >= 0 &&
      trk->flag() <= 1000 &&
      trk->dca(vertInd).perp() >= -0.0001 &&
      trk->dca(vertInd).perp() <= 3. &&
      trk->dcaGlobal(vertInd).perp() >= 0. &&
      trk->dcaGlobal(vertInd).perp() <= 3.) {
    isGood = true;
  }

  return isGood;
}