nus17_fiducial_accounting.C

Go to the documentation of this file.

#include <bitset>
#include <string>
#include <vector>

#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/SpectrumLoaderBase.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "NuXAna/Cuts/NusCuts17.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "OscLib/OscCalcDumb.h"

#include "TFile.h"
#include "TF1.h"
#include "TH1.h"
#include "TH2.h"
#include "TEfficiency.h"
#include "TCanvas.h"


using namespace ana;

const Cut kBoxCutIdeal(
  [](const caf::SRProxy* sr)
  {
    if(sr->mc.nnu == 0) return false;
    std::bitset<14> binary(sr->hdr.dibmask);
    std::pair<int,int> planesA = calcFirstLastLivePlane(sr->slc.firstplane,
                                                        binary);
    std::pair<int,int> planesB = calcFirstLastLivePlane(sr->slc.lastplane,
                                                        binary);
    if((planesA.first != planesB.first) || (planesA.second != planesB.second))
      return false;
    if((planesA.second - planesA.first + 1)/64 < 14) return false;
    if(sr->mc.nu[0].vtx.X() < -100) return false;
    if(sr->mc.nu[0].vtx.X() >  100) return false;
    if(sr->mc.nu[0].vtx.Y() < -100) return false;
    if(sr->mc.nu[0].vtx.Y() >  100) return false;
    if(sr->mc.nu[0].vtx.Z() < 1000) return false;
    if(sr->mc.nu[0].vtx.Z() > 3500) return false;
    return true;
  }); 


const Cut kBoxCut(
  [](const caf::SRProxy* sr)
  {
    if(sr->mc.nnu == 0) return false;
    if(sr->mc.nu[0].vtx.X() < -762) return false;
    if(sr->mc.nu[0].vtx.X() >  763) return false;
    if(sr->mc.nu[0].vtx.Y() < -762) return false;
    if(sr->mc.nu[0].vtx.Y() >  763) return false;
    std::bitset<14> binary(sr->hdr.dibmask);
    std::pair<int,int> planesA = calcFirstLastLivePlane(sr->slc.firstplane,
                                                        binary);
    std::pair<int,int> planesB = calcFirstLastLivePlane(sr->slc.lastplane,
                                                        binary);
    if((planesA.first != planesB.first) || (planesA.second != planesB.second))
      return false;
    if((planesA.second - planesA.first + 1)/64 < 4) return false;
    int first = planesA.first/64;
    int last  = (planesA.second +1)/64;
    if(sr->mc.nu[0].vtx.Z() < ((first)*426.0)) return false;
    if(sr->mc.nu[0].vtx.Z() > ((last)*426.0))  return false;
    return true;
  });
                     
const Var kBlock(
  [](const caf::SRProxy* sr)
  {
    int size = -1;
    std::bitset<14> binary(sr->hdr.dibmask);
    std::pair<int,int> planesA = calcFirstLastLivePlane(sr->slc.firstplane,
                                                        binary);
    std::pair<int,int> planesB = calcFirstLastLivePlane(sr->slc.lastplane,
                                                        binary);
    if((planesA.first != planesB.first) || (planesA.second != planesB.second))
      return size;
    size = (planesA.second - planesA.first + 1)/64;
    return size;
  });

void GetExtrap(PredictionNoExtrap* all,
               PredictionNoExtrap* sel,
               osc::IOscCalc* calc,
               TString title)
{
  TH1::SetDefaultSumw2();
  TH2::SetDefaultSumw2();
  double kPOTFD = all->Predict(calc).POT();
  
  TH1* hSig = all->PredictComponent(calc,
                                    Flavors::kAll,
                                    Current::kNC,
                                    Sign::kBoth).ToTH1(kPOTFD);
  hSig->SetLineColor(kRed-10);

  TH1* hSigSel = sel->PredictComponent(calc,
                                       Flavors::kAll,
                                       Current::kNC,
                                       Sign::kBoth).ToTH1(kPOTFD);
                               

  TH1D* hSigEff = new TH1D("hSigEff_"+title,
                           "Nus NC Selection Efficiency;Diblocks;Efficiency", 11, 3.5, 14.5);
  TH1D* hN = new TH1D("hN_"+title,
                      "Nus NC Selection Efficiency;Diblocks;Efficiency", 11, 3.5, 14.5);
  TH1D* hD = new TH1D("hD_"+title,
                      "Nus NC Selection Efficiency;Diblocks;Efficiency", 11, 3.5, 14.5);
  for(unsigned int i = 5; i < 16; ++i){
    for(unsigned int j = 0; j < hSigSel->GetBinContent(i); ++j){
      hN->Fill(i-1);
    }
    for(unsigned int j = 0; j < hSig->GetBinContent(i); ++j){
      hD->Fill(i-1);
    }
  }
  hN->Write();
  hD->Write();
  hSigEff->Divide(hN, hD, 1.0, 1.0, "b");
  
  std::cout<<"Calculating efficiency for: "<<title<<std::endl;
  for(unsigned int i = 1; i < 12; ++i){
    float sigN = hN->GetBinContent(i);
    float sigD = hD->GetBinContent(i);
    float sigE = hSigEff->GetBinContent(i);
    std::cout << "Signal Eff, DB: " << i+3
              << " Selected: " << sigN
              << " All: " << sigD
              << " Efficiency: " << sigE << std::endl;
  }
  hSigEff->Write();
}

void nus17_fiducial_accounting(std::string outfile)
{
  
  TH1::SetDefaultSumw2();
  TH2::SetDefaultSumw2();

  SpectrumLoader loaderFDSwap("prod_caf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_v1");

  SpectrumLoader loaderFDNonSwap("prod_caf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_v1");

  const Cut kNusEff = kNus17FDPresel && kNHit >= 25;
  
  const HistAxis axis("Diblocks", Binning::Simple(15, -0.5, 14.5), kBlock);

  // Prediction for all true neutrinos in fiducial volume
  PredictionNoExtrap* predAll = new PredictionNoExtrap(
    loaderFDNonSwap, loaderFDSwap,
    axis.GetLabels()[0], axis.GetBinnings()[0], axis.GetVars()[0], kBoxCut
  );
  
  // Prediction for all true neutrinos in fiducial volume that pass preselection
  PredictionNoExtrap* predSel = new PredictionNoExtrap(
    loaderFDNonSwap, loaderFDSwap,
    axis.GetLabels()[0], axis.GetBinnings()[0], axis.GetVars()[0], kNusEff && kBoxCut
  );
  
  // Pprediction for all true neutrinos in small box at center of detector
  PredictionNoExtrap* predIdealAll = new PredictionNoExtrap(
    loaderFDNonSwap, loaderFDSwap,
    axis.GetLabels()[0], axis.GetBinnings()[0], axis.GetVars()[0], kBoxCutIdeal
  );
  // Prediction for all true neutrinos in fiducial volume that pass preselection
  // in small box at center of detector
  PredictionNoExtrap* predIdealSel = new PredictionNoExtrap(
    loaderFDNonSwap, loaderFDSwap,
    axis.GetLabels()[0], axis.GetBinnings()[0], axis.GetVars()[0], kNusEff && kBoxCutIdeal
  );
    
  loaderFDSwap.SetSpillCut(kStandardSpillCuts);
  loaderFDNonSwap.SetSpillCut(kStandardSpillCuts);

  // Load from files
  loaderFDNonSwap.Go();
  loaderFDSwap.Go();

  // Open output file
  TFile file(outfile.c_str(),"RECREATE");

  osc::OscCalcDumb calc;
  
  // Caclulate preselection efficieny as a function of diblock
  GetExtrap(predAll, predSel, &calc, "FidMass");
  // Caclulate preselection efficieny for small window in ideal detector
  GetExtrap(predIdealAll, predIdealSel, &calc, "SmallWindow");

  file.Close();
}