multiverse_macro.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/Var.h"
#include "CAFAna/Core/MultiVar.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Core/SystShifts.h"

#include "CAFAna/Vars/Vars.h"
//#include "CAFAna/Vars/NueVarsExtra.h"

#include "NDAna/numucc_1Pi/NumuCC1PiVars.h"
#include "NDAna/numucc_1Pi/NumuCC1PiCuts.h"

#include "G4RwgtPlots/Geant4WeightVars.h"

#include "CAFAna/XSec/Multiverse.h"
#include "CAFAna/XSec/SystematicDef.h"
//#include "CAFAna/XSec/CutOptimization.h"

// for POT
#include "CAFAna/Analysis/Exposures.h"

// for plotting
#include "CAFAna/Analysis/Plots.h"

// to get datetime in plot names
#include "/nova/app/users/csweeney/xsec/code/Utilities/misc.h"

#include "StandardRecord/Proxy/SRProxy.h"

#include <vector>
#include <string>

#include "TFile.h"
#include "TH1.h"
#include "TCanvas.h"
#include "TGraphAsymmErrors.h"
#include "TString.h"
*/

#include "CAFAna/Core/Binning.h"
#include "CAFAna/Cuts/Cuts.h"

#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"

#include "CAFAna/Core/Var.h"

#include "CAFAna/Vars/Vars.h"

#include "CAFAna/XSec/Multiverse.h"
#include "CAFAna/XSec/SystematicDef.h"

#include "StandardRecord/Proxy/SRProxy.h"

#include <vector>
#include <string>

#include "TFile.h"
#include "TH1.h"
#include "TCanvas.h"
#include "TGraphAsymmErrors.h"                                                                        
#include "TString.h"      

#include "NDAna/numucc_1Pi/NumuCC1PiVars.h"
#include "NDAna/numucc_1Pi/NumuCC1PiCuts.h"

#include "G4RwgtPlots/Geant4WeightVars.h"

// for POT
#include "CAFAna/Analysis/Exposures.h"

// for plotting
#include "CAFAna/Analysis/Plots.h"

// to get datetime in plot names
#include "/nova/app/users/csweeney/xsec-TR/Utilities/misc.h"

// for kVtxX
#include "3FlavorAna/Vars/NueVarsExtra.h"


namespace ana
{   


  UpDownPair<TH1> 
  ToUpDownHist(Multiverse * mv, const TH1 * h_nominal)
  //Scavenged from CAFAna/XSec/CutOptimization.h
  {
    TH1 * h_up   = mv->GetPlusOneSigmaShift (h_nominal);
    TH1 * h_down = mv->GetMinusOneSigmaShift(h_nominal);
    // transform into 1 sided shift as the most conservative shift
    TH1 * ret = (TH1*) h_up->Clone();
    for(auto ibinx = 1; ibinx <= ret->GetNbinsX(); ibinx++) {
      for(auto ibiny = 1; ibiny <= ret->GetNbinsY(); ibiny++) {
        for(auto ibinz = 1; ibinz <= ret->GetNbinsZ(); ibinz++) {
          double up      = h_up      ->GetBinContent(ibinx, ibiny, ibinz);
          double down    = h_down    ->GetBinContent(ibinx, ibiny, ibinz);
          double nominal = h_nominal ->GetBinContent(ibinx, ibiny, ibinz);
          double abs_shift = std::max(std::abs(up - nominal), std::abs(down - nominal));
          ret->SetBinContent(ibinx, ibiny, ibinz, abs_shift + nominal);
        }
      }
    }
    return UpDownPair<TH1>{ret};        
  }


  void PlotDebug(const std::vector<TH1 *> universes,
                 TH1 * h_nominal,
                 std::string title, 
                 std::string name)
  //Scavenged from CAFAna/XSec/CutOptimization.h
  {
    TDirectory * tmp = gDirectory;
    TCanvas * c = new TCanvas();
    c->SetLeftMargin(0.15);
    c->SetRightMargin(0.01);
    c->cd();
    double max_val = -1e9;
    for(auto ihist = 0u; ihist < universes.size(); ihist++) {
      max_val = std::max(max_val, universes[ihist]->GetMaximum());
    }
    h_nominal->GetYaxis()->SetRangeUser(0, max_val * 1.2);
    h_nominal->SetTitle(title.c_str());
    h_nominal->Draw("hist");
    for(auto ihist = 0u; ihist < universes.size(); ihist++) {
      universes[ihist]->Draw("hist same");        
    }
    h_nominal->Draw("hist same");
      
    c->Print(name.c_str());
    tmp->cd();
  }


//  extern const Var kWeight([](const caf::SRProxy* sr){
//      //Just prints out weights for debugging purposes
//      if(sr->mc.nnu != 1 ) return -5.0;
//
//      std::cout << "foo" << std::endl;
//      auto weight = sr->mc.nu[0].geantWeights[0].piplusweight;
//      std::cout << weight << std::endl;
//
//      /*
//      int nUniv = weightTable.size();
//      std::cout << "nUniv: " << nUniv << std::endl;
//      */
//
//      /*
//      for(int i=0; i<nChannels; ++i){
//        std::cout << i << std::endl;
//        std::cout << weightTable[i].piplus << std::endl;
//
//      }
//      */
//      std::cout << "----------------" << std::endl;
//
//      return -5.0;
//  });







  Var kUnity([](const caf::SRProxy *sr){
      // To be used as "common_weight" for Multiverse constructor
      return 1.f;
    });


}


using namespace ana;  

void multiverse_macro()
{

  //std::string inName = "/nova/app/users/csweeney/multirw_files/9file.multirw.caf.root";
  std::string inName = "/pnfs/nova/scratch/users/csweeney/final_grid_prod5_g4rwgt_miniprod/*";

  SpectrumLoader* loader = new SpectrumLoader(inName);

  //std::string outName =  "test.root";
  //TFile *outFile      = new TFile(outName.c_str(), "RECREATE");

  const SystShifts * shift = new SystShifts(kNoShift);
  const Var * common_weight = new Var(kUnity);



  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  std::vector<std::vector<Var>> weights_vec{
    GetkGeantPiplusWeights(), GetkGeantPiplusWeights(), GetkGeantPiplusWeights(),
    GetkGeantPiminusWeights(), GetkGeantPiminusWeights(), GetkGeantPiminusWeights(),
    GetkGeantAllPionsWeights(), GetkGeantAllPionsWeights(), GetkGeantAllPionsWeights()
  };
  std::vector<std::string> myLabel_vec{
    "Vertex in X Direction (cm)", "Vertex in Y Direction (cm)", "Vertex in Z Direction (cm)",
    "Vertex in X Direction (cm)", "Vertex in Y Direction (cm)", "Vertex in Z Direction (cm)",
    "Vertex in X Direction (cm)", "Vertex in Y Direction (cm)", "Vertex in Z Direction (cm)"
  };
  std::vector<std::string> myTitle_vec{
    "Piplus 100 universes",  "Piplus 100 universes",  "Piplus 100 universes",
    "Piminus 100 universes", "Piminus 100 universes", "Piminus 100 universes",
    "Total 100 universes",   "Total 100 universes",   "Total 100 universes"
  };
  std::vector<std::string> myName_vec  = {
    "eff_piplus_vtxX",  "eff_piplus_vtxY",  "eff_piplus_vtxZ", 
    "eff_piminus_vtxX", "eff_piminus_vtxY", "eff_piminus_vtxZ", 
    "eff_total_vtxX",   "eff_total_vtxY",   "eff_total_vtxZ"
  }; 
  std::vector<Binning> myBins_vec{
    Binning::Simple(10, -200, 200), Binning::Simple(10, -200, 200), Binning::Simple(10, 0, 1270),
    Binning::Simple(10, -200, 200), Binning::Simple(10, -200, 200), Binning::Simple(10, 0, 1270),
    Binning::Simple(10, -200, 200), Binning::Simple(10, -200, 200), Binning::Simple(10, 0, 1270)
  };
  std::vector<const Var *> kMyVar_vec{
    new Var(kVtxX), new Var(kVtxY), new Var(kVtxZ),
    new Var(kVtxX), new Var(kVtxY), new Var(kVtxZ),
    new Var(kVtxX), new Var(kVtxY), new Var(kVtxZ)
  };    
  const Cut * kSignal = new Cut(CutFromNuTruthCut(kIsNumuCC1Pi_NT && kTrueDetector_NT));
  //const Cut * kSignal = new Cut(CutFromNuTruthCut(kTrueDetector_NT) );


  //-----

  std::vector<Multiverse *> myMV_vec;
  std::vector<Spectrum *> mySpec_vec;



  for(uint i=0; i<myName_vec.size(); ++i){

    std::string myLabel = myLabel_vec[i];
    Binning myBins = myBins_vec[i];
    const Var * kMyVar = kMyVar_vec[i];
    std::vector<Var> weights = weights_vec[i];

    //    std::string myName = myName_vec[i];



    const HistAxis * axis =  new HistAxis(myLabel,
                                          myBins,
                                          *kMyVar);
  
  
    Multiverse * myMV = new Multiverse(*loader,
                                       *axis,
                                       *kSignal,
                                       *shift,
                                       weights,
                                       *common_weight);
  
  
    Spectrum * sSig = new Spectrum(myLabel, myBins, 
                                   *loader, *kMyVar, 
                                   *kSignal);
    myMV_vec.push_back(myMV);
    mySpec_vec.push_back(sSig);
  }

  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~




  loader->Go();


  for(uint i=0; i<myMV_vec.size(); ++i){
  
    Multiverse * myMV = myMV_vec[i];
    Spectrum * sSig = mySpec_vec[i];

    std::string myName = myName_vec[i];
    std::string myLabel = myLabel_vec[i];
    std::string myTitle = myTitle_vec[i];


    //double pot = kAna2020FHCPOT;
    double pot = sSig->POT();
  
    //TDirectory * dir = outFile->mkdir("vtx_x");
    //TDirectory * nom_dir = dir->mkdir("Nominal");
    //fNominalSignal->SaveTo(nom_dir->mkdir("Signal"));
    //TDirectory * mv_dir = dir->mkdir("mv_geant4");
    //myMV->SaveTo( mv_dir->mkdir("fMVSignal") );
  
    //------------------------------------------
  
    //TH1 * hSig =(TH1*) fNominalSignal->ToTH1(pot);
    //std::vector<UpDownPair<TH1> > hSystSignal;
    //hSystSignal.push_back(ToUpDownHist(myMV,hSig) );
    TH1 * hSig = sSig->ToTH1(pot);
  
    /*
    hSig->SetAxisRange(0., 2.0*hSig->GetMaximum(), "Y" );
    const std::vector<TH1*> mv_signal_hists = myMV->GetUniversesHist(kCyan); 
    PlotDebug(mv_signal_hists, hSig, "this is a plot", "plot.png");
    */
    
    Spectrum * sUp   = myMV->GetPlusOneSigmaShift(sSig);
    TH1 * hUp = sUp->ToTH1(pot);
    std::vector<TH1*> hUp_vec{hUp}; // for some reason syst plotting function wants vector of TH1s
    //hUp->SetLineColor(kRed);
  
    Spectrum * sDown = myMV->GetMinusOneSigmaShift(sSig);
    TH1 * hDown  = sDown->ToTH1(pot);
    std::vector<TH1*> hDown_vec{hDown}; 
    //hDown->SetLineColor(kBlue);
  
    // Get fractional uncertainty (take larger of +/-)
    int nBins = hSig->GetNbinsX();
    double xLow = hSig->GetXaxis()->GetXmin();
    double xHigh = hSig->GetXaxis()->GetXmax();
  
    TH1D * fracUncert = new TH1D("frac", "", nBins, xLow, xHigh);
    for(int i=1; i <= nBins; ++i){
      double up   = hUp->GetBinContent(i);
      double down = hDown->GetBinContent(i);
      double nom  = hSig->GetBinContent(i);
  
      double abs_shift = std::max( std::abs(up - nom), std::abs(down - nom) );
      double frac = abs_shift / nom;
      //std::cout << "{ " << up << ", " << down << " }" << std::endl;
      //std::cout << nom << " : "  << abs_shift << " : " << frac << std::endl;
      fracUncert->SetBinContent(i, frac);
    }
  
    fracUncert->GetXaxis()->SetTitle(myLabel.c_str());
  
  
    /*
    TCanvas * canv2 = new TCanvas("canv2", "canv2", 800, 600);
    fracUncert->Draw();
    */
    
    TGraphAsymmErrors * graph = PlotWithSystErrorBand(hSig, hUp_vec, hDown_vec);  
    graph->SetTitle(myTitle.c_str());
  
    //TCanvas * canv = new TCanvas("canv", "canv", 800, 600);
  
    hSig->SetLineColor(1);
    hSig->SetLineWidth(2);
    /*
    hSig->Draw("hist ][");
    graph->Draw("e2 same");
    hSig->Draw("hist ][ same");
  
    canv->SaveAs("error_band.png");
    */
  
    TString rCname = "rC";
    TCanvas* rC = new TCanvas(rCname, rCname);
  
    rC -> SetBottomMargin(0.);
    double Spl = 0.3;
    TPad* P1 = new TPad( "Temp_1", "", 0.0, Spl, 1.0, 1.0, 0 );
    TPad* P2 = new TPad( "Temp_2", "", 0.0, 0.0, 1.0, Spl, 0 );
    P2 -> SetRightMargin (.03);
    P2 -> SetTopMargin   (.00);
    P2 -> SetBottomMargin(.3);
    P2 -> SetLeftMargin  (.1);
    P2 -> Draw();
    P1 -> SetRightMargin (.03);
    P1 -> SetLeftMargin  (.1);
    P1 -> SetTopMargin   (.10);
    P1 -> SetBottomMargin(.00);
    P1 -> Draw();
    // Set some label sizes.                                                                                
    double Lb1 = 0.07;
    double Lb2 = 0.15;
    // --- First, draw the fracUncert so cd onto Pad2                                                           
    P2 -> cd();
  
    // Set axis ranges etc.                                                                                 
    fracUncert->GetYaxis()->SetTitleSize( Lb2 );
    fracUncert->GetYaxis()->SetTitleOffset(0.4);
    fracUncert->GetYaxis()->SetLabelSize( Lb2 );
    fracUncert->GetXaxis()->SetTitleSize( Lb2 );
    fracUncert->GetXaxis()->SetLabelSize( Lb2 );
    fracUncert->Draw("");
  
    //fracUncert->Draw("axis same");
    P1->cd();
    // graph goes on P1    
    hSig->SetTitle(myTitle.c_str());
    hSig->GetYaxis()->SetTitleSize( Lb1 );
    hSig->GetYaxis()->SetLabelSize( Lb1 );
    hSig->GetYaxis()->SetTitleOffset( 0.7 );
    // Remove the x axis labels                                                                         
    hSig->GetXaxis()->SetLabelSize  (0 );
    hSig->GetXaxis()->SetLabelOffset(99);
    hSig->SetLineColor(4);
    hSig->DrawCopy("hist");
  
    hSig->Draw("hist ][");
    graph->Draw("e2 same");
    hSig->Draw("hist ][ same");
  
    rC->SaveAs(("./plots/multiverse/" + DateTime() + myName + ".png").c_str());
  
  }  
  

}