StMuL3Filter.cxx

/***************************************************************************
 *
 * $Id: StMuL3Filter.cxx,v 1.5 2003/10/20 19:50:13 perev Exp $
 * Author: Frank Laue, BNL, laue@bnl.gov
 *
 ***************************************************************************/

#include "StMuL3Filter.h"
#include "StMuDebug.h"
#include "StMuException.hh"
#include "StEvent/StTrack.h"
#include "StEvent/StTrackGeometry.h"
#include "StEvent/StTrackDetectorInfo.h"
#include "StEvent/StContainers.h"
#include "StEvent/StDedxPidTraits.h"


#include "StarClassLibrary/BetheBloch.h"

ClassImp(StMuL3Filter)

bool StMuL3Filter::accept( const StEvent* e) { cout << "StMuL3Filter::accept( const StEvent* e) not overwritten, returning true" << endl; return true;}
bool StMuL3Filter::accept( const StV0Vertex* v) { cout << "StMuL3Filter::accept(const StV0Vertex* v) not overwritten, returning true" << endl; return true;}
bool StMuL3Filter::accept( const StXiVertex* x) { cout << "StMuL3Filter::accept(const StXiVertex* x) not overwritten, returning true" << endl; return true;}
bool StMuL3Filter::accept( const StKinkVertex* k) { cout << "StMuL3Filter::accept(const StKinkVertex* k) not overwritten, returning true" << endl; return true;}
bool StMuL3Filter::accept( const StV0MuDst* v) { cout << "StMuL3Filter::accept(const StV0MuDst* v) not overwritten, returning true" << endl; return true;}
bool StMuL3Filter::accept( const StXiMuDst* x) { cout << "StMuL3Filter::accept(const StXiMuDst* x) not overwritten, returning true" << endl; return true;}
bool StMuL3Filter::accept( const StKinkMuDst* k) { cout << "StMuL3Filter::accept(const StKinkMuDst* k) not overwritten, returning true" << endl; return true;}
 

StMuL3Filter::StMuL3Filter()  {
  DEBUGMESSAGE3("");
  cerr << "StMuL3Filter::StMuL3Filter(): called. Next BetheBloch instance is made." << endl;
  mBB = new BetheBloch();
  cerr << "StMuL3Filter::StMuL3Filter(): did you see the BetheBloch warning?" << endl;
}

StMuL3Filter::~StMuL3Filter()  {
  delete mBB; mBB = 0;
}

bool StMuL3Filter::accept(const StTrack* track) {
  
  float pCutHigh        = 2.0;    // high momentum cut for RICH/Upsilon candidates 
  int   nHitsCutHighP   = 10;     // nHits cut for all tracks

  // following cuts apply only for tracks with pCutLow < p <pHigh
  float pCutLow             = 0.2;    // low momentum cut
  int   nHitsCutLowP        = 15;    
  int   chargeForLowP       = -1;     // charge for tracks with pCutLow < p < pCutHigh, set to 0 for all tracks
  float dEdxMassCutHigh     = 0.939;  // cut below BetheBloch(p/dEdxMassCutHigh), e.g. proton-band
  float dEdxFractionCutHigh = 0.6;    // cut fraction of dEdx-band, i.e. dEdxFractionCut * BetheBloch(p/dEdxMassCut)
  float dEdxMassCutLow      = 0.494;  // cut above BetheBloch(p/dEdxMassCutLow), e.g. kaon-band
  float dEdxFractionCutLow  = 1.1;

  int iret = 0;

  // next: take all tracks above pCutHigh
  if (track->geometry()->momentum().magnitude() > pCutHigh
      && track->detectorInfo()->numberOfPoints() >= nHitsCutHighP)
        iret = 1;

  // otherwise: take only neg. tracks in a certain dEdx-range
  else {
        if (track->detectorInfo()->numberOfPoints() >= nHitsCutLowP
            && track->geometry()->momentum().magnitude() > pCutLow) {

              int chargeOK = 0;
              int dedxOK = 0;

              // check charge
              if (chargeForLowP==0) 
                    chargeOK = 1;
              else if (track->geometry()->charge() == chargeForLowP) 
                    chargeOK = 1;

              // check dEdx
              //              if (mBB==0) mBB = new BetheBloch();
              float p = track->geometry()->momentum().magnitude();
              float dedxHigh = dEdxFractionCutHigh * mBB->Sirrf(p/dEdxMassCutHigh);
              float dedxLow = dEdxFractionCutLow * mBB->Sirrf(p/dEdxMassCutLow);

              float dedx = 0;
              // get track dEdx
              const StSPtrVecTrackPidTraits& traits = track->pidTraits();
              StDedxPidTraits* dedxPidTr;
              for (unsigned int itrait = 0; itrait < traits.size(); itrait++){
                    dedxPidTr = 0;
                    if (traits[itrait]->detector() == kTpcId) {
                          StTrackPidTraits* thisTrait = traits[itrait];
                          dedxPidTr = dynamic_cast<StDedxPidTraits*>(thisTrait);
                          if (dedxPidTr && dedxPidTr->method() == kTruncatedMeanId) {
                                // adjust L3 dE/dx by a factor of 2 to match offline
                                dedx = 2 * dedxPidTr->mean();
                          }
                    }
              }
              if (dedx > dedxHigh && dedx > dedxLow) 
                    dedxOK = 1;

              // final answer
              iret = chargeOK * dedxOK;
        } // if (pCutLow && nHitsCutLowP)

  }

  return (bool)iret;
}


/***************************************************************************
 *
 * $Log: StMuL3Filter.cxx,v $
 * Revision 1.5  2003/10/20 19:50:13  perev
 * workaround added for TClonesArray::Delete + some cleanup of MuEmc
 *
 * Revision 1.4  2003/01/23 21:59:50  laue
 * Modification to compile on Solaris.
 *
 * Revision 1.3  2002/05/04 23:56:30  laue
 * some documentation added
 *
 * Revision 1.2  2002/03/20 16:04:12  laue
 * minor changes, mostly added access functions
 *
 * Revision 1.1  2002/03/08 17:04:18  laue
 * initial revision
 *
 *
 **************************************************************************/