NCAnalysisCutsNC.cxx

Go to the documentation of this file.

//$Id: NCAnalysisCutsNC.cxx,v 1.13 2009/11/23 16:54:04 rodriges Exp $
//
//NCAnalysisCutsNC.cxx
//
//utility holding analysis cut algorithms
//
//B. Rebel 10/2005

#include "NCUtils/Cuts/NCAnalysisCutsNC.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: NCAnalysisCutsNC.cxx,v 1.13 2009/11/23 16:54:04 rodriges Exp $");
// modified on 10 Nov,2009 by J. de Jong for FD cleaning
//----------------------------------------------------------------------
NCAnalysisCutsNC::NCAnalysisCutsNC() :
  NCAnalysisCuts()
{
}

//----------------------------------------------------------------------
NCAnalysisCutsNC::~NCAnalysisCutsNC()
{
}

//----------------------------------------------------------------------
bool NCAnalysisCutsNC::IsGoodBeamEvent()
{
  // Note: Fiducial cut is done seperatly, but data will not be 
  // clean without it

  // Find out if this is MC or data, and the software versions
  ReleaseType::Release_t softType( NCAnalysisCuts::GetReleaseType() );
  Bool_t isMC( ReleaseType::IsMC(softType) );

  //near detector cleaning is done in the filling of the uDST
  //by a separate call to IsCleanHighMultSnarl
  if ( fEventInfo.header->detector == Detector::kNear ) return true;

  /*
  if ( ReleaseType::IsBirch(softType) ){  
    // R1.18.2 style mainEvt cut done here. 
    // Need to tweak MC to fix the lack of noise problem. Not ideal.
    // (Note: to be consistent we should technically tweak the snarlPE
    // variable as well, to adjust the LI PE cut but the effect on the
    // efficiency there is negligble)
    UInt_t tweakedEvents(0);
    Double_t tweakedPH(0.);
    if ( isMC )  NCAnalysisCuts::RollBirchMCSnarlNoise(tweakedEvents, tweakedPH);
    
    // `overlay' (...after a fashion...)  the event onto the noise.
    tweakedEvents += fEventInfo.header->events;
    tweakedPH += fEventInfo.header->snarlPulseHeight;
    
    // now do mainEvt cut. also remove zero (real) event snarls
    if (fEventInfo.header->events == 0) return false;
    if ( ( tweakedEvents == 2 
           && tweakedPH - fEventInfo.event->pulseHeight > 1e4 )
         || tweakedEvents > 2 ) return false;  
    // [now wash your hands]
    
  }  // if (birch/R1.18.2)
  else*/ {
    //Cedar Style mainevt cut  also remove 0 event snarls.
    if (fEventInfo.header->events == 0) return false;
    //check for div by 0. hould'nt happen with >0 evt snarls, but hey. 
    if (fEventInfo.header->snarlPulseHeight==0.)  return false;
    Float_t evtPHfrac(fEventInfo.event->pulseHeight
                      / fEventInfo.header->snarlPulseHeight);
    if ( ( fEventInfo.header->events == 2 && evtPHfrac < 0.75)
         || fEventInfo.header->events > 2 ) return false;
  }  // if ( !birch )
    
  // Normal background cuts
  if(IsCosmicInSpillOx()) return false;

  if(IsFibreNoiseInSpillOx()) return false;

  if(IsLIInSpillOx()) return false;
  
  // Timing cut. Data only
  if(!isMC){
    double dt_spill = fEventInfo.event->stripTime1st * Munits::s;
    dt_spill -= fEventInfo.beam->nearestNSToSpill * Munits::ns;
    
    if ( dt_spill < -2. * Munits::microsecond 
         || dt_spill > 12. * Munits::microsecond ){
      return false;
    } 
  }

  // Survived all cuts
  return true;
}

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

  bool endZOK = false;
  bool xyOK = false;

  if(fEventInfo.header->detector==(int)Detector::kFar){
    if((fEventInfo.track->endZ > 0.5
        && fEventInfo.track->endZ < 12.)
       || (fEventInfo.track->endZ > 16.5
           && fEventInfo.track->endZ < 28.))
      endZOK = true;

    if(TMath::Sqrt(fEventInfo.track->endX*fEventInfo.track->endX
                   +fEventInfo.track->endY*fEventInfo.track->endY) < 3.5)
      xyOK = true;

    if(TMath::Abs(fEventInfo.track->traceEndZ) > 0.5 && xyOK && endZOK)
      stoppingMuon = true;
  }
  else if(fEventInfo.header->detector==(int)Detector::kNear){

    //see position paper on ND flux by Kordosky and Dorman in DocDB for
    //definition of the fiducial volume

    double u = TMath::Sqrt(0.5)*(fEventInfo.track->endX+fEventInfo.track->endY);
    double v = TMath::Sqrt(0.5)*(fEventInfo.track->endY-fEventInfo.track->endX);

    if(0.3 < u && u < 1.8
       && -1.8 < v && v < -0.3
       && fEventInfo.track->endX < 2.4
       && TMath::Sqrt(u*u + v*v) > 0.8
       && fEventInfo.track->endZ < 15. ) stoppingMuon = true;

  }

  return stoppingMuon;
}

//---------------------------------------------------------------------
bool NCAnalysisCutsNC::InBeamFiducialVolume()
{
  // only in FD FidVol if at FD
  if (fEventInfo.header->detector  == (int)Detector::kFar){

    // Get distance to edge, distance to center & vertex z.
    // Use track if available and longer than the shower
    // RPL 16/11/06: Adjust track vertices to be consistent with event
    // vertices, see MINOS-doc-2409

    Float_t dist_edge(fEventInfo.event->vtxMetersToCloseEdge);
    Float_t dist_coil(fEventInfo.event->vtxMetersToCoil);
    Float_t vtxZ(fEventInfo.event->vtxZ);
    if ( (fEventInfo.event->tracks > 0)
         && ( (fEventInfo.track->planes - fEventInfo.shower->planes)>0)) {
      dist_edge = fEventInfo.track->vtxMetersToCloseEdge;
      dist_coil = fEventInfo.track->vtxMetersToCoil;
      vtxZ = fEventInfo.track->vtxZ - NCType::kTrackVtxAdjustment;
    }

    // RPL 16/11/06 Adjust cut values by -3.92 cm to be constent with
    // vertex definition.

    // Longitudinal cuts: > 4 Scint planes from SM fronts...
    //  \-> vZ > 0.2328m (SM1),  vZ > 16.1408m (SM2)
    // ...and >1m from SM ends...
    //  \-> vZ < 13.6208m (SM1),  vZ < 28.8608m (SM2)
    // ...adjusted (tightened) to be at center of an air gap
    if ( vtxZ < 0.2328 * Munits::m || vtxZ > (28.8608+0.119) * Munits::m ) return false;
    if ( vtxZ > (13.6208+0.119) * Munits::m && vtxZ < 16.1408 * Munits::m ) return false;

    // Transverse cuts: 50cm from edge => 7m (edge2edge) octagon
//    if ( dist_edge < 0.5 * Munits::m ) return false;
    if ( dist_edge < 0.4 * Munits::m ) return false;

    // Coil hole cut: > 0.45m radius around center
//    if ( dist_coil < 0.45 * Munits::m ) return false;
    if ( dist_coil < 0.6 * Munits::m ) return false;
  } //end if far

  //only in ND FidVol if at ND
  else if(fEventInfo.header->detector  == (int)Detector::kNear){

    //Get distance to edge & vertex z.
    // Use track if available and longer than the shower
    // RPL 08/03/07: Adjust track vertices to be consistent with event
    // vertices, as in FD. see MINOS-doc-2409
    Float_t dist_edge(fEventInfo.event->vtxMetersToCloseEdge);
    Float_t vtxZ(fEventInfo.event->vtxZ);
    if ( (fEventInfo.event->tracks > 0)
         && ( (fEventInfo.track->planes - fEventInfo.shower->planes)>0) ) {
      dist_edge = fEventInfo.track->vtxMetersToCloseEdge;
      vtxZ = fEventInfo.track->vtxZ - NCType::kTrackVtxAdjustment;
    }

    // Longitudinal cuts: slightly smaller than golden region,
    // to reduce # of showers getting into the spectrometer
    // Start  1/2 way between reco planes 30 & 31
    // Finish 1/2 way between reco planes 80 & 81
    if ( vtxZ > 4.7368 * Munits::m || vtxZ < 1.7 * Munits::m ) return false;

    // Transverse cuts:
    // Use a Detector-oriented set of cuts, since beam oriented cuts
    // are only apropriate for much smaller fiducial regions
    if ( dist_edge < 0.5 * Munits::m ) return false;

    //No coil hole cut since transvers cuts are based on partial plane
  }// end if near

  // All cuts were passed.
  return true;
}

//---------------------------------------------------------------------
bool NCAnalysisCutsNC::InFiducialVolumeTrue()
{
  Detector::Detector_t det = Detector::kFar;
  PlaneCoverage::PlaneCoverage_t outline = PlaneCoverage::kComplete;
  if(fEventInfo.header->detector == (int)Detector::kNear){
    det = Detector::kNear;
    outline = PlaneCoverage::kNearPartial;
  }

  PlaneOutline po;
  float distU = 0.;
  float distV = 0.;
  float xedge = 0.;
  float yedge = 0.;
  float x = fEventInfo.truth->nuVtxX;
  float y = fEventInfo.truth->nuVtxY;
  float z = fEventInfo.truth->nuVtxZ;

  po.DistanceToOuterEdge(x, y, PlaneView::kU, 
                         outline, distU, xedge, yedge);
  po.DistanceToOuterEdge(x, y, PlaneView::kV, 
                         outline, distV, xedge, yedge);

  bool insideU = false;

  if ( outline != PlaneCoverage::kNearPartial ) {
    // circle of radius 1/sqrt(2) will enclose the coilHole, 
    // but is itself enclosed by the outer edges
    insideU = ( (x*x + y*y) < 0.5 * Munits::m2 );
    //cout << "X: " << x << ";   Y: " 
    //     << y << ";   inside:" << insideU << endl;
  }
    
  bool insideV = insideU;
  if ( !insideU ){ 
    insideU = po.IsInside(x, y, PlaneView::kU, outline);
    insideV = po.IsInside(x, y, PlaneView::kV, outline);
  }
  distU *= ( insideU ? 1. : -1. );
  distV *= ( insideV ? 1. : -1. );

  float vtxMetersToCloseEdge = 0.5*(distU + distV);
  float vtxMetersToCoil = TMath::Sqrt(x*x + y*y);

  //based on oxford fiducial volume.
  // only in FD FidVol if at FD  
  if(det == Detector::kFar){
    
    // Longitudinal cuts: > 4 Scint planes from SM fronts...
    //  \-> vZ > 0.2328m (SM1),  vZ > 16.1408m (SM2)
    // ...and >1m from SM ends...
    //  \-> vZ < 13.6208m (SM1),  vZ < 28.8608m (SM2)
    // ...adjusted (tightened) to be at center of an air gap
    // decrease end of SM by ~2 planes
    if ( z < 0.2328*Munits::m || z > (28.8608+0.119) * Munits::m ) return false;
    if ( z > (13.6208+0.119)*Munits::m && z < 16.1408 * Munits::m ) return false;
    
    // Transverse cuts: 50cm from edge => 7m (edge2edge) octagon 
//    if ( vtxMetersToCloseEdge < 0.5 * Munits::m ) return false;
    if ( vtxMetersToCloseEdge < 0.4 * Munits::m ) return false;
    
    // Coil hole cut: > 0.45m radius around center
//    if ( vtxMetersToCoil < 0.45 * Munits::m ) return false; 
    if ( vtxMetersToCoil < 0.6 * Munits::m ) return false; 
  } //end if far
  
  //only in ND FidVol if at ND  
  else if(det == Detector::kNear){
      
    // Longitudinal cuts: slightly smaller than golden region,
    // to reduce # of showers getting into the spectrometer
    // Start  1/2 way between reco planes 30 & 31
    // Finish 1/2 way between reco planes 80 & 81
    if ( z > 4.7368 * Munits::m || z < 1.728 * Munits::m ) return false;
    
    // Transverse cuts:
    // Use a Detector-oriented set of cuts, since beam oriented cuts 
    // are only apropriate for much smaller fiducial regions
    if ( vtxMetersToCloseEdge < 0.5 * Munits::m ) return false;
    
    //No coil hole cut since transvers cuts are based on partial plane    
  }// end if near
  
  return true;
}

//----------------------------------------------------------------------
bool NCAnalysisCutsNC::PassesFinalSelection()
{
  if(!InBeamFiducialVolume()) return false;
  return true;
}

//----------------------------------------------------------------------
bool NCAnalysisCutsNC::PassesFinalSelection(ANtpRecoInfo *recoInfo)
{
  if(recoInfo->inFiducialVolume < 1)  return false;
  return true;
}

//----------------------------------------------------------------------
bool NCAnalysisCutsNC::IsFibreNoiseInSpillOx()
{
  ReleaseType::Release_t softType(GetReleaseType());

  // New cuts based on dogwood data
  if(ReleaseType::IsCedar(softType))
    {
      
      if(fEventInfo.event->pulseHeight < 2500 &&
         fEventInfo.event->totalStrips <= 8) return true;
      
      //Expand the cut area in chopped data
      //      ReleaseType::Release_t softType(GetReleaseType());
      if(!ReleaseType::IsBirch(softType))
      {
        if(fEventInfo.event->pulseHeight < 5000 &&
           fEventInfo.event->totalStrips <= 4) return true;
      }
    }
  else
    {
      if(fEventInfo.event->pulseHeight < 3750 &&
         fEventInfo.event->totalStrips <= 8) return true;
      if(fEventInfo.event->pulseHeight < 2000 &&
         fEventInfo.event->totalStrips > 8) return true;

    }
  return false;
}


//----------------------------------------------------------------------
bool NCAnalysisCutsNC::IsCosmicInSpillOx()
{
  //Ultra-steep showers
  Float_t strips_f(fEventInfo.event->totalStrips);
  Float_t planes_f(fEventInfo.event->planes);
  if ( planes_f == 0.0 || strips_f/(planes_f*planes_f) >= 1 ) return true;


  //Steep showers
  if ( fEventInfo.event->showers > 0 ){
    Float_t tRMS2 = SQR(fEventInfo.shower->transverseRMSU);
    tRMS2          += SQR(fEventInfo.shower->transverseRMSV);
    Float_t tRMS = TMath::Sqrt( tRMS2 );
    if ( ( 0.3 + 0.1901*TMath::Log10(planes_f) ) < tRMS ) return true;
  } //event.showers > 0


  if ( fEventInfo.event->tracks > 0 ){
    //Track angle
    if ( fEventInfo.track->dcosZVtx < 0.4 ) return true;

    //Stopping muons
    Bool_t upExitTrack(false);
    if ( ( fEventInfo.track->endZ > 28.78
           || fEventInfo.track->endMetersToCloseEdge < 0.5 )
         && fEventInfo.track->endY > -1.657 ) upExitTrack = true;

    if ( upExitTrack
         && fEventInfo.track->dtdz < -0.5
         && -fEventInfo.track->dcosYVtx < -0.4 ) return true;

  } // event.tracks > 0

  // Got here => is not a cosmic.
  return false;
}

//----------------------------------------------------------------------
bool NCAnalysisCutsNC::IsLIInSpillOx()
{
  //RPL 11/07/07:
  // removed snarlPE > 8000 cut as it is not really necessary.
  ReleaseType::Release_t softType(GetReleaseType());

  // New cuts based on dogwood data
  if(ReleaseType::IsCedar(softType))
   {
      return fEventInfo.header->isLIold || fEventInfo.header->triggerPMTTime > 0;
    }
   else // new cuts based on dogwood
    {
       return fEventInfo.header->isLI || fEventInfo.header->triggerPMTTime > 0;
    }
}

---