test_genieweights.C

Go to the documentation of this file.

#include "CAFAna/Core/Binning.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/SpectrumLoaderBase.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Vars/GenieWeights.h"

#include "TCanvas.h"
#include "TH1.h"
#include "TFile.h"
#include "TLegend.h"
#include "TLine.h"
#include "TLatex.h"
#include "THStack.h"
#include "3FlavorAna/Vars/NumuVars.h"

using namespace ana;

void test_genieweights(){
  
  TString filename = "test_genieweight_v0.root";
  TString imgname ="plot_test_genieweight.png";
  //Options from Vars/GenieVars.h: 
  //kTuftsWeightCC CCQE CCnonQE; kUMNWeightCC CCQE CCnonQE; kRPAWeightCCQE 
  const Var kGenieWeight = kTuftsWeightCC;

  SpectrumLoader loaderNDMC 
       ("/nova/prod/mc/S15-05-22/genie/nd/decaf/numu_contain/000108/**/neardet_genie_fhc_nonswap_nogenierw_2000_r000**_s*_c000*");
  
  loaderNDMC.SetSpillCut(kStandardSpillCuts);

  const Cut numuCC([](const caf::SRProxy* sr)
                   { 
                     if(sr->mc.nnu == 0) return false;
                     assert(sr->mc.nnu == 1);
                     return (sr->mc.nu[0].iscc && abs(sr->mc.nu[0].pdg)==14);
                   });
  //*******************
  const HistAxis axisEhad(
                          "E_{had} (GeV)", Binning::Simple(50, 0, 1), kHadE);
  //*******************
  Spectrum sp_QE_noW  (loaderNDMC, axisEhad, numuCC && kIsQE, kNoShift);
  Spectrum sp_RES_noW (loaderNDMC, axisEhad, numuCC && kIsRes,kNoShift);
  Spectrum sp_DIS_noW (loaderNDMC, axisEhad, numuCC && kIsDIS,kNoShift);
  Spectrum sp_Coh_noW (loaderNDMC, axisEhad, numuCC && kIsCoh,kNoShift);
  
  Spectrum sp_QE_Wei  (loaderNDMC, axisEhad, numuCC && kIsQE, kNoShift,kGenieWeight);
  Spectrum sp_RES_Wei (loaderNDMC, axisEhad, numuCC && kIsRes,kNoShift,kGenieWeight);
  Spectrum sp_DIS_Wei (loaderNDMC, axisEhad, numuCC && kIsDIS,kNoShift,kGenieWeight);
  Spectrum sp_Coh_Wei (loaderNDMC, axisEhad, numuCC && kIsCoh,kNoShift,kGenieWeight);
  
  loaderNDMC.Go();

  TFile * file = new TFile (filename,"recreate");
  sp_QE_noW. SaveTo(file, "sp_QE_noweight");
  sp_RES_noW.SaveTo(file, "sp_RES_noweight");
  sp_DIS_noW.SaveTo(file, "sp_DIS_noweight");
  sp_Coh_noW.SaveTo(file, "sp_Coh_noweight");

  sp_QE_Wei. SaveTo(file, "sp_QE_weight");
  sp_RES_Wei.SaveTo(file, "sp_RES_weight");
  sp_DIS_Wei.SaveTo(file, "sp_DIS_weight");
  sp_Coh_Wei.SaveTo(file, "sp_Coh_weight");


  //return;

  //****************************************************
  double thePOT = 1E20;

  auto hQEno = sp_QE_noW.ToTH1(thePOT);
  auto hRESno = sp_RES_noW.ToTH1(thePOT);
  auto hDISno = sp_DIS_noW.ToTH1(thePOT);
  auto hCohno = sp_Coh_noW.ToTH1(thePOT);

  hQEno->SetLineColor(kMagenta);
  hRESno->SetLineColor(kGreen+2);
  hDISno->SetLineColor(kRed);
  hCohno->SetLineColor(kBlue);

  hQEno->SetLineStyle(2);  
  hRESno->SetLineStyle(2);  
  hDISno->SetLineStyle(2);  
  hCohno->SetLineStyle(2);

  THStack * hsno = new THStack ("noW","Stacked histograms; E_{had} (GeV); Events/1E20 POT");
  hsno->Add(hCohno,"histo");
  hsno->Add(hQEno,"histo");
  hsno->Add(hRESno,"histo");
  hsno->Add(hDISno,"histo");  

  auto hQE  = sp_QE_Wei.ToTH1(thePOT);
  auto hRES = sp_RES_Wei.ToTH1(thePOT);
  auto hDIS = sp_DIS_Wei.ToTH1(thePOT);
  auto hCoh = sp_Coh_Wei.ToTH1(thePOT);

  hQE->SetMarkerColor(kMagenta);
  hRES->SetMarkerColor(kGreen+2);
  hDIS->SetMarkerColor(kRed);
  hCoh->SetMarkerColor(kBlue);

  hQE->SetMarkerStyle(29);  
  hRES->SetMarkerStyle(29);  
  hDIS->SetMarkerStyle(29);  
  hCoh->SetMarkerStyle(29);

  hQE->SetMarkerSize(2);  
  hRES->SetMarkerSize(2);  
  hDIS->SetMarkerSize(2);  
  hCoh->SetMarkerSize(2);


  THStack * hsW = new THStack ("W","Stacked histograms; E_{had} (GeV); Events/1E20 POT");
  hsW->Add(hCoh,"hist p");
  hsW->Add(hQE,"hist p");
  hsW->Add(hRES,"hist p");
  hsW->Add(hDIS,"hist p");

  hsno->SetMaximum(max(hsno->GetMaximum(),hsW->GetMaximum()));
  hsno->Draw();
  hsW->Draw("same");

  TLegend* leg = new TLegend(0.6,0.6,0.89,0.89);
  leg->AddEntry(hDIS,"DIS (weighed)", "p");
  leg->AddEntry(hRES,"RES (weighed)", "p");
  leg->AddEntry(hQE, " QE (weighed)", "p");
  leg->AddEntry(hCoh,"Coh (weighed)", "p");
  leg->Draw("same");
  gPad->Print(imgname);

  return;
}