DataMCComp_header.h

Go to the documentation of this file.

#pragma once

#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Prod5Loaders.h"

#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/RunPeriods.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Utilities.h"

#include "CAFAna/Cuts/AnalysisMasks.h"
#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts2020.h"
#include "3FlavorAna/Cuts/QuantileCuts.h"
#include "3FlavorAna/Cuts/QuantileCuts2020.h"
#include "CAFAna/Cuts/PeriodCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "CAFAna/Cuts/TruthCuts.h"

#include "CAFAna/Vars/BPFVars.h"
#include "CAFAna/Vars/CVNProngVars.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/HistAxes.h"
#include "3FlavorAna/Vars/Binnings.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Vars/XsecTunes.h"

#include "CAFAna/Prediction/PredictionInterp.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Prediction/PredictionGenerator.h"

#include "StandardRecord/Proxy/SRProxy.h"

#include "TGraphAsymmErrors.h"
#include "TFeldmanCousins.h"


#include "/nova/app/users/asutton/3Flavor/Ana2020/BoxOpening_br_ana/plotting_header.h"

using namespace ana;



// Function to make data graph with error bars
// same as the one in CAFAna/nue/Ana2018/sideband/fd_plot.h
// but including that caused some conflicts that I didn't want to debug
TGraphAsymmErrors* GraphWithPoissonErrors2(const TH1* h, const TH1* h2, bool noErrorsXaxis, bool drawEmptyBins)
{
  TGraphAsymmErrors* gr = new TGraphAsymmErrors(h);

  TFeldmanCousins fc(0.6827);//1 sigma

  std::vector<int> removePoints; // keep track of the points to remove
  for(int i = 0; i < h->GetNbinsX(); ++i){
    double x, y;
    gr->GetPoint(i, x, y);

    if ( (drawEmptyBins && h2->GetBinContent(i+1) > 0.05) || y>=1 ){
      if ( y < 50 )   gr->SetPointEYlow(i, y-fc.CalculateLowerLimit(y,0));
      else            gr->SetPointEYlow(i,sqrt(y));
      if ( y < 30 )   gr->SetPointEYhigh(i,fc.CalculateUpperLimit(y,0)-y);
      else            gr->SetPointEYhigh(i,sqrt(y));
    } else {
      removePoints.push_back(i);
    }
    
    if(noErrorsXaxis){
      gr->SetPointEXlow(i, 0);
      gr->SetPointEXhigh(i, 0);
    } // Do not use bin width as X-error for points
  }

  // get rid of uninteresting points
  for (uint i = 0; i < removePoints.size(); ++i)
    gr->RemovePoint(removePoints[i]-i);

  return gr;
}

// These are also defined elsewhere but I had some issue with includes
// also only need them to select periods for debugging
int PeriodFromHeader(const caf::SRHeaderProxy& hdr)
  {
    if(hdr.det == caf::kNEARDET){
      return PeriodFromRunND(hdr.run);
    }
    else if(hdr.det == caf::kFARDET){
      return PeriodFromRunFD(hdr.run);
    }
    else{
      std::cout << "PeriodFromHeader(): Unknown detector "
                << hdr.det << ". Abort" << std::endl;
      abort();
    }
  }

const Cut kPeriod10([](const caf::SRProxy* sr) {
    return PeriodFromHeader(sr->hdr) == 10;
  });

// --- Cut for kalbak study.
const Cut kKalBakCut( [](const caf::SRProxy* sr) {
    return (sr->sel.contain.kalbakcell > 40 &&
            sr->sel.contain.kalbakcell < 88);
  });


// ************** Define my additional Vars **************
// --- The Reco - True / True neutrino energy
const Var kERes = ( (kCCE-kTrueE)/kTrueE );

// --- Kalman track E per nHit.
const Var kTrkEPerNHit( [](const caf::SRProxy* sr) {
    return (sr->trk.kalman.tracks[0].calE / sr->trk.kalman.tracks[0].nhit);
  });

// --- The most forward Kalman track hit cell.
const Var kKalFwdCell( [](const caf::SRProxy* sr) { return sr->sel.contain.kalfwdcell; });

// --- The most backward Kalman track hit cell.
const Var kKalBakCell( [](const caf::SRProxy* sr) { return sr->sel.contain.kalbakcell; });

// --- The highest Y vertex of Kalman track.
const Var kMaxKalYPos( [](const caf::SRProxy* sr) {
    if(sr->trk.kalman.ntracks < 1) return -10.0f;
    return std::max( sr->trk.kalman.tracks[0].start.Y()/100, 
                     sr->trk.kalman.tracks[0].stop.Y()/100 );
  });

// --- The highest Z vertex of Kalman track.
const Var kMaxKalZPos( [](const caf::SRProxy* sr) {
    if(sr->trk.kalman.ntracks < 1) return -10.0f;
    return std::max( sr->trk.kalman.tracks[0].start.Z()/100, 
                     sr->trk.kalman.tracks[0].stop.Z()/100 );
  });

// --- Scattering over track length
const Var kScattTrLen( [](const caf::SRProxy* sr) {
    if(sr->trk.kalman.ntracks < 1) return -10.0f;
    return (float)((float)sr->sel.cosrej.scatt / (float)sr->trk.kalman.tracks[0].len);
  });

// --- Min Proj cells to edge
const Var kMinCellEdg( [](const caf::SRProxy* sr) {
    return std::min((sr->sel.contain.kalfwdcell + sr->sel.contain.kalbakcell),
                    (sr->sel.contain.cosfwdcell + sr->sel.contain.cosbakcell));
  });

// --- Ratio of kalman hits in slice.
const Var kRatOfKalHi( [](const caf::SRProxy* sr) {
    if(sr->trk.kalman.ntracks < 1) return -10.0f;
    return (float)((float)sr->trk.kalman.tracks[0].nhit / (float)sr->slc.nhit);
  });

// --- The first hit cell.
const Var kFirstCell( [](const caf::SRProxy* sr) { return sr->slc.firstcell; });  

// --- The last hit cell.
const Var kLastCell ( [](const caf::SRProxy* sr) { return sr->slc.lastcell;  });

// --- The difference in XY Track Length
const Var kTrkLenXY( [](const caf::SRProxy* sr) {
    if(sr->trk.kalman.ntracks < 1) return -10.0f;
    float XLen = ( sr->trk.kalman.tracks[0].start.X() - sr->trk.kalman.tracks[0].stop.X() ) /100;
    float YLen = ( sr->trk.kalman.tracks[0].start.Y() - sr->trk.kalman.tracks[0].stop.Y() ) /100;
    return XLen - YLen;
  });

// --- The True Hadronic Energy
const Var kTrueHadE = ( kTrueE - kTrueMuonE );

// --- The True Hadronic Energy Fraction
const Var kTrueHFrE( [](const caf::SRProxy* sr) {
    float HadE = kTrueHadE(sr);
    float TruE = kTrueE(sr);
    return ( HadE / TruE );
  });

// ******* Define some bin sets which I want to use *******

const Binning NumHitBins   = Binning::Simple( 20,  0.0  , 400  );
const Binning RatioBins    = Binning::Simple( 51, -0.01 , 1.01 );
const Binning CVNLooseBins = Binning::Simple( 31,  0.695, 1.005 );
const Binning CosRejBins   = Binning::Simple( 35,  0.3  , 0.80 );

// Energy bins
const Binning OptEnBins   = kNumuCCEOptimisedBinning;
const Binning EnergyBins  = Binning::Simple( 100, 0   , 5    );
const Binning LowEnBins   = Binning::Simple(  25, 0   , 3    );
const Binning EnPHitBins  = Binning::Simple(  40, 0   , 0.04 );
const Binning EResBins    = Binning::Simple( 100, -0.5, 0.5  );
const Binning HadFracBins = Binning::Simple( 26, -0.02, 1.02 );

// Length / hit bins
const Binning NumTrkBins  = Binning::Simple( 11, -0.5,   10.5 );
const Binning YRangeBins  = Binning::Simple( 24, -8  ,    8   );
const Binning ZRangeBins  = Binning::Simple( 25,  0  ,   60   );
const Binning XYDiffBins  = Binning::Simple( 20,  0  ,   10   );
const Binning LengthBins  = Binning::Simple( 25,  0  ,   20   );
const Binning LgCellsBins = Binning::Simple( 25,  0  , 1000   );
const Binning SmCellsBins = Binning::Simple( 25,  0  ,  400   );

// Time bins
const Binning TimeFBins  = Binning::Simple( 100,   0, 550 );
const Binning TimeZBins  = Binning::Simple(  50, 200, 250 );

// ******* Prediction struct *******
struct PredDef
{
  const IPrediction * pred;
  const TString xp_type;
  const TString syst_name;
  const TString sigma_name;
  const TString edef_name;
};


// ******* Make the variable map *******
// Using pair for shortname in directories and files
// and longname for x-axis title
struct Variable
{
  Var var;
  Binning bin;
};

const std::map< std::pair<std::string, std::string >, Variable > vars =
  {
    // Interesting things
    { std::make_pair("recoE"         , "Reconstructed Neutrino Energy (GeV)"), {kCCE           , kNumuCCEOptimisedBinning} },
    { std::make_pair("sliceTime_Full", "Time of slice - full (us)"          ), {kSliceTime     , TimeFBins               } },
    { std::make_pair("sliceTime_Zoom", "Time of slice - zoom (us)"          ), {kSliceTime     , TimeZBins               } },
    { std::make_pair("muonE"         , "Reconstructed muon energy (GeV)"    ), {kMuE           , LowEnBins               } },
    { std::make_pair("muonEPerHit"   , "Muon energy per hit (GeV/hit)"      ), {kTrkCalEPerNHit, EnPHitBins              } },
    { std::make_pair("hadE"          , "Reconstructed hadronic energy (GeV)"), {kHadE          , LowEnBins               } },
    { std::make_pair("hadEPerHit"    , "Hadronic energy per hit (GeV/hit)"  ), {kHadEPerNHit   , EnPHitBins              } },
    { std::make_pair("hadEFrac"      , "Hadronic energy fraction in slice"  ), {kHadEFrac      , RatioBins               } },
    { std::make_pair("sliceNHits"    , "Number of hits per slice"           ), {kNHit          , NumHitBins              } },
    { std::make_pair("remID"         , "RemID score"                        ), {kRemID         , RatioBins               } },
    { std::make_pair("cvnLoose"      , "CVN loose p_{T}^{sum} score"        ), {kCVNm_looseptp , CVNLooseBins            } },
    { std::make_pair("cosRej"        , "Cosmic rejection score"             ), {kNumuContPID   , CosRejBins              } },
    
    // Track Positions
    { std::make_pair("trackStartX"   , "Track start X position (m)"        ), {kTrkStartX, YRangeBins} },
    { std::make_pair("trackStartY"   , "Track start Y position (m)"        ), {kTrkStartY, YRangeBins} },
    { std::make_pair("trackStartZ"   , "Track start Z position (m)"        ), {kTrkStartZ, ZRangeBins} },
    { std::make_pair("trackEndX"     , "Track stop X position (m)"         ), {kTrkEndX  , YRangeBins} },
    { std::make_pair("trackEndY"     , "Track stop Y position (m)"         ), {kTrkEndY  , YRangeBins} },
    { std::make_pair("trackEndZ"     , "Track stop Z position (m)"         ), {kTrkEndZ  , ZRangeBins} },
    { std::make_pair("trackXYLenDiff", "Difference in XY track lengths (m)"), {kTrkLenXY , XYDiffBins} },

    // Track stuff
    { std::make_pair("kalNTracks"      , "Number of Kalman tracks"            ), {kNKalman   , NumTrkBins  } },
    { std::make_pair("kalTrackNHits"   , "Number of hits per Kal track"       ), {kTrkNhits  , NumHitBins  } },
    { std::make_pair("kalHitSliceRatio", "Fraction of Kal track hits in slice"), {kRatOfKalHi, RatioBins   } },
    { std::make_pair("kalTrackLen"     , "Length of Kalman track"             ), {kTrkLength , LengthBins  } },
    { std::make_pair("kalNumiAngle"    , "Angle of Kalman track to beam"      ), {kCosNumi   , RatioBins   } },
    { std::make_pair("kalFwdCell"      , "Kalman projection forward cells"    ), {kKalFwdCell, LgCellsBins } },
    { std::make_pair("kalBakCell"      , "Kalman projection backward cells"   ), {kKalBakCell, LgCellsBins } },
    { std::make_pair("minCellEdge"     , "Minimum projected cells to edge"    ), {kMinCellEdg, LgCellsBins } },

    // Slice stuff
    { std::make_pair("firstCell"  , "Slice first cell"), {kFirstCell , SmCellsBins } },
    { std::make_pair("lastCell"   , "Slice last cell" ), {kLastCell  , SmCellsBins } },
    { std::make_pair("sliceMaxY"  , "Slice max Y (m)" ), {kSlcMaxY   , YRangeBins  } },
    { std::make_pair("sliceMinY"  , "Slice min Y (m)" ), {kSlcMinY   , YRangeBins  } },
    { std::make_pair("sliceMaxX"  , "Slice max X (m)" ), {kSlcMaxX   , YRangeBins  } },
    { std::make_pair("sliceMinX"  , "Slice min X (m)" ), {kSlcMinX   , YRangeBins  } },
    { std::make_pair("sliceMaxZ"  , "Slice max Z (m)" ), {kSlcMaxZ   , ZRangeBins  } },
    { std::make_pair("sliceMinZ"  , "Slice min Z (m)" ), {kSlcMinZ   , ZRangeBins  } },
        
    // Truth stuff
    { std::make_pair("trueE", "True neutrino energy (GeV)"), {kTrueE, EnergyBins} },
    { std::make_pair("resE" , "(RecoE - TrueE) / TrueE"   ), {kERes , EResBins  } }
  };