plot_NDvsFD_REW.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/SpectrumLoaderBase.h"
#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Decomp/ProportionalDecomp.h"
#include "3FlavorAna/Decomp/NumuDecomp.h"
#include "CAFAna/Extrap/ModularExtrap.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "CAFAna/Vars/Vars.h"

#include "OscLib/OscCalcPMNSOpt.h"
#include "OscLib/OscCalcDumb.h"
#include "CAFAna/Analysis/Style.h"

//#include "CAFAna/nue/SecondAna/draw_plots_util.h"

using namespace ana;

#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TGaxis.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TLine.h"
#include "TSystem.h"
#include <iostream>
#include <iomanip>
#include <fstream>

//make a pretty legend 
void Legend(double x0, double y0, double x1, double y1)
{
  // Doesn't handle log-y well
  //  TLegend* leg = AutoPlaceLegend(x1-x0, y1-y0);
  TLegend* leg = new TLegend(x0, y0, x1, y1);
  leg->SetFillStyle(0);

  TH1* dummy = new TH1F("", "", 1, 0, 1);
  TH1* dummyFill = new TH1F("", "", 1, 0, 1);
  dummyFill->SetLineColor(kTotalMCColor);
  dummyFill->SetFillColor(kTotalMCErrorBandColor);

  dummy->SetMarkerStyle(kFullCircle);

  dummy->SetLineColor(kGray+2);
  dummy->SetMarkerColor(kGray+2);
  dummy->SetLineColor(kTotalMCColor);
  leg->AddEntry(dummy->Clone(), "ND #nu_{#mu} Total MC", "l");
  
  dummy->SetLineColor(kNueSignalColor);
  leg->AddEntry(dummy->Clone(), "FD signal #nu_{e}", "l");

  leg->Draw();
}

void plot_NDvsFD_REW(const std::string weight_name, const std::string fname_ND, const std::string fname_FD)
{
  const double kPOT = 8.09e20;
  
  struct HistDef
  {
    std::string name;
  };

  //which variables do you actually want to look at?
  const int kNumVars = 5;
  const HistDef defs[kNumVars] = {
    "trueQ2",
    "recoQ2",
    "trueW2",
    "PtP",
    "CosNumi"
  };

  //which cut tiers do you want to look at?
  const int kNumSels = 1;
  const std::string selNamesND[kNumSels] = {
        "numuND"
  };
    
  const std::string selNamesFD[kNumSels] = {
        "AllSamples"
  };

  IPrediction* predsND[kNumSels][kNumVars];

  osc::NoOscillations noosc;
  osc::OscCalcPMNSOpt mycalc;
  ResetOscCalcToDefault(&mycalc);
    
  mycalc.SetDmsq32(+fabs(mycalc.GetDmsq32()));
  mycalc.SetdCP(3*TMath::Pi()/2);
  mycalc.SetTh23(asin(sqrt(0.404)));
  
  //pick up spectra/prediction objects from the input root file
  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
    const char* selName = selNamesND[selIdx].c_str();
    const char* wName = weight_name.c_str();
    
    for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
      const char* varName = defs[varIdx].name.c_str();

      std::cout<<TString::Format("%s/spect_%s_%s", selName, wName, varName).Data()<<std::endl;
      
      predsND[selIdx][varIdx] = LoadFromFile<IPrediction>(fname_ND, TString::Format("%s/pred_%s_%s", selName, wName, varName).Data()).release();
      
      if(selIdx == 0 && varIdx == 0)
        std::cout << "ND MC POT: "
                  << predsND[selIdx][varIdx]->Predict(&noosc).POT() << std::endl;

    } // end for varIdx
  } // end for selIdx
    
    
  TFile * file_no = new TFile (fname_FD.c_str(),"READ");
    TDirectory * dpred_no[kNumSels][kNumVars];
    for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
        const std::string selName = selNamesFD[selIdx];
        for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
            const std::string varName = defs[varIdx].name;
            if(varName == "recoQ2") {
                dpred_no[selIdx][varIdx] = file_no->GetDirectory(("prediction/nue_pred_noextrap_"+selName+"_trueQ2").c_str());
            } else if(varName == "CosNumi"){
                dpred_no[selIdx][varIdx] = file_no->GetDirectory(("prediction/nue_pred_noextrap_"+selName+"_CosTheta").c_str());
            } else {
                dpred_no[selIdx][varIdx] = file_no->GetDirectory(("prediction/nue_pred_noextrap_"+selName+"_"+varName).c_str());
            }       
        }
    }
  
  TH1* hFullNDHists[kNumVars-1];
  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
    const std::string selName = selNamesFD[selIdx];
    for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
      
      const std::string varName = defs[varIdx].name;

      //play with the binning choice for given variables:  
      int rebinFactor=1;
    
      //create a canvas   
      TCanvas *c2 =new TCanvas("c2","c2",500,400);
      
      
      TPad *padAbs = new TPad("padAbs", "padAbs", 0, 0.3, 1, 1.0);
      padAbs -> SetBottomMargin(0); 
      padAbs->SetFillStyle(0);
      padAbs -> SetLeftMargin(0.15);
      
      padAbs -> Draw();
      padAbs -> cd();

      double legendxlow=0.55;
      double legendxhigh=0.95;
      
      if(varName=="trueW2"||varName=="CosNumi"){legendxlow=0.2; legendxhigh=0.5;}

      //Pull down the variable we want

      IPrediction* predND = predsND[selIdx][varIdx];
      auto predFD = ana::LoadFrom<IPrediction> (dpred_no[selIdx][varIdx]);

      //Get components of the prediction and store them in histograms
      
      TH1* hFD = predFD->PredictComponent(&mycalc,
                                          Flavors::kNuMuToNuE,
                                          Current::kCC,
                                          Sign::kBoth).ToTH1(kPOT);
                                          
      TH1* hAux = predND->Predict(&noosc).ToTH1(kPOT);                                    
                                                 
      TH1* hND = predND->Predict(&noosc).ToTH1(kPOT*hFD->Integral()/hAux->Integral());
      
      hND->SetLineColor(kTotalMCColor);
      
      hFD->SetLineColor(kNueSignalColor);
      hND->Rebin(rebinFactor);
      
      std::cout << selName << " " << varName << ": " << hND->Integral() << std::endl;
      
      hND->GetXaxis()->CenterTitle();
      hND->GetYaxis()->CenterTitle();
      hND->GetXaxis()->SetDecimals();
      hND->GetYaxis()->SetDecimals();
      hND->GetYaxis()->SetTitleOffset(1.15);
      hND->GetXaxis()->SetLabelSize(.0);
      hND->GetYaxis()->SetLabelSize(.05);
      hND->SetMaximum(1.1*TMath::Max(hND->GetMaximum(), hFD->GetMaximum()));
      hND->GetYaxis()->SetTitle(TString::Format("Events"));
      
      if(varName=="recoQ2" || varName=="trueQ2"){
          hND->GetXaxis()->SetRangeUser(0,2);
      }

      if(varName=="recoE"){
        hND->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
      }
      hND->GetXaxis()->SetNdivisions(405,kFALSE);
      hND->SetMinimum(.001);
      
      hND ->DrawCopy("hist");
  
      hND->Draw("hist same");
      hFD->Draw("hist same");

      //Fix up the axis
      padAbs->RedrawAxis();
      padAbs->Update();
      
      //Add the legend and preliminary
      Legend(legendxlow, .70, legendxhigh, .85);
      std::string wTitle = weight_name+" reweighted";
      
      TLatex* selTitle = new TLatex(.935, .95, wTitle.c_str());
      selTitle->SetNDC();
      selTitle->SetTextSize(2/30.);
      selTitle->SetTextAlign(32);
      selTitle->Draw();
      
      c2->cd();
      TPad *padRatio = new TPad("padRatio", "padRatio", 0, 0.05, 1, 0.3);
      padRatio -> SetTopMargin(0);
      padRatio -> SetBottomMargin(0.25);
      padRatio -> SetFillStyle(0);
      padRatio -> SetLeftMargin(0.15);
      
      padRatio -> Draw();
      padRatio -> cd();
    
      TH1* hRatio = (TH1*)hFD -> Clone("hRatio");
      
      double b=5;
      if(varName=="trueW2") {b = 2.0;}
      if(varName=="trueQ2") {b = 2.0;}
      if(varName=="recoQ2") {b = 2.0;}
      if(varName=="PtP" || varName=="CosNumi") {b = 1.0;}
      
      TLine* lOne = new TLine(0,1.0,b,1.0);
      lOne -> SetLineStyle(2);
      
      hRatio -> Divide(hND);
      
      hRatio->GetXaxis()->CenterTitle();
      hRatio->GetYaxis()->CenterTitle();
      hRatio->GetYaxis()->SetLabelSize(4.0/30.);
      hRatio->GetXaxis()->SetLabelSize(4.0/30.);
      hRatio->GetXaxis()->SetLabelOffset(.04);
      hRatio->GetXaxis()->SetTickSize(.06);
      hRatio->GetYaxis()->SetTitleSize(4.5/30.);
      hRatio->GetXaxis()->SetTitleSize(4.5/30.);
      hRatio->GetYaxis()->SetRangeUser(0.5,1.5);
      hRatio->SetMarkerSize(.1);
      hRatio->GetXaxis()->SetDecimals();
      hRatio->GetYaxis()->SetDecimals();
      hRatio->GetYaxis()->SetTitleOffset(0.4);
      hRatio->GetYaxis()->SetTitle("Sig. FD/#nu_{#mu} ND");

      if(varName=="recoE"){
        hRatio->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
      }
      if(varName=="recoQ2"){
        hRatio->GetXaxis()->SetTitle("ND reco Q^{2} vs. FD true Q^{2} (GeV^{2})");
      }
      
      if(varName=="recoQ2" || varName=="trueQ2"){
          hRatio->GetXaxis()->SetRangeUser(0,2);
      }
      
      hRatio->GetYaxis()->SetNdivisions(502,kFALSE);
      hRatio->GetXaxis()->SetNdivisions(405,kFALSE);
      
      hRatio-> Draw("hist");
      lOne -> Draw("same");
      
      c2->Print(("plots/FDvsND_REW_"+weight_name+"_"+selName+"_"+varName+".png").c_str());
      
    } // end for varIdx
  } // end for selIdx
}