Cuts.cxx

Go to the documentation of this file.

#include "CAFAna/Cuts/Cuts.h"

#include "StandardRecord/Proxy/SRProxy.h"

#include <bitset>
#include <iostream>

namespace ana
{
  //----------------------------------------------------------------------
  /// Helper for \ref kPassesPresel
  inline bool PassesPreselFunc(const caf::SRProxy* sr)
  {
    // TODO HACK this should be cutting those out anyway
    if(sr->sel.lem.pid < -4) return false;
    return (sr->slc.nhit >= 43 && sr->slc.nhit <= 120 &&
            (sr->trk.discrete.ntracks == 0 || sr->trk.discrete.tracks[0].len < 267) &&
            2*sr->slc.calE > 1 && 2*sr->slc.calE < 3);
  }

  const Cut kPassesPresel(PassesPreselFunc);

  //----------------------------------------------------------------------
  const Cut kCrudeCCSel([](const caf::SRProxy* sr)
                        {
                          if(sr->trk.kalman.ntracks == 0) return false;
                          return sr->trk.kalman.tracks[0].len > 200;
                        });

  //----------------------------------------------------------------------
  const Cut kUCCosmicRej([](const caf::SRProxy* sr)
                         {
                           return (sr->sel.cosrej.anglekal > 0.5 &&
                                   sr->sel.cosrej.numuunconttunedpid > 0.55);
                         });

  //----------------------------------------------------------------------
  const Cut kNDContainNue(
    [](const caf::SRProxy* sr)
    {
      if( sr->vtx.elastic.fuzzyk.nshwlid <=0 )
        return false;

      const caf::SRShowerLIDProxy& shwlid = sr->vtx.elastic.fuzzyk.png[0].shwlid;
      const caf::SRVector3DProxy& start = shwlid.start;
      const caf::SRVector3DProxy& stop = shwlid.stop;

      return (std::min(start.X(), stop.X()) > -150 &&
              std::max(start.X(), stop.X()) < +150 &&
              std::min(start.Y(), stop.Y()) > -150 &&
              std::max(start.Y(), stop.Y()) < +150 &&
              std::min(start.Z(), stop.Z()) >   50 &&
              std::max(start.Z(), stop.Z()) < 1230);
    });

  //----------------------------------------------------------------------
  const Cut kFDContainNue(
    [](const caf::SRProxy* sr)
    {
      if( sr->vtx.elastic.fuzzyk.nshwlid <=0 )
        return false;

      const caf::SRShowerLIDProxy& shwlid = sr->vtx.elastic.fuzzyk.png[0].shwlid;
      const caf::SRVector3DProxy& start = shwlid.start;
      const caf::SRVector3DProxy& stop = shwlid.stop;

      return (std::min(start.X(), stop.X()) > -725 &&
              std::max(start.X(), stop.X()) < +725 &&
              std::min(start.Y(), stop.Y()) > -725 &&
              std::max(start.Y(), stop.Y()) < +725 &&
              std::min(start.Z(), stop.Z()) >   50 &&
              std::max(start.Z(), stop.Z()) < 5910);
    });

  //----------------------------------------------------------------------
  const Cut kInTime([](const caf::SRProxy* sr)
                    {
                      return (sr->slc.meantime >= 208000 &&
                              sr->slc.meantime <= 236000);
                    });

  //----------------------------------------------------------------------
  const Cut kAngVar([](const caf::SRProxy* sr)
                    {
                      std::cout << "\n\n\n"
                                << "WARNING:   This cut used sandbox variables which no longer exist. Aborting."
                                << std::endl;
                      abort();
                      /*
                      if(((sr->trk.cosmic.tracks[0].nhit*1.0)/(1.0*sr->slc.nhit))>0.88){
                        if(sr->slc.nhit>50 && sr->slc.nhit<=135 &&
                           (((sr->sand.nue.anglecosmic*sr->sand.nue.anglecosmic)
                             - (sr->trk.cosmic.tracks[0].dir.Y()*sr->trk.cosmic.tracks[0].dir.Y()))
                            <((1.1/85)*sr->slc.nhit-1.65)))
                          return false;

                        if(sr->slc.nhit>135 && sr->slc.nhit<=250 &&
                           (((sr->sand.nue.anglecosmic*sr->sand.nue.anglecosmic)
                             - (sr->trk.cosmic.tracks[0].dir.Y()*sr->trk.cosmic.tracks[0].dir.Y()))
                            <((.6/115)*sr->slc.nhit-0.6)))
                          return false;

                        if(sr->slc.nhit>250 && sr->slc.nhit<=500 &&
                           (((sr->sand.nue.anglecosmic*sr->sand.nue.anglecosmic)
                             - (sr->trk.cosmic.tracks[0].dir.Y()*sr->trk.cosmic.tracks[0].dir.Y()))
                            < ((0.1/250)*sr->slc.nhit+0.6)))
                          return false;

                        if(sr->slc.nhit>500 &&
                           (((sr->sand.nue.anglecosmic*sr->sand.nue.anglecosmic)
                             - (sr->trk.cosmic.tracks[0].dir.Y()*sr->trk.cosmic.tracks[0].dir.Y())) > 0.8))
                          return false;
                      }
                      */
                      return false;
                    });

  //----------------------------------------------------------------------
  const Cut kXYContainment([](const caf::SRProxy* sr)
                           {
                             return (sr->slc.boxmin.X() > -700 &&
                                     sr->slc.boxmin.Y() > -700 &&
                                     sr->slc.boxmax.X() < 700 &&
                                     sr->slc.boxmax.Y() < 660);
                           });

  //----------------------------------------------------------------------
  const Cut kContained([](const caf::SRProxy* sr)
                       {
                         return (sr->slc.boxmin.Z() > 100 &&
                                 sr->slc.boxmax.Z() < 800 );
                       });

  //----------------------------------------------------------------------
  const Cut kThroughGoing([](const caf::SRProxy* sr)
                          {
                            return (sr->slc.boxmin.Z() < 100 &&
                                    sr->slc.boxmax.Z() > 800 );
                          });

  //----------------------------------------------------------------------
  const Cut kEnter([](const caf::SRProxy* sr)
                   {
                     return (sr->slc.boxmin.Z() < 100 &&
                             sr->slc.boxmax.Z() < 800 );
                   });

  //----------------------------------------------------------------------
  const Cut kExit([](const caf::SRProxy* sr)
                  {
                    return (sr->slc.boxmin.Z() > 100 &&
                            sr->slc.boxmax.Z() > 800 );
                  });

  //----------------------------------------------------------------------
  const Cut kCleanup([](const caf::SRProxy* sr)
                     {
                       return (sr->slc.nhit > 20 &&
                               sr->trk.cosmic.tracks[0].view == 2);
                     });

  //----------------------------------------------------------------------
  bool ContainFxn::operator()(const caf::SRProxy* sr) const
  {
    return (sr->slc.boxmin.X() > -745 &&
            sr->slc.boxmin.Y() > -745 &&
            sr->slc.boxmax.X() < 745 &&
            sr->slc.boxmax.Y() < 720 &&
            sr->slc.boxmin.Z() > zmax+12 &&
            sr->slc.boxmax.Z() < 5950 &&
            sr->sel.contain.kalfwdcell > 10 &&
            sr->sel.contain.kalbakcell > 10 &&
            sr->sel.contain.cosfwdcell > 10 &&
            sr->sel.contain.cosbakcell > 10);
  }

  //----------------------------------------------------------------------
  bool NotContainFxn::operator()(const caf::SRProxy* sr) const
  {
    return ((sr->slc.boxmin.X() < -745 ||
             sr->slc.boxmin.Y() < -745 ||
             sr->slc.boxmax.X() > 745 ||
             sr->slc.boxmax.Y() > 720 ||
             sr->slc.boxmin.Z() > 5950 ||
             sr->slc.boxmin.Z() < zmax+12 ||
             sr->sel.contain.kalfwdcell < 11 ||
             sr->sel.contain.kalbakcell < 11 ||
             sr->sel.contain.cosfwdcell < 11 ||
             sr->sel.contain.cosbakcell < 11) && sr->mc.nu[0].vtx.Z() > zmax);
  }

  //----------------------------------------------------------------------
  const Cut kNueNDCont([](const caf::SRProxy* sr)
                       {
                         bool cut = false;

                         //must have a vertex and nonzero numbers of 3d prongs and showers
                         if (!sr->vtx.elastic.IsValid) return cut;
                         if (sr->vtx.elastic.fuzzyk.npng == 0) return cut;
                         if (sr->vtx.elastic.fuzzyk.nshwlid ==0) return cut;

                         //minimum hit requirement
                         if (sr->slc.nhit < 20) return cut;

                         //longest shower containment cut
                         if (std::min(sr->vtx.elastic.fuzzyk.png[0].shwlid.start.X(),sr->vtx.elastic.fuzzyk.png[0].shwlid.stop.X()) < -180.0) return cut;
                         if (std::max(sr->vtx.elastic.fuzzyk.png[0].shwlid.start.X(),sr->vtx.elastic.fuzzyk.png[0].shwlid.stop.X()) > 180.0) return cut;
                         if (std::min(sr->vtx.elastic.fuzzyk.png[0].shwlid.start.Y(),sr->vtx.elastic.fuzzyk.png[0].shwlid.stop.Y()) < -180.0) return cut;
                         if (std::max(sr->vtx.elastic.fuzzyk.png[0].shwlid.start.Y(),sr->vtx.elastic.fuzzyk.png[0].shwlid.stop.Y()) > 180.0) return cut;
                         if (std::min(sr->vtx.elastic.fuzzyk.png[0].shwlid.start.Z(),sr->vtx.elastic.fuzzyk.png[0].shwlid.stop.Z()) < 25.0) return cut;
                         if (std::max(sr->vtx.elastic.fuzzyk.png[0].shwlid.start.Z(),sr->vtx.elastic.fuzzyk.png[0].shwlid.stop.Z()) > 1225.0) return cut;

                         //vertex containment
                         if (sr->vtx.elastic.vtx.X() < -140.0) return cut;
                         if (sr->vtx.elastic.vtx.X() >  140.0) return cut;
                         if (sr->vtx.elastic.vtx.Y() < -140.0) return cut;
                         if (sr->vtx.elastic.vtx.Y() >  140.0) return cut;
                         if (sr->vtx.elastic.vtx.Z() <  100.0) return cut;
                         if (sr->vtx.elastic.vtx.Z() >  550.0) return cut;

                         //cut on longest fuzzyk 3d track length - do we really mean longest here or do we want highest E?
                         if (sr->vtx.elastic.fuzzyk.png[sr->vtx.elastic.fuzzyk.longestidx].len > 600.0) return cut;

                         //cut on maximum slice energy, remove high energy tail
                         if (sr->slc.calE > 5.0) return cut;

                         //cut on gap between longest shower and vertex, remove a class of reco failure
                         if (sr->vtx.elastic.fuzzyk.png[sr->vtx.elastic.fuzzyk.longestidx].shwlid.gap > 150.0) return cut;

                         cut = true;
                         return cut;
                       });

  //----------------------------------------------------------------------
  const Cut kNueIsMuon([](const caf::SRProxy* sr)
                       {
                         if (sr->vtx.elastic.fuzzyk.nshwlid > 0){
                           for (unsigned int i=0; i<sr->vtx.elastic.fuzzyk.nshwlid; ++i){
                             if (sr->vtx.elastic.fuzzyk.png[i].shwlid.lid.ismuon == 1){
                               return false;
                             }
                           }
                           return true;
                         }
                         return false;
                       });

  //----------------------------------------------------------------------
  const Cut kTruthContainedND(
    [](const caf::SRProxy* sr){
      if (sr->mc.nnu == 0)
        return false;
      if (fabs(sr->mc.nu[0].vtx.X()) > 200 ||
          fabs(sr->mc.nu[0].vtx.Y()) > 200 ||
          sr->mc.nu[0].vtx.Z() < 0 ||
          sr->mc.nu[0].vtx.Z() > 1650)
        return false;
      return true;
    });

  //----------------------------------------------------------------------
  const Cut kIsSingleGap([](const caf::SRProxy* sr){
      // Modeled after kRemoveSmallMasks
      // Selects out only spills that have an upstream diblock that is masked
      // off and have all the active diblocks contiguous
      if (sr->hdr.det != caf::kFARDET) return false;
      std::bitset<14> binary(sr->hdr.dibmask);
      // don't look at configs whose first db is active
      if (binary[0] == 1) return false;
      bool gap = false;
      for (int i=0; i<13; ++i){
        // Is this diblock masked, but the next isn't?
        if (binary[i] == 0 && binary[i+1] == 1){
          if (gap) return false;
          gap = true;
        }
      }
      return gap;
    });

  //----------------------------------------------------------------------
  const Cut kIsNeutrinoInUpstreamDeadDBSingleGap([](const caf::SRProxy* sr){
      // Check that we have a neutrino
      if (sr->mc.nnu == 0) return false;

      // To find the upstream start, we need to where the first active db is
      if (sr->hdr.det != caf::kFARDET) return false;
      std::bitset<14> binary(sr->hdr.dibmask);
      // don't look at configs whose first db is active
      if (binary[0] == 1) return false;
      int FirstActiveDB = 0;
      for (int i=0; i<14; ++i){
        if (binary[i] == 0 && binary[i+1] == 1){
          FirstActiveDB = i+1;
        }
      }

      // Calculate the upstream edge of the first active db
      // Got numbers from gdml, checked 14 db's gives you right det length
      int nspace = 1 + 2*FirstActiveDB;
      double start = nspace*0.45 + 32*(6.61+0.0379)*2*FirstActiveDB; // in cm

      bool isUpStreamOfStart = sr->mc.nu[0].vtx.z < start;
      return isUpStreamOfStart;
    });

  //----------------------------------------------------------------------
  const Cut kIsFarDet ([](const caf::SRProxy* sr)
                       {
                         return sr->hdr.det==caf::kFARDET;
                       });
}