Ana2017_sb_opening_macro.C

Go to the documentation of this file.

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

#include "CAFAna/Core/EventList.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"

#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Vars/HistAxes.h"
#include "CAFAna/Vars/GenieWeights.h"

#include "OscLib/OscCalcDumb.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"

#include "TH1D.h"

using namespace ana;

void Ana2017_sb_opening_macro()
{
  //Make time axis for 2017 analysis period
  const Var kSliceTime = SIMPLEVAR(hdr.unixtime);
  const Var kSliceMeanTime = SIMPLEVAR(hdr.subevtmeantime);
  const HistAxis kTimeAxis("unix time (s)",Binning::Simple(10000, 1391000000, 1491000000), kSliceTime);
  const HistAxis kSliceTimeAxis("slice time (ns)",Binning::Simple(600, -50000, 550000), kSliceMeanTime);
  
  const HistAxis kNue2017SBEnergyAxis( "Nue Energy Estimator (GeV)",
                                       Binning::Simple( 20, 0, 10 ), kNueEnergy2017 );
  
  const std::string fnamecos  =
    //"prod_decaf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1";
    "prod_restricteddecaf_R17-03-01-prod3reco.k_fd_numi_fhc_full_v1_goodruns_nue2017";
  
  SpectrumLoader loadercos(fnamecos);
  
  loadercos.SetSpillCut(kStandardSpillCuts);
  
  struct CutDef{
    const Cut cut;
    std::string label;
  };
  
  const int kNCuts = 4;
  const CutDef SBCuts[kNCuts] = 
    {{kNue2017FDPeripheralHELPsideband,"HELP"}, // Peripheral High Energy Low PID
     {kNue2017FDPeripheralHEMPsideband,"HEMP"}, // Peripheral High Energy Mid PID
     {kNue2017FDPeripheralLPsideband,"LPID"},   // Peripheral Low PID
     {kNue2017FDPeripheralMPsideband,"MPID"}    // Peripheral Mid PID
    }; 
  
  Spectrum* spec_cvn_cvn2d_outtime[kNCuts];
  Spectrum* spec_cvn_nueE_outtime[kNCuts];
  Spectrum* spec_cvn_time_outtime[kNCuts];
  Spectrum* spec_cvn_slctime_outtime[kNCuts];
  
  Spectrum* spec_cvn_cvn2d_intime[kNCuts];
  Spectrum* spec_cvn_nueE_intime[kNCuts];
  Spectrum* spec_cvn_time_intime[kNCuts];
  Spectrum* spec_cvn_slctime_intime[kNCuts];
  
  for(unsigned int i = 0 ; i < kNCuts ; ++i){
    // out of time spectra
    spec_cvn_cvn2d_outtime[i] = new Spectrum(loadercos, kNue2017Axis,
                                             SBCuts[i].cut && kInTimingSideband,
                                             kNoShift);
    
    spec_cvn_nueE_outtime[i] = new Spectrum(loadercos, kNue2017SBEnergyAxis,
                                            SBCuts[i].cut && kInTimingSideband,
                                            kNoShift);
    
    spec_cvn_time_outtime[i] = new Spectrum(loadercos, kTimeAxis,
                                            SBCuts[i].cut && kInTimingSideband,
                                            kNoShift);
    
    spec_cvn_slctime_outtime[i] = new Spectrum(loadercos, kSliceTimeAxis,
                                               SBCuts[i].cut && kInTimingSideband,
                                               kNoShift);
    
    // in time spectra
    spec_cvn_cvn2d_intime[i] = new Spectrum(loadercos, kNue2017Axis,
                                            SBCuts[i].cut && kInBeamSpill,
                                            kNoShift);
    
    spec_cvn_nueE_intime[i] = new Spectrum(loadercos, kNue2017SBEnergyAxis,
                                           SBCuts[i].cut && kInBeamSpill,
                                           kNoShift);
    
    spec_cvn_time_intime[i] = new Spectrum(loadercos, kTimeAxis,
                                           SBCuts[i].cut && kInBeamSpill,
                                           kNoShift);
    
    spec_cvn_slctime_intime[i] = new Spectrum(loadercos, kSliceTimeAxis,
                                              SBCuts[i].cut && kInBeamSpill,
                                              kNoShift);
  }
  
  // let's go!
  loadercos.Go();
  
  std::cout<<"Making output file ";
  TFile fout("Ana2017_numi_sb.root","RECREATE");
  std::cout<< fout.GetName() << std::endl;
  
  for(unsigned int i = 0 ; i < kNCuts ; ++i){
    spec_cvn_cvn2d_intime[i]->SaveTo(&fout, ("spec_cvn_cvn2d_numi_"+SBCuts[i].label).c_str());    
    spec_cvn_cvn2d_outtime[i]->SaveTo(&fout, ("spec_cvn_cvn2d_cosmic_"+SBCuts[i].label).c_str());    
  
    spec_cvn_nueE_intime[i]->SaveTo(&fout, ("spec_cvn_nueE_numi_"+SBCuts[i].label).c_str());    
    spec_cvn_nueE_outtime[i]->SaveTo(&fout, ("spec_cvn_nueE_cosmic_"+SBCuts[i].label).c_str());    
    
    spec_cvn_time_intime[i]->SaveTo(&fout, ("spec_cvn_time_numi_"+SBCuts[i].label).c_str());    
    spec_cvn_time_outtime[i]->SaveTo(&fout, ("spec_cvn_time_cosmic_"+SBCuts[i].label).c_str());    
    
    spec_cvn_slctime_intime[i]->SaveTo(&fout, ("spec_cvn_slctime_numi_"+SBCuts[i].label).c_str());    
    spec_cvn_slctime_outtime[i]->SaveTo(&fout, ("spec_cvn_slctime_cosmic_"+SBCuts[i].label).c_str());    
  }
  
  Var kRun = SIMPLEVAR(hdr.run);
  Var kSubrun = SIMPLEVAR(hdr.subrun);
  Var kEvt = SIMPLEVAR(hdr.evt);
  Var kSlice = SIMPLEVAR(hdr.subevt);
  
  std::vector<const Var*> vars = {&kRun, &kSubrun, &kEvt, &kSlice, 
                                  &kSliceMeanTime,
                                  &kCVNe, 
                                  &kNueEnergy2017
  };
  
  std::vector<Cut> cuts;
  std::vector<std:: string> fnames;
  for(unsigned int i = 0 ; i < kNCuts ; ++i){
    cuts.push_back(SBCuts[i].cut && kInBeamSpill);
    cuts.push_back(SBCuts[i].cut && kInTimingSideband);
    
    fnames.push_back(("Ana2017_selected_numi_sb_"+SBCuts[i].label+".txt").c_str());
    fnames.push_back(("Ana2017_selected_cosmic_sb_"+SBCuts[i].label+".txt").c_str());
  }
  
  MakeTextListFile(fnamecos,
                   cuts,
                   fnames,
                   vars,
                   &kStandardSpillCuts );

}