my_star/StRoot/StHbtMaker/Reader/StHbtPythia6DstReader.cxx

#include "StHbtPythia6DstReader.h"
#include "StarClassLibrary/StPhysicalHelixD.hh"
#include "StarClassLibrary/StTimer.hh"

#include "StHbtMaker/Infrastructure/StHbtEvent.hh"
#include "StHbtMaker/Infrastructure/StHbtTrack.hh"
#include "StHbtMaker/Infrastructure/StHbtV0.hh"
#include "StHbtMaker/Infrastructure/StHbtTrackCollection.hh"
#include "StHbtMaker/Infrastructure/StHbtV0Collection.hh"
#include "StHbtMaker/Base/StHbtEventCut.h"

#include "StThreeVector.hh"

#include "TStreamerInfo.h"
#include "TClonesArray.h"
#include "TLorentzVector.h"
#include "TMath.h"
#include "TMCParticle.h"
#include "TDatabasePDG.h"

ClassImp(StHbtPythia6DstReader)

const Float_t PYHIA_BFIELD = -4.989478911781311;    //From AuAu200 y2011

//_________________
StHbtPythia6DstReader::StHbtPythia6DstReader(TPythia6* dst,
					     int genVersion) {

  mReaderStatus = 0;
  mEventCounter = 0;
  mEventTries = 0;
  mEventIsGood = false;
  mPartId = -999;
  mGenVersion = genVersion;
  mSelPartNum = 0;
  mPythiaDst = dst; //Get the pythia instance
}

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

//_________________
void StHbtPythia6DstReader::clear() {
  mEventIsGood = false;
  mSelPartNum = 0;
}

//_________________
int StHbtPythia6DstReader::Init(const char* ReadWrite, StHbtString& Message){
#ifdef STHBTDEBUG
  std::cout << "StHbtPythia6DstReader:Init -- nothing to do here"
	    << std::endl;
#endif
  return kStOK;
}

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

#ifdef STHBTDEBUG
  std::cout << "StHbtPythia6DstReader::ReturnHbtEvent() " << std::endl;
#endif

  //Usage:
  //mGenVersion = 0 - Clean min. bias generated event
  //mGenVersion = 1 - Two-particle selection

  if (!mPythiaDst) return 0;

  if(mGenVersion<0 || mGenVersion>2) {
    mGenVersion = 0;
  }

  //Create empty StHbtEvent
  StHbtEvent* mHbtEvent = new StHbtEvent;

  if(mGenVersion == 0) {
    RandomGeneration(mHbtEvent);
  }
  else if(mGenVersion == 1) {
    PrimaryTrackGeneration(mHbtEvent);
  }
  else if(mGenVersion == 2) {
    V0Generation(mHbtEvent);
  }

  /*
  std::cout << "StHbtPythia6DstReader::ReturnHbtEvent(): " << std::endl
	    << "Tracks pushed to collection: " << mHbtEvent->TrackCollection()->size() << std::endl
	    << "V0s pushed to collection   : " << mHbtEvent->V0Collection()->size() << std::endl;
  */

  return mHbtEvent;
}

//_________________
void StHbtPythia6DstReader::Finish() {
  std::cout << "StHbtPythia6DstReader:Finish" << std::endl
	    << "Number of good events: " << mEventCounter
	    << "\tNumber of generations: " << mEventTries << std::endl;
}

//________________
Double_t StHbtPythia6DstReader::dedxMean(Double_t mass, Double_t momentum) {

  Double_t dedxMean;
  Double_t tpcDedxGain = 0.174325e-06;
  Double_t tpcDedxOffset = -2.71889; 
  Double_t tpcDedxRise = 776.626;
  
  Double_t gamma = TMath::Sqrt(momentum*momentum/(mass*mass)+1.);
  Double_t beta = TMath::Sqrt(1. - 1./(gamma*gamma));
  Double_t rise = tpcDedxRise* beta*beta*gamma*gamma;      
  if ( beta > 0)
    dedxMean = tpcDedxGain/(beta*beta) * (0.5*TMath::Log(rise) - beta*beta 
					  - tpcDedxOffset);
  else
    dedxMean = 1000.;
  return dedxMean;
}

//_________________
void StHbtPythia6DstReader::RandomGeneration(StHbtEvent *mHbtEvent) {

  //Variables for the StHbtEvent
  StThreeVectorF mPrimVertPos = StThreeVectorF(0., 0., 0.);
  Float_t mVpdVz = 0.;
  Float_t mMagneticField = PYHIA_BFIELD;

  //added by gnigmat
#ifdef STHBTDEBUG
  std::cout << "ReturnHbtEvent: start generation process" << std::endl;
#endif

  //Array of generated particles
  TObjArray *pyParticles = NULL;
  TMCParticle *pyParticle = NULL;
  Int_t mNumParticles = 0;
  Int_t pyPosTrkNum = 0;
  Int_t pyNegTrkNum = 0;
  Int_t pyRefMult = 0;

  TLorentzVector pyTrack, pyPosTrack, pyNegTrack;
  float mPseudoRapidity;
  Int_t pyKF, pyKS, pyCharge;
  short pyLocalCharge = 0;
  Float_t pyTrkNSigmaElectron, pyTrkNSigmaPion;
  Float_t pyTrkNSigmaKaon, pyTrkNSigmaProton;
  Float_t pyTrkPidProbElectron, pyTrkPidProbPion;
  Float_t pyTrkPidProbKaon, pyTrkPidProbProton;
  Float_t pyTrkDCAxy, pyTrkDCAz, pyTrkMass;

  //Start event generation here
  do {

    clear();
    mPythiaDst->GenerateEvent();
    pyParticles = mPythiaDst->GetListOfParticles();
    mEventTries++;
    if(pyParticles) {
      mEventIsGood = true;
    }
  } while(!mEventIsGood);

  mEventCounter++;

  mNumParticles = pyParticles->GetEntries();

  //mPythiaDst->Pylist(1);

  //Number of particles in acceptance
  pyRefMult = 0;

  //Main particle loop
  for(Int_t iTrk=0; iTrk<mNumParticles; iTrk++) {

    pyKF = 0;
    pyKS = 0;
    pyCharge = 0;

    pyParticle = (TMCParticle*)pyParticles->At(iTrk);
    //Use final-state particles
    if(!pyParticle || pyParticle->GetKS()!=1) continue;

    pyKF = pyParticle->GetKF();
    pyKS = pyParticle->GetKS();
    pyCharge = TDatabasePDG::Instance()->GetParticle(pyKF)->Charge();

    pyTrack.SetPxPyPzE(pyParticle->GetPx(),
		       pyParticle->GetPy(),
		       pyParticle->GetPz(),
		       pyParticle->GetEnergy());
    
    if(pyTrack.P() < 0.1)  continue;

    mPseudoRapidity = 0.5*TMath::Log( (pyTrack.P() + pyTrack.Pz()) /
				      (pyTrack.P() - pyTrack.Pz()) );

    //Calculate multiplicity
    if( pyTrack.P() >= 0.1 && TMath::Abs(mPseudoRapidity)<=0.5 &&
	pyCharge!=0 ) {

      pyRefMult++;

      if(pyCharge>0) {
	pyPosTrkNum++;
	pyLocalCharge = 1;
      }
      else if( pyCharge<0) {
	pyNegTrkNum++;
	pyLocalCharge = -1;
      }
    }

    //Do not use particles that not pass next cut
    if(pyTrack.P()<0.1 || TMath::Abs(mPseudoRapidity)>1.1) continue;

    //Primary tracks
    if(TMath::Abs(pyKF)!=211 && TMath::Abs(pyKF)!=321 &&
       TMath::Abs(pyKF)!=2212) continue;

    //Factor 10 converts mm to cm
    pyTrkDCAxy = TMath::Sqrt(pyParticle->GetVx()*pyParticle->GetVx() +
			     pyParticle->GetVy()*pyParticle->GetVy()) / 10;
    pyTrkDCAz = pyParticle->GetVz() / 10;

    //Primary track cut (5 cm)
    if(pyTrkDCAxy > 5. ||
       pyTrkDCAz > 5.) continue;

    if(pyKF>0) {
      pyLocalCharge = +1;
    }
    else if(pyKF<0) {
      pyLocalCharge = -1;
    }
    pyTrkMass = pyParticle->GetMass();

    switch(TMath::Abs(pyKF)) {
    case 211:
      pyTrkNSigmaElectron = -999.;
      pyTrkNSigmaPion = 0.;
      pyTrkNSigmaKaon = -999.;
      pyTrkNSigmaProton = -999.;
      pyTrkPidProbElectron = -999.;
      pyTrkPidProbPion = 0.;
      pyTrkPidProbKaon = -999.;
      pyTrkPidProbProton = -999.;
      break;
    case 321:
      pyTrkNSigmaElectron = -999.;
      pyTrkNSigmaPion = -999.;
      pyTrkNSigmaKaon = 0.;
      pyTrkNSigmaProton = -999.;
      pyTrkPidProbElectron = -999.;
      pyTrkPidProbPion = -999.;
      pyTrkPidProbKaon = 0.;
      pyTrkPidProbProton = -999.;
      break;
    case 2212:
      pyTrkNSigmaElectron = -999.;
      pyTrkNSigmaPion = -999.;
      pyTrkNSigmaKaon = -999.;
      pyTrkNSigmaProton = 0.;
      pyTrkPidProbElectron = -999.;
      pyTrkPidProbPion = -999.;
      pyTrkPidProbKaon = -999.;
      pyTrkPidProbProton = 0.;
      break;
    default:
      pyTrkNSigmaElectron = -999.;
      pyTrkNSigmaPion = -999.;
      pyTrkNSigmaKaon = -999.;
      pyTrkNSigmaProton = -999.;
      pyTrkPidProbElectron = -999.;
      pyTrkPidProbPion = -999.;
      pyTrkPidProbKaon = -999.;
      pyTrkPidProbProton = -999.;
    };


    StHbtTrack *track = new StHbtTrack;

    //0-global, 1-primary
    //track->SetTrackType(1);
    track->SetFlag(100);
    //track->SetCharge(pyLocalCharge);
    track->SetNHits(pyLocalCharge * 45);
    //track->SetNHitsFit(45);
    track->SetNHitsPossible(45);
    track->SetNSigmaElectron(pyTrkNSigmaElectron);
    track->SetNSigmaPion(pyTrkNSigmaPion);
    track->SetNSigmaKaon(pyTrkNSigmaKaon);
    track->SetNSigmaProton(pyTrkNSigmaProton);
    track->SetPidProbElectron(pyTrkPidProbElectron);
    track->SetPidProbPion(pyTrkPidProbPion);
    track->SetPidProbKaon(pyTrkPidProbKaon);
    track->SetPidProbProton(pyTrkPidProbProton);
    track->SetdEdx( dedxMean(pyTrkMass, pyTrack.P() ) );
    //track->SetDCAxy(pyTrkDCAxy);
    //track->SetDCAz(pyTrkDCAz);
    track->SetDCAxyGlobal(pyTrkDCAxy);
    track->SetDCAzGlobal(pyTrkDCAz);
    track->SetChi2(0.);
    //track->SetChiSquaredXY(1.);
    //track->SetChiSquaredZ(1.);
    track->SetP(StThreeVectorF(pyTrack.Px(), pyTrack.Py(), pyTrack.Pz()));
    track->SetPGlobal(StThreeVectorF(pyTrack.Px(),pyTrack.Py(),pyTrack.Pz()));
    StPhysicalHelixD mHelix = StPhysicalHelixD(StThreeVectorF(pyTrack.Px(),
							      pyTrack.Py(),
							      pyTrack.Pz()),
					       StThreeVectorF(0,0,0),
					       mMagneticField*kilogauss,
					       pyLocalCharge);
    track->SetHelix(mHelix);
    track->SetHelixGlobal(mHelix);
    track->SetTopologyMap(0, 0.);
    track->SetTopologyMap(1, 0.);
    short pyId=iTrk;
    track->SetTrackId(pyId);
    track->SetTofBeta(pyTrack.P()/pyParticle->GetEnergy());
    track->SetHiddenInfo(0);

    mHbtEvent->TrackCollection()->push_back(track);

  } //for (int iTrk=0; iTrk<mNumParticles; iTrk++)

  mHbtEvent->SetEventNumber(mEventCounter);
  mHbtEvent->SetRunNumber(mEventCounter);
  mHbtEvent->SetCtbMult(0);
  mHbtEvent->SetZdcAdcEast(0);
  mHbtEvent->SetZdcAdcWest(0);
  mHbtEvent->SetNumberOfTpcHits(0);
  mHbtEvent->SetNumberOfTracks(pyRefMult);
  mHbtEvent->SetNumberOfGoodTracks(pyRefMult);
  mHbtEvent->SetReactionPlane(0.);
  mHbtEvent->SetReactionPlaneError(0.);
  mHbtEvent->SetReactionPlaneSubEventDifference(0.);
  mHbtEvent->SetPrimVertPos(mPrimVertPos);
  mHbtEvent->SetUncorrectedNumberOfPositivePrimaries(pyPosTrkNum);
  mHbtEvent->SetUncorrectedNumberOfNegativePrimaries(pyNegTrkNum);
  mHbtEvent->SetUncorrectedNumberOfPrimaries(pyRefMult);
  mHbtEvent->SetMagneticField(mMagneticField);
  mHbtEvent->SetVpdVz(mVpdVz);
}

//_________________
void StHbtPythia6DstReader::PrimaryTrackGeneration(StHbtEvent *mHbtEvent) {

  //Variables for the StHbtEvent
  StThreeVectorF mPrimVertPos = StThreeVectorF(0., 0., 0.);
  Float_t mVpdVz = 0.;
  Float_t mMagneticField = PYHIA_BFIELD;

  //added by gnigmat
#ifdef STHBTDEBUG
  std::cout << "ReturnHbtEvent: start generation process" << std::endl;
#endif

  //Array of generated particles
  TObjArray *pyParticles = NULL;
  TMCParticle *pyParticle = NULL;
  Int_t mNumParticles = 0;
  Int_t pyPosTrkNum = 0;
  Int_t pyNegTrkNum = 0;
  Int_t pyRefMult = 0;

  TLorentzVector pyTrack, pyPosTrack, pyNegTrack;
  float mPseudoRapidity;
  Int_t pyKF, pyKS, pyCharge;
  short pyLocalCharge = 0;
  Float_t pyTrkNSigmaElectron, pyTrkNSigmaPion;
  Float_t pyTrkNSigmaKaon, pyTrkNSigmaProton;
  Float_t pyTrkPidProbElectron, pyTrkPidProbPion;
  Float_t pyTrkPidProbKaon, pyTrkPidProbProton;
  Float_t pyTrkDCAxy, pyTrkDCAz, pyTrkMass;

  //Start event generation here
  do {

    clear();
    mPythiaDst->GenerateEvent();
    pyParticles = mPythiaDst->GetListOfParticles();
    mEventTries++;

    if(!pyParticles) continue;

    mNumParticles = pyParticles->GetEntries();

    //mPythiaDst->Pylist(1);

    //Number of particles in acceptance
    pyRefMult = 0;

    //Main particle loop
    for(Int_t iTrk=0; iTrk<mNumParticles; iTrk++) {

      pyKF = 0;
      pyKS = 0;
      pyCharge = 0;

      pyParticle = (TMCParticle*)pyParticles->At(iTrk);
      //Use final-state particles
      if(!pyParticle || pyParticle->GetKS()!=1) continue;

      pyKF = pyParticle->GetKF();
      pyKS = pyParticle->GetKS();
      pyCharge = TDatabasePDG::Instance()->GetParticle(pyKF)->Charge();

      pyTrack.SetPxPyPzE(pyParticle->GetPx(),
			 pyParticle->GetPy(),
			 pyParticle->GetPz(),
			 pyParticle->GetEnergy());
    
      if(pyTrack.P() < 0.1)  continue;

      mPseudoRapidity = 0.5*TMath::Log( (pyTrack.P() + pyTrack.Pz()) /
					(pyTrack.P() - pyTrack.Pz()) );

      //Calculate multiplicity
      if( pyTrack.P() >= 0.1 && TMath::Abs(mPseudoRapidity)<=0.5 &&
	  pyCharge!=0 ) {

	pyRefMult++;

	if(pyCharge>0) {
	  pyPosTrkNum++;
	  pyLocalCharge = 1;
	}
	else if( pyCharge<0) {
	  pyNegTrkNum++;
	  pyLocalCharge = -1;
	}
      }

      //Do not use particles that not pass next cut
      if(pyTrack.P()<0.1 || TMath::Abs(mPseudoRapidity)>1.1) continue;

      //Primary tracks
      if(TMath::Abs(pyKF)!=211 && TMath::Abs(pyKF)!=321 &&
	 TMath::Abs(pyKF)!=2212) continue;

      if(TMath::Abs(pyKF) != TMath::Abs(mPartId)) continue;

      //Factor 10 converts mm to cm
      pyTrkDCAxy = TMath::Sqrt(pyParticle->GetVx()*pyParticle->GetVx() +
			       pyParticle->GetVy()*pyParticle->GetVy()) / 10;
      pyTrkDCAz = pyParticle->GetVz() / 10;

      //Primary track cut (5 cm)
      if(pyTrkDCAxy > 5. ||
	 pyTrkDCAz > 5.) continue;

      mSelPartNum++;

      pyTrkMass = pyParticle->GetMass();

      switch(TMath::Abs(pyKF)) {
      case 211:
	pyTrkNSigmaElectron = -999.;
	pyTrkNSigmaPion = 0.;
	pyTrkNSigmaKaon = -999.;
	pyTrkNSigmaProton = -999.;
	pyTrkPidProbElectron = -999.;
	pyTrkPidProbPion = 0.;
	pyTrkPidProbKaon = -999.;
	pyTrkPidProbProton = -999.;
	break;
      case 321:
	pyTrkNSigmaElectron = -999.;
	pyTrkNSigmaPion = -999.;
	pyTrkNSigmaKaon = 0.;
	pyTrkNSigmaProton = -999.;
	pyTrkPidProbElectron = -999.;
	pyTrkPidProbPion = -999.;
	pyTrkPidProbKaon = 0.;
	pyTrkPidProbProton = -999.;
	break;
      case 2212:
	pyTrkNSigmaElectron = -999.;
	pyTrkNSigmaPion = -999.;
	pyTrkNSigmaKaon = -999.;
	pyTrkNSigmaProton = 0.;
	pyTrkPidProbElectron = -999.;
	pyTrkPidProbPion = -999.;
	pyTrkPidProbKaon = -999.;
	pyTrkPidProbProton = 0.;
	break;
      default:
	pyTrkNSigmaElectron = -999.;
	pyTrkNSigmaPion = -999.;
	pyTrkNSigmaKaon = -999.;
	pyTrkNSigmaProton = -999.;
	pyTrkPidProbElectron = -999.;
	pyTrkPidProbPion = -999.;
	pyTrkPidProbKaon = -999.;
	pyTrkPidProbProton = -999.;
      };

      StHbtTrack *track = new StHbtTrack;

      //0-global, 1-primary
      //track->SetTrackType(1);
      track->SetFlag(300);
      //track->SetCharge(pyLocalCharge);
      track->SetNHits(pyLocalCharge * 45);
      //track->SetNHitsFit(45);
      track->SetNHitsPossible(45);
      track->SetNSigmaElectron(pyTrkNSigmaElectron);
      track->SetNSigmaPion(pyTrkNSigmaPion);
      track->SetNSigmaKaon(pyTrkNSigmaKaon);
      track->SetNSigmaProton(pyTrkNSigmaProton);
      track->SetPidProbElectron(pyTrkPidProbElectron);
      track->SetPidProbPion(pyTrkPidProbPion);
      track->SetPidProbKaon(pyTrkPidProbKaon);
      track->SetPidProbProton(pyTrkPidProbProton);
      track->SetdEdx( dedxMean(pyTrkMass, pyTrack.P()) );
      //track->SetDCAxy(pyTrkDCAxy);
      //track->SetDCAz(pyTrkDCAz);
      track->SetDCAxyGlobal(pyTrkDCAxy);
      track->SetDCAzGlobal(pyTrkDCAz);
      track->SetChi2(1.);
      //track->SetChiSquaredXY(1.);
      //track->SetChiSquaredZ(1.);
      track->SetP(StThreeVectorF(pyTrack.Px(), pyTrack.Py(), pyTrack.Pz()));
      track->SetPGlobal(StThreeVectorF(pyTrack.Px(),pyTrack.Py(),pyTrack.Pz()));
      StPhysicalHelixD mHelix = StPhysicalHelixD(StThreeVectorF(pyTrack.Px(),
								pyTrack.Py(),
								pyTrack.Pz()),
						 StThreeVectorF(0,0,0),
						 mMagneticField*kilogauss,
						 pyLocalCharge);
      track->SetHelix(mHelix);
      track->SetHelixGlobal(mHelix);
      track->SetTopologyMap(0, 0.);
      track->SetTopologyMap(1, 0.);
      track->SetTrackId(iTrk);
      track->SetTofBeta(pyTrack.P()/pyParticle->GetEnergy());
      track->SetHiddenInfo(0);

      mHbtEvent->TrackCollection()->push_back(track);

    } //for (int iTrk=0; iTrk<mNumParticles; iTrk++)

    if(mSelPartNum>=1) {
      mEventIsGood = true;
    }
  } while(!mEventIsGood);

  mEventCounter++;

  mHbtEvent->SetEventNumber(mEventCounter);
  mHbtEvent->SetRunNumber(mEventCounter);
  mHbtEvent->SetCtbMult(0);
  mHbtEvent->SetZdcAdcEast(0);
  mHbtEvent->SetZdcAdcWest(0);
  mHbtEvent->SetNumberOfTpcHits(0);
  mHbtEvent->SetNumberOfTracks(pyRefMult);
  mHbtEvent->SetNumberOfGoodTracks(pyRefMult);
  mHbtEvent->SetReactionPlane(0.);
  mHbtEvent->SetReactionPlaneError(0.);
  mHbtEvent->SetReactionPlaneSubEventDifference(0.);
  mHbtEvent->SetPrimVertPos(mPrimVertPos);
  mHbtEvent->SetUncorrectedNumberOfPositivePrimaries(pyPosTrkNum);
  mHbtEvent->SetUncorrectedNumberOfNegativePrimaries(pyNegTrkNum);
  mHbtEvent->SetUncorrectedNumberOfPrimaries(pyRefMult);
  mHbtEvent->SetMagneticField(mMagneticField);
  mHbtEvent->SetVpdVz(mVpdVz);
}

//_________________
void StHbtPythia6DstReader::V0Generation(StHbtEvent *mHbtEvent) {

  //Variables for the StHbtEvent
  StThreeVectorF mPrimVertPos = StThreeVectorF(0., 0., 0.);
  Float_t mVpdVz = 0.;
  Float_t mMagneticField = PYHIA_BFIELD;

  //added by gnigmat
#ifdef STHBTDEBUG
  std::cout << "ReturnHbtEvent: start generation process" << std::endl;
#endif

  //Array of generated particles
  TObjArray *pyParticles = NULL;
  TMCParticle *pyParticle = NULL;
  TMCParticle *pyPosParticle, *pyNegParticle;
  Int_t mNumParticles = 0;
  Int_t pyPosTrkNum = 0;
  Int_t pyNegTrkNum = 0;
  Int_t pyRefMult = 0;

  TLorentzVector pyTrack, pyPosTrack, pyNegTrack;
  float mPseudoRapidity;
  Int_t pyKF, pyKS, pyCharge;
  short pyLocalCharge = 0;
  Float_t pyTrkNSigmaElectron, pyTrkNSigmaPion;
  Float_t pyTrkNSigmaKaon, pyTrkNSigmaProton;
  Float_t pyTrkDCAxy, pyTrkDCAz, pyTrkMass;
  Int_t   mV0Num = 0;

  //Start event generation here
  do {

    clear();
    mPythiaDst->GenerateEvent();
    pyParticles = mPythiaDst->GetListOfParticles();
    mEventTries++;

    if(!pyParticles) continue;

    mNumParticles = pyParticles->GetEntries();

    //mPythiaDst->Pylist(1);

    //Initialize counters
    mV0Num = 0;
    pyRefMult = 0;

    //Main particle loop
    for(Int_t iTrk=0; iTrk<mNumParticles; iTrk++) {

      pyKF = 0;
      pyKS = 0;
      pyCharge = 0;

      pyParticle = (TMCParticle*)pyParticles->At(iTrk);

      //Use final-state particles
      if(!pyParticle) continue;

      pyKF = pyParticle->GetKF();
      pyKS = pyParticle->GetKS();
      pyCharge = TDatabasePDG::Instance()->GetParticle(pyKF)->Charge();

      pyTrack.SetPxPyPzE(pyParticle->GetPx(),
			 pyParticle->GetPy(),
			 pyParticle->GetPz(),
			 pyParticle->GetEnergy());

      //Calculate stable particles for refmult
      if(pyParticle->GetKS()==1) {
    
	if(pyTrack.P() < 0.1)  continue;

	mPseudoRapidity = 0.5*TMath::Log( (pyTrack.P() + pyTrack.Pz()) /
					  (pyTrack.P() - pyTrack.Pz()) );

	//Calculate multiplicity
	if( pyTrack.P() >= 0.1 && TMath::Abs(mPseudoRapidity)<=0.5 &&
	    pyCharge!=0 ) {

	  pyRefMult++;

	  if(pyCharge>0) {
	    pyPosTrkNum++;
	    pyLocalCharge = 1;
	  }
	  else if( pyCharge<0) {
	    pyNegTrkNum++;
	    pyLocalCharge = -1;
	  }
	}
      } //if(pyParticle->GetKS()==1)

      //Select V0s here
      if(pyKF==333) {

	//Factor 10 converts mm to cm
	pyTrkDCAxy = TMath::Sqrt(pyParticle->GetVx()*pyParticle->GetVx() +
				 pyParticle->GetVy()*pyParticle->GetVy()) / 10;
	pyTrkDCAz = pyParticle->GetVz() / 10;
	Float_t pyTrkDCA2PrimVtx = TMath::Sqrt(pyParticle->GetVx()*pyParticle->GetVx() +
					       pyParticle->GetVy()*pyParticle->GetVy() +
					       pyParticle->GetVz()*pyParticle->GetVz()) / 10;

	pyTrkMass = pyParticle->GetMass();

	Int_t mFirstDaughterId = pyParticle->GetFirstChild();
	Int_t mSecondDaugherId = pyParticle->GetLastChild();
	Int_t pyPosCharge = 0;
	Int_t pyNegCharge = 0;
	Int_t pyPosKF = 0;
	Int_t pyNegKF = 0;
	Int_t pyPosId;
	Int_t pyNegId;

	pyPosParticle = (TMCParticle*)pyParticles->At(mFirstDaughterId-1);
	pyNegParticle = (TMCParticle*)pyParticles->At(mSecondDaugherId-1);
	pyPosKF = pyPosParticle->GetKF();
	pyNegKF = pyNegParticle->GetKF();
	pyPosTrack.SetPxPyPzE(pyPosParticle->GetPx(), pyPosParticle->GetPy(),
			      pyPosParticle->GetPz(), pyPosParticle->GetEnergy());
	pyNegTrack.SetPxPyPzE(pyNegParticle->GetPx(), pyNegParticle->GetPy(),
			      pyNegParticle->GetPz(), pyNegParticle->GetEnergy());

	//phi->K+K- only.
	if(TMath::Abs(pyPosKF)!=321 || TMath::Abs(pyNegKF)!=321 ||
	   pyPosTrack.P()<0.1 || pyPosTrack.PseudoRapidity()>1.1 ||
	   pyNegTrack.P()<0.1 || pyNegTrack.PseudoRapidity()>1.1) continue;

	pyPosCharge = TDatabasePDG::Instance()->GetParticle(pyPosKF)->Charge();
	pyNegCharge = TDatabasePDG::Instance()->GetParticle(pyNegKF)->Charge();
	if(pyPosCharge>0 && pyNegCharge<0) {

	  pyPosId = mFirstDaughterId;
	  pyNegId = mSecondDaugherId;
	}
	else if(pyPosCharge<0 && pyNegCharge>0) { //Switch positive and negative daughters
	  pyPosParticle  = (TMCParticle*)pyParticles->At(mSecondDaugherId-1);
	  pyNegParticle = (TMCParticle*)pyParticles->At(mFirstDaughterId-1);
	  pyPosKF = pyPosParticle->GetKF();
	  pyNegKF = pyNegParticle->GetKF();
	  pyPosCharge = TDatabasePDG::Instance()->GetParticle(pyPosKF)->Charge();
	  pyNegCharge = TDatabasePDG::Instance()->GetParticle(pyNegKF)->Charge();
	  pyPosTrack.SetPxPyPzE(pyPosParticle->GetPx(), pyPosParticle->GetPy(),
				pyPosParticle->GetPz(), pyPosParticle->GetEnergy());
	  pyNegTrack.SetPxPyPzE(pyNegParticle->GetPx(), pyNegParticle->GetPy(),
				pyNegParticle->GetPz(), pyNegParticle->GetEnergy());
	  pyNegId = mFirstDaughterId;
	  pyPosId = mSecondDaugherId;
	}

#ifdef STHBTDEBUG
	std::cout << "V0 is found: " << std::endl
		  << Form("KF: %5d mass: %5.3f mom: %5.3f PosChildId: %d NegChildId: %d",
			  pyKF, pyParticle->GetMass(), pyTrack.P(), pyPosId, pyNegId) << std::endl
		  << Form("PosChild KF: %5d mass: %5.3f mom: %5.3f charge: %d",
			  pyPosKF, pyPosParticle->GetMass(), pyPosTrack.P(), pyPosCharge) << std::endl
		  << Form("NegChild KF: %5d mass: %5.3f mom: %5.3f charge: %d",
			  pyNegKF, pyNegParticle->GetMass(), pyNegTrack.P(), pyNegCharge) << std::endl;
#endif

	Float_t mDecayLengthV0 = TMath::Sqrt( pyParticle->GetVx()*pyParticle->GetVx() +
					      pyParticle->GetVy()*pyParticle->GetVy() +
					      pyParticle->GetVz()*pyParticle->GetVz() )/10;

	//Increment number of V0s
	mV0Num++;

	//Create an instance ov StHbtV0, fill it and add to mHbtEvent
	StHbtV0 *mHbtV0 = new StHbtV0;

	mHbtV0->SetdecayLengthV0(mDecayLengthV0);
	mHbtV0->SetdecayVertexV0( StHbtThreeVector(pyParticle->GetVx() * 0.1,
						   pyParticle->GetVy() * 0.1,
						   pyParticle->GetVz() * 0.1) );
	mHbtV0->SetdcaV0Daughters(0.);
	mHbtV0->SetmomV0( StHbtThreeVector(pyParticle->GetPx(), pyParticle->GetPy(),
					   pyParticle->GetPz()) );
	mHbtV0->SetprimaryVertex( StHbtThreeVector(0., 0., 0.) );
	mHbtV0->SetHiddenInfo(0);

	//Positively charged daughter
	switch(TMath::Abs(pyPosKF)) {
	case 211:
	  pyTrkNSigmaElectron = -999.;
	  pyTrkNSigmaPion = 0.;
	  pyTrkNSigmaKaon = -999.;
	  pyTrkNSigmaProton = -999.;
	  break;
	case 321:
	  pyTrkNSigmaElectron = -999.;
	  pyTrkNSigmaPion = -999.;
	  pyTrkNSigmaKaon = 0.;
	  pyTrkNSigmaProton = -999.;
	  break;
	case 2212:
	  pyTrkNSigmaElectron = -999.;
	  pyTrkNSigmaPion = -999.;
	  pyTrkNSigmaKaon = -999.;
	  pyTrkNSigmaProton = 0.;
	  break;
	default:
	  pyTrkNSigmaElectron = -999.;
	  pyTrkNSigmaPion = -999.;
	  pyTrkNSigmaKaon = -999.;
	  pyTrkNSigmaProton = -999.;
	};

	mHbtV0->SetdcaPosToPrimVertex(pyTrkDCA2PrimVtx);
	mHbtV0->SetmomPos( StHbtThreeVector(pyPosTrack.Px(), pyPosTrack.Py(), pyPosTrack.Pz()) );
	mHbtV0->SettpcHitsPos(45);
	//mHbtV0->SetnHitsFitPos(45);
	mHbtV0->SetnHitsPossiblePos(45);
	mHbtV0->SetidPos(pyPosId);
	mHbtV0->SetkeyPos(pyPosId);
	mHbtV0->SetdedxPos( dedxMean(pyPosTrack.P(), pyPosParticle->GetMass()) );
	mHbtV0->SetlendedxPos(150.);
	mHbtV0->SetnSigmaElectronPos(pyTrkNSigmaElectron);
	mHbtV0->SetnSigmaPionPos(pyTrkNSigmaPion);
	mHbtV0->SetnSigmaKaonPos(pyTrkNSigmaKaon);
	mHbtV0->SetnSigmaProtonPos(pyTrkNSigmaProton);
	mHbtV0->SettofBetaPos(pyPosTrack.P()/pyPosTrack.Energy());
	StPhysicalHelixD mHelixPos = StPhysicalHelixD(StThreeVectorF(pyPosTrack.Px(),
								     pyPosTrack.Py(),
								     pyPosTrack.Pz()),
						      StThreeVectorF(0.,0.,0.),
						      mMagneticField*kilogauss,
						      pyPosCharge);
	mHbtV0->SetHelixPos(mHelixPos);
	mHbtV0->SetTrackTopologyMapPos(0, 0.);
	mHbtV0->SetTrackTopologyMapPos(1, 0.);

	//Negatively charged daughter
	switch(TMath::Abs(pyNegKF)) {
	case 211:
	  pyTrkNSigmaElectron = -999.;
	  pyTrkNSigmaPion = 0.;
	  pyTrkNSigmaKaon = -999.;
	  pyTrkNSigmaProton = -999.;
	  break;
	case 321:
	  pyTrkNSigmaElectron = -999.;
	  pyTrkNSigmaPion = -999.;
	  pyTrkNSigmaKaon = 0.;
	  pyTrkNSigmaProton = -999.;
	  break;
	case 2212:
	  pyTrkNSigmaElectron = -999.;
	  pyTrkNSigmaPion = -999.;
	  pyTrkNSigmaKaon = -999.;
	  pyTrkNSigmaProton = 0.;
	  break;
	default:
	  pyTrkNSigmaElectron = -999.;
	  pyTrkNSigmaPion = -999.;
	  pyTrkNSigmaKaon = -999.;
	  pyTrkNSigmaProton = -999.;
	};

	mHbtV0->SetdcaNegToPrimVertex(pyTrkDCA2PrimVtx);
	mHbtV0->SetmomNeg( StHbtThreeVector(pyNegTrack.Px(), pyNegTrack.Py(), pyNegTrack.Pz()) );
	mHbtV0->SettpcHitsNeg(45);
	//mHbtV0->SetnHitsFitNeg(45);
	mHbtV0->SetnHitsPossibleNeg(45);
	mHbtV0->SetidNeg(pyNegId);
	mHbtV0->SetkeyNeg(pyNegId);
	mHbtV0->SetdedxNeg( dedxMean(pyNegTrack.P(), pyNegParticle->GetMass()) );
	mHbtV0->SetlendedxNeg(150.);
	mHbtV0->SetnSigmaElectronNeg(pyTrkNSigmaElectron);
	mHbtV0->SetnSigmaPionNeg(pyTrkNSigmaPion);
	mHbtV0->SetnSigmaKaonNeg(pyTrkNSigmaKaon);
	mHbtV0->SetnSigmaProtonNeg(pyTrkNSigmaProton);
	mHbtV0->SettofBetaNeg(pyNegTrack.P()/pyNegTrack.Energy());
	StPhysicalHelixD mHelixNeg = StPhysicalHelixD(StThreeVectorF(pyNegTrack.Px(),
								     pyNegTrack.Py(),
								     pyNegTrack.Pz()),
						      StThreeVectorF(0.,0.,0.),
						      mMagneticField*kilogauss,
						      pyNegCharge);
	mHbtV0->SetHelixNeg(mHelixNeg);
	mHbtV0->SetTrackTopologyMapNeg(0, 0.);
	mHbtV0->SetTrackTopologyMapNeg(1, 0.);

	mHbtV0->UpdateV0();

	mHbtEvent->V0Collection()->push_back(mHbtV0);
      } //if(pyKF==333) 
    } //for (int iTrk=0; iTrk<mNumParticles; iTrk++)
    
    if(mV0Num>=1) {
      mEventIsGood = true;
    }
  } while(!mEventIsGood);
  
  mEventCounter++;

  mHbtEvent->SetEventNumber(mEventCounter);
  mHbtEvent->SetRunNumber(mEventCounter);
  mHbtEvent->SetCtbMult(0);
  mHbtEvent->SetZdcAdcEast(45000);
  mHbtEvent->SetZdcAdcWest(45000);
  mHbtEvent->SetNumberOfTpcHits(450);
  mHbtEvent->SetNumberOfTracks(pyRefMult);
  mHbtEvent->SetNumberOfGoodTracks(pyRefMult);
  mHbtEvent->SetReactionPlane(0.);
  mHbtEvent->SetReactionPlaneError(0.);
  mHbtEvent->SetReactionPlaneSubEventDifference(0.);
  mHbtEvent->SetPrimVertPos(mPrimVertPos);
  mHbtEvent->SetUncorrectedNumberOfPositivePrimaries(pyPosTrkNum);
  mHbtEvent->SetUncorrectedNumberOfNegativePrimaries(pyNegTrkNum);
  mHbtEvent->SetUncorrectedNumberOfPrimaries(pyRefMult);
  mHbtEvent->SetMagneticField(mMagneticField);
  mHbtEvent->SetVpdVz(mVpdVz);
}