NCAnalysisCutsCC.cxx

Go to the documentation of this file.

//$Id: NCAnalysisCutsCC.cxx,v 1.5 2009/03/04 12:07:46 bckhouse Exp $
//
//NCAnalysisCuts.cxx
//
//utility holding analysis cut algorithms
//
//B. Rebel 10/2005

#include "NCUtils/Cuts/NCAnalysisCutsCC.h"
#include "NCUtils/NCType.h"
#include "NCUtils/NCUtility.h"
#include "AnalysisNtuples/ANtpEventInfoNC.h"
#include "AnalysisNtuples/ANtpShowerInfoNC.h"
#include "AnalysisNtuples/ANtpBeamInfo.h"
#include "AnalysisNtuples/ANtpHeaderInfo.h"
#include "AnalysisNtuples/ANtpTrackInfoNC.h"
#include "AnalysisNtuples/ANtpTruthInfoBeam.h"
#include "AnalysisNtuples/ANtpRecoInfo.h"
#include "MessageService/MsgService.h"
#include "Conventions/Detector.h"
#include "Conventions/Munits.h"
#include "Conventions/SimFlag.h"
#include "Conventions/PlaneView.h"
#include "Conventions/PlaneCoverage.h"
#include "DataUtil/PlaneOutline.h"

#include "TMath.h"

using namespace BeamType;
using std::vector;
using NC::Utility::SQR;

CVSID("$Id: NCAnalysisCutsCC.cxx,v 1.5 2009/03/04 12:07:46 bckhouse Exp $");

//----------------------------------------------------------------------
NCAnalysisCutsCC::NCAnalysisCutsCC() :
  NCAnalysisCuts()
{
  return;
}

//----------------------------------------------------------------------
NCAnalysisCutsCC::~NCAnalysisCutsCC()
{
  return;
}

//----------------------------------------------------------------------
bool NCAnalysisCutsCC::IsGoodBeamEvent()
{
  MSG("NCAnalysisCutsCC", Msg::kDebug) << "IsGoodBeamEventCC" << endl;

  bool goodEvent = true;
  bool recovered = false;

  //recover events that fail the passedFit cut
  if(fEventInfo.header->detector==(int)Detector::kNear
     && fEventInfo.track->passedFit < 1
     && TMath::Abs(fEventInfo.track->begPlaneU - fEventInfo.track->begPlaneV) < 6
     && TMath::Abs(fEventInfo.track->endPlaneU - fEventInfo.track->endPlaneV) < 41
     && fEventInfo.track->endPlane < 270
     )
    recovered = true;

  if(fEventInfo.event->tracks < 1
     || (fEventInfo.track->passedFit < 1 && !recovered)
     ) 
    goodEvent = false;

  if(goodEvent && fEventInfo.header->detector==(int)Detector::kFar){

    if(!IsFarDetectorBeamEvent() || !TrackInBeamDirection())
      goodEvent = false;

  }

  return goodEvent;
}

//----------------------------------------------------------------------
bool NCAnalysisCutsCC::IsFarDetectorBeamEvent()
{
  bool goodEvent = true;

  if(fEventInfo.header->dataType == (int)SimFlag::kMC)
    fEventInfo.header->triggerPMTTime = -1.;
  
  bool liEvent = false;
  bool onCrateBoundary = false;
  double eventTime = 0.;
  
  //is the event starting or stopping on a crate boundary?
  if( (fEventInfo.event->begPlane < 3 && fEventInfo.event->endPlane > 62
       && fEventInfo.event->endPlane < 65)
      || (fEventInfo.event->begPlane > 64 && fEventInfo.event->begPlane < 67
          && fEventInfo.event->endPlane > 126 && fEventInfo.event->endPlane < 129)
      || (fEventInfo.event->begPlane > 128 && fEventInfo.event->begPlane < 131
          && fEventInfo.event->endPlane > 190 && fEventInfo.event->endPlane < 193)
      || (fEventInfo.event->begPlane > 192 && fEventInfo.event->begPlane < 195
          && fEventInfo.event->endPlane > 246 && fEventInfo.event->endPlane < 249)
      || (fEventInfo.event->begPlane > 248 && fEventInfo.event->begPlane < 252
          && fEventInfo.event->endPlane > 311 && fEventInfo.event->endPlane < 314)
      || (fEventInfo.event->begPlane > 313 && fEventInfo.event->begPlane < 316
          && fEventInfo.event->endPlane > 375 && fEventInfo.event->endPlane < 378)
      || (fEventInfo.event->begPlane > 377 && fEventInfo.event->begPlane < 380
          && fEventInfo.event->endPlane > 439 && fEventInfo.event->endPlane < 442)
      || (fEventInfo.event->begPlane > 441 && fEventInfo.event->begPlane < 444
          && fEventInfo.event->endPlane > 483)
      ) onCrateBoundary = true;
  
  //*********Should change this to use the isLI in the ANtpHeaderInfo object

  if(fEventInfo.header->triggerPMTTime > 0.
     || (fEventInfo.event->showers > 0 && fEventInfo.event->pulseHeight > 1.e6)
     || onCrateBoundary
     )
    liEvent = true;
    
  //find the difference between the spill server prediction and the
  //time of the event.
  
  //time of the event is defined to be NtpEvenSummarySR->trigtime
  //and is earliest recorded hit in the event
  eventTime = fEventInfo.event->stripTime1st;
  if(fEventInfo.header->dataType == (int)SimFlag::kMC){
    eventTime = 0.;
    fEventInfo.beam->nearestNSToSpill = 0;
  }
  
  //beamInfo->nearestNSToSpill is the difference in ns between
  //the spill time and the start of the snarl
  if(1.e-9*fEventInfo.beam->nearestNSToSpill - eventTime < -20.e-6
     || 1.e-9*fEventInfo.beam->nearestNSToSpill - eventTime > 30.e-6
     || liEvent
     || fEventInfo.event->largestEventInSnarl < 1
     )
    goodEvent = false;
  
  MAXMSG("NCAnalysisCutsCC", Msg::kDebug, 20) 
    << fEventInfo.event->stripTime1st
    << " " << fEventInfo.header->triggerTime
    << " " << fEventInfo.beam->timeToNearestSpill
    << " " << eventTime
    << " " << liEvent
    << " " << fEventInfo.event->largestEventInSnarl
    << " " << goodEvent
    << endl;


  return goodEvent;
}

//----------------------------------------------------------------------
bool NCAnalysisCutsCC::TrackInBeamDirection()
{
  bool goodDirection = true;

  //setup NEAR angle - currently not used for near detector
  double beamDCosY = TMath::Cos(0.0575958653+TMath::Pi()*0.5);
  double beamDCosZ = TMath::Sqrt(1. - beamDCosY*beamDCosY);

  //setup FAR angle
  // use 3 degrees upward in far det since that is what David is using
  // the real number is more like 3.27 degrees
  beamDCosY = TMath::Cos(-0.052359877+TMath::Pi()*0.5);
  beamDCosZ = TMath::Sqrt(1. - beamDCosY*beamDCosY);
  
  if(fEventInfo.track->dcosYVtx*beamDCosY + fEventInfo.track->dcosZVtx*beamDCosZ < 0.6)
    goodDirection = false;

  return goodDirection;
}

//----------------------------------------------------------------------
bool NCAnalysisCutsCC::IsStoppingBeamMuon()
{
  bool stoppingMuon = false;

  if(fEventInfo.header->detector == int(Detector::kNear)
     && fEventInfo.track->endZ < 16.
     && TMath::Sqrt(SQR(fEventInfo.track->endX)
                    +SQR(fEventInfo.track->endY)) > 0.4
     && fEventInfo.track->endX < 2.7 && fEventInfo.track->endX > -1.65
     && fEventInfo.track->endY < 1.65 && fEventInfo.track->endY > -1.65
     && fEventInfo.track->endY > -1.*(fEventInfo.track->endX + 1.65)
     && fEventInfo.track->endY < fEventInfo.track->endX + 1.65
     && fEventInfo.track->endY < 3.55 - fEventInfo.track->endX
     && fEventInfo.track->endY > fEventInfo.track->endX - 3.55
     )
    stoppingMuon = true;
  else if(fEventInfo.header->detector == (int)Detector::kFar
          && fEventInfo.track->endZ > 0.5 && fEventInfo.track->endZ < 29.5
          && fEventInfo.track->endMetersToCloseEdge > 0.5)
    stoppingMuon = true;

  return stoppingMuon;
}

//---------------------------------------------------------------------
bool NCAnalysisCutsCC::InBeamFiducialVolume()
{
  bool inFiducialVolume = true;

  double vtxX = 0.;
  double vtxY = 0.;
  double vtxZ = 0.;

 if(fEventInfo.event->tracks > 0){
  vtxX = fEventInfo.track->vtxX;
  vtxY = fEventInfo.track->vtxY;
  vtxZ = fEventInfo.track->vtxZ;
 }
 else{
   vtxX = fEventInfo.event->vtxX;
   vtxY = fEventInfo.event->vtxY;
   vtxZ = fEventInfo.event->vtxZ;
 }

  //cuts made for selecting far detector CC events based on docdb #1466
  if(fEventInfo.header->detector == (int)Detector::kFar){
    if(SQR(vtxX) + SQR(vtxY) > 14.
       || vtxZ < 0.5
       || (vtxZ > 14.3 && vtxZ < 16.2)
       || vtxZ > 28.)
      inFiducialVolume = false;
  }
  else if(fEventInfo.header->detector == (int)Detector::kNear){

    double x = vtxX;
    double y = vtxY;

    double fBeamCenterX = 1.4885;
    double fBeamCenterY = 0.1397;
    x -= fBeamCenterX;
    y -= fBeamCenterY;

    if(vtxZ <= 1.
       || vtxZ >= 5.
       || x*x + y*y >= 1.)
      inFiducialVolume = false;

  }

  MAXMSG("NCAnalysisCutsCC", Msg::kDebug, 20) << vtxX << " " << vtxY
                                           << " " << vtxZ << endl;

  return inFiducialVolume;
}

//----------------------------------------------------------------------
bool NCAnalysisCutsCC::PassesFinalSelection(ANtpRecoInfo *recoInfo)
{
  //throw out events that are not in the fiducial volume
  return recoInfo->inFiducialVolume >= 1 && fEventInfo.track->fitMomentum <= 0;
}

//----------------------------------------------------------------------
bool NCAnalysisCutsCC::PassesFinalSelection()
{
  return true;
}

---