print_nus17_fd_cut_tables2D.C

Go to the documentation of this file.

//  This macro is used to make the nue energy blessed plots. (see doc 11789)
//  Oscillation parameters are same as those in nue_ana_basic macro

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/NueCutsFirstAna.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/Style.h"
#include "3FlavorAna/Plotting/NuePlotStyle.h"
#include "CAFAna/Core/LoadFromFile.h"

#include "OscLib/OscCalcDumb.h"
#include "OscLib/OscCalcSterile.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "NuXAna/Prediction/PredictionSterile.h"

#include "NuXAna/macros/Nus17/nus17_fd_arrays.h"

#include "TH1D.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TColor.h"
#include "TGraph.h"
#include "THStack.h"
#include "TMultiGraph.h"
#include "TLine.h"
#include "TLegend.h"
#include "TString.h"
#include "Utilities/rootlogon.C"

using namespace ana;


const double pot_val = 9.0e20; // POT
const double livetime = 440.; // seconds


void GetSpectra(IPrediction* pred, osc::IOscCalc* calc,
                TH1*& hSig, TH1*& hNue, TH1*& hCC, TH1*& hTau, TH1*& hTotBkg)
{
  hSig = pred->PredictComponent(calc,
                                Flavors::kAll,
                                Current::kNC,
                                Sign::kBoth).ToTH1(pot_val);
  hSig->SetLineColor(kNCBackgroundColor);


  hNue = pred->PredictComponent(calc,
                                Flavors::kAllNuE,
                                Current::kCC,
                                Sign::kBoth).ToTH1(pot_val);
  hNue->SetLineColor(kNueSignalColor);
  
  hCC = pred->PredictComponent(calc,
                               Flavors::kAllNuMu,
                               Current::kCC,
                               Sign::kBoth).ToTH1(pot_val);
  hCC->SetLineColor(kNumuBackgroundColor);
  
  hTau = pred->PredictComponent(calc,
                                Flavors::kAllNuTau,
                                Current::kCC,
                                Sign::kBoth).ToTH1(pot_val);
  hTau->SetLineColor(kMagenta);
  
  hTotBkg = (TH1*)hNue->Clone(UniqueName().c_str());
  hTotBkg->Add(hCC);
  hTotBkg->Add(hTau);

  hTotBkg->SetLineColor(kTotalMCColor);
}
void GetSpectra2D(IPrediction* pred, osc::IOscCalc* calc,
                  TH2*& hSig, TH2*&hTotBkg)
{
  hSig = pred->PredictComponent(calc,
                                Flavors::kAll,
                                Current::kNC,
                                Sign::kBoth).ToTH2(pot_val);

  TH1* hNue = pred->PredictComponent(calc,
                                     Flavors::kNuMuToNuE,
                                     Current::kCC,
                                    Sign::kBoth).ToTH2(pot_val);
  
  TH1* hCC = pred->PredictComponent(calc,
                                    Flavors::kAllNuMu,
                                    Current::kCC,
                                    Sign::kBoth).ToTH2(pot_val);

  TH1* hBeam = pred->PredictComponent(calc,
                                      Flavors::kNuEToNuE,
                                      Current::kCC,
                                      Sign::kBoth).ToTH2(pot_val);

  hTotBkg = (TH2*)hNue->Clone(UniqueName().c_str());
  hTotBkg->Add(hCC);
  hTotBkg->Add(hBeam);

}


//------------------------------------------------------------

void print_tables(const std::string* cutnames,
                  const std::string pred_basename,
                  const std::string spec_basename,
                  const std::string fname,
                  const double cos_scaling,
                  const int ncuts,
                  bool savePlots=false)
{
  osc::OscCalcSterile* osc = DefaultSterileCalc(3);


  
  for(int ic = 7; ic < ncuts_nminus1; ic++){
    
    char name[50];

    sprintf(name, "%s_%s_%s", 
            pred_basename.c_str(), 
            cutnames_nminus1[ic].c_str(), 
            "2Dtime");
    
    std::cout << "name: "<< name << std::endl;
    IPrediction* pred = 
      LoadFromFile<IPrediction>(fname, name).release();
    
    sprintf(name, "%s_%s_%s", 
            spec_basename.c_str(),
            cutnames_nminus1[ic].c_str(), 
            "2Dtime");
    std::cout << "name: "<< name << std::endl;
    Spectrum* spec = 
      LoadFromFile<Spectrum>(fname, name).release();
    
    TH2 *sighist, *bghist;
    GetSpectra2D(pred, osc, sighist, bghist);
    std::cout << "here: "<< std::endl;
    TH2* coshist = spec->ToTH2(cos_scaling, 
                               ana::EExposureType::kLivetime);
    
    new TCanvas;

    //coshist->RebinX(3);
    //coshist->RebinY(2);
    ///coshist->SetTitle("Cosmic Data");
    //    coshist->GetXaxis()->SetTitle("Shower Maximum Y coordinate (cm)");
    //coshist->GetYaxis()->SetTitle("p_{T}/p");
    coshist->GetXaxis()->CenterTitle();
    coshist->GetYaxis()->CenterTitle();

    coshist->Draw("colz");
    Preliminary();
    gPad->SaveAs(("plots2D/time_"+spec_basename+cutnames_nminus1[ic]+"2Dcos.png").c_str());
    gPad->SaveAs(("plots2D/time_"+spec_basename+cutnames_nminus1[ic]+"2Dcos.eps").c_str());
    gPad->SaveAs(("plots2D/time_"+spec_basename+cutnames_nminus1[ic]+"2Dcos.pdf").c_str());
    gPad->SaveAs(("plots2D/time_"+spec_basename+cutnames_nminus1[ic]+"2Dcos.root").c_str());
    
    new TCanvas;
    ///gStyle->SetTitleAlign(11);
    //sighist->RebinX(3);
    //sighist->RebinY(2);
    //sighist->SetTitle("NC Signal");
    //sighist->GetXaxis()->SetTitle("Shower Maximum Y coordinate (cm)");
    //sighist->GetYaxis()->SetTitle("p_{T}/p");
    sighist->GetXaxis()->CenterTitle();
    sighist->GetYaxis()->CenterTitle();
    sighist->Draw("colz");
    Preliminary();
    gPad->SaveAs(("plots2D/time_"+spec_basename+cutnames_nminus1[ic]+"2Dsig.png").c_str());
    gPad->SaveAs(("plots2D/time_"+spec_basename+cutnames_nminus1[ic]+"2Dsig.eps").c_str());
    gPad->SaveAs(("plots2D/time_"+spec_basename+cutnames_nminus1[ic]+"2Dsig.pdf").c_str());
    gPad->SaveAs(("plots2D/time_"+spec_basename+cutnames_nminus1[ic]+"2Dsig.root").c_str());
    
    
    new TCanvas;
    bghist->Draw("colz");
    gPad->SaveAs(("plots2D/time_"+spec_basename+cutnames_nminus1[ic]+"2Dbg.png").c_str());
    gPad->SaveAs(("plots2D/time_"+spec_basename+cutnames_nminus1[ic]+"2Dbg.eps").c_str());
    gPad->SaveAs(("plots2D/time_"+spec_basename+cutnames_nminus1[ic]+"2Dbg.pdf").c_str());
    gPad->SaveAs(("plots2D/time_"+spec_basename+cutnames_nminus1[ic]+"2Dbg.root").c_str());
    
  }

}

//------------------------------------------------------------

void print_nus17_fd_cut_tables2D(bool savePlots=false)
{
  
  //std::string fname = "/nova/ana/users/gsdavies/nus/nus17/FD_cutTable_extendedE.root";
  std::string fname = "/nova/ana/users/gsdavies/nus/nus17/mcv1_2D_cut_tables_v2.root";
  //std::string fname = "mcv1_cutTable.root";

  // cosmic livetime scaling
  TFile fin(fname.c_str());
  if(fin.IsZombie()){
    std::cerr << "File: " << fname << " does not exist. Check the name." << std::endl;
    exit(0);
  }
    
    
  TH1D* time = (TH1D*)gDirectory->Get("time_cos");
  double scaling = 440.;

  char name[15];
  
  
  std::cout<<"\n--------------------------------------"
           <<"\n---------------N-1 Table -------------\n";
  
  print_tables(cutnames_nminus1,        
               "pred_nminus1",
               "spec_nminus1",
               fname, scaling, ncuts_nminus1, savePlots);
  
  std::cout<<"\n------------------------------------------"
           <<"\n--------------Cut Flow Table -------------\n";
  
  // print_tables(cutnames,
  //           "pred",
  //           "spec",
  //           fname, scaling, ncuts, savePlots);

}