NCAnalysisCuts.cxx

Go to the documentation of this file.

//$Id: NCAnalysisCuts.cxx,v 1.27 2009/12/02 14:21:12 tinti Exp $
//
//NCAnalysisCuts.cxx
//
//utility holding analysis cut algorithms
//
//B. Rebel 10/2005

#include "NCUtils/Cuts/NCAnalysisCuts.h"
#include "NCUtils/NCType.h"
#include "NCUtils/NCUtility.h"
#include "AnalysisNtuples/ANtpDefaultValue.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 "DataUtil/EnergyCorrections.h"

#include "TMath.h"

#include <cassert>

using namespace EnergyCorrections;
using NC::Utility::SQR;

CVSID("$Id: NCAnalysisCuts.cxx,v 1.27 2009/12/02 14:21:12 tinti Exp $");

//----------------------------------------------------------------------
NCAnalysisCuts::NCAnalysisCuts() :
  fBeamType(BeamType::kUnknown),
  fCurrentRun(-1),
  fCurrentSnarl(-1),
  fCurrentSubRun(-1)
{
}

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

//----------------------------------------------------------------------
//use this method to select data from the appropriate beam/horn/intensity
//configurations this method assumes that it is only called for data as
//there are not mixed target positions in the MC files
bool NCAnalysisCuts::IsGoodTarget()
{
  bool goodTarget = false;

  BeamType::BeamType_t bt = (BeamType::BeamType_t)(fEventInfo.beam->beamType);

  bool is185i = (bt == BeamType::kL010z185i
                 || bt == BeamType::kL010z185i_lowintensity);
  bool isL010z = (bt == BeamType::kL010z185i
                  || bt == BeamType::kL010z185i_lowintensity
                  || bt == BeamType::kL010z000i
                  || bt == BeamType::kL010z170i
                  || bt == BeamType::kL010z200i);

  //take any thing if the target position is set to all
  //beams or pME or pHE
  if(fUseAllBeams)
    goodTarget = true;
  else if(fBeamType == BeamType::kL100z200i && bt == fBeamType)
    goodTarget = true;
  else if(fBeamType == BeamType::kL250z200i && bt == fBeamType)
    goodTarget = true;
  else if(fBeamType == BeamType::kL010z185i
          && !fUseAll185i && !fUseAllL010z && bt == fBeamType)
    goodTarget = true;
  else if(fBeamType == BeamType::kL010z170i
          && !fUseAllL010z && bt == fBeamType)
    goodTarget = true;
  else if(fBeamType == BeamType::kL010z200i
          && !fUseAllL010z && bt == fBeamType)
    goodTarget = true;
  else if(fBeamType == BeamType::kL010z185i_lowintensity
          && !fUseAll185i && !fUseAllL010z && bt == fBeamType)
    goodTarget = true;
  else if(fBeamType == BeamType::kL010z000i
          && !fUseAllL010z && bt == fBeamType)
    goodTarget = true;
  else if(fUseAllL010z && isL010z)
    goodTarget = true;
  else if(fUseAll185i && is185i)
    goodTarget = true;


  MAXMSG("NCAnalysisCuts", Msg::kDebug, 10) << fEventInfo.beam->beamType << " "
                                           << BeamType::AsString(bt) << " "
                                           << goodTarget << " "
                                           << (int)fUseAllBeams << " "
                                           << is185i << " " << isL010z << endl;


  return goodTarget;
}

//----------------------------------------------------------------------
bool NCAnalysisCuts::IsNewSnarl()
{
  bool newSnarl = false;

  if(fEventInfo.header->run != fCurrentRun){
    newSnarl       = true;
    fCurrentRun    = fEventInfo.header->run;
    fCurrentSubRun = fEventInfo.header->subRun;
    fCurrentSnarl  = fEventInfo.header->snarl;
  }
  else if(fEventInfo.header->run == fCurrentRun &&
          fEventInfo.header->subRun != fCurrentSubRun){
    newSnarl = true;
    fCurrentRun = fEventInfo.header->run;
    fCurrentSubRun = fEventInfo.header->subRun;
    fCurrentSnarl = fEventInfo.header->snarl;
  }
  else if(fEventInfo.header->run == fCurrentRun &&
          fEventInfo.header->subRun == fCurrentSubRun &&
          fEventInfo.header->snarl != fCurrentSnarl){
    newSnarl = true;
    fCurrentRun = fEventInfo.header->run;
    fCurrentSubRun = fEventInfo.header->subRun;
    fCurrentSnarl = fEventInfo.header->snarl;
  }

  return newSnarl;
}


//----------------------------------------------------------------------
bool NCAnalysisCuts::IsMultiCutsClean()
{
  // 35ns timing cut
  if(TMath::Abs(fEventInfo.event->minTimeSeparation) < 35e-9) return false;

  // tiny events are mostly junk
  if(fEventInfo.event->totalStrips < 5) return false;

  // this cuts very steep showers, leaking in
  if(fEventInfo.event->totalStrips/SQR(fEventInfo.event->planes)>1.15)
    return false;

  // this cuts leakage which leaves activity in partially instrumented
  // region
  if(fEventInfo.event->edgeActivityStrips > 2 &&
     fEventInfo.event->edgeActivityPH > 1000 &&
     fEventInfo.event->energyGeV < 5.0 &&
     fEventInfo.shower->planes > fEventInfo.track->planes) // don't cut out CC
    return false;

  // also cut out events with too much activity in the opposite edge region
  if(fEventInfo.event->oppEdgeStrips > 2 &&
     fEventInfo.event->energyGeV < 5.0 &&
     fEventInfo.shower->planes > fEventInfo.track->planes) // don't cut out CC
    return false;

  // make additional deltaZ cuts if (|minDeltaT|<175ns)
  if(TMath::Abs(fEventInfo.event->closeTimeDeltaZ) < 1.0 &&
     TMath::Abs(fEventInfo.event->minTimeSeparation) < 175e-9)
    return false;

  // Survived all cuts
  return true;
}

//----------------------------------------------------------------------  
bool NCAnalysisCuts::IsSimpleCutsClean()
{
  // 155ns timing cut                                                      
  if(TMath::Abs(fEventInfo.event->minTimeSeparation) < 155e-9) 
    return false;

  //cut on portion of the slice in the event                              
  if(fEventInfo.event->slicePHFraction<0.9) return false;

  // also cut events that share a good fraction of physical strips with
  // other events                                                            
  // and don't cut out CC                                                  
  if(fEventInfo.event->sharedStripFraction>0.3 &&
     ((-fEventInfo.shower->planes + fEventInfo.track->planes)<=5.))
    return false;

// Survived all cuts                                                    
  return true;
}

//----------------------------------------------------------------------
bool NCAnalysisCuts::IsCleanLowMultSnarl()
{
  return fEventInfo.header->events == 1 ||
    (fEventInfo.header->events == 2 &&
     TMath::Abs(fEventInfo.event->minTimeSeparation*1e6) > 0.2);
}

//----------------------------------------------------------------------
void NCAnalysisCuts::SetInfoObjects(ANtpHeaderInfo* headerInfo,
                                    ANtpBeamInfo* beamInfo,
                                    ANtpEventInfoNC* eventInfo,
                                    ANtpTrackInfoNC* trackInfo,
                                    ANtpShowerInfoNC* showerInfo,
                                    ANtpTruthInfoBeam* truthInfo)
{
  fEventInfo.header = headerInfo;
  fEventInfo.beam = beamInfo;
  fEventInfo.event = eventInfo;
  fEventInfo.track = trackInfo;
  fEventInfo.shower = showerInfo;
  fEventInfo.truth = truthInfo;
}

//----------------------------------------------------------------------
void NCAnalysisCuts::SetBeamType(BeamType::BeamType_t type)
{
  fBeamType = type;
}

//----------------------------------------------------------------------
void NCAnalysisCuts::SetUseAllBeams(bool useAll)
{
  fUseAllBeams = useAll;
}

//----------------------------------------------------------------------
void NCAnalysisCuts::SetUseAllL010z(bool useAll)
{
  fUseAllL010z = useAll;
}

//----------------------------------------------------------------------
void NCAnalysisCuts::SetUseAll185i(bool useAll)
{
  fUseAll185i = useAll;
}

//----------------------------------------------------------------------
bool NCAnalysisCuts::IsGoodSpill()
{
  //mc beam is always perfect and only using the
  //correct mc files with each beamType, so no need to
  //check the target
  if(fEventInfo.header->dataType == int(SimFlag::kMC)) return true;

  if(!IsGoodTarget()) return false;

  //check the beam quality cuts, the spill type and
  //error on the gps - error has to be big within 5 minutes of the
  //spill - only check the gps for far detector
  if(fEventInfo.beam->goodSpill < 1 ||
     (fEventInfo.header->detector == int(Detector::kFar) &&
      fEventInfo.beam->deltaSecToSpillGPS < 360 &&
      fEventInfo.beam->gpsError > 1000) ||
     fEventInfo.header->spillType == 3)
    return false;

  return true;
}

//----------------------------------------------------------------------
bool NCAnalysisCuts::IsGoodBeamSnarl()
{
  //MC snarls are always good
  if(fEventInfo.header->dataType == int(SimFlag::kMC)) return true;

  // check that there were no problems with the detector
  // Don't take reversed field data for uDSTs
  if(fEventInfo.header->coilStatus < 1 ||
     (fEventInfo.header->detector == int(Detector::kFar) &&
      fEventInfo.header->isGoodData < 1))
    return false;

  //the cuts in this method come from MINOS-doc-1342
  if(!IsGoodSpill()) return false;

  if(fEventInfo.header->detector == int(Detector::kNear) &&
     (fEventInfo.header->coilCurrent > -1000 ||
      fEventInfo.header->isGoodData < 1))
    return false;

  return true;
}

//----------------------------------------------------------------------
//these cuts were presented in the 2/9/06 reconstruction meeting by
//N. Saoulidou
//this is a little cheesy as the code is similar to what is in NCUtils
bool NCAnalysisCuts::IsGoodShower()
{
  //set the version for the energy corrections.  mike uses R1_18 and cedar, the
  //values in ANtpHeaderInfo are Birch(Cedar) Data(Carrot/Daikon/Etc)
  if( fEventInfo.header->softVersion.Contains("Cedar") ){
    const EnergyCorrections::CorrectionVersion_t ver = EnergyCorrections::kCedar;
    SetCorrectionVersion(ver);
  }
  else if( fEventInfo.header->softVersion.Contains("Birch") ){
    const EnergyCorrections::CorrectionVersion_t ver = EnergyCorrections::kBirch;
    SetCorrectionVersion(ver);
  }


  bool isdata = true;
  if(fEventInfo.header->dataType == (int)SimFlag::kMC) isdata = false;

  Detector::Detector_t det = Detector::kNear;
  if(fEventInfo.header->detector == (int)Detector::kFar)
    det = Detector::kFar;

  double showerEnergy;
  //Protect against events without showers.
  if(fEventInfo.shower->linearCCGeV > 0.){
      showerEnergy = CorrectShowerEnergy(fEventInfo.shower->linearCCGeV,
                                         det,
                                         CandShowerHandle::kCC,
                                         1, isdata);
  }
  else showerEnergy = ANtpDefaultValue::kDouble;

  if(fEventInfo.event->tracks < 1
     || TMath::Abs(fEventInfo.shower->vtxZ - fEventInfo.track->aTrkvtx_ns) < 0.5
     ||  showerEnergy > 2. )
    return true;

  return false;
}

/*
//----------------------------------------------------------------------
void NCAnalysisCuts::RollBirchMCSnarlNoise(UInt_t& extraEvents,
                                             Double_t& extraPH){

  //Seed the generator.  Each snarl is uniquely identified by run,
  // ( 21 001 001 to 21 001 100 ) and snarl ( 0 to 2000ish ) numbers
  // each run is uniquely identifed by index numbers (0 to 5ish)
  // seed is a UInt_t  (0 to 4,294,967,295 )

  UInt_t seed(fEventInfo.header->run % 1000);
  seed *= 2048;
  seed += fEventInfo.header->snarl;
  seed *= 16;
  seed += fEventInfo.event->index;
  fRandom->SetSeed(seed);

  //roll a Poisson to add an event(s) randomly. lambda = 7.741e-4
  const Float_t meanEvt(7.741e-4);

  extraEvents = fRandom->Poisson(meanEvt);

  //Add PH based on how many events we added
  // sqrt(noise) has gaussian distn. Ignore the small probability of
  // going below zero---will square the result anyway.

   Double_t sqrtPH(0);
   if (extraEvents == 0) {
     sqrtPH = fRandom->Gaus(63.33, 8.673);
   }
   else if ( extraEvents == 1 && fRandom->Rndm() > 0.115 ) {
     //11.5% of 1 event fake trigger snarls seem to have something
     // substantial (a cosmic?)
     // Also use this distribution for > 1 events added.
     sqrtPH = fRandom->Gaus(71.02, 7.772);
   }
   else {
     sqrtPH = fRandom->Gaus(259.6, 90.27);
   }

   extraPH = sqrtPH*sqrtPH;

}
*/
//----------------------------------------------------------------------
Int_t NCAnalysisCuts::GetReleaseType()
{
  // Get the ReleaseType from the string.

  //const char* asCharArray= (const char*)fEventInfo.header->softVersion;
  ReleaseType::Release_t softType =
    ReleaseType::StringToType(fEventInfo.header->softVersion);

  assert(softType != -1);

  return softType;
}

---