make_nue_pidtuning.C

Go to the documentation of this file.

#include "OscLib/OscCalcDumb.h"
#include "CAFAna/Analysis/Calcs.h"

#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Core/SystShifts.h"
#include "CAFAna/Core/LoadFromFile.h"



#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "3FlavorAna/Cuts/NueCuts2017.h"
#include "CAFAna/Cuts/SpillCuts.h"


#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Vars/FitVars.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/Vars/HistAxes.h"

#include "CAFAna/Prediction/PredictionNoExtrap.h"

#include "Utilities/func/MathUtil.h"

#include "TCanvas.h"
#include "TH2.h"
#include "TH1.h"
#include "TFile.h"
#include "TLegend.h"

#include <iostream>
#include <cmath>

using namespace ana;



// Quick numbers, POT, livetime, and osc parameters
// livetime/pot from Lyudmila's checks

double pot = 8.0e20;
double lt = 199;

// Change these here to look at cut stability if you want
double dCP = 1.2*TMath::Pi();
double dmsq32 = 0.00244;
double ssth23 = 0.56;
double ss2th13= 0.082;



// Need a way to translate CVN-BDT 1D plot into a TH2
TH2 *TH1ToTH2(TH1 *h)
{
  TH2 *ret = new TH2F(UniqueName().c_str(),"",50,0.5,1,50,0.35,0.85);
  for (int i = 1; i <= 2500; i++){
    int xbin = (i-1)%50 + 1;
    int ybin = (i-1)/50 + 1;
    ret->SetBinContent(xbin,ybin,h->GetBinContent(i));
  }
  return ret;
}

void PeripheralCuts(std::string name)
{
  std::string fname = "out_fdpred.root";
  PredictionNoExtrap *pred =
    LoadFromFile<PredictionNoExtrap>(fname,"pred"+name).release();
  Spectrum *cos = LoadFromFile<Spectrum>(fname,"cos"+name).release();

  // Setup oscillations
  osc::IOscCalcAdjustable* calc = DefaultOscCalc();

  kFitDeltaInPiUnits.SetValue(calc,dCP/TMath::Pi());
  kFitDmSq32.SetValue(calc,dmsq32);
  kFitSinSqTheta23.SetValue(calc,ssth23);
  kFitSinSq2Theta13.SetValue(calc,ss2th13);

  // Make signal, background spectra
  Spectrum sig = pred->PredictComponent(calc,Flavors::kNuMuToNuE,
                                        Current::kCC, Sign::kAntiNu);
  Spectrum bb  = pred->Predict(calc);
  bb -= sig;

  TH2 *hbb  = TH1ToTH2(bb .ToTH1(pot));
  TH2 *hsig = TH1ToTH2(sig.ToTH1(pot));
  TH2 *hcos = TH1ToTH2(cos->ToTH1(lt,kLivetime));

  TH2 *htot = TH1ToTH2(cos->ToTH1(lt,kLivetime));
  htot->Add(hbb);
  htot->Add(hsig);

  std::vector<double> sigs(50,0);
  std::vector<double> tots(50,0);

  // scan over CVN bins
  for (int i = 1; i <= 50; i++){
    double fom = 0;
    double bestY = 0;
    // Calculate a best BDT cut
    for (int j = 1; j <= 50; j++){
      double sig = hsig->Integral(i,i,j,50);
      double tot = htot->Integral(i,i,j,50);
      if (tot == 0) continue;
      double cfom = sig/sqrt(tot);
      if (cfom > fom){
        fom = cfom;
        bestY = 0.35 + 0.01*(j-1);
      }
    }
    if (i >=46)
      std::cout << 0.49 + 0.01*i << "  " << bestY << "  " << fom << std::endl;
    //foms.push_back(fom);
    for (int j = 0; j < i; j++){
      // Accumalate sig / tot events above this CVN cut, and BDT cuts
      sigs[j] += hsig->Integral(i,i,(bestY-0.35)/0.01+1,50);
      tots[j] += htot->Integral(i,i,(bestY-0.35)/0.01+1,50);
    }
  }

  TH1 *h = new TH1F(UniqueName().c_str(),"",50,0.5,1); // FOM vs CVN plot
  for (int i = 50; i >= 1; i--){
    h->SetBinContent(i,sigs[i-1]/sqrt(tots[i-1]));
  }

  h->GetXaxis()->SetTitle("CVN-SS");
  h->GetYaxis()->SetTitle("Total FOM");
  h->GetXaxis()->CenterTitle();
  h->GetYaxis()->CenterTitle();

  TCanvas *can = new TCanvas();
  h->Draw();

}

void CVNCuts(std::string name)
{
  std::string fname = "out_fdpred.root";
  PredictionNoExtrap *pred =
    LoadFromFile<PredictionNoExtrap>(fname,"pred"+name).release();
  Spectrum *cos = LoadFromFile<Spectrum>(fname,"cos"+name).release();

  osc::IOscCalcAdjustable* calc = DefaultOscCalc();

  kFitDeltaInPiUnits.SetValue(calc,dCP/TMath::Pi());
  kFitDmSq32.SetValue(calc,dmsq32);
  kFitSinSqTheta23.SetValue(calc,ssth23);
  kFitSinSq2Theta13.SetValue(calc,ss2th13);


  Spectrum sig = pred->PredictComponent(calc,Flavors::kNuMuToNuE,
                                        Current::kCC, Sign::kAntiNu);
  Spectrum bb  = pred->Predict(calc);
  bb -= sig;

  TH1 *hbb  = bb.ToTH1(pot);
  TH1 *hcos = cos->ToTH1(lt,kLivetime);
  TH1 *hsig = sig.ToTH1(pot);

  TH1 *htot = cos->ToTH1(lt,kLivetime);
  htot->Add(hbb);
  htot->Add(hsig);

  double fom = -1;
  // Bin boundaries for CVN
  double best1 = -1;
  double best2 = -1;
  double best3 = -1;

  std::cout << "Finding optimal FOM for 1 CVN bin analysis" << std::endl;
  for (int i = 1; i <= 50; i++){
    double sig1 = hsig->Integral(i,50);
    double tot1 = htot->Integral(i,50);
    double fom1 = sig1/sqrt(tot1);
    //double fom1 = sig1/sqrt(tot1);
    if (fom1 > fom){
      fom = fom1;
      best1 = 0.5 + 0.01*(i-1);
    }
  }
  std::cout << fom << ":  " << best1 << "-" << "1" << std::endl;


  std::cout << "Finding optimal FOM for 2 CVN bin analysis" << std::endl;
  fom = -1; // reset fom
  for (int i = 1; i <= 49; i++){ // Can't go all the way to 1
    for (int j = i+1; j <= 50; j++){ // Must start above 1 bin max
      double sig1 = hsig->Integral(i,j-1);
      double tot1 = htot->Integral(i,j-1);
      double fom1 = sig1/sqrt(tot1);
      double sig2 = hsig->Integral(j,50);
      double tot2 = htot->Integral(j,50);
      double fom2 = sig2/sqrt(tot2);
      if (sqrt(fom1*fom1+fom2*fom2) > fom){
        fom = sqrt(fom1*fom1+fom2*fom2);
        best1 = 0.5 + 0.01*(i-1);
        best2 = 0.5 + 0.01*(j-1);
      }
    }
  }
  std::cout << fom << ":  " << best1 << "-" << best2 << "-1" << std::endl;


  std::cout << "Finding optimal FOM for 2 CVN bin analysis" << std::endl;
  fom = -1; // Reset fom
  for (int i = 1; i <= 48; i++){
    for (int j = i+1; j <= 49; j++){
      for (int k = j+1; k <= 50; k++){
        double sig1 = hsig->Integral(i,j-1);
        double tot1 = htot->Integral(i,j-1);
        double fom1 = sig1/sqrt(tot1);

        double sig2 = hsig->Integral(j,k-1);
        double tot2 = htot->Integral(j,k-1);
        double fom2 = sig2/sqrt(tot2);

        double sig3 = hsig->Integral(k,50);
        double tot3 = htot->Integral(k,50);
        double fom3 = sig2/sqrt(tot2);

        double tfom = sqrt(fom1*fom1+fom2*fom2+fom3*fom3);
        if (tfom>fom){
          fom = tfom;
          best1 = 0.5 + 0.01*(i-1);
          best2 = 0.5 + 0.01*(j-1);
          best3 = 0.5 + 0.01*(k-1);
        }
      }
    }
  }
  std::cout << fom << ":  " << best1 << "-" << best2 << "-" << best3 << "-1" << std::endl;

}






// TODO: move these cuts to a header, once PID cuts are finalized

const Cut kNue2018CVNVsCosPID([](const caf::SRProxy *sr){
    double cvn = sr->sel.cvnProd3Train.nueid;
    double bdt = sr->sel.nuecosrej.cospidcontain;
    if (cvn > 0.99 || (cvn>0.95 && bdt > 0.57)) return true;
    return false;
  });

const Cut kNue2018FDPeripheral = kNue2017PeripheralPresel &&
                                 kNue2018CVNVsCosPID;

const Var kNue2018SelectionBin([](const caf::SRProxy *sr){
    bool isPeri = kNue2018FDPeripheral(sr);
    double cvn = sr->sel.cvnProd3Train.nueid;
    if (isPeri)
      return float(3.5);
    if (cvn > 0.75 && cvn < 0.87)
      return float(0.5);
    if (cvn > 0.87 && cvn < 0.95)
      return float(1.5);
    if (cvn > 0.95)
      return float(2.5);
    return float(-0.5);
  });


const Var kNue2018AnaBin([](const caf::SRProxy *sr){
    int selBin = kNue2018SelectionBin(sr);

    float nuE  = kNueEnergy2017(sr);
    int nuEBin = nuE/0.5;
    assert(nuEBin <= 8 && "An event with nuE > 4.5 should never happen");

    int anaBin = 9*selBin + nuEBin;

    if (anaBin > 36)
      anaBin = -1;

    return anaBin;
  });


const Var kCVNSSe([](const caf::SRProxy* sr){
    return sr->sel.cvnProd3Train.nueid;
  });

const Cut kNue2018FDCore = kCVNSSe > 0.75 && kNue2017CorePresel;
const Cut kNue2018FD = kNue2018FDCore || kNue2018FDPeripheral;



// Macro for optimizing PID cuts / binning in the FD.  Just works on RHC now
// TODO: incorperate FHC.  Update presel cuts if we decide to change them
void make_nue_pidtuning(bool makePlots=false)
{

  if (!makePlots){

    SpectrumLoader lSwap("dpershey_prod4-cosrej_rhc_mc_fluxswap");
    SpectrumLoader lNonS("dpershey_prod4-cosrej_rhc_mc_nonswap");
    SpectrumLoader lTauS("dpershey_prod4-cosrej_rhc_mc_tauswap");

    SpectrumLoader lCosmic("dpershey_prod4-cosrej_cosmicdata_rhctune");

    lSwap.SetSpillCut(kStandardSpillCuts);
    lNonS.SetSpillCut(kStandardSpillCuts);
    lTauS.SetSpillCut(kStandardSpillCuts);
    lCosmic.SetSpillCut(kStandardSpillCuts);

    Cut selCore = kNue2017CorePresel;
    Cut selPeri = kNue2017PeripheralPresel;
    Var wei = kXSecCVWgt2017*kPPFXFluxCVWgt;



    // Peripheral cuts are 2D PID cuts, set up 2D vars now
    Var kBDTLight = SIMPLEVAR(sel.nuecosrej.cospidlight);
    Var kBDTCont  = SIMPLEVAR(sel.nuecosrej.cospidcontain);

    Binning bdtbins = Binning::Simple(50,0.35,0.85);
    Binning cvnbins = Binning::Simple(50,0.5,1.);

    Var kLightCVN = Var2D(kBDTLight,bdtbins,
                          kCVNSSe,cvnbins);
    Var kContCVN  = Var2D(kBDTCont,bdtbins,
                          kCVNSSe,cvnbins);


    // Prediction as a function of CVN, for main PID tuning
    // Be sure you're using the right core preselection, it may change
    PredictionNoExtrap predCVN(lNonS,lSwap,lTauS,"",
                               Binning::Simple(50,0.5,1),
                               kCVNSSe, selCore, kNoShift, wei);
    Spectrum cosCVN(lCosmic,HistAxis("",Binning::Simple(50,0.5,1),kCVNSSe),
                    kInCosmicTimingWindow&&selCore,kNoShift,wei);

    // 2D predictions for determining peripheral PID cuts
    PredictionNoExtrap predLightPeri(lNonS,lSwap,lTauS,
                                     "",Binning::Simple(2500,0,2500),kLightCVN,
                                     selPeri,kNoShift,wei); // BDTLight
    Spectrum cosLightPeri(lCosmic,
                          HistAxis("",Binning::Simple(2500,0,2500),kLightCVN),
                          selPeri&&kInCosmicTimingWindow,kNoShift,wei);
    PredictionNoExtrap predContPeri (lNonS,lSwap,lTauS,
                                     "",Binning::Simple(2500,0,2500),kContCVN,
                                     selPeri,kNoShift,wei); // BDTContain
    Spectrum cosContPeri (lCosmic,
                          HistAxis("",Binning::Simple(2500,0,2500),kContCVN),
                          selPeri&&kInCosmicTimingWindow,kNoShift,wei);

    // Prediction with 2017 bin boundaries, just for fun
    PredictionNoExtrap pred2017(lNonS,lSwap,lTauS,"",
                                Binning::Simple(36,0,36),
                                kNue2018AnaBin,
                                kNue2018FD,kNoShift,wei);
    Spectrum cos2017(lCosmic,
                     HistAxis("",Binning::Simple(36,0,36),kNue2018AnaBin),
                     kInCosmicTimingWindow&&kNue2018FD,kNoShift,wei);









    lSwap.Go();
    lNonS.Go();
    lTauS.Go();
    lCosmic.Go();

    TFile* outFile = new TFile("out_fdcosrej.root","RECREATE");
    outFile->cd();

    pred2017.SaveTo(outFile, "pred2017");
    cos2017.SaveTo(outFile, "cos2017");

    predCVN.SaveTo(outFile, "predCVN");
    cosCVN.SaveTo(outFile, "cosCVN");

    predLightPeri.SaveTo(outFile, "predLightPeri");
    cosLightPeri .SaveTo(outFile, "cosLightPeri");

    predContPeri.SaveTo(outFile, "predContPeri");
    cosContPeri .SaveTo(outFile, "cosContPeri");


    outFile->Close();
  } // End if !makePlots
  else {

    PeripheralCuts("LightPeri");
    CVNCuts("CVN");

  }
}