fiducial_accounting.C

Go to the documentation of this file.

#include <bitset>
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/SpectrumLoaderBase.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Cut.h"
#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "3FlavorAna/Cuts/NueCuts2017.h"
#include "CAFAna/Decomp/ProportionalDecomp.h"
#include "3FlavorAna/Decomp/NueDecomp.h"
#include "CAFAna/Decomp/NCDecomp.h"
#include "3FlavorAna/Decomp/NumuDecomp.h"
#include "CAFAna/Extrap/ModularExtrap.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "TFile.h"
#include "OscLib/OscCalcDumb.h"
#include "CAFAna/Core/Spectrum.h"
#include "TF1.h"
#include "TH1.h"
#include "TH2.h"
#include "TEfficiency.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "TCanvas.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "CAFAna/Vars/FitVars.h"
#include "CAFAna/Core/Utilities.h"
#include "Utilities/func/MathUtil.h"
#include "CAFAna/Cuts/AnalysisMasks.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;
                     }); 


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::kNuMuToNuE,
                               Current::kCC,
                               Sign::kBoth).ToTH1(kPOTFD );
  hSig->SetLineColor(kRed-10);

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

  TH1D* hN = new TH1D("hN_"+title,"Nue CC selection efficiency;diblocks;efficiency",11,3.5,14.5);
  TH1D* hD = new TH1D("hD_"+title,"Nue CC 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();

}

void fiducial_accounting()
{
  
  TH1::SetDefaultSumw2();
  TH2::SetDefaultSumw2();

  const Cut kSelectionCut = kNue2017CorePresel;

  // SpectrumLoader loaderFDSwap(
  //   "prod_caf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_v1" );
  //
  // NullLoader loaderFDNonSwap;
  Loaders loader_fdmc;
  loader_fdmc.SetLoaderPath("prod_caf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_v1", 
                            caf::kFARDET, Loaders::kMC, kBeam, Loaders::kFluxSwap);

  const HistAxis axis(
    "diblocks", Binning::Simple(15,-0.5,14.5), kBlock );

  // make prediction for all true neutrinos in fiducial volume
  PredictionNoExtrap* predAll = new PredictionNoExtrap(
    loader_fdmc,
    axis, kBoxCut
  );
  //prediction for all true neutrinos in fiducial volume that pass preselection
  PredictionNoExtrap* predSel = new PredictionNoExtrap(
    loader_fdmc,
    axis, kSelectionCut&&kBoxCut
  );
  // make prediction for all true neutrinos in small box at center of detector
  PredictionNoExtrap* predIdealAll = new PredictionNoExtrap(
    loader_fdmc,
    axis, kBoxCutIdeal
  );
  //prediction for all true neutrinos in fiducial volume that pass preselection in small box at center of detector
  PredictionNoExtrap* predIdealSel = new PredictionNoExtrap(
    loader_fdmc,
    axis, kSelectionCut&&kBoxCutIdeal
  );
    
  loader_fdmc.SetSpillCut(kStandardSpillCuts);



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

  // Open output file
  TFile file("fd_fiducial_mass_efficiency.root","RECREATE");
  // Save things in output file

  
  
  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();

}