plot_NDvsFD_weights_RHC.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->SetLineColor(kGray+2);
  dummy->SetMarkerColor(kGray+2);
  dummy->SetLineColor(1);
  leg->AddEntry(dummy->Clone(), "ND #bar{#nu}_{#mu} Total MC", "l");
  dummy->SetLineColor(kNueSignalColor);
  leg->AddEntry(dummy->Clone(), "FD signal #bar{#nu}_{e} + WS", "l");

  leg->Draw();
}

void plot_NDvsFD_weights_RHC(const std::string fname_ND = "datamc_ND_numu_kinematics.root", const std::string fname_FD = "FDprediction_kinematics.root")
{  
  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 = 3;
  const std::string selNamesND[kNumSels] = {
        "numuND",
        "numuND",
        "numuND"
  };
  const std::string selNamesFD[kNumSels] = {
        "CORE",
        "PERI",
        "ALL"
  };
      
  IPrediction* predsND[kNumSels][kNumVars];

  osc::NoOscillations noosc;
  osc::OscCalcDumb mycalc;
  
  //pick up spectra/prediction objects from the input ND root file
  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
    const char* selName = selNamesND[selIdx].c_str();
    for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
      const char* varName = defs[varIdx].name.c_str();

      std::cout<<TString::Format("%s/spect_%s", selName, varName).Data()<<std::endl;
      predsND[selIdx][varIdx] = LoadFromFile<IPrediction>(fname_ND, TString::Format("%s/pred_%s", selName, 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;
            
      // need to pick equivalent variables
      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());
      }     
    }
  }  
  
  //loop over the different selector/cut levels to plot each one
  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
    const std::string selName = selNamesFD[selIdx];
    //TFile * file = new TFile(("Weights_"+selName+"_RHC.root").c_str(),"RECREATE");
    for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
      
      const std::string varName = defs[varIdx].name;
    
      //create a canvas   
      TCanvas *c2 =new TCanvas("c2","c2",500,400);      
      
      //absolute pad
      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.5;
      double legendxhigh=0.9;
      
      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
      
      const double kPOT = predND->Predict(&noosc).POT();
      
      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(1);
      hFD->SetLineColor(kNueSignalColor);
      
      std::cout << selName << " " << varName << ": " << hND->Integral() << std::endl;
      
      hND->GetXaxis()->CenterTitle();
      hND->GetYaxis()->CenterTitle();
      hND->GetXaxis()->SetDecimals();
      hND->GetYaxis()->SetDecimals();
      hND->GetYaxis()->SetTitleOffset(1.0);
      hND->GetXaxis()->SetLabelSize(.0);
      hND->GetYaxis()->SetLabelSize(.05);
      hND->SetMaximum(1.1*hND->GetMaximum());
      hND->GetYaxis()->SetTitle(TString::Format("Events"));
      
      if(varName=="recoQ2" || varName=="trueQ2"){
          hND->GetXaxis()->SetRangeUser(0,2);
      }

      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);
      //Preliminary();
      
      TLatex* selTitle = new TLatex(.935, .95, ("nominal "+selName).c_str());
      selTitle->SetNDC();
      selTitle->SetTextSize(2/30.);
      selTitle->SetTextAlign(32);
      selTitle->Draw();
      
      c2->cd();
      
      //ratio pad
      TPad *padRatio = new TPad("padRatio", "padRatio", 0, 0.02, 1, 0.33);
      padRatio -> SetTopMargin(0.1);
      padRatio -> SetBottomMargin(0.25);
      padRatio -> SetFillStyle(0);
      padRatio -> SetLeftMargin(0.15);
            
      padRatio -> Draw();
      padRatio -> cd();
    
      TH1* hRatio = (TH1*)hND -> 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(hFD);
      
      TH1* hWeight = (TH1*)hFD -> Clone("hWeight");
      hWeight -> 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(3.5/30.);
      hRatio->GetXaxis()->SetTitleSize(3.5/30.);
      hRatio->GetYaxis()->SetRangeUser(0,2.0);
      hRatio->SetMarkerSize(.1);
      hRatio->GetXaxis()->SetDecimals();
      hRatio->GetYaxis()->SetDecimals();
      hRatio->GetYaxis()->SetTitleOffset(0.4);
      hRatio->GetXaxis()->SetTitleOffset(1.05);
      hRatio->GetYaxis()->SetTitle("#bar{#nu}_{#mu} ND/Sig.+WS FD");

      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");
      padRatio -> RedrawAxis();
      padRatio -> Update();
      
      c2->Print(("plots/FDvsND_"+selName+"_"+varName+".png").c_str());
      
      // save further weights in weights file
      /*if(varName=="trueQ2"){
          hWeight->Write(("trueQ2_weight_"+selName).c_str());
      }
      if(varName=="PtP"){
          hWeight->Write(("PtP_weight_"+selName).c_str());
      }
      if(varName=="CosNumi"){
          hWeight->Write(("CosNumi_weight_"+selName).c_str());
      }*/
      
    } // end for varIdx
    //file->Close();
  } // end for selIdx
  
  // close weights file
  
}