StHighPtFilter.cxx

// @(#)root/eg:$Id: StHighPtFilter.cxx,v 1.1 2011/01/25 20:06:53 jwebb Exp $
// Author: Victor Perev  17/03/2009

//______________________________________________________________________________
#include "stdlib.h"
#include "math.h"
#include <iostream>
#include <fstream>

#include "StHighPtFilter.h"
#include "StGenParticle.h"

using namespace std;


static StHighPtFilter particleFilter;


StHighPtFilter::StHighPtFilter():StMCFilter("highpt")
{
  
  mPtCut  = 2.5;
  mEtaCut = 1.2;
  readConfig();
  std::cout << "======================================================" << std::endl;
  std::cout << "highpt hadron filter initialized" << std::endl;
  std::cout << "PtCUT:  " << mPtCut  << std::endl;
  std::cout << "EtaCUT: " << mEtaCut << std::endl;
  std::cout << "======================================================" << std::endl;
};

void StHighPtFilter::readConfig()
{
  ifstream cfile("highpt.cnf");
  while (1) {
    string   key;
    double   value;
    cfile >> key >> value;
    if ( !cfile.good() ) break; // EOF or no file
    parseConfig( key, value );
  }
}

void StHighPtFilter::parseConfig(string key, double value)
{
  if ( key=="PtCUT"  ){ mPtCut  = value; }
  if ( key=="EtaCUT" ){ mEtaCut = value; }
  return;
}


//______________________________________________________________________________
int StHighPtFilter::RejectEG(const StGenParticleMaster &ptl) const
{
  //ptl.Print("************** In RejectEG ************** ");
  //cout<<"********************** In RejectEG ******************"<<endl;
  return 0;
}
//______________________________________________________________________________
int StHighPtFilter::RejectGT(const StGenParticleMaster &ptl) const
{

  const StGenParticle *tk=0;
  int n = ptl.Size();

  //
  // Loop over all particles and accept events with high pT
  // charged hadrons (pi+/-, K+/- and p/pbar) within |y|<1.2
  // and with pT > 2.5 GeV.  The |y|<1.2 cut is consistent
  // with  vertex cuts of |z|<100 cm.
  //
  float maxPt = 0.; 

  for (int i=0;i<n;i++) {

    tk = ptl(i); if (!tk)       continue;
    if (tk->GetStatusCode()!=1) continue;
    if( fabs(tk->Eta())> mEtaCut ) continue;  //Jason's suggestion@Jan13th2011

    if ( tk->Pt() < mPtCut ) continue;

    if( (fabs(tk->GetPdgCode())==211 || 
         fabs(tk->GetPdgCode())==321 || 
         fabs(tk->GetPdgCode())==2212) ) return 0; // accept on any charged hadron with pT > threshold

  }

  return 1; // No charged tracks with pT > threshold


}
//______________________________________________________________________________
int StHighPtFilter::RejectGE(const StGenParticleMaster &ptl) const
{
  //ptl.Print("************** In RejectGE ************** ");
  return 0;
}