NCPi0SemiIncCuts.h

Go to the documentation of this file.

#pragma once
#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Vars/Vars.h"
#include "StandardRecord/Proxy/SRProxy.h"

#include "NDAna/ncpi0_SemiInclusive/NCPi0SemiIncVars.h"


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

namespace ana
{
                   
  const Cut kIs1Vtx  ([](const caf::SRProxy* sr)
                      {
                        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.npng==0) return false;
                        if(sr->vtx.elastic[0].fuzzyk.nshwlid==0) return false;

                        return (sr->vtx.nelastic==1);
                      });

  const Cut kIs2Prongs ([](const caf::SRProxy* sr)
                        {
                    
                          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.npng==0) return false;
                          if(sr->vtx.elastic[0].fuzzyk.nshwlid==0) return false;
                   
                          return (sr->vtx.elastic[0].fuzzyk.npng == 2 );
                        });


  const Cut kVtxInFiducial([](const caf::SRProxy* sr){
                        
      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->vtx.elastic[0].vtx.x < -130) return false;
      if (sr->vtx.elastic[0].vtx.x >  160) return false;
      if (sr->vtx.elastic[0].vtx.y < -150) return false;
      if (sr->vtx.elastic[0].vtx.y >  120) return false;
      if (sr->vtx.elastic[0].vtx.z <  225) return false;
      if (sr->vtx.elastic[0].vtx.z > 950) return false;
                
      return true;
    });


  const Cut kShwStopContainment([](const caf::SRProxy* sr)
                                {
                                  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;
                    
                    
                                  for (unsigned int i=0; i< sr->vtx.elastic[0].fuzzyk.nshwlid; i++){
                     
                                    if ((sr->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.x) < -130.0)  return false;
                                    if ((sr->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.x) > 140.0)  return false;
                                    if ((sr->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.y) < -150.0)  return false;
                                    if ((sr->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.y) > 150.0)  return false;
                                    if ((sr->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.z) < 300.0)  return false;
                                    if ((sr->vtx.elastic[0].fuzzyk.png[i].shwlid.stop.z) > 1025.0)  return false;
                                  }
                                  return true;
                                });




  const SpillTruthCut kIsMyNumuST([](const caf::SRNeutrinoProxy* nu)
                                  {
                                    return(nu->pdg == 14 && nu->pdgorig == 14);
                                  });
  
  const Cut kIsMyNumu  = CutFromSpillTruthCut(kIsMyNumuST);

  const SpillTruthCut kIsMyNueST([](const caf::SRNeutrinoProxy* nu)
                                 {
                                   return(nu->pdg == 12 && nu->pdgorig == 12);
                                 });
  
  const Cut kIsMyNue  = CutFromSpillTruthCut(kIsMyNueST);

  const SpillTruthCut kIsNCurrentST([](const caf::SRNeutrinoProxy* nu)
                                    {
                     
                                      return(nu->iscc == 0);
                                    });

  const Cut kIsNCurrent  = CutFromSpillTruthCut(kIsNCurrentST);


  const SpillTruthCut kIsCCurrentST([](const caf::SRNeutrinoProxy* nu)
                                    {
                     
                                      return(nu->iscc == 1);
                                    });

  const Cut kIsCCurrent  = CutFromSpillTruthCut(kIsCCurrentST);
 
  const SpillTruthCut kIsPi0ST([](const caf::SRNeutrinoProxy* nu)
                               {
                                 float MassOfPi0=0.135;
                                 float kinetic=-10;
                                 float en = -99;

                                 if(nu->prim.size() == 0) return false;
                     
                                 //int countpi=0;                    
                                 //float maximum=0;
                                 int nbofprim=nu->prim.size();

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

                               });

  const Cut kIsPi0  = CutFromSpillTruthCut(kIsPi0ST);

  const SpillTruthCut kIsAllPi0ST([](const caf::SRNeutrinoProxy* nu)
                                  {
                                    float MassOfPi0=0.135;
                                    float kinetic=-10;
                                    float en = -99;

                                    if(nu->prim.size() == 0) return false;
                     
                                    //int countpi=0;
                                    //float maximum=0;
                                    int nbofprim=nu->prim.size();

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

                                  });

  const Cut kIsAllPi0  = CutFromSpillTruthCut(kIsAllPi0ST);

   
  const Cut kIsReMId ([](const caf::SRProxy* sr)
                      {
                        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->sel.remid.pid < 0) return false;
                        return (sr->sel.remid.pid < 0.36);
                      });


  const TVector3 vtxmin(-130,-150,225);
  const TVector3 vtxmax(160,120,950);

  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 Cut kTrueVtxInFiducial  = CutFromSpillTruthCut(kTrueFiducialST);


  const SpillTruthCut kModeDIS ([] (const caf::SRNeutrinoProxy* sr)
                                {
                        
                                  return (sr->mode==caf::kDIS);
                
                                });

  const Cut kModeTypeDIS   = CutFromSpillTruthCut(kModeDIS);

  const SpillTruthCut kModeRES([](const caf::SRNeutrinoProxy* sr)
                               {
  
                                 return (sr->mode==caf::kRes);

                               });

  const Cut kModeTypeRES   = CutFromSpillTruthCut(kModeRES);


  const SpillTruthCut kModeCOH([](const caf::SRNeutrinoProxy* sr)
                               {
                       
                                 return (sr->mode==caf::kCoh);
                               });

  const Cut kModeTypeCOH   = CutFromSpillTruthCut(kModeCOH);


  const SpillTruthCut kModeQE([](const caf::SRNeutrinoProxy* sr)
                              {

                                return (sr->mode==caf::kQE);

                              });

  const Cut kModeTypeQE   = CutFromSpillTruthCut(kModeQE);

  const Cut kcut15hits([](const caf::SRProxy* sr){
                        
      if(sr->mc.nnu==0) return false;
      assert(sr->mc.nnu == 1);
      if(sr->vtx.nelastic == 0) return false;
      if(sr->vtx.elastic[0].fuzzyk.nshwlid==0) return false;
      if(sr->vtx.elastic[0].fuzzyk.npng==0) return false;

      float hit;
      hit=sr->slc.nhit;
      if (hit>15) return true;
                       
      return false;
    });


  const Cut ktMC( [](const caf::SRProxy* sr){
      if(sr->hdr.ismc) return true;
      return false;
    });

  //Since we consider "Neutrino induced" NC interactions so, neutrino can be numu or nue in forward horn current beam
  const SpillTruthCut kIsMyNeutrinoST =  (kIsMyNumuST||kIsMyNueST);
  const Cut kIsMyNeutrino =  (kIsMyNumu||kIsMyNue);


  const Cut kIsPreSel = kIs1Vtx && kIs2Prongs && kVtxInFiducial && kShwStopContainment && kIsReMId && kcut15hits;
  const Cut kIsSigDef = kIsMyNeutrino && kIsNCurrent && kIsPi0;

  //Signal in true fiducial
  const SpillTruthCut kIsSignalST = kIsMyNeutrinoST && kIsNCurrentST && kIsPi0ST && kTrueFiducialST;
  const Cut kIsSignal = kIsMyNeutrino && kIsNCurrent && kIsPi0 && kTrueVtxInFiducial;

  //Signal without kinmetic constraints, to be used in resolution studies where we don't apply anu cut on trueKE of a pi0
  const SpillTruthCut kIsSignalNoEnThrsST = kIsMyNeutrinoST && kIsNCurrentST && kIsAllPi0ST && kTrueFiducialST;
  const Cut kIsSignalNoEnThrs  = kIsMyNeutrino && kIsNCurrent && kIsAllPi0 && kTrueVtxInFiducial;


  //Total background
  const Cut kIsBkgNoEnThrs = !kIsSignalNoEnThrs; 
  const Cut kIsBkg         = !kIsSignal; 

  //Signal with interaction types. PreSel cuts are not included in these defs.
  const Cut kIsRESSignal =  kIsSignal && kModeTypeRES;
  const Cut kIsDISSignal =  kIsSignal && kModeTypeDIS;
  const Cut kIsOtherSignal =  kIsSignal && !(kModeTypeDIS && kModeTypeRES);

  //Bkg with interaction types.  PreSel cuts are not included in these defs.
  const Cut kIsNCBkg = kIsBkg && kIsNCurrent;
  const Cut kIsCCBkg = kIsBkg && kIsCCurrent;
  const Cut kIsRESBkg = kIsBkg && kModeTypeRES;
  const Cut kIsDISBkg = kIsBkg && kModeTypeDIS;
  const Cut kIsOtherBkg = kIsBkg && !(kModeTypeDIS && kModeTypeRES);
  const Cut kIsCCRESBkg = kIsBkg && kIsCCurrent && kModeTypeRES;
  const Cut kIsNCRESBkg = kIsBkg && kIsNCurrent && kModeTypeRES;
  const Cut kIsCCDISBkg = kIsBkg && kIsCCurrent && kModeTypeDIS;
  const Cut kIsNCDISBkg = kIsBkg && kIsNCurrent && kModeTypeDIS;
  const Cut kIsCCOtherBkg = kIsBkg && kIsCCurrent && !(kModeTypeDIS && kModeTypeRES);
  const Cut kIsNCOtherBkg = kIsBkg && kIsNCurrent && !(kModeTypeDIS && kModeTypeRES);


  struct SelDef{
    std::string name;
    Cut cut;
  };

  const int knumcutflow = 5;
  const int knumchns    = 18;

  const SelDef cut_flow[knumcutflow] = { 
    {"Fid",   kIs1Vtx && kVtxInFiducial },
    {"Cont",  kIs1Vtx && kVtxInFiducial && kShwStopContainment},
    {"Pngs",  kIs1Vtx && kVtxInFiducial && kShwStopContainment && kIs2Prongs},
    {"Remid", kIs1Vtx && kVtxInFiducial && kShwStopContainment && kIs2Prongs && kIsReMId},
    {"PreSel", kIsPreSel},
  };

  const SelDef chns[knumchns] = {
    {"TotMC",   ktMC},
    {"Sigall",   ktMC && kIsSignalNoEnThrs},
    {"Sig",     ktMC && kIsSignal},
    {"RESSig",   ktMC &&  kIsRESSignal},
    {"DISSig",   ktMC &&  kIsDISSignal},
    {"OthrSig",  ktMC && kIsOtherSignal},
    {"Bkg",     ktMC && kIsBkg},
    {"CCBkg",   ktMC && kIsCCBkg }, 
    {"NCBkg",   ktMC && kIsNCBkg },
    {"CCRESBkg",   ktMC && kIsCCRESBkg}, 
    {"NCRESBkg",   ktMC && kIsNCRESBkg},
    {"CCDISBkg",   ktMC && kIsCCDISBkg}, 
    {"NCDISBkg",   ktMC && kIsNCDISBkg},
    {"CCOtherBkg",   ktMC && kIsCCOtherBkg}, 
    {"NCOtherBkg",   ktMC && kIsNCOtherBkg },
    {"RESBkg",  ktMC && kIsRESBkg },
    {"DISBkg",  ktMC && kIsDISBkg},
    {"OthrBkg", ktMC && kIsOtherBkg},
  };




}//end of namespace