NueCCIncCuts.h

Go to the documentation of this file.

#pragma once

#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/Utilities.h"
#include "NDAna/nuecc_inc/NueCCIncVars.h"
#include "StandardRecord/Proxy/SRProxy.h"
#include "3FlavorAna/Cuts/NueCutsSecondAna.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Cuts/TruthCuts.h"

#include "NDAna/numucc_inc/NDXSecMuonPID.h"
#include "NDAna/nuecc_inc/NDXSecElectronPID.h"

#include "TMath.h"
#include "TVector3.h"


namespace ana
{
  namespace nueccinc
  {
    // Cut on ND vertex position        
    const Cut kcNueCCIncFiducial(
                                 [](const caf::SRProxy* sr)
                                 {
                                   assert(sr->vtx.nelastic > 0 && 
                                          "Must apply DQ cuts");        
                                   if (sr->vtx.elastic[0].vtx.X() < -130.0) 
                                     return false;
                                   if (sr->vtx.elastic[0].vtx.X() >  150.0) 
                                     return false;
                                   if (sr->vtx.elastic[0].vtx.Y() < -140.0) 
                                     return false;
                                   if (sr->vtx.elastic[0].vtx.Y() >  140.0) 
                                     return false;
                                   if (sr->vtx.elastic[0].vtx.Z() <  150.0) 
                                     return false;
                                   if (sr->vtx.elastic[0].vtx.Z() >  800.0) 
                                     return false;
                                   return true;
                                 });
    
    const Cut 
      kcNueCCIncContainment([](const caf::SRProxy* sr){
          if(sr->vtx.elastic[0].fuzzyk.nshwlid < 1) return false;
          TVector3 start = sr->vtx.elastic[0].fuzzyk.png[0].shwlid.start;
          TVector3 stop  = sr->vtx.elastic[0].fuzzyk.png[0].shwlid.stop;
          
          //Make sure none of the event goes into the muon catcher
          for(uint ix = 0; ix < sr->vtx.elastic[0].fuzzyk.nshwlid; ix++){
            TVector3 start_muoncatcher = 
              sr->vtx.elastic[0].fuzzyk.png[ix].shwlid.start;
            TVector3 stop_muoncatcher  = 
              sr->vtx.elastic[0].fuzzyk.png[ix].shwlid.stop;
            if (std::max(start_muoncatcher.Z(),
                         stop_muoncatcher.Z()) >  1250.0) return false;
          }     
          
          
          if(sr->sel.nuecosrej.distallpngtop < 50) return false;
          if(sr->sel.nuecosrej.distallpngbottom < 30) return false;
          if(sr->sel.nuecosrej.distallpngeast < 50) return false;
          if(sr->sel.nuecosrej.distallpngwest < 30) return false;
          if(sr->sel.nuecosrej.distallpngfront < 150) return false;     
          return true;
        });
    
    
    //Phase Space Cuts
    const NuTruthCut kTrueElectronPhaseSpaceCutST
      ([](const caf::SRNeutrinoProxy* sr){
        float costheta = kTrueElecCosthetaST(sr);
        float elece    = kTrueElecEST(sr);
        
        if( (costheta >= 0.97) && (elece >= 1.4) && (elece < 6.0))
          return true;
        if( (costheta >= 0.94) && (costheta < 0.97) &&
            (elece >= 1.4) && (elece < 4.1))
          return true;
        if( (costheta >= 0.90) && (costheta < 0.94) &&
            (elece >= 1.0) && (elece < 2.5))
          return true;
        if( (costheta >= 0.85) && (costheta < 0.90) &&
            (elece >= 1.0) && (elece < 2.0))
          return true;

        return false;
      });

    const Cut kRecoElectronPhaseSpaceCut
      ([](const caf::SRProxy* sr){
        float costheta = kvRecoCos(sr);
        float elece    = kvElecE(sr);
        
        if( (costheta >= 0.97) && (elece >= 1.4) && (elece < 6.0))
          return true;
        if( (costheta >= 0.94) && (costheta < 0.97) &&
            (elece >= 1.4) && (elece < 4.1))
          return true;
        if( (costheta >= 0.90) && (costheta < 0.94) &&
            (elece >= 1.0) && (elece < 2.5))
          return true;
        if( (costheta >= 0.85) && (costheta < 0.90) &&
            (elece >= 1.0) && (elece < 2.0))
          return true;

        return false;
      });
    
    //Selection Cuts
    //--------------
    const Cut kcDQ           = kNueSecondAnaDQ;
    const Cut kcFiducial     = kNueNDSecondAnaFiducial;
    const Cut kcContainment  = kNueNDSecondAnaContain;
    const Cut kcFrontPlanes  = kNueNDSecondAnaFrontPlanes;
    const Cut kcNHits        = kNueNDSecondAnaNHits;
    const Cut kcProngLength  = kNueNDSecondAnaProngLength;
    const Cut kcRemid        = kRemID < 0.40;//Optimize using uncertainty
                                              // on the cross section
    
    const Cut kcMuonIDCut    = kMuonID < -0.20;//Optimized using uncertainty on 
                                               // the cross section


    const Cut kcCVN          = kvCVN2017 > 0.85; //0.78;  
                                                 //Optimize using uncertainty 
                                                 // on the cross section and
                                                 // statistical FoM 

    const Cut kcElectronIDCut   = kElectronID > 0.2; 
    //Non-optimized but
    // need to test something

    const Cut kcPresel_loose = kcDQ && kcNueCCIncFiducial && 
      kcNueCCIncContainment && kcFrontPlanes && kcNHits;
    
    const Cut kcPresel =  kcPresel_loose && kcMuonIDCut;

    //MC Cuts
    //---------------
    
    const Cut ktMC( [](const caf::SRProxy* sr){
        if(sr->hdr.ismc) return true;
        return false;
      });
    
    const Cut ktSignal( [](const caf::SRProxy* sr){
        if(sr->mc.nnu == 0) return false;
        return (sr->mc.nu[0].iscc && 
                (sr->mc.nu[0].pdg == 12) &&
                (sr->mc.nu[0].pdgorig == 12));
      });
    
    const Cut ktNuebarCC( [](const caf::SRProxy* sr){
        if(sr->mc.nnu == 0) return false;
        return (sr->mc.nu[0].iscc && 
                (sr->mc.nu[0].pdg == -12) &&
                (sr->mc.nu[0].pdgorig == -12));
      });
    
    const Cut ktNumuCC( [](const caf::SRProxy* sr){
        if(sr->mc.nnu == 0) return false;
        return (sr->mc.nu[0].iscc && 
                (sr->mc.nu[0].pdg == 14) &&
                (sr->mc.nu[0].pdgorig == 14));
      });
    
    const Cut ktNumubarCC( [](const caf::SRProxy* sr){
        if(sr->mc.nnu == 0) return false;
        return (sr->mc.nu[0].iscc && 
                (sr->mc.nu[0].pdg == -14) &&
                (sr->mc.nu[0].pdgorig == -14));
      });
    
    const Cut ktNC( [](const caf::SRProxy* sr){
        if(sr->mc.nnu == 0) return false;
        return (!sr->mc.nu[0].iscc);
      });
    
    
    const Cut ktOther = (!ktSignal && !ktNuebarCC &&
                         !ktNumuCC && !ktNumubarCC &&
                         !ktNC);
    
    
    const Cut ktIsPi0([](const caf::SRProxy* sr)
                      {
                        //float MassOfPi0=0.135;
                        //float kinetic=-10;
                        //float en = -99;
                        
                      if(sr->mc.nnu == 0)   return false;
                      //if(sr->vtx.nelastic == 0) return false;
                      assert(sr->mc.nnu == 1);
                      //if(sr->vtx.elastic[0].fuzzyk.nshwlid==0) return false;
                      //if(sr->vtx.elastic[0].fuzzyk.npng==0) return false;
                     
                      if(sr->mc.nu[0].prim.size() == 0) return false;

                      int nbofprim=sr->mc.nu[0].prim.size();
                      int countpi=0;
                      for(int i = 0; i < nbofprim; i++){
                        if(sr->mc.nu[0].prim[i].pdg == 111){
                          //en= sr->mc.nu[0].prim[i].p.E;
                          //kinetic = en - MassOfPi0;
                          //if (kinetic==0) return false;
                          //if(kinetic > 0.5) countpi++; 
                          countpi++;
                        }
                      }
                      if (countpi==0) return false;
                      if(countpi > 0) return true;
                      return false;
 
                    });


    const Cut ktIsNCPi0 = ktNC && ktIsPi0;
    const Cut ktIsNCNotPi0 = ktNC && !ktIsPi0;
    
    const Cut ktIsNumuCCPi0 = ktNumuCC && ktIsPi0;
    const Cut ktIsNumuCCNotPi0 = ktNumuCC && !ktIsPi0;
    
    const Cut ktIsOther2 = (!ktSignal && !ktNuebarCC &&
                            !ktIsNCPi0 && !ktIsNCNotPi0 && !ktNumubarCC &&
                            !ktIsNumuCCPi0 && !ktIsNumuCCNotPi0);
    
    
    //Spill Truth Cuts
    //----------------
    const TVector3 vtxmin(-130,-140,150);//(-100,-160,150);
    const TVector3 vtxmax(150,140,800);//(160,100,900);
    
    const TVector3* vtxmin_ptr = new TVector3(-130,-140,150);//(-100,-160,150);
    const TVector3* vtxmax_ptr = new TVector3(150,140,800);

    const SpillTruthCut 
      kTrueFiducialST([]
                      (const caf::SRNeutrinoProxy* sr){
                        return (sr->vtx.X() < vtxmax.X() &&
                                sr->vtx.X() > vtxmin.X() &&
                                sr->vtx.Y() < vtxmax.Y() &&
                                sr->vtx.Y() > vtxmin.Y() &&
                                sr->vtx.Z() < vtxmax.Z() &&
                                sr->vtx.Z() > vtxmin.Z() );
                      });
    

    
    const SpillTruthCut kIsNueCCST([](const caf::SRNeutrinoProxy* truth)
                                   {
                                     return (truth->iscc &&
                                             truth->pdg == 12 &&
                                             truth->pdgorig==12);
                                   });
    
    const SpillTruthCut kIsNumuCCST([](const caf::SRNeutrinoProxy* truth)
                                    {
                                      return (truth->iscc &&
                                              abs(truth->pdg) ==  14 &&
                                              abs(truth->pdgorig)== 14);
                                    });
    
    const SpillTruthCut kIsCorrectSignal = kTrueEST >= 1.0 && kTrueEST < 6.0 
      && kTrueElecEST >= 1.0 && kTrueElecEST < 6.0 && 
      kTrueElecCosthetaST >= 0.75 && kTrueElecCosthetaST <= 1.01;
    
    
    const SpillTruthCut kIsTrueSigST = 
      kTrueFiducialST &&  kIsNueCCST; //&& kTrueElecECutST && kTrueElecCosCutST
    //;    && kTrueElectronPhaseSpaceCutST;
  
    const SpillTruthCut kIsTrueNotNumuCCST =
      kTrueFiducialST && !kIsNumuCCST;
    
    const Cut kIsTrueSig = CutFromNuTruthCut(kIsTrueSigST);
    
    
    //Definitions
    struct SelDef{
      std::string name;
      Cut cut;
    };
    
    struct SelDefST{
      std::string name;
      SpillTruthCut cut;
    };
    
    const int kcNumSels = 10;
    const int kcNumChns = 8;
    const int kcNumSels2 = 8;
    
    const SelDef cut_flow2[kcNumSels2] = {
      {"slicing",     kNoCut},
      {"dq",          kcDQ},
      {"fiducial",    kcDQ && kcFiducial},
      {"containment", kcDQ && kcFiducial && kcContainment},
      {"frontplanes", kcDQ && kcFiducial && kcContainment && kcFrontPlanes},
      {"nhits",       kcDQ && kcFiducial && kcContainment && 
       kcFrontPlanes && kcNHits},
      {"muonid",       kcDQ && kcFiducial && kcContainment && kcFrontPlanes && 
       kcNHits && kcMuonIDCut},
      {"cvn",         kcDQ && kcFiducial && kcContainment && kcFrontPlanes && 
       kcNHits && kcMuonIDCut && kcCVN},
    };
    
    const SelDef cut_flow[kcNumSels] = {
      {"slicing",     kNoCut},
      {"dq",          kcDQ},
      {"fiducial",    kcDQ && kcNueCCIncFiducial},
      {"containment", kcDQ && kcNueCCIncFiducial && kcNueCCIncContainment},
      {"frontplanes", kcDQ && kcNueCCIncFiducial && kcNueCCIncContainment
       && kcFrontPlanes},
      {"nhits",       kcDQ && kcNueCCIncFiducial && kcNueCCIncContainment && 
       kcFrontPlanes && kcNHits},
      {"muonid",      kcDQ && kcNueCCIncFiducial && kcNueCCIncContainment && 
       kcFrontPlanes && kcNHits && kcMuonIDCut},
      {"phasespace",  kcPresel && kRecoElectronPhaseSpaceCut},
      {"cvn",         kcPresel && kcCVN},
      {"electronid",  kcPresel && kcElectronIDCut},
    };


    const int kcNumSelsTest = 9;
    const SelDef cut_flow_test[kcNumSelsTest] = {
      {"slicing",     kNoCut},
      {"dq",          kcDQ},
      {"fiducial",    kcDQ && kcNueCCIncFiducial},
      {"containment", kcDQ && kcNueCCIncFiducial && kcNueCCIncContainment},
      {"frontplanes", kcDQ && kcNueCCIncFiducial && kcNueCCIncContainment
       && kcFrontPlanes},
      {"nhits",       kcDQ && kcNueCCIncFiducial && kcNueCCIncContainment && 
       kcFrontPlanes && kcNHits},
      {"muonid",      kcDQ && kcNueCCIncFiducial && kcNueCCIncContainment && 
       kcFrontPlanes && kcNHits && kcRemid},
      {"electronid",  kcDQ && kcNueCCIncFiducial && kcNueCCIncContainment && 
       kcFrontPlanes && kcNHits && kcRemid && kcElectronIDCut},
      {"cvn",         kcDQ && kcNueCCIncFiducial && kcNueCCIncContainment && 
       kcFrontPlanes && kcNHits && kcRemid && kcCVN},
    };


    
    
    
    
    const SelDef chns[kcNumChns] = {
      {"mc",        ktMC},
      {"nuecc",     ktMC && ktSignal && CutFromNuTruthCut(kTrueFiducialST)},
      {"nuebarcc",  ktMC && ktNuebarCC},
      {"numucc",    ktMC && ktNumuCC},
      {"numubarcc", ktMC && ktNumubarCC},
      {"nc",        ktMC && ktNC},
      {"other",     ktMC && ktOther},
      {"nuecc_nonfiducial", ktMC && ktSignal && 
       (!CutFromNuTruthCut(kTrueFiducialST))},
    };


    const SelDef chns_data[kcNumChns+1] = {
      {"data",      !ktMC},
      {"mc",        ktMC},
      {"nuecc",     ktMC && ktSignal && CutFromNuTruthCut(kTrueFiducialST)},
      {"nuebarcc",  ktMC && ktNuebarCC},
      {"numucc",    ktMC && ktNumuCC},
      {"numubarcc", ktMC && ktNumubarCC},
      {"nc",        ktMC && ktNC},
      {"other",     ktMC && ktOther},
      {"nuecc_nonfiducial", ktMC && ktSignal && 
       (!CutFromNuTruthCut(kTrueFiducialST))},
    };

    
    
    const int kcNumChnsExtra = 12;
    
    const SelDef chns_extra[kcNumChnsExtra]={
      {"data",      !ktMC},
      {"mc",             ktMC},
      {"nuecc",          ktMC && ktSignal},
      {"nuebarcc",       ktMC && ktNuebarCC},
      {"numucc_pi0",     ktMC && ktIsNumuCCPi0},
      {"numucc_not_pi0", ktMC && ktIsNumuCCNotPi0},
      {"numubarcc",      ktMC && ktNumubarCC},
      {"nc_pi0",         ktMC && ktIsNCPi0},
      {"nc_not_pi0",     ktMC && ktIsNCNotPi0},
      {"other",          ktMC && ktIsOther2},
      {"non_pi0_background", ktMC && (ktIsNumuCCNotPi0 || ktNumubarCC ||
                                      ktIsNCNotPi0 || ktIsOther2)},
      {"background",         ktMC && !ktSignal},
    };
    
    
    const int kcNumIntChnsST = 7;
    const int kcNumIntChns = 8;
    
    const SelDef interaction_chns[kcNumIntChns] = {
      {"nocut", ktMC},
      {"cc",   !kIsNC},
      {"qe",   kIsQE},
      {"res",  kIsRes},
      {"dis",  kIsDIS},
      {"coh",  kIsCoh},
      {"mec",  kIsDytmanMEC},
      {"other", !kIsQE && !kIsRes && !kIsDIS && !kIsCoh && !kIsDytmanMEC},
    };
    
    const SelDefST interaction_chns_st[kcNumIntChnsST] = {
      {"cc", kIsTrueSigST},
      {"qe", kIsTrueSigST && kIsQE_NT},
      {"res", kIsTrueSigST && kIsRes_NT},
      {"dis", kIsTrueSigST && kIsDIS_NT},
      {"coh", kIsTrueSigST && kIsCoh_NT},
      {"mec", kIsTrueSigST && kIsDytmanMEC_NT},
      {"other", kIsTrueSigST && !kIsQE_NT && 
       !kIsRes_NT && !kIsDIS_NT && ! kIsCoh_NT && !kIsDytmanMEC_NT},
    };

    //Cuts For ReMID Cut Optimization
    //-------------------------------
    const int kcReMIDSels   = 3;
    const int kcReMIDSelsST = 1;
    
    const SelDef remid_selection[kcReMIDSels] = {
      {"presel",           kcPresel_loose},
      {"presel_signal",    kcPresel_loose && ktMC && !ktNumuCC && !ktNumubarCC},
      {"presel_background",kcPresel_loose && ktMC && (ktNumuCC || ktNumubarCC)}, 
    };
    
    const SelDefST remid_selection_st[kcReMIDSelsST] = {
      {"spilltruth", kIsTrueNotNumuCCST},
    };
    
    //Cuts For CVN Cut Optimization
    //-----------------------------
    const int kcCvnSels   = 3;
    const int kcCvnSelsST = 1;
    
    const SelDef cvn_selection[kcCvnSels] = {
      {"presel",  kcPresel},
      {"presel_signal", kcPresel && ktMC && ktSignal},
      {"presel_background", kcPresel && ktMC && !ktSignal},
    };
    
    const SelDefST cvn_selection_st[kcCvnSelsST] = {
      {"spilltruth", kIsTrueSigST},
    };

    //Analysis Cuts
    const Cut kPresel2D = 
      kcPresel;
    const Cut kPresel1D = kcPresel && 
      kRecoElectronPhaseSpaceCut;
    
    const NuTruthCut kTrueSignal2DST = kIsTrueSigST;
    const NuTruthCut kTrueSignal1DST = kIsTrueSigST &&
      kTrueElectronPhaseSpaceCutST;
    
    const Cut kTrueSignal2D = CutFromNuTruthCut(kTrueSignal2DST);
    const Cut kTrueSignal1D = CutFromNuTruthCut(kTrueSignal1DST);
    
    //Define Signal Correctly For 1D Result
    const SelDef chns_1d[kcNumChns] = {
      {"mc",        ktMC},
      {"nuecc",     ktMC && ktSignal && CutFromNuTruthCut(kTrueFiducialST) 
       && CutFromNuTruthCut(kTrueSignal1DST)},
      {"nuebarcc",  ktMC && ktNuebarCC},
      {"numucc",    ktMC && ktNumuCC},
      {"numubarcc", ktMC && ktNumubarCC},
      {"nc",        ktMC && ktNC},
      {"other",     ktMC && ktOther},
      {"nuecc_nonfiducial", ktMC && ktSignal && 
       (!CutFromNuTruthCut(kTrueFiducialST) || !CutFromNuTruthCut(kTrueSignal1DST)) },
    };

  }//end namespace nueccinc
} //end namespace ana