signal_count.C

Go to the documentation of this file.

#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Vars/XsecTunes.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/XSec/Multiverse.h"

#include "CAFAna/Analysis/Exposures.h"

#include "TCanvas.h"
#include "TStyle.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "NDAna/nuebarcc_inc/NuebarCCIncCuts.h"
#include "NDAna/nuebarcc_inc/NuebarCCIncVars.h"
#include "NDAna/nuebarcc_inc/NuebarCCIncBins.h"
#include "NDAna/nuebarcc_inc/NuebarCCIncExtra.h"
// temporary measure to use prod4 trained muonid
#include "NDAna/muonid/NDXSecMuonPID.h"
namespace nuebarccinc {
  const ana::Cut kMuonIDProd4Cut = ana::kMuonIDProd4 < -0.2;
}

void load_libs_muonid()
{
  const int ret = gSystem->Load("libNDAnamuonid");
  if(ret != 0) exit(1);
}

using namespace ana;

void signal_count(bool make_plots=false, 
                  std::string input_file_name="/nova/ana/users/ddoyle/NuebarResolutionStudy/signal_count.root",
                  std::string plot_dump="/nova/ana/users/ddoyle/NuebarResolutionStudy/signal_count_plots/") 
{
  if(!make_plots) {
    load_libs_muonid();

    SpectrumLoader lNDMC(nuebarccinc::PROD5_MC_RHC_NOMINAL);
    lNDMC.SetSpillCut(kStandardSpillCuts);

    const Var kCVWeight = kPPFXFluxCVWgt * kXSecCVWgt2020;

    const Cut kSelection = 
      nuebarccinc::kPreselectionLoose &&
      nuebarccinc::kMuonIDProd4Cut;
    //    nuebarccinc::kRemid;

    Spectrum nue(lNDMC,
                 nuebarccinc::kTrueElectronEVsCosStandardAxis,
                 kSelection && nuebarccinc::kNueCC,
                 kNoShift,
                 kCVWeight);


    Spectrum nuebar(lNDMC,
                    nuebarccinc::kTrueElectronEVsCosStandardAxis,
                    kSelection && nuebarccinc::kNuebarCC,
                    kNoShift,
                    kCVWeight);

    lNDMC.Go();

    TFile output("signal_count.root", "recreate");
    nue.SaveTo(output.mkdir("nue"));
    nuebar.SaveTo(output.mkdir("nuebar"));
    output.Close();
  }
  else {
    TFile * input = TFile::Open(input_file_name.c_str());
    TCanvas * c = new TCanvas();

    Spectrum * nue1D    = Spectrum::LoadFrom(input, "nue"   )).release(;
    Spectrum * nuebar1D = Spectrum::LoadFrom(input, "nuebar")).release(;
    Spectrum * nue2D   = new Spectrum(nue1D->ToTH1(kAna2020RHCPOT),
                                      {nuebarccinc::kTrueElectronCosThetaStandardAxis.GetLabels()[0],
                                          nuebarccinc::kTrueElectronEStandardAxis.GetLabels()[0]},
                                      {nuebarccinc::kTrueElectronCosThetaStandardAxis.GetBinnings()[0],
                                          nuebarccinc::kTrueElectronEStandardAxis.GetBinnings()[0]},
                                      kAna2020RHCPOT, 0);
    Spectrum * nuebar2D = new Spectrum(nuebar1D->ToTH1(kAna2020RHCPOT),
                                       {nuebarccinc::kTrueElectronCosThetaStandardAxis.GetLabels()[0],
                                           nuebarccinc::kTrueElectronEStandardAxis.GetLabels()[0]},
                                       {nuebarccinc::kTrueElectronCosThetaStandardAxis.GetBinnings()[0],
                                           nuebarccinc::kTrueElectronEStandardAxis.GetBinnings()[0]},
                                       kAna2020RHCPOT, 0);

    // count plots
    gStyle->SetPalette(kDeepSea);
    gStyle->SetPaintTextFormat(".2f");
    TH2 * nue    = nue2D   ->ToTH2(kAna2020RHCPOT);
    nue->SetTitle("Signal #nu_{e}");
    nue->SetMarkerColor(kWhite);
    nue->GetXaxis()->SetRangeUser(0.5, 1);
    nue->GetYaxis()->SetRangeUser(0, 6);
    nue->Draw("colz text");
    c->Print((plot_dump + "/nue.pdf").c_str());
    TH2 * nuebar = nuebar2D->ToTH2(kAna2020RHCPOT);
    nuebar->SetTitle("Signal #bar{#nu}_{e}");
    nuebar->SetMarkerColor(kWhite);
    nuebar->GetXaxis()->SetRangeUser(0.5, 1);
    nuebar->GetYaxis()->SetRangeUser(0, 6);
    nuebar->Draw("colz text");
    c->Print((plot_dump + "/nuebar.pdf").c_str());

    TH2 * nunubar = (TH2*) nue->Clone("nunubar");
    nunubar->Add(nuebar);
    nunubar->SetTitle("Signal #bar{#nu}_{e} + #nu_{e}");
    nunubar->SetMarkerColor(kWhite);
    nunubar->GetXaxis()->SetRangeUser(0.5, 1);
    nunubar->GetYaxis()->SetRangeUser(0, 6);
    nunubar->Draw("colz text");
    c->Print((plot_dump + "/nunubar.pdf").c_str());


    // fraction plots
    auto kNuebarColor = kGreen;
    auto kNueColor = kMagenta;

    gStyle->SetPaintTextFormat("4.2g");
    TH2 * nue_frac    = (TH2*) nue   ->Clone("nue_frac"   );
    TH2 * nuebar_frac = (TH2*) nuebar->Clone("nuebar_frac");
    nue_frac->SetMarkerColor(kNueColor);
    nuebar_frac->SetMarkerColor(kNuebarColor);

    TLegend * leg = new TLegend(0.2, 0.8, 0.3, 0.7, "", "NDC");

    auto entry = leg->AddEntry((TObject*)0, "#bar{#nu}_{e} Fraction", "");
    entry->SetTextColor(kNuebarColor);
    entry = leg->AddEntry((TObject*)0, "#nu_{e} Fraction"      , "");
    entry->SetTextColor(kNueColor);
    
    nue_frac->Divide(nunubar);
    nuebar_frac->Divide(nunubar);
    
    nunubar    ->Draw("colz");
    // nuebar on top
    nuebar_frac->SetBarOffset(0.2);
    nue_frac   ->SetBarOffset(-0.2);
    nue_frac   ->Draw("text same");
    nuebar_frac->Draw("text same");
    leg->Draw();
    c->Print((plot_dump + "/nunubar_fracs.pdf").c_str());
    


    // uncertainty plots
    TH2 * nuebar_uncert = (TH2*) nuebar->Clone("nuebar_uncert");
    TH2 * nunubar_uncert = (TH2*) nunubar->Clone("nunubar_uncert");
    auto NX = nuebar_uncert->GetNbinsX();
    auto NY = nuebar_uncert->GetNbinsY();
    for(auto x = 1; x <= NX; x++) {
      for(auto y = 1; y <= NY; y++) {
        auto nnuebar = nuebar->GetBinContent(x, y);
        auto nnunubar = nunubar->GetBinContent(x, y);
        
        if(nnuebar) 
          nuebar_uncert->SetBinContent(x, y, 1/std::sqrt(nnuebar));
        else
          nuebar_uncert->SetBinContent(x, y, 0);

        if(nnunubar) 
          nunubar_uncert->SetBinContent(x, y, 1/std::sqrt(nnunubar));
        else
          nunubar_uncert->SetBinContent(x, y, 0);
      }
    }

    TH2 * uncert_diff = (TH2*) nunubar_uncert->Clone("uncert_diff");
    uncert_diff->Add(nuebar_uncert, -1);
    
    nunubar_uncert->SetTitle("Signal #nu_{e} + #bar{#nu}_{e} Stat. Uncert.");
    nuebar_uncert->SetTitle("Signal #bar{#nu}_{e} Stat. Uncert.");
    uncert_diff->SetTitle("#Delta Stat. Uncert. (#bar{#nu} + #nu) - #bar{#nu}");
    
    nunubar_uncert->Draw("colz text");
    c->Print((plot_dump + "/nunubar_uncert.pdf").c_str());

    nuebar_uncert->Draw("colz text");
    c->Print((plot_dump + "/nuebar_uncert.pdf").c_str());

    uncert_diff->Draw("colz text");
    c->Print((plot_dump + "/uncert_diff.pdf").c_str());
    
    
  }
}